Abstract
We present evidence that the twodimensional bulk of monolayer WTe_{2} 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 Vshaped 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 independentelectron 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 chargeordered 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 WTe_{2} exhibits strong correlations over a wide temperature range.
Main
In 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,16,17} as a twodimensional (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}, angleresolved photoemission^{20,21} and density functional theory (DFT) calculations^{20,21,22,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 firstprinciples calculations of exciton dispersion and Bohr radius support this conclusion. The insulating behaviour below 100 K suggests a chargeordered 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_{\rm{p}}^{  1}} \approx {\beta /\sigma} + {R_{\rm{c}}}}\), where β is a geometrical factor considerably larger than one and R_{c} is the contact resistance. The gateinduced 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 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 MIMderived 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 dropoff 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 gatedependent anisotropy, with the lowest conductivity occurring for pdoping 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 gateinduced 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 singleparticle picture, as can be shown by an elementary calculation: in such a picture, μ and n_{g} are related by \({n_{\rm{g}}} = {{{\int}_{  \infty }^{ + \infty }} {D\left( E \right)f\left( E \right){\rm{d}}E}}\), where D(E) 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 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 bestfit 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 singleparticle 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 gateinduced 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 singleparticle 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 gateindependent contribution to the conductance. The vertical shift cannot be accounted for by varying n_{x}; to illustrate this we also plot (green dashdotted 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 singleelectron spectrum, for example with overlapping bands. A more sophisticated treatment of the correlatedelectron 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 spin–orbit effects in a nonperturbative 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 centreofmass frame. The exciton radius is as small as 4 nm. This is comparable with the typical electron separation at the critical doping \({{{n_{\rm{ce}}}^{{  1/2}}} = {4.5\,{\rm{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 manybody 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 temperaturedependent 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 meanfield 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 timereversal 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 timereversal 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 twoband 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 selfconsistent pinning of μ (in the chargeordered 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 manybody gap being computed selfconsistently (Methods). Here the gap has a purely excitonic origin, as the starting noninteracting phase is a semimetal. The step in μ remains clearly visible up to 100 K, in contrast with the smeared profile of the independentelectron 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.
Methods
Groundstate calculations from first principles
We obtained the groundstate 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 normconserving pseudopotentials^{36} to include the spin–orbit interaction. We optimized the lattice parameters and atomic positions using the PBE exchangecorrelation 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 lowestenergy exciton within the framework of manybody perturbation theory^{37,38,39}, by solving the Bethe–Salpeter equation through the YAMBO code^{40,41} and including spin–orbit interaction in a nonperturbative 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 kpoint grid. We extrapolated the excitation energy of the lowest exciton with momentum q = 0 to the limit of a dense kpoint grid, as shown in Supplementary Fig. 7.
Excitonic insulator ground state from indirect excitons and synopsis of theory
From the firstprinciples solution of Bethe–Salpeter equation we find that the lowestenergy exciton with q = k_{Λ} (or −k_{Λ}) is threefold degenerate and separated by 20 meV from a nondegenerate 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 highsymmetry line, both c and v Bloch spinors, each doubly degenerate, may be chosen as the eigenstates of the twofold screwaxis 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 twoband model, which includes both spin and orbital degrees of freedom (Spinful twoband model). In addition to the rotational symmetry, we may classify the excitons according to their triplet or singletlike 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 twofold rotation, even if their actual spin polarization away from Γ is zero. Specifically, the exciton wavefunctions that are even with respect to the screwaxis rotation are
where +/− stands for singlet/tripletlike symmetry, \({{c_\zeta ^{\dagger}} \left( {{{\bf{k}}}} \right)}\) creates and \({v_\zeta \left( {{{\bf{k}}}} \right)}\) destroys an electron of momentum k and screwaxis symmetry ζ in the c and v band, respectively, 0〉 is the noninteracting ground state with all v states filled and c states empty, and ψ is the exciton wavefunction in reciprocal space in the centreofmass 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 electron–hole pair, \({\left\langle {{c_\zeta ^ +} {v_{\zeta \prime }}}\right\rangle}\), has finite magnitude ϕ_{ζζ′} and arbitrary phase θ_{ζζ′}; that is, the complex wavefunction of the condensate is ϕ_{ζζ′}exp(iθ_{ζζ′}) (here \(\left\langle \ldots \right\rangle\) is the average over the manybody 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 timereversal symmetry^{10}. Furthermore, both condensate and exciton wavefunctions share the same screwaxis and spinlike 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 tripletlike symmetry, one has
For the singletlike ground state, one replaces the 2 × 2 Pauli matrix σ_{z} with the identity matrix 1, the righthand 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 lefthand 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 twoband model
The spinful twoband model provides the noninteracting, 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 nonsymmorphic). Within the 4D spin/orbital space, the Hamiltonian is a 4 × 4 Hamiltonian matrix, H_{QSH}(k), whose offdiagonal 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 screwaxis 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 effectivemass 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 offdiagonal 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 firstprinciples calculations. In the spin/orbital space, the inversion operator reads \({{{\bf{I}}}} = {{ \tau _z} \otimes {1_\sigma}}\) and the screwaxis rotation around the x axis is \({{{{\bf{C}}}}_{2x}} = {{i\tau _z} \otimes {\sigma _x\exp (ik_xa/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 fourband tightbinding 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 screwaxis twofold 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 tightbinding calculation of ref. ^{44} and are even with respect to both k_{x} and k_{y} axes. The inversion operator now reads \({{{\bf{I}}}} = {{{\bf{1}}_\tau} \otimes {{\bf{1}}_\sigma}}\) whereas the screwaxis rotation around the x axis is again \({{{{\bf{C}}}}_{2x}} = {i\tau _z} \otimes {\sigma _x\exp (ik_xa/2)}\).
Eigenstates of the screwaxis rotation
The eigenvectors of H_{QSH,1}(k) that were explicitly given in ref. ^{25} are spinpolarized along the direction perpendicular to the W atom chains. Here we use the notation \({{\left {{{{\bf{k}}}},\lambda} \right\rangle} }\) 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,
with
and
as may be checked by direct substitution (followed by rotation of the original frame axis of ref. ^{25}). Importantly, the states ζ = ±i are not spinpolarized (except at k = 0), because the spin and orbital degrees of freedom are now entangled. Note that \({\left {{{{\bf{k}}}},\zeta =  i} \right\rangle}\) and \({\left {  {{{\bf{k}}}},\zeta = + i} \right\rangle}\) are timereversal mates, with \({{\varTheta}\left {{{{\bf{k}}}},\zeta =  i} \right\rangle = i\left {  {{{\bf{k}}}},\zeta = + i} \right\rangle}\), the timereversal operator being \({{\varTheta} = {i{\bf{1}}_\tau \otimes \sigma _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\sqrt 2 )}{[1,1,1, 1]}}\), and the v band state with k = 0 and ζ = +i is (approximately) \({{\left( {1/\sqrt 2 } \right)}{[0,  1,0,1]}}\). The states with ζ = −i as well as those with k = (0, −k_{Λ}) are obtained through timereversal and inversion transformations. Here the first and third (second and fourth) components of the 4D vector, \({[{u_{{\rm{W}}, \uparrow }},{v_{{\rm{Te}}, \uparrow }},{u_{{\rm{W}}, \downarrow }},{v_{{\rm{Te}}, \downarrow }}]}\), 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 screwaxis 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 \({\left\langle {{c_\zeta ^ +} {v_{\zeta \prime }}} \right\rangle}\) 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 kdependent spinless Bogoliubov–Valatinlike 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_k^0}\) and \({v_k^0}\) that occur in equation (21) of ref. ^{10} (these are provided by the selfconsistent gap in equation (3)).
Condensate of direct excitons and simulation of Fig. 5e
We obtain the results shown in Fig. 5e within the spinful twoband 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 noninteracting 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 tightbinding model of ref. ^{44} (table III therein) and partly on ad hoc parameters. In detail, we correct the tightbinding energies by adding the terms \({{2\left( {{t_l^\prime}  {t_l^{\prime\prime}} } \right){\cos 2k_x}} + {{2t_l^{\prime\prime}} \cos 2\left {{{\bf{k}}}} \right}}\), with \({t_g^\prime} = {0.149}\,{\rm{eV}}\), \({t_u^\prime} = {0.075}\,{\rm{eV}}\) (ref. ^{23}), \({t_g^{\prime\prime}} = {0.049}\,{\rm{eV}}\) and \({t_u^{\prime\prime}} = {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 manybody 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 selfconsistent calculation of both the excitonic order parameter, Δ_{X}(k), and μ. Furthermore, we assume that the Coulomb interaction remains longranged 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 noninteracting, 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).
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 noninteracting 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.
Code availability
Manybody perturbation theory calculations were performed using the codes YAMBO (http://www.yambocode.org/) and Quantum ESPRESSO (http://www.quantumespresso.org), which are both opensource software. Results for the twoband model were obtained through custom Fortran codes that are available from M.R. upon reasonable request.
References
Keldysh, L. & Kopaev, Y. Possible instability of semimetallic state toward Coulomb interaction. Sov. Phys. Sol. State 6, 2219–2224 (1964).
Cloizeaux, J. D. Exciton instability and crystallographic anomalies in semiconductors. J. Phys. Chem. Solids 26, 259–266 (1965).
Jerome, D., Rice, T. M. & Kohn, W. Excitonic insulator. Phys. Rev. 158, 462–475 (1967).
Halperin, B. I. & Rice, T. M. Possible anomalies at a semimetal–semiconductor transition. Rev. Mod. Phys. 40, 755–766 (1968).
Rohwer, T. et al. Collapse of longrange charge order tracked by timeresolved photoemission at high momenta. Nature 471, 490–493 (2011).
Kogar, A. et al. Signatures of exciton condensation in a transition metal dichalcogenide. Science 358, 1314–1317 (2017).
Lu, Y. F. et al. Zerogap semiconductor to excitonic insulator transition in Ta_{2}NiSe_{5}. Nat. Commun. 8, 14408 (2017).
Werdehausen, D. et al. Coherent order parameter oscillations in the ground state of the excitonic insulator Ta_{2}NiSe_{5}. Sci. Adv. 4, eaap8652 (2018).
Varsano, D. et al. Carbon nanotubes as excitonic insulators. Nat. Commun. 8, 1461 (2017).
Ataei, S. S., Varsano, D., Molinari, E. & Rontani, M. Evidence of ideal excitonic insulator in bulk MoS_{2} under pressure. Proc. Natl Acad. Sci. 118, e2010110118 (2021).
Eisenstein, J. P. & MacDonald, A. H. Bose–Einstein condensation of excitons in bilayer electron systems. Nature 432, 691–694 (2004).
Nandi, D., Finck, A. D. K., Eisenstein, J. P., Pfeiffer, L. N. & West, K. W. Exciton condensation and perfect Coulomb drag. Nature 488, 481–484 (2012).
Li, J. I. A., Taniguchi, T., Watanabe, K., Hone, J. & Dean, C. R. Excitonic superfluid phase in double bilayer graphene. Nat. Phys. 13, 751–755 (2017).
Liu, X., Watanabe, K., Taniguchi, T., Halperin, B. I. & Kim, P. Quantum Hall drag of exciton condensate in graphene. Nat. Phys. 13, 746–750 (2017).
Qian, X., Liu, J., Fu, L. & Li, J. Quantum spin Hall effect in twodimensional transition metal dichalcogenides. Science 346, 1344–1347 (2014).
Fei, Z. et al. Edge conduction in monolayer WTe_{2}. Nat. Phys. 13, 677–682 (2017).
Wu, S. et al. Observation of the quantum spin Hall effect up to 100 kelvin in a monolayer crystal. Science 359, 76–79 (2018).
Sajadi, E. et al. Gateinduced superconductivity in a monolayer topological insulator. Science 362, 922–925 (2018).
Fatemi, V. et al. Electrically tunable lowdensity superconductivity in a monolayer topological insulator. Science 362, 926–929 (2018).
Tang, S. et al. Quantum spin Hall state in monolayer 1T′WTe_{2}. Nat. Phys. 13, 683–687 (2017).
Cucchi, I. et al. Microfocus laserangleresolved photoemission on encapsulated mono, bi and fewlayer 1T′WTe_{2}. Nano Lett. 19, 554–560 (2019).
Zheng, F. et al. On the quantum spin Hall gap of monolayer 1T′WTe_{2}. Adv. Mater. 28, 4845–4851 (2016).
Ok, S. et al. Custodial glide symmetry of quantum spin Hall edge modes in monolayer WTe_{2}. Phys. Rev. B 99, 121105 (2019).
Song, Y.H. et al. Observation of Coulomb gap in the quantum spin Hall candidate singlelayer 1T′WTe_{2}. Nat. Commun. 9, 4071 (2018).
Varsano, D., Palummo, M., Molinari, E. & Rontani, M. A monolayer transitionmetal dichalcogenide as a topological excitonic insulator. Nat. Nanotechnol. 15, 367–372 (2020).
Jia, Y. et al. Evidence for a monolayer excitonic insulator. Nat. Phys. https://doi.org/10.1038/s4156702101422w (2021).
Lee, P. A. Quantum oscillations in the activated conductivity in excitonic insulators: possible application to monolayer WTe_{2}. Phys. Rev. B 103, L041101 (2021).
Kwan, Y. H., Devakul, T., Sondhi, S. L. & Parameswaran, S. A. Theory of competing excitonic orders in insulating WTe_{2} monolayers. Phys. Rev. B 104, 125133 (2021).
Zhao, W. et al. Determination of the spin axis in quantum spin Hall insulator monolayer WTe_{2}. Preprint at https://arxiv.org/abs/2010.09986 (2021).
Wang, P. et al. Landau quantization and highly mobile fermions in an insulator. Nature 589, 225–229 (2021).
Fei, Z. et al. Ferroelectric switching of a twodimensional metal. Nature 560, 336–339 (2018).
Shi, Y. et al. Imaging quantum spin Hall edges in monolayer WTe_{2}. Sci. Adv. 5, eaat8799 (2019).
Lupke, F. et al. Proximityinduced superconducting gap in the quantum spin Hall edge state of monolayer WTe_{2}. Nat. Phys. 16, 526–530 (2020).
Zittartz, J. Anisotropy effects in the excitonic insulator. Phys. Rev. 162, 752–758 (1967).
Giannozzi, P. et al. Advanced capabilities for materials modelling with Quantum ESPRESSO. J. Phys. Condens. Matter 29, 465901 (2017).
Hamann, D. R. Optimized normconserving Vanderbilt pseudopotentials. Phys. Rev. B 88, 085117 (2013).
Onida, G., Reining, L. & Rubio, A. Electronic excitations: densityfunctional versus manybody Green’sfunction approaches. Rev. Mod. Phys. 74, 601–659 (2002).
Hybertsen, M. S. & Louie, S. G. Electron correlation in semiconductors and insulators: band gaps and quasiparticle energies. Phys. Rev. B 34, 5390–5413 (1986).
Strinati, G. Application of the Green’s functions method to the study of the optical properties of semiconductors. Riv. Nuovo Cim. 11, 1–86 (1988).
Marini, A., Hogan, C., Grüning, M. & Varsano, D. Yambo: an ab initio tool for excited state calculations. Comput. Phys. Commun. 180, 1392–1403 (2009).
Sangalli, D. et al. Manybody perturbation theory calculations using the yambo code. J. Phys. Condens. Matter 31, 325902 (2019).
Elliott, J. D. et al. Surface susceptibility and conductivity of MoS_{2} and WSe_{2} monolayers: a firstprinciples and ellipsometry characterization. Phys. Rev. B 101, 045414 (2020).
Rozzi, C. A., Varsano, D., Marini, A., Gross, E. K. U. & Rubio, A. Exact Coulomb cutoff technique for supercell calculations. Phys. Rev. B 73, 205119 (2006).
Muechler, L., Alexandradinata, A., Neupert, T. & Car, R. Topological nonsymmorphic metals from band inversion. Phys. Rev. X 6, 041069 (2016).
Garcia, J. H. et al. Canted persistent spin texture and quantum spin Hall effect in WTe_{2}. Phys. Rev. Lett. 125, 256603 (2020).
Monney, C. et al. Spontaneous exciton condensation in 1TTiSe_{2}: BCSlike approach. Phys. Rev. B 79, 045116 (2009).
Acknowledgements
The experiments and analysis were supported primarily by NSF DMR awards nos. MRSEC 1719797 and EAGER 1936697. Materials synthesis at the University of Washington was partially supported by the Gordon and Betty Moore Foundation’s EPiQS Initiative, grant no. GBMF6759 to J.H.C. X.H. and Y.T.C. acknowledge support from the NSF under award no. DMR2004701 and a Hellman Fellowship award. The development of the devices was supported as part of Programmable Quantum Materials, an Energy Frontier Research Center funded by the US Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES), under award no. DESC0019443. On the theory side, S.S.A., D.V., E.M. and M.R. were supported in part by the MaX European Centre of Excellence (‘MAterials design at the eXascale’, www.maxcentre.eu) funded by the European Union H2020INFRAEDI20181 programme, grant no. 824143. D.V., E.M. and M.R. were also supported by the Italian national programme PRIN2017 no. 2017BZPKSZ ‘Excitonic insulator in twodimensional longrange interacting systems (EXCINS)’. M.P. was supported by the INFN Time2Quest project and by the European Union H2020MSCARISE project DiSeTCom, grant no. 823728. S.S.A., M.P., D.V., E.M. and M.R. acknowledge access to the Marconi supercomputing system based at CINECA, Italy, through the PRACE and ISCRA programmes.
Author information
Affiliations
Contributions
B.S., W.Z., T.P., Z.F. and E.R. made the devices and carried out the transport measurements. P.M. and J.H.C. provided the crystals. J.C. and X.X. performed the Raman measurements. X.H. and Y.T.C. carried out the MIM measurements. S.S.A., D.V., M.P., E.M. and M.R. developed the theory of excitonic effects and wrote the theoretical sections. D.H.C. led the experiments and their analysis. All authors contributed to all processes and to the manuscript.
Corresponding authors
Ethics declarations
Competing interests
The authors declare no competing interests.
Additional information
Peer review information Nature Physics thanks Vitor Pereira and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary information
Supplementary Information
Supplementary Figs. 1–8 and Discussion.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Sun, B., Zhao, W., Palomaki, T. et al. Evidence for equilibrium exciton condensation in monolayer WTe_{2}. Nat. Phys. 18, 94–99 (2022). https://doi.org/10.1038/s41567021014275
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1038/s41567021014275
Further reading

Topological excitons
Nature Physics (2022)