Nonlocal Coulomb Interaction and Spin Freezing Crossover : 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 oneand three-electron occupations, which leads to almost equal population of oneand the maximum spin three-electron state. Our finding unveils another feature of the Hund’s metal, and has potential implications to the broad range of multiorbital systems as well as the recently discovered charge order in iron-pnictides.

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 to 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 [23,24]. Those which are called Hund's insulator [23] and valenceskipping phase [24][25][26]] -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 [24,27,28]. Interestingly, a purely intraatomic origin, namely the anisotropic orbital-multipole scattering, was suggested to be the key ingredient for such valence-skipping phenomena [24]. Furthermore, this phase has potential implications to the electron pairing mechanisms of unconventional superconductivity [24,29].
The valence-skipping compounds are prevalent in Na-ture most evidently in the form of charge order (CO) [24]. The CO transition has actively been studied in the single-orbital extended Hubbard model presumably in close connection with the superconductivity of cuprates [30]. 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 [31][32][33]. Moreover, relevance of charge fluctuations or CO to the superconductivity of iron-pnictides was reported [34]. 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 multorbital CO transition due to the additional energy scales such as Hund's coupling absent in singleorbital models.
In this Letter, by employing the state-of-theart 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 orbital-multipole scattering mechanism [24].
We first construct a following model for the twodimensional square lattice including both local and nonlocal interaction terms: where c † 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 nearest-neighbor (NN) sites denoted by ij . We use half-bandwidth D = 4t as the unit of energies. n iγσ = c † 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 nonloc is the interaction term between two NN sites coupled via nonlocal Coulomb repulsion V , 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 them -enables us to get analytical solutions, which is found to be a good estimate even under nonzero t and temperature [24,35,36]. In Fig. 1 we plot the obtained phase diagram. 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-and Hund'sinsulator transition emerges [23], as well as the valenceskipping phases cease to exist [24]. The possible existence of d 3 + d 1 phase was previously noticed from the slaveboson mean-field by solving the Kanamori Hamiltonian [23]. 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 charge-ordered Hund's insulator [23].
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 ), 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 [24]. At J/U = 0, V c follows the half-filled single-orbital result of V c = U/4 [35,36].
With insight obtained above, we now turn to our GW +EDMFT results. GW +EDMFT is derivable from the Ψ[G, W ] functional (G: Green's function, W : fully screened Coulomb interaction) [37] where EDMFT is supplemented with nonlocal GW functional [38][39][40][41][42]. 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 [42][43][44][45][46][47][48][49]. We performed calculations within the paramagnetic isotropic phase and inverse temperature of βD = 100 [50]. An impurity model was solved using the ComCTQMC implementation [51] of the hybridization-expansion CTQMC algorithm [52,53]. Both local and nonlocal interaction terms were decoupled via Hubbard-Stratonovich transformation to treat them on an equal footing [42]. In our current implementation, due to the computational complexity we measured only the density-density type of two-particle correlation functions from the impurity; χ imp (τ ) = T τ n γσ (τ )n γ σ (0) (τ : imaginary time). The non-density-density type functions are responsible for the screening of non-monopole terms of charge distribution, making our approximation physically reasonable since these terms are ill-screened [48].
The corresponding phase diagram obtained from GW +EDMFT is shown in Fig. 2(a). 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 Fig. 1 and Fig. 3(b)). One can also confirm that this CO transition is driven by V (compare Fig. 3 Fig. 3(a)) [50].
The actual GW +EDMFT results roughly follow the Mott phase emerges at U = 4, and thus pushing the V c significantly upwards. This behavior is reminiscent of single-orbital model results at half-filling [46]. In the current study, we restrict our discussion to the metallic phase. Even in smaller U regime (Fermi-liquid; see Fig. 4(a)), 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 [54]. 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. This behavior is not captured either by EDMFT or GW approximations (see Fig. 2(b) and (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 result (Fig. 2(c)) and the others is reasonable since this method cannot properly treat the local physics.
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 ) = dke ik·Ri χ(k, iν 0 ) (R i is the position vector of the i-th NN). The magnitude of this quantity is enhanced as V increases as shown in Fig. 3(c) -(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. 3(f)). 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. 3(d)). 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. 3(c) and (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 χ s = β 0 dτ S i (τ )S i (0) with substantial dynamical contribution of ∆χ s = β 0 dτ S i (τ )S i (0) − S i (β/2)S i (0) [22]. 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. Deep inside this crossover where non-FL behavior appears (Γ > 0 and α > 0.5) is called the frozen moment regime [3,8,22]. In Fig. 4, we summarize our analysis of ImΣ loc (iω n ). Note that we used three lowest Matsubara frequencies in fitting ImΣ loc (iωn) to obtain α and Γ. Fig. 4(a) shows the correlation between α and ∆χ s /χ s . Here we identified 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. In FL, ∆χ s /χ s is large since the most part of spin susceptibility stems from the dynamical fluctuations. We can identify a broad region of FL behavior for most of our parameter range. 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. 4(b)). We also note that EDMFT yields qualitatively similar results (not shown) except that in 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 spin-freezing crossover regime. In this range of U and J/U , the increasing V tends to reduce α while maintaining Γ = 0 (see Fig. 4(a) and (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 [1,3,23]. The local population profiles of these eigenstates are plotted in Fig. 5(a) and (b) as approaching the CO boundary.  Fig. 2(a)).
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. 5(c)). 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. 5(a) and (b) with Fig. 5(c). FL regime exhibits non-negligible excursions to every other |N, L, S as expected. We hereafter denote the population of |1, 1, 1/2 and |3, 0, 3/2 by 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. 5(c)). In light of this observation, we construct a phenomenological local wave-function ψ consisting of maximum S states, namely |1, 1, 1/2 , |2, 1, 1 , and |3, 0, 3/2 with corresponding probability p 1 , 1 − 2p 1 , and p 1 , respectively. The re-calculated V c estimate (as is done for Fig. 1) by means of ψ is shown in Fig. 5(d). We can observe the qualitative agreement with the actual behavior obtained from GW +EDMFT at J/U = 0.2 (see stars in Fig. 5(d)). 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. Fig. 5(d) shows, however, deviations of actual GW +EDMFT results at J/U = 0.1 and J/U = 0.15. This can be attributed to the nonnegligible amount smaller S states in N = 3 subspace, and the fundamental inadequacy of this kind of approach for FL regime.
In conclusion, we have shown by employing GW +EDMFT that in 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 orbitalmultipole scattering mechanism to the valence-skipping phenomena [24]. Our study unveils another feature of the Hund's metal, and has potential implications to other multiorbital systems and observed CO in Hund's metal Supplemental material for "Nonlocal Coulomb Interaction and Spin Freezing Crossover: A Route to Valence-skipping Charge Order" Computation detail for GW +EDMFT calculations The impurity action for GW +EDMFT is given by: where c † a (c a ) is the electron creation (annihilation) operator for a composite index a labeling both orbital and spin. G and U is the fermionic and bosonic Weiss field. The latter is the effective impurity interaction which incorporates the dynamical screening effect. In our case, U in matsubara frequency space has the form of which takes into account the dynamical screening D(iν n ) in the density-density terms. U abcd is of the Kanamori form. G and U are determined self-consistently in the spirit of EDMFT using impurity self-energy (Σ imp ) and polarizability (P imp ): where G loc and W loc is the local Green's function and the fully screened Coulomb interaction obtained from Dyson's equations: G 0 is the noninteracting Green's function, and V (k) is the bare Coulomb interaction in k-space. N k is the number of k-points in the first Brillouin zone. In GW +EDMFT, Σ and P read: In our case, χ imp is of the density-density form and so is the resulting P imp . The density-density component of the local GW polarizability, P GW loc (iν n ) = 1 N k k P GW abcd (k, iν n )δ ab δ cd , is subtracted to avoid the double-counting (Eq. (8)). We performed calculations using 32 × 32 k-points in the first Brillouin zone of the square lattice and β = 100. The so-called space-time method [55,56] was applied to calculate the fermionic and bosonic GW quantities, which requires one to evaluate the summation over infinite matsubara frequencies. To this end, we divided the summation interval into two: exact and asymptotic. In the exact interval, we performed the summation exactly. In the asymptotic interval, we replaced G(k, iω n ) and W (k, iν n ) with their high-frequency analytic forms. The cutoff matsubara frequency (ω cutoff ) above which the asymptotic forms of G and W are taken was set to be U -dependent in order to precisely encompass the high-frequency features of nonlocal GW quantities. We used ω cutoff 40, 60, 80, and 100 for U = 1, 2, 3, and 4, respectively. To ensure a stable convergence, we linearly mixed polarizability between iterations: P (k, iν n ) = (1 − R mix )P (k, iν n ) old + R mix P (k, iν n ) new with R mix = 0.1 -0.3 for most cases. We took converged EDMFT solutions as starting points for GW +EDMFT calculations.