Charge diffusion in the one-dimensional Hubbard model

TL;DR

This study investigates charge diffusion in the one-dimensional Hubbard model at high temperatures, revealing robust diffusive behavior across various initial conditions and confirming theoretical predictions with numerical simulations.

Contribution

It demonstrates that charge diffusion in the 1D Hubbard model is robust and occurs under diverse initial conditions, providing numerical evidence supporting linear-response theory.

Findings

Charge diffusion is observed in the strongly interacting regime.

Diffusive behavior is independent of initial state details.

Numerical results agree with linear-response theory.

Abstract

We study the real-time and real-space dynamics of charge in the one-dimensional Hubbard model in the limit of high temperatures. To this end, we prepare pure initial states with sharply peaked density profiles and calculate the time evolution of these nonequilibrium states, by using numerical forward-propagation approaches to chains as long as 20 sites. For a class of typical states, we find excellent agreement with linear-response theory and unveil the existence of remarkably clean charge diffusion in the regime of strong particle-particle interactions. Moreover, we demonstrate that this diffusive behavior does not depend on certain details of our initial conditions, i.e., it occurs for five different realizations with random and nonrandom internal degrees of freedom, single and double occupation of the central site, and displacement of spin-up and spin-down particles.

Equations

$$h_r=-t_{\text{h}}\sum _{s=\downarrow ,\uparrow }(a_{r,s}^{†}{a}_{r+1,s}+\text{H.c.})+U(n_{r,\downarrow }-\frac{1}{2})(n_{r,\uparrow }-\frac{1}{2}),$$

$$|\psi (0)⟩\propto n_{L/2,\uparrow }|\varphi ⟩,|\varphi ⟩=\sum _{k=1}^{4^L}{c}_k|{\psi }_k⟩,$$

$$|{\psi }_{\text{double}}(0)⟩\propto n_{L/2,\downarrow }|\psi (0)⟩.$$

$$p_{r,\uparrow }(t)-p_{\text{eq.}}=2⟨(n_{L/2,\uparrow }-p_{\text{eq.}})(n_{r,\uparrow }(t)-p_{\text{eq.}})⟩,$$

$$\sigma(t)^{2}=\sum_{r=1}^{L}r^{2}\,\delta p_{r,\uparrow}(t)-\Big{[}\sum_{r=1}^{L}r\,\delta p_{r,\uparrow}(t)\Big{]}^{2}\,,$$

$$\frac{\text{d}}{\text{d}t}\sigma (t{)}^2=2D(t),D(t)=\frac{4}{L}{\int }_0^t\text{d}{t}^{\prime }⟨j_{\uparrow }(t^{\prime })j_{\uparrow }⟩,$$

$$p_r(t)-p_{\text{eq.}}=\frac{1}{2}\frac{1}{\sqrt{2\pi }\Sigma (t)}\exp \left[-\frac{(r-L/2{)}^2}{2\Sigma (t{)}^2}\right],$$

$$C_{r,\uparrow }(t)=2⟨(n_{L/2,\uparrow }-p_{\text{eq.}})(n_{r,\uparrow }(t)-p_{\text{eq.}})⟩+p_{\text{eq.}},$$

$$C_{r,\uparrow }(t)=2⟨n_{L/2,\uparrow }{n}_{r,\uparrow }(t)⟩=2\frac{\text{tr}[n_{L/2,\uparrow }{n}_{r,\uparrow }(t)]}{4^L}$$

$$C_{r,\uparrow }(t)=2\frac{\text{tr}[n_{L/2,\uparrow }{n}_{r,\uparrow }(t)n_{L/2,\uparrow }]}{4^L}.$$

$$C_{r,\uparrow }(t)=2\frac{⟨\Phi |n_{L/2,\uparrow }{n}_{r,\uparrow }(t)n_{L/2,\uparrow }|\Phi ⟩}{⟨\Phi |\Phi ⟩}+ϵ,$$

$$C_{r,\uparrow }(t)=2\frac{⟨n_{L/2,\uparrow }\Phi |n_{r,\uparrow }(t)|n_{L/2,\uparrow }\Phi ⟩}{⟨\Phi |\Phi ⟩}$$

$$C_{r,\uparrow }(t)=\frac{⟨e^{-ıHt}{n}_{L/2,\uparrow }\Phi |n_r|e^{-ıHt}{n}_{L/2,\uparrow }\Phi ⟩}{⟨\Phi |\Phi ⟩/2}.$$

$$C_{r,\uparrow }(t)=⟨\psi (t)|n_{r,\uparrow }|\psi (t)⟩=p_{r,\uparrow }(t).$$

$$p_{r,\uparrow }(t)-p_{\text{eq.}}=2⟨(n_{L/2,\uparrow }-p_{\text{eq.}})(n_{r,\uparrow }(t)-p_{\text{eq.}})⟩,$$

$$p_{r,\uparrow }(t)-p_{\text{eq.}}$$

$$⟨n_{L/2,\downarrow }{n}_{L/2,\uparrow }{n}_{L/2,\uparrow }(t)⟩$$

$$=⟨(1-n_{L/2,\downarrow })n_{L/2,\uparrow }{n}_{L/2,\uparrow }(t)⟩$$

$$2⟨n_{L/2,\downarrow }{n}_{L/2,\uparrow }{n}_{L/2,\uparrow }(t)⟩$$

$$=⟨n_{L/2,\uparrow }{n}_{L/2,\uparrow }(t)⟩.$$

$$p_{\text{eq.}}=\frac{L-1}{2L-1}\ne \frac{1}{2}$$

$${ϵ}_r(t)=U(n_{r,\downarrow }-\frac{1}{2})(n_{r,\uparrow }-\frac{1}{2}).$$

$$|{\psi }_{\text{displaced}}(0)⟩=n_{L/2-\delta r,\downarrow }|\psi (0)⟩,$$

Full text

Introduction. Static properties of integrable quantum many-body systems are well understood johnston2000 . In contrast, dynamical questions in these systems continue to be a major challenge in many areas of modern physics and range from fundamental questions in statistical physics to applied questions for a specific class of materials. On the one hand, integrable systems feature a macroscopic number of (quasi-)local conservation laws zotos1997 ; prosen2011 ; prosen2013 and any overlap with these conserved quantities leads to the breakdown of conventional equilibration and thermalization rigol2007 ; ilievski2015 . On the other hand, such overlap is not warranted for all possible initial states, observables, or model parameters and integrability does not rigorously exclude the existence of thermodynamic relaxation such as exponential decay or diffusive transport. This type of relaxation, however, is often traced back to chaos mejiamonasterio2005 ; dalessio2016 , being absent in integrable systems.

In this context, two important and extensively studied examples are (i) the one-dimensional XXZ spin-$1/2$ model and (ii) the (Fermi-)Hubbard chain. As typical for integrable systems, the energy current is (i) strictly or (ii) at least partially conserved zotos1997 such that energy flow is ballistic at any finite temperature kluemper2002 ; langer2011 ; karrasch2014-1 ; karrasch2016-1 , as signaled by a nonzero Drude weight within linear-response theory. However, a much richer dynamical phase diagram develops for other transport quantities. In case (i) of the XXZ spin-$1/2$ chain, the spin current is not strictly conserved. While the partial conservation of this current and a nonzero Drude have been proven analytically below the isotropic point zotos1999 ; benz2005 ; prosen2011 ; prosen2013 , strong numerical evidence for a vanishing Drude weight and nonballistic dynamics has been provided above this point heidrichmeisner2003 ; herbrych2011 ; karrasch2012 ; steinigeweg2014-1 ; steinigeweg2015 . In fact, for the latter regime, clear signatures of diffusion have been reported in various works znidaric2011 ; steinigeweg2011 ; karrasch2014-1 ; steinigeweg2016 . In case (ii) of the Hubbard chain, the situation appears to be similar for charge transport. Even though clarifying the existence of a nonzero Drude weight has turned out to be hard task analytically fujimoto1997 ; kirchner1999 ; peres2000 ; carmelo2013 , numerical studies point to a vanishing Drude weight for strong particle-particle interactions prelovsek2004 ; karrasch2014-2 ; jin2015 . While signatures of diffusion have been observed for such interactions also jin2015 ; karrasch2016-2 ; prosen2012 ; prosen2014 , a direct detection of the characteristic Gaussian broadening is lacking.

The intention of our Letter is to clarify the existence of charge diffusion in the Hubbard chain. For this purpose, we study the nonequilibrium dynamics as resulting for a convenient class of initial states. These initial states are pure and realize density profiles where a peak with the maximum amplitude possible is located in the center of the chain and lies on top of a homogeneous many-body background. First, we focus on a subclass with random internal degrees of freedom and rely on the well-known concept of typicality gemmer2003 ; goldstein2006 ; popescu2006 ; reimann2007 ; bartsch2009 ; bartsch2011 ; sugiura2012 ; sugiura2013 ; elsayed2013 ; hams2000 ; iitaka2003 ; iitaka2004 ; white2009 ; monnai2014 to obtain the real-time broadening of density profiles in the linear-response regime. In this regime, our large-scale numerical simulations for chains as long as $20$ sites allow us to unveil the existence of remarkably clean charge diffusion, as a key result of our Letter. Finally, we extend our analysis to initial states without any randomness and show that the dynamical behavior is stable against varying details of the initial conditions. This stability is another central result of our work and reveals that exactly the same charge diffusion emerges in a far-from-equilibrium setup. These findings clearly demonstrate that thermodynamic relaxation can occur in integrable systems.

Model and Observables. In one spatial dimension and with periodic boundary conditions, the Hamiltonian of the Hubbard model reads $H=\sum_{r=1}^{L}h_{r}$ ,

[TABLE]

where the operator $a_{r,s}^{\dagger}$ ($a_{r,s}$) creates (annihilates) at site $r$ a fermion with spin $s$, $t_{\text{h}}$ is the hopping matrix element, and $L$ is the number of sites. The operator $n_{r,s}=a_{r,s}^{\dagger}a_{r,s}$ is the local occupation number and $U$ is the strength of the on-site interaction. For all values, this model is integrable and the total particle numbers $N_{s}=\sum_{r}n_{r,s}$ and $N=N_{\downarrow}+N_{\uparrow}$ are strictly conserved quantities. We do not restrict ourselves to a particular particle sector, i.e., we study the case $\langle N\rangle=L$ SM . It is worth mentioning that, via the Jordan-Wigner transformation, this model can be mapped onto a spin-$1/2$ two-leg ladder of length $L$, with XY exchange in the legs and Ising exchange in the rungs prosen2012 . In fact, this spin model is used in our numerical simulations.

We are interested in the real-time dynamics of the local occupation numbers $n_{r,s}$ and investigate the expectation values $p_{r,s}(t)=\text{tr}[n_{r,s}\,\rho(t)]$ for the density matrix $\rho(t)$ at time $t$. (It is important to note that $t\neq t_{\text{h}}$.) Doing so, we can follow the broadening of nonequilibrium density profiles, as realized by the preparation of a proper initial state $\rho(0)$.

Initial States. In this Letter, we prepare pure initial states $\rho(0)=|\psi(0)\rangle\langle\psi(0)|$. To specify our $|\psi(0)\rangle$, it is convenient to consider the common eigenbasis of all $n_{r,s}$. Let $|\varphi_{k}\rangle$ be this basis. Then, our initial states read

[TABLE]

where $c_{k}$ are complex and yet arbitrary coefficients. Since $n_{L/2,\uparrow}$ projects only onto states with a spin-up particle in the middle of the chain, $p_{L/2,\uparrow}(0)=1$ has the maximum value possible.

In this Letter, we focus on two particular choices for the coefficients $c_{k}$. First, we choose all $c_{k}$ to be the same number. Second, we choose the $c_{k}$ at random according to the unitary invariant Haar measure bartsch2009 . This choice means in practice that real and imaginary part of the $c_{k}$ are independently drawn from a Gaussian distribution with zero mean. For both, the random and equal choice of the $c_{k}$, all $p_{r\neq L/2,s\neq\uparrow}(0)=p_{\text{eq.}}=1/2$ take on their equilibrium value, and still $p_{L/2,\uparrow}(0)=1$. Therefore, the initial density profile features a central peak on top of a homogeneous many-particle background.

A very similar form for the density profile also results for the state

[TABLE]

Then, due to the additional projection, $p_{L/2,\downarrow}(0)=1$ also for equal and random $c_{k}$. Therefore, the density profiles $p_{r,\uparrow}(t)=p_{r,\downarrow}(t)$ are identical for $t=0$ and all later times $t>0$ as well, see SM for initial displacement.

All initial states introduced have to be considered as far-from-equilibrium states: They are not only pure but also have maximum $p_{L/2,\uparrow}(0)=1$. Remarkably, however, the dynamics of $|\psi(0)\rangle$ in Eq. (2) with random $c_{k}$ can be connected to the linear-response Kubo formula, since the underlying $|\phi\rangle$ is a so-called typical state steinigeweg2016 , see also Refs. gemmer2003 ; goldstein2006 ; popescu2006 ; reimann2007 ; bartsch2009 ; bartsch2011 ; sugiura2012 ; sugiura2013 ; elsayed2013 ; hams2000 ; iitaka2003 ; iitaka2004 ; white2009 ; monnai2014 for the concept of typicality. Exploiting this typicality allows one to derive the relation SM

[TABLE]

where $\langle\bullet\rangle=\text{tr}[\bullet]/4^{L}$ is the thermodynamic average at formally infinite temperature. Thus, for a typical state, the nonequilibrium expectation value is directly related to an equilibrium correlation function. This fact enables a connection to the Kubo formula via the variance

[TABLE]

where $\delta p_{r,\uparrow}(t)=2(p_{r,\uparrow}(t)-p_{\text{eq.}})$ excludes the equilibrium background and is normalized to $\sum_{r}\delta p_{r,\uparrow}(t)=1$. As shown in Ref. steinigeweg2009 , the time derivative of this variance satisfies

[TABLE]

where $j_{\uparrow}=-t_{\text{h}}\sum_{r}(\imath a_{r,\uparrow}^{\dagger}a_{r+1,\uparrow}+\text{H.c.})$ is the total current of the spin-up particles and the quantity $D(t)$ plays the role of a time-dependent diffusion coefficient. For $U=0$, [$j_{\uparrow},H]=0$ necessarily leads to ballistic scaling $D(t)\propto t$ and $\sigma(t)\propto t$. For large $U\gg t_{\text{h}}$, signatures of diffusive scaling $D(t)=\text{const.}$ and $\sigma(t)\propto\sqrt{t}$ have been reported in Refs. jin2015 ; karrasch2016-2 . So far, however, a systematic analysis beyond the mere width of the density profile is lacking and the central issue of our Letter.

Numerical Technique and Results. From a numerical point of view, the Hubbard chain is challenging since the Hilbert-space dimension $\text{dim}=4^{L}$ grows rapidly with $L$, e.g., much faster than the also exponential increase $\text{dim}=2^{L}$ in case of a spin-$1/2$ chain. As a consequence, exact diagonalization of the Hamiltonian is only feasible for a few lattice sites and a real-space experiment like the one done in our Letter would not be reasonable. Hence, we proceed differently and profit from the fact that we only need to deal with pure states. The time evolution of these states can be obtained by forward-propagation methods such as fourth-order Runge-Kutta steinigeweg2014-1 ; steinigeweg2015 ; elsayed2013 or more sophisticated schemes such as Trotter decompositions or Chebyshev polynomials steinigeweg2014-2 ; jin2015 ; weisse2006 . We apply a second-order Trotter formula with a time step $\delta t\,t_{\text{h}}=0.05$, sufficient to reach very good agreement with Chebyshev-polynomial algorithms. A massively parallelized implementation of this formula allows us to treat Hubbard chains as long as $L=20$ sites. For this system size and a maximum time $t\,t_{\text{h}}=8$, the simulation takes about $9$ hours when using $262,144$ double-thread cores. Therefore, apart from the $L=20$ data depicted in Fig. 2, we focus on $L=18$ to reduce computational costs at least a bit.

Now, we turn to our numerical results. We start with the initial state $|\psi(0)\rangle$ in Eq. (2) with a random choice of the coefficients $c_{k}$, i.e., a typical state. It is important to note that we consider a single realization of the $c_{k}$ and do not perform any kind of averaging. Still, we allow the $c_{k}$ to be different for each simulation. In Figs. 1 (a)-(c) we depict our results for the spin-up occupation numbers $p_{r,\uparrow}$(t) for a Hubbard chain of length $L=18$ and with different interactions $U/t_{\text{h}}=16$, $4$, [math], in a 2D time-space density plot. Several comments are in order. First, for the noninteracting case $U=0$ in Fig. 1 (c), the real-time broadening of $p_{r,\uparrow}(t)$ is clearly linear and, as discussed above, has to occur due to the strict conservation of the particle current. The pronounced jets visible are typical for free-particle cases karrasch2014-1 ; steinigeweg2016 and propagate fast without any scattering until they eventually hit the boundary of the chain at short times $t\,t_{\text{h}}\sim 4$. Second, for the interacting cases $U/t_{\text{h}}=16$, $4$ in Figs. 1 (a), (b), these jets and the linear broadening as well disappear, i.e., the dynamics is not ballistic. Note that, in contrast, the flow of energy is ballistic for arbitrary $U$ SM . Third, the broadening is the slower the larger $U$ because scattering becomes stronger as $U$ increases. In particular, for the largest $U/t_{\text{h}}=16$ in Fig. 1 (a) and the maximum time $t\,t_{\text{h}}=8$ calculated, the overall width of $p_{r,\uparrow}(t)$ is still smaller than the chain length. Therefore, we can exclude trivial finite-size effects for such times SM .

To gain insight into the dynamics at large $U/t_{\text{h}}=16$, we show in Fig. 2 (a) the site dependence of the profile $p_{r,\uparrow}(t)$ for various times $t\,t_{\text{h}}=0$, $1$, $2$, $4$, and for an even larger system size $L=20$. We do so by subtracting from $p_{r,\uparrow}(t)$ the equilibrium value $p_{\text{eq.}}$ and using a semi-log plot, to visualize also the outer tails of the profile. It is intriguing to see that, for all times $t$ depicted, the profiles can be described very well by Gaussians

[TABLE]

where the standard deviation $\Sigma(t)$ occurs as the only free parameter and is adjusted by fitting. The excellent fits over several orders of magnitude are a central result of our Letter and already provide strong evidence for the existence of diffusion. Still, however, $\Sigma(t)$ needs to scale as $\Sigma(t)\propto\sqrt{t}$.

For a final conclusion, we show in Fig. 2 (b) the time dependence of $\Sigma(t)$ and compare to the linear-response $\sigma(t)$ in Eq. (6), as resulting from the $D(t)$ calculated in Ref. jin2015 for the same interaction $U/t_{\text{h}}=16$ and length $L=18$. While the perfect agreement illustrates the high accuracy of the typicality relation, this agreement implies that the linear-response result $D(t)=\text{const.}$ jin2015 ; karrasch2016-2 also holds for our nonequilibrium dynamics. Thus, diffusion clearly exists. Note that the same conclusion can be drawn for smaller interactions $U/t_{\text{h}}=8$ also SM , where finite-size effects are still negligibly small. For $U/t_{\text{h}}\ll 8$, however, significant finite-size effects are known to occur jin2015 and a reliable conclusion on the thermodynamic limit $L\to\infty$ is impossible on the basis of $L\sim 20$.

Next, we intend to shed light on the role of the specific initial-state realization, in particular on the influence of randomness. Therefore, in a first step, we investigate the initial state $|\psi(0)\rangle$ in Eq. (2) again but now with equal coefficients $c_{k}$. Recall that, while this nonrandom state has exactly the same initial density profile, the typicality relation does not need to hold any further. For this state, we repeat the calculation in Fig. 1 for $L=18$ sites and different interactions $U/t_{\text{h}}=16$, $4$, [math] and summarize the corresponding results in Fig. 3. In comparison to Fig. 1, jet-like behavior for $U=0$ is enhanced in Fig. 3 (c) and emerges now for $U/t_{\text{h}}=4$ in Fig. 3 (b) in addition. The same observation has been made for the XXZ spin-$1/2$ chain below the isotropic point steinigeweg2016 . Remarkably, however, the diffusive behavior for $U/t_{\text{h}}=16$ in Fig. 3 (a) turns out to be unaltered. In fact, this observation is different to the one found for the XXZ spin-$1/2$ chain above the isotropic point steinigeweg2016 , where the impact of nonrandomness is strong.

The above finding suggests that charge diffusion for strong interactions $U$ is stable against varying details of the initial condition. To substantiate this suggestion, we finally extend our analysis to the initial state $|\psi_{\text{double}}(0)\rangle$ in Eq. (3) and study both, random and equal coefficients $c_{k}$. To repeat, these states have the same initial density profile $p_{r,\uparrow}(0)$ but $p_{r,\downarrow}(0)=p_{r,\uparrow}(0)$ now. For system size $L=18$, interaction $U/t_{\text{h}}=16$, and two different times $t\,t_{\text{h}}=1$, $4$, we compare in Fig. 4 the distribution $p_{r,\uparrow}(t)$ for these states with the one for the others. Apparently, $p_{r,\uparrow}(t)$ is practically indistinguishable for the two cases with random $c_{k}$. Even though not shown explicitly here, these two random cases also coincide for other values of $U$ SM . While the two equal cases in Fig. 4 differ from the two random ones, this difference is minor in view of the semi-log plot used. These observations are another central result of our Letter and clearly show that charge diffusion in the strong-interaction limit does not depend on the specific initial-state preparation, at least for the whole class of nonequilibrium states investigated.

Conclusions. In this Letter, we have investigated the real-time broadening of charge in the Hubbard chain at high temperatures. First, we have introduced a class of pure initial states with density profiles where a sharp peak is located in the middle of the chain and lies on top of a homogeneous many-particle background. Then, we have calculated the dynamics of these nonequilibrium states, by using large-scale numerical simulations. Our results for typical states have unveiled the existence of remarkably clean charge diffusion in the limit of strong particle-particle interactions, in perfect agreement with the Kubo formula. We have additionally shown that this diffusive behavior is stable against varying details of the initial conditions.

Acknowledgments. We sincerely thank T.  Prosen and F. Heidrich-Meisner for fruitful discussions. In addition, we gratefully acknowledge the computing time, granted by the “JARA-HPC Vergabegremium” and provided on the “JARA-HPC Partition” part of the supercomputer “JUQUEEN” stephan2015 at Forschungszentrum Jülich.

I Supplemental Material

I.1 Typicality Approximation

I.1.1 Single Projection

To make our Letter self-contained, we provide details on the calculation yielding the relation in Eq. (4) of the main text. A very similar calculation can be found in Ref. steinigeweg2016 for the XXZ spin-$1/2$ chain.

Starting with the correlation function

[TABLE]

multiplying out the two brackets, and using $p_{\text{eq.}}=1/2$ as well as $\langle n_{r,\uparrow}(t)\rangle=p_{\text{eq.}}$, we get

[TABLE]

and, due to $n_{L/2,\uparrow}^{2}=n_{L/2,\uparrow}$ and a cyclic permutation in the trace, we obtain

[TABLE]

If $|\Phi\rangle$ is a random pure state according to the unitary invariant Haar measure, this correlation function can be expressed as

[TABLE]

where the negligibly small error $\epsilon\propto 4^{L/2}$ is skipped in the following for clarity. Because of $n_{L/2,\uparrow}^{\dagger}=n_{L/2,\uparrow}$, we can rewrite this expression as

[TABLE]

and, applying $n_{r,\uparrow}(t)=e^{\imath Ht}\,n_{r,\uparrow}\,e^{-\imath Ht}$ and moving the factor $2$ to the denominator, this expression becomes

[TABLE]

Finally, since $|\psi(0)\rangle=n_{L/2,\uparrow}\,|\Phi\rangle/\sqrt{\langle\Phi|\Phi\rangle/2}$, we end up with

[TABLE]

Comparing Eqs. (S1) and (S7) leads to

[TABLE]

i.e., the relation in Eq. (4) of the main text.

I.1.2 Double Projection

For the initial state $|\psi_{\text{double}}(0)\rangle=n_{L/2,\downarrow}\,|\psi(0)\rangle$, one can simply repeat the steps in Eqs. (S1) - (S7) to show the relation

[TABLE]

Thus, the resulting relation for $|\psi_{\text{double}}(0)\rangle$ differs from the one for $|\psi(0)\rangle$. However, both relations are identical for $U=0$ (and arbitrary $t$) or for $t=0$ (and arbitrary $U$). For other values of $U$ and $t$, we can use symmetries of our specific Hamiltonian in Eq. (1), i.e., particle-hole symmetry of each particle species. Due to the symmetry for spin-down particles, we get

[TABLE]

or, equivalently,

[TABLE]

Using this identity, Eq. (S9) simplifies to Eq. (4) and the relations for $|\psi(0)\rangle$ and $|\psi_{\text{double}}(0)\rangle$ coincide, at least for the Hubbard chain.

Numerically, we have observed in Fig. 4 of our Letter that the dynamics resulting for the two random initial states $|\psi_{\text{double}}(0)\rangle$ and $|\psi(0)\rangle$ are practically the same, at least for the case of a strong interaction $U/t_{\text{h}}=16$. To demonstrate that no difference is found for other values of $U$ also, we repeat the full calculation in Fig. 1 of the main text for $|\psi_{\text{double}}(0)\rangle$ and depict the corresponding results in Fig. S1. Apparently, the real-time broadening does not change for this initial state, in agreement with the analytical arguments above.

I.2 Intermediate Interactions

In Fig. 2 of our Letter, we have shown the pronounced Gaussian form of the density profile $p_{r,\uparrow}(t)$ for different times $t$ but only for a representative interaction strength $U/t_{\text{h}}=16$. To demonstrate that this form also emerges for other values of $U$, we repeat in Fig. S2 the analysis for the intermediate value $U/t_{\text{h}}=8$. As before, the Gaussian fits in Fig. S2 (a) describe the data very well over several orders of magnitude and the agreement of the resulting standard deviation with linear response theory in Fig. S2 (b) is excellent.

I.3 Chain Length and Particle Sector

All data presented in the main text corresponds to the largest system size available, i.e., $L=18$ (and $L=20$ for a single set of parameters as well). In addition, we have performed a careful finite-size analysis to ensure that these data allow us to draw reliable conclusions on the thermodynamic limit $L\to\infty$. Thus, as an example, we present in Fig. S3 (a) such an analysis for a random initial state $|\psi(0)\rangle$, interaction strength $U/t_{\text{h}}=16$, and the longest time $t\,t_{\text{h}}=4$ considered. Clearly, the density profile $p_{r,\uparrow}(t)$ depicted is practically the same for lengths $L=14$, $18$, $20$ and even the outer tails of this profile do not have a significant finite-size dependence.

Further, it is worth mentioning that we have checked the independence of the specific ensemble used, i.e., our results do not depend on our choice $\langle N\rangle=L$. In Fig. S3 (b) we demonstrate this independence by a comparison of the ensembles $\langle N\rangle=L$ in the average and exactly $N=L$ for the $L=14$ curve in Fig. S3 (a), where the finite-size difference should be strongest. Still, we do not find any significant difference. However, we should stress that we eliminated trivial finite-size effects for the $N=L$ case by taking into account that

[TABLE]

for any finite $L$.

I.4 Energy Flow

In the main text, we have unveiled the emergence of remarkably clean particle diffusion for sufficiently large interaction strengths $U\gg t_{\text{h}}$. But this diffusive behavior should be contrasted against the dynamical behavior of other quantities for exactly the same initial states. Such a quantity occurs in the second term of the Hamiltonian (1), i.e., the local interaction energy

[TABLE]

This quantity becomes the full local energy in the limit of large $U$. Thus, the flow is expected to be ballistic karrasch2016-1 , as typical for integrable systems. In Fig. S4 we confirm this expectation for, e.g., a random initial state $|\psi_{\text{double}}(0)\rangle$, where the central interaction energy $\epsilon_{L/2}(0)=U/4$ takes on the maximum value possible. Apparently, the jet-like broadening of energy observed for $U/t_{\text{h}}=16$ in Fig. S4 (a) is very similar to the one of charge found for $U=0$ in Fig. S1 (c). These findings also illustrate that our initial states do no enforce diffusive dynamics.

I.5 Displacement of Particle Species

For simplicity, our work has focused on the two cases: 1. The spin-down density profiles are $p_{r,\downarrow}(t)=p_{\text{eq.}}$, i.e., they are at equilibrium from the very beginning (initial state $|\psi(0)\rangle$). 2. These profiles are $p_{r,\downarrow}(t)=p_{r,\uparrow}(t)$, i.e., they are always identical to the density profiles of the other particle species (initial state $|\psi_{\text{double}}(0)\rangle$). This is why we have not shown results for $p_{r,\downarrow}(t)$ explicitly in the figures.

Certainly, it is possible to investigate other situations also. For instance, similar to $|\psi_{\text{double}}(0)\rangle$, one may study initial states

[TABLE]

where $p_{r,\downarrow}(0)$ is not peaked in the middle of the chain and displaced to the left by $\delta r$ sites. In such a situation, $p_{r,\downarrow}(t)\neq p_{r,\uparrow}(t)$ and also the dynamical behavior could change in principle. However, first results indicate that charge dynamics is qualitatively the same even for this situation. As an example, we show in Fig. S5 results for a random initial state $|\psi_{\text{displaced}}(0)\rangle$ with a displacement by $\delta r=4$ sites, for a Hubbard chain of length $L=18$ and with interaction $U/t_{\text{h}}=16$. Here, it turns out that the two particle species spread independently from each other, and no significant difference to Fig. 1 (a) of the main text is visible. Once again, this observation can be understood analytically using typicality and arguments analogous to the ones in the context of Eq. (S9).