Dynamical enhancement of symmetries in many-body systems

TL;DR

This paper introduces a novel dynamical decoupling protocol using a self-similar drive to generate and enhance symmetries in many-body quantum systems, maintaining locality and avoiding heating over long times.

Contribution

It presents a new polyfractal drive method for creating multiple symmetries in many-body systems with stable, local Floquet Hamiltonians.

Findings

Effective Floquet Hamiltonian remains local

Symmetries can be generated dynamically in complex systems

Heating is suppressed exponentially long

Abstract

We construct a dynamical decoupling protocol for accurately generating local and global symmetries in general many-body systems. Multiple commuting and non-commuting symmetries can be created by means of a self-similar-in-time ("polyfractal") drive. The result is an effective Floquet Hamiltonian that remains local and avoids heating over exponentially long times. This approach can be used to realize a wide variety of quantum models, and non-equilibrium quantum phases.

Equations

$$t_i=m{2}^{i+n_s(j-1)}{T}_0,i\in [1,n_s],j\in [1,n_f],m\in {\mathcal{Z}}^{+}$$

$$H=\sum _{\mathbf{ϵ}}{A}_{\mathbf{ϵ}}\text{where}\mathbf{ϵ}=({ϵ}_1,...,{ϵ}_{n_s}),{ϵ}_i\in \{0,1\},$$

$$X_j{A}_{\mathbf{ϵ}}{X}_j=(-1{)}^{{ϵ}_j}{A}_{\mathbf{ϵ}}.$$

$$U(T\equiv 2^{n_f{n}_s}{T}_0)=\mathcal{T}\left\{e^{-i{\int }_0^{T}dt{\sum }_{\mathbf{ϵ}}{A}_{\mathbf{ϵ}}{f}_{\mathbf{ϵ}}(t)}\right\},$$

$$\text{where}{f}_{\mathbf{ϵ}}(t)=\pm 1\text{and}{\int }_0^{2^{n_s{n}_f}{T}_0}{f}_{\mathbf{ϵ}}(t)={\delta }_{0,\mathbf{ϵ}}.$$

$$U(T)=e^{-iT{\sum }_{n=0}^{\infty }{T}^n{\Omega }_n}$$

$$H_F^{(n_0)}=\sum _{0\le n\le n_0}{T}^n{\Omega }_n.$$

$${\tau }_X$$

$${\tau }_H$$

$${\tau }_X$$

$${\tau }_H$$

$$H_F^{(n_0)}\approx H_F^{(0)}=A_0$$

$$U(2^n{T}_0)=X\cdot U(2^{n-1}{T}_0)\cdot X\cdot U(2^{n-1}{T}_0),\forall n\ge 1$$

$$U(2T_0)=e^{-i{T}_0(A_0-A_1)}{e}^{-i{T}_0(A_0+A_1)}\equiv e^{-i{T}^{(1)}(A_0^{(1)}+A_1^{(1)})},$$

$$A_0^{(1)}$$

$$A_1^{(1)}$$

$$T^{(n_f)}=$$

$$A_1^{(n_f)}=$$

$$\Vert {\Omega }_n\Vert T^n\le N\frac{(cT\Vert h\Vert {)}^{n}n!}{(n+1{)}^2}\le N(cT\Vert h\Vert n{)}^n$$

$$\Vert \left[H_F^{(n_0)},X\right]\Vert \le 2T^{n_f}\Vert {\Omega }_{n_f}\Vert .$$

$$A_{00}^{(1)}=A_{00},A_{01}^{(1)}=-i{T}_0[A_{00},A_{01}],$$

$$A_{10}^{(1)}=-i\frac{T_0}{2}\left([A_{00},A_{10}]+[A_{01},A_{11}]\right),$$

$$A_{11}^{(1)}=-i{T}_0\left(-i\frac{T_0}{2}\right)\left[A_{00},[A_{00},A_{11}]+[A_{01},A_{10}]\right]$$

$$+i{T}_0\left(-i\frac{T_0}{2}\right)\left[A_{01},[A_{00},A_{10}]+[A_{01},A_{11}]\right].$$

$$H=\sum _{n,k\le 4}{e}^{-k+1}\left(-i{\gamma }_n{\gamma }_{n+k}\right)+V{\gamma }_n{\gamma }_{n+1}{\gamma }_{n+2}{\gamma }_{n+3}$$

$$P_X={\gamma }_2{\gamma }_3{\gamma }_6{\gamma }_7...=\prod _i{\gamma }_{4i-2}{\gamma }_{4i-1}$$

Full text

Introduction.—

Much of the richness of the material universe transpires due to a sequence of spontaneous symmetry breaking events, going from the highest (most symmetric) energy scales down to the lowest. By definition, every equilibrium physical system has already found its energetic low symmetry optimum; can it nevertheless be re-purposed to realize a different spontaneous symmetry breaking pathway?

Sometimes this can be achieved by tuning thermodynamic parameters. For instance, hydrostatic pressure can restore more symmetric crystalline phases at a given temperature, or magnetic fields can suppress superconductivity to reveal other competing instabilities. Such thermodynamic knobs are, unfortunately, quite limited.

A more flexible approach is to re-institute symmetries dynamically. In nuclear magnetic resonance, echo techniques have long been used to improve the coherence of local moments by dynamically suppressing their coupling to the environment. The Hahn echoHahn (1950) for instance reduces inhomogeneous broadening, and the WAHUHA protocolWaugh et al. (1968) can suppress anistropic dipole-dipole interactions. However, the extension of these ideas to creating global symmetries in many-body, interacting systems is less clear. The main concern is that driving concomitantly generates heating in these systems, which in turn suppresses interesting collective phenomena.

In this Letter, we discuss dynamical protocols that can be used to engineer multiple global or local symmetries, while keeping heating at bay for exponentially long times. This paves the way to creating novel symmetry broken and topological phases Senthil (2015a); Wen (2017) out of low symmetry templates.

Our work is informed by recent progress in understanding Floquet dynamics of many-body quantum systems. Crucially, it has been shown that the naive expectation that driving should inevitably lead to heating is not always correct—strong disorder Lazarides et al. (2015); Khemani et al. (2016); Ponte et al. (2015) and/or appropriate drive frequency selection Abanin et al. (2017a, 2016); Mori et al. (2016); Kuwahara et al. (2016); Else et al. (2017) can push heating to exponentially long times. Further, in the context of time-crystals, it is already well appreciated that driving may lead to the creation of $\mathcal{Z}_{2}$ symmetry whose spontaneous rupturing gives rise to time-crystalline phase von Keyserlingk and Sondhi (2016); Else et al. (2017, 2016). In this work, we show how such ideas can be extended to generate multiple local and global symmetries in the effective Floquet Hamiltonian.

We demonstrate our protocol by considering the specific case of how multiple $\mathcal{Z}_{2}$ symmetries can be generated in spin systems. (Generalization to the $\mathcal{Z}_{n}$ case with $n>2$ is straightforward.) The protocol involves injecting a finite set, of say, $n_{s}$ unitary operators $X_{i}$, at specific times corresponding to a fractal pattern, in between regular unitary evolution under the system’s physical Hamiltonian $H$. As we show, such fractal application of $X_{i}$ can be optimized in the number of fractal layers $n_{f}$, to result in an effective Hamiltonian for which $X_{i}$s are symmetries to an accuracy that is nearly exponential in the drive frequency. This sensitivity to the drive frequency allows for accurately implementing global symmetries while requiring a drive frequency that scales merely logarithmically with system size $N$.

The approach is motivated by the following intuition. The periodic application of $X_{i}$ flips the sign of terms in $H$ that anti-commute with $X_{i}$; thus, frequent application of $X_{i}$ leads to their cancellation after even number of drive periods. The resulting effective Hamiltonian commutes with $X_{i}$ up to $\sim\mathcal{O}(T_{0})$ corrections, where $T_{0}$ is the drive period. As we show, these corrections can be further suppressed by applying $X_{i}$ periodically at intervals of $2T_{0},4T_{0},...$ (see Fig. 1). (Note since $X^{2}_{i}=1$, $X_{i}$ will effectively not be applied at some times.) Alternatively, one may apply other generators $X_{j}$s at self-similar intervals to generate additional symmetries. This work follows through with the above intuition, and shows that for local Hamiltonians, there exists an optimal number of fractal layers $n_{f}$ which leads to superpolynomial (in drive-frequency) suppression of symmetry-violating terms, while heating occurs on a stretched-exponentially long timescale.

This work has implications for engineering novel Floquet Hamiltonians and non-equilibrium phases Moessner and Sondhi (2017); Kitagawa et al. (2010); Lindner et al. (2017); Rudner et al. (2013); von Keyserlingk and Sondhi (2016); Bairey et al. (2017); Choi et al. (2018); Yao et al. (2017); Agarwal et al. (2017), symmetry-protected topological phases Senthil (2015b), and stabilization of quantum memories Khodjasteh and Lidar (2007, 2005), among others. When the engineered symmetries anti-commute, they give rise to a robust degeneracy structure in the entire spectrum. This potentially could be exploited to engineer topologically protected qubits manipulated by the same $X_{i}$ to high precision. We explore some of these ideas in a companion paper Martin and Agarwal (2020).These ideas should also easily translate to experiments—particularly in setups exploring non-equilibrium quantum phases, such as trapped ions Zhang et al. (2017), NV centers Choi et al. (2017), among others Bordia et al. (2017); Bernien et al. (2017)—by introducing multiple drives, akin to those already used, but operating at multiples of the base frequency.

Finally, note that there is a precedent for fractal pulsing in finite-sized systems Khodjasteh and Lidar (2007, 2005). However, as opposed to the more general scheme we identify, these are system specific, and rely on operator expansions that have null convergence in the thermodynamic limit Abanin et al. (2017a). One may thus view this work as a formal extension of dynamical decoupling techniques to many-body systems. Below we describe our results, before providing detailed proofs and numerical validation.

Description of protocol and main results.—

We assume a physical system described by a Hamiltonian $H$ comprising of a sum of local terms, with a local norm $\left\lVert h\right\rVert$. We also assume a set of $n_{s}$ unitaries $X_{i}$ that further satisfy the condition $X^{2}_{i}=1$. These unitaries may themselves either commute or anti-commute with one another. The protocol we study involves applying $X_{i}$ at times

[TABLE]

amidst the regular Heisenberg evolution, where $n_{f}$ is the number of “fractal layers” in the composite Floquet unitary; see Fig. 1 for illustration. To illustrate by example, suppose we have $n_{s}=2$ unitaries, and apply these $n_{f}=1$ times. Then, $U(T_{0})=e^{-iHT_{0}}$; $U(2T_{0})=X_{1}U(T_{0})X_{1}U(T_{0})$; $U(4T_{0})=X_{2}U(2T_{0})X_{2}U(2T_{0})$. Subsequent time-evolution at periods of $T\equiv 2^{n_{f}n_{s}}T_{0}=4T_{0}$ is given by the repeated application of the Floquet unitary $U(4T_{0})$. For $n_{f}=2$, the above recursion relations would be repeated for another fractal layer: $U(8T_{0})=X_{1}U(4T_{0})X_{1}U(4T_{0})$; $U(16T_{0})=X_{2}U(8T_{0})X_{2}U(8T_{0})$, and $U(16T_{0})$ would subsequently serve as the Floquet unitary.

We now decompose $H$ into terms which transform differently under $X_{i}$:

[TABLE]

This decomposition is unique if $X_{i}$s commute or anti-commute with one another, which we assume. With this terminology, one may represent the Floquet unitary in time-ordered notation as

[TABLE]

Here $f_{\epsilon}(t)$ tracks times at which $X_{i}$ is applied; this corresponds to a sign change of terms $A_{\bm{\epsilon}}$ for which $\epsilon_{i}=1$. The integral over a complete period is zero except for $f_{\bm{0}}$. Thus, in a time-averaged sense, the effective Hamiltonian is $A_{\bm{0}}$ comprising of only terms even under all $X_{i}$.

Now, one may represent unitary $U(T)$ as an expansion in the exponent

[TABLE]

with operators $\Omega_{n}$ that can be arrived at using the Magnus expansion, or, in this case, a repeated application of the BCH formula; the first term is simply the time-averaged Hamiltonian $\Omega_{0}=A_{\bm{0}}$.

In general, the operator $\Omega_{n}$ involves $n$ nested commutators of the local operators $A_{\bm{\epsilon}}$. Thus, if the local terms comprising $A_{\bm{\epsilon}}$ involve at most $k$ sites, $\Omega_{n}$ can be represented as a sum of terms comprising at most $nk$ sites. Finally, the series expansion is only useful if we can truncate it at some order and effectively approximate the unitary dynamics; we define the approximate Hamiltonian

[TABLE]

Our main results concern the properties of the Floquet unitary $U(T)$, and the associated effective Floquet Hamiltonian $H^{(n_{0})}_{F}$. The first part of our results are directly adapted from the results of Refs. Mori et al. (2016); Kuwahara et al. (2016); Abanin et al. (2017a, b), which state that the difference between the exact reduced density matrix of a region of size $N_{\rho}$ and that obtained by evolving it with $H^{(n_{0})}_{F}$ is bounded in norm by $cN_{\rho}2^{-n_{0}}$, for some finite constant $c$. Here $n_{0}\sim 1/(T\left\lVert h\right\rVert)\equiv\omega/\left\lVert h\right\rVert\gg 1$ scales linearly with the effective drive frequency $\omega=1/\left(2^{n_{f}n_{s}}T_{0}\right)$. The norm of this error sets the inverse of the time scale up to which $H^{(n_{0})}_{F}$ provides a good description of the dynamics of local operators (alternatively, the time scale for heating)—crucially, this timescale grows exponentially with the drive frequency. For global operators $X_{i}$ with finite norm, but also for $H^{(n_{0})}_{F}$ itself, $N_{\rho}$ is the system size $N$, but the error can still be made small by scaling $\omega$ merely logarithmically with system size. Hereon, we will assume such frequency scaling.

Having established the conditions under which $H^{(n_{0})}_{F}$ faithfully describes the time evolution of $X_{i}$, we seek to establish a bound on the norm of terms in $H^{(n_{0})}_{F}$ that do not commute with $X_{i}$. Defining the time scale $\tau_{X}\approx\text{min}_{i}\left\{1/\left\lVert\left[H^{(n_{0})}_{F},X_{i}\right]\right\rVert\right\}$ which sets the shortest timescale at which unitaries $X_{i}$ relax, we find

[TABLE]

where $c_{1},c_{2}$ are $\mathcal{O}\left(1\right)$ combinatorial constants. $\tau_{H}$ is a bound on the time for which global operators such as energy are accurately described by $H^{(n_{0})}_{F}$ Mori et al. (2016); Abanin et al. (2017a). Note that $\tau_{X}$ initially increases with the number of fractal layers $n_{f}$ but eventually begins to decrease again. There is therefore an optimal $n_{f}$ for which $X_{i}$ become effective symmetries. Note also that $n_{f}$ cannot be made arbitrary large since its increase rapidly decreases the thermalization time scale $\tau_{H}$.

We now describe how to optimize $n_{f}$ to maximize $\tau_{X},\tau_{H}$. First, note that to maintain exponential dependence on the reference drive frequency $\omega_{0}\equiv 1/T_{0}$, $n_{f}$ must scale at most logarithmically in the small parameter $T_{0}\left\lVert h\right\rVert$. This implies $n_{f}=\frac{x}{n_{s}}\text{log}_{2}\left(\frac{1}{T_{0}\left\lVert h\right\rVert}\right)$ with $0<x<1$. Plugging this into the result for $\tau_{X}$, we find

[TABLE]

where $c^{\prime}_{2}$ is an $\mathcal{O}(1)$ constant. Thus, we can vary $\omega_{0}$ to control $\tau_{H}$ with (stretched-) exponential sensitivity, and $\tau_{X}$ as a power-law that can be made arbitrarily large. Consequently, a very slow increase of $\omega_{0}$ with system size $N$ is sufficient to cancel the prefactor of $1/N$ in both $\tau_{X}$ and $\tau_{H}$.

Finally, we note that $H^{(n_{0})}_{F}$ is quasi-local in the sense that the amplitude of terms decays exponentially with the spatial range Mori et al. (2016); Abanin et al. (2017a). In general, this operator may be hard to evaluate exactly, but it can be approximated by

[TABLE]

where $A_{\bm{0}}$ commutes with all $X_{i}$ by construction. Since it captures the time-evolution of local operators and, importantly, also reflects the global symmetry properties of $H^{(n_{0})}_{F}$, it is a good approximation to the effective Floquet Hamiltonian for times $t\lesssim\text{min}\left(\tau_{X},\tau_{H}\right)$.

The design of the protocol, which is crucial to the bound obtained in the first part of Eq. (6), its proof, and the result of Eqs. (7) that times $\tau_{X}$ and $\tau_{H}$ are almost exponentially sensitive to the base drive frequency are the central results of this work.

Fractal driving with a single unitary $X$.—

We now derive the bound in the first part of Eqs. (6) for the case ($n_{s}=1$) of a single unitary $X_{1}\equiv X$. The derivation of the result will also help the reader intuit the logic behind fractal driving.

Using the terminology introduced above, the Hamiltonian is composed of two (kinds of) terms: $H=A_{0}+A_{1}$, where $A_{0(1)}$ is even (odd) under $X$. In this case, fractal driving can be described by the simple recurrence relations

[TABLE]

with $U(T_{0})=e^{-iHT_{0}}$. At the first stage, this implies

[TABLE]

where we define $T^{(1)}\equiv 2T_{0}$, and $A^{(1)}_{0}$ and $A^{(1)}_{1}$ are the new effective terms that are even and odd, respectively, under $X$. The BCH formula then yields

[TABLE]

After $n_{f}$ fractal layers, this implies

[TABLE]

Importantly, terms in $H^{(n_{0})}_{F}$ that anti-commute with $X$ appear first at $\mathcal{O}\left(T_{0}^{n_{f}}\right)$. These terms are a subset of all terms that appear at $\mathcal{O}\left(T^{n_{f}}\right)$ in the expansion of the Floquet Hamiltonian, Eq. (5). Their norm is therefore bounded by $T^{n_{f}}\left\lVert\Omega_{n_{f}}\right\rVert$. Further, the norm of all terms that may anti-commute with $X$ can be bounded by $\sum_{n=n_{f}}^{n_{0}}T^{n}\left\lVert\Omega_{n}\right\rVert$. We note from Ref. Mori et al. (2016), that

[TABLE]

for some $\mathcal{O}(1)$ constant $c$. Using the above, we can bound the ratio $\left\lVert\Omega_{n+1}\right\rVert/\left\lVert\Omega_{n}\right\rVert<1/2\;\forall n\leq n_{0}$, if we set $n_{0}=1/(2cT\left\lVert h\right\rVert)$. This finally implies

[TABLE]

Eqs. (14) and (13) give the result in Eq. (6) for $n_{s}=1$.

Polyfractal driving with multiple unitaries.—

We now generalize to the case $n_{s}>1$. As before, we can examine the flow of $A_{\bm{\epsilon}}$ after each fractal layer, that is, at times $2^{nn_{s}}T_{0}$ for integer $n$. (Recall, a fractal layer corresponds to the application of each $X_{i}$ once at progressively doubled periods.)

Let us examine the recursion relations for $n_{s}=2$ to illustrate by example. The Hamiltonian in this case is $H=A_{00}+A_{01}+A_{10}+A_{11}$, where terms $A_{01}$ and $A_{10}$ are odd under $X_{1}$ and $X_{2}$ only, respectively, while $A_{11}$ is odd under both. $U(4T_{0})$ is calculated using the BCH expansion twice. Representing it as $U(4T_{0})=e^{-iT^{(1)}\left(A^{(1)}_{00}+A^{(1)}_{01}+A^{(1)}_{10}+A^{(1)}_{11}\right)}$, where $T^{(1)}=4T_{0}$, we find to $\mathcal{O}\left(T^{2}_{0}\right)$

[TABLE]

Note that terms which are odd under just one symmetry are canceled to $\mathcal{O}\left(T_{0}\right)$, while terms odd under both $X_{1},X_{2}$ are canceled to higher order. Similar conclusions apply for the general case of $n_{s}\geq 2$.

After $n_{f}$ fractal layers, symmetry-violating terms appear at order $\mathcal{O}\left(T^{n_{f}}_{0}\right)$ or higher. (The terms that are odd under just one symmetry appear at the lowest order.) To estimate the norm of these terms, we can apply the same arguments for the case $n_{s}=1$, arriving at the results of Eqs. (6). This completes the proof.

Numerical Results.—

We now provide numerical simulations to illustrate the above results. We consider a short-ranged spin-$1/2$ chain of length $L$, with open boundary conditions. In the majorana representation, the Hamiltonian reads

[TABLE]

where for $n$ odd/even, $\gamma_{n}=\prod_{j<n}\sigma^{z}_{j}\sigma^{x/y}_{n}$. $H$ has parity symmetry $P_{Z}=\prod_{n}\sqrt{i}\gamma_{n}$; we work in the sector $P_{Z}=1$.

One can check that driving with $P_{X}=\prod_{j}\sigma^{x}_{j}=\prod_{n}i\gamma_{4n-2}\gamma_{4n-1}$ suppresses even nearest neighbor majorana bonds and yields a Kitaev chain with terminal majorana zero modes. However, the quality (energy splitting) of the majoranas is not equivalent to the accuracy of $P_{X}$ symmetry; for details see Supplemental Material (SM).

To illustrate the effectiveness of the protocol in creating multiple global symmetries, we drive the system with $P_{X}$ and $P_{Z2}=\prod_{j}\sigma^{z}_{2j}=\prod_{n}i\gamma_{4n-3}\gamma_{4n-2}$. These operators (anti-) commute for $L=(4n+2)\;4n$ for integer $n$. For $L=4n+2$, the operators satisfy the Pauli algebra which leads to a doubly-degenerate spectrum.

To quantify the accuracy of the generated symmetries, we evaluate “decoherences”— $a(t)=1-\text{Tr}\left[P_{a}(t)P_{a}\right]/2^{L-1}$ of relevant operators $P_{a}$. $a(t=0)=0$ and remains zero for perfectly conserved $P_{a}$, while it relaxes to $1$ for non-conserved operators. We study $P_{a}=\{P_{X},P_{Z2},\sigma^{x}_{L/2}\sigma^{x}_{L/2+1}\}$. The first two measure the conservation of $P_{X},P_{Z2}$, and should yield $x(t)=z2(t)=0$ in case they are perfect symmetry generators of the Floquet dynamics, while $x_{l}(t)$ measures the relaxation of a local operator that is not expected to be conserved. Finally, we compute $\epsilon(t)=\left<A_{\bm{0}}(t)-A_{\bm{0}}(0)\right>$ where the average is taken with respect to the ground state of $A_{\bm{0}}$ [the part of $H$ that commutes with $P_{X}$ and $P_{Z2}$, see Eq. (8)]. $\epsilon(t)$ thus characterizes heating in the system.

The numerical results of Fig. 2 (a) generically exhibit rapid initial relaxation on microscopic timescales, before transitioning to a long-lived prethermal state. This is seen via the initial rapid loss of coherence of $P_{Z2}$, $P_{X}$ and increase in the energy $\epsilon(t)$, before plateauing at a fixed value much smaller than $1$. The decoherence $x_{l}(t)$ on the other hand rapidly approaches $1$, as expected. Long-term values of the decoherences, $\overline{z2},\overline{x}$ are seen to improve as fractal layers are increased from $n_{f}=1$, degrading subsequently for larger $n_{f}$—see Fig. 2 (b)—illustrating the existence of an optimal number of fractal layers for symmetry creation. Fig. 2 (c) illustrates the sensitivity of long-time coherences in the optimal protocol (over $n_{f}$) to $T_{0}$, as expected. Finally, in Fig. 2 (d), we confirm our expectations that the eigenspetrum of the Floquet unitary is doubly degenerate for $L=6$ and not for $L=8$.

An interesting aspect of the numerical results is that unlike our expectations, we do not observe heating away from the “prethermal plateau” at times longer than $\tau_{X}$ (Eq. (7). However, the plateau values of $\overline{x},\overline{z2}$ appear to scale with $1/\tau_{X}$. In fact, we observe relaxation (to $1$) only when driving at frequencies smaller than the microscopic scale, or for local operators not designed to be conserved. Whether this is a limitation of the small system sizes, a feature of the particular model that we considered, or an indication that our protocol works generally better than the conservative estimate for heating that we made, deserves further study.

Summary and Outlook.—

We have introduced a novel strong-driving protocol for engineering Floquet Hamiltonians, by creating new local and global symmetries. It may be viewed as an extension of dynamical decoupling techniques to local many-body Hamiltonians. While we describe here creation of $\mathcal{Z}_{2}$ symmetry generators, the results are easily generalized to $Z_{n>2}$ by applying the individual symmetry generators in sets of $n$ instead of twice, as in Eq. (9).

The symmetries can be used to engender a variety of novel Hamiltonians and dynamical phenomena. Creation of topological phases, and quantum memory stabilization using such schemes are explored in Ref. Martin and Agarwal (2020). While we have focused on the quasi-stationary Hamiltonian $H_{F}$, the dynamics of the system inside the Floqut period can be subject to interesting dynamical phenomena and deserves further attention.

Experiments probing non-equilibrium phenomena in driven systems in a variety of setups including Nitrogen-vacancy centers Choi et al. (2017), ion traps Zhang et al. (2017), cold atoms Bordia et al. (2017); Bernien et al. (2017) among others would be the natural setup to explore these ideas.

Acknolwedgements.—

We thank Jonathan Baugh for pointing out the “concatenated driving” protocol proposed in Refs. Khodjasteh and Lidar (2005, 2007); Khodjasteh et al. (2010); Cai et al. (2012) for protecting information in single qubits, which bears resemblance to our approach, and Lorenza Viola for pointing out previous work where dynamical decoupling schemes have been used for Hamiltonian simulation Santos and Viola (2008); Bookatz et al. (2014). We also thank B. Bauer, W. A. Coish, A. Pal, T. Pereg-Barnea, D. Pikulin, and L. Viola for valuable discussions. KA acknowledges support from NSERC Grants RGPIN-2019-06465, and DGECR-2019-00011, and start-up funds from McGill University for support. Work at Argonne was supported by the Department of Energy, Office of Science, Materials Science and Engineering Division.

I Kitaev Chain Example

We study the model as introduced in the main text, with $V=e^{-3}$. We drive the system with the unitary $P_{X}=\prod_{i}\sigma^{i}_{x}$. In the majorana formalism, this corresponds to the operator

[TABLE]

One can easily confirm that the Hamiltonian $H_{0}=(H+P_{X}HP_{X})/2$ which is invariant under $P_{X}$ exactly has weak and strong alternating bonds between nearest-neighbor majorana fermions. It is important to note that an imperfect weakening of the odd set of bonds is sufficient to drive the majorana system into the Kitaev phase (provided interactions are weak enough). In other words, provided the $\mathcal{Z}_{2}$ fermion parity is perfectly conserved, one does not need to introduce additional symmetries to enter the Kitaev phase.

The simplest demonstration of the above fact comes from numerical results of Figs. 3 and 4. In the first figure, we examine the coherence of the majorana at one end of the system. There is initial relaxation that occurs on the microscopic time scales followed by a long plateau where coherence is maintained, before an eventual relaxation. As seen in Fig. 3, the relaxation of the coherence occurs on a time scale that depends exponentially on length $L$ (equal linear displacements of the relaxation time scale can be seen on the log scale as $L$ is increased in fixed steps). In the same plot, dotted lines correspond to the protocol carried out at $n_{f}=2$, while solid lines correspond to $n_{f}=1$. In Fig. 4, we plot the coherence $x(t)$ corresponding to the operator $P_{X}$. Clearly, this coherence is improved as $n_{f}$ is increased. It is important to recall that the protocol engineers a particular symmetry (in this case $P_{X}$) and there is an optimal number of fractal layers $n_{f}$ associated with it. Phases that rely on the preservation of said symmetry will be more robust the better this symmetry is implemented. In this particular instance, the creation of isolated majorana modes relies on the system being in the correct phase (with alternating weak odd and strong even bonds) and not on the engineering of additional global symmetries (besides the fermion parity which we assume to be conserved). Thus, the coherence of the majoranas only weakly depends on the fractal layers $n_{f}$.

II Comments on Numerical Methods

Time evolution has been performed using either exact diagonalization or recursive multiplication to get the unitary for time-evolution at exponentially long times.

In the first method, we diagonalize $U(T_{f})$. Subsequently, all operators are represented in the basis that diagonalizes $U(T_{f})$. In this basis, the phase picked up by each individual Floquet eigenstate can be determined at time easily, and we do so at times that are uniformly spread out on a log scale. Coherences are evaluated as mentioned in the main text, averaging over a complete basis of initial states, except in the case of heating, where the initial state is the ground state of the Hamiltonian at zeroth order of the Magnus expansion, $H_{\vec{0}}$. This method usually works well except in instances where there is massive degeneracy in the spectrum; this makes precise basis rotation hard to compute when diagonalizing $U(T_{f})$. In the example studied in the main text, this unfortunately limits us to studying time-evolution for $L=4,8,12$ as $L=6,10,...$ are degenerate, which further limits our ability to perform finite size scaling in a meaningul way with system sizes accessible.

In the second method, which we use in instances where the spectrum has degeneracies (applicable to the Kitaev chain considered in the previous section), we simply compute unitary matrices for time-evolution at longer times by the recursion relation $U(2t)=U(t)\cdot U(t)$, for $t>T_{f}$. This resolves the issues with exact diagonalization in the presence of degeneracies, but is slower and suffers from the issue of the time-evolution matrix at late times becoming less unitary. This limits the dynamics to shorter times.