Exact results in Floquet coin toss for driven integrable models

TL;DR

This paper derives exact analytical results for the residual energy in a driven integrable quantum model under random binary (coin-toss) driving, showing bounded heating and universal behavior in the asymptotic limit.

Contribution

It introduces an exact formalism using a disorder operator to analytically compute the residual energy in a randomly driven integrable system, revealing universal asymptotic behavior.

Findings

Residual energy saturates to a finite value with random driving.

The formalism simplifies numerical calculations significantly.

Asymptotic residual energy is universal, independent of protocol bias.

Abstract

We study an integrable Hamiltonian reducible to free fermions which is subjected to an imperfect periodic driving with the amplitude of driving (or kicking) randomly chosen from a binary distribution like a coin-toss problem. The randomness present in the driving protocol destabilises the periodic steady state, reached in the limit of perfectly periodic driving, leading to a monotonic rise of the stroboscopic residual energy with the number of periods ($N$). We establish that a minimal deviation from the perfectly periodic driving would always result in a {\it bounded} heating up of the system with $N$ to an asymptotic finite value. Remarkably, exploiting the completely uncorrelated nature of the randomness and the knowledge of the stroboscopic Floquet operator in the perfectly periodic situation, we provide an exact analytical formalism to derive the disorder averaged expectation value of the residual energy through a {\it disorder operator}. This formalism not only leads to an immense numerical simplification, but also enables us to derive an exact analytical form for the residual energy in the asymptotic limit which is universal, i.e, independent of the bias of coin-toss and the protocol chosen. Furthermore, this formalism clearly establishes the nature of the monotonic growth of the residual energy at intermediate $N$ while clearly revealing the possible non-universal behaviour of the same.

Figures

a and 1b of the main text for different protocols and for different values of $p$ for a given protocol.

Equations

$$H=-\sum _{n=1}^L{\tau }_n^x{\tau }_{n+1}^x-h\sum _{n=1}^L{\tau }_n^z.$$

$$f(t)$$

$$=\sum _{n=1}^N{g}_n\alpha \sin \left(\frac{2\pi t}{T}\right)\text{for the sinusoidal drive},$$

$${\mathcal{F}}_k(T)=\exp (-i\alpha {\sigma }_z)\exp (-i{H}_k^{0}T)$$

$$\displaystyle\varepsilon_{res}(NT)=\frac{1}{2\pi}\int dk\biggl{[}\sum_{\alpha=\pm}|r_{k}^{\alpha}|^{2}\langle\phi_{k}^{\alpha}\rvert H_{k}^{0}\lvert\phi_{k}^{\alpha}\rangle$$

$$|{\psi }_k(NT)⟩={\mathbb{U}}_k(g_N){\mathbb{U}}_k(g_{N-1}).......{\mathbb{U}}_k(g_2){\mathbb{U}}_k(g_1)|{\psi }_k(0)⟩$$

$${\mathbb{U}}_k(g_n)=\{\begin{array}{ll}\text{boldmath}{\mathcal{F}}_k(T),& \text{if gn=1}.\\ U_k^0(T),& \text{if gn=0}.\end{array}$$

$$e_k(NT)$$

$$\times H_k^0{\mathbb{U}}_k(g_N){\mathbb{U}}_k(g_{N-1})........{\mathbb{U}}_k(g_2){\mathbb{U}}_k(g_1)|{\psi }_k(0)⟩$$

$$⟨e_k(NT)⟩=/\sum ⟨{\psi }_k(0)|j_k^{(0)}⟩⟨j_k^{(N)}|H_k^0|i_k^{(N)}⟩⟨i_k^{(0)}|{\psi }_k(0)⟩\left[\prod _{m=1}^N\left(\sum _{g_m=1,0}P(g_m)⟨j_k^{(m-1)}|{\mathbb{U}}_k^{†}(g_m)|j_k^{(m)}⟩⟨i_k^{(m)}|{\mathbb{U}}_k(g_m)|i_k^{(m-1)}⟩\right)\right]$$

$$⟨e_k(NT)⟩=/\sum ⟨{\psi }_k(0)|j_k^{(0)}⟩⟨j_k^{(N)}|H_k^0|i_k^{(N)}⟩⟨i_k^{(0)}|{\psi }_k(0)⟩$$

$${\mathbb{D}}_k^{j^{(m-1)},j^{(m)},i^{(m-1)},i^{(m)}}\equiv \left(p{e}^{iT\left({ϵ}_k^{j_{m-1}}-{ϵ}_k^{i_{m-1}}\right)}{\delta }_{j_k^{m-1},j_k^m}{\delta }_{i_k^m,i_k^{m-1}}+(1-p)⟨j_k^{(m-1)}|{\mathbb{U}}_k^{0†}|j_k^{(m)}⟩⟨i_k^{(m)}|{\mathbb{U}}_k^0|i_k^{(m-1)}⟩\right)$$

$$\underset{N\to \infty }{\lim }⟨{\epsilon }_{res}(NT)⟩=\frac{1}{\pi }\left[e_{k=\pi }^g(0)-{\int }_0^{\pi }{e}_k^g(0)dk\right].$$

$$D^{j^{(m-1)},j^{(m)},i^{(m-1)},i^{(m)}}(k)\equiv \left(p{e}^{iT\left({ϵ}_k^{j_{m-1}}-{ϵ}_k^{i_{m-1}}\right)}{\delta }_{j_k^{m-1},j_k^m}{\delta }_{i_k^{m-1},i_k^m}+(1-p)⟨j_k^{(m-1)}|{\mathbb{U}}_k^{0†}|j_k^{(m)}⟩⟨i_k^{(m)}|{\mathbb{U}}_k^0|i_k^{(m-1)}⟩\right)$$

$$U_{k}^{0}(T)\doteq\left(\begin{array}[]{cc}u_{11}(k)-iu_{12}(k)&-iu_{21}(k)\\ -iu_{21}(k)&u_{11}(k)+iu_{12}(k)\\ \end{array}\right)$$

$$u_{11}(k)$$

$$u_{12}(k)$$

$$u_{21}(k)$$

$$\bra{j_{k}^{m}}U_{k}^{0}(T)\ket{j_{k}^{n}}\doteq\left(\begin{array}[]{cc}W_{11}(k)&W_{12}(k)\\ -W_{12}^{*}(k)&W_{11}^{*}(k)\\ \end{array}\right)$$

$$D(k)$$

$$\doteq$$

$$D(k=0,\pi)=\left(\begin{array}[]{cccc}1&0&0&0\\ 0&p\exp{[i\Delta\phi_{k=0,\pi}T]}+(1-p)\exp{[i\Delta e_{k=0,\pi}T]}&0&0\\ 0&0&p\exp{[-i\Delta\phi_{k=0,\pi}+(1-p)\exp{[-i\Delta e_{k=0,\pi}T]}T]}&0\\ 0&0&0&1\\ \end{array}\right)$$

$$\lim_{N\to\infty}D^{N}(k=0)=\left(\begin{array}[]{cccc}1&0&0&0\\ 0&1&0&0\\ 0&0&1&0\\ 0&0&0&1\\ \end{array}\right);\lim_{N\to\infty}D^{N}(k=\pi)=\left(\begin{array}[]{cccc}1&0&0&0\\ 0&0&0&0\\ 0&0&0&0\\ 0&0&0&1\\ \end{array}\right)$$

$$\lim_{N\to\infty}D^{N}(k=0,\pi)=\left(\begin{array}[]{cccc}1&0&0&0\\ 0&0&0&0\\ 0&0&0&0\\ 0&0&0&1\\ \end{array}\right)$$

$$\left[\lambda (k)-s(\lambda ,k)\right]f(\lambda ,k)=0$$

$$\lim_{N\to\infty}D^{N}(k\neq 0,\pi)=\frac{1}{2}\left(\begin{array}[]{cccc}1&0&0&1\\ 0&0&0&0\\ 0&0&0&0\\ 1&0&0&1\\ \end{array}\right)$$

$$\underset{N\to \infty }{\lim }⟨e_{k=0}(NT)⟩$$

$$\underset{N\to \infty }{\lim }⟨e_{k=\pi }(NT)⟩$$

$$\underset{N\to \infty }{\lim }⟨e_{k\ne 0,\pi }(NT)⟩$$

$${\epsilon }_{res}(NT)=\frac{1}{\pi }{\int }_0^{\pi }dk\left[e_k(NT)-e_k^g(0)\right]$$

$$\underset{N\to \infty }{\lim }{\epsilon }_{res}(NT)=\frac{1}{\pi }{\int }_0^{\pi }dk\left[\left\{\sum _{\alpha =1,2}{\left|⟨j_k^{\alpha }|{\psi }_k^0⟩\right|}^2⟨j_k^{\alpha }\left|H_k^0\right|j_k^{\alpha }⟩\right\}-e_k^g(0)\right]$$

$$\underset{N\to \infty }{\lim }⟨{\epsilon }_{res}(NT)⟩=\frac{1}{\pi }\underset{N\to \infty }{\lim }⟨e_{k=\pi }(N)⟩-\frac{1}{\pi }{\int }_0^{\pi }{e}_k^g(0)dk=\frac{1}{\pi }\left[e_{k=\pi }^g(0)-{\int }_0^{\pi }{e}_k^g(0)dk\right]$$

$$a_2$$

$$a_1$$

Full text

The study of non-equilibrium dynamics of closed quantum systems is an exciting as well as a challenging area of recent research both from experimental bloch08 ; lewenstein12 ; jotzu14 ; greiner02 ; kinoshita06 ; gring12 ; trotzky12 ; cheneau12 ; fausti11 ; rechtsman13 ; schreiber15 and theoretical perspectives calabrese06 ; rigol08 ; oka09 ; kitagawa10 ; lindner11 ; bermudez09 ; patel13 ; thakurathi13 ; pal10 ; nandkishore15 ; heyl13 ; budich15 ; sharma16 . One of the prominent areas in this regard is periodically driven closed quantum systems (for review see, dziarmaga10 ; polkovnikov11 ; dutta15 ; eisert15 ; alessio16 ; stat ; bukov16 ) which have an illustrious history dating back to the analysis of the famous Kapitza pendulum kapitza51 and the kicked-rotor model casati79 . The recent interest in periodically driven systems are many fold: e.g., Floquet engineering of materials in their non-trivial phases such as the Floquet graphene oka09 ; kitagawa10 and topological insulators lindner11 (see also cayssol13 ), dynamical generation of edge Majorana thakurathi13 and non-equilibrium phase transitions bastidas11 like recently proposed time crystals else16 ; khemani16 . These studies have received a tremendous boost following experimental studies on light-induced non-equilibrium superconducting and topological systems fausti11 ; rechtsman13 and possibility of realising time crystals zhang17 ; choi17 . The other relevant question deals with fundamental statistical aspects, namely the thermalisation of a closed quantum system under a periodic driving russomanno12 and the possibility of the many-body localisation ponte15 .

Periodically driven closed quantum systems, from a statistical viewpoint, are being extensively studied in the context of defect and residual energy generation mukherjee08 ; russomanno12 , dynamical freezing das10 , many-body energy localization alessio13 dynamical localisationnag14 ; agarwala16 , dynamical fidelity sharma14 , work-statistics anirban_dutta14 ; russomanno15 , and as well as in the context of entanglement entropy russomanno16 and associated dynamical phase transitions sen16 . For periodically driven closed integrable systems, reducible to free fermions, it is usually believed that the system reaches a periodic steady state in the asymptotic limit of driving and hence stops absorbing energy russomanno12 ; gritsev17 : the resulting steady state can be viewed as a periodic Gibbs ensemble with an extensive number of conserved quantities lazarides14 .

However, for a non-integrable model alessio14 or for an aperoidic driving of an integrable model, the system is expected to absorb energy indefinitely; also there exists a possibility of a geometrical generalised Gibbs ensemble in some special situations nandy17 . However, for a driven non-integrable system a MBL state may also arise ponte15 and also a MBL state may get delocalised under a periodic driving lazarides15 . While for a periodically driven non-integrable system it is challenging to prevent the system from heating up, the noninteracting case is fundamentally stable and does not suffer from a “heat death problem”. Considering an integrable model that is reducible to free fermions, we show below that even in this simplest situation, the slightest deviation from the perfect periodicity, which is experimentally inevitable, would always result in heating up of the system with the stroboscopic time to a finite bounded asymptotic value. This bounded nature is in sharp contrast to the situation mentioned in Ref. ott84 ; guarneri84 , where it has been proved, based on an approach of the spectral analysis of a non-random operator guarneri84 , that the expectation value of the kinetic energy operator of a noisy $\delta$-perturbed quantum rotator is unbounded in time.

It should be noted at the outset that the fate of other types of integrable systems (e.g., Bethe-integrable ones) even under a perfectly periodic drives is still an open issue bukov .

Furthermore, this problem of aperiodic driving from a broad scenario raises a plethora of pertinent questions: for example, what would be the fate of an emergent topological phase under such a temporal noise leon13 ? Can the aperiodicity be visualised as an outcome of coupling the system to a bath iwahori16 ? What would happen to the stroboscopic entanglement entropy in the asymptotic limit russomanno16 ? Can the deviation from periodicity in driving also induce localization-delocalization transitions in MBL systems lazarides15 ? Can the problem be connected to a quantum random walk problem and corresponding search algorithms kempe03 ; oka03 ? Although our attention here is limited to the problem of heating up only, we believe that the framework we develop would definitely provide the key foundation to address most of the questions that have been mentioned above.

In this paper, we consider a closed integrable quantum system undergoing an imperfect periodic dynamics. The imperfection or disorder manifests itself in the amplitude of the periodic drive which assumes binary values chosen from a binomial distribution resembling a series of biased coin toss events. The combination of a periodic driving with such a disordered amplitude results in the so-called “Floquet coin toss problem”. The aim of the paper is to provide an exact analytic framework to explore the statistical properties of such a non-equilibrium system observed at $N$-th stroboscopic interval determined by the inverse of the frequency ($\omega$) of the perfectly periodic drive. Our study indeed confirms that the system, even though integrable, never reaches a periodic steady state and rather keeps on absorbing energy till the asymptotic limit.

We establish the above claim considering the one dimensional transverse Ising Model, described by the Hamiltonian,

[TABLE]

where $h$ is the transverse field and $\tau_{n}^{i}$ $\{i=x,y,z\}$ are the Pauli spin matrices at $n^{th}$ site. The Hamiltonian can be decoupled in to $2\times 2$ problems for each Fourier mode via a Jordan-Wigner mapping, such that $H=\sum_{k}{H_{k}}$ with $H_{k}=(h-\mathrm{cos}k)\sigma_{z}+\mathrm{sin}k\sigma_{x}$, where $\sigma$’s are again Pauli matrices. We here use the anti-periodic boundary condition for even $L$ so that $k=\frac{2m\pi}{L}$ with $m=-\frac{L-1}{2},...,-\frac{1}{2},\frac{1}{2},...,\frac{L-1}{2}$.

In our present work we study the effect of aperiodic temporal variation of the external transverse field $h$, considering two different types of driving protocols, namely, the delta kicks and the sinusoidal driving incorporating a binary disorder in the amplitude of driving. In short, we have $h(t)=1+f(t)$ with

[TABLE]

where $\alpha$ is the amplitude of driving and $T$ is the time period. The random variable $g_{n}$ ( where $n$ refers to the $n$-th stroboscopic period) takes the value either $1$ with probability $p$ or [math] with probability $(1-p)$ chosen from a Binomial distribution. Evidently $g_{n}=0$, corresponds to free evolution in the $n^{th}$ time period within the time interval $(n-1)T$ to $nT$ while $g_{n}=1$ corresponds to the periodic perturbation in the form of a kick or sinusoidal driving.

To illuminate the underlying Floquet theory, let us first consider the case of fully periodic situation $(p=1)$ choosing the initial state $|\psi_{k}(0)\rangle$ as the ground-state of the free Hamiltonian in our case; for each $k$ mode, we then have $H_{k}(t+T)=H_{k}(t)$. Using the Floquet formalism, one can define a Floquet evolution operator \boldmath{\cal{F}}_{k}(T)=\cal{T}$$\exp\left(-i\int_{0}^{T}H_{k}(t)dt\right), where $\cal{T}$ denotes the time ordering operator and the solutions of the Schrödinger equation for a time periodic Hamiltonian can be written as, $\lvert\psi_{k}^{(j)}(t)\rangle=\exp(-i\epsilon_{k}^{(j)}T)\lvert\phi_{k}^{(j)}(t)\rangle$. The states $\lvert\phi_{k}^{(j)}(t)\rangle$, the so called Floquet modes satisfying the condition $\lvert\phi_{k}^{(j)}(t+T)\rangle=\lvert\phi_{k}^{(j)}(t)\rangle$ and the real quantities $\epsilon_{k}^{(j)}$ are known as Floquet quasi-energies. In the case of the $\delta$-function kick, $\boldmath{\cal{F}}_{k}$ can be exactly written, in the form of,

[TABLE]

consisting of two pieces; the first one corresponds to the $\delta$-kick at time $t=T$ while the second part represents the free evolution generated by the time independent Hamiltonian $H^{0}_{k}=(1-\cos k)\sigma_{z}+(\sin k)\sigma_{x}$ from time [math] to $T$. However for the sinusoidal driving $\mathbb{\cal{F}}_{k}(T)$ can be numerically diagonalized to obtain the quasi-energies $\epsilon_{k}^{\pm}$ and the corresponding Floquet modes $\rvert\phi_{k}^{\pm}\rangle$. Focussing on the mode $k$, after a time $t=NT$ we have $\lvert\psi_{k}(NT)\rangle=\sum_{j=\pm}r_{k}^{(j)}e^{-i\epsilon_{k}^{(j)}NT}\lvert\phi_{k}^{(j)}(t)\rangle$, where $r_{k}^{\pm}=\langle\phi_{k}^{\pm}\rvert\psi_{k}(0)\rangle$ and hence the residual energy $\varepsilon_{res}(NT)=\frac{1}{L}\sum_{k}\left(e_{k}(NT)-e^{g}_{k}(0)\right)$ with $e_{k}(NT)=\langle\psi_{k}(NT)\rvert H_{k}^{0}\lvert\psi_{k}(NT)\rangle$ and $e^{g}_{k}(0)=\langle\psi_{k}(0)\rvert H_{k}^{0}\lvert\psi_{k}(0)\rangle$. In the thermodynamic limit of $L\rightarrow\infty$, we have

[TABLE]

In the limit $N\rightarrow\infty$, due to the Riemann-Lebesgue lemma the rapidly oscillating off-diagonal terms in Eq. (5) drop off upon integration over all $k$ modes leading to a steady state expression for $\varepsilon_{res}$ russomanno12 .

Deviating from the completely periodic case and considering the situation $0<p<1$, so that we have a probability $(1-p)$ of missing a kick (or a cycle of sinusoidal drive is absent) in every complete period, one can now write the corresponding evolved state after $N$ complete periods as

[TABLE]

with the generic evolution operator given by,

[TABLE]

where $\boldmath{\cal{F}}_{k}(T)$ is the usual Floquet operator and $U^{0}_{k}(T)=\exp(-iH^{0}_{k}T)$ is the time evolution operator for the free Hamiltonian $H^{0}_{k}$. One then readily finds:

[TABLE]

The numerical calculation in Eq. (8) thus involves two steps: (i) the multiplication of $(2N+1)$ matrices corresponding to $N$ complete periods for a given disorder configuration as shown in Eq. (8). (ii) Configuration averaging over the disorder for a fixed value of $N$. Let us refer to the Fig. S2, where numerically obtained residual energy (RE) is plotted as a function of number of complete periods $N$ choosing a high frequency limit ($\omega=2\pi/T>4$) in which the fully periodic situation leads to a periodic steady state both for the $\delta$-kicks and sinusoidal variation and the system synchronises with the external driving. Observing the monotonic rise of RE with $N$ for a given $p\neq 0,1$, we conclude that for any non-zero value of $p$, the periodic steady state gets destabilised and system keeps absorbing energy. We further note, that for a given $N$, the behavior of the RE cannot monotonically rise with $p$ since the RE gets constrained by the fully periodic situation around $p=1$ and the no rise situation around $p=0$. The value of $p$ for which the rate of rise of RE with $N$ will be maximum depends on both $\alpha$ and $\omega$ of the imperfect drive (see supplementary material (SM))  Furthermore, a slightest inclusion of aperiodicity in the drive leads to the difference in the RE to rise linearly with $N$ (up to an appropriate value of $N$) both in the limit $p\to 0$ from the undriven situation or $(1-p)\to 0$ from the perfectly periodic driven situation, as has been elaborated explicitly in the SM.

Having provided the numerical results, we shall proceed to set up the corresponding analytical framework within the space spanned by the complete set of Floquet basis $\{\lvert j_{k}^{\pm}\rangle\}$ states. Introducing $2(N+1)$ identity operators in terms of the Floquet basis states $\sum_{j^{(m)}}\lvert j^{(m)}_{k}\rangle\langle j^{(m)}_{k}\rvert=\hat{1}$ in Eq. (8), where $j^{(m)}_{k}$ can take two possible values corresponding to two quasi-states of the $2\times 2$ Floquet Hamiltonian $\boldmath{\cal{F}}_{k}(T)$ for each mode $k$ and performing the average over disordered configurations and finally upon reorganisation, we find,

[TABLE]

where $\not{\sum}\equiv\sum_{\begin{subarray}{c}j^{(0)},j^{(1)},....,j^{(N)}\\ i^{(0)},i^{(1)},....,i^{(N)}\end{subarray}}$. The uncorrelated nature of $g_{m}$s enables us to perform the configuration average by separately averaging over for each $g_{m}$. Recalling Eq. (7), and the fact that $\boldmath{\cal{F}}_{k}(T)|j^{\pm}_{k}\rangle=\exp(-i\epsilon^{\pm}_{k}T)|j^{\pm}_{k}\rangle$ and $P(g_{m})$ is the probability of being perfectly driven [$P(g_{m}=1)=p$] and free evolution [$P(g_{m}=0)=(1-p)$], respectively, leads us to:

[TABLE]

where the matrix element of ($4\times 4$) matrix $\mathbb{D}_{k}$ is given by,

[TABLE]

It should be noted that in Eq. (10), the $N$ in $\mathbb{D}_{k}^{N}$ is not a label but the matrix $\mathbb{D}_{k}$ that we have defined in Eq. (S1), raised to the power $N$. The above exercise naturally leads to the emergence of $4\times 4$ disorder matrix $\mathbb{D}$ for a imperfectly driven $2\times 2$ system. Given the amplitude, frequency, dimensionality, the form of ${\mathcal{F}}_{k}(T)$ in every stroboscopic intervals and the knowledge of disorder encoded in the probability $p$ of driving, every element of $\mathbb{D}$-matrix can be exactly calculated as shown in the SM where we probe the analytical structure of the disorder matrix emphasising both on the limit $N\to\infty$ and intermediate $N$. The intriguing feature is that one of the eigenvalues of the (non-unitary) disorder matrix becomes unity while the other real eigenvalue and also the modulus of complex eigenvalues are less than unity and thus become vanishingly small in the diagoanlised form of ${\mathbf{D}}^{N}$ in the limit $N\to\infty$. This remarkable property of the $\mathbf{D}$-matrix immediately renders an exact analytical form of the residual energy in the asymptotic limit:

[TABLE]

The asymptotic value obtained in Eq. (S31) is clearly independent of the driving strength, frequency and the protocol implemented, and hence is evidently universal. It is also finite and bounded in contrast to the kicked rotor case in Ref. guarneri84 . Further, the RE as obtained from the exact form given in Eq. (10) confirms the numerically obtained values presented in Fig. S2 for all values of $N$.

For small values of $N$, when the transients have not yet decayed, the slope of the residual energy curves for the two different protocols (as shown in Fig. S2) are different. This difference, however disappears in the long time limit when the initial transients die off. Hence, the curves for residual energy for both the protocols exhibit a universal behaviour at large but intermediate $N$ as depicted in Fig. 2. This coarse-grained (in time) view of $\varepsilon_{res}$ establishes the similar behaviour for two protocols in such large time scale; this can interestingly be understood by noting that at large $N$, the phase of the complex eigenvalues, $\exp[\pm iN\phi(k,\alpha,\omega)]$, of the disorder operator, which are responsible for the initial interferences, oscillate rapidly and becomes vanishingly small for large intermediate $N$ as illustrated in the SM. However, the other real eigenvalue $r=r(k,\alpha,\omega)$, with modulus less than unity survives for large intermediate $N$ and goes to zero only when $N\to\infty$. Thus, out of the four eigenvalues of the disorder matrix only two eigenvalues, $1$ and $r$ (see SM) raised to the power $N$, are non-vanishing and compete against each other. While the unit eigenvalue forces the system to achieve an universal value at $N\to\infty$, $r^{N}$ at large intermediate times is solely responsible for the monotonic rise of the residual energies for any protocol chosen. Of course, $r^{N}\to 0$ eventually at $N\to\infty$. Thus, the monotonic growth of $\varepsilon_{res}$ with increasing $N$, is ensured by the gradual decay of the (positive) $r^{N}$. We note that although at finite $N$ the curves at long time scales are similar, their slopes as they tend towards the asymptotic limit are different for the two protocols. This can be ascribed to the fact that the eigenvalue $r=r(k,\alpha,\omega)$ is not only dependent on the driving strength and frequency but its functional form is dictated by the protocol chosen and hence, its value is non-universal. Thus, we have constructed a complete analytic framework to deal with any “Floquet coin toss” problem in general.

The numerical simulations addressing similar problems involving multiplication of $N$ unitary matrices, followed by averaging over a large number ($N_{c}$) of disorder configurations; as $N$ and $N_{c}\to\infty$, requires huge computational time while being susceptible to numerical inaccuracies. On the contrary, the non-perturbative $\mathbf{D}$-matrix, overcomes such numerical challenges by calculating the exact analytical behavior of $\varepsilon_{res}$ at large $N$ and especially in the asymptotic limit which is otherwise extremely expensive numerically comment .

The sole aim of this paper has been to develop the disorder matrix formalism, given the complete stroboscopic information for a perfectly periodic case, to establish the bounded growth of the residual energy to a asymptotic value following an aperiodic driving keeping free fermion reducible integrable quantum systems in mind. However, we believe it can be generalised to situations where the disorder matrix is higher dimensional. In such scenarios, the asymptotic value of the residual energy and other such operators may still be obtained from our analysis. For non-integrable situations, given uncorrelated binary distribution of disorder, if complete stroboscopic information about the Floquet evolution operator for the $p=1$ situation is completely known, the disorder matrix formalism is expected to hold. We note that the aperiodic kicking situation proposed here has already been experimentally realised for a single rotor by Sarkar $et~{}al.$ sarkar17 ; similar experimental studies for aperiodically driven many body systems are indeed possible. Given a rare possibility of analytical approach to explore a temporally disordered situation and given the wide scope of the validity of our results, we believe that our approach is going to provide a new avenue to a plethora of similar studies. We conclude with a note that while spatial disorder leads to Anderson localisation, the temporal disorder which we study here remarkably leads to a delocalisation (in the Floquet space) with a bound; the importance of our work in this regard has already been noted in Ref. quelle17 .

I The analytical structure of the disorder matrix and the growth of the residual energy

In this supplementary material, we shall probe the analytical structure of the disorder matrix and thus establish the growth of the residual energy as a function of the stroboscopic period $N$, especially emphasising the limit $N\to\infty$ and large intermediate $N$. What is intriguing, as we shall show below, is that one of the eigenvalues of the (non-unitary) disorder matrix is unity while the other eigenvalues ( real and the modulus of the two complex eigenvalues) are less than unity and thus become vanishingly small in the diagonalized form of $D^{N}(k)$ in the limit $N\to\infty$. This remarkable property of the $D$-matrix results in a very simplified form of the disorder matrix which immediately renders an exact analytical form of the residual energy in the asymptotic limit which is independent of the value of the bias $p$ and the protocol used. For finite $N$ on the other hand, the growth of the residual energy shows a non-universal behaviour, i.e., is bias and protocol dependent as in this limit its behaviour is dictated by both the real eigenvalues.

I.1 The exact analytical structure of the disorder matrix

To present the exact analytical structure of the $D$-matrix, let us consider the case of a 1-D Ising model subjected to an imperfect driving protocol referring to Eq. (11) of the main text:

[TABLE]

We recall that $\epsilon_{k}^{j}$ are the quasi-energies and $\ket{j_{k}^{1,2}}$ are corresponding eigenvectors of the $2\times 2$ Floquet Hamiltonian. On the other hand, $U_{k}^{0}(T)$ represents the time evolution operator of the undriven Ising Hamiltonian measured at stroboscopic instant $T$. The matrix representation of $U_{k}^{0}(T)$ is as follows:

[TABLE]

where

[TABLE]

Due to the unitarity, the matrix $U_{k}^{0}$ in the basis of the eigenvectors $\ket{j_{k}^{1,2}}$ can be cast in a general form:

[TABLE]

with $|W_{11}(k)|^{2}+|W_{12}(k)|^{2}=1$. Thus, a general form of the $4\times 4$ disorder matrix can easily be constructed for the above situation to obtain,

[TABLE]

where $\Delta\phi_{k}=\epsilon_{k}^{+}-\epsilon_{k}^{-}$, $r_{i}(k)$ and $c_{i}(k)$ denotes the real and the complex elements of the matrix, respectively. Let us recall that the disorder matrix $D(k)$ emerges due to disorder (classical) averaging over infinite number of configurations and hence evidently is a non-Unitary matrix with the absolute values of all its elements being less than unity. Further, in the case of perfectly periodic driving $p=1$, the off-diagonal terms of the $D$-matrix in (S14) vanishes rendering it in a diagonal form.

I.2 Eigenvalues and eigenvectors of the $D$-matrix in the limit $N\to\infty$

To continue our analysis further, let us first focus on the modes, $k=0$ and $\pi$; the off-diagonal terms of the Hamiltonian $H^{0}_{k}=(1-\cos k)\sigma_{z}+(\sin k)\sigma_{x}$ vanishes for these modes and hence these modes do not evolve with time. The $2\times 2$ matrix in (S6) becomes identity matrix for $k=0$ and a diagonal matrix with two diagonal elements $\exp(\mp 2iT)$ for $k=\pi$. Hence from (S14), we immediately find

[TABLE]

where $\Delta e_{k}=e_{k}^{e}(0)-e_{k}^{g}(0)$ is the energy gap of the Hamiltonian $H^{0}_{k}$. For the sinusoidal drive $\Delta e_{k}=\Delta\phi_{k}$, whereas this is not true for the $\delta$-kicked situation as the integral of $\delta(t-T)$ is not zero over a complete period. Since both the matrices $D(k=0,\pi)$ is diagonal, it is straightforward to estimate $\lim_{N\to\infty}D^{N}$ for these modes. For the sinusoidal drive,

[TABLE]

Whereas for the $\delta$-kicked situation,

[TABLE]

For other modes ($k\neq 0,\pi$), the matrix $D(k)$ indeed has off-diagonal terms and hence to compute the four eigenvalues $\lambda(k)$, one needs to analyse the equation:

[TABLE]

where $s(\lambda,k)=\left\{p+(1-p)\left(\left|W_{11}(k)\right|^{2}+\left|W_{12}(k)\right|^{2}\right)\right\}=1$ as $\left(\left|W_{11}(k)\right|^{2}+\left|W_{12}(k)\right|^{2}\right)=1$, and $f(k,\lambda)$ is a third degree polynomial with all real coefficients. Thus, it is obvious that always one of the eigenvalues $\lambda(k)=\lambda_{1}(k)=1$ with normalised eigenvector, v_{1}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}1\\ 0\\ 0\\ 1\\ \end{array}\right) for $k\neq 0,\pi$. While one eigenvalue sticks to unity, it is straightforward to argue that other eigenvalues (one real and the other two complex conjugates of each other) will have a value (or modulus) less than unity due to the presence of off-diagonal terms in matrix in (S14) and vanish in $D^{N}(k\neq 0,\pi)$ when $N\to\infty$. Given the simple structure of the diagonal form of the $\lim_{N\to\infty}D^{N}(k\neq 0,\pi)=diag(1,0,0,0)$, it is easy to verify:

[TABLE]

I.3 The Residual Energy as $N\to\infty$

In this section we show the behavior of the residual energy $\varepsilon_{res}(NT)$ in the limit of $N\to\infty$ using the matrices given in (S22) and (LABEL:eq_neqkpi); for a transverse Ising chain one readily finds

[TABLE]

Let us recall that $\varepsilon_{res}(NT)=\frac{1}{L}\sum_{k}\left(e_{k}(NT)-e_{k}^{g}(0)\right)$ with $e_{k}(NT)=\langle\psi_{k}(NT)\rvert H_{k}^{0}\lvert\psi_{k}(NT)\rangle$ and $e_{k}^{g}(0)=\langle\psi_{k}(0)\rvert H_{k}^{0}\lvert\psi_{k}(0)\rangle$; in the thermodynamic limit ($L\to\infty$)

[TABLE]

Let us consider the trivial case $p=0$ (no periodic drive), for which the residual energy $\varepsilon_{res}(NT)=0$ as the system always remains in the initial ground state (except for a trivial phase factor) and hence two terms in (S29) identically cancel each other for all values of $N$. On the other hand, in the perfectly periodic situation ($p=1$), all the $k$-modes contribute and the steady state value becomes,

[TABLE]

It should be noted that the steady state value is attained in the asymptotic limit following a partial cancellation of the initial ground state energy .

Let us now immediately contrast this scenario with the case when $p\neq 1,0$ when the disordered average residual energy

[TABLE]

and observe that unlike in the steady state scenario in (S30), as the system gets heated up, only the $k=\pi$ mode contributes in (S31) to the cancellation to finally yield a bias $p$ as well as protocol independent value. Since, only one $k$-mode, the $k=\pi$ mode contributes, the asymptotic value of residual energy for $p\neq 1,0$ is evidently much greater than the steady state value, clearly asserting that the system has heated up to a finite asymptotic value. Equation (S31) yields the maximum energy that the system (after being heated up) can attain in the asymptotic limit which is independent of $p$ as well as the protocol.

However, one remaining question is that how does the system reach this asymptotic value for $p\neq 0,1$ as $N$ increases; this entirely depends on how the eigenvalues of the $D$-matrix (other than the one that sticks to unity) decay in diagonal form of $D^{N}(k)$ as $N$ increases, and hence, as we illustrate below, on the protocol and the value of $p$. This is precisely the reason we observe different initial growth of $\langle\varepsilon_{res}(NT)\rangle$ in Fig.  1a and 1b of the main text for different protocols and for different values of $p$ for a given protocol.

I.4 The Residual Energy at Large Intermediate $N$

In this subsection, we shall show how the growth of the RE depends on the non-universal features in the disorder operator. The D-matrix in Eq. (S14), as has already been mentioned can be diagonalised to obtain four eigenvalues, out of which two eigenvalues $1$ and $r=r(k,\alpha,\omega)$ are real and positive, while the rest two $\lambda_{c}=|\lambda_{c}(k,\alpha,\omega)|e^{\pm i\phi(k,\alpha,\omega)}$ are complex conjugates of each other with modulus less than unity. The three roots other than one are obtained by solving $f(k,\lambda)=0$, where $f(k,\lambda)$ is a third degree polynomial (see Eqs. (S14) and (S24)) with all real coefficients and is of the form $f(k,\lambda)=\lambda^{3}+a_{2}\lambda^{2}+a_{1}\lambda+a_{0}$:

where

[TABLE]

The real root $r=r(k,\alpha,\omega)$ is,

[TABLE]

and the two self-conjugate complex roots are,

[TABLE]

where, $S^{\pm}=(R\pm D)^{\frac{1}{3}}$, $D=Q^{3}+R^{3}$, $R=\frac{9a_{2}a_{1}-27a_{0}-2a_{2}^{3}}{54}$, $Q=\frac{3a_{1}-a_{2}^{2}}{9},$ and $\left(r+\lambda_{c}^{+}+\lambda_{c}^{-}\right)=a_{2}$ shows that $r$ is real.

The eigenvectors of $D(k)$ as a function of the roots $\lambda$ are given as,

[TABLE]

where,

[TABLE]

and

[TABLE]

Imposing the orthonormality condition over the eigenvectors of $D(k)$ one can write down the diagonalising matrix $S(k,\alpha,\omega)$ as:

[TABLE]

At large N, the contribution of $[\lambda_{c}^{\pm}]^{N}=|\lambda_{c}(k,\alpha,\omega)|^{N}e^{\pm iN\phi(k,\alpha,\omega)}$ to the residual energy summed over all the $k$-modes vanish, as the interferences due to the fast oscillating phases $e^{\pm iN\phi(k,\alpha,\omega)}$ cancel each other. This is apparent from Fig. S1, where for both the values of $k$, we see how the contributions of $[\lambda_{c}^{\pm}]^{N}$ to the exact RE$(k)$ (in red) oscillate around an increasing mean (in blue) set by contributions only from $1$ and $r^{N}$. For low $k$ (in Fig. SS1) we observe that the rapid oscillations due to $[\lambda_{c}^{\pm}]^{N}$ in the RE$(k)$ sit on top of the mean, whereas, for high values of $k$ (see Fig. SS1) the contributions of $[\lambda_{c}^{\pm}]^{N}$ to the RE$(k)$ are nearly zero and the mean approximately coincides with the exact RE$(k)$. Therefore, we can easily set the contributions of $[\lambda_{c}^{\pm}]^{N}$ as zero for each mode k throughout the rest of our calculations. However, $1$ and $r^{N}$ survives in such a large $N$ limit to yield,

[TABLE]

where $D_{d}=\operatorname{diag}\left(1,0,0,[r(k,\alpha,\omega)]^{N}\right)$.

Here the first matrix in the R.H.S is a constant matrix independent of the driving frequency, amplitude or protocol chosen and hence, universal. On the other hand, the second term in R.H.S. contains a non-universal matrix $D_{NU}(k,\alpha,\omega)$ which is of the form:

[TABLE]

where $u=u(k,\alpha,\omega)=\left[S^{-1}\right]_{41}$ and $v=v(k,\alpha,\omega)=\left[S^{-1}\right]_{42}$.

Although $D_{NU}(k,\alpha,\omega)$ specifically depends on the driving amplitude, frequency and the protocol implemented, it is essentially independent of the number of stroboscopic periods $N$ in this limit. All the time-dependence of the problem lies in the coefficient $\left[r(k,\alpha,\omega)\right]^{N}$ of the matrix $D_{NU}(k,\alpha,\omega)$. As this coefficient is positive and less than one, with increasing $N$, it gradually goes to zero only at $N\to\infty$ when the system finally attains its universal asymptotic value. But for any large non-zero $N<\infty$, this coefficient competes with the universal part (the first term in the R.H.S.) and generates a turning towards the asymptotic value (see Fig. 2 in the main text). This competition between the universal and the non-universal part is present irrespective of the protocol chosen, and hence, the curves for the residual energy in Fig. 2, are alike. Even though the slopes of the residual energy curves, depend upon both $r^{N}(k,\alpha,\omega)$ and $D^{N}(k,\alpha,\omega)$, the rate at which the curves for any protocol approach the asymptotic value with increasing $N$ is of course governed by the value of $r(k,\alpha,\omega)$ which is non-universal and varies depending upon the applied protocol.

II linear rise in the residual energy with low $N$

In this subsection, we shall use approximate analytical methods to show that for small $p$ (and small $1-p$, where the case with $p=1$ corresponds to the perfectly periodic situation) and low $N$, the residual energy (RE) (and also the difference in residual energy from the periodic steady state value) increases linearly with the number of stroboscopic periods $N$. We shall also investigate the variation of the RE with $p$ for a given number of stroboscopic periods ($N$). For the sake of convenience in explaining the results, we shall restrict our attention only to the aperiodically $\delta$-kicked situation. Our argument below is based on the notion of more probable configurations, for example if $p\to 0$, we shall consider the most probable configuration with either no kick or the second most probable configurations with only one kick being present in the entire process of driving. Similarly, for $(1-p)\to 0$, there are only two more probable configurations that one can probe, namely, the configurations with only one kick missing and the perfectly periodic configuration with $p=1$.

Let us first consider the case with $p\to 0$ and further assume that there is no kick up to the stroboscopic time $mT$ (where $1\leqslant m\leqslant N-1$). A kick is present at time $(m+1)T$; thereafter, the system evolves freely up to the time $NT$. The evolved state at the final time $NT$ is then given by

[TABLE]

where $\boldmath{\cal{F}}_{k}(T)$ is the usual Floquet operator and $U^{0}_{k}(T)=\exp(-iH^{0}_{k}T)$ is the time evolution operator for the free Hamiltonian $H^{0}_{k}$. For the $\delta$-kick situation, as discussed in the main text, we have an exact form of the Floquet operator, $\mathbb{\cal{F}}_{k}(T)=\exp(-i\alpha\sigma_{z})\exp(-iH^{0}_{k}T)$. Using the Baker-Campbell-Hausdorff formula, one can easily find the expectation value

[TABLE]

Let us recall that the ground state of the Hamiltonian $H^{0}_{k}=(1-\cos k)\sigma_{z}+(\sin k)\sigma_{x}$ can be written as $\lvert\psi_{k}(0)\rangle=(-f_{-},f_{+})^{T}$ , where $f_{\pm}=\sqrt{1/2(1\pm\sqrt{(1-\cos k)/2})}$ with the corresponding ground state energy $e^{g}_{k}(0)=-2\sin(k/2)$. Using the identity $\exp(-i\alpha\sigma_{z})=e^{-i\alpha}\lvert+\rangle\langle+\rvert+e^{i\alpha}\lvert-\rangle\langle-\rvert$ , where $\lvert+\rangle=(1,0)^{T}$ and $\lvert-\rangle=(0,1)^{T}$, and the relations $(f_{-}^{2}-f_{+}^{2})=\sin(k/2)$ and $f_{+}f_{-}=(1/2)\cos(k/2)$, we readily arrive at the expression

[TABLE]

Finally, we find the expression for the RE,

[TABLE]

What is important is that in the situation when only one kick is present in the entire driving (i.e., $p\to 0$), the RE is independent of the stroboscopic instant $m$ at which the kick is applied and further it is entirely determined by the strength of the kick. In this situation, the configuration averaged RE can be written as,

[TABLE]

where the first term corresponds to the no-kick situation whose contribution to the RE is zero, whereas the second term provides a non-zero contribution due to the presence of a single kick. Let us note that we have neglected the configurations with a higher number of kicks which occur with a vanishingly small probability in the limit $p\to 0$ and $N\to\infty$. Finally, in the limit $p\rightarrow 0$, we get the RE as

[TABLE]

Remarkably, the RE grows linearly with $N$ as shown in Fig. S2.

In the other limit, $(1-p)\rightarrow 0$, considering the two more probable terms in the configuration averaged RE of the system, one can similarly write,

[TABLE]

Notably, $\epsilon_{res}^{(N)}$ corresponds to the RE of the system for the perfectly periodic driving (as derived in the main text) and $\epsilon_{res}^{(N-1)}$ is the RE when only one kick is missing in the entire driving. The exercise to arrive at an analytic expression in this case is tedious and unilluminating, more because unlike the previous situation the quantity $\epsilon_{res}^{(N-1)}$ involves the sum over the stroboscopic instants at which the kick is missed. To circumvent this problem, we shall average over the all possible permutations i.e. sum over $m$ (position of the missed kick, $1<m<(N-1)$) with an equal weight $1/N$ (see Fig. S2). Henceforth, we see that the difference in the RE for the one missed kick situation and the perfectly periodic situation (i.e., the contribution to the RE that destabilizes the periodic steady state) indeed grows linearly with $N$ as shown in Fig. S2 upto an appropriate value of $N$ for which the approximation of single drive and one drive cycle missing holds.

Finally, we address the question how does variation of the RE (after a given number of stroboscopic periods) depend on the probability $p$; this however depends on the driving frequency $\omega$ and amplitude $\alpha$. In Fig. S3 and Fig. S3, we show that the variation is neither monotonic nor symmetric with $p$.