Competing Coulomb and electron–phonon interactions in NbS2

The interplay of Coulomb and electron–phonon interactions with thermal and quantum fluctuations facilitates rich phase diagrams in two-dimensional electron systems. Layered transition metal dichalcogenides hosting charge, excitonic, spin and superconducting order form an epitomic material class in this respect. Theoretical studies of materials like NbS2 have focused on the electron–phonon coupling, whereas the Coulomb interaction, particularly strong in the monolayer limit, remained essentially untouched. Here, we analyze the interplay of short- and long-range Coulomb as well as electron–phonon interactions in NbS2 monolayers. The combination of these interactions causes electronic correlations that are fundamentally different to what would be expected from the interaction terms separately. The fully interacting electronic spectral function resembles the non-interacting band structure but with appreciable broadening. An unexpected coexistence of strong charge and spin fluctuations puts NbS2 close to spin and charge order, suggesting monolayer NbS2 as a platform for atomic scale engineering of electronic quantum phases. A theoretical study taking account of electron–electron interactions and electron–phonon interactions reveals enhanced quantum fluctuations in monolayer NbS2. Erik van Loon from the Institute for Molecules and Materials at Radboud University and an international team of collaborators perform theoretical studies to treat both electron–phonon interactions and electronic correlations resulting from competing short- and long-range Coulomb interactions in monolayer NbS2. They reveal a simultaneous enhancement of the charge and spin susceptibilities and a sharp transition from tendencies of preferential spin ordering to charge ordering. The results call for caution on the treatment of properties such as magnetic and charge susceptibilities as the electronic correlation effects will come into play beyond the single-particle level, despite a match between photoemission experiments and density functional theory (DFT) calculations. They suggest monolayer NbS2 as a platform for atomic scale engineering of electronic quantum phases.


INTRODUCTION
Layered materials host in many cases pronounced electronic interaction phenomena ranging from eV-scale excitonic binding energies in semiconductors to charge, spin, and superconducting order in metallic systems. Characteristic energy scales and transition temperatures associated with these interaction phenomena change often remarkably when approaching the limit of atomically thin materials. [1][2][3] Two generic contributions to these material-thickness dependencies are quantum-confinement [4][5][6] and enhanced local and long-ranged Coulomb interactions in monolayer thin materials. [7][8][9][10][11][12] In addition, many layered materials feature sizable electron-phonon coupling. [13][14][15][16] The resulting interplay of interactions, which are effective at different length and time scales (see Fig. 1), makes the phase diagrams of twodimensional materials and their response to external stimuli very rich.
The layered metallic transition metal dichalcogenides (TMDC), 17,18 MX 2 , where M denotes one of the transition metals V, Nb, or Ta and X stands for one of the chalcogens S or Se, presents a demonstrative case in this respect where the monolayer limit is becoming experimentally accessible. 2,17,18 Within this material, class a competition of charge-and spinordered, Mott insulating, as well as superconducting states can be found. Here, the V-based compounds show tendencies toward magnetic 19,20 as well as charge order [21][22][23] in their monolayer and bulk phases, respectively, which might partially coexist in the fewlayer limit. 20,24 In contrast, the sub-class of Ta-based compounds 1,25-30 as well as NbSe 2 2,3,31-38 show a competition between charge-density waves, superconducting as well as Mott-insulating states. NbS 2 appears to be a border case. It is superconducting in the bulk 39-41 but does not display any chargedensity wave (CDW) formation there. In the case of few-layer NbS 2 , first experimental studies reported recently metallic transport properties down to three layers, 42 whereas mean-field calculations reveal a tendency to form magnetic states. 43,44 NbS 2 is thus likely on the verge between different instabilities. Whether or how these instabilities are triggered by the interplay of the involved interactions is barely understood, up to now. Although a lot of focus was put on the investigation of electron-phonon coupling effects in the whole class of metallic TMDCs, the effects of the subtle interplay of the local and non-local Coulomb interaction terms have been mostly neglected, so far. It is thus necessary to draw our attention also to these short and long-range electron-electron interactions. Unfortunately, there is no theory that can handle, even qualitatively, the competition of these strong interactions beyond the perturbative regime, yet. To overcome this problem, we combine here the Dual Boson formalism [45][46][47] with first-principles approaches and construct a state-of-the-art material-realistic theory of monolayer NbS 2 , which properly treats electronic correlations as resulting from competing short-and long-range Coulomb interactions. Thereby, we also account for the electron-phonon interactions to gain a universal understanding of all interaction effects.
Our calculations reveal a simultaneous enhancement of the charge and spin susceptibilities owing to the various interactions in monolayers of NbS 2 and a sharp transition from tendencies of preferential spin ordering to charge ordering. Despite these strong interaction effects, the electronic spectral function as measured, e.g., in angularly resolved photoemission (ARPES) experiments largely resembles the non-interacting dispersion in accordance with the available experimental data. We trace this back to a compensation of the different interaction terms, which are partially effective on the single-particle but not on the twoparticle level. From an experimental perspective, this means that finding a match between ARPES results and density functional theory (DFT) bands is not sufficient to rule out strong correlation/ interaction effects, as the competition of the interactions masks correlation effects on the single-particle level, whereas they are still visible on the two-particle level, i.e., in the magnetic and charge susceptibility.

Competing interactions in NbS 2
The non-correlated band structure of NbS 2 monolayers exhibits a half-filled metallic band surrounded by completely filled valence bands 1 eV below and completely empty conduction bands 3 eV above the Fermi level, as shown in Fig. 1a. This motivates a description of the competing interaction effects in terms of an extended Hubbard-Holstein model 48 for the separated metallic band only where c y iσ and c iσ are the creation and annihilation operators of the electrons with spin σ on lattice site i, b y and b are the creation and annihilation operators of a local phonon mode, and n i ¼ c y iσ c iσ is the electron occupation number operator. This model includes on the single-particle level the electron hopping t ij and a local phonon mode with energy ω ph . We include an on-site Coulomb repulsion U and long-range Coulomb interactions V ij , as well as an electron-phonon coupling g which couples the local charge density to the given phonon mode. The latter can actually be integrated out which results in a purely electronic model with an effective dynamic local interaction that is lowered and thereby effectively screened by the phonons. [48][49][50] This treatment of the phonons as simple singlefrequency modes that are coupled locally to the electrons is an assumption necessary to keep the problem tractable. Otherwise, the non-local interaction V ij or V q in momentum space would also become frequency dependent. As we argue in the Methods section, there is a basis for this assumption, however, the simplification might change the exact position at which instabilities occur in the Brillouin zone. Furthermore, we would like to emphasize that, in our treatment, electronically generated phonon anharmonicities are automatically included, whereas bare phonon anharmonicities are not. To realistically describe NbS 2 monolayers, we derive the parameters entering Eq. (1) from first principles. Therefore, we generate a most accurate tight-binding model describing the metallic and the lowest two conduction bands in a first step, and use it afterwards to perform calculations within the constrained Random Phase Approximation (cRPA) 51 to obtain the partially screened Coulomb interaction matrix elements within the same basis. The phonon frequency and the electron-phonon coupling are estimated based on density-functional perturbation theory calculations. The resulting three-band model is subsequently simplified in order to get the final single-band model describing the metallic band only as explained in the method section.
Our ab initio simulations yield an effective local Coulomb interaction U ≈ 1.8 eV, a nearest-neighbor interaction V ≈ 1 eV as well as further long-range interaction terms. The typical bare phonon frequency ω ph and the electron-phonon coupling g for this material are estimated to be 20 and 70 meV, respectively (see methods for further details). Notably, both, the on-site Coulomb repulsion and the effective electron-electron attraction λ = 2g 2 / ω ph = 0.5 eV, are on the order of the electronic band width ≈1.2 eV, as sketched in Fig. 1.
Spectral fingerprints of the interactions Each interaction term on its own can thus trigger strong electronic correlations, which becomes evident from the electronic spectral functions shown in Fig. 2a. These have been calculated using the Dual Boson (DB) method taking into account each interaction term on its own and their combined effects. For local Coulomb interactions only (top left), the half-filled conduction band clearly splits into two Hubbard bands above and below the Fermi level. There is no spectral weight at the Fermi level and the system is insulating. Including the non-local Coulomb interaction terms (top right) markedly changes the spectral function. The lower Hubbard band still retains noticeable spectral weight. However, the upper Hubbard band overlaps now with a broad distribution of spectral weight reaching the Fermi level. That is, the non-local Coulomb interaction drives the system into the state of a correlated metal, similar to what has been shown for graphene. 52 With only electron-phonon interaction (bottom left), the spectrum is also reminiscent of a correlated metal, with again strong spectral weight transfer away from the Fermi level toward polaronic bands at higher energies. Finally, simultaneous inclusion of all interactions yields the spectral function shown in the bottom right of Fig.  2a. We find a single band with a dispersion very similar to the DFT result of Fig. 1 (shown as a red line). Seemingly, the different interaction terms largely compensate each other despite the fact that they are effective at very different length and time scales. The major interaction effect visible in, e.g., ARPES experiments is that the band widens significantly compared to the thermal broadening inherent to any finite temperature measurement. ARPES experiments in other transition metal dichalcogenides monolayers (TaS 2 53 and NbSe 2 37 ) are consistent with this picture: the dispersion follows roughly the DFT band structure with some broadening. A more detailed comparison with our results is, however, not possible, as in the experiments different materials have been used, lower temperatures were applied, and substrates were present.
Our material-specific results can be compared with theoretical findings in model systems. In the Hubbard-Holstein model (U and λ in our language) on the triangular lattice, Mott and polaronic insulating states have been found, 54,55 consistent with our results here. On the other hand, the combination U + V + λ presented here has so far not been studied since there were previously no methods that can deal with these competing interactions in the strongly correlated regime, where vertex corrections beyond GW are important. Therefore, we need to scrutinize this behavior in more detail and examine the local self-energy, which induces all correlation effects. In Fig. 2b, we show these self-energies corresponding to the spectra in Fig. 2a. If we take only local Coulomb interactions into account (stars), we find a strongly enhanced self-energy for small frequencies. By including also long-range Coulomb interactions (circles), the self-energy is reduced around small frequencies. This trend is continued by including electron-phonon interactions (squares), which demonstrates how the long-range Coulomb and the electron-phonon interactions compensate the effects of the local Coulomb interaction. The full self-energy including all interaction terms is thus strongly reduced around small frequencies, but has still sizeable contributions at all energies considered here, which results in the broadened spectral function without significant reshaping. It is interesting to note that when taking only the electron-phonon interactions into account (triangles), the self-energy, and thus the degree of correlation increases by increasing the electron-phonon coupling g. However, in the presence of Coulomb interaction an enhanced electron-phonon coupling necessarily leads to a decrease of the self-energy and hence to a decreased degree of correlation. Thus, the effect of electron-phonon coupling is the exact opposite depending on whether or not the Coulomb interaction is present in the model. It is therefore absolutely crucial to take all interactions simultaneously into account.
Competition of charge and spin fluctuations in NbS 2 monolayers The electronic correlations as resulting from the interplay of the electron-electron and electron-phonon interaction also manifest in the local two-particle correlation functions of the system, which are shown as a function of the electron-phonon coupling g in Fig.  3. These local observables are calculated directly from the DB auxiliary single-site system. The ratio of the static local charge and spin susceptibilities (lower panel; note the logarithmic scale) and the instantaneous double occupancy (upper panel) vary strongly   Fig. 1a. The top panel is calculated at 232 K, the other panels are calculated at 464 K. b The corresponding local self-energies as functions of Matsubara frequencies. Results for two different electron-phonon coupling strengths g 1 = 0.057 eV and g 2 = 0.07 eV are shown in red and blue, respectively as a function of g. Without electron-phonon coupling (g = 0) the system shows typical signs of strong Mott-Hubbard correlation effects: the spin susceptibility is orders of magnitude larger than the charge susceptibility and the probability of finding two electrons at the same site is greatly reduced in comparison to the value of n " n # = 0.25 found in non-interacting half-filled systems. Turning on the electron-phonon interaction screens the local Coulomb interaction according to Eq. (2), and makes the system less correlated. At sufficiently large g ≈ 70-80 meV, the susceptibility ratio and double occupancy even exceed their noninteracting values of 1 and 0.25 (gray lines), respectively. The numerical simulations get unstable close to a transition to the CDW phase. This is why we could perform simulations at 464 K only up to g = 70 meV (red circles). At a higher temperature of T = 2321 K (orange pluses), larger values of λ can be reached. For g = 40-70 meV, the two data sets agree reasonably well, at higher temperatures the double occupancy is less suppressed and the spin susceptibility is substantially smaller, see also the Methods section. Note that 464 K (0.04 eV) is well below the energy scales defined by the band width and the interactions.
Freestanding NbS 2 monolayer, with g ≈ 70 meV as estimated in the methods, thus turn out to be on the verge to form a CDW ground state. The local properties presented in Fig. 3 also show what happens when both, the electron-phonon interaction and the non-local parts of the Coulomb interaction, are ignored. In that case (dashed green lines), the susceptibility ratio goes down another order of magnitude, and the double occupancy decreases to almost zero. These are all characteristics of a Mott-insulating phase. The local Hubbard interaction U is thus in principle strong enough to create an interaction-driven insulator, with a large spin susceptibility, local magnetic moments (small n " n # ) and strongly suppressed charge fluctuations, as we have already anticipated in the discussion of the spectral functions. Only through screening by the non-local Coulomb contributions, and by the electron-phonon coupling, can the system exhibit the large charge fluctuations (local "charge moments", large n " n # ), that are necessary for a CDW. This shows that both, the Hubbard interaction U and the interactions that screen it, are nonperturbatively large in the freestanding monolayer, which casts doubt on approaches that do not explicitly include all interaction terms. Most importantly, the transition from the regime, which is dominated by spin fluctuations to the charge-fluctuation dominated regime is very abrupt as the steep rise of the susceptibility ratio demonstrates. The strong fluctuations in different channels, around g ≈ 70 meV, signal the close proximity of competing charge and spin order and is indeed ubiquitous in correlated electron systems. 56,57 Next, we turn to the static momentum-resolved susceptibilities. The non-interacting susceptibility of the single-band model, χ 0 , shown in Fig. 4a, agrees with previously published data for NbS 2 monolayers. 44 In a non-interacting system the charge and spin susceptibility would be the same and coincide with χ 0 . This is clearly not the case for the charge and magnetic susceptibilities resulting from our DB calculations shown in Fig. 4b-e.
Without electron-phonon coupling (g = 0) the spin susceptibility is enhanced indicating the presence of strong spin fluctuations. The charge susceptibility, on the other hand, is suppressed in the entire Brillouin zone owing to the Coulomb interaction, which is in line with the expectations for a correlated metal. Turning on the electron-phonon interaction (g = 70 meV) reduces the spin susceptibility, which is however still comparable to χ 0 . At the same time, the charge susceptibility is strongly enhanced and is almost divergent at large momenta. At this point it is, however, important to note that the exact position of the ordering vector might change when the ordered phase is actually entered (here, we investigate just its onset based on the susceptibility in the normal phase) and when a more realistic phonon model is used. Nevertheless, these two observations show again one of our main findings that the interactions partially compete and screen each other, leading to a spin susceptibility that is only moderately enhanced. Most importantly, this competition does not lead to a complete cancellation, as is visible in the strong enhancement of the charge susceptibility. Owing to the interplay of these interactions a strong spin and charge response can thus coexist in this system.

DISCUSSION
Using a combination of the Dual Boson approach and ab initio calculations, we investigated the interplay between the Coulomb and electron-phonon interactions in NbS 2 monolayers and the resulting degree of electronic correlations. We found that both, the Coulomb and the phonon-mediated electron-electron interaction energies, are on the same order as the electronic band width allowing both of them to trigger strong electronic correlations. Both types of interactions on their own would drive NbS 2 to the verge of an insulating state, as our analysis of electronic self-energies shows. Remarkably, our simulations with Coulomb and electron-phonon interactions present yield a spectral function, which closely resembles the non-interacting band structure. Yet, in this situation electron correlations have not fully ceased but manifest themselves in a sizeable broadening of the spectral function. In this sense, NbS 2 is very similar to so-called Hund's metals, [58][59][60][61][62] where the exchange coupling drives the electronic system of materials like Fe-based superconductors away from the Mott Hubbard insulating limit into a correlated metallic phase.
For NbS 2 , the observed spectral broadening argues against simple nesting scenarios based on the bare bands for the CDW instabilities. Our findings rather show that the interplay between all interactions is responsible for driving the system in close proximity of a charge-ordered state. The interaction-induced correlations result in strongly modified spin and charge susceptibilities compared with the non-interacting one. Specifically, we found that the competition between the long-range Coulomb and the electron-phonon interactions is responsible for NbS 2 monolayers being on the edge between dominating spin-and charge fluctuations. The transition from a preferential spin order to charge order is thereby abruptly driven by the electron-phonon interaction.
The resulting ground state is thus heavily dependent on the detailed balance between the internal interactions. To study and test this behavior experimentally, there are several points to be aware of. First, from matching ARPES and DFT data one cannot deduce that the mean-field calculation captures the main physics. Many-body effects lead to quasiparticle broadening as well as enhanced magnetic and charge susceptibilities, which can be measured directly in resonant X-ray scattering 63 or electron energy loss spectroscopy. 64 Importantly, the predicted close vicinity of charge order and local magnetic moment formation can be experimentally tested. The electronic system of NbS 2 can be manipulated via environmental screening or strain applied to the monolayer. Increasing the former will mostly reduce the longrange Coulomb interaction V, while the electron-phonon interaction remains largely untouched. Thereby, the effective screening of the local U due to the non-local V is reduced and the spin susceptibility should be enhanced. By applying strain the electron-phonon interaction can be varied without drastic changes to the Coulomb interaction. By increasing the electron-phonon interaction the system would be pushed into a charge-ordered state.
Our calculations show that correlation effects are particularly prominent in the simultaneously enhanced spin and charge susceptibilities of NbS 2 . Hence, we expect a strong response of the material to local perturbations, which can be experimentally realized through charged as well as magnetic adsorbates.  Fig. 4 to be probed in real space but could possibly also exploit NbS 2 as a platform for quantum engineering, where one switches locally between spin and charge order. The combination of the lattice structure of NbS 2 , the sizeable spin-orbit coupling, and the enhancement of spin and charge susceptibilities clearly away from the Brillouin zone center suggests that the competing ordering tendencies is likely subject to frustration effects. For these reasons, monolayer NbS 2 deserves future exploration not only in the light of fundamental interest but possibly also in relation to concepts of miniaturized neuromorphic computing. 67 Finally, our findings allow to speculate about possible superconducting properties. In this context, the enhanced spin and charge susceptibilities point toward interesting unconventional paring mechanisms. At the same time there are the before mentioned striking similarities between the self-energy in NbS 2 and the one found in Hund's metals, yielding a similar scenario as in Fe-based superconductors. In addition, it needs to be pointed out that the appearance of a CDW phase is usually detrimental to superconductivity so that it likely needs to be suppressed to enhance T c . Given the complicated competition between magnetic and charge instabilities in NbS 2 , an analysis including all of these aspects needs to be carried out to gain reliable insights into possible superconducting properties.

Parametrization of the extended Hubbard-Holstein model
All model parameters are derived from first principles based on DFT and cRPA calculations. To do so, we start with a DFT calculation in Fleur 68 for NbS 2 using a lattice constant of a = 3.37 Å, a k mesh of 18 × 18 × 1, a vacuum height of 32 Å, a relaxed sulfur-sulfur distance of Δ = 3.13 Å, and using FLAPW l-expansion cutoffs of 10 (Nb) and 8 (S) and muffin tin radii of 2.58 a 0 (Nb) and 2.01 a 0 (S) to calculate the band structure shown in Fig. 1a. As the spin-orbit coupling leads to severe spin splittings at the K point only, but not around the Fermi level in NbS 2 we neglect it in the following. Afterwards, we construct a three-band tight-binding model by projecting the original DFT wave functions onto the three dominant niobium orbitals (d z 2 , d xy , d x 2 Ày 2 ) using the Wannier90 code, 69 whereby we ensure that the bands are properly disentangled. To preserve the orbital character of the Wannier functions, we do not perform maximal localization. The resulting three-band tight-binding model perfectly interpolates the original DFT band structure and can be used to evaluate the electronic dispersion at arbitrary k points.
The long-range Coulomb interaction is parametrized in a materialrealistic manner using the cRPA method. 51 Therefore, we start with the fully screened dynamic Coulomb interaction W(q, ω) in reciprocal space that is defined by where v(q) ∝ 1/q is the bare interaction in two dimensions and Π(q, ω) is the polarization function rendering all screening processes. According to the cRPA we can reformulate the latter Π(q, ω) ≈ Π mb (q, ω) + Π rest (q) by splitting it into a dynamic part arising from the half-filled metallic band (mb) and a static part resulting from the rest of the band structure. This is appropriate as we are interested in the low-frequency properties of Π(q, ω) and W(q, ω) only, which are completely rendered by the metallic band and thus by Π mb (q, ω). Using this formulation of the full polarization we can rewrite the fully screened interaction as follows with U(q) being the partially screened Coulomb interaction defined by As described in the Supplemental Methods, U(q) needs be evaluated within the same orbital basis as used for the tight-binding dispersions, using 3 × 3 matrices to represent the bare interaction v(q) and the dielectric function ε(q). Importantly, we can fit analytic expressions to all of the involved matrix elements U αβ (q), allowing us to evaluate U(q) at arbitrary q vectors.
In order to derive a single-band model, we neglect the orbital dependencies in a next step. In this case, the dynamic polarization matrix of the metallic band may be approximated via Π mb ðq; ωÞ ¼ 1 9 Π sb ðq; ωÞ where Π sb (q, ω) is the single-band polarization, which is going to be evaluated in the Dual Boson calculations. The factor 1 9 approximates the overlap matrix elements, which are in general orbital and momentum dependent. This is appropriate for small q and as long as all orbital weights are more or less the same. We found that this assumption is indeed valid in the half-filled situation discussed here. Using this, polarization corresponds to a single-band/orbital partially screened Coulomb interaction defined by U sb ðqÞ ¼ 1 9 X αβ U αβ ðqÞ: Thus, under the assumption of vanishing orbital dependencies we can define the partially screened Coulomb interaction of the single-band model as the arithmetic average of all matrix elements of the partially screened interaction matrix U(q) in the orbital basis. This U sb (q) now represents the Fourier transform of the real space Coulomb interactions U and V as used in Eq. (1) and thus serves as the second important ingredient to our extended Hubbard-Holstein model. Finally, we incorporate the phonon frequency and the electron-phonon coupling into our model to describe all important interactions at the same time. To this end, we employ DFPT 70 calculations as implemented in the Quantum Espresso package 71 using local-density approximations potentials, a lattice constant of a = 3.24 Å, a vacuum height of 16 Å, a k mesh of 32 × 32 × 1 for the self-consistent electronic calculation, and a q mesh of 8 × 8 × 1 for the phonons. Within the BZ, i.e., for increased q momenta, the most important electron-phonon couplings arise due to acoustic phonon modes in NbS 2 (the optical modes couple via a Fröhlich interaction which is proportional to 1/q and is thus strongly decreased here). In more detail, the strongest coupling arises owing to the LA mode, which consequently softens and becomes unstable.
To estimate an average bare frequency for this mode in the monolayer, we make use of the other acoustic branches, which are not at all (ZA) or just slightly (TA) renormalized at the Brillouin zone's M point. Thereby, we arrive at an estimation of ω ph = 20 meV for the bare typical phonon frequency, which is comparable to the corresponding modes in bulk NbS 2 . 32 Using this bare frequency, we re-calculate the renormalized phonon frequency ω re ph ðqÞ = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ω 2 ph þ 2ω ph g 2 χ 0 ðqÞ q using the random phase approximation susceptibility χ 0 (q), where we approximate the phonon selfenergy as g 2 χ 0 ðqÞ. From this, we find instabilities starting from g min ≳ 50 meV and similar instabilities as in the full DFPT calculation for g NbS2 % 60 70 meV. This is comparable to the g max = 0.13 eV found by Flicker and van Wezel 35 for bulk NbSe 2 . Accordingly, we use an interval of g = 0.0 … 0.1 eV in order to study the phonon-induced effects, whereas g ≈ 70 meV is supposed to be our material-realistic estimate. At this point it is important to note, that this is clearly just an approximate model. We neglect the phonon dispersions as well as the momentum dependency of the electron-phonon coupling and focus on a single phonon mode only. The latter is, however, well justified as the LA mode is the only mode becoming unstable in DFPT calculations. Furthermore, between the M and K points, the LA mode is rather flat allowing us to describe it as a local Einstein mode. Finally, the Froehlich-like coupling of the those optical modes, which have a finite coupling to the electrons is likely underestimated around Γ and overestimated around the K point in our model. This means, that a full phonon model might lead to changes in the exact position of the arising divergences in the susceptibilities. In more detail, it is likely that in a full model, the charge instability would emerge more within the Brillouin zone and less at its border.

Dual Boson approach
The resulting material-realistic single-band Hubbard-Holstein model is solved using the Dual Boson (DB) method, which is based on the Dynamical Mean-Field Theory (DMFT) 72 philosophy. That is, DB uses an auxiliary single-site problem to take into account strong correlation effects self-consistently. The DB method extends DMFT by also capturing nonlocal interactions via an effective, dynamic local interaction. Here, as in ref. 73 , the impurity model is determined self-consistently on the Extended DMFT level. Then, the DB method calculates the momentum and frequency resolved susceptibilities starting with a DMFT-like interacting Green's function and then adding non-local vertex corrections (in the ladder approximation) to ensure charge conservation. [45][46][47] The auxiliary impurity model was solved using a modified version of the open source CT-HYB solver 74,75 based on the ALPS libraries. 76 The Dual Boson calculations use the single-band dispersion E mb (k) from the tight-binding model and the effective interaction U sb (q) as their input. Both are evaluated on 144 × 144 × 1 k and q meshes. The electron-phonon coupling leads to an additional, retarded, local electron-electron interaction U eÀph ωn = À2g 2 ωph ω 2 ph þω 2 n , where g is the electron-phonon coupling, ω ph is the phonon frequency and ω n is the nth Matsubara frequency.
Unless otherwise noted, all Dual Boson simulations were performed at β = 25 eV −1 (T = 464 K). Calculations without electron-phonon coupling were for temperatures down to β = 150 eV −1 (77 K). These showed few qualitative changes: the system remained in the strongly correlated phase. For the case of (strong) electron-phonon coupling, closer to the chargeordering transition, the temperature is important as ordering is more likely at low temperature, the same holds for spin ordering transitions. However, please note that the model parameters are derived for T = 0 K.
To learn more about the role of temperature, we show the local charge and spin susceptibility in Fig. 5, the ratio of which is plotted in 3, as a function of g and temperature. Near g ≈ 0.07 eV, the low-temperature simulations approach the phase transition and the charge susceptibility sees a large change, whereas the spin susceptibility develops more smoothly. Comparing the susceptibility with that of the non-interacting system at the same temperature (dashed lines), both temperatures show the same trend, although the magnitude of deviations is generally larger for the low-temperature system. At small g, U is the dominant interaction and the spin susceptibility is enhanced and the charge susceptibility reduced with respect to the non-interacting system. For both temperatures, we find a coupling strength g where both susceptibilities are enhanced compared to the non-interacting system (g ≈ 0.07 eV at the lower temperature and g ≈ 0.85 eV at the higher temperature). Thus, we find this simultaneous enhancement of both susceptibilities to be a general feature that does not require a specific temperature. On the other hand, the magnitude and location of the simultaneous enhancement depends on the temperature. In particular, the sharp rise in the charge susceptibility at low-temperature signals the approach to the charge-order transition, which depends on temperature.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.