Ballistic spin transport in a periodically driven integrable quantum system

TL;DR

This paper demonstrates ballistic spin transport in a periodically driven integrable quantum system, revealing a fractal lower bound on the spin Drude weight and providing insights into quantum transport in Floquet spin chains.

Contribution

It introduces an analytic family of quasi-local conservation laws for a Floquet spin chain and establishes a fractal lower bound on the spin Drude weight.

Findings

Ballistic spin transport observed in the model.

Fractal dependence of the lower bound on anisotropy.

Numerical simulations support the tightness of the bound.

Abstract

We demonstrate ballistic spin transport of an integrable unitary quantum circuit, which can be understood either as a paradigm of an integrable periodically driven (Floquet) spin chain, or as a Trotterized anisotropic ($XXZ$) Heisenberg spin-1/2 model. We construct an analytic family of quasi-local conservation laws that break the spin-reversal symmetry and compute a lower bound on the spin Drude weight which is found to be a fractal function of the anisotropy parameter. Extensive numerical simulations of spin transport suggest that this fractal lower bound is in fact tight.

Figures

Equations

$$D=\underset{t\to \infty }{\lim }\underset{N\to \infty }{\lim }\underset{\mu \to 0}{\lim }\frac{⟨J(t){⟩}_{\mu }}{2Nt\mu },$$

$$\begin{array}{rl}D=& \underset{t\to \infty }{\lim }\underset{N\to \infty }{\lim }\frac{1}{2N}\frac{1}{t}\sum _{\tau =1}^t⟨JJ(\tau )⟩\\ & \ge \underset{N\to \infty }{\lim }\frac{1}{2N}\sum _k\frac{|⟨J,Q_k⟩{|}^2}{⟨Q_k,Q_k⟩}.\end{array}$$

$$\begin{array}{r}{\mathcal{U}}_{odd}=\prod _{n=1}^{N/2}{U}_{2n,2n+1},{\mathcal{U}}_{even}=\prod _{n=1}^{N/2}{U}_{2n-1,2n},\end{array}$$

$$U_{n,n+1}=e^{-i\mathcal{J}({\sigma }_n^x{\sigma }_{n+1}^x+{\sigma }_n^y{\sigma }_{n+1}^y)-i{\mathcal{J}}^{\prime }({\sigma }_n^z{\sigma }_{n+1}^z-1)}$$

$$\check{R}(\lambda )=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& \frac{\sin \eta }{\sin (\lambda +\eta )}& \frac{\sin \lambda }{\sin (\lambda +\eta )}& 0\\ 0& \frac{\sin \lambda }{\sin (\lambda +\eta )}& \frac{\sin \eta }{\sin (\lambda +\eta )}& 0\\ 0& 0& 0& 1\end{array}\right)$$

$$e^{2i(\mathcal{J}\pm {\mathcal{J}}^{\prime })}=\frac{\sin \eta -\sin \lambda }{\sin (\eta \pm \lambda )}.$$

$$\displaystyle h_{n,n+1}=\frac{1}{2\sin\eta}\big{(}\sigma_{n}^{x}\sigma_{n+1}^{x}+\sigma_{n}^{y}\sigma_{n+1}^{y}+\Delta\,(\sigma_{n}^{z}\sigma_{n+1}^{z}-\mathbb{1})\big{)},$$

$$\begin{array}{rl}& {\mathcal{U}}^{†}{\sigma }_{2n+1}^z\mathcal{U}-{\sigma }_{2n+1}^z=-j_{2n+2}+j_{2n+1}^{\prime },\\ & {\mathcal{U}}^{†}{\sigma }_{2n}^z\mathcal{U}-{\sigma }_{2n}^z=-j_{2n+1}^{\prime }+j_{2n}\end{array}$$

$$\displaystyle\begin{gathered}j_{2n}=\frac{4\sin\lambda\sin\eta}{\cos 2\eta-\cos 2\lambda}\,\big{(}\sigma^{+}_{2n-1}\,\sigma^{-}_{2n}-\sigma^{-}_{2n-1}\,\sigma^{+}_{2n}\big{)}+\\ +\frac{2\,(\sin\lambda)^{2}}{\cos 2\eta-\cos 2\lambda}\,\big{(}\sigma^{z}_{2n-1}-\sigma^{z}_{2n}\big{)}\end{gathered}$$

$$\displaystyle j_{n}=j^{\prime}_{n}=-\frac{2\lambda}{\sin\eta}\big{(}\sigma^{+}_{n-1}\,\sigma^{-}_{n}-\sigma^{-}_{n-1}\,\sigma^{+}_{n}\big{)}+\mathcal{O}(\lambda^{2}),$$

$$\displaystyle T(\varphi,s)=\operatorname{tr}\!{{}_{a}}\Big{[}\prod_{n=1}^{N}\mathbf{L}_{n,a}\big{(}\varphi-(-1)^{n}\,{\textstyle\frac{\lambda}{2}},s\big{)}\Big{]}.$$

$$\mathbf{L}(\phi ,s)=\frac{1}{\sin \phi }\left(\begin{array}{cc}\sin (\phi +\eta {\mathbf{S}}_s^z)& \sin (\eta ){\mathbf{S}}_s^{-}\\ \sin (\eta ){\mathbf{S}}_s^{+}& \sin (\phi -\eta {\mathbf{S}}_s^z)\end{array}\right)$$

$$\begin{array}{rl}& {\mathbf{S}}_s^z=\sum _{k=0}^{m-1}(s-k)|k⟩⟨k|,\\ & {\mathbf{S}}_s^{+}=\sum _{k=0}^{m-2}\frac{\sin (k+1)\eta }{\sin \eta }|k⟩⟨k+1|,\\ & {\mathbf{S}}_s^{-}=\sum _{k=0}^{m-2}\frac{\sin (2s-k)\eta }{\sin \eta }|k+1⟩⟨k|.\end{array}$$

$$Z(\phi )=\frac{1}{2\eta \sin \eta }{\partial }_{s}T(\phi ,s){\mid }_{s=0},$$

$$j(\phi )=\underset{N\to \infty }{\lim }\frac{1}{N}⟨Z(\bar{\phi }),J⟩=\frac{\sin \lambda }{(\cos \lambda -\cos 2\phi )\sin \eta }.$$

$$\displaystyle\begin{gathered}K(\varphi,\varphi^{\prime})=\frac{\big{(}\cos(\varphi-\varphi^{\prime}+\lambda)+\cos(\varphi-\varphi^{\prime}-\lambda)-2\cos(\varphi+\varphi^{\prime})\big{)}\,\sin[(m-1)(\varphi+\varphi^{\prime})]+(\sin\lambda)^{2}\,\sin[m(\varphi+\varphi^{\prime})]}{4\,(\sin\eta)^{2}\,(\cos 2\varphi-\cos\lambda)\,(\cos\lambda-\cos 2\varphi^{\prime})\,\sin[m(\varphi+\varphi^{\prime})]}.\end{gathered}$$

$$D\ge D_{\mathrm{Mazur}}=\frac{1}{2}\mathrm{Re}\int d{}^2\phi \overline{j(\bar{\phi })}f(\phi ),$$

$$\int d{}^2{\phi }^{\prime }K(\phi ,{\phi }^{\prime })f({\phi }^{\prime })=j(\phi ).$$

$$D_{\mathrm{Mazur}}^{\prime }={\left[\frac{\sin \eta }{\sin (\pi /m)}\right]}^2\left(1-\frac{m}{2\pi }\sin (2\pi /m)\right),$$

$$\underset{m\to \infty }{\lim }{D}_{\mathrm{Mazur}}=2\left(1-\frac{\mathrm{Gd}(|\lambda |)}{\sinh (|\lambda |)}\right),$$

$$D_{\mathrm{Mazur}}=2(1-\mathrm{sech}(|\lambda |)).$$

$${\rho }_{\mu }=\frac{1}{2^N}\left(1-\mu \sum _{n=1}^N\left(n-\frac{N+1}{2}\right){\sigma }_n^z\right).$$

$$\mathcal{U}{\sigma }_n^z{\mathcal{U}}^{†}-{\sigma }_n^z=\mathcal{U}(j_{n+1}-j_n){\mathcal{U}}^{†},$$

$$\begin{array}{c}⟨A(t){⟩}_{\mu }=⟨A⟩-\frac{\mu }{2^N}\sum _{n=1}^N\left(n-\frac{N+1}{2}\right)\mathrm{tr}[{\mathcal{U}}^t{\sigma }_n^z{\mathcal{U}}^{-t}A]=⟨A⟩-\frac{\mu }{2^N}\sum _{n=1}^N\left(n-\frac{N+1}{2}\right)\mathrm{tr}[{\mathcal{U}}^{t-1}{\sigma }_n^z{\mathcal{U}}^{-t+1}A]-\\ -\frac{\mu }{2^N}\sum _{n=1}^N\left(n-\frac{N+1}{2}\right)\mathrm{tr}[{\mathcal{U}}^t(j_{n+1}-j_n){\mathcal{U}}^{-t}A]=⟨A(t-1){⟩}_{\mu }+\frac{\mu }{2^N}\sum _{n=1}^N\mathrm{tr}[j_n{\mathcal{U}}^{-t}A{\mathcal{U}}^t]=⟨A(t-1){⟩}_{\mu }+\mu ⟨JA(t)⟩.\end{array}$$

$$⟨J(t){⟩}_{\mu }=⟨J{⟩}_{\mu }+\mu \sum _{\tau =1}^t⟨JJ(\tau )⟩$$

$$D\equiv \underset{t\to \infty }{\lim }\underset{N\to \infty }{\lim }\underset{\mu \to 0}{\lim }\frac{⟨J(t){⟩}_{\mu }}{2Nt\mu }=\underset{t\to \infty }{\lim }\underset{N\to \infty }{\lim }\frac{1}{2N}\frac{1}{t}\sum _{\tau =1}^t⟨JJ(\tau )⟩.$$

$$\underset{t\to \infty }{\lim }\underset{N\to \infty }{\lim }\underset{\mu \to 0}{\lim }\frac{⟨j_n(t){⟩}_{\mu }}{2t\mu }=\underset{t\to \infty }{\lim }\underset{N\to \infty }{\lim }\frac{1}{2t}\sum _{\tau =1}^t⟨J{j}_n(\tau )⟩.$$

$$D=\underset{t\to \infty }{\lim }\underset{N\to \infty }{\lim }\underset{\mu \to 0}{\lim }\frac{1}{2t\mu }⟨\frac{j_{2n}(t)+j_{2n-1}^{\prime }(t)}{2}{⟩}_{\mu },$$

$$R(\lambda )=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& \frac{\sin \lambda }{\sin (\lambda +\eta )}& \frac{\sin \eta }{\sin (\lambda +\eta )}& 0\\ 0& \frac{\sin \eta }{\sin (\lambda +\eta )}& \frac{\sin \lambda }{\sin (\lambda +\eta )}& 0\\ 0& 0& 0& 1\end{array}\right)$$

$$\displaystyle T(\varphi)=\operatorname{tr}\!{{}_{0}}\Big{[}\prod_{n=1}^{N}R_{n,0}(\varphi-(-1)^{n}\,{\textstyle\frac{\lambda}{2}})\Big{]}.$$

Full text

Introduction.– Understanding transport in various out-of-equilibrium setups, in particular in low dimensions, is one of the main challenges of theoretical condensed matter physics affleck . Experimental evidence corroborates the seemingly controversial proposal castella ; zotos that integrable systems generically exhibit ballistic transport even at high temperatures experimental1 ; experimental2 ; FHM1 . This proposal has received rigorous justification in terms of the existence of an extensive number of quasi-local conserved quantities ip13 ; tomaz11 ; tomaz13 ; tomaz14 ; pereira14 ; review which form a basis for the hydrodynamic theory of interacting integrable systems bertini16 ; doyon16 ; deluca17 .

Recently, periodically driven (Floquet) spin chains with local interactions have attracted considerable attention. This was in particular due to the possibility of exhibiting generalized thermalization towards non-equilibrium steady states ggefloquet and distinct dynamical phases with respect to the spontaneous breaking of time-translation invariance nayak ; sondhi . Still, the possibility of strictly ballistic transport in interacting quantum integrable Floquet systems has never been explored (see katja ; sarang for a classical lattice setting), even though a peculiar robustness of transport to integrability breaking has been observed a while ago tp98 ; tp98ql (also arnab ).

For concreteness, let us consider spin transport. Without resorting to the spectroscopic approach FHM2 , which is harder to justify in Floquet systems, ballistic transport can be defined as a linear growth of the spin current in time, after the system has been prepared in an initial state supporting a small gradient of magnetization. This can be formulated in terms of a nonzero Drude weight

[TABLE]

where $J=\sum_{n}j_{n}$ is the extensive spin current operator on a spin chain of length $N$ and $J(t)$ its time dependence. $\langle\bullet\rangle_{\mu}$ denotes the average in the initial state with a small gradient of magnetization $\mu$, say $\rho_{\mu}\sim\exp\left(\mu\sum_{n}n\,\sigma^{z}_{n}\right)$, where $\sigma^{z}_{n}$ is a local spin variable. A formula similar to (1) holds even if the system is initially prepared in two equilibrated halves at different magnetizations $\mu_{L}$ and $\mu_{R}$ with $\mu\sim(\mu_{L}-\mu_{R})/N$ representing the effective gradient vasseur15 ; karrasch17 ; ilievski17 . This partitioned initial state is easier to simulate using state-of-the-art tensor network simulations.

Expanding to the first order in $\mu$, the Drude weight can be expressed solely in terms of equilibrium auto-correlation functions using the Kubo formula, see Appendix A of the Supplemental material (SM) supp . This can in turn be bounded from below by means of the Mazur inequality M69 ; suzuki ; zotos (see ip13 for a rigorous derivation in extended systems)

[TABLE]

Here $Q_{k}$ are conserved quantities orthogonal with respect to the inner product $\langle A,B\rangle=\operatorname{tr}\!{[}A^{\dagger}B]/2^{N}$, assuming that the reference equilibrium state is the maximum entropy state $\rho_{\mu=0}=2^{-N}\mathbb{1}$ and the local Hilbert space dimension is $2$. In order for the bound to be finite the conserved quantities should be linearly extensive or quasi-local, $\langle Q_{k},Q_{k}\rangle\propto N$, and should have a finite overlap with the spin current, $\langle J,Q_{k}\rangle\neq 0$. For the latter to hold, $Q_{k}$ must not be symmetric, $\mathcal{P}Q_{k}\mathcal{P}\neq Q_{k}$, with respect to the spin-reversal operator $\mathcal{P}$, which flips the spin current ${\cal P}J{\cal P}=-J$. In the easy-plane (gapless) regime of the Heisenberg $XXZ$ model, for example, the integrals of motion with all of the required properties have been shown to exist tomaz11 ; tomaz13 ; tomaz14 ; pereira14 .

In the present paper, we aim to rigorously establish a regime of ballistic transport in a Floquet driven integrable model related to the Trotterized $XXZ$ spin-$1/2$ chain. We introduce the dynamical protocol as a local quantum circuit, establish its connection with the six-vertex $R$-matrix and integrability structure of the $XXZ$ model, and define the spin currents and continuity equations arising from the global $U(1)$ symmetry of the model. Despite its driven nature, we construct a set of quasi-local conservation laws which break the spin-reversal symmetry. We then show how to evaluate the optimized Mazur lower bound on the spin Drude weight. Extensive numerical simulations using the time-evolving block decimation (TEBD) algorithm strongly suggest that this bound, which is a fractal function of parameters, is in fact saturated, similarly as in the continuous-time case ilievski17 ; klumper19 .

The model.– Consider a spin-1/2 chain with $N\in 2\mathbb{Z}$ sites and periodic boundary conditions. The local physical space on each site will be denoted by $\mathcal{V}_{p}\equiv\mathbb{C}^{2}$. We are interested in a discrete-time Liouville-von Neumann equation for a density matrix $\rho_{t+1}=\mathcal{U}\,\rho_{t}\,\mathcal{U}^{\dagger}$.

The propagator $\mathcal{U}=\mathcal{U}_{odd}\,\mathcal{U}_{even}$ acts in two steps

[TABLE]

where

[TABLE]

is a unitary gate acting on two neighbouring sites labelled with $n$ and $n+1$, see Figure 1. Here $\sigma^{\alpha}$ ($\alpha={x},{y},{z}$) are Pauli matrices. By considering infinitesimally small couplings $\mathcal{J}$ and $\mathcal{J}^{\prime}$ and an infinite number of time steps we recover the continuous-time dynamics of the $XXZ$ model, according to the Trotter-Suzuki formula.

The local 2-site unitary gate can be rewritten as $U_{n,n+1}=\check{R}_{n,n+1}(\lambda)$ where

[TABLE]

denotes the braid form of the $R$-matrix of the $XXZ$ model. The new parameters $\eta$ and $\lambda$ can be implicitly expressed as unique functions of $\mathcal{J}$ and $\mathcal{J}^{\prime}$ through the following pair of relations

[TABLE]

The continuous-time limit is recovered as an expansion in small $\lambda$ which gives $U_{n,n+1}=\mathbb{1}+\lambda\,h_{n,n+1}+\mathcal{O}(\lambda^{2})$ with the local Hamiltonian density

[TABLE]

with $\Delta=\cos\eta$ being the anisotropy parameter. Clearly, real $\eta$ and imaginary $\lambda$ correspond to the gapless/easy-plane regime, shown in Figure 2, whereas imaginary $\eta$ and real $\lambda$ correspond to the gapped/easy-axis regime.

In the gapless regime of the continuous-time limit ($|\Delta|<1$) the Drude weight has rigorously been shown to be nonzero for a dense set of anisotropies parametrized by $\eta=l\,\pi/m$, where $l$ and $m$ are coprime integers tomaz13 ; tomaz14 . In this paper we extend this discussion for the same set of anisotropies to a discrete time, i.e., to all imaginary $\lambda$. This will cover the ballistic regime in the phase diagram of our model, shown in red in Figure 3, which is determined by adapting the numerical method of Ref. ljubotina . We stress that we observe ballistic transport for any ratio ${\cal J}^{\prime}/{\cal J}$, even for $|\mathcal{J}^{\prime}|>|\mathcal{J}|$, unlike in the continuous-time case. Note also that the other two transport regimes can be clearly established numerically – the super-diffusive in yellow and the diffusive in blue.

Spin currents.– Due to the $U(1)$ symmetry of the propagator (5) the total magnetization $M=\sum_{n=1}^{N}\sigma^{z}_{n}$ is a conserved quantity. As a result of the discrete time propagation we identify two continuity equations, separately for odd/even sites

[TABLE]

Through them we can define two local current densities, the even current $j_{2n}$ and the odd, $j^{\prime}_{2n+1}$. The even one is

[TABLE]

while the odd current can be computed as $j_{2n+1}^{\prime}=\,\mathcal{U}_{even}^{\dagger}\,j_{2n+1}\,\mathcal{U}_{even}$ and operates on four adjacent sites. In the continuous-time limit both local currents reduce to

[TABLE]

with the pre-factor $(\sin\eta)^{-1}$ coming from the Hamiltonian (7).

The total extensive spin current is now defined as $J=\sum_{n=1}^{N/2}(j_{2n}+j^{\prime}_{2n+1})$. It is clearly anti-symmetric under spin-reversal $\mathcal{P}=\prod_{n=1}^{N}\sigma_{n}^{x}$, i.e., $\mathcal{P}J\mathcal{P}=-J$. We now proceed to construct the relevant conservation laws for all imaginary $\lambda$ and for a dense set of commensurate values of the anisotropy parameter $\eta=l\,\pi/m$.

Quasilocal integrals of motion.– The construction of anti-symmetric conservation laws is similar as in the continuous-time limit tomaz14 . However, they are now generated by the staggered transfer operator

[TABLE]

Here, $\mathbf{L}_{n,a}(\varphi,s)$ denotes the Lax operator acting on the $n$-th physical space $\mathcal{V}_{p}$ in the string $\bigotimes_{n=1}^{N}\mathcal{V}_{p}$ as a $2\times 2$ matrix

[TABLE]

whose elements are themselves matrices in the auxiliary space $\mathcal{V}_{a}$. For $\eta=l\,\pi/m$ the latter is an $m$-dimensional complex spin-$s$ representation of the quantum group $\mathcal{U}_{q}(sl_{2})$ ($q=e^{i\eta}$) traced out in the final expression (13). Its generators have an explicit form reminiscent of the angular momentum generators

[TABLE]

Together with $\check{R}(\lambda)$ given in (5), the Lax operator (14) satisfies the Yang-Baxter equation which implies $[T(\varphi,s),\mathcal{U}]=0$ and $[T(\varphi,s),T(\varphi^{\prime},s)]=0$ (see Appendices B and C of the SM supp ).

For $\lambda=0$ the spin-reversal asymmetric conservation laws of the $XXZ$ spin $1/2$ chain were previously produced tomaz13 ; tomaz14 ; pereira14 as

[TABLE]

and shown to be linearly extensive (quasi-local) inside an analyticity strip $|{\rm Re}\,\varphi-\frac{\pi}{2}|<\frac{\pi}{2m}$. Here we will simply show that this expression can be extended to arbitrary values of parameter $\lambda$ if the staggered form (13) of the transfer matrix is used. Since $\lambda$ is purely imaginary, the region of quasi-locality remains the same. The detailed construction of these conservation laws for finite $\lambda$ is presented in Appendix D of the SM supp .

In order to maximize the Mazur lower bound (2) for the spin Drude weight we need to minimize the norm of the conservation laws $Q_{k}\sim Z(\varphi)$ without reducing the overlap with the spin current $\langle J,Q_{k}\rangle$. Due to the asymmetry of the current operator $\mathcal{P}J\mathcal{P}=-J$, only the spin-reversal anti-symmetric component $Z^{-}(\varphi)=\frac{1}{2}(Z(\varphi)-\mathcal{P}Z(\varphi)\mathcal{P})$ contributes to the lower bound. Furthermore, we can subtract a term proportional to the total magnetization, $Z^{-}_{\perp}(\varphi)=Z^{-}(\varphi)-\frac{1}{N}\langle M,Z(\varphi)\rangle\,M$, since the latter is orthogonal to the spin current, see Appendix E of the SM supp . The overlap between $Z^{-}_{\perp}(\varphi)$ and the current is now given by

[TABLE]

The quasi-locality of $Z^{-}_{\perp}(\varphi)$ is seen from the $N$-independence of $K(\varphi,\varphi^{\prime})\equiv\lim_{N\to\infty}\frac{1}{N}\langle Z_{\perp}^{-}(\bar{\varphi}),Z_{\perp}^{-}(\varphi^{\prime})\rangle$ with the following conjectured analytical result

[TABLE]

Expressions $j(\varphi)$ and $K(\varphi,\varphi^{\prime})$, given by (17) and (19) respectively, are the essential ingredients for the lower bound on the Drude weight, which we discuss next.

Mazur bound.– We now attempt to bound the spin Drude weight by means of the Mazur inequality as elaborated on in Ref. tomaz14 . It can be rewritten in an integral form

[TABLE]

where $f(\varphi)$ solves the following Fredholm equation

[TABLE]

The integrals are formally taken over the area of quasi-locality $|{\rm Re}\,\varphi-\frac{\pi}{2}|<\frac{\pi}{2m}$. However, due to holomorphicity, a single line of integration centred at ${\rm Re}\,\varphi=\frac{\pi}{2}$ is sufficient, which makes for an efficient quasi-exact numerical procedure of computing $D_{\rm Mazur}$. The full Mazur lower bound has a fractal dependence on $\eta$ and a continuous dependence on $|\lambda|$. It has been calculated numerically and compared to TEBD vidal03 ; schollwock11 simulations – see Appendix H of the SM supp . The dependence on $\eta$ is shown in Figure 4, and the dependence on $|\lambda|$ in Figure 5.

In Figure 4 we have re-scaled the Drude weight and the lower bound by a factor of $(\sin\eta)^{2}$ [see Eq. (7)]. This allows us to make a comparison with the established continuous-time result tomaz14

[TABLE]

found by expanding $(\sin\eta)^{2}\,D_{\rm Mazur}=\lambda^{2}\,D^{\prime}_{\rm Mazur}+\mathcal{O}(\lambda^{3})$ around $\lambda=0$. We can see this by noting, that the small-$\lambda$ expansion of Eqs. (20) and (21) reproduces the corresponding equations in the continuous-time case.

The integral equations can be solved analytically for $m\to\infty$, corresponding to an irrational value of $\eta/\pi$. The details are discussed in Appendix G of the SM supp and the result is the enveloping function

[TABLE]

where ${\rm Gd}$ denotes the Gudermannian ${\rm Gd}(x)=2\arctan(e^{x})-\pi/2$. This represents a continuous strict lower bound on top of which an additional fractal structure, shown in Figure 4 for $|\lambda|=1$, emerges.

On the free fermion line ($\mathcal{J}^{\prime}=0$ or $m=2$) exact diagonalization shows the saturation of the lower bound which includes only a single conserved quantity $Z(\pi/2)$. It can easily be computed to give

[TABLE]

For a complete $|\lambda|$ dependence of the Mazur lower bound see Figure 5. Note that the $|\lambda|\to\infty$ limit is always $2$. This can easily be explained, since there, the local propagator reduces to a SWAP gate. As such, transport becomes perfectly ballistic with no scattering at all.

Discussion.– We have demonstrated and proven ballistic transport in a periodically driven interacting quantum spin chain, namely in the Trotterized $XXZ$ spin-1/2 model. We have used the notion of ballistic spin transport referring to a linearly growing extensive spin current after a quench from an inhomogeneous initial state with either a linear gradient or a step bias in the magnetization profile. We argue that this is the most natural definition of ballistic transport in the case of discrete-time propagation. Using the quasi-local conservation laws that we constructed by means of quantum-group theoretic methods, we have calculated the lower bound on the spin Drude weight and explicitly showed its fractal dependence on the anisotropy parameter. Extensive numerical simulations suggest the saturation of the lower bound – see Figure S-2 in Appendix H of the SM. Note, however, that for a fixed commensurate anisotropy $\eta=l\,\pi/m$ the convergence with time seems to become extremely slow with increasing $m$, certainly beyond ultimate verification with state-of-the art numerical methods.

In the continuous-time limit we correctly reproduce the well established results of ballistic spin transport in the $XXZ$ model. However, since the thermodynamic Bethe ansatz has not yet been developed for driven integrable systems vladimir , our results open interesting new avenues for research. The conservation laws that we proposed (see also Ref. zadnik ) can be directly applied for construction of complete generalized Gibbs ensembles and development of generalized hydrodynamics in integrable Floquet systems.

The authors thank M. Žnidarič for helpful discussions on the topic. The authors acknowledge support by the European Research Council (ERC) through the advanced grant 694544 – OMNES and the grant P1-0402 of Slovenian Research Agency (ARRS).

Appendix A Appendix A: Drude weight and Kubo formula

In this appendix we derive the Kubo form of the Drude weight starting from an initial state with a linear gradient of magnetization

[TABLE]

Here $\mu$ is a perturbative parameter which corresponds to a local gradient in magnetization. Our goal is to examine the time dependent average $\langle J(t)\rangle_{\mu}=\operatorname{tr}\!{[\rho_{\mu}\,J(t)]}$ of the extensive current $J=\sum_{n}j_{n}$ for small $\mu$, where $J(t)=\mathcal{U}^{-t}\,J\,\mathcal{U}^{t}$. Note, that in our case the local current densities on even and odd sites are different. In the main text (MT) they are denoted by $j_{2n}$ and $j^{\prime}_{2n-1}$, respectively. Using the continuity equation $\mathcal{U}^{\dagger}\,\sigma_{n}^{z}\,\mathcal{U}-\sigma_{n}^{z}=-j_{n+1}+j_{n}$ (Eq. (8) of the MT), which can be recast as

[TABLE]

we see, for an arbitrary observable $A$,

[TABLE]

Using this recurrence we get

[TABLE]

for the total extensive spin current $J$. Now, our definition of the Drude weight (Eq. (1) of the MT) gives the Kubo formula (Eq. (2) of the MT)

[TABLE]

Taking the local current $j_{n}$ instead of $J$ in (S-6) we see

[TABLE]

Combining (S-7) and (S-8) for even and odd $n$ and using the invariance of the equilibrium average $\langle\bullet\rangle$ under the translation for two sites we now realize

[TABLE]

which shows that it is enough to measure the weight locally at any pair of neighbouring sites.

Appendix B Appendix B: Integrability of the driven model

In this appendix we formulate the integrable structure of our driven model. Let $P$ denote the permutation operator which acts on the tensor square of the local physical space, $\mathcal{V}_{p}\otimes\mathcal{V}_{p}$, and has the explicit form P_{n,n+1}=\frac{1}{2}\big{(}\mathbb{1}+\boldsymbol{\sigma}_{n}\cdot\boldsymbol{\sigma}_{n+1}\big{)}, where $\boldsymbol{\sigma}=(\sigma^{x},\sigma^{y},\sigma^{z})$. Our local propagator, given by Eq. (4) of the MT, is a unitary matrix which can be written as $U=\check{R}(\lambda)=PR(\lambda)$ (see also Eq. (5) of the MT), where

[TABLE]

satisfies the Yang-Baxter equation $R_{1,2}(\lambda-\mu)R_{1,3}(\lambda)R_{2,3}(\mu)=R_{2,3}(\mu)R_{1,3}(\lambda)R_{1,2}(\lambda-\mu)$ defined on the triple $\mathcal{V}_{p}\otimes\mathcal{V}_{p}\otimes\mathcal{V}_{p}$. This implies the commutation of transfer matrices $[T(\varphi),T(\varphi^{\prime})]=0$,

[TABLE]

Index [math] in (S-11) denotes a copy of the physical vector space $\mathcal{V}_{p}$, which is traced out by the partial trace $\operatorname{tr}\!{{}_{0}}$. Elementary calculation yields

[TABLE]

with $\mathcal{U}_{odd}$ and $\mathcal{U}_{even}$ defined in Eq. (3) of the MT. This establishes integrability, since $T(\varphi)$ generates local integrals of motion through logarithmic derivatives. For $\lambda=0$, the transfer matrix (S-11) reduces to the one of the Heisenberg $XXZ$ model with the local Hamiltonian density

[TABLE]

Appendix C Appendix C: Transfer matrix from the complex spin representation

In this appendix we introduce the Lax operator and the transfer matrix which produces the spin-reversal asymmetric quasi-local integrals of motion $Z(\varphi)$. Such integrals exist only for $\eta=l\,\pi/m$, where $l$ and $m$ are co-prime integers. In our model $\lambda$ should then be purely imaginary. The starting point is the Lax operator of the $XXZ$ model,

[TABLE]

which acts on the pair $\mathcal{V}_{p}\otimes\mathcal{V}_{a}$. $\mathcal{V}_{a}$ is an $m$-dimensional complex spin-$s$ representation of the quantum group $\mathcal{U}_{q}(sl_{2})$, where $q=e^{i\eta}$. Its generators satisfy $[\mathbf{S}_{s}^{+},\mathbf{S}_{s}^{-}]=\sin[2\eta\,\mathbf{S}_{s}^{z}]/\sin\eta$, $[\mathbf{S}_{s}^{z},\mathbf{S}_{s}^{\pm}]=\pm\mathbf{S}_{s}^{\pm}$ and can take the following form

[TABLE]

Along with the $R$-matrix (S-10), the Lax operator (S-14) satisfies another Yang-Baxter equation

[TABLE]

on the triple $\mathcal{V}_{p}\otimes\mathcal{V}_{p}\otimes\mathcal{V}_{a}$. As in Appendix B, we can use the Lax operator to construct a family of transfer matrices

[TABLE]

Yang-Baxter equations imply $[T(\varphi,s),T(\varphi^{\prime},s)]=0$ and $[T(\varphi,s),T(\varphi^{\prime})]=0$. As we will see, the spin-reversal asymmetric integrals of motion are

[TABLE]

Appendix D Appendix D: Explicit form of the integrals of motion

Here we describe the explicit form of the conservation laws $Z(\varphi)$. Let us introduce the components of the Lax operator and its derivative at $s=0$ through

[TABLE]

where $\sigma^{0}=\mathbb{1}$ denotes the identity. Components $\mathbf{L}^{\alpha}$, $\widetilde{\mathbf{L}}^{\alpha}$ act on the auxiliary space $\mathcal{V}_{a}$ and explicitly read

[TABLE]

We should now make the following observations:

Of the derivative components, only $\widetilde{\mathbf{L}}^{z}(\varphi)$ preserves the highest-weight state $\ket{0}$. 2. 2.

Of the Lax components, only $\mathbf{L}^{0}(\varphi)$ preserves the highest-weight state $\ket{0}$. 3. 3.

No Lax component $\mathbf{L}^{\alpha}(\varphi)$, where $\alpha\in\{0,+,-,{z}\}$ can lift us from the highest-weight state $\ket{0}$ into the reduced subspace $\mathcal{V}_{a}^{\prime}={\rm lsp}\{\ket{1},\ket{2},...,\ket{m-1}\}$.

Because of the partial trace over the auxiliary space in (S-17), the derivative with respect to $s$ in (S-18) can always be shifted to the rightmost position in the string of Lax operators. The trace can then be separated as $\operatorname{tr}\!{{}_{a}}(\bullet)=\bra{0}\bullet\ket{0}+\sum_{k=1}^{m-1}\bra{k}\bullet\ket{k}$. As we will see, the first term will produce a linearly extensive contribution and the second term a remainder which is exponentially small in the system’s size.

Suppressing the $\varphi$-dependence, let us first focus on the string $\bra{0}\mathbf{L}^{\alpha_{1}}\mathbf{L}^{\alpha_{2}}...\widetilde{\mathbf{L}}^{\alpha_{N}}\ket{0}$. Due to the first two observations, if $\alpha_{N}={z}$ all $\alpha_{n}$ for $n<N$ must be zero. This produces the following magnetization-like terms

[TABLE]

If $\alpha_{N}=+$, the only option is to have $\alpha_{L}=-$ for some fixed $L<N$ and $\alpha_{n}=0$ for all $n<L$. This produces the following terms

[TABLE]

where $S^{n}$ is an automorphism of operator algebra implementing a periodic shift for $n$ lattice sites, defined by $S^{n}(\sigma_{L}^{\alpha})=\sigma^{\alpha}_{L+n}$, and

[TABLE]

The rest of the trace, $\sum_{k=1}^{m-1}\bra{k}\mathbf{L}^{\alpha_{1}}\mathbf{L}^{\alpha_{2}}...\widetilde{\mathbf{L}}^{\alpha_{N}}\ket{k}$, will produce the remainder. Due to the third observation it can be rewritten as

[TABLE]

where $\operatorname{tr}\!{{}^{\prime}_{a}}$ denotes the partial trace over the reduced subspace $\mathcal{V}_{a}^{\prime}={\rm lsp}\{\ket{1},\ket{2},...,\ket{m-1}\}$ onto which the Lax components have been projected by $\sum_{k=1}^{m-1}\ket{k}\bra{k}_{a}$. Gathering all of the terms with explicit matrix product forms provided by equations (S-24) and (S-25) we have

[TABLE]

By $Z_{\infty}(\varphi)$ we have denoted the part of $Z(\varphi)$ that survives the thermodynamic limit $N\to\infty$, as we will see in Appendix F.

Appendix E Appendix E: Spin-reversal symmetric and anti-symmetric component

In this appendix we define the decomposition of the integrals of motion $Z(\varphi)$ with respect to the symmetry under the spin-reversal $\mathcal{P}=\prod_{n=1}^{N}\sigma_{n}^{x}$. Let us assume $\lim_{N\to 0}p_{N}(\varphi)=0$ in (S-29). The spin-reversal symmetric and anti-symmetric components of $Z(\varphi)$ are defined as

[TABLE]

Note, that the Lax operator (S-14) satisfies $\sigma_{n}^{x}\mathbf{L}_{n,a}^{\alpha}(\varphi,s)\sigma_{n}^{x}=\mathbf{L}_{n,a}^{\alpha}(\pi-\varphi,s)^{T}$, where $(\bullet)^{T}$ denotes the partial transposition with respect to the physical space. This immediately implies

[TABLE]

Since the operator $Z(\varphi)$ is not symmetric [see Eq. (S-24)] the spin-reversal anti-symmetric component $Z^{-}(\varphi)$ is nonzero. This is crucial, since it provides the overlap with the spin current.

We can now make these components orthogonal with respect to the total magnetization $M=\sum_{n=1}^{N}\sigma_{n}^{z}$ which does not contribute to the overlap with the spin current. We have

[TABLE]

with

[TABLE]

Appendix F Appendix F: Quasi-locality and the inner product between the integrals of motion

This appendix describes the quasi-locality of $Z(\varphi)$ given in (S-29) and provides explicit formulas for their overlaps. We will refer to operators $Z(\varphi)$ as quasi-local if $\langle Z(\varphi),Z(\varphi)\rangle\propto N$ for large $N$, where $\langle A,B\rangle=\operatorname{tr}\!{(A^{\dagger}B)}/2^{N}$ denotes the Hilbert-Schmidt inner product tomaz14 . We will see that the following holds:

Proposition 1**.**

Let $\lambda$ be purely imaginary and $\eta=\frac{l\,\pi}{m}$ with $l$ and $m$ co-prime integers. Operators $Z(\varphi)$ in (S-29) with terms given by (S-24) and (S-25) are quasi-local if $|{\rm Re}\,\varphi-\frac{\pi}{2}|<\frac{\pi}{2m}$.

We start by recalling the decomposition $Z(\varphi)=Z_{\infty}(\varphi)+p_{N}(\varphi)$ from (S-29) and writing the inner product

[TABLE]

where $(\bar{\bullet})$ denotes the complex conjugation. The first term can be written as

[TABLE]

where

[TABLE]

Using the Cauchy-Schwarz inequality, the absolute value of the last three terms in (S-34) can be bounded by the Hilbert-Schmidt norms of $Z_{\infty}(\varphi)$ and $p_{N}(\varphi)$. Proposition 1 then holds if the following is true:

Lemma 1**.**

The Hilbert-Schmidt norm of the remainder $p_{N}(\varphi)$ is exponentially small in the system size $N$ and $K^{\prime}(\bar{\varphi},\varphi)$ is finite if $|{\rm Re}\,\varphi-\frac{\pi}{2}|<\frac{\pi}{2m}$.

Figure S-1 shows the exponential suppression of the norm $\langle p_{N}(\varphi),p_{N}(\varphi)\rangle$ of the remainder (left) and the difference $K_{(num)}^{\prime}(\bar{\varphi},\varphi)-K^{\prime}(\bar{\varphi},\varphi)$ between the numerically computed $K^{\prime}(\bar{\varphi},\varphi)$ (S-38) and our conjectured formula (S-43) which is finite in the thermodynamic limit (right).

In the following we describe how to calculate $K^{\prime}(\varphi,\varphi^{\prime})$ and conjecture its explicit formula, which we also use in the MT. To this end, let us define a transfer matrix $\mathbb{T}(\varphi,\varphi^{\prime})=\sum_{\alpha\in\{0,+,-,{z}\}}\mathbf{L}^{\alpha}(\varphi)\otimes\mathbf{L}^{\alpha}(\varphi^{\prime})\,\langle\sigma^{\alpha},\sigma^{\alpha}\rangle$ over ${\cal V}_{a}\otimes{\cal V}_{a}$. This enables us to write down the inner products of the local densities in a computationally convenient form

[TABLE]

Calculating the matrix $\mathbb{T}(\varphi,\varphi^{\prime})$ we see, that it

preserves the vector $\ket{0}\otimes\ket{0}$ and 2. 2.

preserves the subspace spanned by the vectors of the form $\ket{k}\otimes\ket{k}$, $k=0,1,...,m-1$.

Due to the first property, any excursion from the reduced auxiliary space $\mathcal{V}_{a}^{\prime}={\rm lsp}\{\ket{1},...,\ket{m-1}\}$ in the string of matrices $\mathbb{T}(\varphi,\varphi^{\prime})$ gives zero. This allows us to project these matrices onto the reduced auxiliary subspace $\mathcal{V}_{a}^{\prime}$ by identifying $\ket{k}\otimes\ket{k}\longleftrightarrow|\sin(k\eta)|\ket{k}$ and $\bra{k}\otimes\bra{k}\longleftrightarrow|\sin(k\eta)|^{-1}\bra{k}$. The matrix $\mathbb{T}(\varphi,\varphi^{\prime})$ becomes

[TABLE]

This matrix has two important properties:

For $|{\rm Re}\,\varphi-\frac{\pi}{2}|<\frac{\pi}{2m}$ its eigenvalues are strictly below $1$ in absolute value. This has been proven in tomaz14 . Obviously this also holds for $\mathbf{T}(\varphi\pm{\textstyle\frac{\lambda}{2}},\varphi^{\prime}\mp{\textstyle\frac{\lambda}{2}})$ if $\lambda$ is strictly imaginary. 2. 2.

Commutation $[\mathbf{T}(\varphi+{\textstyle\frac{\lambda}{2}},\varphi^{\prime}-{\textstyle\frac{\lambda}{2}}),\mathbf{T}(\varphi+{\textstyle\frac{\mu}{2}},\varphi^{\prime}-{\textstyle\frac{\mu}{2}})]=0$, which is actually a result of the Yang-Baxter equation (S-16) and its transpose.

Using these two properties we can sum up the series over $r$ in (S-38) to get

[TABLE]

where $\ket{\psi}=\sum_{k=1}^{m-1}\psi_{k}\ket{k}$ is a solution of a non-homogeneous recurrence relation of the fourth order with non-constant coefficients,

[TABLE]

In order to compute the overlaps explicitly, a solution of the recurrence relation (S-42) is now needed. Up to now we haven’t been able to solve it explicitly, however, we have guessed and extensively numerically checked the following:

Conjecture 1**.**

Let $\lambda$ be purely imaginary and $\eta=\frac{l\,\pi}{m}$, $l,\,m$ co-prime integers. Sum (S-38) of the overlaps between the local densities of $Z(\varphi)$ is

[TABLE]

For small $\lambda$ this expression correctly reproduces the continuous-time result – see tomaz14 . Since the remainder $p_{N}(\varphi)$ in the conservation laws $Z(\varphi)$ vanishes in the thermodynamic limit $N\to\infty$ we have

[TABLE]

and the full overlap between the conservation laws can be computed from (S-36). Using (S-31) we also deduce

[TABLE]

It is now easy to compute the overlaps between the spin-reversal anti-symmetric components $Z^{-}(\varphi)$ given in (S-30) and, separately, the magnetization-orthogonal components $Z_{\perp}^{-}(\varphi)$ given in (S-32):

[TABLE]

The last overlap is essentially Eq. (19) of the MT and is used in the Mazur inequality. Therefore, it deserves its own symbol $K(\varphi,\varphi^{\prime})\equiv\lim_{N\to\infty}\frac{1}{N}\langle Z_{\perp}^{-}(\bar{\varphi}),Z_{\perp}^{-}(\varphi^{\prime})\rangle$. Explicitly it reads

[TABLE]

Appendix G Appendix G: Solving the integral equation

In this appendix we describe how to tackle the integral equation given by Eq. (21) of the MT. Numerics suggests that it is enough to restrict the domain of integration to the line $\rm{Re}\,\varphi=\pi/2$, which can be understood in terms of analyticity of the kernel. Taking $\varphi=\pi/2+i\log x$, $\varphi^{\prime}=\pi/2+i\log x^{\prime}$ and $\lambda=i|\lambda|=i\log\Lambda$ we see, that the integral equation and the lower bound (Eqs. (21) and (20) of the MT, respectively) are equivalent to

[TABLE]

where

[TABLE]

and

[TABLE]

$K(\varphi,\varphi^{\prime})$ is given by (S-47), while Eq. (17) of the MT provides explicit formula for $j(\varphi)$. Function $\widetilde{f}(x)$ has absorbed the Jacobian $1/x$ of the coordinate transformation $\log x\to x$. The goal is now to solve for $\widetilde{f}(x)$ and plug it into the lower bound, provided the explicit form of real functions $\widetilde{K}(x,x^{\prime})$ and $\widetilde{j}(x)$

[TABLE]

Since $\widetilde{j}(x)=\widetilde{j}(x^{-1})$, changing the variables according to $x\to 1/x$ in the second equation of (S-48) suggests $\widetilde{f}(x^{-1})=x^{2}\widetilde{f}(x)$. Using this symmetry, we can rewrite the first equation of (S-48) as

[TABLE]

where

[TABLE]

Taking the limit $m\to\infty$ now gives

[TABLE]

If we define $g(x)=(1+2\Lambda x^{-2}+\Lambda^{2})G(x)$, we have $g(x)+g(x^{-1})=1$. So we can assume the following expansion

[TABLE]

for some coefficients $c_{m}$. If we use this in the second equation of (S-56) we see that $G(x)=x^{4}G(x^{-1})$ must hold for $\widetilde{f}(x^{-1})=x^{2}\widetilde{f}(x)$ to be true. Together with (S-54) this is enough to solve for $G$ and then $\widetilde{f}$, which produces

[TABLE]

It is now easy to compute the lower bound (S-48) which, after substitution $\Lambda=e^{|\lambda|}$ becomes

[TABLE]

This is the enveloping function given in Eq. (23) of the MT. ${\rm Gd}$ denotes the Gudermannian function with explicit formula ${\rm Gd}(x)=2\arctan(e^{x})-\pi/2$.

Appendix H Appendix H: Numerical approaches

Mazur bound calculation.– In order to compute the Drude weight bound numerically (for general $m$), we discretize the region of integration ($|{\rm Re}\,\varphi-\frac{\pi}{2}|<\frac{\pi}{2m}$) into $L$ equally-sized rectangles and choose a cut-off value $|{\rm Im}\,\varphi|\leq\mathcal{C}$. In our experience cut-offs $\mathcal{C}=5$ and $\mathcal{C}=10$ are already indistinguishable and thus the first one ($\mathcal{C}=5$) is sufficient for accurate results. It also appears we can limit ourselves to the line ${\rm Re}\varphi=\frac{\pi}{2}$ as increasing the number of points in the real direction does not affect the value significantly.

Having discretized the region, the integral equation (Eq. (21) of the MT) reduces to a system of linear equations which can be written in matrix form

[TABLE]

and solved efficiently for $\boldsymbol{f}$. Here, $\boldsymbol{K}$ represents a square matrix with elements $K(\varphi,\varphi^{\prime})$ which are computed from (S-47) at the points $\varphi$ and $\varphi^{\prime}$ defining our discretization. These points act as a row and a column index, respectively. Similarly $\boldsymbol{f}$ and $\boldsymbol{j}$ are vectors with components $f(\varphi)$ and $j(\varphi)$, respectively (see Eq. (17) of the MT).

The solution $\boldsymbol{f}$ is now used in a simple scalar product

[TABLE]

which gives the lower bound (compare with Eq. (20) of the MT). Here $\boldsymbol{j^{\prime}}$ is a permuted vector $\boldsymbol{j}$, its components being $j(\bar{\varphi})$ instead of $j(\varphi)$. Symbol $(\bullet)^{\boldsymbol{\dagger}}$ denotes the conjugate transpose.

It appears that even relatively crude discretizations of roughly $L=200$ points along the relevant line provide excellent results, essentially indistinguishable from more precise discretizations. Still, we perform calculations with $L$ up to $5000$ to verify these results.

Tensor network simulations.– In order to check how well our bound describes the actual value we attempt to obtain numerical values for the spin Drude weight through an alternative approach. To this end, we simulate an inhomogeneous quench where the two halves of the system are initially prepared in a slightly polarized product state

[TABLE]

We take $\mu_{R}=-\mu_{L}$ with $\mu_{L}=\mu/2$ typically being small (usually of order $10^{-2}$ or smaller), since this increases the stability of the simulation. In order to obtain the Drude weight, we use the linear response expression

[TABLE]

recently used in vasseur15 ; karrasch17 ; ilievski17 , with $J$ being the total extensive current. We compute its expectation value using a TEBD algorithm. Formula (S-63) differs from (S-7), derived in Appendix A, since this particular protocol gives rise to nontrivial dynamics only inside the light-cone, emerging from the contact between the two halves of the chain.

Typically we simulate chains of length $N=3600$ and $N=7200$ if longer times are needed. The times we can reach are roughly a quarter of the chain’s length, at which point the boundary effects may start affecting the accuracy of our Drude weight estimates. Finally, the size of the bias $\mu$ can be taken to be around $10^{-2}$. Taking a smaller magnetization step appears not to be necessary. For the most part, we take a bond dimension of $\chi=64$ which appears to be sufficient to obtain moderately accurate results with errors in the $\sim 1\%$ range (as compared to a somewhat more precise simulation at $\chi=128$ or even $\chi=256$ which were made occasionally for comparison). Increasing the bond dimension helps increase the Drude weight at points where the less precise calculation might duck under the theoretical lower bound. This can be seen at the peaks at commensurable $\eta/\pi$ (”the fractal spikes”), an example being shown in Figure S-2.

We stress that the numerical simulations clearly suggest convergence to a limiting fractal spin Drude weight as given by the Mazur bound. This convergence is quantitatively illustrated in the inset of Figure S-2, where the full width at half maximum (FWHM) of finite-$t$ DMRG data converges to zero as $\sim t^{-1/2}$.