Diffusion in deterministic interacting lattice systems

TL;DR

This paper rigorously analyzes deterministic lattice systems with charged particles, revealing ballistic, diffusive, and insulating transport behaviors, and provides exact formulas for transport coefficients validated by simulations.

Contribution

It offers the first exact analytical expressions for current autocorrelation and transport coefficients in a deterministic lattice model with interactions.

Findings

Identifies three transport regimes: ballistic, diffusive, insulating.

Derives exact formulas for diffusion constant and Drude weight.

Validates analytical results with Monte-Carlo simulations.

Abstract

We study reversible deterministic dynamics of classical charged particles on a lattice with hard-core interaction. It is rigorously shown that the system exhibits three types of transport phenomena, ranging from ballistic, through diffusive to insulating. By obtaining an exact expressions for the current time-autocorrelation function we are able to calculate the linear response transport coefficients, such as the diffusion constant and the Drude weight. Additionally, we calculate the long- time charge profile after an inhomogeneous quench and obtain diffusive profile with the Green-Kubo diffusion constant. Exact analytical results are corroborated by Monte-Carlo simulations.

Figures

Equations

$$(\varnothing ,\varnothing )\leftrightarrow (\varnothing ,\varnothing ),(\varnothing ,\alpha )\leftrightarrow (\alpha ,\varnothing ),(\alpha ,\beta )\leftrightarrow (\alpha ,\beta ),$$

$${\varphi }^{\text{o}}={\varphi }_{1,2}\circ \cdots \circ {\varphi }_{n-1,n},{\varphi }^{\text{e}}={\varphi }_{2,3}\circ \cdots \circ {\varphi }_{n,1}.$$

$$q_x^t=s_x^t+s_{x+1}^t,$$

$$j_x^{t+1/2}=2(-1{)}^x(s_x^{t+1/2}-s_{x+1}^{t+1/2})(s_x^{t+1/2}+s_{x+1}^{t+1/2}{)}^2,$$

$$q_x^{t+1}-q_x^t+\frac{1}{2}\left(j_{x+1}^{t+1/2}-j_{x-1}^{t+1/2}\right)=0$$

$$[\alpha {]}_x(\underline{s})={\delta }_{{\alpha }_x,s_x},\alpha \in \{\varnothing ,+,-\},$$

$$[{\alpha }_1{\alpha }_2\dots {\alpha }_r{]}_x=[{\alpha }_1{]}_x[{\alpha }_2{]}_{x+1}\cdots [{\alpha }_r{]}_{x+r-1}.$$

$$U=U^{\mathrm{o}}{U}^{\mathrm{e}},U^{\mathrm{o}}=\prod _{x=1}^{n/2}{U}_{2x-1,2x},U^{\mathrm{e}}=\prod _{x=1}^{n/2}{U}_{2x,2x+1}.$$

$$⟨a⟩=\sum _{\underline{s}}a(\underline{s}),$$

$$[0]$$

$$[1]$$

$$[2]$$

$$q_x={[10]}_x+{[01]}_x,$$

$$j_x=2(1-\rho )(-1{)}^x\left({[10]}_x-{[01]}_x+{[12]}_x-{[21]}_x\right).$$

$$\sigma =\frac{1}{2}C(0)+\sum _{t=1}^{\infty }C(t),$$

$$p_x=[0{]}_x^{\prime }+\mu [1{]}_x^{\prime },0\le \rho \le 1,-\rho \le \mu \le \rho .$$

$$⟨[0][1]{⟩}_p=⟨[0][2]{⟩}_p=⟨[1][2]{⟩}_p=0.$$

$$\mathcal{U}=\left[\begin{array}{ccc}1-2\rho & 2\rho & 0\\ 0& 0& 0\\ -2(1-\rho )& 1-2\rho & 0\end{array}\right],$$

$$C(t)$$

$$=\{\begin{array}{ll}8\rho (2-\rho )(1-\rho );& t=0,\\ 16\rho (1-\rho {)}^4(1-2\rho {)}^{2t-2};& t\ge 1,\end{array}$$

$$\sigma =4(1-\rho ),\mathcal{D}={\rho }^{-1}-1,$$

$$z_d^k=\sum _x([0\underset{k}{\underbrace{0\dots 0}}1\underset{d}{\underbrace{0\dots 0}}2{]}_{2x}-[02\underset{d}{\underbrace{0\dots 0}}1\underset{k}{\underbrace{0\dots 0}}{]}_{2x}).$$

$$\frac{C(t)}{8\bar{\rho }}=\{\begin{array}{ll}{\mu }^2+\rho (2-\rho );& t=0,\\ \frac{2{\mu }^2}{\rho }+2{\bar{\rho }}^3(1-2\rho {)}^{2t-2}\left(\rho -\frac{{\mu }^2}{\rho }\right);& t\ge 1,\end{array}$$

$$D=16\left({\rho }^{-1}-1\right){\mu }^2,$$

$$p_{n/2+x}=[0]^{\prime}+\mu_{x}[1]^{\prime},\$$

$$f(x,t)=⟨U^t{q}_{n/2+2x-1}{⟩}_p=⟨q_{n/2+2x-1}{⟩}_{U^{-t}p}.$$

$$\tilde{p}=\sum _{x=-\frac{n}{2}+1}^{\frac{n}{2}}{\mu }_x{e}_x,\text{where}{e}_x=[\underset{\frac{n}{2}+x-1}{\underbrace{0\dots 0}}1\underset{\frac{n}{2}-x}{\underbrace{0\dots 0}}{]}^{\prime }.$$

$$\mathcal{U}=\left[\begin{array}{cccccc}\ddots & \ddots & \ddots & & & \\ & c& a& b& & \\ & & c& a& b& \\ & & & \ddots & \ddots & \ddots \end{array}\right],$$

$$\begin{array}{rl}f(x,t)=\sum _{y=x-t}^{x+t}& {\int }_{-\pi }^{\pi }\mathrm{d}k{e}^{\mathrm{i}k(x-y)}{\mu }_{2y}\cdot \\ & \cdot \left({\lambda }_2^t(k){\tilde{\alpha }}_2(k)-{\lambda }_1^t(k){\tilde{\alpha }}_1(k)\right),\end{array}$$

$$\tilde{f}(\xi )=({\mu }_R+{\mu }_L)+({\mu }_R-{\mu }_L)\text{erf}\left(\frac{\xi }{4\sqrt{\gamma }}\right).$$

$$\mathcal{D}=({\rho }^{-1}-1)\frac{1-\rho \Gamma }{1-\Gamma },$$

$$K(h)=1+h\sum _{x=1}^{n-1}(x-\frac{n}{2})q_x,\text{with}h<\frac{1}{n},$$

$$\tilde{p}=K(h)p,$$

Full text

Introduction.– Understanding out-of-equilibrium phenomena has been at the forefront of condensed matter physics of the last decade. Despite the efforts, we have only recently gained a considerable insight into the microscopic origins of the transport in interacting systems in terms of an emerging field of generalized hydrodynamics (GHD) Spohn (2012); Bernard and Doyon (2016); Castro-Alvaredo et al. (2016); Bertini et al. (2016); Sotiriadis ; De Luca et al. (2016); Ilievski and De Nardis (2017); Bulchandani et al. (2017); Doyon and Spohn (2017); Doyon et al. (2017a, b). While GHD provides a general framework to analytically deal with the ballistic transport, in particular in integrable systems, it lacks an extension which would enable us to study normal, diffusive transport. In this regard, there are only a few results on lower bounding the diffusion constant in terms of local, or bilocal conserved charges Prosen (2014); Medenjak et al. (2017).

Since the integrable systems are characterized by ballistically propagating excitations Bonnes et al. (2014), the question of how the diffusion, which is usually related to microscopic chaos, arises in integrable, locally interacting, clean and deterministic (i.e. non-disordered and non-dissipative) systems is puzzling to say the least Ljubotina et al. (2017); Prosen and Žnidarič (2009); Steinigeweg and Brenig (2011); Steinigeweg (2011); Steinigeweg et al. (2014); Karrasch et al. (2014). A great deal of attention has been devoted to the study of inhomogeneous quench problems, where two chains in equilibrium with distinct temperatures or chemical potentials are joined together and then let to evolve under a homogeneous Hamiltonian Spohn (2012); Bernard and Doyon (2016); Castro-Alvaredo et al. (2016); Bertini et al. (2016); Sotiriadis ; De Luca et al. (2016); Ilievski and De Nardis (2017); Bulchandani et al. (2017); Doyon and Spohn (2017); Spohn and Lebowitz (1977). In such situations, non-equilibrium steady states typically emerge on ballistic lines $x=vt$. However, for systems exhibiting parity-like symmetries, with respect to which the initial quench state is symmetric and the current anti-symmetric, the ballistic transport channel may close and time-dependent DMRG simulations clearly indicate that the steady state arises along diffusive lines $x=\xi\sqrt{t}$ Ljubotina et al. (2017).

In this Letter we expound one possible general mechanism for diffusive behavior in interacting systems, which as we show on an example, is connected to the interplay between freely propagating neutral degrees of freedom and insulating behavior of charge carrying ones. The model in question is a simple, reversible cellular automaton, consisting of three types of particles: freely moving vacancies, and hard-core interacting positive and negative charges. Despite the integrability, obtaining the full time dependence is usually intractable, except for non-interacting systems. Typically this presents an insurmountable barrier for calculation of transport coefficients. In our model, however, we explicitly compute the time dependence of current time autocorrelation functions in separable equilibrium states, and solve the inhomogeneous quench problem with arbitrary initial charge density bias resulting in a universal diffusive error function scaling profile. Depending on the density of vacancies and the imbalance of positive and negative charge three different regimes are identified. The absence of vacancies renders the system insulating, while in a generic case of the charge imbalance the system exhibits ideal transport. The regime with a finite density of vacancies and without the charge imbalance is especially interesting, since it includes the maximum entropy state, and in this regime the model exhibits purely diffusive transport.

The model.– Consider a deterministic, reversible cellular automaton defined on the chain with  an even length of $n$ sites. Each site can be either vacant (state $\emptyset$), or occupied by a positively or negatively charged particle (state $+$ or $-$). The dynamics of the lattice configuration, ${\underline{s}}=(s_{1},\ldots,s_{n}),s_{x}\in\{\emptyset,+,-\}$, can be expressed in terms of a local two site mapping $\phi_{x,x+1}({\underline{s}})=(s_{1},s_{2},\dots,s^{\prime}_{x},s^{\prime}_{x+1},\dots,s_{n})$, where the updated elements $(s_{x}^{\prime},s_{x+1}^{\prime})$ are obtained from the initial ones $(s_{x},s_{x+1})$, by a self-inverse transformation

[TABLE]

with $\alpha,\beta\in\{+,-\}$. The local process describes the elastic scattering of charged particles. The lattice configuration at time $t,$ ${\underline{s}}^{t}=\phi({\underline{s}}^{t-1})$, can be expressed in terms of a two step propagator $\phi=\phi^{\text{o}}\circ\phi^{\text{e}}$ given by sequences of disjoint local mappings (1)

[TABLE]

The dynamics is induced by freely propagating vacancies, while the clustered particles remain completely frozen in time as illustrated in Fig. 1. The cellular automaton admits a mechanical analogy in terms of a two-species, synchronous hard-point gas if particles at occupied sites $x$ are attributed velocities $2(-1)^{x}$ or $-2(-1)^{x}$, at odd or even steps, corresponding respectively to integer or half-integer times 111 In this sense our model belongs to a class of hard-rod systems, where diffusion has been extensively studied Jepsen (1965); Lebowitz et al. (1968, 1969); Dürr et al. (1985); Balakrishnan et al. (2002). However, considering particle positions to take values on the lattice $\mathbb{Z}$ rather than on the line $\mathbb{R}$ has two conceptual advantages: (i) The information (entropy) density in the initial condition is qualitatively smaller and diffusion cannot be attributed to random phases with which different particles collide even if velocities have a discrete distribution. (ii) Dynamics can be analyzed in terms of simple linear algebra and matrix spectral problems..

We are interested in the dynamics of the charge $q^{t}_{x}$ corresponding to sites $x$ and $x+1$ at time $t$,

[TABLE]

with $s^{t}_{x}=\pm 1$, if the site $x$ is occupied by the positive/negative charge, and $s_{x}^{t}=0$ otherwise. The total charge $Q^{t}=\sum_{x}q^{t}_{x}$, is a constant of motion. The corresponding current $J^{t}=\sum_{x}j_{x}^{t}$ can be defined on the intermediate time steps (see Fig. 2) as

[TABLE]

where $j_{x}$ denotes the local density. The only configurations that contribute to the current are $(\pm,\emptyset)$ and $(\emptyset,\pm)$. One can check that the continuity equation

[TABLE]

is satisfied.

To carry out the calculations we introduce a commutative algebra of observables ${\cal A}$ (multiplicative algebra of functions over the configuration space) with a local basis

[TABLE]

The subscript $x$ will be dropped when the position is clear from the context. The dynamics of observables $a^{t}({\underline{s}})\equiv U^{t}a({\underline{s}})=a(\phi^{t}({\underline{s}}))$ is given by a $3^{n}\times 3^{n}$ matrix $U$ 222note that in the canonical basis $[\alpha_{1},\dots,\alpha_{n}]$, $U$ is a permutation matrix given in terms of local propagators $(U_{x,y})_{{\underline{s}}^{\prime},{\underline{s}}}=\delta_{{\underline{s}}^{\prime},\phi_{x,y}({\underline{s}})}$ factorized as

[TABLE]

Note that the local propagator obeys Yang-Baxter equation $U_{x,y}U_{y,z}U_{x,y}=U_{y,z}U_{x,y}U_{y,z}$. Additionally, it proves useful to introduce the unnormalized maximum entropy state

[TABLE]

in terms of which an expectation value in any probability distribution (state) $p$ can be expressed as $\langle a\rangle_{p}=\langle a\,p\rangle$. We introduce an alternative local basis and its dual

[TABLE]

such that $[0]\equiv 1$ and $\langle[\alpha][\beta]^{\prime}\rangle=\delta_{\alpha,\beta}$, $\alpha,\beta\in\{0,1,2\}$. The parameter $\rho$ will be connected to the density of particles later on. The local charge and current now read

[TABLE]

Linear response.– According to Einstein’s relation the diffusion coefficient $\mathcal{D}$ is connected to conductivity as ${\sigma}=\chi\mathcal{D}$, $\chi$ being the static susceptibility (i.e. the second moment of charge $Q$). In the small constant electric field, the conductivity takes the following form (Sect. A of sup )

[TABLE]

where $C(t)=\lim_{n\to\infty}\frac{1}{n}\langle JU^{t}J\rangle_{p}$ is the current correlation function, with $p$ being an equilibrium state, $Up=p$. Another important transport coefficient, corresponding to the rate at which the conductivity diverges Ilievski and Prosen (2013), is the Drude weight, $D=C(\infty)$.

We restrict the discussion to translationally invariant product equilibrium states $p=p(\rho,\mu)=\prod_{x=1}^{n}p_{x}$,

[TABLE]

Note that $[0]^{\prime}_{x}$ depends on $\rho$. The density $\rho$ represents the probability of a lattice site being occupied by a charged particle, while the chemical potential $\mu$ corresponds to the charge imbalance. The static susceptibility in such a state is $\chi=4(\rho-\mu^{2})$.

Diffusive regime.– Initially we consider a balanced equilibrium state with $\mu=0$ and arbitrary $\rho$ for which the following orthogonality relations hold,

[TABLE]

As a consequence only the observables with a single occurrence of ${[1]}$ and at most one ${[2]}$ in the propagated current $J^{t}$ can contribute to the expectation value (13). The local observables that are relevant for the calculation of the diffusion constant are $y_{0}=\sum_{x}\left({[10]}_{2x}-{[01]}_{2x}\right)$, $y_{1}=\sum_{x}\left({[12]}_{2x}-{[21]}_{2x}\right)$ and $y_{2}=\sum_{x}\left({[012]}_{2x}-{[021]}_{2x}\right)$, in terms of which the current reads $J=2(1-\rho)(-2y_{0}-y_{1}+y_{2})$. In addition, we note that two steps of the propagator are conjugated $U^{\text{o}}=S^{-1}U^{\text{e}}S$ by a lattice shift, defined as $S[{\underline{\alpha}}]_{x}=[{\underline{\alpha}}]_{x+1}$. Since $y_{\alpha}$ are translationally invariant, i.e. $S^{2}y_{\alpha}=y_{\alpha}$, the complete propagator on the relevant subspace reads $U=(SU^{\text{e}})^{2}$. Under the half-time propagation $SU^{\text{e}}$ the observables $y_{0}$ and $y_{1}$ map into linear combinations of $y_{0}$ and $y_{2}$. Due to the ballistic propagation of ${[2]}$, $U_{1,2}{[02]}={[20]}$, $U_{1,2}{[20]}={[02]}$, any additional basis operators appearing in $y^{t}_{2}$, $t>0$, are orthogonal to $J$, since there is always at least one occurrence of ${[0]}$ between ${[1]}$ and ${[2]}$. Therefore, in order to compute $C(t)$, we only have to consider the time evolution restricted to the subspace spanned by $\{y_{\alpha}\}$. Half-time step propagator $SU^{\text{e}}$ projected to this basis reads

[TABLE]

and yields the following simple expression for the current-autocorrelation function

[TABLE]

The row vector contains properly normalized overlaps $\langle y_{\alpha}J\rangle_{p}$. One can then easily calculate the conductivity and the diffusion constant

[TABLE]

which agree with Monte-Carlo simulations (see Fig. 3).

Ballistic regime.– By introducing the charge imbalance $\mu\neq 0$, the local observables cease to be orthogonal, $\langle{[2]}{[1]}\rangle_{p(\rho,\mu)}\neq 0$ . However, we may still consider only a subspace of observables, with at most one occurence of ${[2]}$ and single ${[1]}$. In this case the relevant subspace is spanned by $\{y_{0},y_{1},z_{2d}^{0},z_{2d+1}^{1};d\geq 0\}$, with the observables $z_{d}^{k}$ defined as

[TABLE]

Performing a similar calculation as above we obtain the following expression for the time autocorrelation function (Sect. B of sup )

[TABLE]

with $\bar{\rho}=1-\rho$, which immediately yields the exact expression for the Drude weight

[TABLE]

again agreeing excellently with numerics (Fig. 3).

Inhomogeneous quench.– Let us consider dynamics of charges starting from the product initial state $p$ with uniform density $\rho$ and two distinct chemical potentials $\mu_{L},\,\mu_{R}$ on the left and the right half of the chain, which can be expressed locally as

[TABLE]

with $\mu_{x}=\mu_{L}$ for $x\leq 0$ and $\mu_{R}$ for $x>0$. We shall now discuss the time evolution of the centered charge profile

[TABLE]

We are interested in the dynamics inside of the light-cone $|x|\leq t$ with $n\geq 8\,t+2$, so that the boundary values of the profile are constant $f(\mp t,t)=2\mu_{L,R}$.

Due to the conservation of the number of $\alpha_{j}=1$ in the adjoint time evolution $U^{-t}p$ expressed in the dual basis $[{\underline{\alpha}}]^{\prime}$, the initial state $p$ can be replaced by

[TABLE]

Here we use the dual space version of the argument used in the computation of correlation functions. After half of the time step the time propagated charge density (24) consists of the terms including a single occurrence of $[1]^{\prime}$ and combinations of $[1]^{\prime}$ and $[2]^{\prime}$ on neighboring sites on the background of $[0\ldots 0]^{\prime}$. Using the argument of freely propagating $[2]^{\prime}$, namely $U_{1,2}[02]^{\prime}=[20]^{\prime}$, $U_{1,2}[20]^{\prime}=[02]^{\prime}$, we can conclude that the terms containing $[2]^{\prime}$ can be disregarded at all time steps due to the orthogonality to charge densities $q_{x}$. The complete time propagation on the relevant subspace, spanned by $\{e_{x}\}$, is described by a cyclic block three-diagonal matrix

[TABLE]

where $a$, $b$, $c$ are $2\times 2$ blocks. The projected initial state (24) reads ${\underline{p}}=\bigoplus_{x=1}^{n}\mu_{x-n/2}$. Assuming that $n\geq 8\,t+2$, the charge profile (23) inside of the light cone no longer depends on $n$. Thus the limit $n\to\infty$ can be applied and an infinite matrix $\mathcal{U}$ can be diagonalized using the block Fourier transform, yielding two bands of eigenvalues $\lambda_{1,2}(k)$, $k\in[-\pi,\pi]$ and Bloch eigenvectors ${\underline{\omega}}_{1,2}(k)=\bigoplus_{x\in\mathbb{Z}}{\underline{v}}_{1,2}(k)e^{ikx}$ (see Sect. C of sup for details). Additionally note that $\langle q_{x}U^{-t}e_{y}\rangle=0$, if $|x-y|>2t$, so the charge density profile can finally be expressed as a Fourier integral

[TABLE]

where $\tilde{\alpha}_{1,2}(k)=\alpha_{1,2}(k)\left([1,1]\cdot\underline{v}_{1,2}(k)\right)$, and $\alpha_{1,2}(k)$ are the coefficients expressing the vector $[1,1]$ in terms of Bloch vectors $v_{1,2}(k)$ (details in Sect. C of sup ). We are interested in the behavior of $f(x,t)$ on the diffusive lines $2x=\xi\sqrt{t}$ in the large time limit $\tilde{f}(\xi)=\lim_{t\to\infty}f(\xi\sqrt{t}/2,t)$. In this limit, the contribution from the term proportional to $\lambda_{2}(k)^{t}$ can be disregarded, since ${\rm sup}_{k}|\lambda_{2}(k)|<1$. Additionally, one should note that $\lambda_{1}(0)=1$, so large $t$ asymptotics of (26) can be obtained by expanding $\lambda_{1}(k)$ around $k=0$, $\lambda_{1}(k)^{t}\simeq e^{-\gamma k^{2}t}$, $\gamma=\frac{1-\rho}{4\rho}$. Introducing a new variable $h=k\sqrt{t}$ the integral (26) can be calculated exactly in the scalling limit $x,t\to\infty$ with $x/\sqrt{t}$ fixed, yielding

[TABLE]

The agreement between the exact and numerical result can be seen in Fig 4. Since the solution of diffusion equation $\frac{\partial}{\partial t}f(x,t)=\mathcal{D}\frac{\partial^{2}}{\partial x^{2}}f(x,t)$ for a step initial data reads $f(x,t)\sim\text{erf}(\tfrac{x}{\sqrt{4\mathcal{D}t}})$, the diffusion constant is $\mathcal{D}=4\gamma=\rho^{-1}-1$.

Interestingly enough, our inhomogeneous quench does not excite any ballistic transport for the the specific class of initial states considered (22) even for $\mu_{L}+\mu_{R}\neq 0$, unlike the linear response result which has a finite Drude weight for $\mu\neq 0$. The reason is simple: inhomogeneous quench does not excite any imbalance of vacancy momentum, which is a conserved quantity. However, in the perturbative linear-potential quench derivation of linear response coefficients (Sect. A of sup ) one has a natural vacancy momentum bias which generates the Drude weight.

Discussion.– We have studied transport properties of a simple reversible and deterministic cellular automaton. Despite its simplicity the model exhibits a large variety of transport phenomena, including charge diffusion, and offers analytical handle on the calculation of transport coefficients as well as exactly solving interesting initial value problems. Algebraic structure of the model allows to identify microscopic mechanisms behind its various transport regimes. Specifically, the ballistic behavior of certain degrees of freedom can induce the diffusive transport of charge carriers, which are completely frozen in their absence.

Our results open many interesting questions. Firstly, it should be clarified whether the many body deterministic diffusion mechanism disclosed here applies to other integrable models, in particular to quantum lattice models such as $XXZ$. A promising idea in this direction is a formulation of quantum transport in terms of classical-like soliton gas Doyon et al. (2017b). Secondly, our exactly solvable model could serve as a testbed for a precise verification of predictions of GHD. And thirdly, one may imagine various solvable generalizations of our model. For example, we can define a stochastic Markov chain model by introducing a tunnelling probability $\Gamma$ for a particle exchange, i.e. to modify the following matrix elements of the local propagator $U_{(\pm,\mp),(\mp,\pm)}=\Gamma$, $U_{(\pm,\mp),(\pm,\mp)}=1-\Gamma$. Our analysis can be straightforwardly extended, resulting in the same value of the Drude weight as for the deterministic case $\Gamma=0$, while the Green-Kubo diffusion constant has a continuous dependence on $\Gamma$

[TABLE]

and diverges at $\Gamma=1$. This implies that the physics of our model is robust against adding a small noise or dissipation. We conjecture that one could similarly study quantization of the model and demonstrate continuity of the classical limit.

We acknowledge fruitful discussions with Bruno Bertini, Spyros Sotiriadis, Enej Ilievski, Lenart Zadnik and Marko Žnidarič. The work has been supported by grants P1-0044 and N1–0025 of Slovenian Research Agency, and ERC grant OMNES.

I A: Linear response coefficients.

In the first part we revisit the initial value problem, which yields the linear response conductivity and Drude weight. The protocol starts by applying a kick (or a quench) with a linear-potential field $h$

[TABLE]

to the initial translationally invariant equilibrium state $p$

[TABLE]

which induces a constant gradient of charge $q_{x}$, so that $\langle q_{n/2+x}\rangle_{\tilde{p}}=2\mu+\chi hx$, where $\chi$ is the static susceptibility. After the initial quench the system is left to evolve autonomously. The time dependent current, averaged over two neighboring temporal layers, now reads

[TABLE]

Taking into account the continuity equation (5) and discarding boundary terms, we can alternatively write the above expression as

[TABLE]

The last term vanishes as a consequence of a more general property, stating that the current in the bulk is an odd function of time $t^{\prime}$

[TABLE]

which follows from $U^{-{t^{\prime}}}=S^{-1}U^{t^{\prime}}S$. Finally, by noting $\langle{j_{n/2}^{1/2}J^{t^{\prime}+1/2}}\rangle_{p}=\langle{j_{n/2}J^{t^{\prime}}}\rangle_{p}$ and $\langle j_{n/2}J^{t^{\prime}}\rangle_{p}=\langle j_{n/2}^{-t^{\prime}}J\rangle_{p}$, we obtain

[TABLE]

Linear response conductivity $\sigma$ in the main text (13) follows directly from the definition $\langle j\rangle=\sigma h$, after taking first the limit $n\to\infty$ and then the limit $t\to\infty$, while the Drude weight can be interpreted as the rate at which the current in the system increases with time.

Note that this quench (30) at $\mu\neq 0$ also excites the gradient of density of vacancies, which can be understood by explicitly expressing $q_{x}\,p$,

[TABLE]

Here $p_{x}$ are local distributions $p_{x}={[0]}^{\prime}+\mu{[1]}^{\prime}$, as defined in the main text (14). Due to momentum conservation of vacancies, this implies ballistic transport and finite Drude weight unlike the kink shaped inhomogeneous quench.

II B: Autocorrelation functions.

Here we present further details on the calculation of the autocorrelation function. We start by expressing the local two-site propagator in the basis $\{{[0]},{[1]},{[2]}\}$,

[TABLE]

where $\alpha\in\{0,1,2\}$.

As argued bellow, the dynamics can be restricted to the subspace $\mathcal{A}_{J}$ spanned by the basis

[TABLE]

First of all, note that we can only consider the observables with a single ${[1]}$, forming the subspace $\mathcal{A}^{(1)}$. Let $P$ be a linear projector $P:\,\mathcal{A}^{(1)}\to\mathcal{A}_{J}$, $P^{2}=P$ and $Q=1-P$. We then observe that for every $a\in\mathcal{A}^{(1)}$, $\langle JQa\rangle_{p}=0$, and that $Q(\mathcal{A}^{(1)})$ is invariant under $SU^{\text{e}}$. Therefore the current autocorrelation function can be obtained by restricting the dynamics to $\mathcal{A}_{J}$, i.e. considering the reduced operator $\mathcal{U}$,

[TABLE]

To obtain the time propagation of the current density, the matrix $\mathcal{U}$ should be applied to the initial current density $\underline{J}$ expressed in term of basis elements $y_{d}$ and $z_{d}^{k}$

[TABLE]

Vector $\underline{o}$

[TABLE]

incorporates the overlaps between the basis elements $y_{d}$ and $z_{d}^{k}$ and the local current density w.r.t. the state $p$ (14). The correlation function now reads

[TABLE]

Splitting the overlap vector $\underline{o}$ into parts corresponding to $\mu=0$ and $\mu\neq 0$, $\underline{o}=\underline{o}_{D}+\mu^{2}\underline{o}_{B}$, the correlation function decomposes into the $\mu$ independent part and the contribution proportional to $\mu^{2}$. The $\mu$ independent part can be reduced to the $3\times 3$ block of the matrix $\mathcal{U}$ corresponding to the basis elements $y_{0},y_{1},z_{0}^{0}$, as expressed in the main text (17). The contribution from $\underline{o}_{B}$ can be obtained by applying $\mathcal{U}$ to the vector $\underline{o}_{B}$,

[TABLE]

which yields the expression for $C(t)$ from the main text (20).

III C: Inhomogeneous quench.

The complete expressions regarding inhomogeneous initial value problem are presented here. The local propagator takes the following form in the dual basis $\{{[0]}^{\prime},{[1]}^{\prime},{[2]}^{\prime}\}$,

[TABLE]

with $\alpha\in\{0,1,2\}$. Using the local propagator, we can explicitly propagate basis states $e_{x}=[\ldots 010\ldots]^{\prime}$ to obtain the $2\times 2$ blocks of the cyclic block tridiagonal matrix $\mathcal{U}$ (25),

[TABLE]

The solution of the eigenvalue equation $\mathcal{U}{\underline{\omega}}=\lambda{\underline{\omega}}$, within the Bloch ansatz ${\underline{\omega}}=\bigoplus_{r}\exp({\rm i}kr)\underline{v}(k)$ reading explicitly

[TABLE]

yields two bands and the corresponding eigenvectors

[TABLE]

We wish to express the vector $\underline{e}_{2y-1}+\underline{e}_{2y}$, occurring in the expression for $\underline{p}$ (${\underline{e}}_{y}$ is a unit infinite dimensional vector in with a single $1$ at place $y\in\mathbb{Z}$ and zeros elsewhere) in terms of the eigenvectors $\underline{\omega}_{1,2}(k)=\bigoplus_{r}\exp({\rm i}kr)\underline{v}_{1,2}(k)$ by solving the following system of equations

[TABLE]

Using the inverse discrete Fourier transform we obtain

[TABLE]

Expressing $\underline{\omega}_{1,2}(k)$ in terms of local components $\underline{v}_{1,2}(k)$ yields a set of equations

[TABLE]

which admit the following solution:

[TABLE]

The vectors ${\underline{e}}_{2x-1}+{\underline{e}}_{2x}$ expressed in eigenbasis (48) are combined to form the initial state ${\underline{p}}$, which can then be propagated in time. As mentioned in the main text, for a given $x$ and $t$ in the expression for profile $f(x,t)$, only the maximal light cone contributions from ${\underline{p}}$ can be taken into account. To compactly represent the function $f(x,t)$ the time and space independent constants are joined together as $\tilde{\alpha}_{1,2}(k)=\alpha_{1,2}(k)\left([1,1]\cdot\underline{v}_{1,2}(k)\right)$. Inserting the expressions (51) into the main text equation (26), and making an approximation to the eigenvalues discussed in the main text, the following expression is obtained

[TABLE]

Identifying the above sum with the Gaussian integral in the limit $t\to\infty$, we obtain the final result (27).