Introduction

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 conditions1,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 spectroscopy8,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 them10,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 graphene13,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 spectrum16, circumventing the need for any twist engineering.

The theoretical description of many-body effects in 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 excitations17. Also, a strong electron–phonon interaction may result in a charge density wave (CDW) instability predicted recently18. At the same time, monolayer InSe is expected to possess strong many-body electronic effects that have not been accurately studied yet. For instance, the presence of the Mexican-hat-like band in this material might favor 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 based on density functional theory (DFT)19 without systematic consideration of many-body effects. 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 work, 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. 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 electron–phonon 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.

Results

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. 1a). Each layer has D3h point group symmetry. 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. 1b. This shape is characterized by a ring of radius k0 ≈ 0.3 Å−1 formed by the valence band edge, and by a well of E0 = 0.09 eV deep at the Γ point. These parameters are in reasonable agreement with those obtained from photoemission spectroscopy in Ref. 8. This peculiar band structure 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. 1b). The corresponding model is defined in terms of maximally localized Wannier functions on an effective triangular lattice, as shown in Fig. 1a. Each Wannier function is reminiscent of an In-In bonding orbital with tails on Se atoms. The resulting single-band model Hamiltonian on a triangular lattice reads

$$H=\mathop{\sum}\limits_{ij,\sigma }{t}_{ij}^{}{c}_{i\sigma }^{{\dagger} }{c}_{j\sigma }^{}+U\mathop{\sum}\limits_{i}{n}_{i\uparrow }^{}{n}_{i\downarrow }^{}+\frac{1}{2}\mathop{\sum}\limits_{i\ne j,\sigma \sigma ^{\prime} }{V}_{ij}^{}\,{n}_{i\sigma }^{}{n}_{j\sigma ^{\prime} }^{}+{\omega }_{{{{\rm{ph}}}}}\mathop{\sum}\limits_{i}{b}_{i}^{{\dagger} }{b}_{i}^{}+g\mathop{\sum}\limits_{i,\sigma }{n}_{i\sigma }^{}\left({b}_{i}^{}+{b}_{i}^{{\dagger} }\right)$$
(1)

where \({c}_{i\sigma }^{({\dagger} )}\) operator annihilates (creates) an electron on the site i with the spin projection σ = {↑, ↓}. The ab initio electronic dispersion is reproduced by five neighboring hopping amplitudes tij: t01 = 127.9 meV, t02 = − 41.8 meV, t03 = − 45.0 meV, t04 = 13.0 meV, and t05 = −4.4 meV. We neglect spin-orbit coupling because it is not important in truly 2D crystals with the mirror symmetry20, which is also supported by earlier DFT calculations showing that the relevant valence states in monolayer InSe remain unaffected by the presence of spin-orbit coupling21. Figure 1c shows the bare and screened Coulomb interaction between the Wannier orbitals calculated as a function of R within the constrained random phase approximation (cRPA) to avoid a double counting of the screening channels from the correlated subspace (see Methods for more details). R(a) corresponds to the distance between the sites of an effective triangular lattice and is expressed in units of the lattice constant a. The on-site screened Coulomb repulsion U = 1.78 eV greatly exceeds the bandwidth ≈1 eV, which usually indicates strong magnetic fluctuations in the system22. As expected for a 2D material23, the non-local Coulomb interaction Vij in monolayer InSe is weakly screened and long-ranged (see blue line in Fig. 1c). Weak long-range screening in 2D is a direct consequence of the reduced dimensionality. In particular, the softer q-dependence of the bare Coulomb interaction V2D(q) = 2πe2/q (in contrast to ~1/q2 in 3D) ensures that the dielectric constant ε(q, ω) → 1 as q → 0 for any gapped 2D system (see, e.g., Ref. 24). At finite wave vectors this behavior is less obvious, yet the trend is similar23,25,26. Thus, cRPA Coulomb matrix elements for layered materials indeed show as a general trend the described weakly screened and long-ranged behavior with graphene27 and NbS228 being representative metallic examples. Moreover, in monolayer InSe the interaction V01 = 1.04 eV between nearest-neighbor electronic densities \({n}_{i\sigma }={c}_{i\sigma }^{{\dagger} }{c}_{i\sigma }^{}\) is larger than the half of the on-site Coulomb repulsion U. 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 interactions29,30. The full form of the long-range Coulomb interaction is presented in Methods.

Fig. 1: First-principle calculation of model parameters.
figure 1

a Schematic crystal structure of monolayer InSe shown in two projections. Superimposed is an isosurface of the Wannier function describing the valence states in InSe. The black lines depict an effective electronic lattice; (b) Band structure and DOS calculated in DFT (blue), and from tight-binding model parametrization (red); (c) The Coulomb interaction between the Wannier orbitals shown as a function of the distance calculated within the cRPA scheme. The bare (unscreened) Coulomb interaction is shown for comparison.

Full size image

To estimate the phonon properties without double counting screening channels from the correlated subspace we utilize the constrained Density Functional Perturbation Theory (cDFPT)31 at hole doping. We find that the effective electron–phonon coupling \(\lambda =2\int d\omega \,\frac{{\alpha }^{2}F(\omega )}{\omega }\) is dominated by a rather sharp resonance at low phonon frequency, which we can approximate with a local phonon model, i.e α2F(ω) = N0g2δ(ω − ωph). Here, N0 is the DOS at the Fermi level, ωph = 8.5 meV the phonon energy, and g = 34.7 meV its coupling strength to the topmost valence band (see Methods). 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 frequency-dependent attractive interaction \({U}_{\omega }^{{{{\rm{ph}}}}}=2{g}^{2}\frac{{\omega }_{{{{\rm{ph}}}}}}{{\omega }_{{{{\rm{ph}}}}}^{2}-{\omega }^{2}}\) that reduces the repulsive on-site Coulomb potential as \(U\to U-{U}_{\omega }^{{{{\rm{ph}}}}}\)32,33,34 and thus enhances the effect of the non-local Coulomb interaction Vij.

An accurate theoretical investigation of many-body instabilities in InSe monolayer cannot be performed within conventional perturbative methods like the random phase approximation35,36,37 or the GW approach38,39,40. Local correlation effects that are governed by such a large value of the local Coulomb interaction U require a systematic many-body consideration using, e.g., the dynamical mean-field theory (DMFT)41. At the same time, spatial collective electronic fluctuations and the long-range Coulomb interaction cannot be captured by DMFT and require diagrammatic extensions of this theory42. In this work, we solve the considered many-body problem using the dual triply irreducible local expansion (D-TRILEX) method43,44,45 that allows one to account for leading collective electronic fluctuations on equal footing without any limitation on the range46. More details on the theoretical approach are provided in Methods.

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 dispersion8,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 vHS10,11,12. Practically, high concentration of holes of the order of 1014 cm−2 can be achieved in 2D materials by means of electrostatic solid- or liquid-electrolyte gating47 or by surface molecular doping48. First, we solve the many-body problem without considering 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 calculations for the charge and spin susceptibilities \({X}_{{{{\bf{q}}}}\omega }^{\varsigma }\). 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 (see Methods for more details). Fig. 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. Fig. 3a displays the momentum resolved charge susceptibility \({X}_{{{{\bf{q}}}}\omega }^{{{{\rm{c}}}}}\) 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 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 Fig. 3b. It reveals a sharp Bragg peak at the Γ point of the BZ, which indicates the tendency towards FM ordering.

Fig. 2: Phase diagram for the monolayer InSe as a function of temperature and doping.
figure 2

Solid vertical lines correspond to the CDW phase boundaries, dashed lines depict the FM instabilities. Results are obtained in the presence (red line) and in the absence (blue line) of the electron–phonon coupling. The top of the FM dome corresponds to the filling 〈n〉 = 1.70 at which the vHS appears at the Fermi energy. Black arrows with the label “ph” illustrate the effect of phonons that tend to suppress the FM instability and favor the CDW ordering. Orange triangles indicate the CDW transition points 〈n〉 = 1.40 and 〈n〉 = 1.92 obtained at T = 300 K taking into account the electron–phonon coupling in the absence of the non-local Coulomb interaction.

Full size image
Fig. 3: Momentum resolved susceptibilities \({X}_{{{{\bf{q}}}},\omega }^{\varsigma}\) calculated at zero frequency ω = 0 in the absence of the electron–phonon coupling.
figure 3

Results are obtained close to the CDW (T = 300 K, 〈n〉 = 1.29) and the FM (T = 375 K, 〈n〉 = 1.70) instabilities, respectively. Bragg peaks that appear in the charge susceptibility (a) at the K points of the hexagonal BZ correspond to a commensurate CDW ordering. A single peak at the Γ point in the spin susceptibility (b) confirms that the observed instability is FM.

Full size image

Electronic density of states

The obtained phase diagram illustrates that the commensurate CDW represents the main instability in the InSe monolayer contrary to the DFT prediction19. The formation of this 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 quantities49,50 we perform self-consistent D-TRILEX calculations43,44,45. The obtained result is compared to the one of the DMFT that does not take into account spatial correlation effects. Fig. 4 shows the DOS calculated close to the CDW phase boundary (T = 300 K and 〈n〉 = 1.29). We find that the DOS obtained within DMFT (blue) is similar to the one of DFT (black) 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 includes the effect of spatial correlations via the D-TRILEX approach, one observes that at the transition point this peak turns into a pseudogap (red). In this particular case, the pseudogap appears due to strong collective charge fluctuations that lead to an almost diverged charge susceptibility \({X}_{{{{\bf{q}}}}\omega }^{{{{\rm{c}}}}}\) at the q = K point of the BZ close to a phase transition. As a consequence, the renormalized interaction \({W}_{{{{\bf{q}}}}\omega }^{{{{\rm{c}}}}}\) 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 state51,52,53 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 short-range order without long-range order leads to a pseudogap in the electronic spectrum54,55.

Fig. 4: DOS of the InSe monolayer.
figure 4

Results are obtained within DMFT (blue), D-TRILEX (red), and DFT (black) methods close to a CDW phase boundary at T = 300 K and 〈n〉 = 1.29 without taking into account the electron–phonon coupling. Note that the DOS of DFT is divided by a factor of 2 for easier comparison.

Full size image

Effect of the electron–phonon coupling

In order to investigate the effect of phonons on the observed instabilities we repeat the same calculation in the presence of the electron–phonon coupling. As Fig. 2 shows, 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 fact 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 single occupation of lattice sites. On the contrary, the non-local Coulomb interaction promotes the CDW ordering, which upon reducing the local Coulomb interaction becomes energetically preferable.

Remarkably, Fig. 5 demonstrates that taking into account the electron–phonon coupling does not change the position of Bragg peaks in the charge susceptibility calculated close to the CDW phase transition point T = 300 K and 〈n〉 = 1.23. As in the absence of phonons, the momentum-space structure of this susceptibility \({X}_{{{{\bf{q}}}}}^{{{{\rm{c}}}}}(\omega =0)\) consists of six delta-function-like Bragg peaks that appear at K points of the BZ. This fact illustrates that the developed charge ordering also corresponds to a commensurate CDW. Moreover, we do not observe any additional Bragg peaks in the charge susceptibility at twice the Fermi wave vector q* = 2kF, which would indicate a tendency to the charge ordering due to the electron–phonon mechanism.

Fig. 5: Momentum resolved charge susceptibility Xc(q, ω = 0).
figure 5

Result is obtained close to the CDW transition point T = 300 K and 〈n〉 = 1.23 taking into account the electron–phonon coupling.

Full size image

The effect of phonons can be revealed only upon excluding the contribution of the non-local Coulomb interaction Vq from the charge susceptibility. In this case, an additional exclusion of phonons eliminates both sources of the CDW instability, and \({X}_{{{{\bf{q}}}}\omega }^{{{{\rm{c}}}}}\) remains finite at any filling. If the electron–phonon coupling is taken into account, the maximum value of \({X}_{{{{\bf{q}}}}}^{{{{\rm{c}}}}}(\omega =0)\) corresponds to the wave vector q* that changes with doping. For instance, at 〈n〉 = 1.23 this wave vector lies between Γ and K points of the BZ (Fig. 6a). However, the value of this susceptibility is negligibly small compared to the almost diverged one calculated in the presence of the non-local Coulomb interaction (see Fig. 5). Increasing the filling, the maximum of the charge susceptibility calculated for the case of Vq = 0 also increases, and the corresponding wave vector q* shifts towards the K point of the BZ (Fig. 6b). At 〈n〉 ≃ 1.40 (T = 300 K) the q* reaches the K point (Fig. 6c), and the charge susceptibility diverges. This means that the CDW instability mediated purely by phonons is shifted to much larger values of the filling compared to the one driven by strong Coulomb correlations (orange triangles in Fig. 2). Increasing the filling even more, the charge susceptibility becomes finite again at 〈n〉 = 1.92 (Fig. 6d). Remarkably, Figs. 6c and d show that Bragg peaks that correspond to CDW ordering vectors q* are drastically different for two fillings 〈n〉 = 1.40 and 〈n〉 = 1.92. This means that the ordering vector of the phonon-driven CDW also changes inside the CDW phase. These results confirm that the metal-to-CDW phase transition in the monolayer InSe is driven by strong Coulomb correlations rather than by the electron–phonon mechanism as suggested previously18.

Fig. 6: Momentum resolved charge susceptibility Xc (q) calculated at zero frequency ω = 0.
figure 6

Results are obtained at T = 300 K for different values of the filling 〈n〉 = 1.23 (a), 〈n〉 = 1.33 (b), 〈n〉 = 1.40 (c), and 〈n〉 = 1.92 (d) taking into account the electron–phonon coupling and excluding the effect of the non-local Coulomb interaction (Vq = 0).

Full size image

The momentum-structure of the spin susceptibility also changes with doping. At the CDW transition point it stays finite and has a Néel antiferromagnetic (AFM) behavior that manifests itself in the largest value of the susceptibility at K point of the BZ (Fig. 7a). At a larger filling 〈n〉 = 1.42 the Néel AFM form of the spin susceptibility changes to a more complex structure with the largest value at Γ point of the BZ, which indicates the FM instability, and less pronounced intensities at M points corresponding to row-wise AFM fluctuations (Fig. 7b). 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. Fig. 7c shows that close to the magnetic instability the spin susceptibility becomes purely FM, which is indicated by the single delta-function-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 (Fig. 7d), which confirms the dome-like structure of the FM instability.

Fig. 7: Momentum resolved spin susceptibility Xs(q) calculated at zero frequency ω = 0 taking into account the electron–phonon coupling.
figure 7

Results are obtained at T = 300 K for different values of the filling 〈n〉 = 1.23 (a), 〈n〉 = 1.42 (b), 〈n〉 = 1.70 (c), and 〈n〉 = 1.92 (d).

Full size image

For completeness, we also show the electronic DOS calculated close to the CDW phase transition point T = 300 K and 〈n〉 = 1.23 in the presence of the electron–phonon coupling (Fig. 8). We find that in this case the effect of spatial correlations on the DOS is qualitatively the same as in the absence of phonons. Note that in this case DMFT considers the effect of the local electron–phonon coupling via the renormalized on-site Coulomb potential \({U}^{* }=U-{U}_{\omega }^{{{{\rm{ph}}}}}\).

Fig. 8: DOS of the monolayer InSe.
figure 8

Results are obtained within DMFT (blue), D-TRILEX (red), and DFT (black) methods close to a CDW phase boundary at T = 300 K and 〈n〉 = 1.23 taking into account the electron–phonon coupling. Note that the DOS of DFT is divided by a factor of 2 for easier comparison.

Full size image

Discussion

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 non-local 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. The presence of this ordering 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 phonons. 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.

In general, the Coulomb-driven CDW instability in electronic systems appears as the result of a competition between local and non-local electron-electron interactions that favor different distribution of electrons on a lattice. In order to change the balance between these interactions, the free standing InSe monolayer considered here could be exposed to external dielectric screening as, e.g., resulting from a substrate or from encapsulation. This mechanism is especially effective in 2D as the environmental screening is non-local in layered materials, which reduces long-ranged interactions stronger than the on-site interaction25,26,56. This kind of Coulomb engineering paves the way for exploring additional degrees of freedom to tune the correlation effects in this 2D material.

Methods

First-principle calculations

The band structure presented in Fig. 1b was calculated within density functional theory utilizing the projected augmented wave (PAW) formalism57,58 as implemented in the Vienna ab initio simulation package (vasp)59,60 version 5.4. The exchange-correlation effects were considered using the generalized gradient approximation (GGA)61 and the standard In and Se pseudopotentials in version 5.4. A 500 eV energy cutoff for the plane-waves and a convergence threshold of 10−7 eV were used in the calculations together with tetrahedron method for all involved integrals. 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 Å, In-In vertical distance of 2.76 Å, and Se–Se vertical distance of 5.37 Å. 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 localization62,63 using the wannier90 package64 version 1.2.

The Coulomb interaction was evaluated using the maximally localized Wannier functions within the constrained random phase approximation (cRPA)65,66 as Uij = 〈wiwj∣U∣wjwi〉, 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 wi 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 vasp66. To converge the cRPA polarization with respect to the number of empty states with used in total 64 bands. To derive a light-weighted Coulomb model for arbitrary q-grids, we fitted the cRPA Coulomb interaction in momentum space according to:

$$U(q)=v(q)/\varepsilon (q)$$
(2)

with the bare Coulomb interaction of a monolayer

$$v(q)=\frac{h}{2\pi }\int\nolimits_{-\pi /h}^{+\pi /h}\frac{4\pi {e}^{2}}{V{q}^{2}}d{q}_{z}=\frac{4{e}^{2}}{A}\frac{\arctan \left(\frac{\pi }{qh}\right)}{q}$$
(3)

and the dielectric function

$$\varepsilon (q)={\varepsilon }_{0}(q)\frac{{\varepsilon }_{0}(q)+1-({\varepsilon }_{0}(q)-1){e}^{-qd}}{{\varepsilon }_{0}(q)+1+({\varepsilon }_{0}(q)-1){e}^{-qd}}$$
(4)

with

$${\varepsilon }_{0}(q)=\frac{a+{q}^{2}}{\frac{a\sin (qc)}{bqc}+{q}^{2}}.$$
(5)

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. 9 we show the corresponding ab initio values together with their fits using h = 15.35 Å, a = 1.07 Å−2, b = 5.79, c = 0.14 Å, and d = 8.09 Å. In this expression ε0(q) represents the (semi-conducting) cRPA screening in the middle of the film, i.e., z = 0, which is assumed to be the same for the entire film thickness d and which we parameterized following Ref. 67. The surrounding of the film is supposed to be vacuum. This defines a classical electrostatics problem, which can be analytically solved (see e.g., Ref. 68). The full non-local cRPA background screening than takes the form of Eq. (4). More details on this can be found in Ref. 26. The real-space Uij are calculated via a conventional Fourier transform of U(q) from Eq. (2) and using these fits.

Fig. 9: Coulomb interaction from first-principle calculations.
figure 9

Bare Coulomb interaction (a) and cRPA dielectric function (b) within the Wannier basis. Markers represent ab initio data and dashed lines the corresponding fits.

Full size image

The phonon properties are calculated within Quantum Espresso69 version 6.5 using norm-conserving pseudopotentials generated using the “atomic” code by A. Dal Corso (v.5.0.99 svn rev. 10869, 2014), the local density approximation, an energy cutoff of 80 Ry, and a finite hole doping of 0.04 holes per unit cell. For the initial constrained Density Functional Perturbation Theory (cDFPT)31 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 code70 version 5.0.0 to extrapolate the cDFPT results to (32 × 32) k- and 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 α2F(ω) as shown in Fig. 10. From this we calculate the effective electron–phonon coupling \(\lambda =2\int d\omega \,\frac{{\alpha }^{2}F(\omega )}{\omega }\) which we finally use to fit our local phonon model via α2F(ω) ≈ N0g2δ(ω − ωph) with N0 = 3.88 eV−1 per spin per unit cell.

Fig. 10: Phonon properties from first-principle calculations.
figure 10

cDFPT Eliashberg function together with the partial effective electron–phonon coupling \({\lambda }_{p}(\omega )=2\int\nolimits_{0}^{\omega }d\omega ^{\prime} \,\frac{{\alpha }^{2}F(\omega ^{\prime} )}{\omega ^{\prime} }\).

Full size image

Many-body calculations

The introduced many-body electronic model (1) is solved using the finite-temperature diagrammatic D-TRILEX approach43,44,45. In this method, local correlation effects are treated via the self-energy \({{{\Sigma }}}_{\nu }^{{{{\rm{imp}}}}}\) and the polarization operator \({{{\Pi }}}_{\omega }^{\varsigma \,{{{\rm{i}}}}mp}\) of an effective local impurity problem of DMFT. These quantities are written in the fermionic (ν) and bosonic (ω) Matsubara frequency space, respectively. The corresponding impurity problem is solved numerically exactly using the open source CT-HYB solver71,72 based on ALPS libraries73.

Within the self-consistent calculations, the spatial fluctuations are considered in a partially bosonized form of the renormalized charge (ς = c) and spin (ς = s) interactions \({W}_{{{{\bf{q}}}}\omega }^{\varsigma }\) that enter the diagrammatic part of the self-energy \({\overline{{{\Sigma }}}}_{{{{\bf{k}}}}\nu }\) introduced beyond DMFT43,44,45. The dressed Green’s function Gkν of the problem can be found using the standard Dyson equation \({G}_{{{{\bf{k}}}}\nu }^{-1}=i\nu +\mu -{\varepsilon }_{{{{\bf{k}}}}}-{{{\Sigma }}}_{{{{\bf{k}}}}\nu }\) written in momentum k 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 tij, and \({{{\Sigma }}}_{{{{\bf{k}}}}\nu }={{{\Sigma }}}_{\nu }^{{{{\rm{imp}}}}}+{\overline{{{\Sigma }}}}_{{{{\bf{k}}}}\nu }\) is the total self-energy. The renormalized interaction \({W}_{{{{\bf{q}}}}\omega }^{\varsigma }\) can be found via the following Dyson equation \({W}_{{{{\bf{q}}}}\omega }^{\varsigma \,-1}={U}_{{{{\bf{q}}}}\omega }^{\varsigma \,-1}-{{{\Pi }}}_{{{{\bf{q}}}}\omega }^{\varsigma }\), where \({U}_{{{{\bf{q}}}}\omega }^{{{{\rm{c}}}}}=U/2-{U}_{\omega }^{{{{\rm{ph}}}}}+{V}_{{{{\bf{q}}}}}\) and \({U}_{{{{\bf{q}}}}\omega }^{{{{\rm{s}}}}}=-U/2\) are the bare interactions in the charge and spin channels, respectively43,44,45. \({{{\Pi }}}_{{{{\bf{q}}}}\omega }^{\varsigma }={{{\Pi }}}_{\omega }^{\varsigma \,{{{\rm{imp}}}}}+{\overline{{{\Pi }}}}_{{{{\bf{q}}}}\omega }^{\varsigma }\) is the total polarization operator of the problem, where \({\overline{{{\Pi }}}}_{{{{\bf{q}}}}\omega }^{\varsigma }\) is the diagrammatic contribution introduced in the D-TRILEX approach43,44,45. In this way, the dressed lattice Green’s function Gkν takes into account both, local and non-local correlation effects. To obtain the electronic DOS, we take the local part of the lattice Green’s function \({G}_{\nu }^{{{{\rm{loc}}}}}=\frac{1}{N}{\sum }_{{{{\bf{k}}}}}{G}_{{{{\bf{k}}}}\nu }\) and perform an analytical continuation from Matsubara frequency space (ν) to real energies (E) using the stochastic optimization method74.

In the single-shot D-TRILEX calculations for the susceptibility the diagrammatic part of the polarization operator \({\overline{{{\Pi }}}}_{{{{\bf{q}}}}\omega }^{\varsigma }\) is obtained non-self-consistently. In this case, \({\overline{{{\Pi }}}}_{{{{\bf{q}}}}\omega }^{\varsigma }\) accounts only for DMFT Green’s functions \({G}_{{{{\bf{k}}}}\nu }^{{{{\rm{D}}}}}\) that are dressed by the local self-energies: \({G}_{{{{\bf{k}}}}\nu }^{{{{\rm{D\,-1}}}}}=i\nu +\mu -{\varepsilon }_{{{{\bf{k}}}}}-{{{\Sigma }}}_{\nu }^{{{{\rm{imp}}}}}\). Charge and spin susceptibilities \({X}_{{{{\bf{q}}}}\omega }^{\varsigma }\) can then be obtained straightforwardly as \({X}_{{{{\bf{q}}}}\omega }^{\varsigma \,-1}={{{\Pi }}}_{{{{\bf{q}}}}\omega }^{\varsigma \,-1}-{U}_{{{{\bf{q}}}}\omega }^{\varsigma }\)45. We note that this form of the D-TRILEX susceptibility resembles the DMFT susceptibility75,76,77 with a longitudinal dynamical vertex corrections44.