Spin Seebeck coefficient and spin-thermal diffusion in the two-dimensional Hubbard model

We investigate the spin Seebeck coefficient $S_s$ in the square lattice Hubbard model at high temperatures of relevance to cold-atom measurements. We solve the model with the finite-temperature Lanczos and with the dynamical mean-field theory methods and find they give similar results in the considered regime. $S_s$ exceeds the atomic 'Heikes' estimates and the Kelvin entropic estimates drastically. We analyze the behavior in terms of a mapping onto the problem of a doped attractive model and derive an approximate expression that allows relating the enhancement of $S_s$ to distinct scattering of the spin-majority and the spin-minority excitations. Our analysis reveals the limitations of entropic interpretations of Seebeck coefficient even in the high-temperature regime. Large values of $S_s$ could be observed on optical lattices. We also calculate the full diffusion matrix. We quantify the spin-thermal diffusion, that is, the extent of the mixing between the spin and the thermal diffusion and discuss the results in the context of recent measurements of the spin-diffusion constant in cold atoms.


I. INTRODUCTION
Cold-atom systems on optical lattices provide a novel lens on poorly understood transport regimes of correlated electrons.
They can realize the Hubbard model-the standard model of correlated electrons that interact with an onsite repulsion U and move on the lattice (hopping t)-without real world complications, such as lattice vibrations and disorder. Hence one can directly and quantitatively compare the outcome of the experiment to those of the numerical solutions of the Hubbard model [1,2]. Such a cross-verification turned very successful in the measurements of the charge diffusivity [1]. Besides providing an important mutual benchmark of the methods it led to a quantification of the vertex corrections [3]. Intriguingly, related measurements of spin diffusivity revealed a disagreement between the numerical methods and the experiment [2]. An independent numerical investigation [4] confirmed the results of the theory but disagreed with the experiment.
This disagreement thus calls for a close inspection of the underlying assumptions of the experimental analysis. One of the assumptions is that the spin-thermoelectric effect is unimportant. This holds strictly at a vanishing magnetization, but in the actual experiment this condition was only approximately met as some spin imbalance m z = (n ↑ − n ↓ )/2 is seen in the measurements: |m z | 0.05 [2]. The strength of the spinthermoelectric effect is quantified by the spin-Seebeck coefficient S s given by the ratio of the magnetic field and thermal gradient at the condition of vanishing spin current, S s = ∇B/∇T | js=0 . S s is a quantity that is relevant for spintronics applications [5][6][7][8] but has to our knowledge not been discussed for the Hubbard model (surprisingly, as this is the paradigmatic model of corre-lated electrons). The intention of our work is to fill this gap, establish how large S s is in the high-temperature regime and use that knowledge discuss whether the measurements of spin-diffusion in cold atoms could be influenced by the spin-thermoelectric effects.
Some intuition concerning S s could be expected from considerations that relate the ordinary charge Seebeck coefficient S c to thermodynamic quantities, such as the high-temperature Heikes limit S H c = µ/T or the Kelvin formula S K c = dµ/dT that relate the Seebeck coefficient to the temperature dependence of the chemical potential [9][10][11][12]. Namely, it was demonstrated that often at high temperatures the Kelvin formula describes the Seebeck coefficient well [11] (but some exceptions to this were also noted [13]).
To apply this intuition to S s one can use a mapping that relates the spin transport in a magnetized repulsive Hubbard model to the charge transport in a doped attractive Hubbard model [14][15][16][17][18][19][20][21]. This mapping proceeds via a particle-hole transformation on particles of only one spin, e.g, ↓, with c i,↓ → (−1) i c † i,↓ and results in an interchange of spin and charge degrees of freedom, explicitly n ↑ − n ↓ → n ↑ + n ↓ , n ↑ n ↓ → −n ↑ n ↓ . Accordingly, Hubbard repulsion goes to attraction, U n ↑ n ↓ → −U n ↑ n ↓ . Via this mapping, one can relate the spin Heikes estimate for spin Seebeck coefficient S H s = B/T | mz=const (or Kelvin estimate S K s = dB/dT | mz=const ) to the corresponding charge Heikes and Kelvin estimates for a model with an opposite sign of repulsion (that is, an attractive model for the case of interest here). Actually, by exploiting the mapping one can use the results from the literature [22] and obtain S s = 8k B m z at a high temperature (T > U ) and S s = 4k B m z at a lower temperature (T < U ). From these estimates-that one is inclined to trust, especially in the high-temperature regime T > t pertinent to cold atom measurements-one expects only small values of the spin Seebeck coefficient S s k B , since m z 1. In this paper we show that this reasoning is incorrect. We calculate S s for a square lattice Hubbard model at high-T using the finite temperature Lanczos method (FTLM) and the dynamical mean-field theory (DMFT) approaches and find that it strongly exceeds the bounds just discussed. S s violates the thermodynamic expectations even in the high-temperature regime, an unexpected finding based on what was previously known for the ordinary thermoelectric effect. Large values of S s call for a reexamination of the possible spin-thermoelectric effects in cold atom measurements of spin-diffusion. To estimate those, we calculated the full diffusion matrix D. In general, the eigenvalues of D deviate from those found in the absence of spin-thermoelectric effect. Interestingly, we find the deviations are related to the difference between the actual value of spin-Seebeck coefficient and its thermodynamic Kelvin approximate, S s − S K s . We discuss why in spite of this difference being sizable (it exceeds k B at large U ) the final influence on the measured spin-diffusion for moderate magnetization is unimportant. We note that large values of Seebeck coefficient for the attractive model were earlier found in the DMFT [21] but were not compared with the thermodynamic estimates and the importance of those results for the spinthermoelectric response was not discussed. The remainder of the paper is structured as follows. In Sec. II we specify the model, the methods, and the notation. In Sec. III we show our main results for the spin-Seebeck coefficient. In Sec. IV we describe the DMFT calculation of transport and in Sec. V we exploit it in conjunction with the mapping to an attractive model to interpret our results. In Sec. VI we investigate the influence of the spin-thermoelectric effect for the spindiffusion. In Sec. VII we give our conclusions. The Appendix discusses the behavior of spin-Seebeck coefficient for a phenomenological ansatz spectral function.

II. MODEL AND METHOD
We study the square lattice Hubbard model, with t being the hopping between the nearest neighbors. We take = k B = e = gµ B = 1. We likewise take lattice spacing a = 1. We use t as the energy unit. We solve the Hamiltonian with FTLM [23][24][25] on a N = 4 × 4 cluster and in the thermodynamic limit with the DMFT [26] (that is, we solve the Hamiltonian Eq. (1) in a local approximation). The DMFT equations are solved using the numerical-renormalization group (NRG) [27] in the NRG-Ljubljana implementation [28] as the impurity solver.  Figure 1(a) displays S s as a function of temperature for U = 10 evaluated with the FTLM (full, thick) and the DMFT (symbols). At highest temperatures, S s approaches the high-temperature Heikes value, 2 log(1 + 2m z )/(1 − 2m z ) ≈ 8m z , which is the expected behavior. The Kelvin estimate from dB/dT | mz evaluated using the FTLM (thin; DMFT gives similar results) agree with the Kubo evaluation in this regime. On lowering the temperature, surprisingly, instead of diminishing in amplitude as suggested by the Heikes value correspond-

III. SPIN-SEEBECK COEFFICIENT
, S s increases and reaches a maximum value well above the Heikes estimate, and only drops consequently at lower T . This behavior with a substantial increase of the spinthermoelectric coefficient above the high-temperature and thermodynamic estimates becomes even more pronounced for larger U , as displayed in Fig. 1(b). The magnitude of the peak diminishes with decreasing magnetization but increases with increasing U . It is, as we show below, proportional to m z U/T . Throughout the considered regime, the temperatures are high (T > t) and one cannot attribute the deviations from the thermodynamic estimates to the occurrence of coherent transport or to a proximity to a magnetically ordered regime. The observed deviations are in stark contrast to the behavior of the charge thermopower that at temperatures T > t does follow the thermodynamic estimates [11].

IV. DMFT DESCRIPTION OF TRANSPORT
In order to understand this behavior it is convenient to discuss the transport properties within the DMFT approach. The DMFT expresses the transport coefficients in terms of the transport function [29] where v k is the band velocity: v k = d /dk x with k the band energy. A kσ is the spectral function at momentum k and spin σ. The charge Seebeck S c and the spin Seebeck S s coefficients are, respectively, (3) As seen in Fig. 1(a), the DMFT results are similar, although not identical to the FTLM ones. To what extent the differences are technical (it is quite challenging to converge the DMFT calculations in this regime as discussed in Ref. [21]) or physical, such as emanating in non-local fluctuations and/or the vertex corrections neglected in DMFT is a question that goes beyond the scope of the present paper. It is likely that at low T the vertex corrections become more important as the DMFT calculation gives an insulator with a spin gap whereas the actual behavior is that of a spin conductor described by a Heisenberg model. For our purpose we will ignore these differences between the two methods and exploit the more transparent DMFT formulation of transport to interpret the FTLM results.

V. MAPPING TO THE ATTRACTIVE MODEL
It is convenient to analyze the results in terms of the mapping of the spin to the charge degrees of freedom for an attractive model, that is S s (U ) = 2S c (−U ), with the spin polarization 2m z → δ becoming the charge doping with factors of 2 occurring due to the definition of spin. Figure 2(a) presents the local spectral function k A k (ω) of the doped attractive model at T = 6, for two values of attraction −U = 10, 20. As discussed in earlier studies of the attractive model [16,17,21,30], the spectral function consists of two peaks, which are as for the repulsive Hubbard model centered at −µ and −µ + U , with µ ∼ U/2. The important distinction between the doped attractive and the doped repulsive model is in the behavior of the chemical potential with temperature that can at high-T be most simply obtained from a grand-canonical treatment of the atomic problem. There, average electron occupancy can be evaluated from n = 2 exp(βµ) + 2 exp(−β(U − 2µ)) 1 + 2 exp(βµ) + exp(−β(U − 2µ)) In the attractive model terms that include exp(−βU ) grow at low temperatures and should be retained. Hence one obtains The fact that µ ∼ U/2 at low T represents a crucial difference with respect to the repulsive case causes the spectral function to be gapped at a nonvanishing doping. Namely, the spectral function consists of two peaks displaced by δµ from ±U/2 as shown on Fig 2(a) for two values of U . Because at large |U | the gap is well developed, the lower and upper Hubbard bands must have unequal spectral weight to yield a finite doping. Figure 2(d) presents the transport function. One sees that this exhibits the two Hubbard bands and is overall similar to the density of states. There is however an important difference: Because Φ contains A 2 k , the weights of the upper and the lower Hubbard band parts are affected by the amplitude of scattering, as given by self-energy depicted in Figs. 2(b) and 2(c) for |U | = 20, 10, respectively. When the spectral function is a sharply peaked function, A 2 k ∼ A k /(2πΓ k ) with Γ k = −ImΣ(ω k ), and hence the transport function is modulated by the value of the self-energy at the peak frequency ω k . In passing we note that this finding can be used to connect the bubble expression to the Boltzmann calculation, see Refs. [31,32] for a recent discussion. In the Appendix A we take advantage of the simple twopeaked structure of the transport function and find an expression where the first term grows as U/T and is proportional to a coefficient that is expressed in terms of the effective weights of the positive (negative) frequency peaks of the transport functionφ ± . This coefficient (that is found to be approximately given by doping, see Appendix A) grows with the difference in the scattering between the holes and the electrons. When this difference does not vanish, the behavior of the Seebeck coefficient differs from that of the thermodynamic expectations given by the second term in Eq. (6). Let us relate this discussion to the repulsive case. At the particle-hole symmetry, the spectral functions of the repulsive and attractive case for respectively s =↑, ↓ are related by A k↑,↓ (ω)| U >0 = A k↑,↓ (±ω)| U <0 . That is, taking advantage of the mapping, Fig. 2 depicts s =↑ components of the spectral function, the self-energy and the transport function, and the s =↓ components can be obtained by ω → −ω. The different scattering between the electrons and the holes in the attractive model thus relates to a different scattering between the spin majority and the spin minority carriers in the repulsive model, and this in turn leads to S s ∼ m z U/T behavior that explains the enhancement over the thermodynamic estimates seen in Fig. 1.

VI. DIFFUSION MATRIX, DIFFUSION EIGENVALUES, AND RELEVANCE FOR EXPERIMENT
Predicted large values of spin-Seebeck coefficient at large T and the increase of the peak value with increasing U (or the corresponding behavior of the charge-Seebeck coefficient in the attractive case) could be tested in future cold-atom experiments. In the introduction we also raised a possibility that the existing measurements of spin-diffusion would be affected by spinthermoelectric effects, which could account for values of the spin-diffusion and the spin-conductivity that were found to be larger than theoretically expected. Let us try to estimate the influence of thermoelectric effects on spin-conductivity using a hand-waving argument. At finite magnetization, gradients of magnetic field are accompanied by a gradient of energy, and, assuming thermalization, a gradient of temperature, ∇T ∼ ∇E/c = m z ∇B/c, with c as the spe- Importantly, because S s and m z are of equal sign, this estimate anticipates the spin conductivity is actually reduced compared to the case where spin thermoelectricity is neglected. In order to make this discussion more precise one must consider a generalization of the Nernst-Einstein relation to a matrix formulation allowing for a mixed response [33], where the off-diagonal entries involve the mixed transport coefficient L sq and the thermoelectric susceptibility ξ = −∂ 2 f /∂B∂T where f is the free energy density. The diffusion constant in matrix form reads D = −LA −1 , where {j q , j s } = L{∇T, ∇B} defines the conductivity matrix L and the heat-magnetization susceptibility matrix A is defined by {T ∇s, ∇m z } = A{∇T, ∇B} [34]. The diffusion eigenmodes that are involved in general non vanishing spin and heat components are obtained by diagonalizing the matrix D. We denote the diffusion eigenvalue whose mode contains a predominantly spin (heat) component by D − (D + ), respectively. These are shown in Fig. 3(a) as a function of temperature for U = 10, m z = 0.05 and are compared to the bare spin-diffusion constant D s = σ s /χ s and bare heat diffusion constant D q = κ/c, where we use term "bare" to indicate that the spinthermal mixing is neglected. One sees that, consistent with the hand-waving discussion above, the influence of the spin-thermoelectric effects is only moderate, the diffusion eigenvalues are close to the values obtained when there is no mixing. The behavior of D ∓ being smaller (larger) than the corresponding bare diffusion is that of level repulsion.
The relatively small admixing occurs because the geometric mean of the off-diagonal elements D sq D qs is significantly smaller ( 10%) than D qq − D ss . When any of the off-diagonal elements vanish, the spinthermoelectric effect on diffusion vanishes. We inspect more closely D sq . It can be rewritten as D sq = −D s ξ+Ssχs c , which expresses spin-thermal diffusion, i.e. spin-thermoelectric mixing between the spin and the thermal diffusion. Note that S s is multiplied by χ s /c, which becomes large at T → 0. This hints at a rich behavior at low temperatures and should be explored in future work. Using the Kelvin formula S K s = dB/dT | mz , one can further rewrite: D sq = −D s χs c (S s − S K s ). Interestingly, when the spin Seebeck coefficient is given by the Kelvin estimate, D sq = 0, which leads to D − and D + equal D ss and D qq . In this case the spin modulation decays with pure spin diffusion constant D ss = σ s /χ s and has an admixed thermal (heat) component. An initial spin modulation therefore also induces heat currents. On contrary, D + corresponds to pure heat diffusion decaying with the bare heat diffusion constant D qq = κ/c.
In experiment, an initial gradient of magnetization is imposed and D sq . Could larger values of the experimentally inferred spin-diffusion constant occur because there is a significant contribution of the heat eigenmode in the initial state that decays faster (see larger values of D q and D + in Fig. 3)? The components of eigenvectors v ± are shown on Fig. 3(b). At T = 3.1 we find the spin dominated eigenvector v − = {−0.53t, 0.84}, i.e. it contains a significant component of heat current that only increases with temperature. (In this expression we reintroduced t to indicate that the two components of v are in different units; because t is a natural unit for the energy it is meaningful to compare the numerical values of the two components setting t = 1). At the same temperature, the heat dominated v + = {0.999t, −0.02} is mostly single component. One can first assume that the initial state is given by a pure magnetization gradient. Expanding this profile in terms of v ± we find only a small part of magnetization is contained in v + . To be specific, at T = 3.1, v + contains ∼ 2% of the initial spin modulation (the relative weight of v + is sizable, but it carries only a small magnetization). Hence, the faster decay of the v + cannot importantly affect the evolution of the magnetization. What if the heat gradient component is also initially present? Using the estimate δT ∼ mz χs δm z (describing the situation where the magnetic field responsible for δm z is switched off and excess energy is instantly converted into heat), we find that v + is more prominent in the initial state, but still accounts for ∼ 4% of the initial magnetization modulation. In both cases, the majority of spin diffusion is therefore governed by D − ∼ D s .

VII. CONCLUSIONS
In summary, we calculated the spin-Seebeck coefficient in the Hubbard model and discovered a rich behavior with temperature. The spin-Seebeck coefficient significantly exceeds the entropic estimate. This occurs due to the unequal scattering of spin minority and spin majority carriers which gives rise to an ∝ m z U/T dependence that adds up to the Heikes' estimate. This is a striking demonstration of the breakdown of the entropic interpretation of the thermopower in a high-temperature regime where a priori one would trust it the most. Our predictions for a large spin-thermoelectric effect could be tested on optical lattices. We also calculated the diffusion matrix eigenvalues and estimated the influence of spin thermoelectricity in the existing measurements of the spin-diffusion [2]. We found this influence to be moderate and insufficient to explain the discrepancy between the experiment and the theory. Possible directions for future research include simulating explicitly the time dependence in such experiments and the study of possible nonlinear effects.

Appendix A: Phenomenological analysis
Here we take advantage of the known general shape of the transport function to obtain an approximate simple expression for the Seebeck coefficient. The temperatures we consider in our simulations (and that pertain to the cold atom experiments, which motivate our investigation) are large. In Ref. [2] the estimated entropy is 1.1, which pertains to T ≈ 0.3|U | (in this Appendix we will phrase the discussion in terms of the attractive Hubbard model), comparable to the bandwidth. Hence the derivative of the Fermi function df /dω does not change substantially over each of the Hubbard bands and hence in the transport integrals, Eq. (3), one can approximate the transport function by two δ−peaks as with different weights φ − and φ + for the negative and positive frequency peaks, respectively. Note that U < 0 in this case. One can evaluate the Seebeck coefficient using this ansatz transport function. Taking also into account that δµ/T is small, one can expand the derivative of the Fermi function, that is −df /dω(±U/2 − δµ) = −df /dω(U/2)(1 ± t 0 δµ/T ) where we define t 0 = tanh(U/4T ). Introducing effective weights (modified from φ due to df /dω)φ ± = φ ± (1 ∓ t 0 δµ/T ), one ob- tains a simple expression for the Seebeck coefficient The high-T Seebeck coefficient for the attractive model thus has a Heikes' term (second term of this expression, predicted by Chaikin and Beni [22]) but crucially also the first term S 1 , that is proportional to U/2T and the difference between the effective weights of the peaks of the transport function. Whenever this difference (that, as we discuss next, is due to a different scattering of electrons and holes) does not vanish, the Seebeck coefficient cannot be interpreted in terms of the entropic considerations alone. This explains large values of spin-Seebeck coefficient seen in numerical results of Fig. 1.
Is scattering really important? Would not the effective weights differ already due to the different weights of the corresponding Hubbard band weights in the density of states? If this were the case, one would have φ ± ∝ (1 ∓ δ) (simply from the considerations of occupancy). Incidentally, the influence of δµ just cancels at small T . Namely, for small doping δµ/T ≈ δ. As T /|U | becomes small, t 0 → −1. Hence one hasφ ± ≈ φ ± (1 ± δ), and φ + =φ − . In this limit S 1 would vanish. One needs the to take the scattering into account to understand the occurrence of deviation from entropic estimates. We plot the weight φ − , obtained from the integral of the transport function over negative frequencies normalized such that the total integral is 2 on Fig. 4(a). One sees that in most of the studied temperature range φ − is close to the value expected from the dependence (1+δ) 2 .
Only at smaller temperatures the weight φ − decreases and actually approaches a smaller value (1 + δ). In Fig. 4(b) we compare numerical results for S c with the result of Eq. (A2) where we set δµ/T to the hightemperature Heikes value we approximate the transport function weights with φ ± ∝ (1 ∓ δ) z with z = 1 (blue) and with z = 2 (green). At small temperatures these lead to a behavior S 1 = −U/2T (z − 1)δ (where corrections of order δ 2 and higher are ignored). For z = 1, which corresponds to taking into account just the different number of carriers, S 1 vanishes and the strong increase of S c seen in numerical simulations is not reproduced. One needs to take into account also different scattering of carriers (as in this approximation embodied for z = 2), too.