Nonlocal Coulomb interaction and spin-freezing crossover as a route to valence-skipping charge order

Multiorbital systems away from global half-filling host intriguing physical properties promoted by Hund’s coupling. Despite increasing awareness of this regime dubbed Hund’s metal, effect of nonlocal interaction is still elusive. Here we study a three-orbital model with 1/3 filling (two electrons per site) including the intersite Coulomb interaction (V). Using the GW plus extended dynamical mean-field theory, the valence-skipping charge order transition is shown to be driven by V. Most interestingly, the instability to this transition is significantly enhanced in the spin-freezing crossover regime, thereby lowering the critical V to the formation of charge order. This behavior is found to be closely related to the population profile of the atomic multiplet states in the spin-freezing regime. In this regime, maximum spin states are dominant in each total charge subspace with substantial amount of one- and three-electron occupations, which leads to almost equal population of one- and the maximum spin three-electron state. Our finding unveils another feature of the Hund’s metal and has potential implications for the broad range of multiorbital systems as well as the recently discovered charge order in iron pnictides.

In addition to the above-mentioned direct manifestations of Hund's metal regime, its connection and proximity to the symmetry-broken charge-disproportionated phases has recently been highlighted 27,28 . Those that are called Hund's insulator 27 and valence-skipping phase 28-30 -a phase with two different valences while skipping the intermediate one between the two-are prominent examples. One possible route to the valence-skipping is the negative effective Coulomb repulsion, U eff < 0 refs 28,31,32 . Interestingly, a purely intra-atomic origin, namely, the anisotropic orbital-multipole scattering, was suggested to be the key ingredient for such valence-skipping phenomena 28 . Furthermore, this phase has potential implications for the electron pairing mechanisms of unconventional superconductivity 28,33 .
The valence-skipping compounds are prevalent in Nature most evidently in the form of charge order (CO) 28 . The CO transition has actively been studied in the single-orbital extended Hubbard model presumably in close connection with the superconductivity of cuprates 34 . Notably, as in the case of cuprates, recent experiments reported the CO in the vicinity of the superconducting phase of AFe 2 As 2 (A = Rb, K, Cs), archetypal materials of Hund's metal [35][36][37] . Moreover, relevance of charge fluctuations or CO to the superconductivity of iron pnictides was reported 38 . Thus it is tempting to presume that the CO is a common "neighbor" of unconventional high-temperature superconductivity. On the other hand, one can also envisage the more complexity of the multiorbital CO transition due to the additional energy scales such as Hund's coupling absent in single-orbital models.
In this work, by employing the state-of-the-art GW plus extended dynamical mean-field theory (GW+EDMFT) adapted to multiorbital models, we demonstrate that the valence-skipping CO is driven by intersite nonlocal Coulomb repulsion V, and the instability to this phase is significantly enhanced in the spinfreezing crossover regime. This enhancement is shown to be related to the local multiplet population profile. This route to the valence-skipping is distinctive from the anisotropic orbitalmultipole scattering mechanism 28 .

RESULTS AND DISCUSSION
We first construct a following model for the two-dimensional square lattice including both local and nonlocal interaction terms: where c y iγσ (c iγσ ) is the electron creation (annihilation) operator acting on site i with orbital index γ = 1, 2, 3 and spin index σ = ↑, ↓. t (t > 0) is the hopping amplitude between two nearestneighbor (NN) sites denoted by 〈ij〉. We use half-bandwidth D = 4t as the unit of energy. n iγσ ¼ c y iγσ c iγσ is the electron number operator. μ is the chemical potential to be adjusted to obey 1/3 filling per site; ∑ γ,σ 〈n iγσ 〉 = 2. H loc is of the Kanamori form containing the onsite Coulomb repulsion U and the Hund's coupling J, which reads H loc ¼ U P i;γ;σ n iγ" n iγ# þ ðU À 2JÞ P γ≠γ 0 i;γ;γ 0 n iγ" n iγ 0 # þ ðU À 3JÞ P γ<γ 0 i;γ;γ 0 ;σ n iγσ n iγ 0 σ ÀJ P γ≠γ 0 i;γ;γ 0 c y iγ" c iγ# c y iγ 0 # c iγ 0 " þ c y iγ" c y iγ# c iγ 0 " c iγ 0 # : H nonloc is the interaction term between two NN sites coupled via nonlocal Coulomb repulsion V, Vn iγσ n jγ 0 σ 0 : To gain a useful insight for CO transition of the model constructed in Eq. (1), we first investigate a simple case of vanishing t and temperature. This simple atomic limit-a limit where the lattice consists of atoms with zero t among themenables us to get analytical solutions, which is found to be a good estimate even under nonzero t and temperature 28,[39][40][41][42] . In Fig. 1, we plot the obtained phase diagram (see Supplementary Note 1 for more details). Three different phases are classified according to their valence. We used notation d N to denote the N-electron occupation of a site in the primitive cell. Note that the triple point emerges at V/U = 0 and J/U = 1/3, which corresponds to the parameter region where the metal resilient to Mott's and Hund's insulator transition emerges 27 , as well as the valence-skipping phases cease to exist 28 . The possible existence of d 3 + d 1 phase was previously noticed from the slave-boson mean field by solving the Kanamori Hamiltonian 27 . This state, however, is degenerate at J/U = 1/3 with d 2 and 2d 3 + d 0 phases and never the ground state unless V/U > 0. The 2d 3 + d 0 phase is equivalent to the chargeordered Hund's insulator 27 . Note also that other COs such as d 4 + d 0 and 2d 0 + d 6 can be stabilized above the dashed lines depicted in Fig. 1, which are quite irrelevant for the present study.
At 0 < J/U < 1/3, we can observe a transition from the isotropic d 2 to d 3 + d 1 valence-skipping CO with ordering wave-vector (π, π) at the critical V (V c ), V c ¼ U 4 ð1 À 3J=UÞ. It should be noted that this phase is driven by V, not by the anisotropic orbital-multipole scattering since the Kanamori form is free from it by construction 28 . At J/U = 0, V c follows the half-filled single-orbital result of V c = U/4 refs 39,40 .
With insight obtained above, we now turn to our GW+DMFT results. The corresponding phase diagram obtained from GW+EDMFT is shown in Fig. 2a. We identified the CO transition by monitoring the divergence of the static charge susceptibility χ(k, iν 0 ) (ν 0 is the lowest bosonic Matsubara frequency, ν 0 = 0). The divergence actually occurs at the wave vector k = (π, π) indicating the formation of d 3 + d 1 order (see Figs. 1 and 3b). One can also confirm that this CO transition is driven by V (compare Fig. 3b with Fig. 3a; see also Supplementary Note 3 for χ(k, iν 0 ) at V = 0).
The actual GW+EDMFT results roughly follow the atomic limit estimate at J/U ≤ 0.15 and are in fair agreement at large U(U = 3, 4) and J/U = 0.15. Even in smaller U region (Fermi liquid (FL); see Fig.  4a), GW+EDMFT results qualitatively follow the atomic limit estimates. This seemingly unusual behavior is attributed to the leading contribution of interaction energy compared to the kinetic energy in determining the CO transition boundary 43 . Note that the Mott phase emerges at V = 0 for U = 4 (when J/U = 0.05) and U = 5 (when J/U ≤ 0.15). In the current study, we restrict our discussion to the U and J/U region in which the metallic phase is obtained when V = 0.
At J/U = 0.2, GW+EDMFT results exhibit unprecedented behavior at large U (U ≥ 3): CO instability is significantly enhanced, thereby pushing V c further below the atomic limit estimates. Notably, the downturn of V c is most pronounced at U = 4 followed by a rapid upturn of the phase boundary at U = 5. This behavior is not captured either by EDMFT or GW approximations (see Fig. 2b, c). On the other hand, at smaller U (U ≤ 2), V c values obtained from GW+EDMFT are almost identical to those of EDMFT and larger than the atomic limit estimates. The discrepancy between GW results (Fig. 2c) and the others is reasonable since this method cannot properly treat the local physics. We briefly remark that, at further larger J/U, especially near J/U = 1/3 and V = 0 in which the triple point emerges in Fig. 1, a signature of the 2d 3 + d 0 phase is also expected. Notably, the presence of this degeneracy point is claimed to play an important role in stabilizing the metallic phase 27 . The triple degeneracy, however, should be lifted by nonzero V. We expect that an intriguing physics can happen due to this broken degeneracy, which we leave for future study.
To further illustrate the above intriguing result from GW+ EDMFT at large U and J/U regime, we investigate the site-resolved charge susceptibility χðR i ; iν 0 Þ ¼ R dke ikÁRi χðk; iν 0 Þ (R i is the position vector of the ith NN). The magnitude of this quantity is enhanced as V increases as shown in Fig. 3c-f. Near the CO boundary (V ≃ 0.9V c ), the sign of χ(R i , iν 0 ) clearly indicates the CO instability at k = (π, π), which has to be plus (minus) for onsite, second, and third (first and fourth) NNs.
Most interestingly, the large U results exhibit the rapid growth of χ(R i , iν 0 ) as a function of J/U at a finite V (see Fig. 3f). This behavior is in contrast to the smaller U results in which the increase of χ(R i , iν 0 ) is much more gradual (see Fig. 3d). This enhancement of χ(R i , iν 0 ) at J/U = 0.2 is further manifested by the static effective local interaction, Uðiν 0 Þ. The intraorbital elements Uðiν 0 Þ γγ Uðiν 0 Þ γγγγ at U = 4 and V ≃ 0.9V c shows the large screening effect at J/U = 0.2; from Uðiν 0 Þ γγ ¼ 3:31 (3.39) at J/U = 0.1 (0.15) to Uðiν 0 Þ γγ ¼ 2:88 at J/U = 0.2. Note also that the substantial amount of nonlocal χ(R i , iν 0 ) exists even at V = 0 in the larger U and J/U = 0.2 regime (compare Fig. 3c, e and their insets).
Key information for understanding the large enhancement of CO instability is provided by investigating the local self-energy Σ loc (iω n ) (ω n : fermionic Matsubara frequency). Σ loc (iω n ) shows an interesting behavior near spin-freezing crossover regime 3,8 , which is a metal with emerging local moment: large spin susceptibility 25 . S i = (1/2)∑ γ (n iγ↑ − n iγ↓ ) is the local spin operator. In this regime, ImΣ loc ðiω n Þ is claimed to follow the power-law behavior at low frequency: ImΣ loc ðiω n Þ ' ÀΓ þ Aðω n Þ α with α ≃ 0.5 and Γ ≃ 0 8 . Deep inside this crossover where non-FL behavior appears (Γ > 0 and α > 0.5) is called the frozen-moment regime 3,8,25 . In Fig. 4, we summarize our analysis of ImΣ loc ðiω n Þ. Figure 4a shows the correlation between α and Δχ s /χ s . By construction, Δχ s /χ s lies in between 0 and 1. The limiting value of Δχ s /χ s indicates either the FL limit when Δχ s /χ s → 1 or the frozenmoment regime when Δχ s /χ s → 0. Thus we can naturally expect that the spin-freezing regime should lie somewhere in between these two limits. We identify the region of 0.4 ≲ α ≲ 0.5 and Γ ≃ 0 with the spin-freezing crossover regime. In our parameter range, spin-freezing regime appears for 0.25 < Δχ s /χ s < 0.4 (see also Supplementary Note 4 for the correlation of α with χ À1 s and Δχ s ). In FL regime, increasing V drives the system to be less correlated. Interestingly at U = 4 and J/U = 0.1, V drives the system from the (proximity of) frozen moment to FL and eventually to CO. This behavior can be confirmed by vanishing Γ and α > 0.5 near V c (see Fig. 4b). We also note that EDMFT yields qualitatively similar results except that, at U = 4 and J/U = 0.1, increasing V do not show any signal of transition to the FL.
Most notably, the parameter region showing the unusual downturn of V c (U = 3, 4 and J/U = 0.2) corresponds to the spinfreezing crossover regime. As U increases further at J/U = 0.2, an upturn of the phase boundary appears (see Fig. 2a) as entering deeper into the frozen moment regime. In this range of U and J/U, the increasing V tends to reduce α while maintaining Γ = 0 (see Fig. 4a, c). To further clarify the relation between the enhanced CO instability and the spin-freezing crossover, we investigate the local populations (or probabilities) of atomic multiplet states. The U (1) charge × SU(2) spin × SO(3) orbital symmetry of Eq. (2) allows us to have the simultaneous eigenstates of charge N, orbital L, and spin S as N; L; S j i 1,3,27 . The local population profiles of these eigenstates are plotted in Fig. 5a, b as approaching the CO boundary.
One can notice that, in spin-freezing crossover regime, maximum S states are dominant in each total charge subspace with substantial amount of N = 1 and N = 3 populations (contribution of states other than N = 1, 2, 3 subspaces are negligible) (Fig. 5c). It is the effect of J favoring the maximum S. Importantly, these N = 1 and N = 3 charges are directly related to the d 3 + d 1 CO phase, which implies the enhanced CO instability in this regime. On the other hand, in the frozen moment regime, N = 2 population is more dominant with reduced N = 1 and N = 3 portions than the spin-freezing case; compare Fig. 5a, b with Fig.  5c. FL regime exhibits non-negligible excursions to every other N; L; S j i as expected. We hereafter denote the population of 1; 1; 1=2 j iand 3; 0; 3=2 j iby p 1 and p 3 . At U = 4 and J/U = 0.2, only the maximum S is selected in the N = 3 subspace, leading to p 3 ≃ p 1 (see Fig. 5c). At U = 5 and J/U = 0.2 (frozen moment), p 3 ≃ p 1 is also found. This case, however, shows more dominant N = 2 population (~0.79 at V = 0) with reduced N = 1 and N = 3 contributions compared to the U = 4 and J/U = 0.2 case. In light of this observation, we construct a phenomenological local wave-function ψ consisting of maximum S states, namely, ψ ¼ ffiffiffiffi ffi p 1 p 1; 1; 1=2 j i þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 À 2p 1 p 2; 1; 1 j i þ ffiffiffiffi ffi p 1 p 3; 0; 3=2 j i apart from the phase factor, which is irrelevant for the evaluation of energy. The re-calculated V c estimate (as is done for Fig. 1) by means of ψ is shown in Fig. 5d. We can observe the qualitative agreement with the actual behavior obtained from GW+EDMFT at J/U = 0.2 (see stars in Fig. 5d). This result confirms the role of maximum S states in N = 3 subspace in enhancing the CO instability. This type of interpretation should be valid in large U and J/U limit. Figure 5d shows, however, deviations of actual GW+EDMFT results at J/U = 0.1 and J/U = 0.15. This can be attributed to the non-negligible amount of smaller S states in N = 3 subspace and the fundamental inadequacy of this kind of approach for the FL regime.
In conclusion, we have shown by employing GW+EDMFT that, in the spin-freezing regime, significant enhancement of d 3 + d 1 CO instability appears. This enhancement is found to be closely related to the local multiplet population profile: maximum spin states are dominant in each total charge subspace with substantial amount of N = 1 and N = 3 occupations. The observed d 3 + d 1 CO  transition is driven by V and is also a distinctive route from the anisotropic orbital-multipole scattering mechanism to the valence-skipping phenomena 28 . Our study unveils another feature of the Hund's metal and has potential implications for other multiorbital systems and observed CO in Hund's metal AFe 2 As 2 (A = Rb, K, Cs) [35][36][37] .

METHODS
GW+EDMFT is derivable from the Ψ[G, W] functional (G: Green's function, W: fully screened Coulomb interaction) 44 as Ψ GWþEDMFT ½G; W ¼ Ψ EDMFT ½G loc ; W loc þ Ψ GW nonloc ½G; W, where EDMFT is supplemented with nonlocal GW functional [45][46][47][48][49] . This approach allows a nonperturbative solution of the auxiliary impurity model with self-consistently determined local fermionic and bosonic Weiss fields. The bosonic Weiss field Uðiν n Þ (ν n : bosonic Matsubara frequency) is the effective impurity interaction whose value is renormalized by dynamical screening effect. The importance of this effect has recently been highlighted [49][50][51][52][53][54][55][56] . We performed calculations within the paramagnetic isotropic phase and inverse temperature of βD = 100. An impurity model was solved using the COMCTQMC implementation 57 of the hybridization-expansion CTQMC algorithm 58,59 . Both local and nonlocal interaction terms were decoupled via Hubbard-Stratonovich transformation to treat them on an equal footing 49 . In our current implementation, owing to the computational complexity, we measured only the density-density type of two-particle correlation functions from the impurity; χ imp ðτÞ ¼ hT τ n γσ ðτÞn γ 0 σ 0 ð0Þi (τ: imaginary time). The nondensity-density-type functions are responsible for the screening of nonmonopole terms of charge distribution, making our approximation physically reasonable since these terms are ill-screened 55 . All three methods (GW+EDMFT, EDMFT, and GW) were performed selfconsistently. See Supplementary Note 2 for further details of our GW +EDMFT calculations.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding authors upon reasonable request.