Melting in two spatial dimensions, as realized in thin films or at interfaces, represents one of the most fascinating phase transitions in nature, but it remains poorly understood. Even for the fundamental hard-disk model, the melting mechanism has not been agreed on after fifty years of studies. A recent Monte Carlo algorithm allows us to thermalize systems large enough to access the thermodynamic regime. We show that melting in hard disks proceeds in two steps with a liquid phase, a hexatic phase, and a solid. The hexatic-solid transition is continuous while, surprisingly, the liquid-hexatic transition is of first-order. This melting scenario solves one of the fundamental statistical-physics models, which is at the root of a large body of theoretical, computational and experimental research. Generic two-dimensional particle systems cannot crystallize at finite temperature [1][2][3] because of the importance of fluctuations, yet they may form solids [4]. This paradox has provided the motivation for elucidating the fundamental melting transition in two spatial dimensions. A crystal is characterized by particle positions which fluctuate about the sites of an infinite regular lattice. It has long-range positional order. Bond orientations are also the same throughout the lattice. A crystal thus possesses long-range orientational order. The positional correlations of a two-dimensional solid decay to zero as a power law at large distances. Because of the absence of a scale, one speaks of "quasi-long range" order. In a two-dimensional solid, the lattice distortions preserve long-range orientational order [5], while in a liquid, both the positional and the orientational correlations decay exponentially.Besides the solid and the liquid, a third phase, called "hexatic", has been discussed but never clearly identified in particle systems. The hexatic phase is characterized by exponential positional but quasi-long range orientational correlations. It has long been discussed whether the melting transition follows a one-step first-order scenario between the liquid and the solid (without the hexatic) as in three spatial dimensions Two-dimensional melting was discovered [4] in the simplest particle system, the hard-disk model. Hard disks (of radius σ) are structureless and all configurations of nonoverlapping disks have zero potential energy. Two isolated disks only feel the hard-core repulsion, but the other disks mediate an entropic "depletion" interaction (see, e.g., [19]). Phase transitions result from an "order from disorder" phenomenon: At high density, ordered configurations can allow for larger local fluctuations, thus higher entropy, than the disordered liquid. For hard disks, no difference exists between the liquid and the gas. At fixed density η, the phase diagram is independent of temperature T = 1/k B β, and the pressure is proportional to T , as discovered by D. Bernoulli in 1738. Even for this basic model, the nature of the melting transition has not been agreed on.The hard-disk model has been simulated with the...