Aperiodically driven integrable systems and their emergent steady states

TL;DR

This paper investigates how breaking periodic driving with aperiodic noise affects the long-term behavior of integrable quantum systems, revealing new steady states and the conditions under which they form.

Contribution

It introduces the concept of a geometric generalized Gibbs ensemble and analyzes the impact of different types of aperiodic noise on steady state formation in integrable systems.

Findings

Random noise leads to infinite temperature thermalization.

Self-similar noise results in a novel steady state, the geometric GGE.

Large time scales are required for the emergence of these steady states.

Abstract

Does a closed quantum many-body system that is continually driven with a time-dependent Hamiltonian finally reach a steady state? This question has only recently been answered for driving protocols that are periodic in time, where the long time behavior of the local properties synchronize with the drive and can be described by an appropriate periodic ensemble. Here, we explore the consequences of breaking the time-periodic structure of the drive with additional aperiodic noise in a class of integrable systems. We show that the resulting unitary dynamics leads to new emergent steady states in at least two cases. While any typical realization of random noise causes eventual heating to an infinite temperature ensemble for all local properties in spite of the system being integrable, noise which is self-similar in time leads to an entirely different steady state, which we dub as "geometric generalized Gibbs ensemble", that emerges only after an astronomically large time scale. To understand the approach to steady state, we study the temporal behavior of certain coarse-grained quantities in momentum space that fully determine the reduced density matrix for a subsystem with size much smaller than the total system. Such quantities provide a concise description for any drive protocol in integrable systems that are reducible to a free fermion representation.

Figures

Equations

$$H(t)=H(t+nT),$$

$$|\psi (nT)⟩=U(T{)}^n|\psi (0)⟩,$$

$$H(t)=-\sum _{j=1}^L\left[g(t)s_j^x+s_j^z{s}_{j+1}^z\right],$$

$$s_n^x$$

$$s_n^z$$

$$H=P^{+}{H}^{+}{P}^{+}+P^{-}{H}^{-}{P}^{-},$$

$$H^{\pm }$$

$$c_k$$

$$c_n$$

$$H^{+}$$

$$H_k$$

$$H_k=2(g(t)-\cos k){\sigma }_k^z+2\sin k{\sigma }_k^x.$$

$$i\frac{\mathrm{d}}{\mathrm{dt}}|{\psi }_k⟩=H_k(t)|{\psi }_k⟩,$$

$$|\psi (t)⟩$$

$$|{\psi }_k(t)⟩$$

$$g_{\mathrm{ref}}$$

$$U_k(T)=\exp \left(-i{H}_k(g_f)\frac{T}{2}\right)\exp \left(-i{H}_k(g_i)\frac{T}{2}\right),$$

$$|{\psi }_k(n)⟩=U_k(T+{\tau }_{n}dT)U_k(T+{\tau }_{n-1}dT)\cdots U(T+{\tau }_{1}dT)|{\psi }_k(0)⟩=\mathcal{T}\prod _{i=1}^n{U}_k(T+{\tau }_{i}dT)|{\psi }_k(0)⟩,$$

$$m=0:$$

$$m=1:$$

$$m=2:$$

$$m=3:$$

$$m=4:$$

$$⋮$$

$${\rho }_l={\mathrm{Tr}}_{L-l}(\rho ),$$

$$C_{ij}$$

$$F_{ij}$$

$${\mathcal{C}}_n(l)$$

$$S_{\mathrm{ent}}(l)$$

$${\mathcal{D}}_n(l)=\frac{1}{2l}\mathrm{Tr}\sqrt{[{\mathcal{C}}_{\infty }(l)-{\mathcal{C}}_n(l){]}^{†}[{\mathcal{C}}_{\infty }(l)-{\mathcal{C}}_n(l)]}.$$

$${\left(|u_k(n){|}^2\right)}_c$$

$${\left(u_k^{\ast }(n)v_k(n)\right)}_c$$

$$\cos [k(i-j)]$$

$$\sin [k(i-j)]$$

$$C_{ij}$$

$$F_{ij}$$

Full text

I Introduction and Motivation

Statistical mechanics can be derived from Jaynes’ principle of maximum entropy Jaynes1957a ; Jaynes1957b given the constants of motion for any generic many-body system. Remarkably, this principle may also be extended to the level of a single eigenstate in a thermodynamically large, isolated quantum system leading to the eigenstate thermalization hypothesis  Deutsch1991 ; Srednicki1994 ; RigolDO2008 . The eigenstate thermalization hypothesis implies that the expectation values of all operators which are local in space equal those obtained from a maximum entropy statistical description when appropriate conserved quantities are taken into account for constructing the ensemble. The Hamiltonian is the sole such conserved quantity in a generic chaotic system. Exceptions to the eigenstate thermalization hypothesis are provided by disordered systems that show localization, both for non-interacting models Anderson1958 and with interactions BaskoAA2006 (i.e., many-body localized systems NandkishoreH2015 ). Integrable models provide examples of systems that are not localized and yet do not satisfy the eigenstate thermalization hypothesis. However, Jaynes’ principle may again be applied here by considering the extensive number of conservation laws due to integrability that then leads to a generalized Gibbs ensemble (GGE) RigolDYO2007 ; CassidyCR2011 ; CauxE2013 instead of the standard Gibbs ensemble. These statistical descriptions are expected to hold only for local properties (i.e., properties that are determined by the degrees of freedom in a subsystem that is much smaller than the rest of the system, as shown schematically in Fig. 1), and not for the wave function of the entire system which is after all in a pure state PopescuSW2006 .

Tremendous experimental progress in producing and manipulating well-isolated quantum systems such as ultracold quantum gases Blochreview2005 ; Goldman2014 ; Langenetalreview2015 ; Eck17 , has also led to great interest in driven closed quantum systems. Removing the time translational invariance of closed quantum systems via a time-dependent Hamiltonian leads to the much richer possibility of stabilizing non-equilibrium steady states with purely unitary dynamics. Such non-equilibrium steady states may even have properties otherwise forbidden in thermal equilibrium, such as spontaneous time-translational symmetry breaking in many body localized systems ElseBN2016 ; KhemaniLMS2016 and dynamical topological ordering KitagawaBRD2010 ; NathanR2015 . Many-body Floquet systems where the time-dependent Hamiltonian satisfies

[TABLE]

with $T$ being the time period of the drive and $n$ being integer-valued, is a well-studied case where such non-equilibrium steady states have been shown to exist. More precisely, the local properties of Floquet systems eventually synchronize with the period of the drive that then allows the possibility of a “periodic” ensemble LazaridesDM2014a ; LazaridesDM2014b ; Ponte15b ; AlessioR2014 . Propagating the wave function $|\psi\rangle$ of the system in a stroboscopic manner

[TABLE]

where $U(T)$ is the time evolution operator for one time period, generates a discrete quantum map indexed by $n$, and the steady state emerges as $n\to\infty$. This steady state can again be described using Jaynes’ principle LazaridesDM2014a where only quantities that are stroboscopically conserved need to be taken into account. In the non-integrable case, it is expected that no such conserved quantity exists NandkishoreH2015 since the Hamiltonian is no longer conserved; hence the system should locally mimic an infinite temperature ensemble in the steady state by entropy maximization. This is indeed consistent with the results of several recent studies LazaridesDM2014b ; Ponte15b ; AlessioR2014 ; Russomanno12 ; Bukov15 ; Ponte15a ; Eckardt15 , but a complete understanding is still lacking. The situation is markedly different for integrable systems. For systems that are reducible to free fermions e.g., the prototypical one-dimensional (1D) spin-1/2 transverse field Ising model, it has been shown LazaridesDM2014a that an extensive number of quantities that are stroboscopically conserved continue to exist even in the driven case; these prevent the system from approaching the infinite temperature ensemble locally and instead lead to a periodic GGE (p-GGE) description of the synchronized local properties at late times. The fate of other types of integrable systems like Bethe-integrable ones (for instance, the 1D spin-1/2 $XXZ$ model) under periodic drives is still an open issue BurkovH2012 .

I.1 Synopsis

In this work, we consider integrable models that are reducible to free fermions, and our objective is to understand whether well-defined non-equilibrium steady states exist for continually driven protocols that are not periodic functions in time. We restrict ourselves to drives that can be represented by discrete quantum maps similar to Eq. (2) for ease of analysis. For concreteness, we work here with the 1D transverse field Ising model which is a standard many-body integrable quantum system Subirbook ; SIBKCbook ; Duttabook . The time-dependent Hamiltonian is given by

[TABLE]

where $s^{x,y,z}$ are the usual Pauli operators, $L$ denotes the number of spins in the system under consideration, and we assume periodic boundary conditions in space. $g(t)$ represents the time-dependent transverse magnetic field through which the system is continually driven.

We start with a reference function $g_{\mathrm{ref}}(t)$ defined for $0\leq t\leq T$. Repeatedly using the corresponding time evolution operator $U(T)$ for all values of $n>0$ in Eq. (2) results in periodic driving with period $T$. Consider now a different drive protocol where $g(t)$ is constructed by patching together rescaled versions of $g_{\mathrm{ref}}(t)$ in time as a function of $n$ such that the period of $g_{\mathrm{ref}}(t)$ is stretched or reduced by $dT$ depending on $n$. The corresponding time evolution operator is then denoted by $U(T+dT)$ or $U(T-dT)$. The new drive protocol $g(t)$ can now be represented by a sequence in $n$ that takes a value of either $+1$ (representing $U(T+dT)$) or $-1$ (representing $U(T-dT)$) at each $n$. There are an infinite number of ways to choose the sequence of $\pm 1$ as a function of $n$, each of which represents a different purely unitary dynamics, such that the resulting driving protocol $g(t)$ cannot be represented by any periodic function in time. We note that the protocol described above involving two time periods $T+dT$ and $T-dT$ is mathematically equivalent to a protocol in which the time period $T$ is kept fixed but the Hamiltonian is scaled globally by factors of $1+dT/T$ and $1-dT/T$ respectively. The important point about the protocol is that it involves two possible time evolution operators which do not commute with each other.

We will work out the nature of the resulting non-equilibrium steady states as $n\to\infty$ in two cases here. First, we will consider a typical realization of a random sequence of $\pm 1$ where the sign is chosen with probability 1/2 at each $n$ with $dT\ll T$. This mimics a periodic drive protocol perturbed with random noise due to the lack of perfect control over the time variation of $g(t)$ in an experimental set-up. Secondly, we will consider a different sequence of $\pm 1$ that is neither random nor periodic in $n$ but has a fractal structure instead, and thus can be thought of as an example of scale-invariant noise.

For a periodic drive protocol perturbed with random noise where $dT\ll T$, we will show that the system, after initially approaching the p-GGE that corresponds to $dT=0$, eventually heats up to an infinite temperature ensemble as $n\to\infty$ in spite of the underlying integrability of the model. The time spent by the system close to the p-GGE, which can now be thought of as a quasisteady prethermal state, and the eventual manner of heating up to an infinite temperature ensemble can both be understood from a simple random walk argument, provided that $dT\ll T$. We will further show that a coarse-graining procedure in momentum space explains how the reduced density matrix of any subsystem with $l$ adjacent spins approaches an infinite temperature ensemble as long as $l\ll L$ even though the dynamics is completely unitary. It is well known that weakly interacting systems (i.e., systems close to an integrable point) after a sudden perturbation show a prethermalization regime that can be understood by a GGE based on perturbatively constructed constants of motion KollarWE2011 , before eventually thermalizing on a much longer time scale MoeckelK2008 ; Gring2012 with the energy being the only constant of motion. In the perturbed Floquet integrable systems with random noise that we consider, the integrability of the model is never broken but the stroboscopically conserved quantities at $dT=0$ no longer remain conserved when $dT\neq 0$.

However, it is important to note that the lack of conservation of these quantities does not fix the form of the final non-equilibrium steady state for an aperiodic drive protocol. For a periodic drive protocol perturbed with scale-invariant noise, we will show that if the local properties are observed not as a function of $n$, but instead as a function of $2^{n}$ (thus, geometrically instead of stroboscopically), a new steady state emerges when $n\to\infty$ which we call “geometric” GGE (g-GGE). Interestingly, this steady state can also be understood using an effective p-GGE construction even though the drive in not periodic in time because of a fixed point property of products of $2\times 2$ $SU(2)$ matrices. An extensive number of conserved quantities appear when the system is viewed after every $2^{n}$ drives, but only when $n$ is sufficiently large, which shows that these conservations are emergent.

The rest of the paper is arranged in the following manner. In Sec. II, we review some results relevant for our work and set the notations for the rest of the paper. Sec. II.1 details the pseudospin representation for the 1D transverse field Ising model that allows the wave function of the system to be denoted in terms of points on Bloch spheres in momentum space. Sec. II.2 details the aperiodic driving protocols that we study. Sec. II.3 reviews the way to obtain the reduced density matrix of a subsystem of $l$ adjacent spins using appropriate correlation matrices. In Sec. III (including Sec. III.1 and Sec. III.2), we introduce certain coarse-grained quantities in momentum space which fully determine the reduced density matrix of a subsystem when $l\ll L$ and thus all local quantities, and we discuss their behavior for time-periodic drive protocols. In Sec. IV, we present results for perturbed Floquet systems with random noise and use a random walk argument to show how the system heats up to an infinite temperature ensemble in spite of its integrability. In Sec. V, we consider perturbed Floquet systems with noise which is scale-invariant in time and show that for a certain type of scale-invariant noise, the system reaches a steady state which is not an infinite temperature ensemble, but it does so only after an astronomically large number of drives. The final steady state is explicitly constructed using a fixed point property of strings of $2\times 2$ $SU(2)$ matrices. Finally, we summarize our results and conclude with some future directions in Sec. VI.

II Preliminaries

II.1 Pseudospin representation of dynamics for the 1D transverse

field Ising model

To solve the 1D transverse field Ising model (Eq. (3)), we introduce the standard Jordan-Wigner transformation of spin-1/2’s to spinless fermions Lieb61

[TABLE]

Writing $H$ in Eq. (3) in terms of these fermions, we obtain

[TABLE]

where $P^{\pm}$ are projectors on to subspaces with even ($+$) and odd ($-$) number of fermions, and

[TABLE]

with antiperiodic (periodic) boundary conditions for the fermions for $H^{+}$ ($H^{-}$). We now focus on the case of even $L$ with antiperiodic boundary conditions for the fermions (which corresponds to periodic boundary conditions for the spins). To go to $k$ space, we define

[TABLE]

where $k=2\pi m/L$ with $m=-(L-1)/2,\cdots,-1/2$, $1/2,\cdots,(L-1)/2$ for even $L$. Rewriting $H^{+}$ in terms of $c_{k}$ and $c_{k}^{\dagger}$, we get

[TABLE]

This Hamiltonian connects the vacuum $|0\rangle$ of the $c$ fermions with the two-particle state $|k,-k\rangle=c_{k}^{\dagger}c^{\dagger}_{-k}|0\rangle$, and $|k\rangle=c_{k}^{\dagger}|0\rangle$ with $|-k\rangle=c_{-k}^{\dagger}|0\rangle$. It is more convenient to write Eq. (8) as $H=\sum_{k>0}H_{k}$ where

[TABLE]

We now introduce a “pseudospin representation” $\sigma_{k}$ (which is different from the original $s$ matrices in Eq. (3)) where $|\uparrow\rangle_{k}=|k,-k\rangle=c_{k}^{\dagger}c^{\dagger}_{-k}|0\rangle$ and $|\downarrow\rangle_{k}=|0\rangle$ KolodrubetzCH2012 . Writing $H_{k}$ in Eq. (9) in this basis, we see that

[TABLE]

Then the pseudospin state $|\psi_{k}\rangle$ for each $k$ mode evolves independently as

[TABLE]

where

[TABLE]

Thus, specifying $(u_{k}(t),v_{k}(t))^{T}$ (where the superscript $T$ denotes the transpose of the row $(u_{k}(t),v_{k}(t))$) for each allowed $k>0$ for a finite $L$ specifies the complete wave function $|\psi(t)\rangle$ through Eq. (12). The wave function of the system can be equivalently characterized by a two-dimensional state space at each allowed positive $k$, and $(u_{k}(t),v_{k}(t))^{T}$ can then be represented as a point that evolves with time according to Eq. (11) on the corresponding Bloch sphere.

In the rest of the paper, we will take $(u_{k},v_{k})^{T}=(0,1)^{T}$ to be the starting state at each $k$ from which the unitary dynamics proceeds. This corresponds to the system being initially prepared in the pure state where all the spins have $\sigma^{x}=+1$ (which corresponds to the ground state for $g\to\infty$).

II.2 Details of the driving protocol

While our basic conclusions are protocol independent and do not depend on the exact nature of $g_{\mathrm{ref}}(t)$, we take it to have a square pulse variation in time for mathematical convenience, namely,

[TABLE]

The corresponding unitary time evolution operator $U_{k}(T)$ for the pseudospin mode with momentum $k$ equals

[TABLE]

where $H_{k}$ has the form given in Eq. (10). For the perturbed Floquet dynamics that we consider in this work, we need to define $U_{k}(T\pm dT)$; these follow in a straightforward manner from Eq. (14) by replacing $T\to T\pm dT$.

Rewriting Eq. (2) for the perturbed Floquet dynamics for the mode with momentum $k$, we get

[TABLE]

where the sequence ${\tau_{i}}=\tau_{1},\tau_{2},\tau_{3},\cdots$ (with each $\tau_{i}$ being equal to either $+1$ or $-1$) is the same for all the allowed $k$ modes at a finite $L$. (In the last expression in Eq. (15), ${\cal T}$ denotes time-ordering).

For Floquet systems perturbed with random noise, the sequence ${\tau_{i}}$ can be taken to be any typical realization of a random process where each element is chosen to be $\pm 1$ randomly and independently with probability 1/2; we will discuss the results for such a realization later in the paper. For Floquet systems perturbed with scale-invariant noise, we take the well-known Thue-Morse sequence  Thue1906 ; Morse1921a ; Morse1921b which is an infinite sequence of $\tau_{i}=\pm 1$ that is obtained by starting with $\tau_{1}=-1$ and successively appending the negative of the sequence obtained thus far. The first few steps of this recursive procedure yield

[TABLE]

At each recursion level $m$, we thus obtain the first $2^{m}$ elements of this infinite sequence. The self-similar structure of the Thue-Morse sequence can be verified by removing every second term from it, which then results in the same infinite sequence. This sequence is thus an example of a quasiperiodic sequence in time which is neither periodic nor random. Note that the average time period of this perturbed Floquet system equals $T$ (the time period of the unperturbed Floquet system), just like in the random noise case.

II.3 Reduced density matrices and the distance measure

Given the wave function $|\psi(n)\rangle$ of the entire system, which can be described by the corresponding density matrix $\rho=|\psi(n)\rangle\langle\psi(n)|$, all local properties within a length scale $l$ can be understood by considering the reduced density matrix of a subsystem of $l$ adjacent spins (Fig. 1); we denote this by $\rho_{l}$ which is given by

[TABLE]

where all the other $L-l$ degrees of freedom have been integrated out. Even though $\rho$ is a pure density matrix at any $n$, the reduced density matrix $\rho_{l}$ is mixed since the pure state $|\psi(n)\rangle$ typically gets more entangled as $n$ increases.

The reduced density matrix of $l$ adjacent spins when $|\psi(n)\rangle$ has the form given in Eq. (12) is determined in terms of the correlation functions of the $c$ fermions at these $l$ sites VidalLRK2003 . For free fermions, since all higher-point correlation functions are expressible using two-point correlations from Wick’s theorem, the reduced density matrix can be determined from the knowledge of two $l\times l$ matrices Peschel2003 , $\mathbf{C}$ and $\mathbf{F}$, whose elements are constructed from $u_{k}(n)$ and $v_{k}(n)$ as follows:

[TABLE]

where $i,j$ refer to sites in the subsystem. Using these expressions, we construct the $2l\times 2l$ matrix ${\mathcal{C}}_{n}(l)$ given by

[TABLE]

The reduced density matrix $\rho_{l}$ is then completely fixed by the eigenvalues and eigenvectors of $\mathcal{C}_{n}(l)$. For example, the entanglement entropy $S_{\mathrm{ent}}(l)$ equals:

[TABLE]

where $p_{k}$ denotes the $k$-th eigenvalue of $\mathcal{C}_{n}(l)$.

If a non-equilibrium steady state exists as $n\to\infty$, then the convergence of the local properties to their final steady state values as $n$ increases can be characterized by defining a distance measure $\mathcal{D}_{n}(l)$ as

[TABLE]

We note that $0\leq\mathcal{D}_{n}(l)\leq 1$ and is identically zero only if $\mathcal{C}_{n}(l)=\mathcal{C}_{\infty}(l)$. For a periodically driven protocol with $dT=0$, $\mathcal{C}_{\infty}(l)$ can be calculated using the appropriate p-GGE; the distance measure decays to zero as $\mathcal{D}_{n}(l)\sim n^{-\alpha(T)}$ for large $n$ (and $l\ll L$), with $\alpha(T)=1/2$ or $3/2$ depending on the time period $T$ of the drive SenNS2016 . We note that $\mathcal{D}_{n}(l)$ does not decay to zero as $n\to\infty$ if $l/L$ is finite; this explicitly shows that only local quantities can be described by the p-GGE.

III Coarse-graining in momentum space

Even though the entire wave function $|\psi(n)\rangle$ needs the specification of $(u_{k}(n),v_{k}(n))^{T}$ at $k=2\pi m/L$ with $m=1/2,3/2,\cdots,(L-1)/2$, $\mathcal{C}_{n}(l)$ and consequently the reduced density matrix for any $l\ll L$ depends only on certain coarse-grained variables defined as follows LaiY2015 ; NandySDD2016

[TABLE]

These variables are defined using $N_{c}(\gg 1)$ consecutive $k$ modes that lie within a cell ($k_{\mathrm{cell}}$) which has an average momentum $k_{c}$ in momentum space and a size $\delta k$ where $1/L\ll\delta k\ll 1/l$. With this condition on $\delta k$, the sinusoidal factors in Eq. (18) can be replaced to a very good approximation by

[TABLE]

for all the $k$ modes within the cell. The sum over momentum in Eq. (18) can thus be carried out in two steps, first summing over the consecutive momenta in a single cell, and then summing over the different momentum cells. This gives

[TABLE]

where $\mathcal{N}_{\mathrm{cell}}\ll L/2$ represents the number of momentum cells after the coarse-graining procedure.

Note that even though quantities like $|u_{k}(n)|^{2}$ and $u^{*}_{n}(k)v_{n}(k)$ continue to display Rabi oscillations and thus have no well-defined $n\to\infty$ limit, the coarse-grained variables $(|u_{k}(n)|^{2})_{c}$ and $(u^{*}_{k}(n)v_{k}(n))_{c}$ must reach final steady state values for a non-equilibrium steady state to exist.

If the system is locally described by an infinite temperature ensemble, then $\rho_{l}=I$, where $I$ is the identity matrix, for any $l\ll L$. This immediately implies that $(|u_{k}(n)|^{2})_{c}\to 1/2$ (also $(|v_{k}(n)|^{2})_{c}\to 1/2$) and $(u^{*}_{n}(k)v_{n}(k))_{c}\to 0$ for all the coarse-grained momentum cells as $n\to\infty$, since $C_{ij}=(1/2)\delta_{i,j}$ and $F_{ij}=0$ for $\rho_{l}=I$. (A similar behavior has been found for repeated quenching of the transverse field in the 1D spin-1/2 $XY$ model Mukherjee08 ). We calculate the $n$ dependence of such coarse-grained quantities numerically by taking a sufficiently large system size $L$ (typically $L=131072$, but $L=524288$ in some cases), so that the allowed momentum modes are sufficiently dense in the interval $k\in[0,\pi]$, and then dividing the Brillouin zone into $\mathcal{N}_{\mathrm{cell}}=32$ cells, each of which has $N_{c}\gg 1$ consecutive momenta. For the periodic drive case where $dT=0$, we see that although $(|u_{k}(n)|^{2})_{c}$ indeed converges to a finite value (Fig. 2) and thus has a well-defined non-equilibrium steady state, this constant depends on the value of the average momentum of the coarse-grained cell ($k_{c}$) and is different from 1/2, showing that the local description is different from that of an infinite temperature ensemble.

III.1 Probability distribution on the Bloch sphere

At each $k$, the pseudospin state $|\psi_{k}(n)\rangle$ obtained after $n$ drives can be parameterized in terms of three angles $(\theta_{k}(n),\phi_{k}(n),\beta_{k}(n))$,

[TABLE]

The overall phase $\beta_{k}(n)$ is not important and we will ignore it henceforth. We will concentrate on the angles $\theta_{k}(n)$ and $\phi_{k}(n)$ which define the motion of a point on the unit Bloch sphere for each momentum $k$. From Eq. (14), it is clear that $U_{k}(T)$ is an element of the group $SU(2)$ (namely, $2\times 2$ unitary matrices with determinant equal to 1), and it can be written as

[TABLE]

where $0\leq\gamma_{k}\leq\pi$ and ${\hat{e}}_{k}=(e_{1k},e_{2k},e_{3k})$ is a unit vector. For a drive protocol that is perfectly periodic, $(\theta_{k}(n),\phi_{k}(n))$ traces out a circle on the Bloch sphere formed by the intersection of this unit sphere with the plane that contains the point $(\theta_{k}(0),\phi_{k}(0))$ and whose normal vector equals ${\hat{e}}_{k}$. It is then clear that the quantities $\mathcal{I}_{k}$ defined as

[TABLE]

are all independent of $n$ for a periodic drive protocol and thus define the extensive number of stroboscopically conserved quantities.

The behaviors of the coarse-grained quantities in Eq. (24) depend on the simultaneous positions of $(\theta_{k}(n),\phi_{k}(n))$ on the unit sphere of all the momentum modes that lie within a coarse-grained cell. Since their number $N_{c}\gg 1$ when $L\gg 1$, it is useful to instead look at the probability distribution of these $N_{c}$ points on the unit sphere. In particular, we will study the probability distribution of $|u_{k}(n)|^{2}-|v_{k}(n)|^{2}=\cos\theta_{k}$ for each such coarse-grained momentum cell, which we denote as $P_{n}(\cos\theta_{k})$, as a function of $n$. We can also define another probability distribution from the motion of a single $k$ mode for a large number of drive cycles and denote it by $\overline{P}(\cos\theta_{k})$. Remarkably, we find that

[TABLE]

for large $n$, where $k_{c}$ denotes the average momentum of a coarse-grained cell that has $N_{c}$ momenta on the left hand side of Eq. (30), and the right hand side is obtained from the motion of a single momentum mode at $k_{c}$.

For an infinite temperature ensemble, we expect that the $N_{c}$ points will cover the unit sphere uniformly at large $n$ even though all these points started from the south pole at $n=0$ due to the choice of $|\psi(0)\rangle$. Thus, $P_{n}(\cos\theta_{k_{c}})$ at large $n$ should be independent of the value of the average momentum of the cell, $k_{c}$, and thus equal 1/2, if the local description is that of an infinite temperature ensemble. We will derive $\overline{P}(\cos\theta_{k_{c}})$ for the perfectly periodic case in the next subsection and see that numerically $P_{n}(\cos\theta_{k_{c}})$ approaches $\overline{P}(\cos\theta_{k_{c}})$ for large $n$. It turns out that $\overline{P}(\cos\theta_{k_{c}})$ depends on $k_{c}$; this again shows that the description is not locally that of an infinite temperature ensemble for periodic driving.

III.2 Probability distribution for periodic driving: $dT=0$

We will first derive the form of $\overline{P}(\cos\theta_{k})$ when there is no randomness, i.e., when $dT=0$. In this case, the time evolution operator after one time period $T$ for the momentum mode $k$ (where $0\leq k\leq\pi$) is given by

[TABLE]

Since $U_{k}(T)$ is an element of the group $SU(2)$, it can be written in the form of Eq. (28), where $\gamma_{k}$ and the unit vector ${\hat{e}}_{k}=(e_{1k},e_{2k},e_{3k})$ can be derived by using the expression given in Eq. (LABEL:uk). Then the two-component spinor obtained after $n$ drives is given by

[TABLE]

Hence

[TABLE]

where we have used the identity $e_{1k}^{2}+e_{2k}^{2}+e_{3k}^{2}=1$. We now assume that $\gamma_{k}/\pi$ is an irrational number; this is justified since rational numbers lying in the range $[0,1]$ form a set of measure zero. Hence the set of numbers $2n\gamma_{k}~{}{\rm mod}~{}\pi$ covers the region $[0,\pi]$ uniformly as $n\to\infty$. We now use the fact that if $2n\gamma_{k}$ is uniformly distributed from 0 to $\pi$, the probability distribution of the variable $v=\cos(2n\gamma_{k})$ is given by

[TABLE]

for $-1\leq v\leq 1$. [This is because $P(v)=(1/\pi)\int_{0}^{\pi}~{}d\theta~{}\delta(v-\cos\theta)=1/(\pi\sqrt{1-v^{2}})$.] Using Eq. (41) in Eq. (40), we find that the probability distribution of $\cos\theta_{k}$ is given by

[TABLE]

for $-1\leq\cos\theta_{k}\leq 1-2e_{3k}^{2}$, while $\overline{P}(\cos\theta_{k})=0$ for $1-2e_{3k}^{2}<\cos\theta_{k}\leq 1$. We observe that this distribution has square root divergences at two points given by $\cos\theta_{k}=-1$ and $1-2e_{3k}^{2}$, and the range of $\cos\theta_{k}$ is given by $2(1-e_{3k}^{2})$. Further, the distribution is correctly normalized so that $\int_{-1}^{1}d(\cos\theta_{k})\overline{P}(\cos\theta_{k})=1$.

Given the values of the parameters $g_{i},~{}g_{f},~{}T$ and $k$, one can calculate ${\hat{e}}_{k}$ and therefore $e_{3k}$ from Eqs. (LABEL:uk) and (28). Eq. (42) then gives the probability distribution of $\cos\theta_{k}$ which is generated by a large number of drives, as shown in Fig. 3. We also calculate $P_{n}(\cos\theta_{k})$ for the 1D transverse field Ising model with $L=131072$ where the positive momenta in the Brillouin zone are divided into $\mathcal{N}_{\mathrm{cell}}=32$ cells, each thus containing $N_{c}=2048$ consecutive momenta; we display the result for a large $n$ ($=5\times 10^{5}$) in Fig. 3 which validates Eq. (30).

Eq. 30 can be analytically understood as follows. For large $n$, we have argued above that the right hand side of Eq. 30 is given by Eq. (42) with $k=k_{c}$. Now we look at the left hand side of Eq. 30; here we have to consider a large but fixed value of $n$ and average over a range of momenta $\delta k$ around $k=k_{c}$. We begin with the expressions

[TABLE]

For $\delta k\ll 1$, $\gamma_{k}$ will be close to $\gamma_{k+\delta k}$ and ${\hat{e}}_{k}$ will be close to ${\hat{e}}_{k+\delta k}$. However,

[TABLE]

and it is clear that if $n$ is very large, $n\gamma_{k}$ and $n\gamma_{k+\delta k}$ will not be close to each other. In other words, a small value of $\gamma_{k+\delta k}-\gamma_{k}$ gets amplified as $n$ increases but a small value of ${\hat{e}}_{k+\delta k}-{\hat{e}}_{k}$ does not get amplified. For $\delta k\ll 1$, we can therefore make the approximation of setting ${\hat{e}}_{k}={\hat{e}}_{k+\delta k}$, but we cannot set $\gamma_{k}=\gamma_{k+\delta k}$ if we are interested in the $n\to\infty$ limit. Hence,

[TABLE]

imply, using Eq. (40), that

[TABLE]

Now, we can write $\gamma_{k+\delta k}-\gamma_{k}=g_{k}\delta k$, where $g_{k}$ is a number of order 1. Suppose that $n$ is large enough that $2ng_{k}\delta k\gg\pi$. This implies that as $k$ goes over all the momenta in a cell lying in the range $k_{c}-\delta k/2$ to $k_{c}+\delta k/2$, $2n\gamma_{k}$ mod $\pi$ will cover the range $[0,\pi]$ uniformly. Hence in this cell with an average momentum $k_{c}$, the distribution of the variable $v=\cos(2n\gamma_{k})$ will again be given by Eq. (41). Arguments similar to the ones leading from Eq. (40) to Eq. (42) will then show that the left hand side of Eq. (30) is also given by Eq. (42).

We will now derive the value to which $(|u_{k}(n)|^{2})_{c}$ converges as $n\rightarrow\infty$ (Fig. 2). It is straightforward to see that $1/2-(|u_{k}(n)|^{2})_{c}=-(1/2)(\cos\theta_{k}(n))_{c}$, where the coarse-grained quantity $(\cos\theta_{k}(n))_{c}$ is given by

[TABLE]

As $n\rightarrow\infty$, using Eq. 30 and Eq. 42, we get

[TABLE]

This is in agreement with the results shown in Fig. 2 for $1/2-(|u_{k}(n)|^{2})_{c}$ which, in the limit $n\to\infty$, is equal to $(1/2)e_{3k_{c}}^{2}$.

IV Results for perturbed Floquet system with random noise

We will now study the unitary dynamics when $dT\ll T$ and the $\tau_{i}$’s in Eq. (15) are taken to be any typical realization of a random sequence of $\pm 1$’s which are chosen with probability 1/2 each. We again look at the behavior of the coarse-grained quantities to understand the nature of the resulting non-equilibrium steady state. We will see that in spite of the integrability of the model, the resulting non-equilibrium steady state is locally described by an infinite temperature ensemble as $n\to\infty$, unlike the case with $dT=0$.

The first quantity we look at is the overlap of a reference $k$ mode $|\psi_{k}(n)\rangle$ with other $k+\delta k$ modes that lie very close to it, i.e., where $\delta k\ll 1$, as a function of $n$ when the unitary dynamics is of the form given in Eq. (15) with the $\tau_{i}$’s taken from a typical realization of a random sequence of $\pm 1$’s. Fig. 4 shows a plot of $|\langle\psi_{k}(n)|\psi_{k+\delta k}(n)\rangle|^{2}$ versus $n$. At $n=0$, all the $|\psi_{n}(k)\rangle$ start from the same state $(0,1)^{T}$. We see that while $|\langle\psi_{k}(n)|\psi_{k+\delta k}(n)\rangle|^{2}$ stays localized between 1 and another value (which is derived in the next paragraph) for the perfectly periodic drive with $dT=0$ (see Fig. 4), the $dT\neq 0$ case is extremely sensitive to initial conditions; we then find that $|\langle\psi_{k}(n)|\psi_{k+\delta k}(n)\rangle|^{2}$ covers the entire range $[0,1]$ at large $n$ when the overlap with different $k+\delta k$ modes is considered where all $\delta k\ll 1$ (see Fig. 4). Note that if $\delta k\ll 1$, then each individual $U_{k+\delta k}(T+\tau_{n}dT)$ is only slightly different from $U_{k}(T+\tau_{n}dT)$ of the reference $k$ mode. This behavior of the overlaps already suggests that the coarse-grained quantities behave completely differently in the $dT\neq 0$ case at large $n$.

For the perfectly periodic drive with $dT=0$, we can understand the behavior of $|\langle\psi_{k}(n)|\psi_{k+\delta k}(n)\rangle|^{2}$ versus $n$ as follows. Following the arguments given in Eqs. (43-50), we find that

[TABLE]

For $n\gg 1/|\gamma_{k}-\gamma_{k+\delta k}|$, $\sin^{2}[n(\gamma_{k}-\gamma_{k+\delta k})]$ will oscillate between 0 and 1. Eq. (54) then implies that $|\langle\psi_{k}(n)|\psi_{k+\delta k}(n)\rangle|^{2}$ will oscillate between 1 and $e_{3k}^{2}$, and the range of values is given by $1-e_{3k}^{2}$. We note that this is exactly half the range of values of $\cos\theta_{k}$ which is given by $2(1-e_{3k}^{2})$ according to Eq. (42).

We will now study what happens to the probability distribution $P_{n}(\cos\theta_{k})$ for the unitary process with $dT\neq 0$. To evaluate this quantity, we use a chain with $L=524288$ spins where each coarse-grained cell in momentum space contains $N_{c}=8192$ consecutive momenta (hence, $\mathcal{N}_{\mathrm{cell}}=32$) and the spinor at each $k$ equals $(0,1)^{T}$ at $n=0$. The results for one such cell (with average momentum $k_{c}=31\pi/64$) are shown in Fig. 5. For small $n$, we see that $P_{n}(\cos\theta_{k})$ approaches the distribution for the perfectly periodic drive with the same $T$ but with $dT=0$ which has characteristic square-root divergences as shown in Eq. (42). However, as $n$ increases further, the probability distribution starts deviating strongly from this form and approaches a constant distribution eventually at large $n$ (see Fig. 5); this shows that the $N_{c}$ points that constitute the coarse-grained momentum cell are uniformly spread over the unit sphere. The lower panel of Fig. 5 shows the positions of the $N_{c}$ points on the unit sphere, from which $P_{n}(\cos\theta_{k})$ is constructed, for $n=100$, $n=10000$ and $n=200000$. We will now provide an understanding of this spreading based on the idea of random walks.

We define two unitary matrices $U_{k}(T-dT)$ and $U_{k}(T+dT)$ given by

[TABLE]

These are the time evolution operators for the two possible time periods $T-dT$ and $T+dT$ respectively. The time evolution operator $U_{nk}$ for $n$ drives is obtained by multiplying $n$ matrices each of which is randomly chosen to be either $U_{k}(T-dT)$ or $U_{k}(T+dT)$ with probability 1/2 each to obtain a particular realization of this random process as shown in Eq. (15).

Next, we parameterize the two matrices as

[TABLE]

For $dT=0$, these matrices reduce to $U_{k}(T)$ given in Eq. (LABEL:uk); in that case, we have seen that the evolution operator $U_{k}(T)$ corresponds to a rotation about a unit vector ${\hat{e}}_{k}$ given in Eq. (28). If $dT$ is non-zero but small, the unit vector corresponding to each drive will have a small random component. To quantify this, we note that if $dT$ is small, the unit vectors ${\hat{a}}_{k}$ and ${\hat{b}}_{k}$ given in Eq. (56) differ by small angles $\Delta\phi_{1k}=\arccos({\hat{e}}_{k}\cdot{\hat{a}}_{k})$ and $\Delta\phi_{2k}=\arccos({\hat{e}}_{k}\cdot{\hat{b}}_{k})$ from the unit vector ${\hat{e}}_{k}$. These two small angles are of the same order and we can denote their average as $\Delta\phi_{k}$. For $dT=0$ or $\sin k=0$, it is clear from Eqs. (LABEL:uk) and (55) that the unit vectors ${\hat{e}}_{k},~{}{\hat{a}}_{k}$ and ${\hat{b}}_{k}$ are identical, and therefore $\Delta\phi_{k}=0$. Hence, we see that

[TABLE]

when either $dT$ or $\sin k\to 0$, and Eq. (57) holds independently of the exact form of the drive protocol.

If the drives have time periods given by $T+dT$ and $T-dT$ randomly, we expect from the theory of random walks that after $n$ drives, the unit vector corresponding to $U_{nk}$ will differ from ${\hat{e}}_{k}$ by an angle proportional to $\sqrt{n}\Delta\phi_{k}$, assuming that $\Delta\phi_{k}$ is small (this will be true if $dT$ is small). However, since the unit vectors lie on a compact space (given by the points on the surface of the unit sphere), the deviation of the unit vector of $U_{nk}$ from ${\hat{e}}_{k}$ cannot go to $\infty$ as $n\to\infty$; rather, the unit vectors will cover the unit sphere uniformly giving rise to a uniform probability distribution $\overline{P}(\cos\theta_{k})=1/2$. The behavior of the probability distribution $P_{n}(\cos\theta_{k})$ shown in Fig. 5 for large $n$ is thus consistent with Eq. (30).

Furthermore, from the above argument, the spreading of the overlaps $|\langle\psi_{k}(n)|\psi_{k+\delta k}(n)\rangle|^{2}$ as a function of $n$ (as shown in Fig. 4) when many neighboring momentum modes very close to the reference mode $k$ are considered is controlled by a time scale $nT$, where $n(\Delta\phi_{k})^{2}\simeq n(dT)^{2}(\sin k)^{2}\sim\mathcal{O}(1)$; hence the time scale $nT$ varies with the momentum $k$. We explicitly verify this by showing the behavior of the corresponding coarse-grained variables $(|u_{k}(n)|^{2})_{c}$ (defined in Eq. (24)) in Fig. 6. The data is consistent with an exponential decay to 1/2; however, the rate of decay is clearly sensitive to the values of both $dT$ and $k_{c}$ (the average momentum of the coarse-grained cell) as we see in Fig. 6 (a). However, when the data is plotted versus $n(dT)^{2}(\sin k_{c})^{2}$ (see Fig. 6 (b)) instead of $n$, we clearly see that the time scale of the exponential decay to 1/2 is controlled by $\tau_{k_{c},dT}=1/\left[(dT)^{2}(\sin k_{c})^{2}\right]$ as predicted by the random walk argument. We have further checked that the coarse-grained quantity $(u^{*}_{k}(n)v_{k}(n))_{c}$ (Eq. (24)) relaxes to zero for large $n$ for any value of $k_{c}$.

The behaviors of both $P_{n}(\cos\theta_{k})$ in Fig. 5 and $(|u_{k}(n)|^{2})_{c}$ in Fig. 6 as functions of $n$ show the irreversible approach to an infinite temperature ensemble when local quantities are probed. Finally, we consider the implication of the dependence of the time scale $\tau_{k,dT}=1/\left((dT)^{2}(\sin k_{c})^{2}\right)$ (that controls the relaxation of $(|u_{k}(n)|^{2})_{c}$) on $dT$ and $k_{c}$ for the rate at which the trace distance of the reduced density matrix of a subsystem of $l$ consecutive spins from an infinite temperature ensemble (using the definition $\mathcal{D}_{n}(l)$ in Eq. (23)) approaches zero after a large number of drives $n$. Since the trace distance receives a contribution from all momenta $k$, it should approach zero as

[TABLE]

where $\mathcal{C}_{\infty}(l)$ in Eq. (23) is calculated using an infinite temperature ensemble.

For large $n$, the integral in (58) is dominated by the regions in $k$ where $\tau_{k,dT}$ is large, namely, $k=0$ and $\pi$ where $\sin k\to 0$. Absorbing $n$ in $k$ (or $\pi-k$) in the exponential, we see that the trace distance will scale with $n$ as $1/(\sqrt{n}dT)$ when $dT$ is small. Thus, $\mathcal{D}_{n,ITE}(l)\sim\mathcal{F}_{l}(n(dT)^{2})$, where $\mathcal{F}_{l}(x)\sim 1/\sqrt{x}$ for $x\gg 1$. Secondly, $\mathcal{F}_{l}(x)\sim\mathcal{O}(1)$ when $x\ll 1$ since in this limit the system must relax to the p-GGE that emerges for $dT=0$. Thus, the prethermal regime where the system initially approaches the p-GGE exists for time scales that depend on $dT$ as $n\sim 1/(dT)^{2}$ for $dT\ll T$. The numerical data for $dT\neq 0$ for the behavior of $\mathcal{D}_{n,ITE}(l)$ as a function of $x=n(dT)^{2}$ is consistent with both these expectations, i.e., $\mathcal{F}_{l}\sim 1/\sqrt{x}$ for $x\gg 1$ and $\mathcal{F}_{l}\sim\mathcal{O}(1)$ for $x\ll 1$, as shown in Fig. 7 for a fixed $T=0.2$ and different values of $dT$. It is also instructive to look at $\mathcal{D}_{n,GGE}$ for non-zero values of $dT$ (Fig. 7, inset) where $\mathcal{C}_{\infty}(l)$ is calculated using the $dT=0$ p-GGE as a function of $n$. For time scales $n\sim 1/(dT)^{2}$, $\mathcal{D}_{n,GGE}$ decreases as a function of $n$ which shows that the system approaches the p-GGE in the prethermal regime. However, at larger values of $n$, $\mathcal{D}_{n,GGE}$ starts increasing as a function of $n$ which explicitly shows that the large-$n$ steady state is not described by a p-GGE. The initial relaxation to a p-GGE in the prethermal regime scales either as $\mathcal{D}_{n,GGE}\sim n^{-3/2}$ (the case shown in Fig. 7, inset) or as $\mathcal{D}_{n,GGE}\sim n^{-1/2}$ depending on which dynamical phase the periodic drive protocol ($dT=0$) belongs to SenNS2016 .

For completeness, we point out that Ref. MarinoS2012, (see also Ref. RooszJI2016, ) also considered a noisy 1D transverse field Ising model though not periodically driven on average, and showed that the asymptotic steady state is an infinite temperature ensemble. However, they considered the behavior of quantities averaged over different realizations of noise which resembles an open quantum system with non-unitary dynamics. Here, we show that if one restricts to understanding only local quantities as functions of $n$ in a perturbed Floquet integrable system, even the unitary dynamics of a single typical realization of a random sequence leads to an infinite temperature ensemble.

V Results for perturbed Floquet system with scale-invariant noise

We will now study the unitary dynamics when $dT\neq 0$ but the $\tau_{i}$’s in Eq. (15) are not a random sequence of $\pm 1$ but are instead given by the scale-invariant Thue-Morse sequence in Eq. (16) that is neither periodic nor random. The average time-period of such a drive protocol continues to equal $T$. To this end, we first show the coarse-grained quantity $(|u_{k}(n)|^{2})_{c}$ in Eq. (24) as a function of $n$ where the $\tau_{i}$’s are chosen to be $\pm 1$ according to the Thue-Morse, with $T=0.2$ and $dT=0.05$ (Fig. 8). We see that these quantities behave entirely differently from either the periodic case ($dT=0$) in Fig. 2 or the randomly perturbed case in Fig. 6. From Fig. 8, it is not even clear whether a well-defined non-equilibrium steady state exists as $n\to\infty$. We will show below that in fact a non-equilibrium steady state does exist in the large-$n$ limit if the local quantities are observed not as a function of $n$, but instead as $2^{n}$ (hence, geometrically instead of stroboscopically). New emergent conserved quantities, similar to Eq. (29) when $2^{n}$ is used instead of $n$, appear in this unitary dynamics but only when $n$ becomes sufficiently large, which then allow for the construction of a “geometric” GGE for the resulting steady state.

We begin by noting that the evolution operators for drives corresponding to $T-dT$ and $T+dT$ at a particular $k$ are respectively given by the matrices $U_{k}(T-dT)$ and $U_{k}(T+dT)$ given in Eqs. (55); for simplicity of notation, we will not write the subscript $k$ in this section henceforth. Furthermore, we will denote $U(T-dT)$ by $A_{0}$ and $U(T+dT)$ by $B_{0}$ here. Then the operators corresponding to the pair of drives $(\tau_{i},\tau_{i+1})$ = $-1,+1$ and $(\tau_{i},\tau_{i+1})$=$+1,-1$ are given by $A_{1}=B_{0}A_{0}$ and $B_{1}=A_{0}B_{0}$ respectively. We now focus on the case where $n=2^{m}$ where $m=0,1,2,\cdots$. It is then easy to see that the unitary dynamics can be entirely expressed in terms of the new matrices $A_{1}$ and $B_{1}$. We illustrate this in the second line of Eq. (59) for the case of $m=4$ (thus $n=16$).

[TABLE]

This motivates us to define matrices $A_{m}$ and $B_{m}$ recursively as

[TABLE]

for all $m\geq 0$. We can then show that the evolution operator after exactly $2^{m}$ drives is given by $A_{m}$. (Eq. (59) illustrates this for $m=4$, i.e., $n=16$). Note that this provides a very efficient method to calculate the dynamics after $n=2^{m}$ which only requires $\mathcal{O}(m)$ matrix multiplications at each $k$ instead of $\mathcal{O}(2^{m})$. We will show that remarkably $A_{m}=B_{m}$ when $m$ becomes sufficiently large and thus there is an emergent periodicity when the system is viewed geometrically (i.e., after every $n=2^{m}$).

We will now study how the matrices $A_{m},~{}B_{m}$ defined in Eqs. (60) behave in the limit $m\to\infty$. To this end, we parameterize these matrices as

[TABLE]

We also define the angle $\phi_{m}$ between the unit vectors ${\hat{a}}_{m}$ and ${\hat{b}}_{m}$, namely,

[TABLE]

with $0\leq\phi_{m}\leq\pi$. Using Eqs. (60-62), we find that

[TABLE]

From Eqs. (63), we can prove the following results. First, we find that $\alpha_{m}=\beta_{m}$ for all $m\geq 1$. We will therefore use only the variable $\alpha_{m}$ henceforth. Next, we find the recursion relations

[TABLE]

Since $\alpha_{m}=\beta_{m}$ for all $m\geq 1$, Eqs. (63) imply that

[TABLE]

This means that the unit vector which points in the direction of ${\hat{a}}_{m}~{}+~{}{\hat{b}}_{m}$ is the same for $m$ and $m+1$, up to a sign; hence, up to a sign, the direction of this unit vector does not change at all with $m$, all the way from $m=1$ to $\infty$. These arguments will break down if either ${\hat{a}}_{m}$ or ${\hat{b}}_{m}$ is ill-defined; according to Eq. (61), this can only happen if $\alpha_{m}$ or $\beta_{m}$ becomes exactly equal to 0 or $\pi$. However these cases require very special choices of $A_{0}$ and $B_{0}$ which form a set of measure zero; we will therefore ignore such special cases.

To summarize, we find that the dynamics given by Eq. (60) conserves two quantities at every iteration: $\alpha_{m}=\beta_{m}$, and the direction of the unit vector ${\hat{a}}_{m}+{\hat{b}}_{m}$ is conserved for all $m\geq 1$. The existence of these conserved quantities considerably simplifies the analysis as we will see below.

For sufficiently large $m$, we find numerically that $\phi_{m}\to 0$ although $\phi_{m}$ does not approach zero monotonically (see Fig. 9). We can understand this as follows. If $\phi_{m}\to 0$, Eq. (64) implies that $\alpha_{m+1}\to 2\alpha_{m}~{}{\rm mod}~{}\pi$. Assuming that $\alpha_{m}/\pi$ is irrational, we see that $\alpha_{m}$ will cover the region $[0,\pi]$ uniformly as $m\to\infty$. Next, assuming $\phi_{m},~{}\phi_{m+1}$ to be small, we obtain

[TABLE]

from the lowest order expansion of the terms appearing in Eq. (65). Since $\alpha_{m}$ covers $[0,\pi]$ uniformly, $|\tan(\alpha_{m})|$ varies all the way from zero to $\infty$; this explains why $\phi_{m}$ does not monotonically approach zero. In fact,

[TABLE]

and for large $N$, a uniform distribution of $\alpha_{m}$ implies that $\sum_{j=0}^{N-1}\log(|\tan\alpha_{m+j}|)\to(N/\pi)\int_{0}^{\pi}d\alpha\log(|\tan\alpha|)$ which is equal to zero. Thus we would expect that as $m\to\infty$, $\phi_{m}$ should stay at a constant value if we stopped at the first order expression in Eq. (67). We therefore have to go to the next order in the expansion of Eq. (65). We find that

[TABLE]

The fact that the last factor on the right hand side is always less than 1 explains why $\phi_{m}$ eventually goes to zero as $m\to\infty$; this explains the numerical behavior shown in Fig. 9.

As $m\to\infty$, the fact that $\phi_{m}\to 0$ means that ${\hat{a}}_{m}-{\hat{b}}_{m}\to 0$. The discussion following Eq. (66) then implies that ${\hat{a}}_{\infty}={\hat{b}}_{\infty}$ is equal, up to a sign, to the unit vector which points in the direction of ${\hat{a}}_{1}+{\hat{b}}_{1}$. Hence the value of ${\hat{a}}_{\infty}={\hat{b}}_{\infty}$ can be found right in the beginning when we calculate ${\hat{a}}_{1}$ and ${\hat{b}}_{1}$.

We now examine what happens if we view the system geometrically after $n=2^{m}$ drives where $m=0,1,2,\cdots$. In Fig. 9 (inset), we show the components of the unit vector ${\hat{a}}_{m}$ versus $m$ for $g_{i}=4,~{}g_{f}=2,~{}T=0.2,~{}dT=0.05$, and $k=\pi/2$. (In that figure, we have taken the third component of ${\hat{a}}_{m}$, called $a_{3m}$, to be non-negative for all values of $m$. According to Eq. (61), this can always be ensured, by flipping $\alpha_{m}\to-\alpha_{m}$ and ${\hat{a}}_{m}\to-{\hat{a}}_{m}$ if necessary.) We see that for $m$ up to about 1300, there are intermittent large fluctuations in both ${\hat{a}}_{m}$ and ${\hat{b}}_{m}$. As Eq. (67) shows, these fluctuations occur when $|\tan\alpha_{m}|$ becomes large, namely, when $\alpha_{m}$ comes close to $\pi/2$. For very large values of $m$, ${\hat{a}}_{m}$ and ${\hat{b}}_{m}$ settle down to the same value (Fig. 9, main panel). We therefore conclude that after an extraordinarily large number of drives given by $2^{m}$, where $m\gtrsim 1300$ (we have checked this for other values of $k$ also), ${\hat{a}}_{m}$ and ${\hat{b}}_{m}$ become identical and independent of $m$.

Since the unitary dynamics for $n=2^{m}$ can be represented by a single unitary matrix $A_{m}$ for the Thue-Morse sequence, this immediately shows that $\mathcal{J}_{k}$ defined as

[TABLE]

equals $\mathcal{J}_{k}\equiv\langle\psi_{k}(0)|\hat{a}_{\infty}\cdot\vec{\sigma}|\psi_{k}(0)\rangle$ and becomes independent of $m$ for sufficiently large $m$ at each momentum $k$. We should stress here, that unlike $\mathcal{I}_{k}$ defined in Eq. (29) for the $dT=0$ case which is conserved as a function of $n$, the quantities $\mathcal{J}_{k}$ are conserved geometrically, i.e., when $n=2^{m}$. Furthermore, one requires $m$ to be larger than a threshold value ($m\gtrsim 1300$ for the drive parameters shown in Fig. 9) for these conserved quantities to emerge at all momenta $k$, and these conserved quantities are not manifest in the unitary dynamics below this threshold value of $m$. We have also verified numerically that $\hat{a}_{m}$ and $\hat{b}_{m}$ become identical and independent of $m$ at large $m$ not only for $dT\ll T$ but for any value of $dT/T$.

The motion of the $N_{c}(\gg 1)$ consecutive momenta contained in any of the coarse-grained momentum cells on the unit sphere as a function of $n$ provides another way to geometrically see the emergent conserved quantities of this unitary dynamics (with the $\tau_{i}$’s chosen according to the Thue-Morse sequence). As shown in Fig. 10, when the system is viewed geometrically, i.e., with $n=2^{m}$, the $N_{c}$ points that belong to a momentum cell with average momentum $k_{c}$ initially perform a complicated dynamics on the unit sphere but eventually settle down to a very simple motion along a circle in the plane defined by $\hat{a}_{\infty}$ (at $k_{c}$) for large enough values of $m$ beyond a particular threshold.

Finally, we show in Fig. 11 the components of $\delta n_{k}={\hat{a}}_{k}(T,dT)-{\hat{a}}_{k}(T,0)$ versus $k$ when both ${\hat{a}}_{m}$ and ${\hat{b}}_{m}$ have settled down to the same value for a unitary process with the $\tau_{i}$’s chosen according to the Thue-Morse sequence for a given $T$ and $dT$. Since ${\hat{a}}_{k}(T,0)$ represents the unit vector that corresponds to the unitary matrix $U_{k}(T)$ for the perfectly periodic case ($dT=0$) given in Eqs. (28) and (LABEL:uk), we must have $\delta n_{k}\to 0$ as $dT\to 0$. In Fig. 11 (a), we show the variation of $\delta n_{k}$ versus $k$ for $dT=0.05,~{}T=0.2$, and the corresponding results for $dT=0.8,~{}T=0.2$ (thus, $dT>T$) in Fig. 11 (b). We see that indeed ${\hat{a}}_{\infty}(={\hat{b}}_{\infty})$ equals (up to a sign) the unit vector in the direction of ${\hat{a}}_{1}+{\hat{b}}_{1}$. The nature of the variation with $k$ is rather different for $dT\ll T$ (Fig. 11 (a)) and $dT\gg T$ (Fig. 11 (b)). Also, we note that $\delta n_{k}$ equals [math] for $k=0$ and $k=\pi$ as expected since $|a_{3}|=1$ and $a_{1}=a_{2}=0$ independently of the value of $dT$ at these two momenta.

Knowing the value of ${\hat{a}}_{\infty}={\hat{b}}_{\infty}$ as a function of $k$ (Fig. 11) allows us to construct the local description of the final non-equilibrium steady state. The construction is completely analogous to the perfectly periodic driven case ($dT=0$) where the relevant integrals of motion are now $\mathcal{J}_{k}$ (in Eq. (LABEL:conservedJk)) at each momentum instead of $\mathcal{I}_{k}$ (in Eq. (29)). The density matrix of the g-GGE then equals

[TABLE]

where the Lagrange multipliers $\lambda_{k}$ are fixed by the condition

[TABLE]

and the normalization constant $Z$ ensures that $\mathrm{Tr}[\rho_{\mathrm{g-GGE}}]=1$. The final steady state value of any local operator can then be evaluated by using this ensemble.

We show the results for one such local quantity $\langle\psi(n)|s^{x}|\psi(n)\rangle$ in Fig. 12. Fig. 12 (a) shows the convergence of this quantity to the final steady state value calculated using $\rho_{\mathrm{g-GGE}}$ for two different cases ($dT\ll T$ and $dT\gg T$) as a function of $n=2^{m}$; the convergence to the final g-GGE is extremely slow as we have already discussed in this section. Fig. 12 (b) shows the data for the case $dT\ll T$ as a function of $n$ from which it is clear that even after a large number of drives $n\sim 10^{5}$, the system does not seem to relax to the final non-equilibrium steady state in any simple manner if $dT\neq 0$ (whereas the relaxation of local quantities at $dT=0$ is much faster as is clear from the inset of Fig. 12 (a)), and the relaxation to the non-equilibrium steady state takes place over a much longer time scale, as we see in Fig. 12 (a). We note here that this extremely slow relaxation of local quantities to their final steady state appears to be quite insensitive to the value of $dT$ as long as it is non-zero. Fig. 12 (b) also shows the perfect agreement of the data calculated at $n=2^{m}$ using the unitary matrix $A_{m}$ at each $k$ to the results of the more direct calculation at each $n$ (using Eq. (15)) without using the recursive structure of the Thue-Morse sequence.

VI Conclusions and outlook

In this work, we have considered a prototypical integrable model, the one-dimensional transverse field Ising model, which is continually driven with a time-dependent Hamiltonian that results in a unitary dynamics. Since the model is translation invariant, the modes can be labeled by their momentum $k$. For a subsystem whose size is much smaller than that of the total system (shown schematically in Fig. 1), the reduced density matrix of the subsystem is controlled by certain coarse-grained quantities in momentum space as defined in Eq. (24). Even though the corresponding quantities for a single momentum continually fluctuate as $n\to\infty$ and do not have a well-defined large-$n$ limit, these coarse-grained quantities have a well-defined $n\to\infty$ limit if a non-equilibrium steady state exists. Furthermore, we define a probability distribution on the unit sphere at each average momentum $k_{c}$ by employing a pseudospin representation at each $k$ and using the large number of momenta lying in a small cell in momentum space centered around a momentum $k_{c}$ during the coarse-graining procedure. This probability distribution at each $k_{c}$ changes in an irreversible manner to the final $n\to\infty$ distribution as a function of $n$. At large $n$, this distribution generated from the points in a cell (centered at $k_{c}$) at that $n$ has the property stated in Eq. (30) that it equals the probability distribution generated from the motion of the pseudospin at a single momentum $k_{c}$ (thus, with no coarse-graining in momentum) for a large number of drives.

We first recapitulate the results known for the well-studied case in which the system is driven perfectly periodically with a time period $T$. Since the model is completely integrable with an infinite number of conserved quantities $\mathcal{I}_{k}$ (defined in Eq. (29)), the system does not locally heat up to infinite temperature, and instead this driving leads to the periodic generalized Gibbs ensemble. We also calculate the behavior of the coarse-grained quantities in momentum space and the probability distributions on the unit sphere as a function of $k_{c}$ which immediately shows that the non-equilibrium steady state is not an infinite temperature ensemble.

We then consider two new driving protocols both of which involve driving with two possible periods $T+dT$ and $T-dT$, such that the unitary time evolution operators corresponding to the two periods do not commute with each other; these driving protocols can be thought of as perturbations of the perfectly periodic driving if $dT\ll T$. In the first driving protocol, the successive drive periods are randomly chosen to be $T+dT$ and $T-dT$ with probability 1/2 each. In this case, we find that after a large number of drives the system eventually approaches an infinite temperature ensemble. Remarkably, this ensemble also has the property that the probability distribution of a coarse-grained quantity averaged over many momenta in a cell is equal to the probability distribution of the same quantity generated by many drives for a single momentum. The second driving protocol that we consider is one in which the time periods $T+dT$ and $T-dT$ are chosen according to the Thue-Morse sequence. This sequence is neither periodic nor random, but it has the property of being scale-invariant if the number of drives is doubled. In this case, we find that the system approaches a non-equilibrium steady state with conserved quantities $\mathcal{J}_{k}$ (defined in Eq. (LABEL:conservedJk)) which are smooth functions of $k$ if we view it after $n=2^{m}$ drives rather than after $n$ drives. We call this steady state a geometric generalized Gibbs ensemble; to the best of our knowledge, this is the first known example of such an ensemble. We note however that this ensemble emerges only when $m$ is large (of the order of 1300 for one particular set of parameters) which implies that the number of drives $n$ required is astronomically large.

An interesting problem for future studies may be to understand how the local properties relax in time to the final geometric generalized Gibbs ensemble that we have discussed here. Such a study would lead to a deeper understanding of the fixed point of a sequence of $SU(2)$ matrices which is generated from two arbitrary initial matrices $A_{0}$ and $B_{0}$ (which do not commute with each other) by the recursion relations given in Eq. (60). The numerical results that we have presented in Sec. V provide only a glimpse of this general problem which may possibly involve a very rich mathematical structure.

More generally, it will be interesting to understand the conditions which determine whether a continually driven integrable system locally heats up to an infinite temperature ensemble or not ongoing . For the latter case, it will be instructive to find other examples of quasiperiodic drive protocols Sutherland1986 , apart from the one we have pointed out here, where well-defined non-equilibrium steady states may emerge without the restriction of a time-periodic drive.

Acknowledgments

A.S. is grateful to Kedar Damle and Pushan Majumdar for useful discussions. The work of A.S. is partly supported through the Partner Group program between the Indian Association for the Cultivation of Science (Kolkata) and the Max Planck Institute for the Physics of Complex Systems (Dresden). D.S. thanks Department of Science and Technology, India for Project No. SR/S2/JCB-44/2010 for financial support.