Spin subdiffusion in disordered Hubbard chain
Maciej Kozarzewski Peter Prelovsek Marcin Mierzejewski




Institute of Physics, University of Silesia, 40-007 Katowice, Poland

J. Stefan Institute, SI-1000 Ljubljana, Slovenia 
Faculty of Mathematics and Physics, University of Ljubljana, SI-1000 Ljubljana, Slovenia 

Department of Theoretical Physics, Faculty of Fundamental Problems of Technology, Wrocław University of Science and Technology, 50-370 Wrocław, Poland



We derive and study the effective spin model that explains the anomalous spin dynamics in the one-dimensional Hubbard model with strong potential disorder. 
Assuming that charges are localized, we show that spins are delocalized and their subdiffusive transport originates from a singular 
random distribution of spin exchange interactions.  The exponent relevant for the subdiffusion  is determined by the Anderson localization length 
and the density of electrons. While the analytical derivations are valid for low particle density, numerical results for the full model reveal a qualitative
agreement up to half–filling.


 

Spin subdiffusion in  disordered Hubbard chain
    Marcin Mierzejewski
    Received: 24 October 2018 / Accepted: TBD
==============================================



Introduction.–
The many–body localization  (MBL) <cit.> has recently  been intensively studied.  Vast amount of numerical data allowed to identify the main properties of the MBL systems:
 vanishing steady transport, <cit.>
 absence of thermalization <cit.>  and  logarithmic growth of the entanglement entropy <cit.>.
 It has also been found that MBL prevents a driven system from heating <cit.>.
 These unusual properties  can be  explained via  the existence of a macroscopic number of local integrals of motion <cit.>.

While most of theoretical studies so far concentrated on the one-dimensional (1D) disordered model of interacting spinless fermions,
the experiments on MBL are performed  on cold–fermion lattices <cit.> where the relevant model is
the Hubbard model with spin–1/2 fermions, whereby the disorder enters only via a random (or quasi-periodic) charge potential. 
Recent  numerical studies of such a  model  <cit.> 
reveal that even at strong disorder, localization and nonergodicity occurs only in the charge subsystem, implying a partial MBL. 
Unless one introduces also an additional random magnetic field <cit.>, the spin remain delocalized,
<cit.>, although the spin transport is anomalously slow and subdiffusive  
<cit.>.

In the present work we focus on the explanation of the slow spin dynamics and subdiffusion within the disordered 1D Hubbard model.
We first demonstrate that in the case of potential disorder and low particle density  the spin dynamics  can be described by a squeezed 
isotropic Heisenberg model, whereby the distribution of the random exchange interactions is singular. Such an effective model
can be studied  numerically quite in detail, but also analytically taking into account that the 1D spin dynamics is dominated by weak links.   




In this manner we show that spin excitations spread over distance M in a characteristic time t such  that M ∝ t^α with α≃λ /(d+λ), where 
d is the average distance between singly–occupied sites and λ is determined by the Anderson localization length  in the noninteracting system. 
While the mapping on the Heisenberg model is valid for dilute systems d ≫ 1, numerical results for strongly disordered Hubbard model
reveal that the same qualitative spin dynamics remains valid for all densities even up to half–filling. 

Model.–
Our aim is to establish the spin dynamics in  the disordered Hubbard chain,

    H   =    -t_ h∑_i σ   (c^†_i+1σ c_iσ+ H.c.)+∑_iε_i  n_i 
           + U ∑_i  n_i↑ n_i↓.
 
where n_i=n_i ↑ +n_i↓ and  n_i σ=c^†_iσ c_iσ. We study the model with L sites
and N electrons, fixing also total spin projection S^z_ tot=0. We assume  a uniform distribution of random charge potentials, ε_i ∈ [-W,W] and  set the hopping integral t_h=1.
 

Two electrons. In order to gain a preliminary insight to the spin dynamics, we first study two electrons.
The dynamics of a few interacting spinless particles has been studied  previously <cit.>. 
Here, we study for illustration N=2 electrons with opposite spin projections which  propagate on the chain with L=16 sites.
Assuming that particles are initially at sites j and l,  | ψ (0) ⟩ = c^†_j↓ c^†_j↑ |0⟩,   the propagation of  | ψ (t) ⟩ is obtained via exact diagonalization.  Figs. <ref>a and  <ref>b show, respectively,  time–dependence of the local spin ⟨ S^z_i⟩(t) =  1/2⟨ψ (t) | n_i ↑ - 
n_i↓    | ψ (t) ⟩  and   density  ⟨ n_i⟩(t) =   ⟨ψ (t) | n_i  | ψ (t) ⟩
 for one configuration of ε_i corresponding to  W=8.  While for such strong disorder,  the charge degrees appear to be fully
localized,  spins  undergo oscillations. 


Effective spin model.– The coexistence of almost frozen charges and oscillating spins suggests that one can derive an effective  spin model.
To this end, we use   the Anderson states as the basis, i.e., we use the single–particle eigenstates, ϕ_i a=⟨ i |a ⟩, of the noninteracting (U=0) model.
Then, 

    H   =   ∑_a σϵ_a c^†_aσ c_aσ +U/2∑_a a' b b' σχ^ab_a'b' c^†_aσ  c^†_b σ̅c_b' σ̅c_a' σ, 
    χ^ab_a'b'   =   ∑_i  ϕ^*_i aϕ^*_i bϕ_i b'ϕ_i a' .

Assuming that the charge-dynamics is frozen, the main effect arising from the presence of the Coulomb interaction comes from terms with either 
a=a', b=b' or a=b', b=a',  and in both cases a  b.  Then, the Hubbard  term in Eq. (<ref>) can be written in a SU(2)–invariant form

    H_U=1/2∑_a  b   J_ab( 1/4 n_a n_b - S⃗_a ·S⃗_b ),

where we use standard density and spin operators:  n_a=n_a ↑+n_a ↓, 
S^z_a=1/2 (n_a ↑-n_a ↓), S^+_a= c^†_a↑ c_a ↓, S^-_a=c^†_a ↓ c_a↑.
The effective exchange interaction is ferromagnetic

    J_ab=2 U χ^ab_ab=2 U ∑_i   |ϕ_ia |^2  |ϕ_ib |^2.


In order to test the approximation Eq. (<ref>), we consider  N=2  electrons with opposite spins, located 
at sites j and j+d. We evaluate the spin–oscillation frequency, ω_2, directly from results for the Hubbard model, see Fig. <ref>a. 
For the same set of ε_i we then identify Anderson states a, b which maximize |ϕ_ja|^2  and  |ϕ_j+d b|^2,  respectively.
This enables the calculation of J_ab from Eq. (<ref>) that should lead to  spin oscillations ⟨ S^z_j,l⟩ = ±1/2cos(J_ab t ).
Figs. <ref>c and  <ref>d  show correlations between ω_2 and J_ab for various realizations of disorder and various distances d between the electrons.
One finds that indeed ω_2 ≃ J_ab for strong disorder W ≫ 1 and d ≫ 1.

For low density of carriers and larger disorder,  the maxima of the occupied Anderson states a and b are typically 
separated  by  x_ab > ξ, exceeding the single-particle localization length, ξ. Then, one  obtains an approximate relation

    J_ab≃ 2U exp(-x_ab/λ),     λ∼ξ


The squeezed spin model.– Assuming that  charges are frozen to the initial occupations n_i= 0,1,2 it is evident that the effective spin model 
Eq. (<ref>) acts only on singly occupied sites with n_i=1. Spin dynamics of the Hubbard model at high temperatures T ≫ W,U can be then studied 
by first  randomly positioning  N electrons on L sites. This allows us to establish the distribution of distances between the  singly occupied sites as well as the distribution 
of the effective J_ab, using Eq. (<ref>) or its simplified version, Eq. (<ref>).  Due to the exponential decay of J_ab  we  consider only 
couplings between the neighboring singly occupied sites.  The effective Heisenberg model on a squeezed chain then reads

    H_H=-∑_i  J_iS⃗_i ·S⃗_i+1,
 
where the summation is carried out over  singly occupied sites n_i=1, i ∈{1,...,Ñ} with Ñ≤ N.  
Note that  at infinite temperature, the  average lattice–spacing in the effective model equals d= L/ Ñ =(n̅-n̅^2/2)^-1, where n̅=N/L is the average filling in the original Hubbard model.
  
 In order to establish the probability distribution of J_i, we first consider a section of length L ≫ 1, 
where we randomly choose the continuous positions of Ñ=L/d points and study the regime L ≫ 1.  The probability density 
for the distances between the neighboring points is 
f_d(x)=1/dexp (-x/d).  While the latter result has formally been obtained for continuous positions of points,  it should hold true also for discrete positions of singly occupied sites  provided that d ≫ 1. Using this result one may find the probability density for the random exchange interaction f_J(J). 
To this end, we use Eq. (<ref>) and  compare the cumulative distribution functions

    ∫_0^y dx  1/dexp(- x/d) =∫_2U exp(-y/λ)^2Ud J  f_J(J).

Taking the derivative of  Eq. (<ref>) with respect to y and introducing  the dimensionless interaction J̃=J/2U one gets

    f_J̃ (J̃)=λ̃J̃^λ̃-1,     λ̃= λ/d.

It is clear  that the interaction U  sets the energy scale (and the time scale)  of the effective model, 
whereas the qualitative spin dynamics depends on the ratio between effective localization length, λ, and inter–particle distance d. 
The important message is that for low doping (d≫ 1) and strong disorder (λ∼ 1)  one obtains λ̃≪ 1 with
 the distribution of J̃ being singular at J̃=0. Still, lim_δ→ 0^+∫_0^δdJ̃  f_J̃ (J̃)=0, 
hence the probability for cutting the Heisenberg chain into disconnected subchains is vanishingly small. 

 

 

In order to test feasibility and accuracy of Eq. (<ref>) we have numerically generated the distribution of 
J̃= J_ab/2U also  directly from  Eq. (<ref>),  in the same way as discussed for N=2 case.
The positions of two electrons l  and j have been randomly chosen in such a way that the distance x=| l-j | 
is distributed according to  f_d(x) for various d.  Numerical results for W=4 and 8 are shown in 
Figs. <ref>a and  <ref>b, respectively. These results have been fitted by Eq.(<ref>),  whereby we adjusted a single parameter
λ for all values of d.  We have obtained λ≃ 0.75 and  λ≃ 0.4  for W=4 and W=8, respectively. 
Although Eq. (<ref>) has been derived for d ≫ 1,  it turns out to remain qualitatively valid also for d=2, i.e. for the
average distance between singly–occupied sites in the half–filled Hubbard model.
We conclude that Eq. (<ref>) accurately describes
f_J̃ (J̃)  at least for small J̃, i.e.,  in the regime which is essential for the long–time spin dynamics.  

Local spin correlations.– 

We first calculate the time–dependent local spin correlations at infinite temperature,

    S_L(t)= ⟨ S^z_i S^z_i(t) ⟩ = 1/ Tr  1⟨ Tr [  S^z_i (t) S^z_i ] ⟩_ dis
 
where the spin evolution is determined  by the effective H_H. We take the random interaction J_i=J̃ as given by Eq.  (<ref>), i.e.,  
we express time in units of  1/2U. ⟨ ... ⟩_ dis means averaging over various realizations of J_i and  we use at least 2000 disorder samples.
  
Figs.  <ref>c and  <ref>d show S_L(t).   For longer times and  λ̃ < 1 one observes
power-law decay S_L(t) ∝ (2Ut)^-α with α  < 1/2, hence the spin dynamics is clearly subdiffusive.  
In Fig. <ref>a we demonstrate α obtained from fitting numerical results  by S_L(t) ∝ t^-α in the  time–window t ∈ [10,50]. 
The main message coming from these studies is that α >0 for arbitrary nonzero λ̃ > 0, i.e., for arbitrary nonzero filling.
For λ̃≪ 1 we obtain the exponent α≃λ̃. Still, it should be noted that the distribution Eq. (<ref>) is singular only for λ̃<1 which should be the regime of subdiffusion. For λ̃≪ 1 the finite size effects are negligible (Fig. <ref>c) but become significant for larger λ̃, Fig. <ref>d. 
Nevertheless, for the regime with λ̃≃ 1 (relevant for larger filling n̅∼ 1 and/or weaker disorder)  our results shown in Fig. <ref>d are consistent  with normal spin diffusion with α =1/2, which is also expected in the weakly disordered Hubbard model.  Results in  Figs. <ref>c, <ref>d and <ref>a support the scenario, that the spin excitations spread subdiffusively  due to the singular distribution of random exchange interactions, Eq. (<ref>).

Single weak–link scenario.–
To explain the relation of the subdiffusive  dynamics and the singular distribution of J̃_i,  we consider a single spin excitation and 
estimate the time, t,  in which the excitation spreads over M sites in the effective Heisenberg chain. We assume the weak–link scenario, where the long–time dynamics
is govern by rare regions with the smallest J_i. Similar approach has been used to describe the subdiffusive transport of spinless particles in the vicinity of the MBL transition <cit.>. 
Here, we assume that t ∼ 1/(2U J̃_m), where J̃_m is the weakest exchange out of J̃_i for i=1,...,M.

The probability that each random J̃_i is larger than J̃_0 is

    [ ∫_J̃_0^1 dJ̃  f_J̃ (J̃) ]^M  
    
    =∫_J̃_0^1 dJ̃_m  f_m(J̃_m),
 
where  f_m(J̃_m) is the probability density for  the smallest interaction.



 Using Eq. (<ref>) and calculating derivative of Eq. (<ref>) with respect to  J̃_0 one finds  the distribution function
 
 f_m(J̃_0)= λ̃ M  J̃_0^λ̃-1  (1-J̃^λ̃_0 )^M-1.
 
Then, the expectation value of the smallest exchange interaction out of M  random J̃_i reads
 
    ⟨J̃_m⟩   =    ∫_0^1 dJ̃_m  f_m(J̃_m)  J̃_m 
    
     
      ≃Γ(1+ 1/λ̃) M^-1/λ̃.
 
In the latter equation we have used  

formulas for asymptotics  at M ≫ 1/λ̃. 


So we find the relation between the spread of the spin excitations Λ and t as,

    Λ∼ M d ∝( 2U  t   )^λ̃,      S_L(t) ∝Λ^-1∝( 2U  t  )^-λ̃.


 The exponent α=λ̃  is the same as previously obtained from numerical studies of the effective Heisenberg model for λ̃≪ 1. 
 

Multiple weak–link scenario.–  The single–weak link scenario breaks down for λ̃∼ 1, where α≃λ̃/2 instead of  λ̃, as shown 
in  Fig. <ref>d.
As an alternative explanation for the subdiffusive transport we consider 
distribution of effective hopping times between neighboring sites in the squeezed spin model. The relevant dimensionless quantity is  τ=1/J̃. Using Eq. (<ref>) one finds
the probability density f_τ(τ)=λ̃/τ^λ̃+1.  For such broad distribution of hopping times in a classical model of random traps <cit.>, 
 one gets subdiffusive transport Λ∝( 2U  t   )^α where α=λ̃/(1+λ̃) for λ̃  < 1.  In the latter model, a classical particle
may hop between neighboring traps in time τ. The hopping time, τ, is randomly chosen for each site but remains the same for each visit of the same site. This simple model quite accurately reproduces the dynamical exponent  α for arbitrary λ̃, as shown in Fig. <ref>a, whereas for λ̃≪ 1 it gives the same relation as
the single weak–link scenario.

Comparison with the Hubbard model.– 
Finally,  we compare our analytical predictions with numerical results obtained directly for the Hubbard model. 

The time-dependent local spin correlation function S_L(t) at infinite temperature has been obtained using the microcanonical Lanczos method (MCLM) <cit.> (in analogy to the imbalance correlations presented previously <cit.>) for the Hubbard model with L=18 and N=6, i.e. d ≃ 3.5.
Results are shown in Figs. <ref>b,  <ref>c  together with analytical prediction, S_L(t) ∝ (2Ut)^-α with α=λ/(λ+d) and λ obtained from fits in Fig. <ref>b.  
Despite significant finite–size effects, one observes  that the latter estimate  correctly describes the subdiffusive spin dynamics in the Hubbard model at low filling n̅≪ 1. 
In particular, the exponent α obtained directly from  the Hubbard model weakly depends on U.   

Moreover, one can consider the validity of the subdiffusion scenario in the full Hubbard model beyond the limit of low filling.
It has previously been  found  that spins reveal a subdiffusive dynamics even for n̅ =1 <cit.>. We therefore analyze the MCLM  results for 
S_L(t) ∝  t^-α  considering various system sizes L=14,16,18 and various numbers of electrons N. 
For the time–window t ∈ [1,10]  we extract  the dynamical exponent α  and compare with  α=λ/(λ+d), as shown in Fig. <ref>d.

Our approach works even up to  n̅=1 since the average distance between singly occupied sites  d ≥  2,
while λ<1 provided that the disorder is sufficiently strong. In particular, for W=8  and   n̅≤ 1 we have estimated  that λ̃≲ 0.2.


 Conclusions.– In this paper we presented an explanation for the anomalous spin dynamics  in 1D Hubbard model with large potential disorder
 in the regime of partial MBL, where the charge dynamics appears to be frozen whereas  spins exhibit ergodic but subdiffusive transport <cit.>.     
 We have derived an effective isotropic Heisenberg model
with random exchange interactions between neighboring singly–occupied  sites J_i.  Our derivation is formally best applicable to the regime of low filling,
n̅≪ 1,   and strong disorder. 
The main origin of the subdiffusive behavior then appears to be the singular distribution of the effective exchange interaction, J_i, with the crucial parameter  λ̃= λ/d  
representing the ratio of the single-particle Anderson localization length
and the  average distance between singly occupied sites.  Results for the Hubbard model reveal that such scenario seems to remain qualitatively  
valid beyond the considered limits  of low filling,  even at n̅=1, provided that the disorder is sufficiently strong.  
It appears that there is no  threshold filling n̅_c, below which also spins  would become localized and full MBL would prevail.

There are still questions concerning the  dynamics within disordered Hubbard model, being relevant also to cold-atom experiments on MBL <cit.>.

Our derivation of the effective model, remains on the level of spin dynamics, while  charge degrees are assumed to be frozen. It is evident that higher order terms in the Anderson basis, following from Eq. (<ref>), would lead also to the dynamical coupling between charge and spin degrees of freedom. Since the spin dynamics is ergodic, it is not excluded that also charges would eventually delocalize, but then on much larger time scales.









 
 
  
This work is supported by the National Science Centre, Poland via project 2016/23/B/ST3/00647 (MM and MK) and P.P. acknowledges the support by the program P1-0044 of the
Slovenian Research Agency.






