Evidence of ideal excitonic insulator in bulk MoS2 under pressure

Spontaneous condensation of excitons is a long sought phenomenon analogous to the condensation of Cooper pairs in a superconductor. It is expected to occur in a semiconductor at thermodynamic equilibrium if the binding energy of the excitons---electron (e) and hole (h) pairs interacting by Coulomb force---overcomes the band gap, giving rise to a new phase: the 'excitonic insulator' (EI). Transition metal dichalcogenides are excellent candidates for the EI realization because of reduced Coulomb screening, and indeed a structural phase transition was observed in few-layer systems. However, previous work could not disentangle to which extent the origin of the transition was in the formation of bound excitons or in the softening of a phonon. Here we focus on bulk MoS2 and demonstrate theoretically that at high pressure it is prone to the condensation of genuine excitons of finite momentum, whereas the phonon dispersion remains regular. Starting from first-principles many-body perturbation theory, we also predict that the self-consistent electronic charge density of the EI sustains an out-of-plane permanent electric dipole moment with an antiferroelectric texture in the layer plane: At the onset of the EI phase, those optical phonons that share the exciton momentum provide a unique Raman fingerprint for the EI formation. Finally, we identify such fingerprint in a Raman feature that was previously observed experimentally, thus providing direct spectroscopic confirmation of an ideal excitonic insulator phase in bulk MoS2 above 30 GPa.

print in a Raman feature that was previously observed experimentally, thus providing direct spectroscopic confirmation of an ideal excitonic insulator phase in bulk MoS 2 above 30 GPa.
The long-sought excitonic insulator (EI) is a permanent Bose-Einstein condensate of excitons in the absence of optical excitation, hosted in a narrow-gap semiconductor or a semimetal [1][2][3][4] . As the exciton condensate shares similarities with the superconductor ground state 5 , it may exhibit macroscopic quantum coherence and exotic low-energy excitations [6][7][8][9][10][11][12] . These intriguing features are linked to the arbitrariness of the phase of the condensate wave function, ϕ (defined in Eq. 2 below): whereas in the superconductor this phase degeneracy is protected by the conservation of electronic charge, in the EI it is contingent on the preservation of excitons 7,11 , and hence lifted by those terms in the Hamiltonian that annihilate or create e-h pairs. This is the case of e-phonon 13 and spin-orbit 14 interactions, which pin ϕ while hybridizing conduction and valence bands [remarkably, spin-orbit coupling provides excitonic insulators with topological properties 14 ]. So far, the most accomplished EIs were realized in bilayer heterostructures in the presence of a magnetic field, requiring both low temperature and complex engineering to maximize the impact of e-h correlations as well as the degeneracy of ϕ 10,15 . A related concept aims to achieve the temporary condensation of indirect excitons, made of spatially separated e and h, through the optical pumping of artificial bilayers designed to maximize the exciton lifetime [16][17][18] . electronic, with only small adjustments of the lattice 4,28 .
Here, we follow an early suggestion by Hromadová et al. 29 and focus on bulk MoS 2 under hydrostatic pressure [29][30][31][32] . We use many-body perturbation theory from first principles 33,34 to demonstrate that MoS 2 is unstable against exciton condensation but stable against lattice distortion. Bulding a self-consistent effective-mass model on top of ab initio calculations, we show that the true ground state is an ideal, anti-ferrolectric EI with a distinctive Raman fingerprint that has already been observed 35 .
In bulk MoS 2 , the pressure (P ) closes the indirect gap, G, between the top of the filled valence band-located at the center of the Brillouin zone (Γ point), and the bottom of the sixdegenerate valleys of the empty conduction band-placed at Λ points (approximately midway between Γ and K, see Fig. 1C for P = 34 GPa). The energy landscape along one of the ΓΛ cuts (sketched in Fig. 1A) favours the Coulomb binding of an e, located at Λ, with a h, placed at Γ, creating an exciton of finite momentum |q| = ΓΛ and binding energy E b . Whereas ordinarily E b < G, it may occur that E b > G above a critical pressure, a condition that makes the semiconductor unstable against the condensation of excitons. This is actually the case, as we show below from first principles.
So far, ultra high pressure has been used as a handle to make MoS 2 superconducting 32 (at P ∼ 90 GPa), though the pairing mechanism remains unclear [36][37][38] . The putative EI must be searched at lower pressure (P ∼ 25 GPa), close to the semiconductor-semimetal transition that was observed by several groups 30,31,[39][40][41] . Near this boundary, theory 29,42 -including our own  calculations (SI Appendix, Fig. S1)-predicts an isostructural transition from the 2H c (Fig. 1B) to the 2H a (SI Appendix, Fig. S2) phase, which does not affect the crystal space group D 4 6h , as the two structures transform into each other through the sliding of the layers in the unit cell (the layer unit is made of one Mo and two S atoms, represented respectively by violet and yellow balls in the sketch of Fig. 1B). Raman and x-ray spectroscopic observations 30,32,35,40,43 suggest that 2H c and 2H a phases coexist in diamond-anvil cells, in a range that varies between 25 and 50 GPa in powders but has narrower extension (∼ 4 GPa) in single crystals. Importantly, we find that both 2H c and 2H a polytypes experience a similar excitonic instability-unrelated to the structural transition, as the electronic bands of the two phases are basically identical close to the Fermi energy. Below, we discuss the 2H c stacking and leave the analysis of 2H a to the SI Appendix, Figs. S3 and S4.

Results
The indirect gap of 2H c -MoS 2 is sensitive to pressure, as its value drops from 1.31 eV at P = 0 ( Fig. 2A) to only 9 meV at P = 34 GPa (Fig. 2C), close to the semimetal limit. With respect to the accurate band structure calculated within the GW approximation (circles in Fig. 2A to C, see Methods), density functional theory (triangles) underestimates the gap of about 0.4 eV at P = 0.
However, as pressure reduces the out-of-plane lattice parameter c (SI Appendix, Fig. S1), forcing sulfur orbitals belonging to adjacent layers to overlap 44 , virtual e-h pairs start tunnelling among layers, screening effectively Coulomb interaction at long wavelength. This reduces the GW energy correction to DFT bandgap, as evident in Fig. 2C. Consistently, the conduction band increases its dispersion along the k z direction (Fig. 2F), as well as the other axes of the effective mass tensor (D to F) Band dispersion of conduction and valence bands close to Λ and Γ points, respectively, for P = 0 GPa (black colour), 25 GPa (red), and 34 GPa (blue). In panels E and F the conduction band has been rigidly translated by the wave vector − ΓΛ. GW predictions (dots) are shown together with effective-mass fits (curves). The directions of the cuts (shown as red arrows in the Brillouin zone) are the principal axes of the effective-mass tensor, two being in the k z = 0 plane (panels D and E) and one being parallel to the k z axis (panel F).
semiconductor becomes progressively more isotropic as it turns into a semimetal, loosing its twodimensional character.

Exciton binding and instability
The exciton candidate for the instability has a finite center-of-mass momentum q, i.e., it travels in space. We compute its excitation energy-the difference between the GW bandgap and the binding energy-by solving the Bethe-Salpeter equation from first principles (Methods). The dispersion exhibits a dip for q = Λ, whose energy is first positive at P = 0 (1.26 eV, black dots in Fig. 3A) but then quickly lowers with P , eventually changing sign close to the semimetal threshold (−27 meV at P = 34 GPa, blue dots). This negative value signals that excitons spontaneously form, which leads to a reconstructed many-body phase of lower energy.
The softening of the exciton shown in Fig. 3A validates from first principles the seminal prediction by des Cloizeaux 2 : the binding energy remains finite even if the gap vanishes, as explicitly shown in Fig. 3B (black dots). The reason is that conduction and valence band profiles are almost unaffected by P (Fig. 2), as the band edges are displaced in k space, which prevents the macroscopic dielectric constant from diverging (red dots in Fig. 3B). Were the closing gap direct, metal-like screening would dissociate the exciton.
The square modulus of the exciton wave function is illustrated in Figs. 3C and D, as the conditional probability density to locate the bound electron (green contour map), provided the hole is fixed (black dot). Note that the center-of-mass motion does not appear in this frame.
The probability extends tens of Angstroms-the feature of Wannier excitons familiar from bulk The violet (yellow) colour in the stick-and-ball skeleton points to Mo (S) atoms.
semiconductors-both in-and out-of-plane, as apparent in panels C and D, respectively (the Bohr radius is 50Å at 34 GPa, as shown in SI Appendix, Fig. S5). The exciton becomes lighter and more isotropic with pressure, i.e., more delocalized in real space (here shown at P = 0).

Two-band model
The major source of numerical error is the finite sampling of the Brillouin zone 14 , since the exciton is significantly localized in k space while the computational load prevents us from refining the mesh

The excitonic insulator phase
Close to the semiconductor-semimetal boundary, the ground state undergoes a reconstruction from the 'normal' phase, |Φ 0 , which is either insulating or semimetallic, to the excitonic insulator, |Ψ EI . In the following, we highlight the essential features of |Ψ EI within the simpler two-band with the valence band. The new band structure at Γ is replicated at Λ, since both Γ and Λ points belong to the EI reciprocal lattice. Apart from spin degeneracy, EI renormalized bands exhibit an additional orbital degeneracy reminescent of the pristine multivalley structure: bands (solid curves) from top to bottom are respectively one-, three-, two-, and one-fold degenerate, respectively.
The two-fold degenerate band, which overlaps with the pristine conduction band (dashed curve), is actually split due to the tiny anisotropy of ∆ k in the kx, ky plane (splitting hardly visible in the plot). Circles point to EI conduction and valence bands obtained within the two-band model. model (as a mnemonic, we adopt the apex '0' to identify quantities of interest defined within this model). Then, we take into account the EI multivalley nature by adapting the theory first proposed for the candidate material TiSe 2 45 .
Within the two-band model 3 , |Ψ 0 EI is formally analogous to the superconductor wave func- provided the Cooper pairs of the metal are replaced with the e-h pair excitations of the normal state,b + kâ k |Φ 0 . Hereb + k creates an electron with momentum k + ΓΛ and energy ε b (k) in the conduction band,â k annihilates an electron with momentum k and energy ε a (k) in the valence The value of ϕ is-ideally-arbitrary and solely fixed by the spontaneous breaking of the conservation law for e-h pairs, as The EI band structure is obtained by solving the pseudo Bethe-Salpeter equation for ζ 0 k selfconsistently, where W (q) is the screened Coulomb interaction and the minimum value of 2E k is the bandgap (Methods). Reassuringly, Eq. 3 turns into the Bethe-Salpeter equation for the zero-energy exciton at the onset of the EI phase (∆ 0 k → 0+). As a consequence of the condensation energy gain, the EI conduction and valence bands (circles in Fig. 4A) are flattened and distorted with respect to those of the pristine semiconductor (dashed curves), the gap widening by ≈ 15 meV at P = 34 GPa.

Multivalley effects
As e-h pairs may be formed by exciting an electron from the valence band to any one of the six conduction band valleys, Λ i , the condensate wave function is multi-component 45 , ik creating an electron with momentum k + ΓΛ i and energy ε ib (k) in the ith valley (i = 1, . . . , 6). In principle, one must solve up to six coupled equations for ζ ik to account for the distortion of the condensate in k space, due to intervalley coupling. Nevertheless, we note that ∆ 0 k has hardly any angular dependence in the k x , k y plane (the maximum amplitude of the azimuthal modulation is smaller than 0.07 meV, see SI Appendix, Fig. S6), whereas ε ib (k) depends on the angle between ΓΛ i and (k x ,k y ) due to mass anisotropy. As Coulomb interaction protects the cylindrical symmetry of ζ ik , and since the bare-band anisotropy has negligible effect at valley bottom k ≈ 0 (where the value of ζ ik is largest), we neglect the azimuthal dependence of ζ ik and obtain (Methods): Here only the magnitude of ζ ik is fixed (from the self-consistent solution of equation 3), whereas the six phases ϕ i remain undetermined. This is sufficient to compute the band structure of the EI (Fig. 4A), as the ground state energy is independent from ϕ i (Methods).
There are now one valence and six conduction bands (solid thin lines in Fig. 4A), in place of the two bands (circles) of the superconductor-like model. Some of the conduction bands are degenerate, the degeneracy being respectively one, three, two, and one, from the topmost conduc-tion to the valence band. Importantly, the band structure at Γ is replicated at Λ, as the electronic charge exhibits a super-modulation in real space that we discuss below, the corresponding unit cell (solid frame in Fig. 4B) being larger than the cell of the crystal lattice (dashed frame). As a consequence, bands are folded into the smaller Brillouin zone (SI Appendix, Fig. S6), changing the gap from indirect to direct. Only the valence and topmost conduction bands repel each other, in agreement with the two-band model (circles), whereas the remaining bands, which are unaffected by the presence of the exciton condensate, replicate at Γ the bare valleys and hence reduce the direct gap. Since the location of the valence band top is slightly displaced from Γ along the k z axis (SI Appendix, Fig. S7), by ∼ 0.2 Bohr −1 , the actual EI gap is indirect and around ∼ 5 meV, smaller than the direct gap at Γ. Note that in Fig. 4A the two-fold degenerate band, which almost overlaps with the bare conduction band (dashed curve), splits due to the tiny anisotropy of ∆ k in the k x , k y plane (the splitting is hardly visible in the plot).

Anti-ferroelectric excitonic insulator
The EI ground state is invariant under time reversal, hence the phases of the condensate components that live in two antipodal valleys must have opposite sign (modulus a multiple integer of 2π), i.e., ϕ 1 = −ϕ 4 , ϕ 3 = −ϕ 6 , and ϕ 5 = −ϕ 2 (see SI Appendix, Fig. S6 and Methods). This constraint leads to the formation of a purely electronic, self-sustained charge density wave, ∆ (r), which breaks the inversion symmetry of the pristine crystal (the proof is given in the Methods). The total wave ∆ is the coherent superposition of three contributions, ∆ = ∆ 1,4 + ∆ 3,6 + ∆ 5,2 , each one originating from a couple of antipodal valleys. For example, exhibits the new periodicity 2π/| ΓΛ 1 | given by the momentum of those excitons that condense in valleys 1 and 4, and similarly ∆ 3,6 and ∆ 5,2 display an analogous modulation along directions ΓΛ 3 and ΓΛ 5 with phase shifts ϕ 3 and ϕ 5 , respectively. Here ψ Γ and ψ Λ 1 are the periodic envelopes of Bloch states respectively at Γ and Λ 1 , ψ Λ 4 = ψ * Λ 1 , and the spin has been factored out, since the lattice space group contains a center of inversion and a unique z axis 46 . It is clear that the total amount of charge displaced from the pristine background, as well as the amplitude of the charge modulation, are both driven by the condensate through k u 0 Importantly, the arbitrariness of the phases ϕ 1 , ϕ 3 , and ϕ 5 points to a huge, continuous degeneracy of the ground state. Since the effect of any given two arbitrary phases is merely to rigidly shift the charge pattern ∆ with respect to the frame origin (Methods), in the following we take ϕ 1 = ϕ 3 = ϕ 5 = 0. The resulting density wave is slightly distorted in the generic case, in which all three phases take arbitrary values (see discussion below). Figure 4C shows the overlap charge density of the envelopes obtained from first principles, This dipole may be regarded as the polarization of the excitons coherently built in the condensate 8 . Since the contributions to the dipole due to the remaining valleys, P 3,6 and P 5,2 , are obtained by rotating P 1,4 by respectively 2π/3 and −2π/3 along the z axis, the total dipole P = P 1,4 + P 3,6 + P 5,2 is parallel to the z axis. We evaluate this parallel component, P z , through direct integration over the unit cell ( Fig. 4D and Methods).
The overall charge pattern, P z (x, y), exhibits an anti-ferroelectric texture that breaks inversion symmetry. This is shown in Fig. 4B, where local dipoles, which point out of the plane, are depicted as red arrows having length proportional to |P z |. The electric dipole averages to zero over the unitary cell of the superstructure (solid frame), which contains 72 atoms [with ΓΛ ≈ 2π/(3a)] against 6 of the original cell (dashed frame). The reconstructed Brillouin zone, which is again hexagonal in the plane but rotated by π/6 (SI Appendix, Fig. S6D), is spanned by any two independent vectors chosen among the ΓΛ i 's. In the generic, degenerate case that ϕ 1 , ϕ 3 , and ϕ 5 take arbitrary values, we expect a reduction of the maximum local value of |P z | up to 2/3, together with a variable tilt of the dipole in the plane.

Semiconductor-semimetal crossover
The formation of a Fermi surface, made of six e pockets in the Λ valleys and one h pocket at Γ, signals the transition from the semiconductor (Fig. 5B) to the semimetal (Fig. 5D) occurring in the absence of excitonic effects. Figures 5B and 5D show one of the conduction valleys, displaced by − ΓΛ in k space, and the valence band, the filled states being shadowed by gray colour. As the free  5D) and the e-h excitation spectrum would be gapless.

Raman fingerprint
Were ion displacements responsible for the building of electric dipoles in place of excitons, the frequency of the phonon of momentum q = ΓΛ and consistent symmetry would soften (or at list exhibit a dip) at the onset of the new phase 49 . The phonon dispersion obtained from first principles, respectively at P = 0 (Fig. 6A) and 34 GPa (Fig. 6B), shows the opposite behaviour, with all lowenergy modes hardening with P (Methods and SI Appendix, Fig. S8 for the 2H a phase). Therefore, the anti-ferroelectricity has a purely electronic origin. This prediction is consistent with recent diffraction measurements, which ruled out any periodic lattice distortion above 40 Kelvin 43 .
The evolution of Raman spectrum with pressure, as obtained from first principles in Fig. 6C (structure 2H c ) and SI Appendix, Fig. S9 (2H a ), compares with observed data with the exception of the E peak at 174 cm −1 [ Fig. 4(b) of Ref. 35 ], which appears below 150 K and above 30 GPa but is missed by the theory for the normal phase. Cao and coworkers proposed 35  a transverse acoustic phonon of finite momentum, which becomes bright at the onset of a charge density wave, due to the reconstruction of the Brillouin zone. Whereas the first-principles spectrum for the excitonic phase is presently out of reach, below we confirm the essence of Cao's explanation by identifying E as the lowest optical phonon at Λ. This is the fingerprint of the anti-ferroelectric charge density wave associated with exciton condensation.
The symmetry group of the anti-ferroelectric ground state depicted in Fig. 4B only includes the identity operation. Therefore, all 216 vibrational modes are in principle infrared and/or Raman active. However, since the EI critical temperature is relatively low and the E peak is extremely bright, we expect that the new mode is an optical phonon of momentum Λ, which is Raman active through the folding into the zone center and strongly couples with P. Since P(x, y) originates everywhere in the cell from the inter-layer vertical displacement of the charge between two neighbour Mo atoms, it will mainly couple with those optical oscillations of Mo atoms that occur along the z axis. In fact, these vibrations linearly change the Mo-Mo distance and hence the local dipole strength, whereas the amount of displaced charge, which is ruled by the long-range part of Coulomb interaction, changes weakly with the oscillation. From direct inspection of phonon eigenvectors, there is one candidate only below 400 cm −1 , i.e., the lowest optical mode of frequency 164 cm −1 located at Λ, which is highlighted by a red dot in Fig. 6B. As shown by the displacement vectors in the EI reconstructed cell displayed in Fig. 6D, the Mo atoms oscillate out of phase along the z direction with an in-plane modulation of period 3a along the ΓΛ direction (parallel to the vertical axis of Fig. 4B), hence matching the periodicity of P z (x, y) in the plane.
This superlattice vibration is twice degenerate, due to the additional folding of the phonon with independent wave vector ΓΛ . Note that the observed intensity of the E mode is constant up to 60 K, which compares with the EI critical temperature. In summary, the E mode points to the excitonic insulator in the P − T space.

Discussion
Both 2H c and 2H a phases coexist [30][31][32]40 in the region of visibility of the E mode, which extends between 30 and 50 GPa at 5 Kelvin 35 . The lower bound agrees with our prediction, since in the 2H a structure the EI sets in at P ∼ 28 GPa (SI Appendix, Fig. S4) with a mode frequency of 166 Fig. S8). The upper bound of 50 GPa is larger than our expectation of ∼ 34 GPa for the 2H c phase. However, recent diffraction measurements on single crystals 43 , though only available at temperatures higher than 40 Kelvin, suggest that the critical upper pressure could be actually much lower, being artificially enhanced in powders due to the deviatoric stress field applied to randomly oriented crystallites.
In addition, other Raman features unexplained so far 35 point to the EI scenario: (i) the observation of modes supposedly forbidden or silent (ii) the anomalous frequency variation of the out-of-plane A 1g mode accompanying the onset of the E mode. Since the understanding of the available electrical transport measurements 30,35 is complicated by the mixture of phases in the high-pressure cell, we do not speculate on the origin of the resistivity peak that was tentatively assigned 30 to the EI.
The huge degeneracy of the EI ground state, associated with condensate phases ϕ 1 , ϕ 3 , and we estimate c exciton ∼ 2 · 10 4 m/s, which is much higher than the sound velocity of the stiffest acoustic phonon branch (Fig. 6B), c phonon ∼ 8 · 10 3 m/s. Therefore, the phase mode of the exciton condensate should be experimentally accessible.

Conclusion
In summary, we have demonstrated that a real excitonic insulator phase sets in between the semiconducting and semimetallic phases of MoS 2 , building on calculations from first principles and available spectroscopic data. These findings call for further investigation of some fascinating possibilities. A first question is the manifestation of the macroscopic quantum coherence of the exciton condensate, which might occur through the observation of low-lying collective modes associated with the oscillation of the condensate phase ϕ(r, t). Another issue is whether the superconductivity observed above 90 GPa is related to the excitonic phase, as the overscreening action of surviving exciton-plasmons might act as unconventional glue for Cooper pairs. We hope our study may stimulate further work along these paths.

Methods
Computational details of ground-state calculation from first principles The lattice parameters and the ground-state electronic structure for the three values of pressure were obtained within density functional theory (DFT), with a plane wave basis set as implemented in the Quantum ESPRESSO package 51, 52 , using the generalized gradient approximation Perdew-Burke-Ernzerhof (PBE) parametrization 53 . A kinetic energy cutoff of 100 Ry was adopted for the wave functions, and fully relativistic norm-conserving pseudopotentials 54 were used to take into account spin-orbit interaction. Van der Waals interactions, included by using the Grimme approximation method, were found to be relevant only at zero pressure, as already shown in Ref. 29 .
Phonons Phonon dispersions were calculated by using a Density Functional Perturbation Theory approach 55 . We used a 10 × 10 × 3 Monkhorst-Pack grid for the integration in the Brillouin zone; the dynamical matrix at a given point of the Brillouin zone was obtained from a Fourier interpolation of the dynamical matrices computed on a 5 × 5 × 1 q-point mesh.
Quasiparticles and excitons Many-body calculations 33,56,57 were performed by using the Yambo code 58,59 . Quasiparticle corrections to the Kohn-Sham energies were evaluated using the G 0 W 0 approximation for the self-energy, the dynamical dielectric screening been accounted for within the plasmon-pole approximation 60 . To speed-up the convergence of quasiparticle energies with respect to the number of empty bands in the sum over states occurring in the calculation of the polarizability and self energy, we have adopted the scheme proposed in Ref. 61 Fig. S12).
Computational details of the two-band model The effective-mass framework builds on the knowledge of conduction and valence energy bands. Here G > 0 (G < 0) is the indirect bandgap (band overlap) for pressures below (above) the semiconductor-semimetal threshold-in the absence of excitonic effects-and the momentum components, k , k ⊥ , k z , are projected along the principal axes of the effective mass tensor 64 , the corresponding masses being m i , m i⊥ , m iz , with i = a, b. These axes are respectively parallel (k ) and perpendicular [in-(k ⊥ ) and out-of-plane (k z )] to the ΓΛ direction, the axis origin being placed at the band edge. We emphasize that all parameters of the two-band model, for a given pressure, are fixed and obtained from first principles. In particular, the bandgap and the effective masses are extracted from GW bands, as illustrated in Figs. 2D to F, and hence include the mean-field renormalization due to e-e interactions. The (modulus) of the screened e-h Coulomb attraction in momentum space, depends on the static dielectric constant, κ r , which is obtained as the inverse of the first-principles dielectric tensor, 1/[ −1 (q = 0)] G=G =0 , in the long-wavelength, macroscopic limit, as illustrated in Fig. 3B (here Ω is the crystal volume and G the reciprocal lattice vector).
In the semimetal, the P -dependent values of G, m i , m i⊥ , m iz , and κ r are derived as linear extrapolations of first-principles data at P = 25 and 34 GPa, respectively. Since free e and h carriers effectively screen the interaction by adding a metal-like, intraband contribution to the polarizability, we modify the dressed Coulomb potential as Here the Thomas-Fermi term, proportional to the density of states, D(ε), evaluated at the Fermi energy, ε F , removes the long-wavelength divergence of W . We obtain numerically D through the summation of localized Gaussian functions over a fine grid in k space, as well as ε F by imposing overall charge neutrality (we take into account the six-fold degeneracy of conduction band).

Two-band Bethe-Salpeter equation
In the semiconductor, the exciton wave function is where φ k is the probability amplitude of a bound e-h pair in momentum space. The Bethe-Salpeter equation of motion for φ k is where ε exc is the excitation energy of the exciton, whose negative value signals the instability. We solve this equation by numerical discretization in k space and assess convergence by refining the mesh as well as varying the momentum cutoff. Note that the singularity of Coulomb potential for |q| → 0 is harmless, as we integrate W over a small parallelepiped, in a semi-analytical, accurate manner. We have benchmarked the convergence of our calculations against known analytical or high-precision results, as shown for bulk Wannier excitons in SI Appendix, Fig. S13 and for anisotropic excitons with a well-defined azimuthal quantum number 65 in SI Appendix, Fig. S14.
In the semimetal ground state, a small area of k space around the origin is populated by electrons in band b and holes in band a. In addition, due to band anisotropy 66 , in narrow regions nearby there are either electrons or holes only, which prevents from exciting e-h pairs due to Pauli exclusion principle. Therefore, the Bethe-Salpeter equation of motion must be modified as 4 where n i (k) is the occupancy factor of the ith band in the normal ground state, which takes value either 0 or 1. The 'counting' prefactor of W , [n a − n b ], removes scattering channels forbidden by Pauli blocking and is responsible of the plasmon-like features shown Fig. 5E. Note that in the semiconductor, n a (k) = 1 and n b (k) = 0, hence one regains the standard form of equation 11.
Self-consistent theory of the excitonic insulator within the two-band model The EI bands (circles in Fig. 4A) Coulomb interaction W is renormalized by a vertex correction associated with the EI ground state 47 , as the opening of the many-body gap significantly enhances the e-h attraction-by suppressing screening-with respect to the gapless normal phase. Therefore, following Kozlov and Maksimov 47 , for small momentum transfer q the dressed interaction W appearing in equation 3 of main text takes the self-consistent form where the gap function at the Fermi surface, ∆ 0 k F , which is determined recursively, removes the long-wavelength divergence as one approaches the EI-semimetal boundary. Here ∆ 0 k F is an average value defined as ∆ 0 , with k xF given implicitly by ε F = ε b (k xF , 0, 0), and similarly for k yF and k zF . The constant α, for given band overlap G < 0, is α = [|G 0 | (ε F − G/2) 3/2 ] 1/2 , where |G 0 | = 9.38 meV is the maximum magnitude of the band overlap at which e-h pairing takes place. We neglect the modification of Eq. 13 for large momentum transfer, as it turns out to be irrelevant numerically. Whereas the vertex form 13 was originally proposed 47 for the case of spherically symmetric e and h pockets, we notice that, at the semiconductor-semimetal threshold, the exciton responsible for the instability is essentially isotropic (SI Appendix, Fig. S5).
At finite temperature, T , the gap equation takes the form where f F (x) = 1/[1 + exp (βx)] is Fermi distribution function, with β = 1/k B T and k B being Boltzmann constant, and we neglect the small renormalization of the chemical potential due to the presence of the exciton condensate.
Multivalley band structure The calculation of the EI band structure relies on the theory by Monney and coworkers 45 to include valley degeneracy. This approach, based on Green functions, generalizes to multiple bands the original theory by Jérome and coworkers 3 . For every k point, the EI band energies (solid lines in Fig. 4A and SI Appendix, S7) are found as the seven roots of the equation [cf. Eq. (8) of Ref. 45 ], after the magnitudes of the excitonic gap components, ∆ i (k), are obtained as follows. The gap function is defined as with ζ ik , apart from a phase factor, being the equal-time interband excitonic coherence F † i (k, t, t) defined in Eq. (4) of Ref. 45 , and δ → 0 + being a positive infinitesimal quantity. The integral 17 is evaluated through contour integration, the Fourier transform F † i (k, ω) being derived from the equations of motion of Green functions 45 as Whereas this expression would generically lead to an intractable system of six coupled equations for the ∆ i 's, we exploit the high symmetry of the problem to simplify the form of F † i (k, ω) and recover a single gap equation. As discussed in the main text and SI Appendix, Fig. S6, the symmetrizing effect of e − h attraction makes ∆ 0 k almost independent from the azimuthal angle ϕ k , with k ≡ (k, ϕ k , k z ) being expressed in cylindrical coordinates (k is the in-plane radial distance and k z the component along the z axis). Therefore, it is natural to assume that ∆ i has cylindrical symmetry, ∆ i (k) = ∆(k, k z )e iϕ i . Since we are mainly interested in the region k ≈ 0, we also neglect the azimuthal dependence of ε ib (k) in the denominator of F † i , obtaining where we omitted the dependence of terms on k in the notation. Equation 19 is now easily integrated, giving a single self-consistent gap equation. This has the same form of the equation 3 of the two-band model, provided that ∆ 0 k is replaced with Ground state wave function The contour integration of equal-time Green functions provides us with all interband coherences and band populations, i.e., Ψ EI |b + jkb ik |Ψ EI = ∆ * i ∆ j /2E(E+ε b /2− ε a /2), Ψ EI |b + ikb ik |Ψ EI = (v 0 ) 2 /6, Ψ EI |â + kâ k |Ψ EI = (u 0 ) 2 , where we omitted the dependence of right-hand-side terms on k to ease the notation, neglected the in-plane anisotropy of valence band, ε ib = ε b , and put E = {[ε b − ε a ] 2 /4 + |∆ 0 | 2 } 1/2 . This allows us to write explicitly the ground state wave function, in terms of Bogoliubov-Valatin-like creation operators,γ + , which are defined aŝ As discussed in the main text, time reversal symmetry limits the number of independent condensate phases to three: ϕ 1 , ϕ 3 , and ϕ 5 (recall that u 0 k = u 0 −k and v 0 k = v 0 −k are real positive quantities; see SI Appendix, Fig. S6C).
Inversion symmetry breaking The ground state wave function allows us to understand the symmetry breaking associated with exciton condensation. The inversion operator,Î, acts differently on b i and a Bloch states, since the envelope function at Γ is odd:Îâ + k = −â + −k ,Îb + 1k =b + 4−k , etc.
The magnitude of the expression enclosed in curly brackets is less than one (unless ϕ 1 = ϕ 3 = ϕ 5 = ±π/2, i.e.,Î |Ψ EI = − |Ψ EI ), hence, in the thermodynamic limit, the overlap between I |Ψ EI and |Ψ EI tends to zero as the two states become orthogonal. Since the ground state has a lower symmetry than the Hamiltonian, inversion symmetry is broken.
Charge density wave The form 5 of the purely electronic charge density wave, ∆ = ∆ 1,4 + ∆ 3,6 + ∆ 5,2 , is derived in a straightforward manner by averaging the density operator,ˆ (r) = 2ψ † (r)ψ(r), over |Ψ EI , with the Fermi field operator,ψ(r), being defined aŝ Cross-terms proportional to ψ * Λ i ψ Λ j average out to zero, once summed together, as the various ψ Λ i 's are obtained one from the other by either rotation by ±2π/3 along the z axis or complex conjugation. Apart from the envelope functions, which have the lattice periodicity, ∆ depends on r through a sum over three exponentials, whose imaginary arguments are respectively (times the prefactor i) ΓΛ 1 · r − ϕ 1 , ΓΛ 3 ·r − ϕ 3 , and ΓΛ 5 ·r − ϕ 5 , as illustrated in the main text.
Let us construct explicitly R shift as R shift = −R − R ⊥ , where R = n t 2 and R ⊥ = n ⊥ (2t 1 + t 2 ) are respectively parallel and perpendicular to ΓΛ 1 (SI Appendix, Fig. S6C), n and n ⊥ are integers to be determined, and t 1 , t 2 are the primitive vectors that generate the hexagonal lattice in Mattheiss' coordinate frame 46 . Since ΓΛ 1 is generically not commensurable with the reciprocal lattice vectors, there exists an integer n such that ΓΛ 1 ·R = ϕ 1 with arbitrary accuracy 4 , modulus an integer multiple of 2π. Similarly, we may fix n ⊥ such that ΓΛ 3 ·R ⊥ = − ΓΛ 5 ·R ⊥ = ϕ 3 − ΓΛ 3 ·R .
This theorem implies that the set of charge density waves [∆ (r)] 0,0,ϕ 5 labeled by the continuous parameter ϕ 5 spans all possible modulations of the electronic charge density of the EI, each realization having in turn a huge translational degeneracy, which is parametrized by the two continuous variables ϕ 1 and ϕ 3 .
Anti-ferroelectric order The electronic charge density wave of the EI ground state (Eq. 5 of main text) induces an out-of-plane electric dipole, P z (R i ), in the ith cell of the pristine 2H phase located at R i , with i = 1, . . . , N (N is the total number of cells). This is illustrated in Fig. 4B, where the dipoles P z (R i ) are depicted as red arrows. The local dipole P z (R i ) is given by the coherent superposition of three density waves, whose characteristic wave vectors are q i = ΓΛ i , with i = 1, 3, 5, 3,5 cos (q j · R i ).
The maximum value, P z (0), is shown in Fig. 4D. Here P z (R i ) is evaluated within the envelope function approximation, the factor P z0 being derived from first principles through the overlap charge density of the periodic part of conduction and valence Bloch states at Γ and Λ, respectively, which is shown in Fig. 4C. The latter is numerically integrated over the pristine unit cell volume, Ω cell : the frame origin being placed at the inversion center-the midpoint between the two Mo atoms of the 2H cell. As the charge displacement that gives rise to the dipole is essentially localized on Mo atoms (Fig. 4C), we expect |P z0 | to be well defined. We obtain P z0 /e = 15.1 Bohr at P = 34 GPa.