Coexisting charge density wave and ferromagnetic instabilities in monolayer InSe

Recently fabricated InSe monolayers exhibit remarkable characteristics that indicate the potential of this material to host a number of many-body phenomena. Here, we consistently describe collective electronic effects in hole-doped InSe monolayers using advanced many-body techniques. To this end, we derive a realistic electronic-structure model from first principles that takes into account the most important characteristics of this material, including a flat band with prominent van Hove singularities in the electronic spectrum, strong electron-phonon coupling, and weakly-screened long-ranged Coulomb interactions. We calculate the temperature-dependent phase diagram as a function of band filling and observe that this system is in a regime with coexisting charge density wave and ferromagnetic instabilities that are driven by strong electronic Coulomb correlations. This regime can be achieved at realistic doping levels and high enough temperatures, and can be verified experimentally. We find that the electron-phonon interaction does not play a crucial role in these effects, effectively suppressing the local Coulomb interaction without changing the qualitative physical picture.

Recently fabricated InSe monolayers exhibit remarkable characteristics that indicate the potential of this material to host a number of many-body phenomena. Here, we consistently describe collective electronic effects in hole-doped InSe monolayers using advanced many-body techniques. To this end, we derive a realistic electronicstructure model from first principles that takes into account the most important characteristics of this material, including a flat band with prominent van Hove singularities in the electronic spectrum, strong electron-phonon coupling, and weakly-screened long-ranged Coulomb interactions. We calculate the temperature-dependent phase diagram as a function of band filling and observe that this system is in a regime with coexisting charge density wave and ferromagnetic instabilities that are driven by strong electronic Coulomb correlations. This regime can be achieved at realistic doping levels and high enough temperatures, and can be verified experimentally. We find that the electron-phonon interaction does not play a crucial role in these effects, effectively suppressing the local Coulomb interaction without changing the qualitative physical picture.
Two-dimensional (2D) group III-VI metal chalcogenides have recently attracted great interest because of their appealing characteristics. Among them are high charge carrier mobility, controllable energy gaps, excellent thermoelectric and optical properties, as well as excellent stability at ambient conditions [1][2][3][4][5][6]. The electronic structure of ultrathin InSe features flat regions in the valence band dispersion leading to prominent van Hove singularities (vHS) in the hole density of states (DOS) [7]. Importantly, this kind of electronic structure is only observed in the monolayer limit of these materials, as has been experimentally demonstrated by means of angular resolved photoemission spectroscopy [8,9]. If the vHS appears at the Fermi energy, it may result in numerous competing channels of instabilities such as magnetic, charge, or superconducting order with a very non-trivial interplay between them [10][11][12]. Scientific interest to flat-band materials has been triggered in recent years by the discovery of unconventional superconductivity and related exotic phenomena in twisted bilayer graphene [13][14][15]. Therefore, thin films group III-VI materials turn out to be prospective candidates for studying many-body correlation effects, which are closely related to the above-mentioned features in the electronic spectrum [16], circumventing the need for any twist engineering.
The theoretical description of many-body effects in the monolayer InSe is challenging and is limited to a few works mainly focusing on the electron-phonon coupling. In particular, it has been shown that hole states in this material undergo significant renormalization due to the interaction with acoustic phonons, which gives rise to the appearance of unusual temperature-dependent optical excitations [17]. Also, a strong electron-phonon interaction may result in a charge density wave (CDW) instability predicted recently [18]. At the same time, monolayer InSe is expected to possess strong electronic Coulomb correlations that have not been accurately studied yet. For instance, the presence of the Mexican-hat-like band in this material might favors a magnetic instability that could lead to the formation of a magnetically ordered state at low temperatures. Up to now the existence of a magnetic solution for InSe monolayer has been demonstrated only at the level of density functional theory (DFT) [19]. In addition, the weakly screened Coulomb interaction in 2D may result in a Coulomb-driven CDW instability. However, a systematic many-body consideration of these effects is still missing in the literature.
In this Letter, we address the problem of collective electronic effects in monolayer InSe. For this purpose, we derive a realistic model that considers both, long-range Coulomb interactions and the electron-phonon coupling. The introduced interacting electronic problem is further solved using an advanced many-body approach that explicitly takes into account non-local collective electronic fluctuations. In the regime of hole-doping, we find that charge ordering represents the main instability in the monolayer InSe. It is formed in a broad range of doping levels and corresponds to a commensurate CDW, which indicates that this instability is rather driven by strong electronic Coulomb correlations than by an electronphonon mechanism as discussed previously. Inside the CDW phase, we detect another collective effect that drives the system towards a ferromagnetic (FM) ordering. This instability is formed in close proximity to the vHS in the DOS. Finally, we observe that the electron-phonon coupling tends to suppress the FM ordering, enlarging the CDW phase. However, the presence of the electron-phonon coupling does not qualitatively affect the observed effects.
Model -InSe is a layered van der Waals material, where each layer consists of two vertically displaced In-Se honeycombs, giving rise to four Se-In-In-Se atomic planes [ Fig. 1(a)]. In the monolayer limit DFT calculations predict InSe to be a semiconductor with an indirect energy gap of ∼ 2 eV. The electronic dispersion shows a single well-separated valence band, which has the shape of a Mexican hat, as depicted in Fig. 1(b). This peculiarity is of great advantage for many-body considerations as it allows us to reduce the correlated subspace to a single band. To this end, we construct a tractable tight-binding model that accurately reproduces this highest valence band [ Fig. 1(b)]. The corresponding model is defined in terms of maximally-localized Wannier functions on an effective triangular lattice, as shown in Fig. 1(a). Each Wannier function is reminiscent of an In-In bonding orbital with some tails on Se atoms. The resulting single-band model Hamiltonian on a triangular lattice reads shows the bare and screened Coulomb interaction between the Wannier orbitals calculated as a function of the distance within the constrained random phase approximation (cRPA). The on-site screened Coulomb repulsion U = 1.78 eV greatly exceeds the bandwidth ≈ 1 eV, which usually indicates well-developed magnetic fluctuations in the system. As expected for a 2D material, the non-local Coulomb interaction V i j in monolayer InSe is weakly-screened and long-ranged [see blue line in Fig. 1(c)]. Moreover, the interaction V 01 = 1.04 eV between nearest-neighbor electronic densities n iσ = c † iσ c iσ is larger than the half of the on-site Coulomb repulsion, V 01 > U/2. This suggests that the considered system may have a tendency to form a CDW phase due to the competition between local and non-local Coulomb interactions. The full form of the long-range Coulomb interaction is presented in the Supplemental Material (SM) [20].
To estimate the phonon properties we utilize the constrained Density Functional Perturbation Theory (cDFPT) [21] at hole doping. We find that the effective electron-phonon coupling λ = 2 dω α 2 F(ω) ω is dominated by a rather sharp resonance at low phonon frequency, which we can approximate with a local phonon model, i.e α 2 F(ω) = N 0 g 2 δ(ω − ω ph ). Here, N 0 is the DOS at the Fermi level, ω ph = 8.5 meV the phonon energy, and g = 34.7 meV its coupling strength to the upmost valence band [20]. The strong local coupling of electrons to phonons renders the CDW formation even more favorable. Indeed, upon integrating out bosonic operators b ( †) that correspond to phonon degrees of freedom one gets an effective local frequencydependent attractive interaction U ph ω = 2g 2 ω ph ω 2 ph −ω 2 [22][23][24]. This interaction reduces the repulsive on-site Coulomb potential as U → U − U ph ω and thus enhances the effect of the non-local Coulomb interaction V i j . Method -An accurate theoretical investigation of manybody instabilities in InSe monolayer cannot be performed within conventional perturbative methods like the random phase approximation [25][26][27] or the GW approach [28][29][30]. Local correlation effects that are governed by such a large value of the local Coulomb interaction U should be addressed using at least the dynamical mean-field theory (DMFT) [31]. At the same time, spatial collective electronic fluctuations and the long-range Coulomb interaction cannot be taken into account in the framework of DMFT and require to use diagrammatic extensions of this theory [32]. In this work, we solve the considered many-body problem using the dual triply irreducible local expansion (D-TRILEX) method [33][34][35]. The D-TRILEX approach is one of the simplest consistent diagrammatic extensions of DMFT that allows one to account for leading collective electronic fluctuations on equal footing without any limitation on the range. Within this method local correlation effects are treated via the self-energy Σ imp ν and the polarization operator Π ς imp ω of an effective local impurity problem of DMFT. The corresponding impurity problem is solved numerically exactly using the open source CT-HYB solver [36,37] based on ALPS libraries [38]. The spatial fluctuations are considered in a partially-bosonized form of the renormalized charge (ς = c) and spin (ς = s) interactions W ς qω that enter the diagrammatic part of the self-energy Σ kν introduced beyond DMFT. Thus, the D-TRILEX method can be seen as the GW extension of the DMFT that, however, additionally considers exact local three-point vertex corrections in diagrams for the self-energy and the polarization operator [33,34]. Such vertices are of a crucial importance for a correct description of magnetic, optical and transport proper- ties of the system [35,[39][40][41][42][43][44][45][46][47][48][49]. The dressed Green's function G kν of the problem can be found using the standard Dyson equation G −1 kν = iν + µ − ε k − Σ kν written in momentum k and fermionic Matsubara frequency ν space. In this expression µ is the chemical potential, ε k is the electronic dispersion that can be obtained as a Fourier transform of the hopping amplitudes t i j , and Σ kν = Σ imp ν + Σ kν is the total self-energy. The renormalized interaction W ς qω can be found via the following Dyson ph ω + V q and U s qω = −U/2 are the bare interactions in the charge and spin channels, respectively [33,34].
qω is the total polarization operator of the problem, where Π ς qω is the diagrammatic contribution introduced in the D-TRILEX approach [33,34]. Charge and spin susceptibilities X ς qω can then be obtained straightforwardly as X ς −1 qω = Π ς −1 qω − U ς qω [35]. Collective electronic instabilities -One of the most remarkable features of the InSe monolayer is the presence of a Mexican-hat-like valence band in the electronic dispersion [8,9]. The top of this band exhibits flat regions that lead to a sharp vHS in the DOS. However, under normal conditions this valence band is fully filled, making the material an indirect semiconductor. To enhance correlation effects in the system we consider realistic hole dopings with the Fermi level close to the vHS. Practically, high concentration of holes of the order of 10 14 cm −2 can be achieved in 2D materials by means of electrostatic solid-or liquid-electrolyte gating [50] or by surface molecular doping [51]. First, we solve the many-body problem without taking into account the electron-phonon coupling in order to investigate purely Coulomb correlation effects. For detecting main instabilities in the system we perform single-shot D-TRILEX calcu- lations for the charge and spin susceptibilities X ς qω . This way the diagrammatic part of the polarization operator Π ς qω is obtained non-self-consistently in terms of DMFT Green's functions G D kν , which are dressed only by the local self-energies: ν . This form of the D-TRILEX susceptibility resembles the DMFT susceptibility [52][53][54] with a longitudinal dynamical vertex corrections [34]. In this case, divergences of charge and spin susceptibilities do not affect each other through the self-energy, which allows one to detect instabilities inside broken-symmetry phases. Figure 2 shows the obtained phase diagram for the InSe monolayer, where n is the filling of the considered valence band ( n = 2 in the fully filled band that corresponds to the undoped case). Phase boundaries indicate points in the temperature (T ) vs. doping space, where corresponding susceptibilities diverge. We find that the charge susceptibility diverges in a broad range of hole dopings n ≥ 1.29, and the corresponding phase boundary is independent of temperature. The left panel of Fig. 3 displays the momentum resolved charge susceptibility X c qω obtained at the zero frequency ω = 0 near the transition point (T = 300 K, n = 1.29). It shows that the corresponding Bragg peaks in the charge susceptibility appear at the K points of the Brillouin zone (BZ), which indicates the formation of a commensurate CDW ordering. In turn, the spin susceptibility remains finite at the CDW transition point (see SM [20]) and diverges only inside the CDW phase. The corresponding instability has a dome shape as depicted in Fig. 2 by a blue dashed line. Remarkably, we find that the top of the dome corresponds to the filling n = 1.70 at which the vHS appears exactly at the Fermi level. The momentum resolved spin susceptibility obtained close to the top of the dome (T = 375 K, n = 1.70) is shown in right panel of Fig. 3. It reveals a sharp Bragg peak at the Γ point of the BZ, which indicates the tendency of the system towards FM ordering.
The obtained phase diagram illustrates that the commensurate CDW represents the main instability in the InSe monolayer contrary to the DFT prediction [19]. The formation of this insulating phase is associated with strong long-range collective charge fluctuations that are expected to renormalize the electronic dispersion. Therefore, the development of the CDW ordering should also be reflected in single-particle observables such as the electronic DOS. In order to account for the feedback effects of long-range collective electronic fluctuations to single-particle quantities we perform self-consistent D-TRILEX calculations [33,34] and obtain the full lattice Green's function G kν of the considered many-body problem that takes into account both, local and non-local correlation effects. Further, we take the local part of the lattice Green's function G loc ν = 1 N k G kν and perform an analytical continuation from Matsubara frequency space ν to real energies using the stochastic optimization method [55] to get the DOS. The obtained result is compared to the one of the DMFT that does not take into account spatial correlation effects. Figure 4 shows the DOS calculated via both methods close to the CDW phase boundary (T = 300 K and n = 1.29). We find that the DOS obtained within DMFT (blue line) is similar to the one of DFT [ Fig. 1(b)] and shows a sharp peak in the vicinity of the Fermi energy, which reflects the presence of the vHS in the electronic spectrum. However, if one additionally accounts for the effect of spatial correlations via the D-TRILEX approach, one observes that at the transition point this peak turns into a pseudogap (red line). In this particular case, the pseudogap appears due to strong collective charge fluctuations that lead to an almost diverged charge susceptibility X c qω at the q = K point of the BZ close to a phase transition. As a consequence, the renormalized interaction W c qω that enters the self-energy also becomes nearly divergent, which causes the formation of a pseudogap according to the tendency towards electron-hole pairing. This mechanism is similar to the formation of the excitonic insulator state [56][57][58] with the only difference that in our case electrons and holes belong to the same band and that excitons in the condensate have non-zero momentum corresponding to the CDW wave vector. If the tendency towards electron-hole pairing results in long-range order, one gets a true gap whereas strong shortrange order without long-range order leads to a pseudogap in the electronic spectrum [59,60]. This result illustrates the importance of a self-consistent consideration of long-range collective electronic fluctuations for capturing the formation of the insulating CDW phase in monolayer InSe.
In order to investigate the effect of phonon degrees of freedom on the observed instabilities we repeat the same calculation in the presence of the electron-phonon coupling. We find that in this case the CDW phase boundary is shifted to smaller values of the filling n = 1.23. At the same time, the FM instability is pushed down to lower temperatures, but the top of the FM dome remains at the vHS filling n = 1.70 as in the absence of phonons. This result is consistent with the presented argument above that the electron-phonon coupling effectively reduces the on-site Coulomb potential, which consequently decreases the critical temperature for the magnetic instability. The observed shift of the CDW phase boundary can also be explained by the same argument. Indeed, the local Coulomb repulsion favors the single occupation of lattice sites. On the contrary, the non-local Coulomb interaction promotes CDW ordering, which upon reducing the local Coulomb interaction becomes energetically preferable. Remarkably, we find that the position of the Bragg peaks in the charge and spin susceptibilities calculated close to corresponding phase boundaries remain unchanged [20]. In addition, we do not find any additional Bragg peak in the charge susceptibility at twice the Fermi wave vector 2k F , which would indicate a charge instability due to the electron-phonon mechanism. This fact suggests that many-body effects in the monolayer InSe are driven by strong Coulomb correlations rather than by the electron-phonon mechanism as suggested previously [18].
Conclusions -We have systematically studied many-body effects in the hole-doped InSe monolayer. We have found that this material displays coexisting instabilities that are mainly driven by nonlocal Coulomb correlations. The commensurate CDW ordering represents the main instability in the system and is revealed in a broad range of doping levels and temperatures. This ordering is accompanied by the formation of an insulating phase that is confirmed by the appearance of a pseudogap in the DOS close to the transition point, which illustrates the importance of considering spatial electronic fluctuations. We also observed the tendency to a FM ordering that manifests itself only inside the CDW phase and is related to a vHS in the electronic spectrum. The inclusion of the electron-phonon coupling results in a shift of the CDW and FM ordering phases on the phase diagram, which can be explained by an effective reduction of the local Coulomb interaction. However, the qualitative physical picture does not change in the presence of phonon degrees of freedom. Our results suggest that monolayer InSe can serve as an attractive playground for investigation of coexisting many-body correlation effects and, in particular, of 2D magnetism, although in a bulk phase this material is non-magnetic.
We thank Andy J. Millis for fruitful discussions and Jan Berges for sharing his cDFPT-related changes to Quantum Espresso with us. The work of E.A.S. was supported by the European Union's Horizon 2020 Research

DETAILS OF FIRST-PRINCIPLES CALCULATIONS
The band structure presented in Fig. 1(b) was calculated within density functional theory utilizing the projected augmented wave (PAW) formalism [61,62] as implemented in the Vienna ab initio simulation package (vasp) [63,64]. The exchange-correlation effects were considered using the generalized gradient approximation (GGA) [65]. A 500 eV energy cut-off for the plane-waves and a convergence threshold of 10 −7 eV were used in the calculations. The Brillouin zone was sampled by a (36×36) k-point mesh. We adopted fully relaxed atomic structure with a lattice constant of 3.94 Å. A ∼30 Å-thick vacuum layer was added in the direction perpendicular to the 2D plane in order to avoid spurious interactions between supercell images. The Wannier functions and the tight-binding Hamiltonian were calculated within the scheme of maximal localization [66,67] using the wannier90 package [68].
The Coulomb interaction was evaluated using the maximally localized Wannier functions within the con- strained random phase approximation (cRPA) [69,70] as U i j = w i w j |U|w j w i , where U is the partially screened Coulomb interaction defined by U = v + vΠU with v being the bare Coulomb interaction, Π the cRPA polarization, and w i is the Wannier function at the lattice site i. The polarization operator Π describes screening from all electronic states except those given by the tight-binding Hamiltonian obtained in the Wannier basis. For these calculations, we used a recent cRPA implementation by Kaltak within vasp [70]. To converge the cRPA polarization with respect to the number of empty states with used in total 64 bands. To derive a lightweighted Coulomb model for arbitrary q grids, we fitted the cRPA Coulomb interaction in momentum space according to: with the bare Coulomb interaction of a monolayer v(q) = h 2π Momentum resolved charge susceptibility X c (q, ω = 0) calculated taking into account the electron-phonon coupling. Result is obtained close to the CDW transition point T = 300 K and n = 1.23. and the dielectric function with ε 0 (q) = a + q 2 a sin(qc) Here, e, A, and h are the elementary electron charge, the InSe unit cell area and its effective height, respectively, and a, b, c, and d are fitting parameter. In Fig. 5 we show the corresponding ab initio values together with their fits using: The real-space U i j are calculated via a coventional Fourier transform of U(q) from Eq. (1) and using these fits. The phonon properties are calculated within Quantum Espresso [71] using norm-conserving pseudopotentials, the local density approximation, an energy cutoff of 80 Ryd, and a finite hole doping of 0.04 holes per unit cell. For the intial constrained Density Functional Perturbation Theory (cDFPT) [21] calculation we use (16×16) k and (8×8) q-point meshes and exclude the highest (hole doped) valence band within the evaluation of the Sternheimer equation. Afterwards we use the EPW code [72] to extrapolate the cDFPT results to (32×32) kand q−point meshes using a Wannier interpolation based on a single Wannier projection for the highest valence band (in the very same way as we constructed them in vasp for the cRPA calculations). This allows us to accurately calculate the cDFPT Eliashberg function α 2 F(ω) as shown in Fig. 6. From this we calculate the effective electronphonon coupling λ = 2 dω α 2 F(ω) ω which we finally use to fit our local phonon model via α 2 F(ω) ≈ N 0 g 2 δ(ω − ω ph ) with N 0 = 3.88 states/spin/eV/unit cell.

SUSCEPTIBILITY
In this section we show the results for the charge and spin susceptibilities obtained within single-shot D-TRILEX calculations [33,34] in the presence of the electron-phonon coupling. As discussed in the main text, taking into account phonon degrees of freedom shifts the charge density wave (CDW) phase boundary to a smaller value of the filling n = 1.23 compared to the one n = 1.29 obtained in the absence of the electron-phonon coupling. The charge susceptibility calculated close to the CDW phase transition point T = 300 K and n = 1.23 is presented in Fig. 7. The structure of the susceptibility in momentum space consists of six delta-function-like Bragg peaks that appear at K points of the Brillouin zone (BZ). This fact illustrates that the developed charge ordering corresponds to a commensurate CDW. Remarkably, the form of the charge susceptibility calculated with the electron-phonon coupling is similar to the one obtained neglecting phonon degrees of freedom (see main text) and does not show any additional features, like Bragg peaks at twice the Fermi wavevector 2k F , that could be associated with the effect of phonons. This result suggests that the CDW instability in the monolayer InSe is driven by strong electronic Coulomb correlations.
At the CDW transition point the spin susceptibility stays finite and has an antiferromagnetic (AFM) behavior that manifests itself in the largest value of the susceptibility at K point of the BZ (see top left panel of Fig. 8). At a larger filling n = 1.42 the AFM form of the spin susceptibility changes to FIG. 9.
DOS of the monolayer InSe obtained within DMFT (blue line) and D-TRILEX (red line) methods close to a CDW phase boundary at T = 300 K and n = 1.23 taking into account the electron-phonon coupling. a more complex structure with the largest value at Γ point of the BZ, which indicates the ferromagnetic instability (FM), and less pronounced intensities at M points (see top right panel of Fig. 8). Increasing the filling to n = 1.70 the van Hove singularity (vHS) in the electronic spectral function appears at the Fermi energy, which strongly enhances collective electronic fluctuations. Bottom left panel of Fig. 8 shows that close to the magnetic instability the spin susceptibility becomes purely FM, which is indicated by the single deltafunction-like Bragg peak that at the Γ point of the BZ. At even larger filling of the band n = 1.92 the spin susceptibility remains FM, but the value of the susceptibility at Γ point decreases compared to the case of the vHS filling, which confirms the dome-like structure of the FM instability.

LOCAL DENSITY OF STATES
In this section we show the density of states (DOS) obtained via the self-consistent D-TRILEX calculation at the CDW phase transition point T = 300 K and n = 1.23 in the presence of the electron-phonon coupling. The corresponding result is shown in Fig. 9 (red line), and is compared to the one of the dynamical mean-field theory (DMFT) (blue line). Note that DMFT does not take into account spatial correlations, but considers the effect of the local electron-phonon coupling via the renormalized on-site Coulomb potential U * = U − U ph ω . We find that the DOS predicted by DMFT has a peak close to the Fermi energy, which corresponds to the vHS of the Mexican-hat-like band. However, at the CDW transition point collective charge fluctuations are strong. Taking them into account via D-TRILEX method turns this peak into a pseudogap, which confirms the formation of the insulating ordered state in the system. Remarkably, we find that considering phonon degrees of freedom does not change a qualitative structure of the DOS obtained close to the CDW phase transition point.