A potential all-electronic route to the charge-density-wave phase in monolayer vanadium diselenide

The transition metal dichalcogenides offer significant promise for the tunable realisation and application of correlated electronic phases. However, tuning their properties requires an understanding of the physical mechanisms underlying their experimentally observed ordered phases, and in particular the extent to which lattice vibrations are a necessary ingredient. Here we present a potential mechanism for charge-density-wave formation in monolayers of vanadium diselenide in which the key role at low energies is played by a combination of electron–electron interactions and nesting. There is a competition between superconducting and density-wave fluctuations as sections of the Fermi surface are tuned to perfect nesting. This competition leads to charge-density-wave order when the effective Heisenberg exchange interaction is comparable to the effective Coulomb repulsion. When all effective interactions are purely repulsive, it results instead in d-wave superconductivity. We discuss the possible role of lattice vibrations in enhancing the effective Heisenberg exchange during the earlier stages of the renormalisation group flow. Charge density waves are periodic modulations of the electron density which occur at low temperatures in many transition metal dichalcogenide monolayers, though the underlying mechanism is still a matter of debate. Here the authors use renormalisation group analysis to demonstrate that electron–electron interactions and Fermi surface nesting may play a prominent role in mediating the competition between the charge density-wave and superconducting states.

S ince the isolation and characterisation of graphene in 2004 1 , the field of two-dimensional materials has seen an explosion in research activity 2 , and a search has begun for twodimensional materials that can be tuned to exhibit a wider range of properties than graphene. Of particular interest in this regard are monolayers of the transition metal chalcogenides FeX (X = Se, Te) and transition metal dichalcogenides MX 2 (M = Ti, V, Nb, Mo, Ta, W; X = S, Se, Te) 3 . The transition metal dichalcogenides (TMDs) display an especially wide range of behaviours, including Mott-insulating, semi-metallic, charge-density-wave (CDW), excitonic, and superconducting phases. The development of van der Waals heterostructures made from two or more TMDs 4 is expected to further increase the range of strongly correlated physics that can be realised in this family of materials.
However, tuning the properties of TMDs requires an understanding of the way in which variations in microscopic parameters affect their phase diagrams. This, in turn, necessitates an understanding of the physical mechanisms that underlie the experimentally observed ordered phases. For several of the ordered states of monolayer TMDs, especially the CDW phases, the mechanism remains the subject of debate.
Many TMDs exhibit CDW phases with rather high critical temperatures, which are often further enhanced in the monolayer limit 5 . One well known route to CDW formation is via Fermi surface nesting: here sections of the Fermi surface lie parallel to each other, giving an enhanced particle-hole susceptibility at a non-zero wavevector Q 6-10 . This is an inherently electronic mechanism. However, there are other candidate mechanisms for the CDW phases in the TMDs, including the Rice-Scott saddle band mechanism 11 , the softening of phonon modes 12 , and a mechanism based on the transition to an exciton insulator 13 . It is often assumed that there must be significant influence of the phonon sector in driving the CDW transitions in monolayer TMDs; however, recent evidence of electronically driven densitywave transitions in other correlated materials 14 motivate us to ask whether a mechanism dominated by the electronic sector might be operative here as well.
We focus on the 1T structural isomer of vanadium diselenide, VSe 2 , in the monolayer limit. Theoretical and experimental attempts to determine the low-temperature Fermi surface of this material do not all agree. Several studies show column-like Fermi surface pockets protruding from the edge of the Brillouin zone [15][16][17][18] ; others show a Fermi surface with large triangular pockets around the K and K 0 points of the Brillouin zone with an additional small Fermi surface pocket at the Γ point 19,20 . Which of these Fermi surfaces is realised appears to depend on the exact position of the chemical potential with respect to a van Hove singularity in the band structure 20 . Such singularities are usually associated with an enhancement of the susceptibilities to various forms of ordered phase, with superconductivity typically dominant 21,22 .
This variation in the predicted Fermi surface leads to a disagreement over the predicted Q-vector of any CDW, and thus also over the reconstructed unit cell. Some studies propose a Qvector perpendicular to the Brillouin zone edge 19,23,24 , in the k y direction; however, others propose alternative nesting vectors parallel to the Brillouin zone edges 17,25 . These studies agree on a renormalisation to flat Fermi surface sections in the lowtemperature and low-dimensional limit.
In this article we consider an idealised model of monolayer 1T-VSe 2 . For definiteness, we assume column-like Fermi surfaces 17,25 , though the patch scheme we employ should also be applicable to the triangular Fermi surface case with appropriate modifications to intra-and inter-pocket scattering and the definitions of superconducting symmetries. Starting from an effective model for the intermediate-energy (~10 meV) physics in the conduction-electron sector, we implement a renormalisation group (RG) analysis retaining both particle-particle and particlehole channels. We are thus able to capture the interplay of superconducting and density-wave fluctuations, and the effect of Fermi surface nesting on both 26,27 , as the eventual ordered state is approached.

Results
Renormalisation group flow equations. (For the meaning of the coupling constants and other details of our model, please see the 'Methods' section below.) We define the RG flow parameter y ¼ which diverges to infinity as Ω → 0. Introducing the rescaled interaction parameters g i ! k c 2π 2 g i we perform a oneloop RG analysis including terms that contribute with a divergent susceptibility at low energies 26,27 . We find the following RG flow equations: _ g 4 ¼ Àð1 À βÞg 2 4 À g 2 5 À g 2 6 ; ð4Þ _ g 6 ¼ Àg 3 g 5 À g 4 g 6 þ 2d ε ðyÞðg 1 g 5 þ g 2 g 6 À 2g 1 g 6 Þ; ð6Þ where _ g i denotes the derivative dg i dy . The y-dependence of the couplings g i has been suppressed for brevity. The function d ε (y) is defined as i.e. it measures the relative strengths of the Q 1 particle-hole susceptibility and the q = 0 particle-particle susceptibility in terms of the flow parameter y and nesting parameter ε. Solving these differential equations, we find that the couplings diverge at a critical value of y. In order to allow a numerical solution we stop the flow when the largest of the couplings g i becomes equal to 1; this defines a critical value y = y c . At this point a subset of the couplings have already become several orders of magnitude larger than their initial values, signalling the breakdown of our perturbation theory and the onset of order. The finite critical value y c is an artefact of the one-loop RG; higherloop corrections should shift the divergence to y c → ∞. We consider that no phase transition occurs if no coupling has reached 1 by the time y = 1/(N 0 U eff ), where N 0 = k c /(2π 2 ) is the density of states.
In the limit y → 0 with ε small, one may show from (7) that d ε % 1 À ε 2 3log 2 . In the large-y limit d ε (y) takes the form d ε (y → y c ) = ((1 − ε 2 )e y )/((1 + e y )(1 + ε 2 e y )). We therefore use the following approximation to d ε (y): which interpolates between the y → 0 and y → y c limits. The exact form of this interpolation function is selected for calculational convenience and is not unique; however, the properties of the RG flow should not depend strongly on the details of the choice we make here.
Initial conditions. The initial conditions for the couplings are approximated by We find the initial conditions at y = 0 for the couplings to be g 0 1 % g 0 3 % g 0 6 % V 1 and g 0 2 % g 0 4 % g 0 5 % V 2 in our approximation. The effect of this approximation is to split the solutions into three regions: , all couplings attractive. In a more general microscopic model the values of the couplings g 0 i would be independent. The mapping between the intermediate-scale effective couplings J eff and V eff and the RG couplings V 1 and V 2 is illustrated in Fig. 1.

.
We refer to the possible superconducting symmetries using their continuum analogues, despite the fact that our system is on a lattice and we furthermore only utilise a discrete set of patches. To make the meanings of these order parameters clear, we note that the s-wave eigenvector predicts an isotropic gap, while the dwave eigenvector leads to four nodes on each Fermi surface pocket. The p-wave and f-wave eigenvectors each give two nodes per pocket; however, the p-wave order parameter naïvely changes sign twice as a function of angle in the Brillouin zone, whereas the f-wave order parameter changes sign six times.
Due to our patch approximation we can predict neither the relative phases of the superconducting order parameter between pockets nor which vector(s) Q i will form the CDW. To calculate the latter, a multi-component order parameter theory is required 25 .
Given the divergence of the couplings at y c we introduce the asymptotic form g i = G i /(y c − y). As y → y c we can express the divergences of order parameter susceptibilities in the power-law CDW }. The exponents are given by the following equations: Due to the nature of our patch scheme, ferromagnetic instabilities cannot be investigated: they require the full Fermi surface to calculate the susceptibilities. Ferromagnetic phases have been observed experimentally in monolayer VSe 2 28 . However, there is evidence to suggest that ferromagnetism is suppressed near the CDW phase as our nested approximation would suggest 29 .
Considering the alternative triangular Fermi surface case, the definitions of intra-vs. inter-pocket scattering have to be altered. This does not change the CDW nesting vectors; however, the fwave superconductivity would be replaced by an s ± -like order parameter.
Phase diagrams. Solving (1-6) numerically with the initial conditions g i ðy ¼ 0Þ ¼ g 0 i , and utilising the definitions of the divergent susceptibilities, we can investigate the phase diagram of the model. When the effective interaction is of pure contact form, for which V eff = J eff = 0, only two instabilities are predicted: s-wave superconductivity for an initially attractive interaction and dwave superconductivity for an initially repulsive one. For V eff and J eff non-zero, the phase diagrams for a range of nesting strengths (ε = 10 −1 , 10 −2 , 10 −3 ) are plotted in Fig. 1. When all effective interactions are initially repulsive the predicted instability is again Fig. 1 Phase diagrams of our extended Hubbard model for monolayer 1T-VSe 2 . a An illustration of the mapping between the effective nearest-neighbour Coulomb repulsion V eff and the effective Heisenberg exchange interaction J eff and the coupling constants V 1 and V 2 , defined in (9). b-d Calculated phase diagrams for our model of monolayer 1T-VSe 2 . The assumed degree of Fermi surface nesting increases from left to right (ε = 10 −1 , 10 −2 , 10 −3 ), and the parameter β that represents the finite length of the nested sections is set to 1/2. As well as the metallic phase, we find regions of s-wave superconductor (s-SC), p-wave superconductor (p-SC), d-wave superconductor (d-SC), and a charge density-wave with wavevector Q 1 (Q 1 −CDW).
to d-wave superconductivity. As some of the initial effective interactions become attractive, regions of s-wave and p-wave superconductivity arise. As the nesting strength is increased, these regions become occupied by a CDW phase.

Discussion
The tuning of the Fermi surface nesting is a control parameter in our analysis. Our analysis is therefore complementary to that of Jang et al. 25 , and predicts CDW formation without any meanfield assumption, taking into account the competition of superconducting and density-wave fluctuations.
Our CDW phase emerges in the limit of strong nesting. There is evidence in the literature that supports this limit: for example, taking into account self-energy corrections to the patch dispersions, a flow to perfect nesting is predicted by some previous RG calculations 25,30,31 . Other analyses, however, show 10,32 a Fermi surface that retains a power-law curvature at the hot-spots and is not perfectly nested. However, given the relatively robust nature of our predicted CDW phase, we do not expect the incorporation of such curvature effects to alter our resulting phase diagrams qualitatively.
It is also worth noting that in our analysis, following Jang et al. 25 , we have-for technical reasons-considered hot-spots displaced from the boundary of the Brillouin zone. It is of course important to ask what effect this has. Moving our four hot-spots on pocket 1 towards the zone edge would suppress those fluctuation channels in which the converging hot-spots have different phases, i.e., the d-and f-wave superconducting channels. However, since these are not the channels in direct competition with the Q 1 CDW phase, we do not expect major alterations in our resulting phase diagrams in the range where the CDW phase occurs.
It might seem surprising that a strong effective exchange coupling J eff should lead to CDW formation, since one might have expected that it would favour a magnetically ordered state. However, this intuition is too 'local moment': in our itinerant model, where the kinetic energy scales dominate over the bare interaction energies, the development of order is determined by the interplay between the different channels of fluctuation as the temperature is lowered rather than by the mean-field ground state of any single interaction term.
The fact that some of the intermediate-scale effective interactions must be attractive for a CDW phase to be favoured is an interesting result in the context of a purely electronic calculation. It is well known that electron-phonon interactions lead to an effective attractive interaction between electrons. Whether the inclusion of such electron-phonon interactions in the high-energy flow is necessary to arrive at the values of J eff that provoke a charge-density-wave transition, or whether these could be reached by a (theoretical) purely electronic high-energy flow, is beyond the scope of this calculation to determine.
The Q 1 CDW wavevector is favoured as the chosen instability, even with the artificial enhancement of the q 1 channel due to lack of curvature corrections to the dispersions. Thus this behaviour again agrees with that seen by Jang et al. 25 , and gives a potential all-electronic mechanism for CDW formation in the monolayer TMDs.
To analyse the effect of Fermi surface nesting on CDW formation in the TMDs, we have performed an RG analysis of an effective extended Hubbard model for monolayer 1T-VSe 2 , retaining both particle-particle (superconducting) and particle-hole (density-wave) channels. In the region of intermediate-energy effective coupling strengths where some or all of the effective two-particle interactions are attractive, regions of superconductivity give way to CDW order as the strength of Fermi surface nesting is increased.

Methods
Since the Fermi surface of 1T-VSe 2 is derived from one of the vanadium dbands 20,33 , we adopt a single-band model to describe the physics of monolayer 1T-VSe 2 17,25 . We separate the development of correlations as the temperature is lowered into two different regimes: high-energy flow (from scales of~1 eV down to~10 meV), and low-energy flow (below scales of~10 meV). We know from experiment that the high-energy flow tends to increase the degree of nesting of opposite sides of the Fermi surface pockets; it will also renormalise the on-site Coulomb repulsion U, the nearest-neighbour Coulomb repulsion V, and the nearest-neighbour magnetic exchange J from their microscopic values. Since we cannot describe the high-energy flow quantitatively, we take the nesting parameter ε and the effective interactions at a scale of~10 meV -U eff , V eff , and J effas input parameters of our electronic model. The resulting Hamiltonian is Here the operator c y iσ creates an electron with spin-projection σ on vanadium site i; n iσ ¼ c y iσ c iσ is the number operator for electrons on vanadium site i with spinprojection σ, while S i ¼ 1 2 P σσ 0 c y iσ τ σσ 0 c iσ 0 is the operator for the spin on site i, where τ ¼ ðτ x ; τ y ; τ z Þ T is the vector of Pauli matrices. t eff,ij denotes the effective hopping matrix elements for our single-band model of VSe 2 , which already incorporate the increase in nesting resulting from the high-energy flow. U eff and V eff are the effective strengths of the on-site and nearest-neighbour parts of the Coulomb repulsion respectively, and J eff is the effective Heisenberg exchange coupling. As noted above, these effective couplings, applicable at scales~10 meV, are renormalised from those in the microscopic model due to the high-energy stage of the flow. 〈i, j〉 indicates that the sum runs over all pairs of nearestneighbour sites.
At low energies it is sufficient to linearise the non-interacting dispersion relation around the Fermi level, and thus we do not need to know the precise form of the matrix elements t eff,ij . A schematic non-interacting Fermi surface is shown in Fig. 2. As discussed above, nested sections of the Fermi surface arise at lower temperatures. Since these dominate the relevant susceptibilities, we can safely use a simplified form of the dispersion relation that agrees with the true dispersion in these nested regions.
We utilise the patch scheme of Jang et al. 25 . This scheme consists of twelve patches that lie on sections of the Fermi surface that become nested at low temperatures, as shown in Fig. 2. The absolute wavevector of the centre of patch 1 + is denoted K 1+ , and similarly for the other patches. In each patch we linearise the dispersion relation, i.e., we write the single-electron energy (measured with respect to the Fermi energy) as a linear function of the components of k, the wavevector measured relative to the centre of the patch. For the four patches labelled '1', this gives ξ ± 1 ¼ ±k x þ εk y and ξ ± 1 ¼ Àξ ± 1 , in units where both ℏ and the Fermi velocity v F are set to 1. The parameter ε controls the nesting of the Fermi surface, with the limit ε → 0 corresponding to perfect nesting. We use only the bare effective dispersions in our calculations as fermion self-energy corrections are independent of the renormalisation of interactions at one-loop 34 . a Schematic Fermi surface of monolayer 1T-VSe 2 with non-nested Fermi surfaces. In the lowtemperature limit sections of the Fermi surface become nested. We adopt the notation of Jang et al. 25 to describe the patch scheme and nesting vectors. b Schematic change in one pocket of the Fermi surface of monolayer 1T-VSe 2 as the temperature T is lowered. c A zoomed view of the left-hand side of that pocket in our linearised approximation.
The effective dispersions on the second and third pockets n = 2, 3 may be obtained from a similar expression, ξ ± n ξ ± ðk ðnÞ x ; k ðnÞ y Þ ¼ ±k ðnÞ x þ εk ðnÞ y , where the wavevector k (n) is obtained by an appropriate rotation: together with the relation ξ ± n ¼ Àξ ± n . We can then use these dispersions to calculate the particle-particle and particlehole susceptibilities for all possible nesting vectors between patches with Gðω; kÞ ¼ ðiω À ξ k þ μÞ À1 . The range of integration is ω ∈ (−∞, ∞) and k x , k y ∈ (− k c , k c ), where k c is an ultraviolet momentum cutoff 27 . In the presence of near-perfect nesting, these susceptibilities both have similarly strong divergences. The nature of the eventual low-temperature ordered state is thus determined by a rather intricate competition between various superconducting and density-wave fluctuation channels, which must be tracked carefully as the energy scale is lowered from~10 meV where we take the low-energy flow to begin. The complete particle-hole susceptibility at wavevector The Fermi surface nesting parameter ε cuts off the Ω → 0 divergence of the logarithm in this channel, and the height of the Ω = 0 peak in the susceptibility reduces as ε is increased. By contrast, the particle-particle susceptibility at zero momentum in the low-energy limit has the usual logarithmic dependence, independent of ε, Π 0 pp ðΩÞ % k c 2π 2 log k c Ω . Here we have discarded contributions from non-divergent arctan terms as they are negligible as Π 0 pp ðΩÞ becomes large at low energies.
The particle-particle susceptibility Π q 1 pp ðΩÞ with q 1 = K 1+ + K 1− is logarithmically divergent and dependent on the nesting parameter ε; indeed, Π q 1 pp ðΩÞ ¼ Π Q 1 ph ðΩÞ. The particle-hole susceptibility Π 2K 1þ ph ðΩÞ is always perfectly nested for the case of linear dispersion. However, the nested sections of the VSe 2 Fermi surface are finite in length and there will be curvature corrections to the dispersion which will cutoff the divergence of the integral. We therefore introduce an additional parameter β with 0 ⩽ β ⩽ 1 to reduce the magnitude of this susceptibility and emulate the effect of finite length nested sections: Π 2K 1þ ph ðΩÞ ¼ βΠ 0 pp ðΩÞ. Interactions between Fermi surface patches belonging to different pockets do not give divergent contributions, since the particle-particle bubble has non-zero q and there is no particle-hole nesting between patches on separate pockets. Therefore in our low-energy model we retain only one of the Fermi surface pockets, thereby reducing the number of patches to four. This greatly simplifies our effective Lagrangian; however, we lose information about the relative phase of the superconducting order parameter between different Fermi surface pockets and the competition of particle-hole nesting vectors.