Density dynamics in the mass-imbalanced Hubbard chain

We consider two mutually interacting fermionic particle species on a one-dimensional lattice and study how the mass ratio $\eta$ between the two species affects the (equilibration) dynamics of the particles. Focussing on the regime of strong interactions and high temperatures, two well-studied points of reference are given by (i) the case of equal masses ${\eta = 1}$, i.e., the standard Fermi-Hubbard chain, where initial non-equilibrium density distributions are known to decay, and (ii) the case of one particle species being infinitely heavy, ${\eta = 0}$, leading to a localization of the lighter particles in an effective disorder potential. Given these two opposing cases, the dynamics in the case of intermediate mass ratios ${0<\eta<1}$ is of particular interest. To this end, we study the real-time dynamics of pure states featuring a sharp initial non-equilibrium density profile. Relying on the concept of dynamical quantum typicality, the resulting non-equilibrium dynamics can be related to equilibrium correlation functions. Summarizing our main results, we observe that diffusive transport occurs for moderate values of the mass imbalance, and manifests itself in a Gaussian spreading of real-space density profiles and an exponential decay of density modes in momentum space. For stronger imbalances, we provide evidence that transport becomes anomalous on intermediate time scales and, in particular, our results are consistent with the absence of strict localization in the long-time limit for any ${\eta>0}$. Based on our numerical analysis, we provide an estimate for the ''lifetime'' of the effective localization as a function of $\eta$.


I. INTRODUCTION
Understanding the dynamics of quantum many-body systems is a central objective of modern physics which has been reignited by experimental advancements featuring, e.g., cold atoms or trapped ions [1,2], and has experienced an upsurge of interest also from the theoretical side [3][4][5][6][7]. In this context, an intriguing and fundamental direction of research is to explain if and how thermodynamic behavior can emerge from the unitary time evolution of isolated quantum systems. One notable explanation for this occurrence of thermalization is the eigenstate thermalization hypothesis (ETH) [8][9][10], which has been numerically verified in numerous instances [5].
However, despite thermalization certainly being a rather common observation, there are also classes of systems which generically evade to reach thermal equilibrium even at indefinitely long times. In particular, it has been realized early on by Anderson that non-interacting particles in one or two spatial dimensions localize for an arbitrarily weak disorder potential [11,12] (for experimental confirmations see, e.g., [13,14]). Moreover, it is now widely believed that for sufficiently strong disorder, localization is also possible in the presence of interactions [15,16], which is supported by experimental results as well [17].
While the majority of studies on many-body localization (MBL) typically focus on one-dimensional and shortranged models composed of, e.g., spin-1/2 degrees of freedom, there has been much effort recently to generalize the * tjark.heitmann@uos.de † rsteinig@uos.de notion of MBL to a wider class of models [18]. This includes, e.g., systems which are weakly coupled to a thermal bath [19], models with long-range interactions [20] or degrees of freedom with higher spin S > 1/2 [21,22], as well as Hubbard models where the disorder only couples to either one of the charge or spin degrees of freedom [23,24].
A particularly interesting question is whether MBL can also occur in systems which are translational invariant, i.e., without any explicit disorder [25][26][27][28][29][30][31][32][33][34][35][36][37][38][39]. A convenient model to investigate this question is given by the mass-imbalanced Hubbard chain [29-31, 39, 40]. In this model, two mutually interacting particle species are defined on a one-dimensional lattice and exhibit different hopping amplitudes. Here, the imbalance is parametrized by the ratio η between the two hopping strengths, ranging from η = 0, where the heavy particles are entirely static, to η = 1, where the hopping amplitudes are the same. On the one hand, in the balanced limit η = 1, numerical evidence for diffusive [41][42][43] (or superdiffusive [41,44]) charge transport has been found in the regime of high temperatures and strong interactions. On the other hand, for η = 0, the static particle species creates an effective disorder potential which induces localization of the lighter particles [40,[45][46][47][48]. In view of these two opposing cases, it is intriguing to study the dynamics in the regime of intermediate imbalances 0 < η < 1. While genuine localization (i.e. on indefinite time scales) is most likely absent for any η > 0 [35,39], e.g., due to slow anomalous diffusion which ultimately leads to thermalization [35], this does not exclude the possibility of interesting dynamical properties such as a "quasi-MBL phase" at short to intermediate times [35].
In this paper, we scrutinize the impact of a finite arXiv:2004.13604v1 [cond-mat.str-el] 28 Apr 2020 Figure 1. (Color online) Illustration of the imbalanced Fermi-Hubbard chain. Spin-↑ and -↓ particles with on-site interaction of strength U and different hopping amplitudes t ↑ and t ↓ . Diffusive broadening of the initially peaked spin-↑ density profile is sketched as a possible scenario depending on the imbalance ratio η = t ↓ /t ↑ .
mass imbalance 0 < η < 1 from a different perspective by studying the real-time dynamics of pure states featuring a sharp initial non-equilibrium density profile. Relying on the concept of dynamical quantum typicality, the resulting non-equilibrium dynamics can be related to equilibrium correlation functions. Summarizing our main results, we observe that diffusive transport occurs for moderate values of the mass imbalance, and that it manifests itself in a Gaussian spreading of real-space density profiles and an exponential decay of density modes in momentum space. Moreover, for stronger imbalances, we find evidence that on the time and length scales numerically accessible, transport properties become anomalous, albeit we cannot rule out that normal diffusion eventually prevails at even longer times. Furthermore, our results are consistent with the absence of genuine localization for any η > 0. In particular, we find that for smaller and smaller values of η > 0, the resulting dynamics resembles the localized η = 0 case for longer and longer time scales. However, we conjecture that this "lifetime" of effective localization always remains finite for a finite η. This paper is structured as follows. After introducing the model in Sec. II, we give an introduction to the employed typicality approach, and the initial states and observables in Sec. III. We then present our results in Sec. IV and conclude with a discussion in Sec. V.

II. MODEL
We study the Hubbard chain describing interacting spin-↑ and -↓ fermions on a one-dimensional lattice. The Hamiltonian for L lattice sites with periodic boundary conditions (L + 1 ≡ 1) reads with local terms where the creation (annihilation) operator c † r,σ (c r,σ ) creates (annihilates) a fermion with spin σ at site r, and n r,σ = c † r,σ c r,σ is the particle-number operator. The first term on the right hand side of Eq. (2) describes the siteto-site hopping of each particle species with amplitude t σ . The second term is the on-site interaction between the particle species with strength U , see Fig. 1. The imbalance between t ↑ and t ↓ ≤ t ↑ is parametrized by the ratio ranging from η = 0 for t ↓ = 0 to η = 1 in the case of t ↓ = t ↑ .
While the Hamiltonian H in Eqs. (1) and (2) is integrable in terms of the Bethe Ansatz for η = 1 (i.e. in the case of the standard Fermi-Hubbard chain, see, e.g., Ref. [49]), this integrability is broken for any finite imbalance 0 < η < 1. Moreover, despite its integrability, there has been numerical evidence that, in the regime of high temperatures and strong interactions, charge transport in the one-dimensional Fermi-Hubbard is diffusive [41][42][43] (or superdiffusive [41,44]). In order to have this well-controlled point of reference for our analysis of finite imbalances η ≤ 1, we here fix the interaction strength to the large value U/t ↑ = 16.
Using this symmetry, the Hamiltonian (1) can be decoupled into 2 L independent subspaces, effectively describing non-interacting spin-↑ particles on a onedimensional lattice with random (binary) on-site potentials ±1/2 U (n r,↑ − 1/2), which implies the onset of Anderson localization [11]. It is worth mentioning that by means of a Jordan-Wigner transformation, the fermionic model in Eqs. (1) and (2) can be mapped to a spin-1/2 model with ladder geometry [35]. This spin model is described by the local terms + J ⊥ s z r,1 s z r,2 with J 1 = t ↑ , J 2 = t ↓ and J ⊥ = U . Here, the different particle species ↑ and ↓ are represented as local magnetizations on the two separate legs k = 1, 2 of the ladder. The hopping term and the interaction term in the Hubbard formulation correspond to the XY interaction along the legs and the Ising interaction on the rungs of the ladder, respectively. The particle number conservation N σ = r n r,σ = const. for both particle types translates into magnetization conservation M k = r s z r,k = const. on each leg.

A. Initial states and observables
We investigate the real-time dynamics of local particle densities given by the expectation values with the density matrix for pure initial states |ψ(0) , such that In order to realize inhomogeneous particle densities, we prepare the initial states via the projection The reference pure state |φ is constructed as a random superposition, where the |ϕ k denote the common eigenbasis of the local occupation number operators n r,σ , and the sum runs over the full Hilbert space with finite dimension d = 4 L . (In spin language, this simply is the Ising basis.) Moreover, the complex coefficients c k are randomly drawn from a distribution which is invariant under all unitary transformations in the Hilbert space (Haar measure) [52,53], i.e., real and imaginary parts of these coefficients are normally distributed with zero mean. As a consequence, the initial density profile exhibits a sharp delta peak for the spin-↑ particles in the middle of the chain on top of a homogeneous many-particle background [43,54], Rather than taking the full Hilbert space into account, one could also consider the half-filling sector (respectively the zero-magnetization sector).

B. Dynamical quantum typicality
Given the specific construction of the pure state |φ in Eq. (10), the concept of dynamical quantum typicality (DQT) provides a direct connection between the nonequilibrium expectation value p r,↑ (t) and an equilibrium correlation function (see Ref. [43] and also Appendix A), where the thermodynamic average • = Tr [ • ] /d is carried out at formally infinite temperature. As a consequence, the dynamics of the non-equilibrium expectation value p r,↑ (t) can be used to study transport properties within the framework of linear response theory. Importantly, the variance of the statistical error = (|φ ) of Eq. (12) is bounded from above by i.e., the accuracy of the typicality approximation improves exponentially upon increasing the size of the system. In principle, this error can be further reduced by averaging over multiple realizations of the random state |φ [55,56]. However, for the system sizes studied here, the DQT approach is already very accurate, and this additional sampling becomes unnecessary. More details on the concept of dynamical quantum typicality (and on error bounds) can be found in Refs. [57][58][59][60][61][62][63][64][65][66][67][68][69][70][71][72].

C. Time evolution via pure-state propagation
For the time evolution of the pure state we can bypass the exact diagonalization of the Hamiltonian and rather solve the time-dependent Schrödinger equation directly via iterative forward propagation in small time steps δt. Aside from the many numerical methods such as Trotter decompositions [73,74], Chebyshev polynomials [75][76][77] or Krylov-space methods [78], the action of the time-evolution operator in each step can be calculated by a fourth-order Runge-Kutta scheme [66,67], Crucially, the matrix-vector multiplications in Eq. (15) can be implemented very memory efficiently due to the sparse matrix representation of the given Hamiltonian. While the action of H on |ψ can also be calculated onthe-fly, we save the sparse Hamiltonian matrix for the sake of run time. Moreover, symmetries of the system can be exploited in order to split the problem into smaller sub-problems and to further reduce the computational effort [79]. In this paper, we exploit the particle number (magnetization) conservation for both particle species (legs) separately. As a consequence, the maximum memory consumption for the largest symmetry sector in a system of length L = 15, with full Hilbert-space dimension d ∼ 10 9 , amounts to about 20 GB. While L = 14 or L = 15 are already comparatively large (especially in view of the extensive parameter screening and the long simulation times considered here), let us note that even larger system sizes can be treated by the usage of largescale supercomputing (see, e.g., Refs. [40,43] with the diffusion constant D. For the δ-peak initial conditions (11), the solution of Eq. (16) can be well approximated by the Gaussian function where the spatial variance scales as Σ 2 (t) = 2Dt and is given by with δp r,↑ (t) ∝ p r,↑ (t) − p eq. fulfilling r δp r,↑ (t) = 1 for all times t. More generally, a scaling of the variance according to Σ(t) ∝ t α is called ballistic for α = 1, superdiffusive for 1/2 < α < 1, diffusive for α = 1/2, subdiffusive for 0 < α < 1/2, and localized for α = 0. Moreover, away from the case α = 1/2, the density profiles p r,↑ (t) are not expected to take on a Gaussian shape.

Connection to current-current correlation functions
Due to the typicality relation (12), the spatial variance in Eq. (18) can be related to the dynamics of currentcurrent correlation functions via [81] where the time-dependent diffusion constant D(t) is given by and j ↑ denotes the total current operator of the spin-↑ particles, (Note that the relation (19) requires δp r,↑ (t) to vanish at the boundaries of the chain [81].) We therefore can compare the spatial variance of density profiles calculated according to Eq. (18) to the one already obtained from current-current correlation functions [42,81,82], A detailed analysis of transport in the mass-imbalanced Hubbard chain extracted from current-current correlation functions can be found in [40].

Momentum space
In addition to the real-space perspective, it is also instructive to look at momentum-space observables as given by the lattice Fourier transform of the density profiles, with the momentum q = 2πk/L and wave numbers k = 0, 1, . . . , L − 1. In particular, the Fourier transformation of the diffusion equation (16) yields the corresponding diffusion equation for the p q,↑ (t), with q 2 = 2(1 − cos q). From Eq. (24), it becomes clear that diffusion manifests itself in momentum space by exponentially decaying modes

IV. RESULTS
We now turn to our numerical results. To begin with, the two limiting cases η = 1 and η = 0 are presented in Sec. IV A. Intermediate imbalances 0 < η < 1 are discussed in Secs. IV B and IV C.  The initial density peak in the center of the chain spreads rather quickly over the system for η = 1, whereas it appears to be frozen for η = 0.

A. Limiting cases
In order to mark out the two completely different behaviors of the density dynamics in the limiting cases of the model, we first discuss the limit of equal particle masses (η = 1) and contrast it with the limit of infinite mass-imbalance (η = 0). Recall that the interaction strength is set to U/t ↑ = 16 in the following.
First, Fig. 2 shows the real-time broadening of the initially peaked density profiles p r,↑ (t) for both limits in a time-space density plot. While the particle density for   17) and (18). In the case η = 0, an overall triangular shape survives even for long times, with some local fluctuations.   For a more detailed analysis, the spatial dependence of the profiles p r,↑ (t) − p eq. is shown in Fig. 3 for fixed times t in a semi-logarithmic plot. Remarkably, the profiles for η = 1 in Fig. 3 (a) can be very well described by Gaussians [see Eq. (17)] over three orders of magnitude. These Gaussian profiles indicate that charge transport in the integrable Fermi-Hubbard chain is diffusive [41][42][43], at least in this parameter regime (strong interactions and high temperatures) and for the time scales depicted, see also Refs. [41,44] for the possibility of superdiffusive transport. Note that the Gaussians in Fig. 3 (a) are no fit, since the width Σ(t) has been calculated exactly according to Eq. (18), i.e., there is no free parameter involved. In contrast, the profiles for η = 0 in Fig. 3 (b) are clearly non-Gaussian and remain, even for the long times shown, in an overall triangular shape with variations on short length scales.

Real space
Next, let us study a finite imbalance between the particle masses. In analogy to Fig. 2, time-space density plots are shown in Fig. 4 for η = 0.8 and η = 0.6. For these ratios the broadening of the initial density peak apparently happens on a time scale comparable to the one observed for η = 1 in Fig. 2, with a barely noticeable slowdown with the increasing imbalance. Similar observations can be made for the density profiles at fixed times, as shown in Fig. 5. At weak imbalance η = 0.8 [ Fig. 5 (a) 17) and (18). At moderate imbalance (a) η = 0.8, the profiles can be very well described by Gaussians. At slightly smaller (b) η = 0.6, the density profiles are still in good agreement with Gaussians, although little deviations become apparent at t t ↑ = 4. η = 1. Even for stronger imbalance η = 0.6 [ Fig. 5 (b)], the profiles appear to be of Gaussian shape, although little deviations start to appear at t t ↑ = 4, which might be seen as the onset of a drift from normal to anomalous diffusion, see also the discussion below.

Spatial width
In order to analyze the broadening of the density profiles further, Fig. 6 shows the time-dependence of the spatial width Σ(t) obtained by Eq. (18) for moderate imbalance η = 0.8 and different system sizes L = 10, . . . , 14. Necessarily, there is an initial linear increase Σ(t) ∝ t for t t ↑ 1, indicating ballistic transport, as it is expected for short times below the mean-free time. Subsequently,   Σ(t) shows a scaling ∝ √ t, consistent with diffusion. However, for later times, we find that Σ(t) approaches a saturation value which increases with increasing L. This behavior of Σ(t) can be easily understood since the width of a density profile on a finite lattice with L sites is obviously bounded from above. Namely, assuming equilibration, i.e., a perfectly homogeneous distribution of the p r,↑ for t → ∞ with δp r,↑ = 1/L at each site, we obtain the saturation value This L-dependent saturation value is reached quickly for the weakly imbalanced case η = 0.8 in Fig. 6, e.g., Σ ≈ 4 for L = 14. Moreover, for the biggest size L = 14, Fig. 6 also shows Σ(t) calculated from current-current correlation functions via Eq. (22). Overall, the behavior of this Σ(t) is in good agreement with the one described above. Note that the small deviations between the two widths setting in at t t ↑ ∼ 6 presumably arise when the tails of the density distribution reach the boundaries of the system (cf. Fig. 5). Additionally, we note that the finite-size saturation value (26) does not apply to Eq. (22), which, by definition, is not bounded. Rather, for times t t ↑ 6, we find an accelerated increase of Σ(t). This is caused by the fact that the current-current correlation function j ↑ (t)j ↑ does not completely decay to zero in a system of finite size, see also Refs. [40,80].

Momentum space
Complementary to the real-space data for η = 0.8, 0.6 shown in Fig. 5, the corresponding Fourier modes p q,↑ (t) with momentum q = 2πk/L are shown in Fig. 7 for the four longest wavelengths available, i.e., k = 1, . . . , 4. While p q,↑ (t) decays rather quickly for k ≥ 2 (with the decay rate increasing with k), we find that at least for k = 1, p q,↑ (t) is to good quality described by an exponential decay [see Eq. (25)], consistent with the onset of diffusion on the corresponding length scales.

C. Strong mass-imbalance
Now, let us study how the equilibration dynamics alter for stronger imbalances and also discuss the possibility of localization for η > 0.

Real-space dynamics
Before discussing the full density profile in detail, let us for simplicity focus on the decay of the central peak p L/2,↑ (t), as shown in Fig. 8 for imbalance ratios η = 0, . . . , 1 and two system sizes L = 12 and L = 14. While p L/2,↑ (t) ∝ t −1/2 to good quality for η = 1, consistent with diffusive transport, this decay is slowed down with decreasing η. At small but finite η = 0.1, we find that p L/2,↑ (t) approximately coincides with the η = 0 curve up to times t t ↑ ≈ 40, until it eventually starts to decay towards the equilibrium value p L/2,↑ (t → ∞) − p eq. = (1 − p eq. )/L. Note that the two curves for L = 12, 14 agree very well with each other before the equilibration value is reached. On these time scales, the behavior of the density dynamics thus appears to be independent of the system size. For additional data with smaller η and longer time scales, see Appendix B. Moreover, a more detailed finite-size analysis can be found in Appendix C.
Next, let us come back to a discussion of the full density profile. To this end, Fig. 9 shows time-space density   plots for the two η = 0.4 and η = 0.2. We find that the broadening of the density profiles visibly slows down with decreasing η, until no substantial spreading of the density can be observed for η = 0.2 up to the maximum time t t ↑ = 40 shown here, consistent with Fig. 8 discussed before.
The corresponding cuts of the density profiles at fixed times are shown in Figs. 10 (a) and (b). Note that, owing to the slow broadening of the profiles, we show cuts at later times compared to Fig. 5. One clearly sees that the profiles are not Gaussian anymore, but rather exhibit a pronounced triangular shape in the semi-logarithmic plot used. In particular, they can be well described by the function  with the time-dependent fit parameters α(t), Σ f (t), and β(t). In particular, the exponent α(t) ∈ [1, 2] is introduced to capture the triangular shape. This shape indicates a crossover to anomalous diffusion for small ratios η 0.4 [83]. This is another central result of this paper.

Spatial width
Additionally, Fig. 11 shows the width Σ(t) of the density profiles for η = 0.4 and η = 0.2, as calculated by Eqs. (18) and (22). Compared to the weakly imbalanced case shown in Fig. 6, Σ(t) now grows much slower, and Eqs. (18) and (22) are in better agreement, since the distribution is still well concentrated in the center of the chain. For η = 0.2, Σ(t) appears to remain at a constant plateau up to the maximum time t t ↑ = 20 shown.
To analyze the η-dependence of the width in more detail, Fig. 12 shows Σ(t) in Eq. (18) on a longer time scale t t ↑ ≤ 150 for various values of η and a fixed system size L = 14. While the growth of Σ(t) towards the saturation value becomes slower and slower with decreasing η, we find that even for the smallest value of η = 0.05 shown here, Σ(t) clearly increases at long times. In contrast, the width in the η = 0 case fluctuates around a constant and lower value, which might be interpreted as the Anderson localization length.

Momentum-space dynamics
Let us now turn to momentum-space dynamics again. To this end, Figs. 13 (a) and (b) show the discrete Fourier modes p q,↑ (t) for imbalance ratios η = 0.2 and η = 0.1. Note that the data is obtained for an even larger system with L = 15 lattice sites and for momenta q = 2πk/L with k = 1, 2, 3. Compared to Fig. 7, we find that the p q,↑ (t) now decay visibly slower for all wave numbers k. Moreover, in contrast to the scaling of decay rates in the case of normal diffusion [cf. Eq. (25)], the density modes now seem to decay at a similar rate for all k. Furthermore, even for small η = 0.1 and k = 1, we find that p q,↑ (t) is clearly non-constant, which suggests that genuine localization is absent for η > 0.
To analyze the dependence on system size, Fig. 13 also shows the Fourier mode p q,↑ (t) for L = 10 and wave number k = 2. This mode has the same momentum q = 2π/5 as the mode k = 3 for L = 15. We find that for both η = 0.2 and η = 0.1 the decay of significant differences up to the maximum time t t ↑ = 120 shown. Finally, Fig. 14 (a) shows the relaxation of the Fourier mode p q,↑ (t) with the smallest wave number k = 1 for various 0 ≤ η ≤ 1. The decay appears to be exponential for η 0.2, albeit very slow for strong imbalances. While for sufficiently small times, all η > 0 curves agree with the η = 0 curve, they start to deviate at a certain point in time. In order to analyze this separation time from the η = 0 curve in more detail, we define using the running averages of the density modes It measures the maximum time up to which η = 0 and η > 0 curves do not deviate up to a distance δ. (Note that this maximum time can not exceed the maximum simulation time, here t max t ↑ = 1000. Moreover, the running averages are used to mitigate the fluctuations of the p q,↑ (t), which complicate the extraction of precise separation times.) The physical picture for this analysis can be understood as follows. For very small but nonzero η, the heavy particles still appear as a quasi-static disorder potential for the lighter particles, which induces localization analogous to η = 0. At some point in time, however, the residual hopping of the heavy particles becomes relevant, which can be seen as an η-dependent "lifetime" of the Anderson insulator. The corresponding data for different distances δ is shown in Fig. 14 (b). For every δ, the lifetime grows fast with decreasing η, but apparently is always finite for all η considered. A complementary analysis of τ η,δ , based on the spatial width Σ(t) (cf. Fig. 12), can be found in Appendix D and provides a similar picture.

V. CONCLUSION
In this paper, we have studied the real-time dynamics of local charge densities in the Fermi-Hubbard chain with a mass-imbalance between the spin-↑ and -↓ particles. To this end, we have prepared a certain class of pure states featuring a sharp initial peak of the density profile for the (lighter) spin-↑ particles in the middle of the chain and investigated the resulting non-equilibrium dynamics. Relying on dynamical quantum typicality, this dynamics can be related to time-dependent correlation functions at equilibrium.
In the regime of weak and moderate imbalance, η 0.6, we have provided evidence for the emergence of diffusive dynamics, manifesting in (i) Gaussian shape of density profiles, (ii) square-root scaling of the spatial variance in time, and (iii) exponentially decaying modes for small momenta.
In contrast, in the regime of strong imbalance, η 0.6, we have observed signatures of anomalous transport, emerging as an exponential rather than a Gaussian shape of density profiles and subdiffusive scaling of spatial variance and density modes in time, consistent with other works [35,40]. However, we cannot rule out that this anomalous transport is just a transient effect which crosses over to normal diffusion at even longer times, e.g., at time scales much longer than the "lifetime" of the Anderson insulator.
For very small but nonzero η, our results are consistent with the absence of genuine localization and support long but finite equilibration times.
Promising future research directions include extensions of the model such as nearest-neighbor interactions and the study of lower temperatures.  to (28), but based on the spatial width Σ(t) (cf. Fig. 12), In comparison to Fig. 14 (b), Fig. 17 provides a very similar picture for the η-dependent lifetime.