We present a general, rigorous theory of Lee-Yang zeros for models with first-order phase transitions that admit convergent contour expansions. We derive formulas for the positions and the density of the zeros. In particular, we show that for models without symmetry, the curves on which the zeros lie are generically not circles, and can have topologically nontrivial features, such as bifurcation. Our results are illustrated in three models in a complex field: the low-temperature Ising and Blume-Capel models, and the q-state Potts model for q large enough. PACS numbers: 05.50.+q, 05.70.Fh, 64.60.Cn, 75.10.Hk Almost half a century ago, in two classic papers [1], Lee and Yang studied the zeros of the Ising partition function in the complex magnetic field plane, showed rigorously that the zeros lie on the unit circle, and proposed a program to analyze phase transitions in terms of these zeros. A decade later, Fisher [2] extended the study of the Ising partition function zeros to the complex temperature plane. Since that time, there have been numerous studies, both exact and numerical, of the Lee-Yang and Fisher zeros in a wide variety of models [3]. However, with a few notable exceptions [4], remarkably little progress has been made in extending the rigorous LeeYang program. This is due to the fact that rigorous statistical mechanics has relied almost exclusively on probabilistic techniques which fail in a complex parameter space. In this Letter, we adapt complex extensions [5,6,7] of Pirogov-Sinai theory [8] to realize the Lee-Yang program in a general class of models with first-order phase transitions.The purpose of this work is threefold. First, it is of interest to establish the mathematical foundation of a program that has been so central to statistical physics. Second, our theory gives a novel physical interpretation of the existence and position of partition function zeros by relating them to the phase coexistence lines in the complex plane. Finally, from a practical viewpoint, our theory provides a framework for the interpretation of numerical data by allowing explicit, rigorous computation of the position of the zeros. Indeed, we find rigorous results which clarify many ambiguities in published data. Specifically, in models without an underlying symmetry, we prove that the zeros generically do not lie on circles, even in the thermodynamic limit. This applies, in particular, to the Blume-Capel and Potts models in complex magnetic fields; see [9,10] for heuristic studies of these models. We also prove that the curves defined by the asymptotic positions of the zeros can have topologically nontrivial features, such as bifurcation and coalescence, and show that these features correspond to triple (or higher) points in the complex phase diagram.The results to be stated next are rather technical. Roughly speaking, they say that, for models with a convergent contour expansion, the partition function can be written in the form (1) and its zeros are given as solutions to equations (2) Main Result: Consider a d-di...