Evidence for equilibrium exciton condensation in monolayer WTe2

We present evidence that the two-dimensional bulk of monolayer WTe2 contains electrons and holes bound by Coulomb attraction—excitons—that spontaneously form in thermal equilibrium. On cooling from room temperature to 100 K, the conductivity develops a V-shaped dependence on electrostatic doping, while the chemical potential develops a step at the neutral point. These features are much sharper than is possible in an independent-electron picture, but they can be accounted for if electrons and holes interact strongly and are paired in equilibrium. Our calculations from first principles show that the exciton binding energy is larger than 100 meV and the radius as small as 4 nm, explaining their formation at high temperature and doping levels. Below 100 K, more strongly insulating behaviour is seen, suggesting that a charge-ordered state forms. The observed absence of charge density waves in this state is surprising within an excitonic insulator picture, but we show that it can be explained by the symmetries of the exciton wavefunction. Therefore, in addition to being a topological insulator, monolayer WTe2 exhibits strong correlations over a wide temperature range. Exciton condensation has been observed in various three-dimensional (3D) materials. Now, monolayer WTe2—a 2D topological insulator—also shows the phenomenon. Strong electronic interactions allow the excitons to form and condense at high temperature.

I n a monolayer semimetal, low carrier densities combined with reduced dimensionality provide the conditions for strong correlation effects. One possible form of such correlations is pairing of electrons and holes in the equilibrium state to form excitons. At low temperatures, such excitons could condense to form an excitonic insulator [1][2][3][4] . However, exciton formation is expected to be easily disrupted by free carriers (which screen the binding interaction) and thus occur only at low temperatures and near charge neutrality. A number of materials have been mooted as excitonic insulator candidates [5][6][7][8][9][10] , but there is no consensus as to whether any of them truly contains excitons in equilibrium, either as an incoherent gas or as a coherent condensate, except in the case of bilayer heterostructures at high magnetic fields [11][12][13][14] .
WTe 2 is a layered semimetal, but an exfoliated WTe 2 monolayer behaves 15-17 as a two-dimensional (2D) topological insulator, exhibiting helical conducting edge modes, and becomes superconducting when electrostatically doped 18,19 . The monolayer has the 1T′ structure shown in Fig. 1a. Its bands are spin degenerate due to inversion symmetry and, near the Fermi energy E F , there is a valence (v) band maximum at Γ flanked by two conduction (c) band minima located at k x = ±k Λ , as sketched in Fig. 1b. Some tunnelling spectroscopy measurements 20 , angle-resolved photoemission 20,21 and density functional theory (DFT) calculations 20-23 point to a positive bandgap, E g , of the order of 50 meV, while others suggest overlapping bands 15,24 . Noting that, in the photoemission spectra, the v and c band photoemission features are broad enough that they overlap at least somewhat, we use thick lines in the sketch to signify this uncertainty in E g .
This band structure immediately invokes the possibility that excitons could occur in equilibrium in monolayer WTe 2 , and therefore that the insulating state might not be a simple band insulator [25][26][27][28] . In this Article we argue that the behaviour of the conductivity and the electron chemical potential, even well above 100 K, is impossible to reconcile with an independent particle picture and strongly indicates the presence of excitons in the equilibrium state. Our first-principles calculations of exciton dispersion and Bohr radius support this conclusion. The insulating behaviour below 100 K suggests a charge-ordered state, but in an excitonic insulator one would normally expect charge density waves, no signs of which are seen in scanning tunnelling microscopy or Raman spectroscopy. To explain this, we show that the entanglement of spin, orbital and valley degrees of freedom hides the charge order, as the contributions to the density wave paired through time reversal cancel out.
The measurements were made on exfoliated monolayer WTe 2 flakes with platinum contacts, encapsulated by hexagonal boron nitride (hBN), with graphite gates either below or above (Supplementary Section 1). To study the sheet conductivity while excluding edge conduction, we used the approach illustrated in Fig.  1c. As indicated in the insets, a bias V is applied to one contact and the current I flowing to ground through an opposite contact is measured. When the intervening side contacts are grounded, this current must flow through the bulk. The edge current, which produces the plateau at low V g , is thereby eliminated and the 'partial conductance' G p ≡ I/V reflects the sheet conductivity σ via G −1 p ≈ β/σ + Rc, where β is a geometrical factor considerably larger than one and R c is the contact resistance. The gate-induced areal number density n g is deduced from the voltages applied to the graphite gate(s) and the geometric capacitances. Figure 1d shows measurements of G p versus n g and temperature T. These characteristics are not measurably affected by either a normal displacement field or a normal magnetic field of 14 T, confirming the rejection of edge conduction, which is highly sensitive to magnetic field at low temperatures 29 . (A recent paper 30 reports surprising quantum oscillations at low n g , but we have not seen these in any device.) On cooling from room temperature to 100 K, G p versus n p develops a sharp 'V' shape centred close to n g = 0. We have seen consistent behaviour across a dozen monolayer devices, although the sharpness of the V varies, probably as a result of variable sample homogeneity. We also note that a similar sharp V occurs in bilayer WTe 2 but at temperatures about five times lower 31 (Supplementary Section 7). As shown in the inset, for positive n g smaller than a value n ce ≈ +5 × 10 12 cm −2 , G p decreases monotonically on cooling, whereas for n g > n ce it initially increases. For negative n g (hole doping), a similar but less clearcut transition occurs around n cp ≈ −10 × 10 12 cm −2 . As indicated in the band at the bottom of Fig. 1d, these values of n ce and n cp are consistent with the thresholds for metallic behaviour reported in previous work, where it was found that the metallic state at n g > n ce becomes superconducting below ~0.8 K.
Below 100 K, as T decreases G p collapses over an increasingly wide range of n g . This insulating behaviour, which allows the edge conduction to dominate in normal geometries, is consistent with the findings in ref. 26 . We used microwave impedance microscopy 32 (MIM; Supplementary Section 2) on devices with no top gate to confirm that this is not a contact effect, as well as to detect any cracks in  a, 1T′ structure of monolayer WTe 2 . The x axis is taken to be along the zigzag W chains. b, Schematic Brillouin zone (above) and bands near E F (below). c, The technique used to exclude edge conduction: when the side contacts are grounded, the measured current between top and bottom contacts must flow through the interior and the edge conduction plateau disappears. (G p ≡ I/V). These measurements were taken on device MW2 at T = 10 K. d, G p versus gate-induced density n g at a series of temperatures T on the same device. Inset: temperature dependence at positive values of n g . Bottom: regimes of insulating, metallic and superconducting behaviour identified in previous work.  Here, the uncalibrated imaginary MIM signal is plotted for the purpose of identifying conducting features. Multiple conducting cracks are seen running from upper left to lower right. Inset: polar plot of the Raman peak P11 (210 cm −1 ) intensity, used to determine the crystal axes. b, Conductivity versus n g and T at 2 GHz deduced from MIM measurements in the centre of the white dotted square in a. c, Partial conductance parallel to the x axis (left) and y axis (right), measured in the white dotted square using the configurations shown in the insets. the monolayer WTe 2 that could invalidate the measurements. Figure  2a is a MIM image of device MW10 at 11 K. Red dashed lines mark the edges of the monolayer WTe 2 flake, and the wiggly bright lines are cracks. Figure 2b shows the MIM-derived conductivity σ MIM , measured in the centre of the white dashed square. Like G p , it collapses over a range of n g that grows as T falls, indicated by the dotted white contour which is drawn at σ MIM ≈ 0.1 μS. The drop-off of G p at low temperatures, even at large n g , as is evident in Fig. 1d, can be explained by the fact that the monolayer adjacent to the metal contacts is partially screened from the gates and so is less doped and remains insulating. In addition, the contact pattern in MW10 was aligned with the crystal axes, as determined by Raman spectroscopy (Fig. 2a inset and Supplementary Section 3), allowing us to compare conductivities along the x and y axes (Fig. 2c,d and Supplementary Section 4). We see that there is substantial gate-dependent anisotropy, with the lowest conductivity occurring for p-doping parallel to the x axis. This is consistent with the direction in which the valence band edge has a large effective mass.
To measure the chemical potential, we study a device including a separately contacted graphene sheet in parallel with the WTe 2 monolayer, as shown in Fig. 3a. Briefly, the WTe 2 is almost an equipotential because it has finite conductivity and carries no current. With both the graphene and bottom gate grounded, a voltage V g is applied to the top gate relative to the WTe 2 and the voltage V W on the WTe 2 is adjusted to bring the graphene conductance to a minimum. This keeps the graphene neutral and maintains zero electric field beneath the WTe 2 . The electrostatic potential in the WTe 2 is thus in effect fixed to that of the graphene, so the change in V W is due to the change in chemical potential, Δμ = −eΔV W , associated with the gate-induced charge density −en g = ε r ε 0 V g /d, where d is the distance to the gate ε r is the hBN dielectric constant and ε 0 the vacuum permitivity. From V W versus V g we thereby obtain μ(n g ), choosing the zero of μ at each temperature for convenience. Figure 3b shows measurements of both μ (black) and G p (red) versus n g made on device MW12. As usual, G p forms a sharp V as a function of n g at 100 K. Meanwhile, μ exhibits a step at the centre of the V that is still discernible at room temperature and which grows on cooling, saturating at ~40 meV in height below ~50 K. The same behaviour was seen in two devices (Supplementary Section 5).
The variations of G p and μ with n g are impossible to reconcile with a single-particle picture, as can be shown by an elementary calculation: in such a picture, μ and n g are related by is the total electron density of states and f(E) = [1 + exp{(E − μ)/kT}] −1 . To match the variation of μ with n g at low temperatures, for energies in the relevant range, D(E) must have roughly the form shown in Fig. 3c: a constant value, D 0 , in the conduction band to give a uniform slope for n g > 0; a gap, E g ; and a large spike at the valence band edge (E = 0). The latter is needed to pin μ to the valence band edge and thus make it independent of n g when n g < 0, to be consistent with the data at 25 K

Fig. 3 |
Chemical potential measurements and comparison with the single-particle model. a, Schematic of the device structure used to measure the chemical potential versus doping. In parallel with (below) the monolayer WTe 2 is a graphene sheet, which is maintained at its Dirac point so that the electric field beneath the WTe 2 is zero. The doping n g and chemical potential μ are then obtained from the voltages as shown (see text). The hBN dielectric between the layers is not shown. b, Measurements of μ (black) and conductance G p (red) versus n g on monolayer WTe 2 device MW12. The length of the green bars indicates the thermal energy, kT. c, Single-particle density of states D(E) that reproduces the low-temperature μ-n g behaviour. d, Calculated μ (black line), electron density n (red dotted line) and hole density p (blue dashed line) using this D(E) at the same temperatures as the measurements. and 10 K. In Fig. 3d we plot the chemical potential calculated at the same temperatures as the measurements in Fig. 3b, using this D(E) with best-fit parameters D 0 = 3.7 × 10 11 cm −2 meV −1 and E g = 43 meV. At 150 K the step is washed out; in fact, no choice of D(E) can yield a distinct step in μ whose height is less than kT as is needed to match the measurements at T ≥ 150 K. In particular, states in the gap will only smear the step more. We also plot (dotted lines) the calculated populations of electrons in the conduction band, n, and holes in the valence band, p = n g − n, to contrast their thermally smeared dependence on n g with the sharp V shape seen in the conductance when T ≥ 100 K. The contradictions between the single-particle picture and the observed dependence of μ and G on n g and T can be largely resolved simply by positing that some electrons and holes are bound as neutral excitons with density n x so that the conductivity is σ = μ e e(n − n x ) + μ h e(p − n x ), the sum of the contributions of n − n x free electrons and p − n x free holes with respective mobilities μ e and μ h . In Fig. 4 we compare predictions based on this equation with data from Fig. 3 (Supplementary Section 6 provides more details). When n x = 0 we just have σ = μ e en + μ h ep, which is thermally smeared for any choice of D(E). To illustrate this, in Fig. 4a we plot p and n calculated using the same D(E) shown in Fig. 3c, at 100 K, and in Fig. 4b we plot the calculated σ for n x = 0 (blue dotted), using a mobility ratio chosen to obtain the best match to the measured conductance at 100 K (black line). However, when n x takes its maximal value, determined by the number of minority carriers n x = min(n, p), then for n g > 0 we have σ = μ e e(n − p) = μ e en g while for n g < 0 we have σ = μ h e(p − n) = μ h en g . The result is a sharp, asymmetric V shape (red dashed line) that matches the measurements much better. Note that, in this limit, where only the unbalanced gate-induced charge is free to move, the detailed form of D(E) becomes immaterial.
As T is increased from 100 K, the conductance at n g = 0 rises and the sides of the V become shallower. This is illustrated in Fig. 4c, where we replot the conductance at 150 K (black line). The behaviour remains highly incongruous with the single-particle model (n x = 0, blue dotted line), and is again more similar to the calculation for the case of maximal n x with slightly decreased mobilities (red dashed line). However, there is now a discrepancy in the form of a vertical shift, equivalent to an extra gate-independent contribution to the conductance. The vertical shift cannot be accounted for by varying n x ; to illustrate this we also plot (green dash-dotted line) the result of assuming that n x varies with gate voltage according to a chemical equilibrium condition, n x = K(n − n x )(p − n x ), where K is an equilibrium constant. It also cannot be reproduced by using a single-electron spectrum, for example with overlapping bands. A more sophisticated treatment of the correlated-electron system may therefore be needed to understand this aspect of the behaviour.
Excitons that persist in equilibrium at 100 K and at doping levels above 1 × 10 12 cm −2 must have binding energy much larger than the thermal energy of ~10 meV and small size to survive screening by free charges. To see whether this is plausible, we solved the exciton (Bethe-Salpeter) equation of motion from first principles, building on the DFT band structure (Fig. 5a) and including spinorbit effects in a non-perturbative way (Methods). The resulting excitation energy versus momentum q is shown in Fig. 5b. As the binding energy only weakly depends on E g because the gap is indirect 2,10 , and the value of E g is uncertain, we tuned the DFT hybrid functional to make E g vanish. The dielectric function was evaluated in the random phase approximation, and the uncertainty induced by the numerical discretization of k space is shown by the error bars. The excitation energy is negative for all q, ranging from −100 meV for direct excitons at q = 0 to a minimum of −330 meV for indirect excitons made of a hole at Γ and an electron at Λ. In Fig. 5c,d we plot the spatial profile of an exciton in the centre-of-mass frame. The exciton radius is as small as 4 nm. This is comparable with the typical electron separation at the critical doping nce −1/2 = 4.5 nm, suggesting that excitons could play a role in the insulator-metal transition at n ce .
In interpreting the situation below 100 K in terms of an excitonic insulator, formed by condensation of the excitons, we face two problems. First, at low temperatures, insulating behaviour sets in over a wide doping range, approaching n cp < n g < n ce . Conventional neutral excitonic insulator theory provides no mechanism for localizing the unbalanced charge. Nevertheless, it seems possible that the Coulomb interaction, which has a long range for these doping values, could stabilize both the excitonic phase 10 and the Wigner crystallization of unbound carriers, whose effective mass is enhanced by the opening of a many-body gap. Second, the exciton condensate is naively expected to exhibit a charge density wave (CDW) with wavevector k Λ , but CDWs are not seen in tunnelling microscopy 20,24,33 and our detailed temperature-dependent Raman spectra show no evidence of any CDW transition (Supplementary Section 3). One possible explanation is that the condensate is made of direct excitons having q = 0, as predicted 25 for the T′ phase of monolayer MoS 2 . We have checked that this leads to no substantial symmetry breaking, due to the anisotropic character of the WTe 2 band structure (Methods and Supplementary Fig. 8). However, this state would have higher energy than a condensate made of indirect excitons with finite q. A more likely possibility is that the peculiar symmetries of excitons with q = ±k Λ prevent the condensate from exhibiting charge order. A mean-field theory of this condensate (Methods), which builds on the knowledge of indirect excitons obtained from first principles and takes into account both spin and valley degrees of freedom, shows that there is no CDW with momentum q = k Λ . We find that the charge order is hidden by the combined effect of time-reversal symmetry and spin-orbit interaction, which is ultimately related to the topological properties of WTe 2 (ref. 15 ). The CDW may be unveiled by breaking time-reversal symmetry; practically, one may split the hole states that sustain the density wave through Zeeman coupling with the magnetic field, as the hole spin is polarized along x close to Γ. Finally, the ground state exhibits a spin density wave of momentum q = k Λ , plus a weak charge modulation of period q = 2k Λ . These features were also found in ref. 28 by adding an intervalley scattering term to a two-band model, in the absence of which different ordered states would have been degenerate. In contrast, the hidden order we predict here is inherent to all possible spin symmetries of the condensate, and hence robust.
We also consider the behaviour of μ versus n g . We start by putting forward the following heuristic argument. At low temperatures (T ≲ 50 K), where every minority carrier is paired, for n g < 0 when an electron is added to the system it pairs with a hole, reducing the addition energy by the exciton binding energy and leading to diverging compressibility and self-consistent pinning of μ (in the charge-ordered state). For n g > 0, the added electron does not pair and so μ is not affected by the binding energy; the result is a step in μ at n g = 0 whose height is related to the binding energy. The decrease in the height of the step for T ≳ 100 K could be because the exciton binding weakens due to screening by the free carriers. This scenario is supported by simulating the behaviour of an excitonic insulator in the presence of free charge carriers (Fig. 5e), with both μ and the many-body gap being computed self-consistently (Methods). Here the gap has a purely excitonic origin, as the starting non-interacting phase is a semimetal. The step in μ remains clearly visible up to 100 K, in contrast with the smeared profile of the independent-electron model (Fig. 3d). This is a peculiar consequence of exciton condensation, as electrons injected for n g > 0 fill in the lowest c states blockading the formation of e-h pairs-an effect due to the Pauli exclusion principle and suppressed with temperature 34 . At even higher T the step is smeared anyway, because the condensate is depopulated by thermal excitations and the excitonic gap melts. The simulation of Fig. 5e was performed for direct excitons for the sake of illustration 25 , but indirect excitons will exhibit the same qualitative features.
In conclusion, examination of the conductivity and thermodynamics of monolayer WTe 2 provides evidence that neutral excitons are present in thermal equilibrium, not only at low temperatures where they probably form a collective excitonic insulator state, but also even at temperatures above 100 K when they coexist with free carriers. The size of the step in chemical potential and theoretical calculations suggest that the exciton binding energy is at least 40 meV.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41567-021-01427-5.

Methods
Ground-state calculations from first principles. We obtained the ground-state electronic structure within DFT with a plane wave basis set, as implemented in the Quantum ESPRESSO package 35 . We fixed a kinetic energy cutoff of 80 Ry for the wavefunctions and used fully relativistic norm-conserving pseudopotentials 36 to include the spin-orbit interaction. We optimized the lattice parameters and atomic positions using the PBE exchange-correlation functional, the final cell parameters being a = 3.52 Å and b = 6.29 Å. We set the cell side along z to 15 Å. We obtained the band structure using a PBE0 pseudopotential, for which we considered a small fraction of exact exchange, 2%.
Excitation energies and exciton wavefunctions from first principles. We calculated the excitation energies of excitons as well as the dispersion of the lowest-energy exciton within the framework of many-body perturbation theory [37][38][39] , by solving the Bethe-Salpeter equation through the YAMBO code 40,41 and including spin-orbit interaction in a non-perturbative way 42 . We considered the PBE0 electronic structure as a starting point and calculated the static screening in the direct term within the random phase approximation, with inclusion of local field effects, and we employed the Tamm-Dancoff approximation for the Bethe-Salpeter Hamiltonian. To avoid spurious interactions among layers, we employed a truncated Coulomb cutoff technique 43 . We obtained converged excitation energies considering, respectively, two valence and two conduction bands in the Bethe-Salpeter matrix, the irreducible Brillouin zone being sampled up to a 48 × 24 × 1 k-point grid. We extrapolated the excitation energy of the lowest exciton with momentum q = 0 to the limit of a dense k-point grid, as shown in Supplementary Fig. 7.

Excitonic insulator ground state from indirect excitons and synopsis of theory.
From the first-principles solution of Bethe-Salpeter equation we find that the lowest-energy exciton with q = k Λ (or −k Λ ) is three-fold degenerate and separated by 20 meV from a non-degenerate first excited state. This energy splitting is due to the residual exchange interaction, present despite the strong spin-orbit coupling. These excitons are made of electrons and holes that populate, respectively, the lowest c and highest v band, with momentum k lying near the ΓΛ line. Along this high-symmetry line, both c and v Bloch spinors, each doubly degenerate, may be chosen as the eigenstates of the two-fold screw-axis rotation, which has complex conjugated values as irreducible representations. We label these Bloch spinors as ζ = ±i and use them to write explicitly the exciton wavefunctions within a minimal two-band model, which includes both spin and orbital degrees of freedom (Spinful two-band model). In addition to the rotational symmetry, we may classify the excitons according to their triplet-or singlet-like character, that is, whether they respectively maximize or minimize the electron-hole spatial overlap. In fact, Bloch spinors labelled by ζ transform like spins polarized along the x axis under the two-fold rotation, even if their actual spin polarization away from Γ is zero. Specifically, the exciton wavefunctions that are even with respect to the screw-axis rotation are where +/− stands for singlet/triplet-like symmetry, c † ζ (k) creates and v ζ (k) destroys an electron of momentum k and screw-axis symmetry ζ in the c and v band, respectively, |0〉 is the non-interacting ground state with all v states filled and c states empty, and ψ is the exciton wavefunction in reciprocal space in the centre-of-mass frame, which is nodeless and even in k y (and, approximately, in k x ): ψ(k x , k y ) = ψ(k x , −k y ).
In the excitonic insulator phase, the condensate is macroscopically occupied by excitons, hence the expectation value of the operator that creates an electronhole pair, ⟨ c + ζ v ζ′ ⟩ , has finite magnitude ϕ ζζ′ and arbitrary phase θ ζζ′ ; that is, the complex wavefunction of the condensate is ϕ ζζ′ exp(iθ ζζ′ ) (here ⟨…⟩ is the average over the many-body ground state). As excitons with q = (k Λ , 0) and (−k Λ , 0) are degenerate, they may separately condense, so the wavefunction of the condensate has two valley components of equal magnitude, whose respective phases are related by time-reversal symmetry 10 . Furthermore, both condensate and exciton wavefunctions share the same screw-axis and spin-like symmetries 3 , as well as the same parity in k space. These fundamental constraints relate all components of the condensate to a unique wavefunction magnitude, ϕ(k), and phase, θ. If the excitonic ground state has triplet-like symmetry, one has For the singlet-like ground state, one replaces the 2 × 2 Pauli matrix σ z with the identity matrix 1, the right-hand side of the equation then reading 1exp(±iθ)ϕ + (k). The charge/spin density wave of the excitonic insulator is dictated by the interband contribution to the expectation value of the corresponding density operator, which is proportional to the left-hand side of the above equation. Therefore, we may assess the occurrence of charge order in the excitonic phase without actually computing ϕ ± (k), which is given by a gap equation 10 that depends on E g . Indeed, ϕ only provides the intensity of the density modulation, whereas θ rigidly shifts the density wave with respect to the frame origin.
Spinful two-band model. The spinful two-band model provides the non-interacting, doubly degenerate c and v Bloch states of crystal momentum k that comply with the symmetry group of monolayer WTe 2 (the T′ structure is centrosymmetric and non-symmorphic). Within the 4D spin/orbital space, the Hamiltonian is a 4 × 4 Hamiltonian matrix, H QSH (k), whose off-diagonal elements are the spin-orbit interaction terms. Here the orbital degree of freedom identifies the c and v Bloch states at Γ, whose energies are 'inverted' with respect to the usual order of bulk semiconductors. These two states have been variously identified in the literature as a pair of orbitals having either opposite 15,25 or like 29,44 parities under spatial inversion, leading to two different forms of the spin-orbit interaction: we label the corresponding model Hamiltonians, respectively, as H QSH,1 (k) and H QSH,2 (k). We find that the charge order is hidden in the condensate of indirect excitons, regardless of the model. The reason is that the screw-axis rotation (around the W atom chain direction) maintains the same form within the spin/orbital space, entangling the degrees of freedom of its eigenstates, labelled by ζ = ±i. This is pivotal to the discussion of the main text.
In the following we detail the Hamiltonians H QSH,1 (k) and H QSH,2 (k). The first model, H QSH,1 (k), is taken from refs. 15,25 with minor adjustments. The Hamiltonian, written in the present reference frame with the W atom chain parallel to the x axis, reads where v 2 is the spin-orbit coupling parameter, τ x , τ y , τ z and σ x , σ y , σ z are 2 × 2 Pauli matrices in orbital and spin space, respectively, and the 2 × 2 unit matrices are 1 τ and 1 σ . The diagonal matrix elements, ε l (k) with l = u, g, are the band energies in the absence of spin-orbit interaction, which are inverted at Γ and cross at the points ±k Λ of the ΓX line in the Brillouin zone. The energies ε l (k x , k y ) are even with respect to both k x and k y axes. The corresponding Bloch states transform like p x and d xz orbitals at Γ. The functional dependence of ε l on k differs from the effective-mass expression given in ref. 25 , but its precise form is irrelevant to the discussion of the main text, solely based on symmetry arguments. We discard the off-diagonal spin-orbit term linear in σ x considered in ref. 25 , the Hamiltonian being now even in k y , H QSH,1 (k x , k y ) = H QSH,1 (k x , −k y ). This choice agrees with the evidence that the spin-orbit field lies in the x-z plane 29 , as also predicted by ref. 45 . Furthermore, we find that the spin-orbit term proportional to σ y provides a good matching with the strongly anisotropic DFT bands, as we checked through comparison with our own first-principles calculations. In the spin/orbital space, the inversion operator reads I = −τz ⊗ 1σ and the screw-axis rotation around the x axis is C2x = iτz ⊗ σx exp(ikxa/2), with a being the lattice constant in the direction parallel to the W chains. The second model, H QSH,2 (k), is taken from ref. 29 and builds on the four-band tight-binding Hamiltonian proposed in ref. 44 to improve the matching between model and DFT bands. The c and v orbital states are Wannier functions, respectively an antibonding combination of d x 2 −y 2-type orbitals localized on W atoms (energy ε W (k)) and a bonding superposition of p x -type orbitals localized on Te atoms (ε Te (k)). These Wannier functions have the same parities under inversion but opposite parities under the screw-axis two-fold rotation, like the Bloch states at Γ of our own DFT calculations. The Hamiltonian is where λ SO > 0 is the spin-orbit coupling parameter, which is independent of k. For the sake of simplicity, here we have neglected the spin-orbit term proportional to σ z proposed by ref. 29 . Within the envelope function approximation, the band energies ε W (k) and ε Te (k) are provided by the tight-binding calculation of ref. 44 and are even with respect to both k x and k y axes. The inversion operator now reads I = 1τ ⊗ 1σ whereas the screw-axis rotation around the x axis is again Eigenstates of the screw-axis rotation. The eigenvectors of H QSH,1 (k) that were explicitly given in ref. 25 are spin-polarized along the direction perpendicular to the W atom chains. Here we use the notation |k, λ⟩ to identify these eigenvectors, which belong to either the c or v band and have spin polarization λ = ↑,↓. In the main text we introduced an alternative, equally legitimate set of eigenvectors, which are simultaneous eigenstates of H QSH,1 (k) and C 2x (the latter with eigenvalues ζ = ±i), by applying a unitary rotation of the basis for any wavevector k. Explicitly, |k, ζ = −i⟩ = (i |k, ↑⟩ + |k, ↓⟩)/ √ 2, C2x|k, ζ = +i⟩ = i exp(ikxa/2) |k, ζ = +i⟩, as may be checked by direct substitution (followed by rotation of the original frame axis of ref. 25 ). Importantly, the states ζ = ±i are not spin-polarized (except at k = 0), because the spin and orbital degrees of freedom are now entangled. Note that |k, ζ = −i⟩ and |−k, ζ = +i⟩ are time-reversal mates, with Θ |k, ζ = −i⟩ = i |−k, ζ = +i⟩, the time-reversal operator being Θ = i1τ ⊗ σyK (K is the complex conjugation operator). The simultaneous eigenvectors of H QSH,2 (k) and C 2x that we use to build the excitonic insulator ground state (within the envelope function approximation) are as follows. The c band state with k = (0, k Λ ) and ζ = +i is (−1 + i)/(2 Here the first and third (second and fourth) components of the 4D vector, [u W,↑ , v Te,↑ , u W,↓ , v Te,↓ ], correspond to a spinor whose orbital part is a Wannier function localized on W (Te) atoms in the crystal unit cell.

Condensate of indirect excitons and charge/spin order.
We assess the charge (spin) order of a permanent condensate of indirect excitons of momentum q = (±k Λ , 0) within the multivalley framework developed in ref. 10 . This approach, in turn, relies on the scheme to decouple the equations of motion for Green functions of ref. 46 . The theory, which deals with spinless electrons, may be straightforwardly generalized to spinors labelled by the screw-axis symmetry =±i. Indeed, for any given electron and hole species ζ and ζ′, the structure of the equations of motion for Green functions that govern the condensate component ⟨ c + ζ v ζ′ ⟩ remains the same, because ζ electrons pair with ζ′ holes only. As far as pairing is concerned, ζ spinors behave as if they were spinless fermions. Therefore, the k-dependent spinless Bogoliubov-Valatin-like creation operator that defines the excitonic insulator ground state in equations (20) and (21) of ref. 10 is simply replicated for ζ = ±i, provided one specializes equation (21) to the present case of two valley components and chooses the condensate phases as shown above (Excitonic insulator ground state from indirect excitons and synopsis of theory). Finally, we compute the expectation value of the charge (spin) density operator over the excitonic ground state, after making the spin/orbital structure of ζ Bloch states given in the previous section explicit, the derivation being lengthy but straightforward. As discussed above, to simply assess the occurrence of charge (spin) order without computing the density wave modulation intensity, it is not necessary to evaluate explicitly the coherence coefficients u 0 k and v 0 k that occur in equation (21) of ref. 10 (these are provided by the self-consistent gap in equation (3)). Fig. 5e. We obtain the results shown in Fig. 5e within the spinful two-band model H QSH,1 (k) described above, the band energies ε l (k) being parametrized through comparison with our own DFT calculations (l = u, g). The non-interacting ground state, in the absence of spin-orbit interaction, is taken to be a semimetal with a band overlap of 38 meV. The energy parametrization relies partly on the tight-binding model of ref. 44 (table III therein) and partly on ad hoc parameters. In detail, we correct the tight-binding energies by adding the terms 2 ( t ′ l − t ′′ l ) cos 2kx + 2t ′′ l cos 2 |k|, with t ′ g = 0.149 eV, t ′ u = 0.075 eV (ref. 23 ), t ′′ g = 0.049 eV and t ′′ u = 0.055. The spin-orbit parameter is v 2 = 1 × 10 14 Å s −1 , and the strength of the Coulomb interaction is fixed by the 2D polarizability, α 2D = 5.5 (the notation of ref. 25 ). To compute the chemical potential μ versus charge density n g in the many-body excitonic phase, we adapt the theory of ref. 25 , which deals with a condensate of direct excitons in an intrinsic semiconductor, to the case of a doped system, by means of a fully self-consistent calculation of both the excitonic order parameter, Δ X (k), and μ. Furthermore, we assume that the Coulomb interaction remains long-ranged for any doping value, which is supported by the evidence that charge carriers localize in a wide doping interval at low T. The free carriers populating the renormalized bands of the excitonic insulator are conventionally taken as non-interacting, which is the origin of the unphysical behaviour of μ for small, positive values of n g at T = 25 K (dμ/dn g is negative close to the axis origin in Fig. 5e).

Condensate of direct excitons and simulation of
We have checked that the observable effects related to the breaking of inversion symmetry, due to the condensation of excitons with q = 0, are negligible for WTe 2 , in contrast to the case of T′-MoS 2 . This is related to the limited extent of the excitonic order parameter, Δ X (k), along the Brillouin zone direction perpendicular to the ΓΛ line, which is in turn caused by the strong anisotropy of the non-interacting c and v bands close to Λ (compare Fig. 5a with Supplementary  Fig. 8). In particular, the real part of Δ X (k) is negligible, so the c and v bands, renormalized by electron-hole pairing, each remain doubly degenerate.

Data availability
The datasets generated during and/or analysed during the current study are available from the corresponding authors on reasonable request. Source data are provided with this paper.