Strain-induced stripe phase in charge ordered single layer NbSe$_2$

Charge density waves are ubiquitous phenomena in metallic transition metal dichalcogenides. In NbSe$_2$, a triangular $3\times3$ structural modulation is coupled to a charge modulation. Recent experiments reported evidence for a triangular-stripe transition at the surface, due to strain or accidental doping and associated to a $4\times4$ modulation. We employ \textit{ab-initio} calculations to investigate the strain-induced structural instabilities in a pristine single layer and analyse the energy hierarchy of the structural and charge modulations. Our results support the observation of phase separation between triangular and stripe phases in 1H-NbSe$_2$, relating the stripe phase to compressive isotropic strain, favouring the $4\times4$ modulation. The observed wavelength of the charge modulation is also reproduced with good accuracy.

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 P 6 3 /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 D 3h 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 freedom [32][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 superconductivity [37,42]. In the case of 2H-NbSe 2 , the CDW is recognised to arise from a phonon instability [43], inducing an electronic reconstruction in a sizable energy range [23] in the bulk and at the surface, although single layers show a CDW gap centred at the Fermi level [37]. Electronic and magnetic texture may be layer-resolved [44,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 bulk [46], may be synthesised [47] or induced by external sources [48]; the resulting charge order is not the usual star of David phase, unlike TaSe 2 [49,50] and TaS 2 [51,52]. In the 2H-type, 4 × 4 [53,54] and/or 2 × 2 [54] modulations, linked to a stripe phase at the surface, are reported. Distinct CDW structures, supporting the widely known CDW incommensurate character in 2H-NbSe 2 [55] 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 NbSe 2 [20,56], but an early scanning tunnelling microscopy/spectroscopy (STM/STS) study [37] confirms that the 3 × 3 modulation dominates in single layers deposited on bilayer graphene; on the other hand, a recent theoretical study [57] showed that within a 3 × 3 modulation, uniaxial strain leads to a triangular-stripe CDW transition in 2H-NbSe 2 . 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-NbSe 2 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.

II. 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) pseudopotentials [58,59], as implemented in the QUANTUM ESPRESSO suite [60,61] and the Vienna Abinitio 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 parametrisation [62,63]. The pseudopotentials used in the simulations provide explicit treatment of the following valence electrons: 4s 2 4p 6 5s 1 5d 4 for Nb and 4s 2 4p 4 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. 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 Zunger [64]. STM simulations were also performed, by means of the BSKAN code [65,66]. The revised Chen method [67] with an electronically flat and spatially spherical tip orbital has been used, which is equivalent to the Tersoff-Hamann model of electron tunnelling [68]. 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 Appendix.

III. 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 inplane stress is 7.2 × 10 −3 eVÅ −3 , 1.1 × 10 −3 eVÅ −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 study [19] reports a lattice constant of 3.458Å, in agreement with the residual strain found for |a| = 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 √ n × √ n [70] and n × m (suggested by discussion of results from STM [55]).
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, (a) to (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 √ 13 × √ 13 modulation is expected from the lowest values of the phonon dispersion, see especially Fig. 1 (b), 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 early [71] and recent [72] STM observations. In particular, a 2H-1T structural phase separation is re-  The energy of each structure is given as the difference, expressed in meV/f.u., with the energy of the symmetric structure, and it is negative for favoured structures. The names of the structures are explained as follows: MM and MW in the 2 × 3 cells are two stripe phases with a different Nb-Nb pattern; HX, CC and HC in the 3 × 3 cells are the known hexagonal, triangular chalcogen centred and triangular hollow centred structures, respectively; CCccws, HXHC and HC in the √ 13 × √ 13 cells are phases with triangular chalcogen centred plus counterclockwise vortex, hexagonal with triangular chalcogen centred merged and just triangular chalcogen centred, respectively; finally, the 1q , 3q and 1q in the 4 × 4 cells are a single ordering vector CDW structure, a triple ordering vector CDW structure and a single ordering vector structure, respectively, as also illustrated in Fig. 2 Further, we note similarities between the phonon dispersion in the 1H-type and that in the 1T-type [70]. 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).
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 findings [54,55]. We study a representative number of structures for each periodicity (in the 3 × 3 periodicity we study three known CDW structures [30,55,56,73] anew, complementing previous works by showing the energy hierarchy upon strain); here we refer to them as hexagonal (HX), chalcogen centered triangular (CC) and hollow centered triangular (HC). Energy differences, shown in Table I, are used to test whether and how each periodicity can give rise to a CDW structure.
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, figure 3, where HC and CC merge. The energy hierarchy between CDW structures with 3 × 3 periodicity is in agreement with existing ab-initio studies [30,56]. The √ 13× √ 13 and the 4×4 periodicities allow the relaxation of three structures each, energetically close to the known 3×3 general, a plethora of structural modulations are available in single layer NbSe 2 , 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 (Figs. 2 (a) and (c), respectively); a third one, referred to as 3q, is characterised by three equivalent ordering vectors, see Fig. 2 (b); 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.43 × a 0 , where a 0 is the lattice constant of the unstrained unit cell), 11.96Å (3.47 × a 0 ) and 12.09Å (3.50 × a 0 ), 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 surface [53,54]. 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, 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 I). Further details on these calculations are presented in the Appendix. 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 applied [54] and/or accidental strain [53] favour such transition. A recent paper by Kim et al. [74] discussed this point in terms of dielectric screening (comparing two different sub-  strates), 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 paper [57], 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 degen- erates 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 (c,f) and at the equilibrium lattice constant (b,e), compare the three-spikes ring in the CC CDW with the three-fold fidgetspinner-like star in the HC CDW, but the HC degenerates into the CC for compressive strain (a,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 shown [16,19] that it is not the case of NbSe 2 , 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 Appendix.
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. 6 (e), the 3q CDW has a 3q symmetry which is maintained down to −0.25 eV, see Fig. 6 (c); 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. ?? in the Appendix, with reference to the 3 × 3 CDW ground state, and compare it with ref. 19; around −2.6 eV 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).

IV. CONCLUSIONS
In conclusion, we analysed the lattice distortions and their accompanying charge modulations in single layer 1H-NbSe 2 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 NbSe 2 , which can be useful for future experimental investigations on the manipulation of collective excitations in NbSe 2 and related transition metal dichalcogenides. Moreover, we observe that the stripe phase observed in NbSe 2 thin films is very likely associated to 4 × 4 modulations, and favoured by compressive strain. In the main text, we presented the results for integrated charge density outside the energy range shown in Fig. 5 (main text). For completeness, we report the DOS of the 3 × 3 HC CDW and the undistorted structure over an extended energy range in Fig. S1. This plot provides evidence for the presence of electronic states at various ranges and for the regions where major CDW-driven changes occur. We also analyse the phase shifts from Fig. 6 (main text) in more details. As discussed in the main text, the symmetry of this CDW is a 3q around the Fermi level, and down to −0.25 eV. Within the Nb-derived band, although maintaining a triple q ordering vector, a π inversion - Fig. 6 (c)-(d) (main text) -and a rigid shift of the charge patches - Fig. 6 (e)-(f) (main text) -occur; compare this trend with recent measurements [S1]. Conversely, the 1q phases keep the same symmetry, see Fig. 7 (a)-(f) (main text). In the next section, we illustrate the STM maps at different tunnelling voltages. Furthermore, we compare the behaviour of the 3q CDW in the 4 × 4 superlattice with the CDWs in the 3 × 3 superlattice. Further consideration is deserved by the strain-induced variations in the DOS, illustrated in Fig. S2. Important changes at the Fermi level are seen especially in the 4 × 4 CDWs, where tensile (compressive) strain tends to increase (decrease) the spectral weight.

VII. STM/STS
Simulations of STM were carried out by employing the BSKAN code [S2, S3]. The revised Chen method [S4] with an electronically flat and spatially spherical tip orbital is used, which is equivalent to the Tersoff-Hamann model of electron tunnelling [S5]. The reported STM images are in constant current mode (with the maxima of the current contours at 5.8 Angstroms for all structures for a better comparability).
The STM maps of the 3 × 3 and 4 × 4 CDWs are shown in Figs. S3 and 3 (main text), respectively. For the 3 × 3 CDWs, we find a good agreement with existing experimental STM data [S1, S6-S11]. First, we point out that it is well known that the tip-sample distance can have a major influence on the observed STM contrasts, and a complete contrast reversal is also possible, see e.g. ref. S12. This effect is crucial for CDWs in hetero-layered structures such as transition metal dichalcogenides as well, where the more charged regions can be imaged as protrusions or depressions depending on the tip and tunneling conditions. In our case, modulations can result from a superposition of contribution from the transition metal layer and the chalcogen layer. As illustrated in figures 9 and 1 of ref. [S13], the Se-Se bond patterns of CC and HX have the same symmetry as the Nb-Nb bond patterns of HX and CC, respectively. Therefore, analysing the STM maps and the charge distribution, some care has to be taken. As our set of examples show, the observed STM contrasts may depend dramatically on the bias voltage and on the tip-sample distance (not shown), e.g., all observed STM contrasts are inverted at -1.25 V with respect to the other considered voltages in both Figs. S3 and 3 (main text). In some cases, the contrast, see e.g. Fig. S3, is completely inverted with respect to what it is expected from the visualised isosurface charge, shown e.g. in Fig. 4 (main text) or in previous ab-initio studies [S13]. The reasons for this puzzling trend are the different electronic states contributing to the tunneling current at different bias voltages, and the variation of the decay of the electron wave  functions of different symmetries in the vacuum. Nevertheless, a good agreement with a measured phase shift [S1] is found. The STM maps for two of the 4 × 4 CDWs, 1q and 3q are shown in Fig. 3 (main text). As opposed to the integrated charge density around specific energy values, where the symmetry of the 3q is not maintained below a certain energy, these STM maps show that the three-fold symmetry is maintained for integrated states from the Fermi level below −0.4 V, pointing to a stronger influence of the Nb-derived band in determining the symmetry of the charge modulation. A phase shift of the 1q phase, more clear than in the charge plots shown in Fig. 7 (main text) is observed; the same phase shows a clearer stripe character as a consistent voltage is considered. Regarding the comparison between periodicities, the STM maps for the 3q resemble some characteristic feature of the 3 × 3 CDWs.

VIII. BAND STRUCTURES
In view of further exploration of the stripe phase and the 4 × 4 vs 3 × 3 periodicity, we analyse the band structures for the 3 × 3 and 4 × 4 CDWs, with reference to Fig. S5, comparing it with the band structure of a unit cell, shown in   S4. The calculated band structure of the supercells were unfolded back to the primitive cell Brillouin zone to get an effective band structure of the systems [S14]. This is done following the method of Popescu and Zunger [S15], which is based on calculating the spectral weights of the supercell eigenvalues and constructing the spectral function. The orange lines in Figs. S5 correspond to the eigenvalues within the Brillouin zone of the supercell (folded), whereas the green curves represent the eigenvalues projected onto the unfolded Brillouin zone with a broadening proportional to their spectral weight. The Nb-derived band crossing the Fermi level features band splittings around M and along KΓ, but the latter is much less pronounced for the 4 × 4 CDWs. The valence band below features a splitting along KΓ, MK and ΓM, accompanied to a plethora of valence bands arising from the lattice distortions. Compared to angle resolved photoemission spectroscopic measurements [S7, S11, S16, S17], these bands look more detailed. The conduction band splitting along ΓM and MK are in agreement; the lower valence band features a two-fold bifurcation at ∼ −1 eV along ΓM and especially ΓK; recent experiments [S11] detect such feature along ΓK only, as a result of spin-orbit interaction -which in our calculations is not taken into account.

IX. A HINT TO THE PROPERTIES OF A BILAYER
Throughout this manuscript, we presented a model for a stripe phase in single layer 1H-NbSe 2 . Although the stripe phase was not observed for single layers, it was observed for thin films at the surface. As mentioned in the main text, the strain a substrate can apply on a single layer may not be sufficient, or it may be opposite in sign, for the stripe phase to occur. In order to strengthen our argument, we computed the phononic spectra of a bilayer 2H-NbSe 2 , using a k-mesh of 40 × 40 × 2 for the wavefunctions and a q-mesh of 8 × 8 × 1 for the force constants. Due to the centrosymmetry of the 2Hbilayer, instabilities occur along the ΓM line (around 2/3). We modelled bilayers with 3 × 3 and 4 × 4 periodicities, taking few examples from the single layer case, finding that the energy difference between 4 × 4 and 3 × 3 lattice distortions (0.3 meV per unit cell in favour of the 4 × 4 periodicity) is in line with the above reported results for the single layer. Although we have not tested all possible configurations, these results suggest that the ordering between 3 × 3 and 4 × 4 is not highly altered by the symmetry; however, defining a clear trend is not within the scope of the present work. The relaxed structures for the 3×3 periodicity maintain the same geometry they have in the single layer, whereas the 3q interfaced with a 1q relaxes to a 1q structure; note that the ordering vectors of the two interfaced layers are twisted by π/3 with respect to each other, see Fig. S6. The 3 × 3 bilayers are not shown because there is no new information about the structure with respect to the single layers.