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

TL;DR

This study uses ab-initio calculations to explore how strain influences phase transitions between triangular and stripe charge density wave patterns in single-layer NbSe2, revealing strain-induced phase separation and modulation characteristics.

Contribution

It provides a detailed theoretical analysis of strain-induced structural instabilities and phase separation in pristine single-layer NbSe2, linking strain to specific charge modulation patterns.

Findings

Stripe phase is favored by compressive isotropic strain.

The wavelength of charge modulation matches experimental observations.

Phase separation occurs between triangular and stripe phases.

Abstract

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.

Tables

Table 1. Table 1: Total energy of various CDW structures, grouped per periodicity, at three different values of the in-plane lattice constant. The energy of each structure is given as the difference, expressed in meV millielectronvolt \mathrm{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 2 3 2\times 3 cells are two stripe phases with a different Nb-Nb pattern; HX, CC and HC in the 3 × 3 3 3 3\times 3 cells are the known hexagonal, triangular chalcogen centred and triangular hollow centred structures, respectively; CCccws, HXHC and HC in the 13 × 13 13 13 \sqrt{13}\times\sqrt{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 1 ​ 𝐪 ′ 1 superscript 𝐪 ′ 1{\mathbf{q}}^{\prime} , 3 ​ 𝐪 3 𝐪 3{\mathbf{q}} and 1 ​ 𝐪 ′′ 1 superscript 𝐪 ′′ 1{\mathbf{q}}^{\prime\prime} in the 4 × 4 4 4 4\times 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 .

	$2\times 3$		$3\times 3$			$\sqrt{13}\times \sqrt{13}$			$4\times 4$
	MM	MW	HX	CC	HC	CCccws	HXHC	HC	$1{𝐪}^{\prime }$	$3𝐪$	$1{𝐪}^{\prime \prime }$
$\left|𝐚\right|$ = 3.49	$-1.9$	$-0.9$	$-3.0$	$-3.5$	$-4.4$	$-0.3$	$-0.4$	$-0.4$	$-2.9$	$--$	$-3.0$
$\left|𝐚\right|$ = 3.45	$-1.9$	$-0.0$	$-2.8$	$-3.4$	$-4.0$	$-1.4$	$-0.6$	$-1.4$	$-4.4$	$-2.5$	$-4.3$
$\left|𝐚\right|$ = 3.41	$-2.3$	$-0.0$	$-2.8$	$-3.6$	$-3.6$	$-3.4$	$-2.4$	$-3.6$	$-6.2$	$--$	$-6.2$

Full text

I Introduction

Transition metal dichalcogenides are among the most exciting hosts for ordered phases, such as superconductivity and charge/spin modulations Saito et al. (2016); Wilson et al. (1974, 2001); Moncton et al. (1975); Revolinsky et al. (1965) because of their cooperative interactions Kiss et al. (2007). 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 works Rice and Scott (1975); Liu et al. (1998); Straub et al. (1999), 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) Wilson et al. (1974); Moncton et al. (1975); Revolinsky et al. (1965); Rahn et al. (2012); Castro Neto (2001); Cho et al. (2018). Among the proposed explanations, the electronically driven mechanism was ruled out Inosov et al. (2008); Rossnagel et al. (2001); Rossnagel and Smith (2007); Johannes et al. (2006); Johannes and Mazin (2008); Zhu et al. (2015); Ángel Silva-Guillén et al. (2016) in favor of a momentum-assisted mechanism Calandra et al. (2009); Valla et al. (2004); Rahn et al. (2012); Zheng et al. (2018); Arguello et al. (2014); Flicker and van Wezel (2015a); Zhu et al. (2015), which is likely influenced by external conditions, such as applied fields Galvis et al. (2018); Cho et al. (2018) or epitaxially induced strain Fu et al. (2017); Wei et al. (2017), and chemical Di Salvo et al. (1975); Chatterjee et al. (2015); Cossu et al. (2018) or gate Wei et al. (2017); Shao et al. (2016) 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 $P6_{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 Kotov et al. (2012); Guinea et al. (2014); Mak et al. (2012); Frindt (1972); Xi et al. (2015); Ugeda et al. (2016); Butler et al. (2013); Geim and Grigorieva (2013); Novoselov et al. (2016); Ge and Liu (2012). 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 Menard et al. (2015); Ugeda et al. (2016). In the case of 2H-NbSe2, the CDW is recognised to arise from a phonon instability Weber et al. (2011), inducing an electronic reconstruction in a sizable energy range Arguello et al. (2014) in the bulk and at the surface, although single layers show a CDW gap centred at the Fermi level Ugeda et al. (2016). Electronic and magnetic texture may be layer-resolved Riley et al. (2014); Bawden et al. (2016), 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 Kadijk and Jellinek (1971), may be synthesised Nakata et al. (2016) or induced by external sources Bischoff et al. (2017); the resulting charge order is not the usual star of David phase, unlike TaSe2 Liu et al. (2016); Börner et al. (2018) and TaS2 Law and Lee (2017); Fazekas and Tosatti (1980). In the 2H-type, $4\times 4$ Soumyanarayanan et al. (2013); Gao et al. (2018) and/or $2\times 2$ Gao et al. (2018) modulations, linked to a stripe phase at the surface, are reported. Distinct CDW structures, supporting the widely known CDW incommensurate character in 2H-NbSe2 Gye et al. (2019) are observed, further speculating on their topological connection. A change from the $3\times 3$ to the $4\times 4$ modulation was previously predicted by ab-initio studies on 2H/1H-type NbSe2 Calandra et al. (2009); Lian et al. (2018), but an early scanning tunnelling microscopy/spectroscopy (STM/STS) study Ugeda et al. (2016) confirms that the $3\times 3$ modulation dominates in single layers deposited on bilayer graphene; on the other hand, a recent theoretical study Flicker and van Wezel (2015b) showed that within a $3\times 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 $3{\mathbf{q}}$-$1{\mathbf{q}}$ (triangular-stripe) transition occurs in the $4\times 4$ modulation, whereas it is incipient in the $3\times 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 Blöchl (1994); Kresse and Joubert (1999), as implemented in the Quantum ESPRESSO suite Giannozzi et al. (2009, 2017) 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 parametrisation Perdew et al. (1996, 1997). 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 $\sim$ $1142\text{\,}\mathrm{eV}$ and that for the supercell calculations is $500\text{\,}\mathrm{eV}$; prior to the phonon calculations, an energy cutoff is also applied for the charge density and potential, of $\sim$ $4354\text{\,}\mathrm{eV}$. The energy tolerance on the electronic loops for the calculations in the unit cell is $1.4\text{\times}{10}^{-9}\text{\,}\mathrm{eV}$, whereas it is ${10}^{-6}\text{\,}\mathrm{eV}$ for the calculations in the supercells (${10}^{-7}\text{\,}\mathrm{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}\text{\,}\mathrm{eV}$ per unit cell and a residual force of $\sim$ ${10}^{-3}\text{\,}\mathrm{eV}\text{\,}{\mathrm{\SIUnitSymbolAngstrom}}^{-1}$ on each atom. The calculations of phonons are performed sampling the first Brillouin Zone with a $45\times 45\times 1$ k-grid (for the wavefunctions) and a $15\times 15\times 1$ q-grid (for the force constants). The formation of CDWs is analysed sampling the first Brillouin Zone by means of the $30\times 30\times 1$, $23\times 23\times 1$, $20\times 20\times 1$, $17\times 17\times 1$ and $15\times 15\times 1$ grids of ${\mathbf{k}}$-points for the $2\times 1$, $\sqrt{7}\times\sqrt{7}\times 1$, $3\times 3\times 1$, $\sqrt{13}\times\sqrt{13}\times 1$ and $4\times 4\times 1$ supercells, respectively, keeping approximately an equal distance between ${\mathbf{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 Popescu and Zunger (2012). STM simulations were also performed, by means of the BSKAN code Hofer (2003); Palotás and Hofer (2005). The revised Chen method Mándi and Palotás (2015) with an electronically flat and spatially spherical tip orbital has been used, which is equivalent to the Tersoff-Hamann model of electron tunnelling Tersoff and Hamann (1985). 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\text{\,}\mathrm{\SIUnitSymbolAngstrom}$, $3.45\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ and $3.49\text{\,}\mathrm{\SIUnitSymbolAngstrom}$, which aim at modelling the system under compressive strain ($\sim 1\%$), no strain and tensile strain ($\sim 1\%$), respectively. The precise value of the computed in-plane stress is $7.2\text{\times}{10}^{-3}\text{\,}\mathrm{eV}\text{\,}{\mathrm{\SIUnitSymbolAngstrom}}^{-3}$, $1.1\text{\times}{10}^{-3}\text{\,}\mathrm{eV}\text{\,}{\mathrm{\SIUnitSymbolAngstrom}}^{-3}$ and $-4.5\text{\times}{10}^{-3}\text{\,}\mathrm{eV}\text{\,}{\mathrm{\SIUnitSymbolAngstrom}}^{-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. Fang et al. (2018). Since the discrepancy is very small, this detail does not affect the conclusions of the present work. A recent ab-initio study Ángel Silva-Guillén et al. (2016) reports a lattice constant of $3.458\text{\,}\mathrm{\SIUnitSymbolAngstrom}$, in agreement with the residual strain found for $\left|{\mathbf{a}}\right|=$ $3.45\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ 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}$ Kamil et al. (2018) and $n\times m$ (suggested by discussion of results from STM Gye et al. (2019)).

The phonon dispersions are represented as two-dimensional plots, Fig. 1, with positive (negative) values of $(\omega({\mathbf{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 $\Gamma$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. 1 (b), where the troughs appear in couples at symmetric positions with respect to the $\Gamma$M high symmetry line. Each couple of instability represent two periodic lattice distortions with opposite chirality. This result supports early Komori et al. (1997) and recent Wang et al. (2009) STM observations. In particular, a 2H-1T structural phase separation is reported in refs. Wang et al., 2009 and Bischoff et al., 2017; Nakata et al., 2016. Further, we note similarities between the phonon dispersion in the 1H-type and that in the 1T-type Kamil et al. (2018). Nevertheless, these instabilities may couple to give rise to a $3\times 3$ modulation, which is the most observed. In a similar fashion, instability regions around M may couple to give rise to a $2\times 2$ (or even a $2\times 3$) modulation. Instabilities are expected to dominate along $\Gamma$M in a 2H-type symmetry (because the crystal is overall centrosymmetric).

Compressive strain reduces the instabilities along the MK line and the minima of $(\omega({\mathbf{q}}))^{2}$ dominate around 2/4 $\Gamma$M (off the high symmetry line), see Fig. 1 (a), pointing towards a $4\times 4$ modulation. Conversely, a zone of instability centered around, but not including, M becomes dominant under tensile strain, pointing towards a $2\times 2$ modulation; further local minima of $(\omega({\mathbf{q}}))^{2}$ appear at 2/4 $\Gamma$M, Fig. 1 (c), justifying recent reports Gao et al. (2018) for which both $2\times 2$ and $4\times 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 findings Gao et al. (2018); Gye et al. (2019). We study a representative number of structures for each periodicity (in the $3\times 3$ periodicity we study three known CDW structures Cossu et al. (2018); Gye et al. (2019); Malliakas and Kanatzidis (2013); Lian et al. (2018) 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 1, are used to test whether and how each periodicity can give rise to a CDW structure.

All $2\times 2$ structures tested in our calculations relaxed to the symmetric structure, for all values of the strain; therefore, the $2\times 2$ periodicity in ref. Gao et al., 2018 does not arise from an intrinsic character of a single layer. The $2\times 3$ periodicity accounts for the occurrence of instabilities at 2/3 $\Gamma$M and at M; compare with ref. Gye et al., 2019, figure 3, where HC and CC merge. The energy hierarchy between CDW structures with $3\times 3$ periodicity is in agreement with existing ab-initio studies Lian et al. (2018); Cossu et al. (2018). The $\sqrt{13}\times\sqrt{13}$ and the $4\times 4$ periodicities allow the relaxation of three structures each, energetically close to the known $3\times 3$ CDWs. In particular, the $4\times 4$ CDW structures are favoured over the $3\times 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\times 4$ periodicity; two of them are characterised by a single ordering vector and referred to as $1\mathbf{q}^{\prime}$ and $1\mathbf{q}^{\prime\prime}$ (Figs. 2 (a) and (c), respectively); a third one, referred to as $3\mathbf{q}$, is characterised by three equivalent ordering vectors, see Fig. 2 (b); it is energetically unfavoured at the equilibrium lattice constant and degenerates into the $1\mathbf{q}^{\prime}$ ($1\mathbf{q}^{\prime\prime}$) for compressive (tensile) strain. The wavelength of the $1\mathbf{q}^{\prime}$ and the $1\mathbf{q}^{\prime\prime}$ charge distributions is $11.82\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ ($3.43\times\mathbf{a}_{0}$, where $\mathbf{a}_{0}$ is the lattice constant of the unstrained unit cell), $11.96\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ ($3.47\times\mathbf{a}_{0}$) and $12.09\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ ($3.50\times\mathbf{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 Soumyanarayanan et al. (2013); Gao et al. (2018). In fact, two features of the modelled $4\times 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 figure 1F in ref. Gao et al., 2018 – and the wavelength of the modulation – compare our result with the measured wavelength from STM data, \sim$$12\text{\,}\mathrm{\SIUnitSymbolAngstrom} Soumyanarayanan et al. (2013), $3.5\times\mathbf{a}_{0}$. Further comparison to these experimental features is provided by our simulated STM plots of the $1{\mathbf{q}}^{\prime}$ phase, shown in Fig. 3, left column. Additional STM plots, as well as a more extensive discussion, are reported in the Appendix.

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\times 3$ and $4\times 4$ periodicities were simulated, and the latter were found to be more favourable of $0.3\text{\,}\mathrm{eV}$ with respect to the former (compare with Table 1). 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. Ugeda et al., 2016. This is likely due to unfavourable or insufficient epitaxial strain of the substrates where the single layers are deposited, whereas applied Gao et al. (2018) and/or accidental strain Soumyanarayanan et al. (2013) favour such transition. A recent paper by Kim et al. Kim and Son (2017) 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\times 4$ periodicity and remains latent in the $3\times 3$ periodicity. In a recent paper Flicker and van Wezel (2015b), 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 $3{\mathbf{q}}$ phase, which degenerates into a $1{\mathbf{q}}$ phase upon uniaxial strain. The authors restrained their study to a $3\times 3$ CDW. Within $3\times 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\times 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 fidget-spinner-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\times 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 ${\mathbf{q}}$ and triple ${\mathbf{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 Johannes et al. (2006); Ángel Silva-Guillén et al. (2016) 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. Ángel Silva-Guillén et al., 2016). Figure 5 shows the density of the electronic states (DOS) of the single-vector CDWs (denoted both as $1{\mathbf{q}}$ because their DOS curves have no substantial difference), the $3{\mathbf{q}}$ CDW structure and the undistorted structure. A steep curve, showing increase (depletion) of states below (above) the Fermi level, characterises the $3{\mathbf{q}}$ CDW, whereas the $1{\mathbf{q}}$ CDW features a depletion of spectral weight all around the Fermi level. The $3{\mathbf{q}}$ CDW displays also a clear and narrow dip at $-0.08\text{\,}\mathrm{eV}$. A more detailed analyis of the spectral properties, complemented with unfolded band structures, is reported in the Appendix.

Charge density distributions of the $3{\mathbf{q}}$ and $1{\mathbf{q}}$ 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\text{\,}\mathrm{eV}$, in analogy with ref. Arguello et al., 2014. Starting from the Fermi level, see Fig. 6 (e), the $3{\mathbf{q}}$ CDW has a $3{\mathbf{q}}$ symmetry which is maintained down to $-0.25\text{\,}\mathrm{eV}$, see Fig. 6 (c); below this energy, the symmetry is reduced to a single ${\mathbf{q}}$. On the other hand, the $1{\mathbf{q}}$ CDW survives down to $-2.6\text{\,}\mathrm{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. LABEL:SI-totdos-fig in the Appendix, with reference to the $3\times 3$ CDW ground state, and compare it with ref. Ángel Silva-Guillén et al., 2016; around $-2.6\text{\,}\mathrm{eV}$ and $-1.2\text{\,}\mathrm{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 $3{\mathbf{q}}$ 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-NbSe2 under biaxial isotropic strain, finding that a compressive strain of $\sim 1\%$ favours $4\times 4$ modulations with a single ${\mathbf{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\times 3$ modulations allow only an incipient $3{\mathbf{q}}$-$1{\mathbf{q}}$ 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\times 4$ modulations, and favoured by compressive strain.

V Acknowledgements

We are grateful to Yunkyu Bang, Poonam Kumari, Dhani Milind Nafday, Erio Tosatti, Stefano de Gironcoli, Ali G. Moghaddam, Tristan Cren, Hermann Suderow and Han-Woong Yeom for insightful discussions. F. C. and A. A. acknowledge financial support from the National Research Foundation (NRF) funded by the Ministry of Science of Korea (Grants No. 2016K1A4A01922028, No. 2017R1D1A1B03033465, & No. 2019R1H1A2039733); K. P. acknowledges support from NRDIO-Hungary project No. FK124100. This research was supported by appointments to the JRG program at the APCTP through the Science and Technology Promotion Fund and Lottery Fund of the Korean Government, and by the Korean Local Governments – Gyeongsangbuk-do province and Pohang City.

**Supplemental material for

Triangular-stripe transition in charge ordered single layer NbSe2**

by

Fabrizio Cossu, Krisztián Palotás, Sagar Sarkar, Igor Di Marco, and Alireza Akbari

VI Extended DOS and Phase shifts

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\times 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 $3{\mathbf{q}}$ around the Fermi level, and down to $-0.25\text{\,}\mathrm{eV}$. Within the Nb-derived band, although maintaining a triple ${\mathbf{q}}$ ordering vector, a $\pi$ 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 SI-arguello-PRB.89.235115 . Conversely, the $1{\mathbf{q}}$ 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 $3{\mathbf{q}}$ CDW in the $4\times 4$ superlattice with the CDWs in the $3\times 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\times 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 SI-hofer-PSS2003 ; SI-palotas-JPCM2005 . The revised Chen method SI-mandi-PRB.91.165406 with an electronically flat and spatially spherical tip orbital is used, which is equivalent to the Tersoff-Hamann model of electron tunnelling SI-tersoff-PRB.31.805 . 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\times 3$ and $4\times 4$ CDWs are shown in Figs. S3 and 3 (main text), respectively. For the $3\times 3$ CDWs, we find a good agreement with existing experimental STM data SI-soumyanarayanan-PNAS2013 ; SI-ugeda-nphys2016 ; SI-arguello-PRB.89.235115 ; SI-chatterjee-ncomm2015 ; SI-gye_gc-PRL.122.016403 ; SI-gao-PNAS2018 ; SI-nakata-NPJ2DMA2018 . 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. SI-palotas-PRB.86.235415, . 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. SI-cossu-PRB.98.195419 , 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 SI-cossu-PRB.98.195419 . 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 SI-arguello-PRB.89.235115 is found. The STM maps for two of the $4\times 4$ CDWs, $1{\mathbf{q}}^{\prime}$ and $3{\mathbf{q}}$ are shown in Fig. 3 (main text). As opposed to the integrated charge density around specific energy values, where the symmetry of the $3{\mathbf{q}}$ 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\text{\,}\mathrm{V}$, pointing to a stronger influence of the Nb-derived band in determining the symmetry of the charge modulation. A phase shift of the $1{\mathbf{q}}^{\prime}$ 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 $3{\mathbf{q}}$ resemble some characteristic feature of the $3\times 3$ CDWs.

VIII Band structures

In view of further exploration of the stripe phase and the $4\times 4$ vs $3\times 3$ periodicity, we analyse the band structures for the $3\times 3$ and $4\times 4$ CDWs, with reference to Fig. S5, comparing it with the band structure of a unit cell, shown in Fig. 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 SI-popescu-PRL.104.236403 . This is done following the method of Popescu and Zunger SI-popescu-PRB.85.085201 , 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$\Gamma$, but the latter is much less pronounced for the $4\times 4$ CDWs. The valence band below features a splitting along K$\Gamma$, MK and $\Gamma$M, accompanied to a plethora of valence bands arising from the lattice distortions. Compared to angle resolved photoemission spectroscopic measurements SI-bawden-ncomm2016 ; SI-nakata-NPGAM2016 ; SI-ugeda-nphys2016 ; SI-nakata-NPJ2DMA2018 , these bands look more detailed. The conduction band splitting along $\Gamma$M and MK are in agreement; the lower valence band features a two-fold bifurcation at $\sim$ $-1\text{\,}\mathrm{eV}$ along $\Gamma$M and especially $\Gamma$K; recent experiments SI-nakata-NPJ2DMA2018 detect such feature along $\Gamma$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-NbSe2. 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-NbSe2, using a k-mesh of $40\times 40\times 2$ for the wavefunctions and a q-mesh of $8\times 8\times 1$ for the force constants. Due to the centrosymmetry of the 2H-bilayer, instabilities occur along the $\Gamma$M line (around 2/3). We modelled bilayers with $3\times 3$ and $4\times 4$ periodicities, taking few examples from the single layer case, finding that the energy difference between $4\times 4$ and $3\times 3$ lattice distortions ($0.3\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ per unit cell in favour of the $4\times 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\times 3$ and $4\times 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\times 3$ periodicity maintain the same geometry they have in the single layer, whereas the $3{\mathbf{q}}$ interfaced with a $1{\mathbf{q}}$ relaxes to a $1{\mathbf{q}}$ structure; note that the ordering vectors of the two interfaced layers are twisted by $\pi/3$ with respect to each other, see Fig. S6. The $3\times 3$ bilayers are not shown because there is no new information about the structure with respect to the single layers.