Introduction

Transition metal dichalcogenides are among the most exciting hosts for ordered phases, such as superconductivity and charge/spin modulations1,2,3,4,5 because of their cooperative interactions6. In fact, their coexistence is counter-intuitive, being collective electronic modes which build an excitation gap, hence requiring a large number of electrons available at the Fermi level. Since early works7,8,9, the transition to a charge/spin ordered state driven by the electron gas at the Fermi level has been questioned, partly because of the difficulties in explaining the intricate coexistence of superconductivity with charge density waves (CDWs)2,4,5,10,11,12. Among the proposed explanations, the electronically driven mechanism was ruled out13,14,15,16,17,18,19 in favour of a momentum-assisted mechanism10,18,20,21,22,23,24, which is likely influenced by external conditions, such as applied fields12,25 or epitaxially induced strain26,27, and chemical28,29,30 or gate27,31 doping.

Further intricacies in the occurrence of these collective electronic excitations arise at the surface and single layers due to a reduced symmetry. While the bulk, 2H type, is characterised by a P63∕mmc symmetry, with centrosymmetry and inversion symmetry operations, the layer at the surface has no inversion symmetry and lacks a van der Waals bonded layer; its structural type remains 2H, although high reactivity to (accidental) dopants and diffusion into interstitial regions affect experimental investigations, modifying its structure and symmetry operations; the single layer is characterised by a D3h symmetry, with an inversion symmetry operation but lack of centrosymmetry and no van der Waals bonded layers on either sides. In passing from the bulk to thin layers and down to a single layer, enhanced effects from fluctuations and lack of van der Waals bonded layers may determine dramatic changes in the structural and electronic degrees of freedom32,33,34,35,36,37,38,39,40,41. Due to the fact that CDW consists of periodic lattice distortions coupled to charge modulations, structural properties may modify such collective excitation dramatically, and in turn its interplay with superconductivity37,42. In the case of 2H-NbSe2, the CDW is recognised to arise from a phonon instability43, inducing an electronic reconstruction in a sizable energy range23 in the bulk and at the surface, although single layers show a CDW gap centred at the Fermi level37. Electronic and magnetic texture may be layer-resolved44,45, contributing to a complex scenario to explain the thickness-dependent properties of these materials. Reduced dimensionality allows for different ways of tuning the crystal and the electronic structure. The competition between the trigonal prismatic (2H) and the octahedral (1T) phases becomes stronger; as a consequence, the 1T crystal phase, unstable in the bulk46, may be synthesised47 or induced by external sources48; the resulting charge order is not the usual star of David phase, unlike TaSe249,50 and TaS251,52. In the 2H-type, 4 × 453,54 and/or 2 × 254 modulations, linked to a stripe phase at the surface, are reported. Distinct CDW structures, supporting the widely known CDW incommensurate character in 2H-NbSe255 are observed, further speculating on their topological connection. A change from the 3 × 3 to the 4 × 4 modulation was previously predicted by ab initio studies on 2H/1H-type NbSe220,56, but an early scanning tunnelling microscopy/spectroscopy (STM/STS) study37 confirms that the 3 × 3 modulation dominates in single layers deposited on bilayer graphene; on the other hand, a recent theoretical study57 showed that within a 3 × 3 modulation, uniaxial strain leads to a triangular-stripe CDW transition in 2H-NbSe2. However, phase transitions remain often latent because they are suppressed by growth/synthesis conditions, choice of the substrate and/or quantum fluctuations. In this respect, discerning intrinsic from extrinsic characteristics is crucial to understand these phases and finding possible routes to their manipulation.

The present work focusses on 1H-NbSe2 single layers under isotropic biaxial strain. In order to shed light on the competition and coexistence of various periodicities and structures, we first perform phonon calculations, suggesting the wave-vectors of the structural instabilities. Thereafter, the suggested periodicities are investigated with a full set of total energy calculations, yielding both the energy hierarchy of the possible CDW structures and their charge distributions. A 3q−1q (triangular-stripe) transition occurs in the 4 × 4 modulation, whereas it is incipient in the 3 × 3 modulation. Although our model represents only an approximation of detailed experimental conditions, e.g. due to the fact that we do not include a substrate, we believe it identifies the most important intrinsic characteristics of single layers under isotropic biaxial strain. Moreover, our results are likely applicable also to thin films, due to the fact that CDW properties depend on the structural properties more strongly than on the electronic properties, as well as due to the above-mentioned high reactivity to accidental doping which may affect the structural features at the surface.

Methods

Our results are obtained by ab initio calculations within the formalism of the density-functional theory (DFT). The projected augmented wave (PAW) method with Perdew–Burke–Ernzerhof (PBE) pseudopotentials58,59, as implemented in the Quantum ESPRESSO suite60,61 and the Vienna Ab-initio Simulation Package (VASP), is used. The relaxation of the unit cell and the computation of its relative phononic spectra are run in the former code, while structural relaxation, total energy calculations and charge distributions are performed with the latter. The exchange-correlation functional is treated within the generalised gradient approximation using the PBE parametrisation62,63. The pseudopotentials used in the simulations provide explicit treatment of the following valence electrons: 4s24p65s15d4 for Nb and 4s24p4 for Se. The cutoff energy of the plane waves for the unit cell (QUANTUM ESPRESSO) is ~1142 eV and that for the supercell calculations is 500 eV; prior to the phonon calculations, an energy cutoff is also applied for the charge density and potential, of ~4354 eV. The energy tolerance on the electronic loops for the calculations in the unit cell is 1.4 × 10−9 eV, whereas it is 10−6 eV for the calculations in the supercells (10−7 eV for accurate determination of the electronic properties, such as density of electronic states and charge distributions). A conjugate gradient algorithm is employed for structural relaxation in all cases, with a tolerance in energy of 10−4 eV per unit cell and a residual force of ~10−3 eV Å−1 on each atom. The calculations of phonons are performed sampling the first Brillouin Zone with a 45 × 45 × 1 k-grid (for the wavefunctions) and a 15 × 15 × 1 q-grid (for the force constants). The formation of CDWs is analysed sampling the first Brillouin Zone by means of the 30 × 30 × 1, 23 × 23 × 1, 20 × 20 × 1, 17 × 17 × 1 and 15 × 15 × 1 grids of k-points for the 2 × 1, \(\sqrt{7}\times \sqrt{7}\times 1\), 3 × 3 × 1, \(\sqrt{13}\times \sqrt{13}\times 1\) and 4 × 4 × 1 supercells, respectively, keeping approximately an equal distance between k-points.

Finally, we also computed some quantities offering a more direct comparison to experimental data. Band structures were calculated using method of Popescu and Zunger64. STM simulations were also performed, by means of the BSKAN code65,66. The revised Chen method67 with an electronically flat and spatially spherical tip orbital has been used, which is equivalent to the Tersoff–Hamann model of electron tunnelling68. The reported STM images are in constant current mode (with the maxima of the current contours at 5.8 Å for all structures for a better comparability). Additional technical details on band structures and STM spectra are provided in the Supplementary Information.

Results and discussion

Structural relaxation of the unit cell was performed for three values of the lattice constant, 3.41, 3.45 and 3.49 Å, which aim at modelling the system under compressive strain (~1%), no strain and tensile strain (~1%), respectively. The precise value of the computed in-plane stress is 7.2 × 10−3, 1.1 × 10−3 and −4.5 × 10−3 eV Å−3, respectively. The small residual compressive strain present in the no-strain case is due to the choice of a standard literature value for the equilibrium lattice constant of the bulk, in line with the work of Fang et al.69. Since the discrepancy is very small, this detail does not affect the conclusions of the present work. A recent ab initio study19 reports a lattice constant of 3.458 Å, in agreement with the residual strain found for \(\left| \bf a \right| =\) 3.45 Å in our calculations. The relaxed structures were used to compute phonon spectra, whose imaginary values identify regions of structural instability, thus speculating on the periodicity of the supercells describing the stable crystal. These supercells are then investigated by total energy calculations to confirm (or dismiss) the relevance of a periodicity and to evaluate the energy hierarchy of a set of symmetries within each periodicity. The lattice reconstructions considered on the basis of the analysis of phonon spectra are of the type \(\sqrt{n}\times \sqrt{n}\)70 and n × m (suggested by discussion of results from STM55).

The phonon dispersions are represented as two-dimensional plots (Fig. 1), with positive (negative) values of (ω(q))2) shown in blue (red). Compressive strain, equilibrium, and tensile strain are shown from left to right (Fig. 1a–c). The deepest troughs of the single layer phonon dispersion appear off the ΓM line, as a consequence of the non-centrosymmetric character of the single layer (1H). The instabilities have to be interpreted as perturbations which drive the system away from its most symmetric configuration, providing only an indication, rather than a determination, of the resulting periodicities. A \(\sqrt{13}\times \sqrt{13}\) modulation is expected from the lowest values of the phonon dispersion, see especially Fig. 1b, where the troughs appear in couples at symmetric positions with respect to the ΓM high symmetry line. Each couple of instability represent two periodic lattice distortions with opposite chirality. This result supports early71 and recent72 STM observations. In particular, a 2H-1T structural phase separation is reported in refs. 72,47,48. Further, we note similarities between the phonon dispersion in the 1H-type and that in the 1T-type70. Nevertheless, these instabilities may couple to give rise to a 3 × 3 modulation, which is the most observed. In a similar fashion, instability regions around M may couple to give rise to a 2 × 2 (or even a 2 × 3) modulation. Instabilities are expected to dominate along ΓM in a 2H-type symmetry (because the crystal is overall centrosymmetric).

Fig. 1: Phonon dispersion of NbSe2 single layer under strain.
figure 1

Two-dimensional phonon dispersion for the NbSe2 single layer unit cell as two-dimensional density plots, under compressive strain a, at the equilibrium lattice constant b and under tensile strain c. The plots, as function of kx and ky, in units of 2π∕a, show only the branch with the lowest values of ω(q)2; red (blue) represent regions of negative (positive) values of ω(q)2, see the colour bar, where the scale is given in tens of cm−1.

Full size image

Compressive strain reduces the instabilities along the MK line and the minima of (ω(q))2 dominate around 2/4 ΓM (off the high symmetry line) (see Fig. 1a), pointing towards a 4 × 4 modulation. Conversely, a zone of instability centred around, but not including, M becomes dominant under tensile strain, pointing towards a 2 × 2 modulation; further local minima of (ω(q))2 appear at 2/4 ΓM (Fig. 1c), justifying recent reports54 for which both 2 × 2 and 4 × 4 periodicities arise under tensile strain.

After having found the most likely periodicities for the occurrence of periodic lattice distortions by the calculation of phonon spectra, we investigate a set of structures for each periodicity by total energy calculation; in particular, the distorted structures we start with have fewer point symmetry operations with respect to the full symmetric structure, and model a periodic lattice distortion with three-fold rotational symmetry. In order to keep a close connection to existing materials, we also explore periodicities compatible with experimental findings54,55. We study a representative number of structures for each periodicity (in the 3 × 3 periodicity we study three known CDW structures30,55,56,73 anew, complementing previous works by showing the energy hierarchy upon strain); here we refer to them as hexagonal (HX), chalcogen-centred triangular (CC) and hollow centred triangular (HC). Energy differences, shown in Table 1, are used to test whether and how each periodicity can give rise to a CDW structure.

Table 1 Total energy of various CDW structures, grouped per periodicity, at three different values of the in-plane lattice constant.
Full size table

All 2 × 2 structures tested in our calculations relaxed to the symmetric structure, for all values of the strain; therefore, the 2 × 2 periodicity in ref. 54 does not arise from an intrinsic character of a single layer. The 2 × 3 periodicity accounts for the occurrence of instabilities at 2/3 ΓM and at M; compare with ref. 55 (Fig. 3), where HC and CC merge. The energy hierarchy between CDW structures with 3 × 3 periodicity is in agreement with existing ab initio studies30,56. The \(\sqrt{13}\times \sqrt{13}\) and the 4 × 4 periodicities allow the relaxation of three structures each, energetically close to the known 3 × 3 CDWs. In particular, the 4 × 4 CDW structures are favoured over the 3 × 3 CDW ground state under compressive strain. In general, a plethora of structural modulations are available in single layer NbSe2, whose competition depends and may be tuned by epitaxial strain, charge transfer, or proximity effects.

All starting structures for all the supercells considered here have been given a symmetry compatible with a triple ordering vector for the lattice distortion. The relaxed structures are found in accordance; however, we find a surprising result in one supercell. A total of three relaxed structures are identified within the 4 × 4 periodicity; two of them are characterised by a single ordering vector and referred to as 1q′ and 1q″ (Fig. 2a, c, respectively); a third one, referred to as 3q, is characterised by three equivalent ordering vectors (see Fig. 2b); it is energetically unfavoured at the equilibrium lattice constant and degenerates into the 1q′ (1q″) for compressive (tensile) strain. The wavelength of the 1q′ and the 1q″ charge distributions is 11.82 Å (3.43a0, where a0 is the lattice constant of the unstrained unit cell), 11.96 Å (3.47a0) and 12.09 Å (3.50a0), for compressive strain, no strain and tensile strain, respectively, for both phases. These results show a strong connection between the phases shown in Fig. 2 and the stripe phase observed at the surface53,54. In fact, two features of the modelled 4 × 4 CDWs are shown to be in remarkable agreement with experiments: the Fourier transform of the charge distribution—compare the inset of Fig. 2 with that resulting from STM in Fig. 1F in ref. 54—and the wavelength of the modulation—compare our result with the measured wavelength from STM data, ~12 Å53, 3.5a0. Further comparison to these experimental features is provided by our simulated STM plots of the 1q′ phase, shown in Fig. 3, left column. Additional STM plots, as well as a more extensive discussion, are reported in the Supplementary Information. The previous analysis of the stripe phase is based on the comparison between our calculations for the single layer (1H) and experimental results for the surface (2H). As pointed out in the “Introduction” section, these systems have different characteristics and therefore a certain care is needed when comparing them. To support our conclusions, we also performed selected calculations for the bilayer, whose symmetries and physical properties are a better model of the surface. Relevant structures with 3 × 3 and 4 × 4 periodicities were simulated, and the latter were found to be more favourable of 0.3 eV with respect to the former (compare with Table 1). Further details on these calculations are presented in the Supplementary Information. Finally, it is also important to comment on the fact that the stripe phase remains unobserved in single layers (1H) (see e.g. ref. 37). This is likely due to unfavourable or insufficient epitaxial strain of the substrates where the single layers are deposited, whereas applied54 and/or accidental strain53 favour such transition. A recent paper by Kim et al.74 discussed this point in terms of dielectric screening (comparing two different substrates), elucidating the role of spin–orbit coupling in the electronic properties and in the CDW formation. The emergence of a stripe phase under isotropic-simulated strain reveals an intrinsic tendency for the CDW towards a lower symmetry, which manifests in the 4 × 4 periodicity and remains latent in the 3 × 3 periodicity. In a recent paper57, calculations based on a Ginzburg–Landau formalism show that phonon fluctuations suppress long-range order in favour of a short-ranged pseudogap phase with a 3q phase, which degenerates into a 1q phase upon uniaxial strain. The authors restrained their study to a 3 × 3 CDW. Within 3 × 3 superlattices, with (biaxial) isotropic strain, we cannot find a stripe phase either, and even when a stripe structure is given as an input, it degenerates into a triangular CDW (CC or HC, depending on the starting input). All in all, isotropic strain induces only an incipient stripe phase with a 3 × 3 periodicity (see Fig. 4). The HC and CC CDW patches are distinct under tensile strain (Fig. 4c, f) and at the equilibrium lattice constant (Fig. 4b, e), compare the three-spikes ring in the CC CDW with the three-fold fidget-spinner-like star in the HC CDW, but the HC degenerates into the CC for compressive strain (Fig. 4a, d), where the rings and stars patches are merged. On the other hand, allowing the 4 × 4 periodicity, in agreement with the imaginary phonon dispersion, enhances the degrees of freedom and widens the ordering parameters space. Such periodicity hosts both single q and triple q phases, and therefore is suitable to analyse a relation between the electronic reconstruction and the symmetry of the ordering vectors. We consider the electron–phonon coupling as the driving mechanism for the CDW, as widely reported in the literature. In a Peierls instability scenario, where the electrons dominate over the phonons in driving the transition, a CDW results in a spectral weight loss at the Fermi level, which is also consistent with a weak-coupling. It has been already shown16,19 that it is not the case of NbSe2, where instead the electronic spectral weight shifts are more important at a lower energy range, implying that the transition is governed by phonons (note also the presence of electronic states below the Nb-derived band and its dependence on the CDW, in ref. 19). Figure 5 shows the density of the electronic states (DOS) of the single-vector CDWs (denoted both as 1q because their DOS curves have no substantial difference), the 3q CDW structure and the undistorted structure. A steep curve, showing increase (depletion) of states below (above) the Fermi level, characterises the 3q CDW, whereas the 1q CDW features a depletion of spectral weight all around the Fermi level. The 3q CDW displays also a clear and narrow dip at −0.08 eV. A more detailed analyis of the spectral properties, complemented with unfolded band structures, is reported in the Supplementary Information.

Fig. 2: Emerging structures in single layer NbSe2 under strain.
figure 2

Structures a–c and corresponding charge d–f modulations resulting from a 4 × 4 periodicity. The wavelength of the stripe phases 1q′ a and d and 1q″ c and f is 11.8 Å; the central column b and e refers to the 3q phase. Dark green (grey) spheres represent Nb (Se) atoms; Nb–Nb bonds shorter than the equilibrium distance (the respective lattice constant) are represented by solid lines, in order to help visualising the CDW structure pattern. Dashed lines mark the supercells borders. The inset in d shows the Fourier transform of the charge density of the 1q′ phase as computed in a \(4\times 4\sqrt{3}\) rectangular cell; the black (green) circle denotes one of the Bragg (CDW) peaks.

Full size image
Fig. 3: Simulated STM maps in the emergent CDW.
figure 3

STM maps for 1q′ and 3q 4 × 4 CDWs, left and right columns, respectively; rows 1, 2 and 3 stand for 0.00, −0.20 and −1.25 V, respectively.

Full size image
Fig. 4: A missed stripe character in the 3 × 3 CDW.
figure 4

Charge distribution of the CC a–c and the HC d–f CDW structures under compressive strain a, d, no strain b, e and tensile strain c, f. At the bottom, a sketch shows how the charge patches in the HC and CC merge in the HC CDW under compressive strain, as appears in d. Each pair of green, yellow and red bars connect two structures with equivalent, similar and different charge patches, respectively. At the bottom, a sketch shows what happens in figure d, as two different structures merge under compressive strain.

Full size image
Fig. 5: Density of states of 4 × 4 CDW in comparison.
figure 5

Density of the electronic states (DOS) of the non-modulated structure (in black, denoted by “sym”) and two representative CDW structures with a 4 × 4 periodicity, 1q and 3q (DOS curves for CDW structures with a single ordering vector do not differ sensibly). All curves are rescaled in formula units for a meaningful comparison between different supercells.

Full size image

Charge density distributions of the 3q and 1q CDWs in direct space are reported in Figs. 6 and 7, respectively. These data show a phase change of the electronic modulation across −0.10 eV, in analogy with ref. 23. Starting from the Fermi level (see Fig. 6e), the 3q CDW has a 3q symmetry which is maintained down to −0.25 eV (see Fig. 6c); below this energy, the symmetry is reduced to a single q. On the other hand, the 1q CDW survives down to −2.6 eV (see Fig. 7), pointing to an even stronger phononic involvement in the transition for this phase. The lattice distortions induce variations of the DOS at several energy levels, see Fig. S1 in the Supplementary Information, with reference to the 3 × 3 CDW ground state, and compare it with ref. 19; around −2.6 and −1.2 eV large variations occur. Note that the symmetry of the charge distribution has particular features depending on the energy range of the electrons. As CDW-induced spectral weight shifts may differ in intensity at different energy ranges, STS maps may differ from STM maps; for example, Fig. 3 shows that the 3q phase still exhibits a three-fold symmetry for all indicated bias voltages in constant-current STM images, where the electronic states are integrated within a corresponding energy window (from the Fermi level to the bias voltage).

Fig. 6: Integrated charge density for the triple q CDW.
figure 6

Charge distribution of the 3q CDW integrated over the following energy range"s: (−2.65, −2.55) eV a, (−0.45, −0.35) eV b, (−0.25, −0.15) eV c, (−0.15, −0.05) eV d, (−0.05, +0.05) eV e, (+0.05, + 0.15) eV f.

Full size image
Fig. 7: Integrated charge density of the single q CDW.
figure 7

Charge distribution of the 1q′ CDW integrated over the following energy ranges: (−2.65, −2.55) eV a, (−0.45, −0.35) eV b, (−0.25, − 0.15) eV c, (−0.15, −0.05) eV d, (−0.05, +0.05) eV e, (+0.05, +0.15) eV f.

Full size image

Conclusions

In conclusion, we analysed the lattice distortions and their accompanying charge modulations in single layer 1H-NbSe2 under biaxial isotropic strain, finding that a compressive strain of ~ 1% favours 4 × 4 modulations with a single q ordering vector. Together with two phases characterised by a single ordering vector, a metastable triple ordering vector phase exists. On the other hand, 3 × 3 modulations allow only an incipient 3q–1q transition under strain, but the charge patches of the ground state CDW structure assume features from the next most favourable structure. Although the strain in our model may not have been applicable by suitable (e.g. chemically non-reactive) substrates so far, the present study offers insights on the intrinsic properties of the single layer NbSe2, which can be useful for future experimental investigations on the manipulation of collective excitations in NbSe2 and related transition metal dichalcogenides. Moreover, we observe that the stripe phase observed in NbSe2 thin films is very likely associated to 4 × 4 modulations, and favoured by compressive strain.