Effective Heisenberg model and exchange interaction for strongly correlated systems
E. A. Stepanov S. Brener F. Krien M. Harland A. I. Lichtenstein M. I. Katsnelson




Radboud University, Institute for Molecules and Materials, 6525AJ Nijmegen, The Netherlands



Institute of Theoretical Physics, University of Hamburg, 20355 Hamburg, Germany


Institute of Theoretical Physics, University of Hamburg, 20355 Hamburg, Germany


Institute of Theoretical Physics, University of Hamburg, 20355 Hamburg, Germany


Institute of Theoretical Physics, University of Hamburg, 20355 Hamburg, Germany



Radboud University, Institute for Molecules and Materials, 6525AJ Nijmegen, The Netherlands





We consider the extended Hubbard model and introduce a corresponding Heisenberg-like problem written in terms of spin operators. The derived formalism is reminiscent of Anderson's idea of the effective exchange interaction and takes into account nonlocal correlation effects. The results for the exchange interaction and spin susceptibility in the magnetic phase are expressed in terms of single-particle quantities. This fact can be used not only for realistic calculations of multiband systems, but also allows to reconsider a general description of many-body effects in the most interesting physical regimes where the physical properties of the system are dominated by collective (bosonic) fluctuations. In the strongly spin-polarized limit, when the local magnetic moment is well-defined, the exchange interaction reduces to a standard expression of the density functional theory that has been successfully used in practical calculations of magnetic properties of real materials. 



Effective Heisenberg model and exchange interaction for strongly correlated systems
    M. I. Katsnelson
    December 30, 2023
====================================================================================

The theory of magnetism is one of the most attractive and discussed areas of physics. An additional interest to this topic is heated up by the theoretical prediction <cit.> and experimental observation <cit.> of topologically stable skyrmionic spin textures that are intensively studied nowadays in the context of spintronics and magnetic data storing <cit.>. Also, a correct account for spin excitations is important for realization of Kitaev spin model <cit.> and its practical application in the Majorana quantum computers <cit.>. A quantitative description of the mentioned effects requires the knowledge of the exchange interaction between two spins. However, this problem is challenging when applied to many magnetic materials that are by definition strongly correlated quantum systems. 

Originally, the development of the theory of exchange interactions in solids and molecules was based on the Heitler-London theory of the hydrogen molecule <cit.>. It has been demonstrated, however, in the early 60's by Freeman and Watson <cit.> that this theory, being applied to ferromagnetic transition metals, gives a completely wrong order of magnitude and even an incorrect sign of the exchange parameters. For magnetic insulators, a semi-empirical theory of exchange interactions has been developed in the 50's, known as Goodenough-Kanamori-Anderson rules <cit.>, however, it was not quantitative. An analysis of “superexchange” in particular compounds always assumed some model considerations, that is, the importance and non-importance of specific intermediate states. When the density functional theory (DFT) became the base of microscopic quantum theory of molecules and crystals <cit.> the most straightforward way to estimate the exchange interactions was simply the calculation of the total energy difference between ferromagnetic and antiferromagnetic phases. This assumes the applicability of the Heisenberg model, which is frequently not the case, especially for itinerant electron systems <cit.>.

A general, model-independent and parameter-free method to calculate exchange interactions within DFT was suggested in Refs. <cit.> based on the “magnetic local force theorem”. It is based on the consideration of second-order variations of the total energy with respect to small rotations of magnetic moments starting from equilibrium  ground states. Later this approach was generalized to strongly correlated systems <cit.> (within the framework of dynamical mean-field theory (DMFT) <cit.>), magnetic systems out of equilibrium <cit.>, and relativistic magnetic interactions, such as the Dzyaloshinskii-Moriya interaction <cit.>. This theory was successfully used for many calculations of real systems, such as magnetic semiconductors <cit.>, molecular magnets <cit.>, ferromagnetic transition metals <cit.> and half-metallic ferromagnets <cit.>.

Despite the success of this approach its conceptual status remains unclear. Indeed, a mapping from DFT or from a  Hubbard model to the Heisenberg model is in general impossible; exchange interactions obtained from the magnetic force theorem are classical and dependent on the magnetic configuration (see, e.g. <cit.>). Their relation to observables is not very clear; strictly speaking, only the spin-wave stiffness constant in ferromagnets is a well-defined quantity since we can be sure that in the limit of slow times and large spatial scales the phenomenological Landau-Lifshitz equations are correct. This was emphasized already in the first paper <cit.>. 
Observables are directly related to the dynamic magnetic susceptibility, but to establish relations between the magnetic local force approach and the standard language of response functions is not an easy problem. It was solved only within the local spin-density approximation in DFT <cit.> and within the time-dependent mean-field approach in the Hubbard model <cit.>. However, most of the interesting magnetic materials are strongly correlated systems, and these approximations seem to be insufficient (or, at least, not completely justified) to describe spin dynamics.

In this Letter we show that the extended Hubbard Hamiltonian can be mapped onto an effective Heisenberg model. Inspired by the Dual Boson (DB) formalism <cit.> we construct a bosonic model, whose interaction is reminiscent of Anderson's superexchange mechanism <cit.>. Importantly, the derived formalism remains applicable not only in the strongly localized regime and allows the description of every magnetic system with a well-defined local magnetic moment. Moreover, the presence of the latter
allows to reveal a general way of the description of a complicated quantum many-body problem in terms of single-particle quantities with the use of Ward identities <cit.>.

Effective s-d model — 
We consider the action of the extended Hubbard model for correlated electrons,

    S =    -∑_,ν,σ c^*_νσ[iν+μ-ε^†_]
      c^*_νσ
       + U∑_,ω n^*_ω↑ n^*_ω↓ +1/2∑_,ω,ςρ^* ς_ω[V^*_]_ςςρ^ ς_ω.

Here c^*_νσ (c_νσ) are Grassmann variables corresponding to creation (annihilation) of an electron with momentum ,  fermionic Matsubara frequency ν and spin σ labels.
The label ς={c,  s} depicts charge c and spin s = {x,y,z} degrees of freedom, so that U corresponds to local Coulomb interaction, [V_]_cc=V_ and [V_]_ss=-J^ d_/2 describe nonlocal Coulomb and direct ferromagnetic exchange interactions, respectively. Here, we also introduce bosonic variables:  ρ^ ς_ω=n^ ς_ω-⟨ n^ ς_ω⟩, where n^ς_ω = ∑_νσσ'c^*_νσσ^ς_σσ' c^*_+,ν+ω,σ' is the charge (ς=c) and spin (ς=s) density of electrons with the momentum , bosonic frequency ω and Pauli matrices σ^ς={1,σ^s}.  

Expressing the effective exchange interaction in terms of correlation functions is a nontrivial task, since it is not an observable. Furthermore, in the strongly correlated regime charge and spin fluctuations are entangled in a complicated way.
Both challenges can be approached within the Dual Boson formalism <cit.>, since it naturally separates charge and spin degrees of freedom by representing them in terms of bosonic fields entering an effective action. To this aim one splits the lattice action (<ref>) into the local impurity problem 

of the extended dynamical mean-field theory (EDMFT, <cit.>) and the remaining non-local part, 

which is a bilinear function of c^*(c) and ρ variables. 
Within the DB approach this remaining part is decoupled by two Hubbard-Stratonovich transformations, thus introducing dual fermionic f^* (f) and bosonic ϕ fields. Then, the initial fermionic degrees of freedom c^* (c) can be integrated out, leading to the interaction part W̃[f,ϕ] of the resulting dual action being expressed in terms of the full vertex functions of the local impurity problem (for details see Suppl. Mat. <cit.>). Thus, by construction, local correlations are already embedded into the bare propagators and interactions of the DB problem, which is very convenient for practical calculations. In the following we restrict ourselves to the lowest order terms in W̃[f,ϕ] stemming from the four-point γ_νν'ω and three-point γ_νω vertices <cit.>.

Dual fields f^* (f) and ϕ have no direct physical interpretation, but this fact does not represent a significant obstacle for the calculation of physical observables, since there is an exact connection between dual and lattice quantities <cit.>. However, for our goal of deriving an effective bosonic model that describes initial (lattice) degrees of freedom it is crucial to formulate the problem in terms of bosonic fields that have a clear physical meaning. To remedy this problem, we perform the reverse Hubbard-Stratonovich transformation for the bosonic variables ϕ introducing fields ρ̅. In this we were inspired by works of Dupuis <cit.>, where a similar trick was performed for fermionic degrees of freedom. After integrating over dual bosonic fields ϕ one gets the following action reminiscent of the s-d model <cit.>

    S_s-d =    -∑_,ν,σ f^*_νσG̃^-1_0f^*_νσ
    -1/2∑_,ω,ς(')ρ̅^* ς_ω[X^*_ E]^-1_ςς'ρ̅^ ς'_ω + W.

Here, X_ E is the EDMFT susceptibility and G̃_0 is the nonlocal part of the EDMFT Green's function. Importantly, after all transformations the field ρ̅ indeed has the same physical meaning as original composite bosonic field ρ of the lattice problem (<ref>) as shown in <cit.>. The decisive advantage of the variable ρ̅ is that it can now be treated as the elementary bosonic field that has a well-defined propagator and is independent of fermionic degrees of freedom c^* (c). Remarkably, W[f,ρ̅] keeps the practical form of the dual interaction W̃[f,ϕ] with the replacement of bosonic variable ϕ→ρ̅, although the four-fermionic term is modified under these transformations. As we argue in <cit.> and numerically check below, in the case of well-developed bosonic fluctuations this modification results in the corresponding contribution to the interaction W[f,ρ] becoming negligibly small and the latter takes the simple form 
W[f,ρ]
≃∑_,∑_ν,ω,ςγ^ ς_νω ρ^* ς_ω 
 f^*_νσ f^*_+,ν+ω,σ'. 
At last we mention, that the fermionic degrees of freedom are kept in the dual space, which will prove to be useful to discriminate between local and nonlocal contributions to the lattice susceptibility.


Magnetic susceptibility —
In order to design an effective Heisenberg model for spin degrees of freedom, one has to assume that the local magnetization m=2S^z is described well at the dynamical mean-field level and fluctuations revealed by the system beyond EDMFT are mostly bosonic. In order to have well-defined local magnetic moment, the effective impurity model has to be considered for the spin polarized state. For easier description, one can transform spin variables from s={x,y,z} to s={+,-,z} basis with S^±=(ρ^x±iρ^y)/2. In the spin-polarized case charge and spin z channels are yet entangled, but the ± spin channel can be separated in the collinear case <cit.>. Thus, for the correct description of the spin fluctuations, one may consider correlations only in the ± spin channel and the contribution of the z channel to the exchange interaction can be later restored from symmetry arguments. For simplicity, ± spin labels are omitted wherever they are not crucial for understanding.  

Now, one can integrate out fermionic degrees of freedom in the effective action (<ref>) and get the following spin model

    S_ spin   =-1/2∑_,ω S^ -_ω[X_ω^ -+]^-1  S^ +_-,-ω +  h.c.

A first approximation for the magnetic susceptibility X_ω can be obtained for the case when the interaction W[f,ρ̅] contains only the three-point vertex γ^±_νω, as discussed above. Therefore,
the expansion of the partition function of the action (<ref>) up to the second order with respect to bosonic fields gives <cit.>

    [X^(2)_ω]^-1   = J^ d_ q + Λ^1_ω + χ^-1_ω - Π̃^(2)_ω.

Here, Λ_ω and χ_ω are the bosonic hybridization function and susceptibility of the impurity problem, respectively. Also,

    Π̃^(2)_ω=∑_,νγ^ -_ν+ω,-ω G̃^2_+,ν+ω↑G̃^2_ν↓ γ^ +_ν,ω=
    
        < g r a p h i c s >

is the second order polarization function <cit.>. 
Note that a conserving description of spin fluctuations is given by the two-particle ladder approximation of the magnetic susceptibility provided by the ladder DB approach <cit.> that accounts for the four-fermionic contribution in W[f,ρ̅] and treats bosonic hybridization Λ as a constant <cit.>

    [X^ ladd_ω]^-1 = J^ d_ + Λ + [X^ DMFT_ω]^-1.

Here, X^ DMFT_ω=χ^2_ω+χ^2_ωΠ̃^ ladd_ωχ^2_ω is the DMFT- <cit.>, or DΓA-like <cit.> susceptibility written in terms of local two-particle irreducible four-point vertices and lattice Green's functions. Π̃^ ladd_ω is the dual polarization in the ladder form <cit.> that contains Π̃^(2)_ω as the lowest order term. 
Therefore, the hybridization Λ plays the role of the Moriyaesque λ correction that was introduced in DΓA <cit.> by hand similarly to the Moriya and Kawabata theory of weak itinerant magnets <cit.> and now is derived analytically.



Importantly, the expressions for the magnetic susceptibility (<ref>) and (<ref>) can be drastically simplified to be applicable for realistic multiband calculations, for which the two-particle quantities can hardly be obtained. As it was discussed above, the system with a well-defined local magnetic moment exhibits mostly bosonic fluctuations. Therefore, one can expect that local vertex functions are mostly described by the bosonic frequency ω, while the dependence on fermionic frequencies ν,ν' is negligible and can be averaged out. In order to perform this averaging consistently, it is carried out using the local Ward identities <cit.>, which leads to the following approximation of three-point vertex <cit.>

    γ^+_νω = γ^-_ν+ω,-ω≃χ^-1_ω + δΣ^1_νω ≃χ^0 -1_ω

Here, χ^0_ω=∑_νg_ν+ω↑g_ν↓ is the bare spin susceptibility, g_νσ and Σ_νσ are the full Green's function and self-energy of the impurity problem and δΣ_νω = (Σ_ν+ω↑-Σ_ν↓)/m. Therefore, exploiting the system being in the magnetic phase allows to rewrite the complicated many-body problem (<ref>) in a much simpler form of Eq. <ref> introducing bosonic fields that correspond to the collective magnetic fluctuations. In this case, the expression for the corresponding fermion-boson coupling γ^±_νω can be in drastically simplified (<ref>), leading to a similar expression that was recently postulated in <cit.> and numerically checked using brute force calculations <cit.>.

Exact numerical solution — 
In order to exemplify the above approximations we consider the half-filled Hubbard model (<ref>) (V_, J^ d_, Λ=0) on the hypercubic lattice in infinite dimensions. In this case, the exact result for the magnetic susceptibility is known to be given by the DMFT expression (<ref>) and can be compared to the simplified result of Eq. <ref>. At low temperatures this system favors antiferromagnetic (AFM) order over paramagnetism as shown in the phase diagram in Fig. <ref>.

 

The local four-point vertex γ_νν'ω is measured at U=5 for the three temperatures marked in Fig. <ref>, roughly below the maximum of the AFM dome, where T_N≈0.186 is obtained using DMFT <cit.>. As the temperature is lowered from β=6 to 10, the magnetization increases from m≃0.42 to 0.84. We validate in Fig. <ref> that at large magnetization the dependence of the four-point vertex γ_νν'ω on fermionic frequencies ν,ν' is small. Consequently, one may indeed use the approximated form of the vertex γ_νν'ω≃γ_νν'ω, which leads to Eq. <ref>.

We evaluate Eq. (<ref>) in the AFM phase at the =0 point of the reduced Brillouin zone. The transversal susceptibility is a 2×2 matrix with the homogeneous susceptibility X^hom(ω) as a diagonal element <cit.>. Fig. <ref> shows X^hom(ω), which is real, as well as the off-diagonal element X^ off(ω). Remarkably, despite the approximation of the vertex functions, X^ hom(ω≠0)=0 and X^ off(ω≠0)=-2im/ω hold to very good accuracy, which are exact constraints due to global spin conservation <cit.>. 

At U=5 the eigenvalue of the ladder Eq. (<ref>) corresponding to X^hom(ω=0)
is large (≃0.715). Therefore, one can not approximate the polarization Π̃^ ladd_ω by the second order expression Π̃^(2)_ω in Eq. <ref>. The corresponding approximation for X^hom(ω=0) and X^off(ω=2πβ) is marked in Fig. <ref> with open triangles and indeed clearly distinguishable from Eq. (<ref>).

Nonetheless, the simplified expression for magnetic susceptibility X^(2) (<ref>) with the vertex approximation (<ref>) shows a good agreement with X^ ladd (<ref>).
Importantly, the approximation for the magnetic susceptibility obtained in Eq. <ref> should not be confused with the truncation of the ladder equation, even though it formally uses the same quantity Π̃^(2)_ω. The good agreement of the simplified result X^(2) with the much more advanced ladder approximation (<ref>) shows that the bosonic fluctuations indeed dominate in the polarized regime of the impurity model, which was assumed while deriving Eq. <ref>. 



Classical Heisenberg Hamiltonian — Although the action (<ref>) is general and can be used for the description of quantum effects in terms of susceptibilities, at low temperatures it can be mapped onto an effective classical Heisenberg Hamiltonian H_ spin = - ∑_ J_  S_  S_- that describes small spin fluctuations around the AFM ground state <cit.>. To this aim, spin variables S^±_ω in (<ref>) are replaced by classical vectors S_ of the length ⟨ S_z⟩ and the contribution from the z spin channel is restored from the requirement of rotational invariance. Then, an effective exchange interaction J_ can be defined as a nonlocal part of the inverse spin susceptibility at the zero bosonic frequency <cit.>. 
Thus, the effective exchange interaction that corresponds to the simplified form of magnetic susceptibility (<ref>) reads

    J^d_ = J^ d_ - ∑_,νγ^ -_ν,ω=0 G̃^2_+,ν↑G̃^2_ν↓ γ^ +_ν,ω=0,

while the exchange interaction in the ladder approximation is detailed in <cit.>. 
This result reminds of Anderson's idea of the superexchange interaction <cit.>. Indeed, the first and the second term in Eq. <ref> describe the direct ferromagnetic and kinetic antiferromagnetic exchange interactions, respectively. As a result, in the strongly localized regime and in the case of antiferromagnetic dimer the kinetic part of the exchange interaction takes the well-known form J=-2t^2/U <cit.>.

It is worth mentioning that the three-point vertex γ_ν,ω that enters the kinetic part of the exchange interaction describes the total spin splitting. In the spin polarized case one can again use the simplified form of the vertex function (first approximation in Eq. <ref>). In the strongly polarized regime the potential contribution to the spin splitting δΣ_νω is much larger than the kinetic one χ^-1_ω. Therefore, the latter can be neglected and the result for the exchange interaction (<ref>) reduces to the expression obtained in <cit.> that was successfully applied to the description of many realistic systems <cit.>. Note that in <cit.> the exchange interaction was derived assuming the existence of the collinear spin ground state, while here we show that the limit of applicability of the derived expression is much broader. If the dependence of the three-point vertex on the fermionic frequencies is fully disregarded (second approximation in Eq. <ref>), the exchange interaction reduces to the “Hartree-Fock” approximation J_ = χ^0 -1_ω=0 X^ 0_,ω=0 χ^0 -1_ω=0 <cit.> derived in <cit.>.

Conclusion — 
To conclude, here we derived the action for effective s-d and Heisenberg-like problems for the extended Hubbard model. We observed that by virtue of a local Ward identity the vertex functions of the impurity model can be well approximated, provided its weak dependence on the fermionic frequencies. Our results show that this criterion is indeed satisfied in the AFM phase of the Hubbard model in infinite dimensions when the staggered magnetization is sufficiently large. As a consequence, it is possible to obtain the magnetic susceptibility without a costly measurement of the impurity vertex functions, which is very useful for the realistic multiband calculations. 
For the considered parameters this approximation becomes accurate enough to reach an agreement with the global spin conservation. In finite dimensions this is of importance for a sound description of magnon spectra in accord with Goldstone's theorem.
In the classical limit, the derived spin action reduces to an effective Heisenberg Hamiltonian. In the spin-polarized case the result for the kinetic part of the effective exchange interaction simplifies to the expression derived in <cit.>, which is argued to be a good approximation for the case of many real materials. 
We believe that this approximation can be applied in different and, in particular, more realistic contexts. We further speculate that similar approximations could prove valuable in any physical regime where it can be argued that the behavior of the vertex functions is strongly dominated by the transferred momentum.


The authors thank Hugo Strand and Vladimir Antropov for inspiring discussions. We kindly thank Igor Krivenko for providing the solver for stochastic optimization method for analytical continuation <cit.> used to obtain the data shown in Fig. <ref>. We also would like to thank Junya Otsuki for providing the CTQMC solver <cit.> that was used to obtain the data shown in Fig. <ref>. The work of E.A.S. was supported by the Russian Science Foundation, Grant 18-12-00185. The work of M.I.K. was supported by NWO via Spinoza Prize and by ERC Advanced Grant 338957 FEMTO/NANO. M.H., F.K., and A.I.L. acknowledge support from the excellence  cluster “The Hamburg Centre for Ultrafast Imaging - Structure, Dynamics and Control of Matter at the Atomic Scale” and North-German Supercomputing Alliance (HLRN) under the Project No. hhp00040. Also, the work was partially supported by the Stichting voor Fundamenteel Onderzoek der Materie (FOM), which is financially supported by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO).


apsrev4-1






Supplemental Material for 

“Effective Heisenberg model and exchange interaction for strongly correlated systems”






§ EFFECTIVE SPIN PROBLEM FOR THE EXTENDED HUBBARD MODEL


Here we explicitly derive a spin problem for the extended Hubbard model and obtain magnetic susceptibility. For this reason, let us consider the following action written in momentum space

    S =    -∑_,ν,σ c^*_νσ[iν+μ-ε^*_]
     c^†_νσ
    + U∑_,ωn^*_ω↑n^*_-,-ω↓
    +1/2∑_,ω,ςρ^* ς_ω[V^*_]_ςςρ^ ς_ω.

Here, c^*_νσ (c_νσ) are Grassmann variables corresponding to creation (annihilation) of an electron with momentum , fermionic Matsubara frequency ν and spin σ. Quantities ε_ and [V_]_ςς are the Fourier transforms of the hopping amplitude and nonlocal part of an interaction written in the matrix form, respectively. The label ς={c,  s} depicts the charge c and spin s={x,y,z} degrees of freedom, so that U and [V_]_cc=V_ describe the local and nonlocal parts of the Coulomb interaction respectively, and [V_]_ss=-J^ d_/2 is the nonlocal direct ferromagnetic exchange interaction. The latter ensures the following form of the Heisenberg Hamiltonian H=-∑_J^ d_ S^d_S^d_-. Here, we also introduce bosonic variables  ρ^ ς_ω=n^ ς_ω-n^ ς_ω, where n^ς_ω = ∑_νσσ'c^*_νσσ^ς_σσ' c^*_+,ν+ω,σ' is the charge (ς=c) and spin (ς=s) density of electrons with the momentum , bosonic frequency ω and Pauli matrices σ^ς={1,σ^s}.  

Following the standard procedure of the Dual Boson theory <cit.>, the lattice action is divided into the local impurity S_ imp and nonlocal S_ rem parts as

    S_ imp   = 
    -∑_ν,σc^*_νσ[iν+μ-Δ^*_ν]c^*_νσ 
    +U∑_ωn^*_ω↑ n^*_-ω↓ 
    +1/2∑_ω,ςρ^* ς_ω[Λ^*_ω]_ςςρ^ ς_ω, 
     S_ rem   = 
    -∑_k,ν,σc^*_kνσ[Δ^*_ν-ε^*_]
    c^*_kνσ +
    1/2∑_,ω,ςρ^* ς_ω[V^*_-Λ^*_ω]_ςςρ^ ς_ω + ∑_,ω,ς
    j^ * ς_ω ρ^ ς_ω,

where we introduced fermionic Δ_ν and bosonic [Λ_ω]_ςς hybridization functions and sources j^ ς_ω for bosonic variables.
Since here we consider a spin-polarized case of local impurity model, the fermionic hybridization function Δ_νσ becomes spin-dependent. The partition function of our problem is given by the following relation

    Z=∫ D[c^*,c]   e^- S,

where S is the lattice action introduced in Eq. <ref>.
Using a matrix form of the Hubbard–Stratonovich transformation of the remainder term S_ rem (<ref>) one can introduce dual fermionic f^*,f and bosonic variables ϕ^ ς 

    exp{∑_,ν,σ c^*_νσ[Δ^*_νσ-ε^*_]c^*_νσ}   = D_f∫ D[f^*,f]exp{-∑_,ν,σ( f^*_νσ[Δ^*_νσ-ε^*_]^-1f^*_νσ + c^*_νσf^*_νσ + f^*_νσc^*_νσ)},
    exp{∑_,ω,ς(')1/2 ρ^* ς_ω[Λ^*_ω-V^*_]_ςς'ρ^ ς'_ω}   = D_ϕ∫ D[ϕ]exp{-∑_,ω,ς(')( 1/2 ϕ^* ς_ω[Λ^*_ω-V^*_]_ςς'^-1ϕ^ ς'_ω + 
    
    ϕ^* ς_ω ρ^ ς_ω)},

where terms D_f =  det(Δ_νσ-ε_) and D^-1_ϕ = √( det[Λ_ω-V_]) can be neglected when calculating expectation values. 
Rescaling fermionic fields on the Green's function g_νσ of impurity problem (<ref>) as f^*_νσ→f^*_νσg^-1_νσ and f_νσ→g^-1_νσf_νσ, and bosonic field on the susceptibility [χ_ω]_ςς' as ϕ^* ς_ω→ϕ^* ς'_ω [χ^*_ω]^-1_ς'ς,
and shifting bosonic variables, the nonlocal part (<ref>) of the lattice action (<ref>) transforms to

    S_ DB
    =    -∑_,ν,σ f^*_νσg^-1_νσ [ε^*_-Δ^*_νσ]^-1 g^-1_νσf^*_νσ 
    +∑_,ν,σ[ c^*_νσ g^-1_νσf^*_νσ + 
    f^*_νσ g^-1_νσc^*_νσ] +∑_,ω,ς(')ϕ^* ς_ω[χ^*_ω]^-1_ςς'ρ^ ς'_ω
    
       -1/2∑_,ω,ς(')( ϕ^* ς_ω - j^* ς'_ω[χ^*_ω]_ς'ς)
    [χ^*_ω]^-1_ςς”[V^*_ - Λ^*_ω]^-1_ς”ς”'[χ^*_ω]^-1_ς”'ς””( ϕ^ ς””_ω - [χ^*_ω]_ς””ς””'
    j^ ς””'_ω) .

Now, the initial degrees of freedom can be integrated out with respect to the impurity action (<ref>) in the following way

    ∫ D[c^*,c] exp{-∑_i S^ i_ imp - ∑_,ν,σ[ c^*_νσ g^-1_νσf^*_νσ + 
    f^*_νσ g^-1_νσc^*_νσ] 
    -∑_,ω,ς(')ϕ^* ς_ω[χ^*_ω]^-1_ςς'ρ^ ς'_ω
    } = 
    
        Z_ imp×exp{ -∑_,ν,σf^*_νσ g^-1_νσf^*_νσ 
    -1/2∑_,ω,ς(')ϕ^* ς_ω[χ^*_ω]^-1_ςς'ϕ^ ς'_ω - W̃[f,ϕ]},

where Z_ imp is a partition function of the impurity problem.
Here, the interaction part of the action W̃[f,ϕ] is presented as an infinite series of full vertex functions of impurity problem (<ref>) as discussed in <cit.>. The lowest order interaction terms are following 

    W̃[f,ϕ]
       ≃∑_,',∑_ν,ν',ω∑_σ('),ς(')(ϕ^* ς_ωγ^ ς_νω f^*_νσ f^*_+,ν+ω,σ'
    - 1/4 γ^ σσ'σ”σ”'_νν'ω  f^*_νσf^*_+,ν+ω,σ'f^*_'+,ν'+ω,σ” f^*_'ν'σ”'),

where the full three-point vertex function (and its Hermitian conjugate) is defined as

    γ^ ς_νω   = ∑_ς'[χ_ω]^-1_ςς'ρ^ ς'_ωc^*_νσ c^*_ν+ω,σ'_ impg^-1_νσ g^-1_ν+ω,σ' 
    = ∑_ς'c^*_νσ c^*_ν+ω,σ' ρ^ ς'_ω_ imp[χ_ω]^-1_ςς'g^-1_νσ g^-1_ν+ω,σ', 
    
    [γ^ ς_νω]^*   = γ^ ς^*_ν+ω,-ω = ∑_ς'c^*_ν+ω,σ' c^*_νσ ρ^* ς'_ω_ imp[χ_ω]^-1_ς'ς g^-1_ν+ω,σ' g^-1_νσ.

The full four-point vertex determined in the particle-hole channel is equal to

    γ^ σσ'σ”σ”'_νν'ω =
    c^*_νσ c^*_ν+ω,σ' c^*_ν'+ω,σ” c^*_ν'σ”'_ c  imp g^-1_νσ g^-1_ν+ω,σ' g^-1_ν'+ω,σ” g^-1_ν'σ”'.

Therefore, the initial lattice problem (<ref>) transforms to the following dual action 

    S̃
    =    -∑_,ν,σ f^*_νσg^-1_νσ [ε^*_-Δ^*_νσ]^-1 g^-1_νσf^*_νσ 
    +∑_,ν,σf^*_νσ g^-1_νσf^*_νσ 
    +1/2∑_,ω,ς(')ϕ^* ς_ω[χ^*_ω]^-1_ςς'ϕ^ ς'_ω + W̃[f,ϕ] 
       -1/2∑_,ω,ς(')( ϕ^* ς_ω - j^* ς'_ω[χ^*_ω]_ς'ς)
    [χ^*_ω]^-1_ςς”[V^*_ - Λ^*_ω]^-1_ς”ς”'[χ^*_ω]^-1_ς”'ς””( ϕ^ ς””_ω - [χ^*_ω]_ς””ς””'
    j^ ς””'_ω).


In order to come back to the original bosonic variables, one can perform the third Hubbard-Stratonovich transformation as 

    exp{1/2∑_,ω,ς(')( ϕ^* ς_ω - j^* ς'_ω[χ^*_ω]_ς'ς)
    [χ^*_ω]^-1_ςς”[V^*_ - Λ^*_ω]^-1_ς”ς”'[χ^*_ω]^-1_ς”'ς””( ϕ^ ς””_ω - [χ^*_ω]_ς””ς””'
    j^ ς””'_ω)} = 
    
    
    
       D_ρ̅∫ D[ρ̅]exp{-∑_,ω,ς(')(1/2 ρ̅^* ς_ω [V^*_ - Λ^*_ω]_ςς'ρ̅^ ς'_ω
    - ϕ^* ς_ω[χ^*_ω]^-1_ςς'ρ̅^ ς'_ω
    + j^ * ς_ωρ̅^ ς_ω)}.

Comparing this expression to the Eq. <ref>, one can see that sources j^* ς_ω introduced for the initial degrees of freedom ρ^ ς_ω are also the sources for new bosonic fields ρ̅^ ς_ω. Therefore, fields ρ̅^ ς_ω indeed represent initial degrees of freedom and have the same physical meaning as original composite bosonic variables ρ^ ς_ω=∑_νσσ'c^*_νσσ^ς_σσ' c^*_+,ν+ω,σ'-n^ ς_ω of the lattice problem (<ref>). Nevertheless, ρ̅^ ς_ω can now be treated as elementary bosonic fields that have a well-defined propagator, since they are introduced as a decoupling fields of dual degrees of freedom ϕ^ ς_ω and therefore, independent on fermionic variables c^*_νσ (c_νσ). Taking sources to zero and replacing ρ̅^ ς_ω by ρ^ ς_ω, dual bosonic fields can be integrated out with respect to the Gaussian bosonic part of the dual action as

    ∫ D[ϕ^*,ϕ] exp{-∑_,ω,ς(')(1/2 ϕ^* ς_ω[χ^*_ω]^-1_ςς'ϕ^ ς'_ω
    - ϕ^* ς_ω[χ^*_ω]^-1_ςς'ρ̅^ ς'_ω) - W̃[f,ϕ]}
    = 
     Z_ϕ×exp{1/2∑_,ω,ς(')ρ^* ς_ω[χ^*_ω]^-1_ςς'ρ^ ς'_ω - W[f,ρ] },

where Z_ϕ is a partition function of the Gaussian part of the bosonic action.
Here we restrict ourselves to the lowest order interaction terms of W̃[f,ϕ] shown in Eq. <ref>. Then, the integration of dual bosonic fields in Eq. <ref> simplifies and W[f,ρ] keeps an efficient dual form of W̃[f,ϕ] (<ref>) with replacement of bosonic variables ϕ^ ς→ρ̅^ ς. Also the four-point vertex becomes irreducible with respect to the full local bosonic propagator χ_ω, as can be seen from the works of <cit.>, while the three-point vertex γ_νω remains invariant

    W[f,ρ]
    =∑_,',∑_ν,ν',ω∑_ς(')(ρ^* ς_ωγ^ ς_νω f^*_νσ f^*_+,ν+ω,σ'
    - [γ - θ]^ ςς'_νν'ω f^*_νσf^*_+,ν+ω,σ'f^*_'+,ν'+ω,σ” f^*_'ν'σ”').

Here,

    θ^ ςς'_νν'ω = - γ^ ς_νω [χ_ω]_ςς'[γ^ ς'_ν'ω]^*

is the full reducible bosonic contribution to the full local four-point vertex γ^ ςς'_νν'ω introduced in <cit.> and spin labels σ, σ', σ”, σ”' are fixed by the channel indices ς, ς'.
Therefore, the problem transforms to the following action of an effective s-d model 

    S_s-d =    -∑_,ν,σ f^*_ kνσG̃^-1_0f^*_ kνσ 
    -1/2∑_,ω,ς(')ρ^* ς_ω[X^*_ E]^-1_ςς'ρ^ ς'_ω + W[f,ρ],

where [X_ E]_ςς' = [χ^-1_ω + Λ^1_ω - V^1_]^-1_ςς' is the susceptibility of the extended dynamical mean-field theory.
As it is shown below, when the three-point vertex function γ'_νω of impurity problem that connects two fermionic propagators and interaction is close to unity (<ref>), the main contribution to the local four-point vertex is given by the full reducible bosonic contribution, i.e. γ≃θ, or diagrammatically

    < g r a p h i c s >
     .

Here, the dotted wave line depicts full local bosonic propagator and the minus sign in Eq. <ref> appears due to Feinman rules <cit.>. Then, the interaction part of the action (<ref>) takes the most simple form that contains only three-point vertex functions

    W'[f,ρ]
    ≃∑_,∑_ν,ω∑_σ('),ς(')
    ρ^* ς_ωγ^ ς_νω f^*_νσ f^*_+,ν+ω,σ'
    .




 §.§ Transformation of spin basis

Let us consider an effective impurity model in the spin-polarized case. For easier description, one can transform spin variables from the s={x,y,z} to the s={+,-,z} basis as S^±=(ρ^x±iρ^y)/2. In the spin-polarized case fluctuations in the charge and spin z channels are yet entangled, but the ± spin channel can be separated in the collinear case. Thus, for a correct account for spin fluctuations, one may consider correlations only in the ± spin channel and the contribution of the z channel to the exchange interaction can be later restored from the symmetry arguments. In is worth mentioning that the transformation {x,y}→{+,-} is very useful for calculation of physical observables, since it diagonalizes the spin susceptibility. Nevertheless, one has to remember that operators S^+ and S^- are not Hermitian. Therefore, components of bosonic operator in matrix representation in the old and new basis are defined as

    ρ̂_ω = 
    [ ρ^ x_ω; ρ^ y_ω;        ];   Ŝ_ω = 
    [ S^+_ω; S^-_ω;       ];   ρ̂^*_ω = 
    (ρ^ x_-,-ω, ρ^ y_-,-ω);   Ŝ^*_ω = 
    (S^-_-,-ω,  S^+_-,-ω).

Connection between these bases can be obtained using the following matrix transformation

    Ŝ^*_ω = ρ̂^*_ω×Â,
        or    (S^-_-,-ω,  S^+_-,-ω) = 
    (ρ^ x_-,-ω, ρ^ y_-,-ω) ×[  1/2  1/2; -i/2  i/2;      ]

and

    Ŝ_ω = B̂×ρ̂_ω,
        or    [ S^+_ω; S^-_ω;       ] = 
    [  1/2  i/2;  1/2 -i/2;      ]×[ ρ^ x_ω; ρ^ y_ω;        ].

Then, all matrices M̂_xy involved in above derivations can also be transformed to the new basis M̂_± as

    M̂_± = Â^-1×M̂_xy×B̂^-1.

In particular, the matrix form of the nonlocal interaction [V_]_ςς' remains diagonal

    [V_]_± = 
    Â^-1×[ -J^ d_/2        0;        0 -J^ d_/2;          ]×B̂^-1 = 
    [ -J^ d_      0;      0 -J^ d_;        ]

and inverse susceptibility is transformed to a diagonal form as

    [χ_ω]_±^-1 = 
    Â^-1×[ χ^xx_ω χ^xy_ω; χ^yx_ω χ^yy_ω;        ]
    ^-1×B̂^-1 = 
    1/χ^xx_ωχ^yy_ω - χ^xy_ωχ^yx_ω[  χ^xx_ω+χ^yy_ω + iχ^xy_ω - iχ^yx_ω  χ^xx_ω-χ^yy_ω + iχ^xy_ω + iχ^yx_ω; -χ^xx_ω+χ^yy_ω + iχ^xy_ω + iχ^yx_ω  χ^xx_ω+χ^yy_ω - iχ^xy_ω + iχ^yx_ω;                                    ].

Defining χ^+-_ω = -S^+_ω S^-_-ω = 1/4(χ^xx_ω+χ^yy_ω - iχ^xy_ω + iχ^yx_ω) and χ^-+_ω = -S^-_ω S^+_-ω = 1/4(χ^xx_ω+χ^yy_ω + iχ^xy_ω - iχ^yx_ω), and taking into account that χ^xx_ω = χ^yy_ω and χ^xy_ω = -χ^yx_ω, one gets that
χ^xx_ωχ^yy_ω - χ^xy_ωχ^yx_ω = 4 χ^+-_ωχ^-+_ω
and 

    [χ_ω]_±^-1  = 
    
    
    
    
    
    [ [χ^+-_ω]^-1           0;           0 [χ^-+_ω]^-1;             ].




 §.§ Magnetic susceptibility


In order to obtain the effective problem written in terms of bosonic degrees of freedom only, one can integrate out dual fermionic degrees of freedom from the Eq. <ref>.
Taking into account transformation of the spin basis presented above, the spin ± part of the effective action reads

    S_ spin   = -1/2∑_,ω S^ -_ω[X_ω^ -+]^-1  S^ +_-,-ω
    -1/2∑_,ω S^ +_ω[X_ω^ +-]^-1  S^ -_-,-ω

The first approximation for the spin susceptibility X^-+_ω can be obtained after expanding the simplified form of interaction W[f,ρ] given by Eq. <ref> up to the second order with respect to bosonic fields ρ in the expression for the partition function of the action (<ref>). This results in 

    [X^(2)_ω]^-1   = J^ d_ q + Λ^1_ω + χ^-1_ω - Π̃^(2)_ω,

where

    Π̃^(2)_ω=∑_,νγ^ -_ν+ω,-ω G̃^2_+,ν+ω↑G̃^2_ν↓ γ^ +_ν,ω

is the second order polarization function and χ^+_ω=χ^-+_ω and Λ^+_ω=Λ^-+_ω are the spin susceptibility and bosonic hybridization function of impurity problem, respectively. Hereinafter, ± spin labels are omitted for simplicity wherever they are not crucial for understanding. The three-point vertex functions in the spin channel are defined as in Eq. <ref>, or explicitly as

    γ^ +_ν,ω   =   c^*_ν↓ c^*_ν+ω↑ S^-_ω_ imp[χ^-+_ω]^-1g^-1_ν↓ g^-1_ν+ω↑, 
    γ^ -_ν+ω,-ω   =  
    S^+_-ω c^*_ν+ω↑ c^*_ν↓_ imp[χ^-+_ω]^-1 g^-1_ν+ω↑ g^-1_ν↓.

The more accurate approximation for the spin susceptibility can be found when expanding the full form of interaction  W[f,ρ] given by Eq. <ref> up to the second order with respect to bosonic fields ρ as previously. Using the ladder approximation, one gets 

    [X^ ladd_ω]^-1 =
    J^ d_ q + Λ^1_ω + χ^-1_ω - Π^ ladd_ω,

where the polarization function Π^ ladd_ω expressed in the matrix form in the space of fermionic frequencies ν,ν' reads

    Π^ ladd_ω = {γ̂^-_ω X̂̃̂^0_ω[ I + (γ̂^+_ω-θ̂^+_ω) X̂̃̂^0_ω]^-1γ̂^+_ω}.
 
Here, I is the identity matrix in the same space. Multiplication and inversion should be understood as a standard matrix operations. For simplicity, we omit the fermionic indices  wherever they are not crucial for understanding. 
The trace is taken over the external fermionic indices.  
Matrix elements of the bare dual spin susceptibility X̃^0_ω and three-point vertex function γ_ω are defined as X̃^0_ω; νν' = ∑_G̃_+,ν+ω↑G̃_ν↓ δ_νν' and γ^±_ω; νν' = γ^±_νω δ_νν', where γ^±_νω are defined in Eq. <ref>. The four-point vertex functions γ_νν'ω and θ_νν'ω in the ± spin channel are defined above in Eqs. <ref> and (<ref>), or explicitly as

    γ^1_νν'ω   = γ^ ↓↑↑↓_νν'ω =
    c^*_ν↓ c^*_ν+ω↑ c^*_ν'+ω↑ c^*_ν'↓_ c  imp g^-1_ν↓ g^-1_ν+ω↑ g^-1_ν'+ω↑ g^-1_ν'↓, 
    θ^1_νν'ω   = - γ^ +_νω χ_ωγ^ -_ν'+ω,-ω.

Substituting the above expressions to the Eq. <ref>, one recovers conserving result for the spin susceptibility provided by the ladder DB approach <cit.> in the case of the constant bosonic hybridization function Λ <cit.>

    [X^ ladd_ω]^-1 = J^ d_ + Λ + [X^ DMFT_ω]^-1.

Here, 

    X^ DMFT_ω=χ^2_ω+χ^2_ωΠ̃^ ladd_ωχ^2_ω

and Π̃^ ladd_ω is the dual polarization function in the ladder form <cit.> given by the following matrix form in the space of fermionic frequencies ν,ν'

    Π̃^ ladd_ω = {γ̂^-_ω X̂̃̂^0_ω[ I + γ̂^+_ωX̂̃̂^0_ω]^-1γ̂^+_ω}.
 
As it was already noted in <cit.>, the difference between the lattice (<ref>) and dual (<ref>) polarization functions is that the first one is irreducible with respect to the (local and nonlocal parts of) EDMFT susceptibility X_ E, while the dual one is irreducible only with respect to the bare dual susceptibility, which is identically equal to the nonlocal part of X_ E.

Expression for the spin susceptibility (<ref>) can be rewritten in the more convenient way. 
For this reason one can define the two-particle irreducible (2PI) vertex function in the ± spin channel as 

    γ̂^ 2PI_ω = γ̂^1_ω[I - χ̂^0_ωγ̂^1_ω]^-1,

where the matrix elements of the bare local spin susceptibility are χ^0_ω; νν' = g_ν+ω↑g_ν↓ δ_νν'. 
Then, the spin susceptibility of the impurity problem can be expressed as

    χ_ω = -S^-_ω S^+_-ω = {χ̂^0_ω - χ̂^0_ω γ̂_ω χ̂^0_ω} = {χ̂^0_ω[I+γ̂^ 2PI_ωχ̂^0_ω]^-1}.
 
Rewriting the relation for the dual polarization function Π̃^ ladd_ω (<ref>) through the 2PI vertex function and using the exact relation between the three- and four-point vertex functions of impurity problem

    γ^+_νω   =  
    c^*_ν↓ c^*_ν+ω↑ S^-_ω_ impχ^-1_ωg^-1_ν↓ g^-1_ν+ω↑ 
    = ∑_ν'c^*_ν↓ c^*_ν+ω↑ 
    c^*_ν'↓ c^*_ν'+ω,↑_ impχ^-1_ωg^-1_ν↓ g^-1_ν+ω↑
       = ∑_ν'{δ_ν,ν' -  γ^1_ν,ν',ω g_ν'+ω↑ g_ν'↓}χ^-1_ω
    = ∑_ν'[δ_νν' + γ^ 2PI_ν,ν',ω  g_ν'+ω↑g_ν'↓]^-1χ^-1_ω,

and the fact that in the case of zero dual self energy Σ̃_ν=0 the following relation holds 

    X̃^0_ω; νν' + χ^0_ω; νν' = X^0_ω; νν' = ∑_G_+,ν+ω↑G_ν↓ δ_νν',

one finds that

    X^ DMFT_ω = χ^1_ω + χ^1_ωΠ̃^ ladd_ωχ^1_ω = 
    {X̂^0_ω[I + γ̂^  2PI_ωX̂^0_ω]^-1}

is the DMFT-like <cit.>
susceptibility written  in terms of the 2PI vertex functions of impurity model and lattice Green's functions. 
Therefore, the spin susceptibility (<ref>) derived within the ladder Dual Boson approach <cit.> can be rewritten as

    X^ ladd_ω = {X̂^0_ω[ I + (γ̂^ 2PI_ω + I[J^ D_ q + Λ]) X̂^0_ω]^-1}.




 §.§ Classical Heisenberg Hamiltonian
 
In order to map the initial problem onto a classical Heisenberg Hamiltonian the
spin variables S^±_ω in Eq. (<ref>) have to be replaced by the classical vectors S_ of the length S_z.

Then, an effective exchange interaction J_ can be defined as a nonlocal part of the inverse susceptibility at zero bosonic frequency <cit.>.
After all, the action (<ref>) maps on an effective Heisenberg Hamiltonian 

    H_ spin = - ∑_ J_  S_  S_-,

where the contribution from the z spin channel is restored from the requirement of rotational invariance. Here, the effective exchange interaction obtained from the simplified form of magnetic susceptibility (<ref>) is

    J_ = J^ d_ - ∑_,νγ^ -_ν,ω=0 G̃^2_+,ν↑G̃^2_ν↓ γ^ +_ν,ω=0.

and the exchange interaction in the ladder approximation obtained from the Eq. <ref> reads

    J^d_ = J^ d_ - Π̃^ ladd_,ω=0[1+χ^d_ω=0 Π̃^ ladd_,ω=0]^-1.




 §.§ Ward identity for the vertex function of impurity model


When the system exhibits mostly bosonic fluctuation, one can expect that local vertex functions of impurity problem are mostly described by the bosonic frequency ω, while the dependence on fermionic frequencies ν,ν' can be averaged. In order to account for single electronic degrees of freedom correctly, the averaging procedure over the fermionic frequencies is carried out using Ward identity for the two-particle irreducible four-point vertex function of the impurity problem <cit.> as

    Σ_ν+ω↑-Σ_ν↓ = -∑_ν”γ^ 2PI_ν,ν”,ω(g_ν”+ω↑-g_ν”↓) ≃ - γ^ 2PI_ν,ν”,ω{∑_ν”(g_ν”+ω↑-g_ν”↓)} 
    = - 2 γ^ 2PI_ν,ν”,ω S^z.

Then, one can approximate the two-particle irreducible vertex function as

    γ^ 2PI_ν,ν”,ω≃γ^ 2PI_ν,ν”,ω = - Σ_ν+ω↑-Σ_ν↓/2S^z = 
    -δΣ_νω.

The three-point vertex function (<ref>) is then simplified as

    γ^+_νω   = ∑_ν'ν”{δ_ν,ν' -  γ^ 2PI_ν,ν”,ω[δ_ν”ν' + g_ν”+ω↑g_ν”↓γ^ 2PI_ν”,ν',ω]^-1
    g_ν'+ω↑ g_ν'↓}χ^-1_ω
       ≃{1 - γ^ 2PI_ν,ν”,ω∑_ν'ν”{[δ_ν”ν' + g_ν”+ω↑g_ν”↓γ^ 2PI_ν”,ν',ω]^-1
    g_ν'+ω↑ g_ν'↓}}χ^-1_ω
       =  {1 - γ^ 2PI_ν,ν”,ω χ^1_ω}χ^-1_ω
    = χ^-1_ω +δΣ_νω≃χ^ 0 -1_ω,

where χ^ 0_ω = ∑_νg_ν+ω↑ g_ν↓.
Similarly, one can show that 

    γ^-_ν+ω,-ω≃χ^-1_ω + δΣ_νω≃χ^ 0 -1_ω

and the magnetic susceptibility (<ref>) can be written as

    [X^(2)_ω]^-1   = J^ d_ q + Λ^1_ω + χ^-1_ω - ∑_,ν(χ^-1_ω + δΣ_νω) G̃_+,ν+ω,↑ G̃_,ν,↓ (χ^-1_ω + δΣ_νω)
       = J^ d_ q + Λ^1_ω + χ^-1_ω - χ^ 0 -1_ω X̃^0_ω χ^ 0 -1_ω.

Here, we also introduce X̃^0_ω = ∑_,νG̃_+,ν+ω↑ G̃_,ν↓.

The ladder form of the magnetic susceptibility (<ref>) can also be simplified. Taking into account that the last approximation in Eq. <ref> is nothing else than averaging of the 2PI four-point vertex function (<ref>) over the second fermionic frequency 

    γ^ 2PI_νν'ω≃γ^ 2PI_νν'ω=χ^-1_ω-χ^0 -1_ω,

one gets the following RPA-like approximation for the magnetic susceptibility in the ladder approximation (<ref>)

    X^ ladd_ω = X^0_ω[1+(γ^ 2PI_νν'ω + Λ + J^ d_ q)X^0_ω]^-1,

where the bare lattice magnetic susceptibility X^0_ω = ∑_νG_+,ν+ω↑G_ν↓ was introduced.

As it is shown below, in the strongly polarized regime the three-point vertex function of impurity problem γ'_νω that connects two fermionic propagators and interaction is close to unity (<ref>). Then, the local polarization function of impurity can be approximated as Π_ω≃χ^0_ω and the full local susceptibility in the spin channel reads

    χ^-1_ω = χ^0 -1_ω -  U^±_ω.

Here, U^±_ω=-U+Λ is the bare interaction of impurity problem in the spin channel.
Then, the averaged 2PI four-point vertex function (<ref>) reads γ^ 2PI_νν'ω≃U-Λ and one finally gets the following simple expression for the magnetic susceptibility (<ref>)

    X^ ladd_ω = X^0_ω[1+(U + J^ d_ q)X^0_ω]^-1.




§ APPLICATION: THE HUBBARD MODEL ON THE HYPERCUBIC LATTICE IN INFINITE DIMENSIONS


We consider the half-filled Hubbard model 

    H=-(2d)^-1/2∑_⟨ ij⟩σ c_iσ^† c_jσ+U∑_i n_i↑n_i↓

on the hypercubic lattice in infinite dimensions d→∞, where the summation over ij runs over nearest neighbors. In this limit the non-interacting density of states becomes a Gaussian, D(ϵ)=(2π)^-1/2e^-ϵ^2/2 <cit.>. At low temperatures this system favors antiferromagnetic order over paramagnetism. Within the symmetry-broken phase one has to consider two sublattices A and B of the bipartite hypercubic lattice with a staggered magnetization, ⟨ m^A⟩=-⟨ m^B⟩=m. In a bipartite ordered state the volume of the Brillouin zone (BZ) is halved, such that Fourier transforms may only be performed up
to the magnetic unit cell, see, for example, <cit.>. 
In the reduced Brillouin zone (RBZ) the noninteracting Hamiltonian reads,

    H_0=∑_σ(a^*_σ,b^*_σ)
        [  0 ε_; ε_  0 ][ a_σ; b_σ ],

where a^(*)_σ and b^(*)_σ annihilate (create) a σ-electron with momentum  in sublattice A and B, respectively.
ε_ is the dispersion of the hypercubic lattice and  a vector of the RBZ. Therefore, the Green's function becomes a two-by-two matrix Ĝ=(G^AA, G^AB; G^BA, G^BB) in sublattice space. Since the Hubbard model in infinite dimensions is an exact limit of DMFT, the lattice model (<ref>) is mapped exactly to a single-site Anderson impurity model (AIM). Therefore, the self-energy Σ_νσ of Ĝ_νσ is local and it reads

    Ĝ_νσ=
    [  ζ_ν  -ε_;  -ε_ ζ_νσ ]^-1,

where ζ_νσ=iν+μ-Σ_νσ. The impurity ĝ_νσ and local part of the lattice Green's function Ĝ_loc=∑_Ĝ_ are tied via the following prescription

    Ĝ_loc,νσ=ĝ_νσ=∫_-∞^∞D(ϵ)dϵ/ζ_νσζ_ν-ϵ^2[  ζ_ν    0;    0 ζ_νσ ].

The momentum summation was rewritten as an integral over the density of states D(ϵ) of the hypercubic lattice. It was used that the off-diagonal elements of Green's function are an odd function of ϵ and thus vanish upon integration. By symmetry, an exchange of the sublattice indices A↔ B is equivalent to a flip of the spin label σ↔. The prescription is satisfied by fixing the dynamical Weiss field 𝒢^-1_νσ=G^-1_loc,νσ+Σ_νσ of the AIM self-consistently.




 §.§ DMFT susceptibility of the ordered phase


In order to calculate the transversal spin susceptibility of the Hubbard model in the antiferromagnetically ordered phase we introduce the bare susceptibility.
On a bipartite lattice it is in general necessary to consider two-particle quantities with four indices a,b,c,d. The bubble is then a 4×4 matrix given by the tensor product Ĝ_⊗Ĝ^T_. The locality of the irreducible vertex in DMFT allows to consider the Bethe-Salpeter equation only in a 2×2 subspace, where the bare susceptibility is given by the following point-wise product Ĝ_∘Ĝ^T_

    X̂^0_ω;ν=   ∑_Ĝ_+,ν+ω∘Ĝ^T_ν
    
    =   ∑_1/ζ_ν+ω,ζ_ν+ω,-ε^2_+1/ζ_νζ_ν-ε^2_[ ζ_ν+ω,ζ_ν,      ε_+ε_;      ε_+ε_ ζ_ν+ω,ζ_ν, ]
    
    =   ∬_-∞^∞ dϵ_1dϵ_2
    D_(ϵ_1,ϵ_2)/(ζ_ν+ω,ζ_ν+ω,-ϵ^2_1)(ζ_νζ_ν-ϵ^2_2)[ ζ_ν+ω,ζ_ν,     ϵ_1ϵ_2;     ϵ_1ϵ_2 ζ_ν+ω,ζ_ν, ],

where  and  are vectors of the reduced Brillouin zone (RBZ). Here, the momentum summation leads to a double integral involving the expression D_(ϵ_1,ϵ_2). This reduces to D(ϵ_1)D(ϵ_1) for any generic wave vector . The term “generic” may be understood such that  is a vector of the RBZ with an infinite number of random entries (see <cit.> and references therein).
As a consequence, the integrals in Eq. (<ref>) factorize and the bubble is given as ĝ_ν+ω∘ĝ_ν.
Hence, the nonlocal bubble

    X̃^0,ab_ω;ν=∑_ G^ab_+,ν+ω,G^ba_ν-g^a_ν+ω,g^a_νδ_ab,

vanishes identically at generic . In the following, we consider the non-generic vector _0=0 of the RBZ, where D__0(ϵ_1,ϵ_2)=δ(ϵ_1-ϵ_2)D(ϵ_1). This may be used to eliminate one of the integrals in Eq. (<ref>), the remaining integral is solved numerically. Vectors  of the Brillouin zone (BZ) are marked with a tilde.

From the nonlocal bubble in Eq. (<ref>) and from the full local four-point vertex of impurity model
γ^ab_νν'ω=γ^ a_νν'ω δ_ab one obtains the
T-matrix F^ab_ω; νν' via the Bethe-Salpeter equation (BSE),

    F̂^-1_ω = γ̂^-1_ω +X̂̃̂^0_ω,

where X̃^0,ab_ω; νν'=X̃^0,ab_ω;ν δ_νν' and Ô^-1 denotes a super-matrix inversion
with respect to the indices (a,ν) and (b,ν').
One further calculates the dual polarization in the ladder approximation (<ref>) as

    Π̂̃̂^ ladd_ω=Tr_νν'[γ̂_ωX̂̃̂^0_ω(I - V.C.)γ̂_ω],

where Tr_νν' denotes a trace over fermionic frequencies and 
V.C. indicates vertex corrections given by F̂_ωX̂̃̂^0_ω. Finally, the lattice susceptibility is obtained using the relation (<ref>). 
Further, we consider an approximation for the magnetic susceptibility given by equation (<ref>) in the case of Hubbard model

    [X^(2)_ω]^-1   = χ^-1_ω - Π̃^(2)_ω,

where Π̃^(2)_ω is obtained when neglecting vertex corrections in Eq. (<ref>). The case of magnetic susceptibility (<ref>) where the polarization function Π̃^ ladd_ω is approximated by the second-order correction Π̃^(2)_ω is also considered. However, is does not provide a good approximation for the exact result of Eq. <ref> as shown in the main text.



 §.§ Numerical calculations



The numerical calculations are performed using 10^8 measurements with 50 Monte-Carlo moves between them. Aside from the segment insertion and removal we also use the shift and the double move as well as the spin-flip, a global move. We measure the Green's function in the Legendre basis with 35 coefficients. The Hilbert-transform for the local Green's function is done on an energy mesh of ω∈ [-20, 20] with 4000 mesh points. The initial DMFT-cycle is performed with an external magnetic field, that is switched off for the following cycles. For the DMFT updates we use a mixing parameter of 0.5.

We fit the model S^z = √(T_N - T) to the DMFT results to estimate T_N≈ 0.186, see Fig. <ref>. This refined scan is done for U=5 with a more dense temperature mesh. Only data points in proximity to the transition were taken into account and the domain of the paramagnetic region is treated by a Heaviside step-function.

For simplicity, we calculate an approximated versions (<ref>) and (<ref>) of the magnetic susceptibilities (<ref>) and (<ref>) for the Hubbard model (V_, J^ d_, Λ=0), respectively. We note that the wave vector _0 of the RBZ maps to two vectors _0=_0 and _π=(π,...,π) of the BZ. In the paramagnet this mapping diagonalizes the susceptibility matrix X̂(_0)=(X^AA, X^AB; X^BA, X^BB), where the diagonal elements are X(_π)=X^AA+X^BB-X^AB-X^BA and X(_0)=X^AA+X^BB+X^AB+X^BA. In the ordered phase the same mapping does not diagonalize X̂, since the offdiagonal element X^±(_0)=X^AA-X^BB+X^AB-X^BA does not vanish. Approaching T_N from above, X(_π) diverges and it remains divergent in the ordered phase, signaling that the crystal is prone to a spontaneous tilt of its magnetization axis. We verified in our calculations that at U=5 and T=0.1<T_N one of the two eigenvalues of the BSE (<ref>) is very close to unity, |λ__π|≈0.993, and that this channel corresponds to X(_π). We account the slight deviation of this eigenvalue from unity to our approximation of the impurity vertex γ. The second eigenvalue, however, remains smaller than one, |λ__0|≈0.715, and belongs to the homogenous susceptibility X(_0). In the main text we show the real and imaginary parts of the lattice susceptibility X̂(_0), which corresponds to X̂(_0) =  X(_0) and X̂(_0) =  X^±(_0).



 §.§ Ward identity

We deduce two exact statements about the dynamical homogenous susceptibility X̂(_0,ω) from the Ward identity of the two-particle correlation function G^abcd_kk'q=-a^*_kb^*_k+q,c^*_k'+q,d^*_k', where each of the operators labeled with a,b,c,d denotes either a^(*) or b^(*), respectively, and k=(,ν), q=(,ω) abbreviate momenta from the RBZ and Matsubara frequencies. From the equation of motion ∂_τρ^a_=[H,ρ^a_] of the density operator
ρ^a_=∑_ a^*_a^*_+ one obtains the Ward identity (see, for example, <cit.>),

    -iω∑_k'(G^aaaa_kk'q+G^aabb_kk'q)
    +∑_k'[ε_k'+q-ε_k'](G^aaba_kk'q+G^aaab_kk'q)
    =G^aa_k+q-G^aa_k,  (a≠ b),

where ∑_k implies a summation over the RBZ and Matsubara frequencies.
Evaluating Eq. (<ref>) at q_0^+=(_0=0,ω^+>0) the term in the second line vanishes.
Upon summation over k and using that ∑_kk'G^aabb_kk'q_0^+=X^ab(_0,ω^+) it follows,

    -iω^+[X^aa(_0,ω^+)+X^ab(_0,ω^+)]=m^a,

where it was also used that ∑_k G^aa_kσ=n^a_σ and m^a=⟨ n^a_⟩-⟨ n^a_⟩.
Adding up above relation for a=A, b=B and a=B, b=A it follows that

    X(_0,ω^+)=X^AA(_0,ω^+)+X^AB(_0,ω^+)+X^BB(_0,ω^+)+X^BA(_0,ω^+)=0,

since m^A=-m^B=m. Subtraction likewise leads to

    X^±(_0,ω^+)=X^AA(_0,ω^+)+X^AB(_0,ω^+)-X^BB(_0,ω^+)-X^BA(_0,ω^+)=2im/ω^+.

Eqs. (<ref>) and (<ref>) follow from the equation of motion of the total spin density, ρ^A__0+ρ^B__0,
and are therefore necessary criteria for global spin conservation.



§ SPIN POLARIZED SOLUTION OF ATOMIC PROBLEM

One can perform an exact diagonlization of a magnetically polarized single orbital Hubbard atom at the half-filling. The thermodynamic potential operator of the atom is given by

    Ĥ-μN̂=
    ∑_σΔ^*_σ
    c^*_σc^*_σ
    +Un_↑n_↓,

Here Δ_↑,↓=± B-μ and the magnetic field B is considered much larger than the temperature T≡1/β.    
The system has four eigenstates |0>, |↑⟩, |↓⟩ and |↑↓⟩ with the corresponding energies E_0=0, E_↑,↓=± B-μ and E_↑↓=U-2μ. Half-filling corresponds to μ=U/2, so that E_↑↓=0. Indeed, the partition function for μ=U/2 is

    Z=2+e^β(μ-B)+e^β(μ+B)≈ e^β(B+μ),

and the average filling is given by (1×e^β(μ-B)+e^β(μ+B)+2×1)/Z=1, where we used that Bβ≫1.

The non-zero matrix elements of the creation and annihilation operators are

    <↑|c^*_↑|0>=1;  <↓|c^*_↓|0>=1;  <↑↓|c^*_↑|↓>=1;  <↑↓|c^*_↓|↑>=-1.

    
Now we use the Lehmann representation to obtain the Green's functions of the system

    g_νσ=1/Z∑_ij|<i|c_σ|j>|^2e^-βE_i+e^-βE_j/iν+E_i-E_j.

This yields

    g_ν↑   =1/Z[|<0|c_↑|↑>|^21+e^β(μ-B)/iν+μ-B+|<↓|c_↑|↑↓>|^21+e^β(μ+B)/iν-μ-B],
    
    g_ν↓   =1/Z[|<0|c_↓|↓>|^21+e^β(μ+B)/iν+μ+B+|<↑|c_↓|↑↓>|^21+e^β(μ-B)/iν-μ+B],

or, using Bβ≫1,

    g_ν↑=1/iν-μ-B≃1/iν-U/2,    
    g_ν↓=1/iν+μ+B≃1/iν+U/2.

Now we calculate the magnetic susceptibility

    χ^-+_ω=-1/Z∫_0^βdτ e^iωτS^-(τ) S^+(0),

where S^±(τ) are Heisenberg representations of S^± operators.
The non-zero matrix elements of the latter are

    ↑|S^+|↓=↓|S^-|↑=1.

Lehmann representation reads

    χ^-+_ω=-1/Z∑_ij|i|S^-|j|^2[(e^-βE_j-e^-βE_i)(1-δ_E_i,E_j)/iω+E_i-E_j+βδ_E_i,E_jδ_ω0] =1/iω-2B.


Finally we turn to calculating of χ^↓↑+(τ_1,τ_2)≡T_τc^*_↓(τ_1)
c^*_↑(τ_2)S^+(0).
Unlike the previous cases, here we have to explicitly consider the time-ordering operator.

    χ^σσ'ς(τ_1,τ_2)=c^*_σ(τ_1) c^*_σ'(τ_2) S^ς(0)θ(τ_1-τ_2) - c^*_σ'(τ_2)c^*_σ(τ_1) S^ς(0)θ(τ_2-τ_1).

The usual trick here is to split the integration region 0<τ_1,τ_2<β in the Fourier transform integral into two parts: ∫_0^β∫_0^βdτ_1dτ_2…=∫_0^βdτ_1∫_0^τ_1dτ_2…+∫_0^βdτ_2∫_0^τ_2dτ_1… and swapping the integration variables in the second term. This immediately gives for

    χ^σσ'ς(ν_1,ν_2)≡∫_0^β∫_0^βdτ_1dτ_2
    e^iν_1τ_1+iν_2τ_2χ^σσ'ς(τ_1,τ_2)

the following Lehmann representation

    χ^σσ'ς(ν_1,ν_2)=1/Z[∑_ijki|c^*_σ|jj|c^*_σ'|kk| S^ς|if_ijk(ν_1,ν_2)
    -∑_ijki|c^*_σ'|jj|c^*_σ|kk| S^ς|if_ijk(ν_2,ν_1)],

where

    f_ijk(ν_1,ν_2)
       = ∫_0^β dτ_1 ∫_0^τ_1 dτ_2 e^-βE_i e^iν_1τ_1+iν_2τ_2 e^τ_1(E_i-E_j)+τ_2(E_j-E_k)
       =(e^-βE_k-e^-βE_i)(1-δ_E_i,E_k)/(iν_1+iν_2+E_i-E_k)(iν_2+E_j-E_k)+e^-βE_j+e^-βE_i/(iν_1+E_i-E_j)(iν_2+E_j-E_k)
    +βδ_E_i,E_kδ_ν_1+ν_2,0/iν_2+E_j-E_k.

For our particular case σ=↓, σ'=↑ and ς=+, so i=|↓>, k=|↑> and j can be either |0> for the first term in (<ref>) or |↑↓> for the second one. Thus

    χ^↓↑+(ν_1,ν_2)=1/Z(f_↓,0,↑(ν_1,ν_2)+f_↓,↑↓,↑(ν_2,ν_1)).

Using (<ref>) and β B≫1 we obtain

    χ^↓↑+(ν_1,ν_2)   =(1-2μ/iν_1+iν_2-2B)1/(iν_1-B-μ)(iν_2-B-μ)
       =-g_ν_1↓g_-ν_2↑(1-Uχ^-+(ν_1+ν_2)).

Let us define the three-point vertex γ'_νω for the spin channel that connects two fermionic propagators and interaction in the same way as in <cit.> with the cut-off on the renormalization parameter α^-+_ω = (1+U^-+χ^-+(ν_1+ν_2))= W^-+_ω/U^-+ instead of χ^-+_ω. The difference between these two definitions is that in the case of γ'_νω the full bosonic propagator of the impurity problem that is attached to the vertex is the full local susceptibility χ_ω, while in the case of γ'_νω vertex function it is equal to the renormalized interaction of impurity problem W^-+_ω in the spin channel. 

Remarkably, the three-point vertex function γ'_νω in the spin-polarized case is equal to unity

    γ'^ -(ν_1,ν_2)= -c^*_↓(ν_1) c^*_↑(-ν_2) S^ +(ν_1+ν_2)/g_ν_1↓g_-ν_2↑ α^-+(ν_1+ν_2)=1,

because in the spin channel the bare interaction is equal to U^+-=-U.
Using the relation between the three- and four-point vertices derived in <cit.>, one gets

    γ'^ ς_νω   =
    α^ ς -1_ω∑_ν'[1 - γ^ ς_νν'ωg_ν'σ g_ν'+ω,σ'] 
    γ'^ ς_νω(1+ W^ ς_ω Π^ ς_ω)    =
    ∑_ν'[1 - γ^ ς_νν'ω g_ν'σ g_ν'+ω,σ'] 
    γ'^ ς_νω   =
    ∑_ν'[1 - (γ^ ς_νν'ω +γ'^ ς_νω W^ ς_ωγ'^ ς_ν'+ω,-ω)g_ν'σ g_ν'+ω,σ'],

where the Hedin expression for the polarization function of impurity problem Π_ω = ∑_νγ'_ν+ω,-ω g_ν,σg_ν+ω,σ' is used. Therefore, when the three-point vertex function γ'_νω is close to unity, the main contribution to the four-point vertex function is given by the following expression 

    γ^ ς_νν'ω≃-γ'^ ς_νω W^ ς_ωγ'^ ς_ν'+ω,-ω.

Transforming back to the definition of the three-point vertex function used in this Letter γ'_νω→γ_νω, one also has to replace the full local bosonic propagator as W_ω→χ_ω. Then, the final expression for the four-point vertex reads

    γ^ ς_νν'ω≃-γ^ ς_νω χ^ ς_ω γ^ ς_ν'+ω,-ω.





§ APPLICATION: EXCHANGE INTERACTION IN THE STRONGLY LOCALIZED REGIME


Here, we calculate the exchange interaction for the Hubbard model in the strongly localized regime t≪U. For this reason, let us find the nonlocal Green's function as the first order correction to the atomic limit solution with respect to the hopping amplitude. Then, using the definition of the Green's function, one gets

    G_ijσ = 1/Z∫ D[c^*,c]  c^*_iσ c^*_jσ  e^-∑_iS^i_ at - ∑_i'j'σ't_i'j'c^*_i'σ' c^*_j'σ' = 
    1/Z∫ D[c^*,c]  c^*_iσ c^*_jσ  e^-∑_iS^i_ at - 
    1/Z∫ D[c^*,c] ∑_i'j't_i'j'c^*_iσ c^*_jσc^*_i'σ' c^*_j'σ'  e^-∑_iS^i_ at.
 
Since the atomic action is purely local, the contribution to the nonlocal Green's function for i≠j is given only by the second term when i=j' and j=i', so σ=σ'. Then, the nonlocal Green's function can be rewritten as   

    G_ijσ = t_ji 1/Z^i_ at∫ D[c^*,c]  c^*_iσ c^*_iσ  e^-S^i_ at×1/Z^j_ at∫ D[c^*,c]  c^*_jσ c^*_jσ  e^-S^j_ at = t g^2_νσ,
    
where g_νσ is the local Green's function of atomic problem. Taking into account the result of Eq. <ref> and that
the difference of the self-energies is equal to
Σ_ν↑ - Σ_ν↓ = U+2B and that 2S^z=1, the exchange interaction reads

    J_ij   = -∑_ν(χ^-1_ω=0 + Σ_ν↑-Σ_ν↓/2S^z)
    G^+_ij,ν↑G^+_ji,ν↓(χ^-1_ω=0 + Σ_ν↑-Σ_ν↓/2S^z)
    = -∑_νtU/(iν-U/2)^2tU/(iν+U/2)^2 
    =  -2t^2/U.
 



§ APPLICATION: ANTIFERROMAGNETIC DIMER

One can also perform an exact diagonlization of a two-site model with the antiferromagnetic ground state 

    Ĥ-μN̂= t c^*_1σc^*_2σ + t c^*_2σc^*_1σ + 
    ∑_i=1,2; σ(Δ^*_iσ
    c^*_iσc^*_iσ
    +Un_i↑n_i↓),

where Δ_↑,↓=± B-μ and the magnetic field B is again considered much larger than the temperature T≡1/β. 
Using the Lehmann representation, one can obtain the nonlocal Green's functions at the low temperatures, i.e. βU≫1 in the strongly-correlated regime t≪U as

    G_ab,σ = 1/Z∑_ij<Ψ_i |c^*_aσ| Ψ_j><Ψ_j |c^*_bσ| Ψ_i>e^-βE_i+e^-βE_j/iν+E_i-E_j
 
Since only the low-lying energy states contribute to the Green's function at low temperatures, because the contribution of higher energy states is exponentially suppressed, we give only relevant energies and (unnormalized) eigenstates below




 E_5 = -U/2+√(B^2+t^2)      
 Ψ_5 = - -B - √(B^2 + t^2)/t |↑↓.↑ > + |↑.↑↓ > 

 E_4 = -U/2+√(B^2+t^2)     
 Ψ_4 = - -B - √(B^2 + t^2)/t |↓. 0> + | 0 .↓ >  

 E_31 = -U/2-√(B^2+t^2)      
 Ψ_31 = - -B + √(B^2 + t^2)/t |↑↓.↑ > + |↑.↑↓ >  

 E_21 = -U/2-√(B^2+t^2)      
 Ψ_21 = - -B + √(B^2 + t^2)/t |↓. 0> + | 0 .↓ >  

 E_11 = -U-2B      
 Ψ_11 = |↑↓. 0> - 2t/U+2B |↓.↑ > - U+2B/t |↑.↓ > + | 0 .↑↓ >  




Then, one gets

    G_ab↑   = 1/Z<Ψ_3 |c^*_a↑| Ψ_1><Ψ_1 |c^*_b↑| Ψ_3> e^-βE_3+e^-βE_1/iν+E_3-E_1 + 1/Z<Ψ_1 |c^*_a↑| Ψ_2><Ψ_2 |c^*_b↑| Ψ_1> e^-βE_1+e^-βE_2/iν+E_1-E_2
        + 1/Z<Ψ_5 |c^*_a↑| Ψ_1><Ψ_1 |c^*_b↑| Ψ_5> e^-βE_5+e^-βE_1/iν+E_5-E_1 + 1/Z<Ψ_1 |c^*_a↑| Ψ_4><Ψ_4 |c^*_b↑| Ψ_1> e^-βE_1+e^-βE_4/iν+E_1-E_4
       = 1/Z1/N^2_1N^2_3(5)(- -B ±√(B^2 + t^2)/t < ↑↓.↑| + < ↑.↑↓|)(- 2t/U+2B |↑↓.↑ > + |↑.↑↓ >)×
       ×(- <↑↓.↑| + U+2B/t <↑.↑↓|)(- -B ±√(B^2 + t^2)/t |↑↓.↑ > + |↑.↑↓ >)
    e^β(U/2±√(B^2+t^2))+e^β(U+2B)/iν-(U/2±√(B^2+t^2))+(U+2B)
        + 1/Z1/N^2_1N^2_2(4)(< ↓. 0 | - U+2B/t <0 .↓|)(- -B ±√(B^2 + t^2)/t |↓. 0> + | 0 .↓ >)×
       ×(- -B ±√(B^2 + t^2)/t <↓. 0 |  +<0 .↓|) 
    (- 2t/U+2B |↓. 0> + | 0 .↓ >)
    e^β(U/2±√(B^2+t^2))+e^β(U+2B)/iν+(U/2±√(B^2+t^2))-(U+2B).
 
Note that all eigenstates Ψ_i were normalized as 1/N_iΨ_i. Simplifying the previous equation one gets

    G_ab↑   = 1/Z1/N^2_1N^2_3(5)(-2B ± 2√(B^2 + t^2)/U+2B+1)(-B ±√(B^2 + t^2)/t +U+2B/t)
    e^β(U/2±√(B^2+t^2))+e^β(U+2B)/iν+U/2+2B∓√(B^2+t^2)
        + 1/Z1/N^2_1N^2_2(4)(- -B ±√(B^2 + t^2)/t- U+2B/t)(-B ±√(B^2 + t^2)/t2t/U+2B + 1)
    e^β(U/2±√(B^2+t^2))+e^β(U+2B)/iν-U/2-2B±√(B^2+t^2).
 
Taking into account that

    Z    = 4(1+e^βU/2coshβB)≃e^β(U+2B)
    
    N^2_1   = 2 + (2t/U+2B)^2 + (U+2B/t)^2 ≃(U+2B/t)^2 
    
    N^2_2(4)   = N^2_3(5) = 1 + (-B±√(B^2+t^2)/t)^2

one can finally get

    G_ab↑   = t^2/(U+2B)^21/1 + (-B±√(B^2+t^2)/t)^2U ± 2√(B^2 + t^2)/U+2BU+B ±√(B^2 + t^2)/t1/iν+U/2+2B∓√(B^2+t^2)
        - t^2/(U+2B)^21/1 + (-B±√(B^2+t^2)/t)^2U+B ±√(B^2 + t^2)/tU ± 2√(B^2 + t^2)/U+2B1/iν-U/2-2B±√(B^2+t^2).
 
If U≫B and U≫t, we get

    G_ab↑   = ∑_±t/U1/1 + (-B±√(B^2+t^2)/t)^2[
    1/iν+U/2 - 1/iν-U/2] = t/U[
    1/iν+U/2 - 1/iν-U/2] .
 
The same result can be found for 

    G_ba↓   = t/U[
    1/iν+U/2 - 1/iν-U/2] .
 
Therefore, the exchange interaction reads

    J_ab= 
    
    -U^2/β∑_ν G_ab↑G_ba↓ = t^2/2π∫_-∞^+∞2 dx/(x-iU/2)(x+iU/2) = - 2t^2/U.
 






