Quantum Monte Carlo simulation for the spin-drag conductance of the Hubbard model

In the situation of two electro-statically coupled conductors a current in one conductor may induce a current in the other one. We will study this phenomenon, called Coulomb drag, in the Hubbard chain where the two ``conductors'' are given by fermions with different spin orientation. With the aid of a Monte Carlo (MC) approach which we presented in a recent paper we calculate the Transconductance in different variants of the Hubbard chain (with/without impurity and additional [long-ranged] interactions) for different fillings.


I. INTRODUCTION
The Coulomb-drag effect describes how two conductors (only coupled by the Coulomb force) may influence each other. Since the Coulomb repulsion is relatively small, a sizeable effect will only arise when the two conductors are very close to each other. This condition can be met in mesoscopic systems-where with the advent of new technologies (e.g., carbon nanotubes) the problem of Coulomb drag attracted more and more attention (e.g., Refs. 1,2,3,4)-or in the spin-drag problem. For instead of considering two conductors one may look at drag effects between different fermion species, e.g., fermions with different spin orientation. Since fermions with different spins are not spatially separated, there is a large Coulomb force between them which can lead to all kinds of correlation effects, e.g., a drag effect. In the last years, the interest in spin-dependent transport increased. One key problem is the generation of a spin-polarized current, i.e., a current where only fermions with one of the two spin orientations flow. In this context it is important to keep in mind that the spin-polarized current may affect the fermions with the opposite spin orientation. Hence, the drag effect may play here a crucial rôle even though it is experimentally not directly accessible. (This is because the driving potentials are in general not spin dependent.) For this "spin-drag" problem the trans-resistivity of higher dimensional systems has been investigated in previous publications, e.g., using the Boltzmann equation 5 or the random-phase approximation. 6 In this paper we focus on the Transconductance for the Hubbard model. While most authors used a bosonization approach to compute the conductance 7,8 we will use here, for the study of the Transconductance, a Monte Carlo method which we introduced in a recent paper. 9 The strategy we followed there was to map our fermion system via the Jordan-Wigner transformation to a spin system which can be analyzed by efficient though standard Monte Carlo techniques. We will now extend this method to the onedimensional Hubbard model concentrating on the question how a spin-polarized current (driven by a voltage drop which is assumed to be spin dependent) affects the fermions with opposite spin orientation.
In the Section II we present the model and give some central definitions for the spin-drag problem. The Sec. III contains the technical details on the subject of the MC simulations. The MC method of our choice 9,10 was a variant of the Stochastic Series Expansion (SSE) as introduced in Refs. 11,12,13,14. This method allows an investigation of the one-dimensional Hubbard model. 15,16 In the following Sec. IV A we present our results for the standard Hubbard model and compare with analytical predictions from bosonization theory. To obtain a spin-polarized current we add in Sec. IV B an impurity to the system, which acts like a combination of a onesite chemical potential and a one-site magnetic field. We show that such a "magnetic" impurity can produce the desired spin-polarized current.
The Hubbard model can also be mapped to a spinlessfermion ladder. Hence, our results may also describe that situation, but there one might argue that the very specific modeling of the Coulomb (on-site) interaction in the Hubbard model is unrealistic (and may differ from other approaches, e.g., Refs. 1,2,3,4). We therefore discuss two variants (with additional interaction terms) of our model. First, we will discuss in Sec. IV C a situation, where fermions with different spins live on different sites. The full system has the geometry of a zig-zag chain. Second, we show in the appendix that a spin-polarized interaction leads to equal Cis-and Transconductance. This is similar to the "absolute" drag result found, e.g., in Ref.

II. DEFINITION OF THE SPIN DRAG IN THE HUBBARD MODEL
Our model Hamiltonian is the standard Hubbard model (with N sites or atoms) where n n,σ is the occupation number of fermions with spin σ at site n, and c ( †) n,σ is the corresponding annihilation (creation) operator. To perform our transport calculations, we will use the approach from Ref. 9. Since the hopping term does not connect fermions with different spins, it is natural to consider current and potential operators for each spin orientation separately. Following Ref. 9 the potential operators read then: (e being the charge unit and x being the position of the voltage drop) The conductance (of a spinless-fermion chain) is the linear response of one potential operator to another; therefore the explicit form of the current operators is not needed here. 9 As we have two potentials, we can define four transport quantities (conductances) g ij which describe the (linear) response of P i to P j where i, j ∈ {↑, ↓ }. Further details on how to evaluate the g ij 's are to be found in section III. For the moment we will discuss only symmetric models, i.e., we have spin-rotational invariance (the only asymmetric model that we will discuss appears in Sec. IV B); thus, we end up with only two distinct quantities. We call g c := g ↓↓ = g ↑↑ the Cisconductance and g t := g ↓↑ = g ↑↓ the Transconductance.
The naming conventions come from the physical interpretation of these coefficients which is the following: If we switch on at a certain time a (supposedly) spin-polarized potential which acts only on spin-up fermions (i.e., we add a time-dependent perturbation of the form with V being the voltage amplitude and θ the Heavysidestep function) we will find a current of spin-up fermions This is the drive current governed by the Cisconductance, but there will also be a current of spin-down fermions the drag current (governed by the Transconductance). The latter may be nonzero, even though the spin-down fermions do not feel the applied potential.
The situation of a nonvanishing Transconductance (or drag current) is called Coulomb drag. This problem has been studied, e.g., for coupled spinless-fermion systems by bosonization in Refs. 1,2,3 and to second-order perturbation theory in Ref. 4.
Normally, (since spin-polarized potentials are not available) one is only interested in the (full) conductance of the Hubbard model, where both fermion species feel the same potential, and the full current is the sum of the currents of the spin-up and spin-down fermions. As is straightforward to see, the Cis-and Transconductance give us directly the conductance of the Hubbard model via the relation g Hubbard = 2(g c + g t ).

III. DESCRIPTION OF THE MONTE CARLO METHOD
Before starting with Monte Carlo we have to cast our Hamiltonian in a convenient form.
Using the Jordan-Wigner transform the Hubbard model can be mapped to a spin ladder. To each occupation operator we introduce a spin operator, i.e., we replace n n,↑ → S z 2n + 1/2 and n n,↓ → S z 2n+1 + 1/2. If we express the Hubbard Hamiltonian with those new operators we obtain the following spin ladder (with 2N sites): where the sites with even number represent spin-up fermions and the sites with odd numbers, spin-down fermions. (For the Hubbard model one has to put J z = 0 = U ′ . These parameters are used to model the spin-polarized interaction from the appendix and the zigzag chain from Sec. IV C; see bottom half of Fig. 1. The sites 2n and 2n + 1 in the ladder represent therefore one atom of the Hubbard model and interact via an Ising interaction representing the Coulomb force (see Fig. 1). The hopping amplitudes satisfy J x = 2t, and the strength of the magnetic field B in Eq. (2) is obtained from the chemical potential via the relation Half filling corresponds therefore to B = 0. The two potential operators P ↑,↓ x introduced in the previous section can be related to potential operators for the two chains (e being the charge unit) Then we may obtain the four conductances g ij introduced in the previous section by computing (i, j ∈ {↓, ↑}) at the Matsubara frequencies ω M = 2πM (β ) −1 , M ∈ N, and then extrapolating to ω = 0. (The extrapolated value should not depend on x or y. 9 We chose x = N/2 = y − 1.) For the extrapolation from g(ω M ) to g(ω = 0) we will use a quadratic fit from the first three Matsubara frequencies. [  plaquettes (see Fig. 1). The following operators belong to the plaquette P (containing interactions between the sites 2n, 2n + 1, 2n + 2, and 2n + 3): The set h consists of all h (t) P for all plaquettes P and all t = 1, . . . , 5.
The heart of the SCSE program is the so called loop update, where a spin flip of a subset (loop) of all spin variables is proposed. Since the sites with even and odd numbers form two chains, which are only coupled by a z-z interaction term, we find that the set of spin variables that will be flipped in the loop update belongs entirely to one of the chains. Therefore, we can view the new algorithm as making loop updates for each chain separately. During a loop update for one chain the spin variables of the other chain remain fixed. The consequence of this is that, if we update, e.g., the even chain, then operators with superscript i = 3, 4 can be neglected (are irrelevant for the loop construction), and the coupling terms (between the chains) reduce to magnetic field terms (for the even chain).
It is however advantageous to consider another variant of the loop update. The construction is similar to the first variant, but now we propose spin flips for both chains, i.e., the spin variables belonging to sites 2n and 2n + 1 are flipped simultaneously. This may be viewed as a construction of two parallel loops-one for each chain. Since the two loops must be parallel, the number of possible transitions between different plaquette states is reduced. This may lead to a less efficient algorithm, 13 but one should note that this parallel-loop update becomes deterministic for the case of B = 0 and hence enhances the efficiency of the algorithm (at least for this situation) considerably.

IV. NUMERICAL RESULTS
A. Transconductance in the Hubbard model

1.
Comparison with the Hubbard model In bosonization theory the Hubbard model is described by two boson fields Φ ↑,↓ representing the degrees of freedom of different spin orientations. The current operators for the spin sectors are then given by The conductance can be written in the form of a currentcurrent correlator 7 and may be evaluated in terms of the Luttinger-liquid parameters K ρ,σ of the charge and spin The result is (using the linearity of the correlator and the results from Ref. 17) When µ = U/2 we are at half filling, where Umklapp processes are responsible for a gap in the system. 18 The (charge) gap ∆(U ) depends on the Hubbard repulsion U and is finite for all U .

Numerical simulations
We present now Monte Carlo results for the Transconductance in the Hubbard model Eq.
We performed simulations for two different chemical potentials: First, µ = U/2 corresponding to half filling and second µ = 0. In the latter case the system is no longer at half filling, but has a U -dependent filling, which is shown in Fig. 11 of the appendix (in the large-U limit the system reaches quarter filling).
We show g c and g t as a function of U for the two different µ's in Fig. 2. The figure shows that the Coulomb drag is very sensitive to a change in chemical potential.
Let us first look at the half-filled case (Fig. 2). If U is very large the Coulomb repulsion acts as an effective projection to those configurations satisfying P ↓ = −P ↑ . This implies that g c + g t → 0 as U → ∞. This contemplation is in accordance with Fig. 2. We should actually expect from Eq. (4) that for T = 0 we have g c + g t = K ρ = 0 because of the charge gap ∆(U ) > 0 for all U > 0 (cf. Ref. 18). This should lead to a discontinuous jump at U = 0, because without the Coulomb force evidently g t = 0 and g c is the conductance of uncoupled spinless-fermion chains from Ref. 17. Here we emphasize that our method is a finite-temperature method, which means that the conductances calculated by us interpolate smoothly between the values for U = 0 and U = ∞. The crossover is expected to take place at that interaction value U T which satisfies ∆(U T ) = k B T . It is therefore interesting to see how g c,t scale with temperature. However, our method gives only access to the low-T regime, 9 such that we will compare here results for only two different temperatures, k B T = 0.01t and k B T = 0.02t. In the two subsequent simulations for the Cis-/Transconductance we did not find any difference at all. This implies only a weak temperature dependence (at low T ) for the interaction U T which governs the crossover. Since ∆(U ) is known from analytical results, 19 we can calculate the two crossover interactions-defined by ∆(U 0.02 ) = 0.02t and ∆(U 0.01 ) = 0.01t-finding that U 0.02 ≈ 1.25t and U 0.01 ≈ 1.12t do not differ much, as expected (they are also both indicated by arrows in Fig.  2). Another important consequence of g c + g t → 0 is that the signs of the Cis-and Transconductance are opposite or-in terms of the spin-up and spin-down currents-that the induced drag current flows in the opposite direction of the drive current. 20 Now we will turn to the spin sector. We have K σ = 1 by spin-rotational invariance of the Hubbard model 18 implying [see Eq. (4)] g c − g t = 1 for all U which is very well satisfied by Fig. 2.
Putting the two results for g c ± g t together, we obtain g c = 0.5 e 2 /h = −g t valid at high U . This large-U limit of g c may also be computed in second-order perturbation theory. In this approximation the Hubbard model can be mapped to a Heisenberg model. The operator [on the Hilbert space of the original Hamiltonian Eq. (1)] P ↑ x , which is effectively equal to −P ↓ x , is identified with the operator (on the Hilbert space of the effective Heisenberg model) P x = n>x T z n . [ T z n = (n n,↑ − n n,↓ )/2 is the spin operator for fermions; here denoted by T to avoid confusion with the spin operators appearing in Sec. III.] Applying the results from Ref. 9 the computation of g c reduces then to the computation of the spin conductance of the Heisenberg model, which equals one half in units e 2 /h. In the case of zero µ (again Fig. 2) where the system is away from half filling there is no charge gap. Hence, K ρ is finite, and so Eq. (4) tells us that g c +g t = K ρ does not decay with U . [Here we note that g c + g t agrees (within error bars) with the values for K ρ available in Ref. 18.] We still have K σ = 1, which leads to g c − g t = 1 for all U (again very well satisfied by the figure).
Finally, we consider the large-U limit. Inserting K ρ (U = ∞) = 0.5 and K σ ≡ 1 from Ref. 18 in Eq. (4) yields g c = 0.75 and g t = −0.25 (units e 2 /h). These results are in accordance with the figure. (Note that the statistical error increases with U such that we cannot compute g(U ) for sufficiently high U in order to extract the large-U limit accurately.)

B. Magnetic Impurity
In this subsection we will study the influence of an impurity. The modeling of the impurity follows Ref. 7, but for the spin-drag problem it is natural to consider a spin-dependent impurity, as we will do here. We extend therefore our Hamiltonian in the following way H = H Hubb + B Imp n N/2,↑ , i.e., we introduce an (impurity) potential at exactly one central site (which acts only on one spin orientation, see Fig. 3). Fig. 4 shows Cisconductances and Transconductances as a function of the impurity potential B Imp at half filling. (The example is chosen such that the Transconductance in the unperturbed system is relatively large.) The conductance of the Heisenberg chain with one impurity, which is the large-U limit of g c , is given for comparison.
Although the two Cisconductances, g ↑↑ and g ↓↓ , could in principle differ (the model is now asymmetric) they do not in the case of half filling-at least within error bars. Both Cisconductance and Transconductance gomore or less linearly-to zero as the impurity strength increases.
We note that within error bars g c = −g t such that the full conductance of the system g = 2(g c + g t )  remains zero after insertion of the impurity. Furthermore, investigations with our method at different temperatures find no sizeable T -dependence.
In the case of zero chemical potential, µ = 0, we find a splitting of the two Cisconductances (see Fig. 5). This is particularly interesting since this finding implies a spinpolarized current. If we assume a (spin-independent) driving potential of the form V (P ↑ + P ↓ ) (V being the voltage amplitude), then the current is The average spin of a fermion in the current is therefore (using g ↓↑ = g ↑↓ = g t ) different from zero (see Fig. 6). An interesting problem is the question whether for µ = 0 the Cisconductance in the pure (spin-down) sector g ↓↓ survives or not when we increase B Imp to infinity. Note that a finite g ↓↓ would imply a total spin polarization, i.e., S = −1/2. The limit B Imp → ∞ of g ↓↓ cannot be taken directly (because of problems with the MC simulation), but here we note that there is another way to model the impurity. Instead of applying a local magnetic field on one site, one can introduce a weak link, i.e., decrease the hopping amplitude for spin-up fermions between the sites N/2 and N/2+1 from the initial value t to t ↑ Imp . These two variants of impurities behave similarly. 7 We computed the Cisconductance for the unaffected spin orientation for the model with t ↑ Imp = 0 (corresponds to B Imp = ∞) at µ = 0, U = 4t, N = 200 sites. We find a value of about [0.79 ± 0.03]e 2 /h.
The different behavior of the Cisconductance in the (unaffected) spin-down sector at half filling and away from half filling may be explained as follows: Suppose B Imp and U are large. The effect of the B Imp term on the fermions is that it forbids occupation of the impurity site for one of the two fermion species (in our case spin-up fermions). At half filling a spin-up fermion can hop only from one site to another by exchanging the site with a spin-down fermion (there are no empty sites), i.e., simultaneously with the spin-up a spin-down fermion must hop in the opposite direction et vice versa (implying g c = −g t ). A fermion of a certain spin index can then only pass the impurity site if accompanied by a fermion of opposite spin (which moves in the opposite direction). Since the impurity site is forbidden for one of the two fermion species, no fermion can pass the impurity site, and both Cisconductances must go to zero.
Away from half filling the hopping of an spin-up fermion does not necessarily require the hopping of a spin-down fermion (the spin-up fermion can hop to an empty site) and hence the impurity affects only one of the two Cisconductances.

C. Zig-zag chain
So far we have dealt with a system of two fermion species, where the two species reside on the same set of sites. In contrast to this, the bosonization approaches considered mostly systems of two coupled spinless-fermion conductors. We can compare our results with that situation, if we interpret the Hubbard model as a spinless-fermion ladder. Here one assumes that each fermion species lives on a different conductor [i.e., the two indices (n, ↑) and (n, ↓) are supposed to label (spatially) different sites; compare Sec. III and upper half of Fig.  1]. But one should note that for this case the parameter U should be small (since the distance between separated conductors is large) and the Coulomb interaction should be long-ranged (not on-site as in the Hubbard model). Hence we are led to the question how a variation in the interaction term modifies our results.
To address this question we add a new interaction term to the Hubbard model H = H Hubb + U n n n,↓ n n+1,↑ .
This Hamiltonian corresponds to Eq. (2) with U ′ = U . One can justify introducing this new interaction term if the fermion species live on different sites where each spindown site (n, ↓) lies between two spin-up sites, (n, ↑) and (n + 1, ↑). This model has therefore the geometry of a frustrated zig-zag chain as depicted in the lower half of Fig. 1.
One should note that this system has a total of 2N sites, N sites for each fermion species. [Although it would be useful to adopt the notion of a system with two coupled (spinless-fermion) chains, we will keep here the notation of a system of spinfull fermions.] The results for the spin drag in this model are shown in Fig. 7. We discuss again two chemical potentials: one is µ = U implying half filling, the other is again µ = 0. In the latter case the (mean) occupation number per site n ρ is different from one half (the occupation at half filling) and depends on U . It is shown in Fig. 8.
One sees in Fig. 7 that |g c,t | grow with the strength of the interaction. This may be explained as follows: First, the Coulomb interaction mediates an attractive nearestneighbor interaction for fermions with equal spin orientation (this is a consequence of the frustration). Therefore, in a simple approximation the only effect of the Coulomb interaction is to renormalize the Luttinger-liquid parameters for the two spin sectors K ↑,↓ . Since the Luttinger parameter for a spinless-fermion chain increases with the strength of the attraction, 8 we expect that K ↑,↓ increases as U increases. Since K ↑,↓ gives the conductance of one spin sector 17 (which is essentially the Cisconductance) we have that g c increases with U . One may also infer from the figure that the dependence of g c,t on a chemical-potential shift is weak. Within error bars g c decreases only slightly upon shifting µ away from half filling.
One should note that in the limit U = ∞ the ground state is a spin-polarized configuration (see Fig. 9). For µ = U this means that all conductances are zero in the large-U limit, for one spin sector is empty and the other, completely filled. In contrast to this, for µ = 0 one of the two spin sectors may remain conducting. The crossover to the ordered state occurs at values of U larger than 3t which may be seen by simulating and comparing the occupation number for different states. For µ = 0 the difference in occupation (between the two spin sectors) n σ = | n (n n,↑ − n n,↓ ) | /(2N ) is shown in Fig. 8; for µ = U it is zero within error bars as long as U ≤ 3t. We conclude that for the values of U considered in Fig. 7 the two spin sectors have approximately the same filling.

V. CONCLUSION
In this paper we discussed the spin drag for the Hubbard model at zero temperature. We found that the Transconductance is negative-at half filling the Umklapp even enforces g c = −g t . In that respect our situation is different from two coupled Tomonaga-Luttinger mod- If we assume that a given potential is in general not spin dependent, the only relevant quantity is the full conductance g = 2(g c + g t ), which is only nonzero away from half filling. Here both spin orientations contribute equally to the current. However, the situation changes when we add a magnetic impurity. Even if the driving potential is still spin independent, the resulting current will be (partially) spin polarized, if we are away from half filling. In the limit of a large impurity potential the current will be fully spin polarized. restriction of the original Hamiltonian to the assumed ground state configurations, i.e., zeroth order in U ) which is an xx chain in magnetic field, where the "spin" operator R z n (this time denoted by R to avoid confusion with previous spin operators) is defined by R z n = −1/2 if there is a particle on site n and R z n = +1/2 if site n is empty. Since the effective Hamiltonian describes the charge part of the Hamiltonian the full conductance 2(g c + g t ) for the original model should coincide with the conductance of the new Hamiltonian which is e 2 /h as the system is noninteracting. 8,17 Since we have g c = g t by the assumption of a spin gap and Eq. (4), the relation g = 2(g c + g t ) yields g c = 0.25 e 2 /h = g t .
In principle the model Eq. (A1) can also be analyzed with the Monte Carlo method developed in this paper, but we found that the simulation for this case is problematic: We measured large autocorrelation times for finite J z and µ (e.g., for the computation of the compressibility). We therefore must restrict ourselves to J z ≤ 0.8t.
For J z = 0.8t we present results for the Cis-and Transconductance in Fig. 10. In the large-U limit we find good agreement with our prediction that g c = g t = 0.25 e 2 /h which gives credit to the simulation data despite the large autocorrelation times.
Here we want to stress once again the remarkable fact that the sign of the Transconductance (the direction of the induced current) changes when we switch the magnetic field and the spin-polarized interaction on. (The Transconductance is for all U negative in Fig. 2 whereas in the present situation we expect g t = g c = K ρ /2 > 0 for T = 0, U = ∞.) Occupation in the ground state-Since we identified the ground state of the Hamiltonian Eq. (A1) H(µ = 0, U → ∞) with the ground state of the xx chain in magnetic field, we can calculate the occupation per state of this Hamiltonian in the large U -limit, the result being: This prediction may be tested against a Monte Carlo simulation. We find good agreement (see Fig. 11).