# Strong and almost strong modes of Floquet spin chains in Krylov subspaces

Daniel J. Yates  Aditi Mitra  Center for Quantum Phenomena, Department of
Physics, New York University, 726 Broadway, New York, NY, 10003, USA

###### Abstract

Integrable Floquet spin chains are known to host strong zero and $\pi$ modes
which are boundary operators that respectively commute and anticommute with
the Floquet unitary generating stroboscopic time-evolution, in addition to
anticommuting with a discrete symmetry of the Floquet unitary. Thus the
existence of strong modes imply a characteristic pairing structure of the full
spectrum. Weak interactions modify the strong modes to almost strong modes
that almost commute or anticommute with the Floquet unitary. Manifestations of
strong and almost strong modes are presented in two different Krylov
subspaces. One is a Krylov subspace obtained from a Lanczos iteration that
maps the time-evolution generated by the Floquet Hamiltonian onto dynamics of
a single particle on a fictitious chain with nearest neighbor hopping. The
second is a Krylov subspace obtained from the Arnoldi iteration that maps the
time-evolution generated directly by the Floquet unitary onto dynamics of a
single particle on a fictitious chain with longer range hopping. While the
former Krylov subspace is sensitive to the branch of the logarithm of the
Floquet unitary, the latter obtained from the Arnoldi scheme is not. The
effective single particle models in the Krylov subspace are discussed, and the
topological properties of the Krylov chain that ensure stable $0$ and $\pi$
modes at the boundaries are highlighted. The role of interactions is
discussed. Expressions for the lifetime of the almost strong modes are derived
in terms of the parameters of the Krylov subspace, and are compared with exact
diagonalization.

##  I Introduction

The Kitaev chain [1], which after a Jordan-Wigner transformation maps to the
transverse field Ising model (TFIM) [2], despite its apparent simplicity, is a
toy model for understanding diverse phenomena such as quantum phase
transitions [2] and topological systems [3]. It also forms a building block
for realizing non-Abelian braiding [4, 5, 6, 7]. In recent years, periodic or
Floquet driving of this model has helped conceptualize new phenomena such as
boundary and bulk discrete time crystals [8, 9, 10].

A key feature of the TFIM model, and its anisotropic generalization, is that
it hosts robust edge modes known as strong zero modes (SZMs) [1, 11]. An
operator $\Psi_{0}$ is a SZM if it obeys certain properties. Firstly it should
commute with the Hamiltonian $H$ in the thermodynamic limit i.e,
$\left[\Psi_{0},H\right]\approx 0$, secondly it should anticommute with a
discrete (say ${\mathbb{Z}}_{2}$) symmetry of the Hamiltonian ${\mathcal{D}}$,
$\\{\Psi_{0},{\mathcal{D}}\\}=0$, and thirdly it should be a local operator
with the property $\Psi_{0}^{2}=O(1)$. The presence of a SZM immediately
implies that the entire spectrum of $H$ is doubly degenerate, where the
degenerate pairs are $\\{|n\rangle,\Psi_{0}|n\rangle\\}$, with each member of
the pair being eigenstates of ${\mathcal{D}}$, but in the opposite symmetry
sector. This eigenstate-phase makes the edge modes extremely stable to adding
symmetry preserving perturbations, such as exchange interactions between
spins. In particular, the edge modes acquire a finite lifetime in the presence
of interactions, but with the lifetime being non-perturbative in the strength
of the interactions [12, 13, 14, 15, 16, 17]. These long-lived quasi-stable
modes in the presence of perturbations are referred to as almost strong zero
modes (ASZMs) [13].

Under Floquet driving, besides the SZMs, new edge modes arise [18, 19, 20, 21,
22], which are called strong $\pi$ modes (SPMs) [22], with the possibility of
having phases where SZMs and SPMs co-exist [23, 24, 25, 26, 22]. In order to
define the strong modes in a Floquet system, one needs to adapt the previous
definition of the SZM to that of a Floquet unitary over one drive cycle, $U$.
The SZM $\Psi_{0}$ and the SPM $\Psi_{\pi}$ are now defined as follows. Both
these operators anticommute with the discrete symmetry of $U$,
$\\{\Psi_{0},{\mathcal{D}}\\}=0,\\{\Psi_{\pi},{\mathcal{D}}\\}=0$. But while
the SZM $\Psi_{0}$ commutes with the Floquet unitary in the thermodynamic
limit $\left[\Psi_{0},U\right]\approx 0$, the SPM $\Psi_{\pi}$ anticommutes
with it $\\{\Psi_{\pi},U\\}\approx 0$ in the same limit. Moreover, as for the
static case, both operators are local, with the property
$\Psi_{0,\pi}^{2}=O(1)$. Thus existence of $\Psi_{0}$ implies that the entire
spectrum of $U$ is doubly degenerate $\\{|n\rangle,\Psi_{0}|n\rangle\\}$, with
the two eigenstates of a pair also being eigenstates of opposite symmetry. The
existence of $\Psi_{\pi}$ also implies that the spectrum is paired, but with
each pair $|n\rangle$, $\Psi_{\pi}|n\rangle$ not only being eigenstates of
opposite symmetry, but also being separated by the quasi-energy of $\pi/T$,
with $T$ being the period of the drive. In particular, the eigenstate pair
$|n\rangle,\Psi_{\pi}|n\rangle$ pick up a relative phase of $-1$ after odd
stroboscopic time-periods. Thus the SPM is an example of a boundary time-
crystal.

The Majorana mode of the Kitaev chain is an example of a SZM discussed above,
but written in the Majorana basis. From the discussion above it follows that
with open boundary conditions, not just the ground state, but all the
eigenstates of the Kitaev chain are doubly degenerate in the topologically
non-trivial phase. In a similar manner, Floquet driving the Kitaev chain gives
rise to $\pi$ Majorana modes. These are examples of the SPMs discussed above,
with their existence implying a $\pi$-pairing for the entire spectrum of the
Floquet unitary.

The strong modes are a property of the entire Hilbert space. When interactions
are included in the Kitaev chain, the SMs reduce to ASMs, i.e, the exact
pairing structure of the entire Hilbert space reduces to an approximate
pairing. However even with interactions, there is an integrable limit,
corresponding to the XYZ chain after a Jordan-Wigner transformation, where SMs
survive [11]. For a generic non-integrable chain which does not have SMs, when
the interactions are weak and hence irrelevant in a renormalization group
sense, one can still have zero modes that ensure a degeneracy only in the
ground state sector. These zero modes are known as weak zero modes [27, 28,
29].

In the Floquet setting, SZMs and SPMs have been identified in free fermion
(integrable) Floquet chains [18, 20, 22]. In addition, the effect of weak
interactions was explored using exact diagonalization (ED) where it was found
that quasi-stable and long lived edge modes exist despite bulk heating [22].
Drawing inspiration from the static case [13], these quasi-stable modes are
called ASZMs and almost strong $\pi$ modes (ASPMs) [22]. The current paper
proposes a route to understanding these quasi-stable modes by mapping their
dynamics to a Krylov subspace. A similar study was performed for ASZMs in
static chains [16, 17]. In addition, a recent paper showed how the Krylov
method can be generalized to Floquet systems, where the method used was a
Lanczos iteration scheme based on the Floquet Hamiltonian [30]. Here we
further expand on this method, and also introduce the Arnoldi method, which is
based on the Floquet unitary rather than the Floquet Hamiltonian, as a
promising alternate route. Since the Krylov subspace dynamics, whether
generated by the Lanczos method or the Arnoldi method, is the dynamics of a
free particle on a fictitous one-dimensional (1$d$) lattice, the Krylov chain,
this mapping paves the way for developing analytic methods to extract the
lifetimes for general interacting and driven settings.

The paper is organized as follows. The integrable model is introduced in
Section II, and its phases outlined. In Section III the Heisenberg time-
evolution of the SZM and SPM are mapped to a Krylov subspace via a Lanczos
iteration scheme where the generator of the dynamics is the logarithm of the
Floquet unitary over one drive cycle, i.e., the Floquet Hamiltonian. In
Section IV, the dynamics generated by the Floquet unitary is mapped to a
Krylov subspace via the Arnoldi iteration [31]. In section V analytic
expressions of the Krylov subspace arising from the Arnoldi method are derived
in the free limit. In Section VI, the Arnoldi iteration is used to arrive at a
compact expression for the lifetime of the edge modes that holds for finite
size systems as well as non-zero interactions. In Section VII, the effect of
interactions is discussed, and the lifetime of the edge modes obtained from
Krylov space methods are compared with ED. We present our conclusions in
Section VIII.

Figure 1: The autocorrelation function $A_{\rm\infty}$ at logarithmically
separated stroboscopic times for a chain of length $L=8$ and for the free
binary drive. The parameters are $g=0.3$ and $T=2.0,3.5,5.5,8.25$
corresponding respectively to SZM phase, SZM-SPM phase, trivial phase, and a
SPM phase.

##  II Model

In this section we introduce the integrable model of the Floquet chain. We
consider an open chain of length $L$ where the stroboscopic time evolution is
generated by the following Floquet unitary [20, 23, 25, 32, 33, 22]

$U=e^{-i\frac{T}{2}J_{x}H_{xx}}e^{-i\frac{T}{2}gH_{z}},$ (1)

where

$\displaystyle H_{z}$ $\displaystyle=\sum_{i=1}^{L}\sigma^{z}_{i},$ (2)
$\displaystyle H_{xx}$
$\displaystyle=\sum_{i=1}^{L-1}\sigma^{x}_{i}\sigma^{x}_{i+1}.$ (3)

In what follows we set $J_{x}=1$. Eq. (1) is a binary drive where the first
part of the drive involves time-evolution purely by the local magnetic field
of strength $gT$, while the second part of the drive involves time-evolution
with respect to the nearest-neighbor Ising interaction of strength $T$. The
Floquet unitary has a discrete symmetry as it commutes with

$\displaystyle\mathcal{D}=\sigma_{1}^{z}\dots\sigma_{L}^{z}.$ (4)

When we add interactions in Section VII, $\mathcal{D}$ will continue to be a
symmetry of the problem.

When $T\ll 1$, a high frequency expansion gives the Floquet Hamiltonian
$H_{F}=H_{xx}/2+gH_{z}/2$. Moreover, in this high frequency limit, $T$ can be
interpreted as the period of the drive, with the stroboscopic time-evolution
being $U=e^{-iH_{F}T}$. The phases of the Kitaev chain are recovered in the
high frequency limit, where as $g$ is tuned from $g<1$ to $g>1$, one
encounters a quantum phase transition from a topologically non-trivial phase
to a topologically trivial phase. We are interested in general $T$ where the
phase diagram is much richer. Eventually we are interested in how these phases
manifest in the Krylov subspace.

It is convenient to map the problem to Majorana fermions as follows

$\displaystyle
a_{2\ell-1}=\prod_{j=1}^{\ell-1}\sigma^{z}_{j}\sigma^{x}_{\ell};\,\,a_{2\ell}=\prod_{j=1}^{\ell-1}\sigma^{z}_{j}\sigma^{y}_{\ell},$
(5)

where $\ell=1\ldots L$. Denoting the vector

$\displaystyle\vec{a}=\begin{pmatrix}a_{1}\\\ a_{2}\\\ a_{3}\\\ \vdots\\\
a_{2L}\end{pmatrix},$ (6)

the stroboscopic time-evolution of $\vec{a}$ is as follows

$\displaystyle U^{\dagger}\vec{a}U=K\vec{a},$ (7)

where on defining

$\displaystyle c_{1}$ $\displaystyle=\cos(Tg),$ (8a) $\displaystyle c_{2}$
$\displaystyle=\cos(T),$ (8b) $\displaystyle s_{1}$ $\displaystyle=\sin(Tg),$
(8c) $\displaystyle s_{2}$ $\displaystyle=\sin(T),$ (8d)

we find that $K$ takes the following form [22] for $L=3$

$\displaystyle K$ $\displaystyle=\begin{pmatrix}c_{1}&-s_{1}&0&0&0&0\\\
s_{1}c_{2}&c_{1}c_{2}&-c_{1}s_{2}&s_{1}s_{2}&0&0\\\
s_{1}s_{2}&c_{1}s_{2}&c_{1}c_{2}&-s_{1}c_{2}&0&0\\\
0&0&s_{1}c_{2}&c_{1}c_{2}&-c_{1}s_{2}&s_{1}s_{2}\\\
0&0&s_{1}s_{2}&c_{1}s_{2}&c_{1}c_{2}&-s_{1}c_{2}\\\
0&0&0&0&s_{1}&c_{1}\end{pmatrix}.$ (9)

$K$ is an orthogonal matrix which, depending on the parameters can admit
eigenvalues at $\pm 1$. While $K$ has been shown explicitly for $L=3$, the
form for general $L$ can be easily guessed as the bulk consisting of the
second and third rows are repeated. For the example of $L=3$ given above, we
see that the bulk structure is repeated twice.

A key quantity that we will analyze is the autocorrelation function of
$\sigma^{x}_{1}=a_{1}$, evaluated at stroboscopic times

$\displaystyle A_{\infty}(nT)$ $\displaystyle=\frac{1}{2^{L}}{\rm
Tr}\biggl{[}\sigma^{x}_{1}(nT)\sigma^{x}_{1}(0)\biggr{]},$ (10a)
$\displaystyle\sigma^{x}(nT)$
$\displaystyle=\left[U^{\dagger}\right]^{n}\sigma^{x}_{1}U^{n}.$ (10b)

When strong modes exist, this quantity has a non-zero overlap with them, and
therefore its stroboscopic time-evolution is a diagnostic for whether strong
modes exist or not. This point will be further clarified below when we discuss
some exactly solvable limits. Moreover, from Eq. (7) it is clear that the
above autocorrelation function is completely determined by the action of $K$
on an initial vector that is localized on the first site. Thus $\pm 1$
eigenvalues of $K$ imply SZMs and SPMs respectively.

We now discuss some exactly solvable limits where $\sigma^{x}_{1}$ is the
strong mode, or has an $O(1)$ overlap with it. Below $n$ is an integer.

  * •

$Tg=(2n+1)\pi$, $T$ arbitrary: We have $e^{-iTgH_{z}/2}\propto{\mathcal{D}}$,
thus
$U^{\dagger}\sigma_{1}^{x}U=\mathcal{D}\sigma_{1}^{x}\mathcal{D}=-\sigma_{1}^{x},$
and a SPM exists. There is however no SZM.

  * •

$Tg=(2n)\pi$, $T$ arbitrary: Now $e^{-iTgH_{z}/2}\propto 1$, thus
$U^{\dagger}\sigma_{1}^{x}U=\sigma_{1}^{x}$. We now have a SZM while there is
no SPM.

  * •

$T=(2n+1)\pi$, $Tg$ arbitrary: Now we have
$e^{-iTH_{xx}/2}\propto\sigma_{1}^{x}\sigma_{L}^{x}$. With
$e^{-iTgH_{z}/2}=\prod_{l=1}^{L}\left[\cos(Tg/2)-i\sigma_{l}^{z}\sin(Tg/2)\right]$.
It is straightforward to check that

$\displaystyle U^{\dagger}\sigma_{1}^{x}U$
$\displaystyle=\cos(Tg)\sigma_{1}^{x}-\sin(Tg)\sigma_{1}^{y},$ $\displaystyle
U^{\dagger}\sigma_{1}^{y}U$
$\displaystyle=-\sin(Tg)\sigma_{1}^{x}-\cos(Tg)\sigma_{1}^{y}.$ (11)

The above linear combinations of $\sigma^{x}_{1}$ and $\sigma^{y}_{1}$ yield a
SZM and a SPM.

  * •

$T=2n\pi$, $Tg$ arbitrary: We have $e^{-iTH_{xx}/2}\propto 1$. Then,

$\displaystyle U^{\dagger}\sigma_{1}^{x}U$
$\displaystyle=\cos(Tg)\sigma_{1}^{x}-\sin(Tg)\sigma_{1}^{y},$ $\displaystyle
U^{\dagger}\sigma_{1}^{y}U$
$\displaystyle=\sin(Tg)\sigma_{1}^{x}+\cos(Tg)\sigma_{1}^{y}.$ (12)

The above linear combinations do not yield any strong modes unless $Tg=m\pi$
where $m$ is any integer. Here two cases arise depending on whether $m$ is
even or odd. If $Tg=2m\pi$ then $\sigma_{1}^{x}$ is trivially a SZM because
$U\propto 1$. If $Tg=(2m+1)\pi$, then $\sigma_{1}^{x}$ is a SPM because
$U\propto{\mathcal{D}}$.

The above examples show that whenever a SPM exists, the Floquet unitary has a
characteristic non-local structure by either being of the form
$U\propto{\mathcal{D}}$ (SPM phase) or being of the form
$U\propto\sigma^{x}_{1}\sigma^{x}_{L}$ (SZM-SPM phase) [23, 24, 25, 26]. We
will show later that this non-local structure will have implications on the
Krylov subspace.

The above phases are stable to perturbing around it. Fig. 1 shows $A_{\infty}$
for some special cases for stroboscopic and logarithmically separated times.
The parameters chosen are $g=0.3$, and $T=2.0,3.5,5.5,8.25$. We will work with
these parameters for the rest of the paper when presenting the numerical
results. Among the chosen parameters, all of them except $T=5.5$ host strong
modes. Thus for all parameters except $T=5.5$, $A_{\infty}$ is long-lived with
its lifetime only set by the system size $L$. Among these different cases,
$T=2.0$ is a phase that has a SZM, $T=3.5$ is a phase that simultaneously
hosts a SZM and a SPM, while $T=8.25$ hosts a SPM. For a SPM, the stroboscopic
time-evolution flips between $O(1)$ positive and negative values, consistent
with period-doubled dynamics.

When both SZM and SPM exist, $\sigma^{x}_{1}$ has an overlap with both of
them. Thus we can write, $\Psi_{0}=a\sigma^{x}_{1}+\ldots$ and
$\Psi_{\pi}=b\sigma_{1}^{x}(-1)^{n}+\ldots$ where $a,b=O(1)$ and positive
$a,b>0$. After the initial transients, we can write the form of
$\sigma^{x}_{1}(t)$ as $\sigma^{x}_{1}(t)\approx
c_{0}\Psi_{0}+c_{\pi}\Psi_{\pi}+\ldots$, where $c_{0,\pi}\leq 1/2$ for
normalization. Thus the autocorrelation function becomes
$A_{\infty}(nT)\approx c_{0}a+c_{\pi}b(-1)^{n}$. This shows that the signal
will oscillate between $c_{0}a+c_{\pi}b$ and $c_{0}a-c_{\pi}b$. There are
several indications that $c_{0}>c_{\pi}$. Firstly, note from Fig. 1 that when
only a SZM is present, $c_{0}\approx 1$, while when only a SPM is present
$c_{\pi}\approx 0.75$. This is further corroborated at the exactly solvable
limit for the SZM-SPM phase discussed in Eq. (11). Here one finds that at odd
stroboscopic times $A_{\infty}(n=1)=\cos(gT)$, while at even stroboscopic
times $A_{\infty}(n=2)=1$. Since $c_{0}>c_{\pi}$, the signal oscillates
between two positive numbers, as shown in Fig. 1.

Let us define the Floquet Hamiltonian $H_{F}$ as

$H_{F}=\frac{i}{T}\ln(U).$ (13)

Note that for the binary drive in Eq. (1), since the problem is free, the
matrix $K$ in Eq. (7) and $U$ are simultaneously diagonalized. Thus $i\ln{K}$
directly gives us the Floquet Hamiltonian $H_{F}$ in the Majorana basis.

The structure of $\ln{K}$ in the vicinity of some exactly solvable points was
discussed in Ref. 30. In particular it was shown that when $T\ll 1$, and
$|gT-\pi|\ll 1$, i.e, parameters for which the system hosts a SPM, $i\ln{K}$
has the form of a topologically non-trivial Su-Schrieffer-Heeger (SSH) model
[34, 35], upto an overall diagonal matrix $\pi$. Thus the energy of the zero
mode of the SSH model is shifted by $\pi$, giving a SPM. Moreover, such an
interpretation survived even when one tuned away from this exactly solvable
limit. Ref. 30 also showed that an interpretation of the strong modes in terms
of boundary modes of generalized SSH models, also survive in the Krylov
subspace constructed from the Lanczos method. We will further expand on this
point below, and also discuss the appearance of SSH-type chains in the Krylov
subspace constructed from the Arnoldi method.

##  III Krylov chain from the Lanczos iteration

The time evolution of an operator in a generic integrable or non-integrable
system can be mapped to single particle dynamics on a semi-infinite chain
employing a recursive Lanczos scheme [36]. This method has made a reapparance
recently as a way to identify chaotic dynamics [37, 38, 39, 40]. In this
section we outline this method. The exponential complexity of solving the
dynamics enters into the calculation of the hopping parameters on this chain
which we denote by $b_{n}$.

Note the definition of the Floquet Hamiltonian $H_{F}$ in Eq. (13). The
stroboscopic time-evolution after $m$ periods can be written in terms of
$H_{F}$ as follows

$\left[U^{\dagger}\right]^{m}OU^{m}=e^{iH_{F}mT}Oe^{-iH_{F}mT}=\sum_{n=0}^{\infty}\frac{(imT)^{n}}{n!}\mathcal{L}^{n}O,$
(14)

where we define

$\mathcal{L}O=[H_{F},O].$ (15)

To employ the Lanczos algorithm, we recast the operator dynamics into vector
dynamics by defining $|O)=O$. Since we are concerned with infinite temperature
quantities, we have an unambiguous choice for an inner product on the level of
the operators,

$(A|B)=\frac{1}{2^{L}}\text{Tr}\left[A^{\dagger}B\right].$ (16)

The Lanczos algorithm iteratively finds the operator basis that tri-
diagonalizes $\mathcal{L}$. We begin with the seed “state”, $|O_{1})$, and let
$\mathcal{L}|O_{1})=b_{1}|O_{2})$, where
$b_{1}=\sqrt{|\mathcal{L}|O_{1})|^{2}}$. The recursive definition for the
basis operators $|O_{n\geq 2})$ is,

$\mathcal{L}|O_{n})=b_{n}|O_{n+1})+b_{n-1}|O_{n-1}),$ (17)

where we define $b_{n\geq
2}=\sqrt{|\mathcal{L}|O_{n})-b_{n-1}|O_{n-1})|^{2}}$. It is straightforward
algebra to check that the above procedure will iteratively find basis
operators that yield a $\mathcal{L}$ which is tri-diagonal, and of the
following form,

$\mathcal{L}=\begin{pmatrix}&b_{1}&&&\\\ b_{1}&&b_{2}&&\\\ &b_{2}&&\ddots&\\\
&&\ddots&&\end{pmatrix}.$ (18)

This basis spanned by $|O_{n})$ lies within the Krylov subspace of
$\mathcal{L}$ and $|O_{1})$. We refer to this tri-diagonal matrix as the
Krylov Hamiltonian $H_{K}$,

$H_{K}=\sum_{n}b_{n}(c_{n}^{\dagger}c_{n+1}+c_{n+1}^{\dagger}c_{n}),$ (19)

and the 1$d$ lattice it represents, the Krylov chain.

An important aspect of this technique, often overlooked when discussing chaos,
is that the values of $b_{n}$ are highly dependent on the choice of seed
operator. Further, outside of special cases, namely a Hamiltonian that is
free, the exact solution to the operation $\mathcal{L}|O_{n})$ will require
ED, or similar methods with equivalent costs. This method does not escape the
rapidly growing exponential wall of complexity. In cases where the calculation
of all $b_{n}$ are possible, the above algorithm will return a value of
$b_{\text{end}}=0$ for the very last $b_{n}$. If the iteration is continued
further, all the previous $b_{n}$s will be repeated until $b_{\rm end}=0$ is
again reached.

Figure 2: Upper panel: The $b_{n}$s for the binary drive with $g=0.3$ and
system size $L=20$. Lower panel: The corresponding spectra with $x$-axis label
“index” denoting the location of the eigenvalues arranged in ascending order.
For both upper and lower panels, from left to right $T=2.0,3.5,5.5,8.25$.
These parameters correspond to SZM, SZM-SPM, trivial, and SPM phases
respectively. Horizontal red lines in lower panels correspond to energies $0$
and $\pm\pi$ in units of $T^{-1}$.

For free systems, the operation $\mathcal{L}|O_{n})$ can be efficiently solved
when in the Majorana basis. If the starting operator is a single Majorana then
the dimension of the Krylov subspace of that operator will scale as $2L$, as
free system dynamics can only mix the individual Majoranas among themselves.
Outside of free problems, the size of the full set of $|O_{n})$ will be large.
For example, a system size of $L$ will have $\sim 2^{2L}$ possible basis
operators. For all intents and purposes we treat $\mathcal{L}$ as a semi-
infinite chain. If the number of solved $b_{n}$ is insufficient for the
quantity of interest, an approach that works well is to supplement the known
set with approximate $b_{n}$ that are calculated based off of trends
established among the known hoppings [16, 17].

Starting with the seed state $|O_{1})=|\sigma_{1}^{x})$, we can recast
$A_{\infty}$ into the following form

$A_{\infty}(nT)=(e^{i\mathcal{L}nT})_{1,1}.$ (20)

Now, following the above discussion, the dynamics of $A_{\infty}$ has been
transformed into that of a semi-infinite single-particle problem. The details
of the semi-infinite chain will be discussed in subsequent sections. Ref. 16
showed that the slow dynamics of $A_{\infty}$ for static spin chains was a
result of topological modes residing at the left boundary (origin) of the
Krylov chain. We discuss the analogous situation for the Floquet problem.

Since for the binary drive, the problem is free, the matrix $K$ in Eq. (7) and
$U$ are simultaneously diagonalized. Thus we construct the Krylov Hamiltonian
by taking the Floquet Hamiltonian to be $TH_{F}=i\ln(K)$ and performing the
iterative steps outlined above. The parameters $b_{n}$ of the Krylov
Hamiltonian, and the corresponding spectra are shown in Fig. 2. The
eigenfunctions with eigenvalues corresponding to $0$ and $\pm\pi$ are shown in
Fig. 3.

Fig. 2 clearly shows that the Krylov chain has a dimerized structure like that
of a SSH chain. The case of the SZM phase at $T=2.0$ and the trivial phase at
$T=5.5$ are the easiest to understand as these correspond to SSH chains of
opposite sign of the dimerization, one sign being topologically trivial
($T=5.5$, with $b_{\rm odd}>b_{\rm even}$) and the other being topologically
non-trivial ($T=2.0$ with $b_{\rm odd}<b_{\rm even}$). For the case of the SPM
at $T=8.25$ (discussed in detail in Ref. 30), one finds that the very first
hopping is pinned at a large value of $O(\pi)$, while the rest of the chain,
from $n\geq 4$ onwards has a topologically non-trivial dimerization, implying
an edge mode. The local strong hopping at the first site results in an edge
mode that is pinned at $\pm\pi$. Thus $\pi$ modes of the Floquet chain appear
as regular zero modes of the topologically non-trivial Krylov chain with an
overall shift in energy by $\pi$ due to boundary conditions imposed by the
hoppings on the first few sites. For the more complicated SZM-SPM phase at
$T=3.5$, the first two sites at the edges act as a two-level system that can
host a pair of edge modes. These modes act as the boundary of an SSH chain
with non-trivial sign of the dimerization. This point will be demonstrated
further when we discuss the Krylov Hamiltonian obtained from the Arnoldi
method. Also note that all the eigenvalues of $H_{K}$ (lower panels of Fig. 2)
lie within the Floquet Brillouin Zone (FBZ) defined as $\epsilon
T\in\left[-\pi,\pi\right]$. However, this outcome is not guaranteed, and as we
will discuss in detail later, a different basis choice can lead to an unfolded
spectrum for the Krylov chain.

Fig. 3 plots the mod-square of modes at $0$ and $\pm\pi$ energies
($|\psi_{0,\pi}(n)|^{2}$) for all the phases. Note that these modes appear in
pairs in that, when a zero edge mode exists, there are two of them (orange and
blue lines in the upper panels of the plot). While when a $\pi$ mode exists,
there are a pair of them located at $\pi$ and $-\pi$ (orange and blue lines in
the lower panels of the plot). The periodicity of the spectrum under Floquet
driving implies that the $\pi$ modes are doubly degenerate in the same way as
the $0$ modes are because $-\pi=\pi+2\pi n$, where $n$ is any integer.

Fig. 3 shows that the degenerate pairs of edge modes are not symmetrically
localized on the two ends of the Krylov chain as the chain is constructed for
the operator $\sigma^{x}_{1}$, and is therefore heavily weighted on the left
edge. If the seed operator would have been a symmetric linear combination of
operators on the first and the last site of the physical chain, then the edge
modes of the Krylov chain would also have been symmetrically located at the
two ends of the chain. However it is interesting to note that the $\pi$ edge
mode pairs are a lot more asymmetric than the $0$ edge mode pairs, where for
the latter one does see two peaks, one each on the left and right edges of the
Krylov chain. In contrast the $\pi$ edge mode pairs are strongly peaked at the
left edge. This is because of the highly non-local structure of the Floquet
unitary when a $\pi$ mode exists. Recall that in the previous section we
showed that when a SPM exists, the Floquet unitary is proportional to the non-
local operators ${\mathcal{D}}$ (SPM phase) or $\sigma^{x}_{1}\sigma^{x}_{L}$
(SZM-SPM phase), and therefore has weight on both ends of the physical chain.
Thus the very first Lanczos step, even when it starts out with an operator
localized at the left edge, does end up “seeing” the right end of the chain,
resulting in the manifestation of the right edge mode on the left edge of the
Krylov chain.

Figure 3: The mod-square of the eigenfunctions of the Krylov chain for the
same parameters as Fig. 2. Upper (lower) panels correspond to modes at zero
($\pm\pi$) energy. The orange and blue data sets in each panel reflect the
fact that these modes appear in pairs. In particular, there are two zero edge
modes and two $\pi$ edge modes, the latter occurring at $\pm\pi$ energy. The
plots show that $T=5.5$ hosts no edge modes, consistent with a trivial phase.
On the other hand $T=2.0$ hosts a pair of edge modes at zero and no edge modes
at $\pi$ energy. $T=8.25$ hosts a pair of edge modes at $\pm\pi$ but no edge
modes at $0$ energy. $T=3.5$ hosts edge modes at both $0$ and $\pm\pi$
energies.

###  III.1 Majorana basis vs spin basis

Since the Floquet Hamiltonian can always be modified by shifting the quasi-
energies by integer multiples of $2\pi$, such shifts lead to different Krylov
Hamiltonians, but ultimately the same physics. In this section we discuss this
aspect in detail using the binary drive as an example. This issue was also
discussed in Ref. 30, both for the binary drive and the interacting ternary
drive, and it was pointed out that a particular choice of branchcut for
$\ln(U)$ gives a Krylov Hamiltonian that is more easy to interpret. This
particular choice is one where all the quasi-energies $\epsilon$ lie within
the FBZ defined by $\epsilon T\in\left[-\pi,\pi\right]$.

We now discuss different ways to construct the Krylov Hamiltonian for the
binary drive. One is working in the Majorana basis, and the second is working
in the spin basis. Moreover, working in the spin basis, we will discuss two
different Krylov chains, one where the spectrum is unfolded, while the other
where the spectrum is folded into the FBZ. Note that in the Majorana basis,
the spectrum is already naturally folded in the FBZ. This is because we have
direct access to the orthogonal matrix $K$ whose eigenvalues are pure phases.
Moreover, the eigenvalues of the Krylov Hamiltonian in the Majorana basis
equal those of $i\ln{K}$. With a choice of branchcut of the logarithm along
$(-\infty,0)$, this gives a spectrum for the Krylov Hamiltonian that is
bounded in the FBZ.

Fig. 4 constructs the Krylov Hamiltonian for the binary drive in a few
different ways. The first method is to work in the Majorana basis and to use
the Lanczos method where the generator of stroboscopic dynamics is $i\ln{K}$,
with $K$ given in Eq. (9). For this case, the Lanczos method is a simple basis
rotation and the spectrum of the Krylov Hamiltonian $H_{K}$, will be exactly
the same as that of $i\ln{K}$. That is, we begin with a single-particle
propagator $i\ln{K}$, and we end with a single-particle propagator
$\mathcal{L}$. The orange data labeled as ”free” in Fig. 4 corresponds to this
case where all the computations are performed in the Majorana basis. The top
panels of Fig. 4 shows the hopping parameters obtained from this procedure,
and these are identical to the top panels of Fig. 2. The middle panels show
the corresponding spectra and are identical to the lower panels of Fig. 2.
Note that for a system size of $L$, the Hilbert space is exhausted at $2L$ in
the Majorana basis. We are working with a system size of $L=4$. Hence top
panels show orange data that correspond to a total of $2L-1=7$ hoppings, while
middle panels show $2L=8$ eigenvalues. The lower panels show $A_{\infty}$ as
defined in Eq. (20) for this Krylov Hamiltonian. They agree perfectly with the
ED calculation (the orange data is masked by the ED data in the lower panels).

The second way to perform the Lanczos iteration is to work in the many-body or
spin basis, even though the problem is free. In this case, we are concerned
with the branch of the many-body Floquet Hamiltonian $i\ln{U}$, and not of the
single-particle Hamiltonian $i\ln{K}$. In Fig. 4, the blue data set labeled as
“spin1” corresponds to Lanczos performed with $TH_{F}=i\ln{U}$ without any
special treatment for the single-particle status of the problem. The seed
state is $\sigma_{1}^{x}$. We see that the resulting spectrum of $H_{K}$ in
the middle panel extends beyond the FBZ as indicated by the horizontal red
lines which pass through energies $\pm\pi$. The $b_{n}$ are shown in the top
row and share little in common with the $b_{n}$ from the Majorana or free
basis (orange crosses). The main differences between the blue data and the
orange data is as follows: the spectrum is unfolded for the former and folded
for the latter. Secondly, the former shows 3 gaps (at $\epsilon
T=0,\pi,-\pi$), while the latter shows only two gaps ($\epsilon T=0,\pi$ with
$\pi$ continuously connected to $-\pi$). This difference is due to the fact
that the blue data does not have information regarding the periodicity of $U$,
and so the gaps at $\pi$ and $-\pi$ appear to be distinct. The orange data has
information about the periodicity of the spectrum as it is obtained from the
exact $K$. Despite these differences, both Krylov Hamiltonians generate
identical autocorrelation functions, as can be seen in the lower panels of
Fig. 4.

We now wrap the blue spectrum of $H_{K}$ such that it fits within the FBZ.
This is now given by the green data set labeled “spin2” in the middle panels.
This folding requires transforming the Krylov Hamiltonian as follows [30]
$H_{K}\rightarrow U_{K}\hat{\epsilon}_{\rm FBZ}U_{K}^{\dagger}$, where
$\hat{\epsilon}_{\rm FBZ}$ is a diagonal matrix where all the energies lie in
the FBZ, and $U_{K}$ is the unitary matrix that diagonalizes $H_{K}$ before
the folding. After this folding, the Krylov Hamiltonian is no longer tri-
diagonal, and one needs to perform a second Lanczos iteration to return to a
tri-diagonal Krylov subspace. The resulting $b_{n}$ of this new $H_{K}$ are
shown in the top panel in green, and are labeled as “spin2”. We note that this
wrapping of the spectrum requires us to fully diagonalize the problem, which
will typically not be possible for us when dealing with large system sizes. We
see that the green $b_{n}$ of the top row match the $b_{n}$ obtained from the
single-particle Majorana basis upto $b_{n=2L-1}$. At $n=2L$, $b_{2L}=0$ so
that all the $b_{n}$ appearing after that represent a second chain decoupled
from the part that lies between $n=1\ldots 2L-1$. The significance of these
decoupled chains on the spectrum is shown in the middle panels (green data).
Here we see that the spectrum is repeated twice for our example. This
repetition does not influence the physics, as can be seen in the lower panel
which plots $A_{\infty}$. This is because once $b_{m}=0$, the remaining
$b_{n>m}$ do not affect the dynamics. The bottom panel in Fig. 4 plots
$A_{\infty}$ for the Krylov Hamiltonian obtained after the folding. The
results agree perfectly with ED, and the two other Krylov Hamiltonians
discussed previously.

Thus while the $b_{n}$ and the spectra (top and center rows of Fig. 4) are
sensitive to whether the Lanczos is performed on the level of the many-body or
spin basis versus the single-particle or Majorana basis, the dynamics is not
sensitive to it. As touched upon in the previous paragraph, to align the
Majorana picture and the spin picture, we needed to wrap the spectrum of
$H_{K}$, requiring the full diagonalization of $H_{K}$. This is only possible
in the free case where we can easily find the full Krylov subspace of the seed
operator. Moreover, the reason for the very different $b_{n}$s in the spin and
in the Majorana bases, arises due to the ambiguity associated with the branch
of $H_{F}=i\ln(U)/T$. To avoid this ambiguity, it is better to set up the
problem in an alternate Krylov subspace that involves working directly with
the Floquet unitary $U$, rather than with the Floquet Hamiltonian $H_{F}$. We
discuss this alternate approach in the next section.

Figure 4: Exploring the difference between the spin basis and the Majorana
basis for the binary drive with $g=0.3$ and $L=4$. All rows, from left to
right $T=2.0,3.5,5.5,8.25$. These correspond respectively to SZM, SZM-SPM,
trivial, and SPM phases. Top rows show the $b_{n}$s for the Majorana basis
(labeled as “free”), and the two spin bases (labeled as spin1,2), where spin1
is the basis that gives the unfolded spectrum, while spin2 is the basis that
gives the spectrum folded into the FBZ. Middle rows show the spectra of the
corresponding Krylov Hamiltonians, with the $x$-axis label “index” labeling
the eigenvalues. The bottom rows plot $A_{\infty}$ obtained from performing
the time-evolution using the three different Krylov Hamiltonians. The red data
labeled as “ED” in the bottom rows is $A_{\infty}$ from ED.

##  IV Krylov chain from the Arnoldi iteration

We now discuss a different Krylov subspace, one that arises from the action of
$U$ on the seed operator and is known as the Arnoldi iteration [31]. This
differs from the previous section where the generator of the dynamics was a
Hermitian operator, the Floquet Hamiltonian $TH_{F}=i\ln{U}$. Instead, we now
work with a unitary operator, and below we outline the steps for obtaining the
corresponding Krylov subspace.

Instead of $\mathcal{L}$, we now have the quantity $W$ which is defined as
follows

$\displaystyle W|O)$ $\displaystyle=U^{\dagger}\hat{O}U$ (21) $\displaystyle
W^{n}|O)$
$\displaystyle=\left[U^{\dagger}\right]^{n}\hat{O}\left[U\right]^{n}.$ (22)

$W$ is unitary because

$\displaystyle W^{\dagger}W|O)$
$\displaystyle=U(U^{\dagger}\hat{O}U)U^{\dagger}$ (23)
$\displaystyle=\hat{O}=WW^{\dagger}|O).$ (24)

We now outline the Arnoldi method. This reduces to the Lanczos scheme outlined
in the previous section when $W$ is Hermitian. In particular we will see that
the unitarity of $W$ will no longer produce a simple tri-diagonal matrix.
Hermiticity is needed for the appearance of a tri-diagonal form as we saw in
the previous section on the Lanczos iteration scheme.

Let $|1)=\sigma_{1}^{x}$, then assuming we have found the orthonormal basis
vectors $|n),|n-1),\dots$, we find $|n+1)$ by time evolving it, and projecting
out the overlaps with the known basis vectors. Thus we have

$\displaystyle|n+1^{\prime})$ $\displaystyle=W|n)-\sum_{l=1}^{n}|l)(l|W|n)$
$\displaystyle=W|n)-\sum_{l=1}^{n}w_{l,n}|l)$ (25)
$\displaystyle=\left[1-\sum_{l=1}^{n}|l)(l|\right]W|n)=P_{n}W|n).$ (26)

Above $w_{l,n}=(l|W|n)$ and $P_{n}=1-\sum_{l=1}^{n}|l)(l|$ projects out
overlaps with the previously calculated basis vectors. Following this we
normalize $|n+1^{\prime})$,

$|n+1)=\frac{|n+1^{\prime})}{\sqrt{(n+1^{\prime}|n+1^{\prime})}}.$ (27)

We note that

$\displaystyle(n+1^{\prime}|n+1^{\prime})$
$\displaystyle=(n+1^{\prime}|P_{n}W|n)=(n+1^{\prime}|W|n),$ (28)

because $(n+1^{\prime}|l)=0$ for $l=1\ldots n$ due to the projector $P_{n}$.
Using Eq. (27) the above becomes

$\displaystyle(n+1^{\prime}|n+1^{\prime})=\sqrt{(n+1^{\prime}|n+1^{\prime})}(n+1|W|n),$
(29)

implying that

$\displaystyle\sqrt{(n+1^{\prime}|n+1^{\prime})}$
$\displaystyle=(n+1|W|n)=w_{n+1,n}.$ (30)

Thus we start with $|1)$, find $W|1)$ and $(1|W|1)=w_{1,1}$. Following this we
find the unnormalized $|2^{\prime})=W|1)-w_{1,1}|1)$. Then we determine the
coefficient $w_{2,1}$ from $(2^{\prime}|2^{\prime})=w_{2,1}^{2}$. This first
step fills in the first two rows of the first column of $W$ in the Arnoldi
basis.

Flipping Eq. (25) around, we have

$W|n)=w_{n+1,n}|n+1)+\sum_{l=1}^{n}w_{l,n}|l).$ (31)

Thus as we iterate through this Arnoldi algorithm, we will find that $W$ has
the following upper Hessenberg form

$W=\begin{pmatrix}w_{1,1}&w_{1,2}&\dots\\\ w_{2,1}&w_{2,2}&\dots\\\
0&w_{3,2}&\dots\\\ 0&0&\ddots\end{pmatrix}.$ (32)

An upper Hessenberg matrix is a square matrix whose elements $w_{i,j}=0$ for
$i>j+1$, i.e, it is a matrix whose elements are zero below the first sub-
diagonal [41]. In addition, when the seed operator is Hermitian, all the
elements of $W$ are real because under time-evolution a Hermitian operator
stays Hermitian, with the matrix elements of $W$ simply denoting the weight of
different Hermitian operators at a particular step in the iterative procedure.
Typically the structure of $W$ is such that the elements on the first sub-
diagonal are dominant. These elements besides measuring the norm of the
operators (see Eq. (30)), also measure how much the operator is spreading into
new parts of the Krylov subspace. In contrast, all other elements of $W$ are
simply the overlap of the new operator generated at each iteration, with
previous elements of the Krylov subspace. Moreover, the further an element is
from the first sub-diagonal, the smaller it is. This observation will be handy
later when we derive analytic expressions for $W$.

The autocorrelation function in terms of $W$ now has the form

$A_{\infty}(nT)=(1|W^{n}|1).$ (33)

Just as for $\mathcal{L}$, the iteration must be cut short after $N$ steps.
Due to the truncation, the resulting $W$ will not be unitary. As we will show
later in the context of ASZMs and ASPMs, as long as $N$ is sufficiently large,
the truncated $W$ does well in reproducing the dynamics. The success of the
approximation is a good indication that if the edge modes are sufficiently
localized at the edge of the Krylov chain, the truncation scheme does not
affect the physics.

We note that for the last state $|N)$ we have,

$W|N)=\sum_{l=1}^{N}w_{l,N}|l)+w_{N+1,N}|N+1),$ (34)

and by stopping at $|N)$, we are making some form of a Markov approximation at
a rate of $w_{N+1,N}$. Thus we expect almost all operator dynamics produced by
a truncated $W$ to eventually cause a decay in time.

###  IV.1 Strong modes in limiting cases

We now discuss the form of the $W$ matrix for some analytically tractable
limiting cases. These limits were already discussed in previous sections in
the context of the matrix $K$ in Eq. (9).

  * •

$Tg=(2n+1)\pi$, $T$ arbitrary: For this case
$e^{-iTgH_{z}/2}\propto\mathcal{D}$, thus
$U^{\dagger}\sigma_{1}^{x}U=\mathcal{D}\sigma_{1}^{x}\mathcal{D}=-\sigma_{1}^{x},$
and a SPM exists, whereas there is no SZM. The first step of the Arnoldi
iteration is $|1)=\sigma_{1}^{x}$ and the next state is

$|2^{\prime})=W|1)=-|1).$ (35)

Thus the Krylov subspace terminates after the first basis vector
$|1)=\sigma_{1}^{x}$, and $W$ is a ”matrix” with a single element, $-1$.

  * •

$Tg=(2n)\pi$, $T$ arbitrary: For this case $e^{-iTgH_{z}/2}\propto 1$ and
$U^{\dagger}\sigma_{1}^{x}U=\sigma_{1}^{x}$. We have a SZM, but no SPM.
Following the same steps as in the previous case, it is clear that $W=1$.

  * •

$T=(2n+1)\pi$, $Tg$ arbitrary: We now have
$e^{-iTH_{xx}/2}\propto\sigma_{1}^{x}\sigma_{L}^{x}$, with
$e^{-iTgH_{z}/2}=\prod_{l}\left[\cos(Tg/2)-i\sigma_{l}^{z}\sin(Tg/2)\right]$.
For this case Eq. (11) holds. We now construct the Krylov subspace for
$\sigma_{1}^{x}$. Beginning with $|1)=\sigma_{1}^{x}$, we have using Eq. (11)

$W|1)=\cos(Tg)\sigma_{1}^{x}-\sin(Tg)\sigma_{1}^{y},$ (36)

giving,

$|2^{\prime})=W|1)-\cos(Tg)|1)=-\sin(Tg)\sigma_{1}^{y}.$ (37)

After normalization we obtain

$|2)=-\sigma_{1}^{y}.$ (38)

Moreover, again using Eq. (11), we can see that $W|2)$ will be a linear
combination of $|1),|2)$, and in particular

$W=\begin{pmatrix}\cos(Tg)&|\sin(Tg)|\\\ |\sin(Tg)|&-\cos(Tg)\end{pmatrix},$
(39)

which corresponds to having both a SZM and a SPM because the eigenvalues of
$W$ are $\pm 1$.

  * •

$T=2n\pi$ and $Tg$ arbitrary: For this case we have
$e^{-iTJ_{x}H_{xx}/2}\propto 1$ and the time-evolution is given by Eq. (12).
Carrying out the same steps as before one finds that the Krylov subspace for
$\sigma_{1}^{x}$ is two dimensional with

$W=\begin{pmatrix}\cos(Tg)&-|\sin(Tg)|\\\ |\sin(Tg)|&\cos(Tg)\end{pmatrix}.$
(40)

Thus the eigenvalues are pure phases that do not equal $\pm 1$, unless
$Tg=n\pi$, where $n$ is an integer. When $n\in$ even, we have a SZM, and when
$n\in$ odd, we have a SPM.

##  V Arnoldi iteration for the binary drive

In this section we present the analytic form of $W$ for the binary drive, and
discuss it for some limiting cases which are more general than those discussed
in the previous section. This discussion is aimed at showing the connection
between strong modes and the topological properties of the Hamiltonian
$i\ln{W}$.

Recall that for the binary drive, we can work in the Majorana basis where a
chain of length $L$ requires diagonalizing the problem in a reduced Hilbert
space of size $2L$. It is straightforward, but tedious to verify that for a
chain of length $L=4$ the Arnoldi alogrithm produces the following $8\times 8$
$W$ matrix,

$W=\begin{pmatrix}c_{1}&-s_{1}c_{2}&s_{2}s_{1}c_{1}&-s_{2}s_{1}^{2}c_{2}&s_{2}^{2}s_{1}^{2}c_{1}&-s_{2}^{2}s_{1}^{3}c_{2}&s_{2}^{3}s_{1}^{3}c_{1}&-s_{2}^{3}s_{1}^{4}\\\
s_{1}&c_{2}c_{1}&-s_{2}c_{1}^{2}&s_{1}c_{1}s_{2}c_{2}&-s_{2}^{2}s_{1}c_{1}^{2}&s_{2}^{2}s_{1}^{2}c_{2}c_{1}&-s_{2}^{3}s_{1}^{2}c_{1}^{2}&s_{2}^{3}s_{1}^{3}c_{1}\\\
0&s_{2}&c_{2}c_{1}&-s_{1}c_{2}^{2}&s_{1}c_{1}s_{2}c_{2}&-s_{2}s_{1}^{2}c_{2}^{2}&s_{2}^{2}s_{1}^{2}c_{2}c_{1}&-s_{2}^{2}s_{1}^{3}c_{2}\\\
0&0&s_{1}&c_{2}c_{1}&-s_{2}c_{1}^{2}&s_{1}c_{1}s_{2}c_{2}&-s_{2}^{2}s_{1}c_{1}^{2}&s_{2}^{2}s_{1}^{2}c_{1}\\\
0&0&0&s_{2}&c_{2}c_{1}&-s_{1}c_{2}^{2}&s_{1}c_{1}s_{2}c_{2}&-s_{2}s_{1}^{2}c_{2}\\\
0&0&0&0&s_{1}&c_{2}c_{1}&-s_{2}c_{1}^{2}&s_{1}c_{1}s_{2}\\\
0&0&0&0&0&s_{2}&c_{2}c_{1}&-s_{1}c_{2}\\\
0&0&0&0&0&0&s_{1}&c_{1}\end{pmatrix},$ (41)

where $c_{1,2},s_{1,2}$ were defined in Eq. (8). We discuss the structure of
$W$ as it will help us make further approximations. Firstly since $W$ is
unitary, its rows and columns form an orthonormal basis. Secondly, for even
rows, starting on the column corresponding to the sub-diagonal, we have

$\\{s_{1},\ c_{1}c_{2},\ -s_{2}c_{1}^{2},\ (s_{1}s_{2})c_{1}c_{2},\
-(s_{1}s_{2})s_{2}c_{1}^{2},\ \dots,\ (s_{1}s_{2})^{n}(c_{1}c_{2}),\
(s_{1}s_{2})^{n}(-s_{2}c_{1}^{2}),\ \dots\\}.$ (42)

Thirdly, for odd rows, starting on the column corresponding to the sub-
diagonal, we have

$\\{s_{2},\ c_{1}c_{2},\ -s_{1}c_{2}^{2},\ (s_{1}s_{2})c_{1}c_{2},\
-(s_{1}s_{2})s_{1}c_{2}^{2},\ \dots,\ (s_{1}s_{2})^{n}(c_{1}c_{2}),\
(s_{1}s_{2})^{n}(-s_{1}c_{2}^{2}),\ \dots\\}.$ (43)

The rows and columns at the edge of the matrix can be determined from using
the above rules and then dividing by $c_{2}$. This also ensures orthonormality
of the rows and columns. In what follows we will further explore $W$, and in
particular $\ln{W}$ for limiting cases that give SZMs and SPMs. The reason for
exploring $\ln{W}$ is that it has a Hamiltonian description, and therefore the
topological properties can be more easily discerned.

###  V.1 $\ln{W}$ for $gT\ll 1$, $T\ll 1$: SZM phase

First let us consider the case when $gT\ll 1$. For this case we have
$s_{2},c_{2}$ arbitrary, and $s_{1}\ll 1,c_{1}\approx 1$. Truncating the $W$
matrix in Eq. (41) to terms of $O(s_{1})$ we have

$\displaystyle W=\begin{pmatrix}1&-s_{1}c_{2}&s_{2}s_{1}&0&0&0&0&0\\\
s_{1}&c_{2}&-s_{2}&s_{1}s_{2}c_{2}&-s_{2}^{2}s_{1}&0&0&0\\\
0&s_{2}&c_{2}&-s_{1}c_{2}^{2}&s_{1}s_{2}c_{2}&0&0&0\\\
0&0&s_{1}&c_{2}&-s_{2}&s_{1}s_{2}c_{2}&-s_{2}^{2}s_{1}&0\\\
0&0&0&s_{2}&c_{2}&-s_{1}c_{2}^{2}&s_{1}s_{2}c_{2}&0\\\
0&0&0&0&s_{1}&c_{2}&-s_{2}&s_{1}s_{2}\\\ 0&0&0&0&0&s_{2}&c_{2}&-s_{1}c_{2}\\\
0&0&0&0&0&0&s_{1}&1\end{pmatrix}+O(s_{1}^{2})=W_{0}+\delta W,$ (44)
$\displaystyle W_{0}=\begin{pmatrix}1&0&0&0&0&0&0&0\\\
0&c_{2}&-s_{2}&0&0&0&0&0\\\ 0&s_{2}&c_{2}&0&0&0&0&0\\\
0&0&0&c_{2}&-s_{2}&0&0&0\\\ 0&0&0&s_{2}&c_{2}&0&0&0\\\
0&0&0&0&0&c_{2}&-s_{2}&0\\\ 0&0&0&0&0&s_{2}&c_{2}&0\\\
0&0&0&0&0&0&0&1\end{pmatrix}=W(gT=0),$ (45) $\displaystyle\delta
W=\begin{pmatrix}0&-s_{1}c_{2}&s_{2}s_{1}&0&0&0&0&0\\\
s_{1}&0&0&s_{1}s_{2}c_{2}&-s_{2}^{2}s_{1}&0&0&0\\\
0&0&0&-s_{1}c_{2}^{2}&s_{1}s_{2}c_{2}&0&0&0\\\
0&0&s_{1}&0&0&s_{1}s_{2}c_{2}&-s_{2}^{2}s_{1}&0\\\
0&0&0&0&0&-s_{1}c_{2}^{2}&s_{1}s_{2}c_{2}&0\\\ 0&0&0&0&s_{1}&0&0&s_{1}s_{2}\\\
0&0&0&0&0&0&0&-s_{1}c_{2}\\\ 0&0&0&0&0&0&s_{1}&0\end{pmatrix}.$ (46)

Above $W_{0}=W(gT=0)$, while $\delta W$ is the leading correction about this
point and is $O(s_{1})$. If in addition we impose $T\ll 1$ which is equivalent
to $s_{2}\ll 1$, the matrices $W_{0}$ and $\delta W$ commute if terms of
$O(s_{1}s_{2}),O(s_{1}^{2}),O(s_{2}^{2})$ and higher are dropped. Recall that
if $A$ and $B$ are two commuting matrices $\left[A,B\right]=0$, then,

$\displaystyle\ln{(A+B)}\approx\ln\left[A\right]+\left[A\right]^{-1}B+\ldots.$
(47)

Using, Eq. (47) and setting $A=W_{0},B=\delta W$, we find the following
expressing to first order in $s_{1}$ and $T$

$\displaystyle
i\ln\left[W\right]\approx\left(\begin{array}[]{cccccccc}0&-is_{1}&0&0&0&0&0&0\\\
is_{1}&0&-iT&0&0&0&0&0\\\ 0&iT&0&-is_{1}&0&0&0&0\\\ 0&0&is_{1}&0&-iT&0&0&0\\\
0&0&0&iT&0&-is_{1}&0&0\\\ 0&0&0&0&is_{1}&0&-iT&0\\\ 0&0&0&0&0&iT&0&-is_{1}\\\
0&0&0&0&0&0&is_{1}&0\\\ \end{array}\right).$ (56)

Above $i\ln{W}$ represents a Floquet Hamiltonian that is effectively a SSH
model with topologically non-trivial sign of the dimerization for $|s_{1}|<T$.
Thus a zero mode is guaranteed. This shows that the SZM of the original
Hamiltonian manifests as an edge mode of a topologically non-trivial
Hamiltonian $i\ln{W}$ in the single particle Krylov subspace.

###  V.2 $\ln{W}$ for $gT\approx\pi$, $T\ll 1$: SPM phase

Let us first consider $gT\approx\pi$ where $s_{1}\ll 1$ and $c_{1}=-1$.
Truncating $W$ in Eq. (41) to $O(s_{1})$ we obtain

$\displaystyle W=\begin{pmatrix}-1&-s_{1}c_{2}&-s_{2}s_{1}&0&0&0&0&0\\\
s_{1}&-c_{2}&-s_{2}&-s_{1}s_{2}c_{2}&-s_{2}^{2}s_{1}&0&0&0\\\
0&s_{2}&-c_{2}&-s_{1}c_{2}^{2}&-s_{1}s_{2}c_{2}&0&0&0\\\
0&0&s_{1}&-c_{2}&-s_{2}&-s_{1}s_{2}c_{2}&-s_{2}^{2}s_{1}&0\\\
0&0&0&s_{2}&-c_{2}&-s_{1}c_{2}^{2}&-s_{1}s_{2}c_{2}&0\\\
0&0&0&0&s_{1}&-c_{2}&-s_{2}&-s_{1}s_{2}\\\
0&0&0&0&0&s_{2}&-c_{2}&-s_{1}c_{2}\\\
0&0&0&0&0&0&s_{1}&-1\end{pmatrix}+O(s_{1}^{2})=W_{0}+\delta W,$ (57)
$\displaystyle W_{0}=W(gT=\pi)=-\begin{pmatrix}1&0&0&0&0&0&0&0\\\
0&c_{2}&s_{2}&0&0&0&0&0\\\ 0&-s_{2}&c_{2}&0&0&0&0&0\\\
0&0&0&c_{2}&s_{2}&0&0&0\\\ 0&0&0&-s_{2}&c_{2}&0&0&0\\\
0&0&0&0&0&c_{2}&s_{2}&0\\\ 0&0&0&0&0&-s_{2}&c_{2}&0\\\
0&0&0&0&0&0&0&1\end{pmatrix},$ $\displaystyle\delta
W=\begin{pmatrix}0&-s_{1}c_{2}&-s_{2}s_{1}&0&0&0&0&0\\\
s_{1}&0&0&-s_{1}s_{2}c_{2}&-s_{2}^{2}s_{1}&0&0&0\\\
0&0&0&-s_{1}c_{2}^{2}&-s_{1}s_{2}c_{2}&0&0&0\\\
0&0&s_{1}&0&0&-s_{1}s_{2}c_{2}&-s_{2}^{2}s_{1}&0\\\
0&0&0&0&0&-s_{1}c_{2}^{2}&-s_{1}s_{2}c_{2}&0\\\
0&0&0&0&s_{1}&0&0&-s_{1}s_{2}\\\ 0&0&0&0&0&0&0&-s_{1}c_{2}\\\
0&0&0&0&0&0&s_{1}&0\end{pmatrix}.$ (58)

Above $\delta W$ is $O(s_{1})$, and represents the leading correction to $W$
around the point $gT=\pi$. If we now also assume $T\ll 1$ which is equivalent
to $s_{2}\ll 1$, then $W_{0}$ and $\delta W$ commute if terms of
$O(s_{1}s_{2}),O(s_{1}^{2}),O(s_{2}^{2})$ and higher are dropped. Performing
the expansion in Eq. (47) we find, to first order in $s_{1}$ and $T$

$\displaystyle
i\ln\left[W\right]\approx\pm\pi+\left(\begin{array}[]{cccccccc}0&is_{1}&0&0&0&0&0&0\\\
-is_{1}&0&iT&0&0&0&0&0\\\ 0&-iT&0&is_{1}&0&0&0&0\\\ 0&0&-is_{1}&0&iT&0&0&0\\\
0&0&0&-iT&0&is_{1}&0&0\\\ 0&0&0&0&-is_{1}&0&iT&0\\\ 0&0&0&0&0&-iT&0&is_{1}\\\
0&0&0&0&0&0&-is_{1}&0\\\ \end{array}\right).$ (67)

Thus we find that the Floquet Hamiltonian $i\ln{W}$ now represents a
topologically non-trivial SSH model (for $|s_{1}|<T$), with a constant shift
of $\pi$, the latter ensuring that the zero mode of the SSH model is shifted
in energy by $\pi$. A similar picture was discussed in Ref. 30, where the
discussion was presented for the Floquet Hamiltonian $i\ln{K}$.

###  V.3  $\ln{W}$ for $T\approx\pi$: SZM-SPM phase

Let us consider the case $T\approx\pi$. For this case we have
$c_{2}\approx-1,s_{2}\ll 1$ while $c_{1},s_{1}$ are arbitrary. Truncating Eq.
(41) to $O(s_{2})$ we obtain

$\displaystyle W$
$\displaystyle=\begin{pmatrix}c_{1}&s_{1}&s_{1}s_{2}c_{1}&s_{1}^{2}s_{2}&0&0&0&0\\\
s_{1}&-c_{1}&-s_{2}c_{1}^{2}&-s_{1}s_{2}c_{1}&0&0&0&0\\\
0&s_{2}&-c_{1}&-s_{1}&-s_{1}s_{2}c_{1}&-s_{1}^{2}s_{2}&0&0\\\
0&0&s_{1}&-c_{1}&-s_{2}c_{1}^{2}&-s_{1}s_{2}c_{1}&0&0\\\
0&0&0&s_{2}&-c_{1}&-s_{1}&-s_{1}s_{2}c_{1}&s_{1}^{2}s_{2}\\\
0&0&0&0&s_{1}&-c_{1}&-s_{2}c_{1}^{2}&s_{1}s_{2}c_{1}\\\
0&0&0&0&0&s_{2}&-c_{1}&s_{1}\\\
0&0&0&0&0&0&s_{1}&c_{1}\end{pmatrix}+O(s_{2}^{2})=W_{0}+\delta W,$ (68)
$\displaystyle W_{0}$
$\displaystyle=W(T=\pi)=\begin{pmatrix}c_{1}&s_{1}&0&0&0&0&0&0\\\
s_{1}&-c_{1}&0&0&0&0&0&0\\\ 0&0&-c_{1}&-s_{1}&0&0&0&0\\\
0&0&s_{1}&-c_{1}&0&0&0&0\\\ 0&0&0&0&-c_{1}&-s_{1}&0&0\\\
0&0&0&0&s_{1}&-c_{1}&0&0\\\ 0&0&0&0&0&0&-c_{1}&s_{1}\\\
0&0&0&0&0&0&s_{1}&c_{1}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}&0&0&0&0&0&0\\\ \pm ie^{\mp
i\frac{\pi}{2}(c_{1}\sigma_{z}+s_{1}\sigma_{x})}&0&0&0&0&0&0\\\
0&0&&0&0&0&0\\\ 0&0&-e^{iT\sigma_{y}}&0&0&0&0\\\ 0&0&0&0&&0&0\\\
0&0&0&0&-e^{iT\sigma_{y}}&0&0\\\ 0&0&0&0&0&0&\\\ 0&0&0&0&0&0&\pm ie^{\mp
i\frac{\pi}{2}(-c_{1}\sigma_{z}+s_{1}\sigma_{x})}\end{pmatrix},$
$\displaystyle\delta W$
$\displaystyle=\begin{pmatrix}0&0&s_{1}s_{2}c_{1}&s_{1}^{2}s_{2}&0&0&0&0\\\
0&0&-s_{2}c_{1}^{2}&-s_{1}s_{2}c_{1}&0&0&0&0\\\
0&s_{2}&0&0&-s_{1}s_{2}c_{1}&-s_{1}^{2}s_{2}&0&0\\\
0&0&0&0&-s_{2}c_{1}^{2}&-s_{1}s_{2}c_{1}&0&0\\\
0&0&0&s_{2}&0&0&-s_{1}s_{2}c_{1}&s_{1}^{2}s_{2}\\\
0&0&0&0&0&0&-s_{2}c_{1}^{2}&s_{1}s_{2}c_{1}\\\ 0&0&0&0&0&s_{2}&0&0\\\
0&0&0&0&0&0&0&0\end{pmatrix}.$ (69)

Above $\delta W$ is the leading correction to $W_{0}$ and is $O(s_{2})$.

At this stage we simply discuss the Floquet Hamiltonian $i\ln(W_{0})$,

$\displaystyle\ln\left[W_{0}\right]\approx\begin{pmatrix}&0&0&0&0&0&0\\\
\ln\biggl{[}\pm ie^{\mp
i\frac{\pi}{2}(c_{1}\sigma_{z}+s_{1}\sigma_{x})}\biggr{]}&0&0&0&0&0&0\\\
0&0&&0&0&0&0\\\ 0&0&\ln\biggl{[}-e^{iT\sigma_{y}}\biggr{]}&0&0&0&0\\\
0&0&0&0&&0&0\\\ 0&0&0&0&\ln\biggl{[}-e^{iT\sigma_{y}}\biggr{]}&0&0\\\
0&0&0&0&0&0&\\\ 0&0&0&0&0&0&\ln\biggl{[}\pm ie^{\mp
i\frac{\pi}{2}(-c_{1}\sigma_{z}+s_{1}\sigma_{x})}\biggr{]}\end{pmatrix}.$ (70)

Thus we see that $i\ln{(W_{0})}$, for these parameters, has a complex edge
structure represented by the $2\times 2$ block in the upper and lower
diagonals. The eigenvalues of these edge modes are $0,\pm\pi$. Moreover, these
edge modes are completely decoupled from a bulk which has a dimerization of a
topologically non-trivial SSH model represented by the repeated sub-block
$\ln\biggl{[}-e^{iT\sigma_{y}}\biggr{]}$. Deviating from this exactly solvable
limit will couple the edge modes weakly to the bulk states, and will also
generate longer range hopping. However, these edge modes are protected as long
as the bulk gap remains non-zero. Thus we see that the SZM-SPM phase also can
be interpreted as topologically protected edge modes of a single-particle
Hamiltonian in Krylov subspace. The advantage of the discussion in this
section is that, even when interactions are non-zero, one may construct an
effectively single-particle $W$. Thus, topological features, or lack thereof,
of $W$ can shed light on the lifetime of the almost strong modes.

##  VI Lifetime of edge modes from the Arnoldi Algorithm

In this section we construct the approximate edge mode and derive its lifetime
for the most general $W$ matrix whose upper Hessenberg form is highlighted in
Eq. (32). The results derived here hold for the free and the interacting
cases. Recall that for the free system we have a SZM and/or a SPM. These
acquire a lifetime due to finite size effects, i.e, their lifetime scales
exponentially in the system size $L$. On the other hand, when the interactions
are non-zero, the edge modes acquire a system-size independent lifetime
provided the system size is large enough. For small system sizes, there is not
much of a difference between strong modes and almost strong modes as both have
lifetimes that grow exponentially with system size.

We denote edge modes of $W$, both exact and approximate, as
$|\psi_{0,\pi}\rangle$, where the subscript indicates whether it is a zero
mode or a $\pi$-mode. We then argue that the lifetime $\gamma_{0,\pi}^{-1}$ of
the edge mode is well approximated by the expression

$\displaystyle
e^{-\gamma_{0,\pi}}=|\langle\psi_{0,\pi}|W|\psi_{0,\pi}\rangle|,$ (71)

where $\gamma_{0,\pi}$ is in dimensionless units. We will present numerical
evidence for this. It is straightforward to see that when
$|\psi_{0,\pi}\rangle$ are exact edge modes,
$|\langle\psi_{0,\pi}|W|\psi_{0,\pi}\rangle|=1$ and the decay rate
$\gamma_{0,\pi}=0$. We will apply the formula to the case where we have quasi-
stable edge modes, so that $|\langle\psi_{0,\pi}|W|\psi_{0,\pi}\rangle|<1$,
resulting in a non-zero decay rate.

We solve for the edge modes by determining the left eigenvectors of $W$,

$\displaystyle\langle\psi_{0,\pi}|W=\tilde{\lambda}\langle\psi_{0,\pi}|,$ (72)

with $\tilde{\lambda}=1$ for $\langle\psi_{0}|$ and $\tilde{\lambda}=-1$ for
$\langle\psi_{\pi}|$. Essentially this allows us to work from the first site
and iterate into the bulk. We denote the elements of $\langle\psi_{0,\pi}|$ by
$\psi_{i}$ where the subscript $i$ is the site index. The first step involves
the inner product of the left eigenvector with the first column of
$W-\tilde{\lambda}$ which gives

$\displaystyle\psi_{1}$ $\displaystyle=1,$ (73) $\displaystyle\psi_{2}$
$\displaystyle=-\biggl{(}\frac{w_{1,1}-\tilde{\lambda}}{w_{2,1}}\biggr{)}.$
(74)

Above we set $\psi_{1}=1$ as we are interested in an edge mode. The inner
product of the left eigenvector onto the second column of the matrix
$W-\tilde{\lambda}$ gives

$\displaystyle\psi_{3}$
$\displaystyle=-\frac{1}{w_{3,2}}\left[(w_{2,2}-\tilde{\lambda})\psi_{2}+w_{1,2}\psi_{1}\right].$
(75)

The pattern is clear for the rest of the eigenvector components,

$\psi_{n+1}=-\frac{1}{w_{n+1,n}}\left[(w_{n,n}-\tilde{\lambda})\psi_{n}+\sum_{i=1}^{n-1}w_{i,n}\psi_{i}\right],$
(76)

with a final normalization step at the end of the calculation.

We now rewrite Eq. (76) as follows as it will be helpful later

$\sum_{i=1}^{n-1}\psi_{i}w_{i,n}=-\psi_{n+1}w_{n+1,n}-(w_{n,n}-\tilde{\lambda})\psi_{n}.$
(77)

Since all the elements of $W$ are real, the components $\psi_{i}$ of the edge
mode are also real. We now set out to compute the right hand side of Eq. (71).
It is staightforward to see that

$\displaystyle\langle\psi_{0,\pi}|W|\psi_{0,\pi}\rangle$
$\displaystyle=\sum_{k=1}^{N-1}\left[\psi_{k+1}w_{k+1,k}+\psi_{k}w_{k,k}+\sum_{m=1}^{k-1}\psi_{m}w_{m,k}\right]\psi_{k}+\sum_{m=1}^{N}\psi_{m}w_{m,N}\psi_{N}$
(78)
$\displaystyle=\sum_{k=1}^{N-1}\left[\psi_{k+1}w_{k+1,k}+\psi_{k}w_{k,k}+\left(-\psi_{k+1}w_{k+1,k}-(w_{k,k}-\tilde{\lambda})\psi_{k}\right)\right]\psi_{k}+\sum_{m=1}^{N}\psi_{m}w_{m,N}\psi_{N}$
(79)
$\displaystyle=\sum_{k=1}^{N-1}\tilde{\lambda}\psi_{k}^{2}+\sum_{m=1}^{N}\psi_{m}w_{m,N}\psi_{N}$
(80)
$\displaystyle=\tilde{\lambda}\left[\sum_{k=1}^{N}\psi_{k}^{2}+\tilde{\lambda}\sum_{m=1}^{N}\psi_{m}w_{m,N}\psi_{N}-\psi_{N}^{2}\right].$
(81)

In the second equality above we have used Eq. (77) to simplify the expression.

We normalize $\psi_{0,\pi}\rightarrow\psi_{0,\pi}/N_{0,\pi}$ and use that
$|\tilde{\lambda}|=1$ to obtain

$|\langle\psi_{0,\pi}|W|\psi_{0,\pi}\rangle|=1+\frac{\tilde{\lambda}}{N_{0,\pi}^{2}}\sum_{m=1}^{N}\psi_{m}w_{m,N}\psi_{N}-\frac{\psi_{N}^{2}}{N_{0,\pi}^{2}}.$
(82)

Using Eq. (71), we identify the decay rate from Taylor expanding the left hand
side of the same equation to obtain

$\gamma_{0,\pi}\approx\frac{\psi_{N}^{2}}{N_{0,\pi}^{2}}-\frac{\tilde{\lambda}}{N_{0,\pi}^{2}}\sum_{m=1}^{N}\psi_{m}w_{m,N}\psi_{N}.$
(83)

We cannot proceed further without making approximations to $W$ or to $\psi$.
As it stands, Eq. (83) agrees well with the quantity
$\ln{|\langle\psi_{0,\pi}|W|\psi_{0,\pi}\rangle|}$ which was how the decay-
rate was defined in Eq. (71). We give evidence of this in the next section.
Empirically, we also find that the second term in Eq. (83) is often much
smaller than the first, thus we can approximate the expression for the decay
rate even further as

$\gamma_{0,\pi}\approx\psi_{N}^{2}.$ (84)

Above we are explicitly working with the normalized wavefunction and therefore
we have dropped the factor of $N_{0,\pi}^{2}$.

Thus the decay rate is the probability of the edge mode to be located at the
site $N$. Recall that $N$ is the dimension of the $W$ matrix. For the free
case, if $N=2L$, this exhausts the Hilbert space. For the interacting problem
$N$ has to be very large $\approx e^{2L}$ in order to exhaust the Hilbert
space. Thus when working with interacting systems, we will always be
truncating $W$. However as long as this truncation occurs for an $N$ which is
large enough, and the quasi-stable edge mode is sufficiently localized at the
edge, the decay rate $\gamma_{0,\pi}$, and therefore the probability of
finding the particle at site $N$, will be independent of $N$ (for a given
system size $L$).

##  VII Krylov subspace with Interactions

We now consider the interacting example by modifying the binary drive in Eq.
(1) into a ternary drive where the third part of the drive breaks the free
fermion nature of the problem. Specifically we study

$U=e^{-i\frac{T}{2}J_{z}H_{zz}}e^{-i\frac{T}{2}J_{x}H_{xx}}e^{-i\frac{T}{2}gH_{z}},$
(85)

where the interacting part is

$\displaystyle H_{zz}=\sum_{i=1}^{L-1}\sigma^{z}_{i}\sigma^{z}_{i+1}.$ (86)

The above drive was studied in detail in [22, 30], and it was found that the
SZMs and SPMs get modified to ASZMs and ASPMs where the latter are
characterized by a system size independent lifetime, for sufficiently large
system sizes. Moreover, this lifetime is still very long as compared to the
time needed for the bulk of the system to heat to infinite temperature. These
studies were carried out employing the Lanczos iteration scheme. Therefore, we
discuss the role of interactions only within the Arnoldi iteration scheme.

###  VII.1 Krylov subspace from the Arnoldi iteration

Section IV outlined the Arnoldi procedure. The key quantity is the unitary
matrix $W$ in Eq. (32) which has a characteristic upper Hessenberg form. The
analytic expression for $W$ and $\ln{W}$ for the binary drive was presented in
Section V, and the lifetime of the edge mode operators for the general case
including interactions was derived in Section VI. In this section we present
the $W$ matrix for the interacting problem, and apply the results for the
lifetime, derived in Section VI.

Fig. 5 shows the spectrum of $i\ln{W}$ for the free case (top panel) and for
the interacting case (bottom panel) where $J_{z}=0.05$ for the interacting
case. All other parameters are common between the top and bottom panels. In
particular, $g=0.3$, the system size is $L=10$, and $W$ is a $20\times 20$
matrix. Thus for the free case, $W$ is exact. However for the interacting
case, $W$ is not exact because the Hilbert space is larger than $2L$, and
construction of the $W$ matrix leads to truncation and loss of unitarity. A
consequence of this is the appearance of unphysical zero modes. This is
clearly seen in the lower left-most panel of Fig. 5 which shows three zero
modes, while the upper panel (free case) shows only two zero modes. Thus the
truncation for the interacting case has lead to the appearance of an
additional spurious zero mode.

To understand why spurious zero modes can appear on truncation, let us
consider a $W$ with a simple $3\times 3$ structure. The $W$ have the property
that the lower sub-diagonal (c.f. Eq. (30)) is the strongest as it measures
the part of the operator that explores new regions of the Krylov subspace. For
an ergodic system, it is natural that this element will be largest. However,
an artificial truncation will cause a $W$ of the form

$\displaystyle W=\begin{pmatrix}0&0&1\\\ 1&0&0\\\ 0&1&0\end{pmatrix}$ (87)

to become

$\displaystyle W_{\rm trunc}=\begin{pmatrix}0&0&0\\\ 1&0&0\\\
0&1&0\end{pmatrix}$ (88)

Now $W_{\rm trunc}$ has a null vector

$\displaystyle W_{\rm trunc}\begin{pmatrix}0\\\ 0\\\ 1\end{pmatrix}=0$ (89)

This null vector is not shared by the untruncated $W$. Thus quite generally,
the fact that the lower sub-diagonal is dominant in $W$ can lead to the
appearance of null vectors when the $W$ matrix is truncated, where the null
vectors have a large weight at the lower end.

Figure 5: Spectrum of $i\ln{W}$ for free example (upper row) and an
interacting example with $J_{z}=0.05$ (lower row). $x$-axis label “index”
denotes the location of the eigenvalues arranged in descending order of their
magnitude. For all data $L=10$ and $W$ is a $20\times 20$ matrix. In addition
$g=0.3$. The four different values of $T$ are from left to right
$T=2.0,3.5,5.5,8.25$. These correspond respectively to (A)SZM, (A)SZM-(A)SPM,
trivial, and (A)SPM phases. The horizonal red lines indicate $0,\pm\pi$.

We now discuss the (A)SPMs. The $\pm\pi$ modes are clearly visible on the top
panel, and they are found to persist in the presence of interactions, although
the edge modes are not so well separated in energy from the bulk states when
interactions are present. This is expected as we now have ASPMs which will now
decay into the bulk. Sometimes the mode at $-\pi$ can appear at $\pi$ as seen
in the lower fourth panel. This is not a spurious effect because $-\pi$ and
$\pi$ states are degenerate states due to the periodicity of the spectrum.
However the truncation can lead to spurious effects such as the disappearance
of one of the $\pi$ modes as can be seen in the lower second panel.

Figure 6: Time evolution of approximate edge states of $W$ (orange) for an
interacting example ($J_{z}=0.05$) and compared to $A_{\infty}$ obtained from
ED (blue). The top row is for a very small system size of $L=4$ and the bottom
row is for $L=8$. The four columns compare four different phases. These are
from left to right: ASZM phase ($T=2.0$), ASZM-ASPM phase ($T=3.5$), trivial
phase ($T=5.5$), and ASPM phase ($T=8.25$). $W$ is truncated to a 18 $\times$
18 matrix, and the time-evolution generated by it (orange) is according to Eq.
(91). The data labeled by $d_{1,\lambda},d_{2,\lambda}$ with $\lambda=0,\pi$
are two different approximations corresponding to Eq. (92) for $d_{1,\lambda}$
and Eq. (93) for $d_{2,\lambda}$.

We now turn to the computation of $A_{\infty}$ and extracting the lifetime.
Fig. 6 shows the autocorrelation functions for $g=0.3$, $J_{z}=0.05$ and for
the same set of $T$ as Fig. 5. The top row is for $L=4$ while the bottom row
is for $L=8$. The blue dots are $A_{\infty}$ from ED. The orange data are the
$A_{\infty}$ obtained from the $W$ matrix, with the $W$ matrix truncated to
$N=18$. We have checked that the truncation at $N=18$ does not influence the
results for the system sizes chosen ($L=4,8$).

The precise quantity being plotted in orange (labeled $W$) is Eq. (33), but
with some approximations made to it. In particular,

$\displaystyle A_{\infty}(nT)$ $\displaystyle={\rm
Tr}\biggl{[}W^{n}|1\rangle\langle 1\biggr{]}$ $\displaystyle\approx{\rm
Tr}\biggl{[}W^{n}\biggl{(}|\psi_{0}\rangle\langle\psi_{0}|+|\psi_{\pi}\rangle\langle\psi_{\pi}|\biggr{)}|1\rangle\langle
1\biggr{]}$ (90) $\displaystyle=\sum_{\lambda=0,\pi}\psi_{\lambda,1}\langle
1|W^{n}|\psi_{\lambda}\rangle,\,\,\psi_{\lambda,1}=\langle\psi_{\lambda}|1\rangle.$
(91)

Above $|1\rangle$ is a state which is completely localized on the first site
of the Krylov subspace. The approximation made in the second line involves
replacing the complete set of states of $W$ by only its edge modes
$|\psi_{\lambda=0,\pi}\rangle$. Thus we are dropping all the bulk modes of
$W$. Also note that, $\psi_{\lambda,1}$ is the normalized amplitude of the
edge mode on the first site.

The results for $A_{\infty}$ in Fig. 6 are compared to two different
approximations. One is the data labeled as $d_{1,\lambda=0,\pi}$ where

$\displaystyle d_{1,\lambda}(nT)$
$\displaystyle=|\psi_{\lambda,1}|^{2}e^{-\gamma_{\lambda}n}$
$\displaystyle=|\psi_{\lambda,1}|^{2}|\langle\psi_{\lambda}|W|\psi_{\lambda}\rangle|^{n},\,\,\lambda=0,\pi.$
(92)

The above expression uses the definition of the decay-rate in Eq. (71)
accounting for the amplitude $\psi_{\lambda,1}$ of finding the particle on the
first site, before the time-evolution. Since $|\psi_{0,\pi}\rangle$ edge modes
are only approximate edge modes, Eq. (92) does lead to a decay. The data set
labeled by $d_{2,\lambda=0,\pi}$ corresponds to approximating $\gamma_{0,\pi}$
by Eq. (84), thus

$\displaystyle
d_{2,\lambda}(nT)=|\psi_{\lambda,1}|^{2}e^{-|\psi_{\lambda,N}|^{2}n},\,\,\lambda=0,\pi,$
(93)

where $\psi_{\lambda,N}$ is the normalized amplitude of the $\lambda=0,\pi$
approximate edge mode at the last site $N$.

We find that $W$, despite its truncation, agrees very well with ED. In
addition, the two approximations to $W$ given by Eq. (92) and Eq. (93) also
agree well as far as capturing the decay rates are concerned. The
approximations appear to not work quite so well for the second panel
corresponding to the ASZM-ASPM phase ($T=3.5$). This may be because when both
almost strong modes are present, as they decay they influence each other
through the bulk states in a way that our simple approximation has not
accounted for. However, the agreement does improve with increasing system size
(compare upper and lower panels for $T=3.5$). In particular, when both ASZM
and ASPM are present, the edge structure is more complex and has a $2\times 2$
structure at each end (see Section V). This can make the effective system size
of $L=4$ (upper panels) look effectively shorter for an ASZM-ASPM phase as
compared to an ASPM or an ASZM phase.

##  VIII Conclusions

Almost strong edge modes are quasi-stable edge modes that have unusually long
lifetimes such that they can coexist with a metallic bulk for many drive
cycles. These edge modes include both zero modes as well as $\pi$ modes, where
the latter show period doubled dynamics. It is notoriously hard to develop
analytic methods for interacting and driven systems in general, and
determining lifetimes of quasi-stable modes in particular. However in this
paper, building on previous work [16, 17, 30], we showed a promising route
which involves mapping the dynamics of the edge mode operator to single
particle dynamics in Krylov subspace. While the detailed modeling of the
Krylov subspace, such as the precise values of the hopping parameters, still
requires the same computational costs as ED, yet when the operator of interest
has some universal features, the Krylov subspace can have some general
properties. For example, the Krylov subspace of maximally chaotic systems have
certain universal features much discussed in the literature [37, 38, 39, 40].

In this paper we showed that when the operators are strong modes or almost
strong modes, they appear as stable or quasi-stable edge modes of Krylov
subspaces with topological properties. For the examples studied in this paper,
the Krylov subspaces are given by generalized SSH models with long range and
spatially inhomogeneous hopping. Exploiting these topological structures can
lead to better understanding of the long lifetimes of the almost strong modes.

Most studies on Krylov subspace dynamics use the Lanczos method where the
generator of dynamics is a static Hamiltonian [36, 37, 38, 39, 40]. In
studying Floquet systems, as we do in this paper, one has to adapt the Krylov
subspace because the generator of dynamics is a unitary operator rather than a
Hermitian operator. This lead us to derive the Krylov subspace using an
alternate approach, known as the the Arnoldi iteration [31]. Topological
features of the resulting chain were highlighted. A compact expression (c.f.
Eq. (93)) for the lifetime of the edge modes was derived and compared with ED.
In particular, we showed that the decay rate is simply determined by the
probability of finding the particle on the last site $N$ of the $W$ matrix
that generates the time-evolution in this Krylov subspace. This observation
opens up the possibility of developing analytic tools for calculating the
decay rate, such as by employing the WKB approximation.

Our methods can be generalized to study other kinds of slow dynamics such as
dynamics of scar states [42, 43, 44, 45]. It is also promising to perform a
topological classification of the Krylov subspaces of edge modes of strongly
interacting topological insulators, both static and Floquet [46, 47, 48].

Acknowledgements: The authors thank Sasha Abanov for helpful discussions. This
work was supported by the US Department of Energy, Office of Science, Basic
Energy Sciences, under Award No. DE-SC0010821.

## References

  * Kitaev [2001] A. Y. Kitaev, Unpaired majorana fermions in quantum wires, Physics-Uspekhi 44, 131 (2001).
  * Sachdev [2011] S. Sachdev, Quantum Phase Transitions, Cambridge University Press (2011).
  * Bernevig and [with T. L. Hughes] B. A. Bernevig and (with T. L. Hughes), Topological Insulator and Topological Superconductors, Princeton University Press, Princeton (2013).
  * Kitaev [2006] A. Kitaev, Anyons in an exactly solved model and beyond, Annals of Physics 321, 2 (2006), january Special Issue.
  * Nayak  _et al._ [2008] C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, Non-abelian anyons and topological quantum computation, Rev. Mod. Phys. 80, 1083 (2008).
  * Alicea [2012] J. Alicea, New directions in the pursuit of majorana fermions in solid state systems, Reports on Progress in Physics 75, 076501 (2012).
  * Sarma  _et al._ [2015] S. D. Sarma, M. Freedman, and C. Nayak, Majorana zero modes and topological quantum computation, njp Quantum Information 1, 15001 (2015).
  * Sacha and Zakrzewski [2017] K. Sacha and J. Zakrzewski, Time crystals: a review, Reports on Progress in Physics 81, 016401 (2017).
  * Else  _et al._ [2020] D. V. Else, C. Monroe, C. Nayak, and N. Y. Yao, Discrete time crystals, Annual Review of Condensed Matter Physics 11, 467 (2020).
  * Khemani  _et al._ [2019] V. Khemani, R. Moessner, and S. Sondhi, A brief history of time crystals, arXiv:1910.10745 (2019).
  * Fendley [2016] P. Fendley, Strong zero modes and eigenstate phase transitions in the xyz/interacting majorana chain, Journal of Physics A: Mathematical and Theoretical 49, 30LT01 (2016).
  * Else  _et al._ [2017] D. V. Else, P. Fendley, J. Kemp, and C. Nayak, Prethermal strong zero modes and topological qubits, Phys. Rev. X 7, 041062 (2017).
  * Kemp  _et al._ [2017] J. Kemp, N. Y. Yao, C. R. Laumann, and P. Fendley, Long coherence times for edge spins, Journal of Statistical Mechanics: Theory and Experiment 2017, 063105 (2017).
  * Parker  _et al._ [2019a] D. E. Parker, R. Vasseur, and T. Scaffidi, Topologically protected long edge coherence times in symmetry-broken phases, Phys. Rev. Lett. 122, 240605 (2019a).
  * Kemp  _et al._ [2020] J. Kemp, N. Y. Yao, and C. R. Laumann, Symmetry-enhanced boundary qubits at infinite temperature, Phys. Rev. Lett. 125, 200506 (2020).
  * Yates  _et al._ [2020a] D. J. Yates, A. G. Abanov, and A. Mitra, Lifetime of almost strong edge-mode operators in one-dimensional, interacting, symmetry protected topological phases, Phys. Rev. Lett. 124, 206803 (2020a).
  * Yates  _et al._ [2020b] D. J. Yates, A. G. Abanov, and A. Mitra, Dynamics of almost strong edge modes in spin chains away from integrability, Phys. Rev. B 102, 195419 (2020b).
  * Jiang  _et al._ [2011] L. Jiang, T. Kitagawa, J. Alicea, A. R. Akhmerov, D. Pekker, G. Refael, J. I. Cirac, E. Demler, M. D. Lukin, and P. Zoller, Majorana fermions in equilibrium and in driven cold-atom quantum wires, Phys. Rev. Lett. 106, 220402 (2011).
  * Bastidas  _et al._ [2012] V. M. Bastidas, C. Emary, G. Schaller, and T. Brandes, Nonequilibrium quantum phase transitions in the ising model, Phys. Rev. A 86, 063627 (2012).
  * Thakurathi  _et al._ [2013] M. Thakurathi, A. A. Patel, D. Sen, and A. Dutta, Floquet generation of majorana end modes and topological invariants, Phys. Rev. B 88, 155133 (2013).
  * Asbóth  _et al._ [2014] J. K. Asbóth, B. Tarasinski, and P. Delplace, Chiral symmetry and bulk-boundary correspondence in periodically driven one-dimensional systems, Phys. Rev. B 90, 125143 (2014).
  * Yates  _et al._ [2019] D. J. Yates, F. H. L. Essler, and A. Mitra, Almost strong ($0,\pi$) edge modes in clean interacting one-dimensional floquet systems, Phys. Rev. B 99, 205419 (2019).
  * Khemani  _et al._ [2016] V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi, Phase structure of driven quantum systems, Phys. Rev. Lett. 116, 250401 (2016).
  * von Keyserlingk  _et al._ [2016] C. W. von Keyserlingk, V. Khemani, and S. L. Sondhi, Absolute stability and spatiotemporal long-range order in floquet systems, Phys. Rev. B 94, 085112 (2016).
  * von Keyserlingk and Sondhi [2016a] C. W. von Keyserlingk and S. L. Sondhi, Phase structure of one-dimensional interacting floquet systems. i. abelian symmetry-protected topological phases, Phys. Rev. B 93, 245145 (2016a).
  * von Keyserlingk and Sondhi [2016b] C. W. von Keyserlingk and S. L. Sondhi, Phase structure of one-dimensional interacting floquet systems. ii. symmetry-broken phases, Phys. Rev. B 93, 245146 (2016b).
  * Katsura  _et al._ [2015] H. Katsura, D. Schuricht, and M. Takahashi, Exact ground states and topological order in interacting kitaev/majorana chains, Phys. Rev. B 92, 115137 (2015).
  * Wouters  _et al._ [2018] J. Wouters, H. Katsura, and D. Schuricht, Exact ground states for interacting kitaev chains, Phys. Rev. B 98, 155119 (2018).
  * Wada  _et al._ [2021] K. Wada, T. Sugimoto, and T. Tohyama, Coexistence of strong and weak majorana zero modes in an anisotropic xy spin chain with second-neighbor interactions, Phys. Rev. B 104, 075119 (2021).
  * Yates  _et al._ [2021] D. J. Yates, A. G. Abanov, and A. Mitra, Long-lived $\pi$ edge modes of interacting and disorder-free floquet spin chains, arXiv:2021.13766 (2021).
  * Arnoldi [1951] W. E. Arnoldi, The principle of minimized iterations in the solution of the matrix eigenvalue problem, Quart. Appl. Math. 9, 17 (1951).
  * Prosen [2000] T. Prosen, Exact time-correlation functions of quantum ising chain in a kicking transversal magnetic field: Spectral analysis of the adjoint propagator in heisenberg picture, Progress of Theoretical Physics Supplement 139 (2000).
  * Gritsev and Polkovnikov [2017] V. Gritsev and A. Polkovnikov, Integrable Floquet dynamics, SciPost Phys. 2, 021 (2017).
  * Su  _et al._ [1979] W. P. Su, J. R. Schrieffer, and A. J. Heeger, Solitons in polyacetylene, Phys. Rev. Lett. 42, 1698 (1979).
  * Su  _et al._ [1980] W. P. Su, J. R. Schrieffer, and A. J. Heeger, Soliton excitations in polyacetylene, Phys. Rev. B 22, 2099 (1980).
  * Vishwanath and Müller [2008] V. Vishwanath and G. Müller, The Recursion Method: Applications to Many-Body Dynamics, Springer, New York (2008).
  * Parker  _et al._ [2019b] D. E. Parker, X. Cao, A. Avdoshkin, T. Scaffidi, and E. Altman, A universal operator growth hypothesis, Phys. Rev. X 9, 041017 (2019b).
  * Dymarsky and Gorsky [2020] A. Dymarsky and A. Gorsky, Quantum chaos as delocalization in krylov space, Phys. Rev. B 102, 085137 (2020).
  * Barbón  _et al._ [2019] J. Barbón, E. Rabinovici, R. Shir, and R. Sinha, On the evolution of operator complexity beyond scrambling, Journal of High Energy Physics 2019, 264 (2019).
  * Avdoshkin and Dymarsky [2020] A. Avdoshkin and A. Dymarsky, Euclidean operator growth and quantum chaos, Phys. Rev. Research 2, 043234 (2020).
  * Horn and Johnson [2013] R. A. Horn and C. R. Johnson, Matrix Analysis, Cambridge University Press, 2nd Edition (2013).
  * Moudgalya  _et al._ [2020] S. Moudgalya, E. O’Brien, B. A. Bernevig, P. Fendley, and N. Regnault, Large classes of quantum scarred hamiltonians from matrix product states, Phys. Rev. B 102, 085120 (2020).
  * Turner  _et al._ [2018] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić, Weak ergodicity breaking from quantum many-body scars, Nature Physics 14, 745 (2018).
  * Ho  _et al._ [2019] W. W. Ho, S. Choi, H. Pichler, and M. D. Lukin, Periodic orbits, entanglement, and quantum many-body scars in constrained models: Matrix product state approach, Phys. Rev. Lett. 122, 040603 (2019).
  * Lin and Motrunich [2019] C.-J. Lin and O. I. Motrunich, Exact quantum many-body scar states in the rydberg-blockaded atom chain, Phys. Rev. Lett. 122, 173401 (2019).
  * Fidkowski and Kitaev [2010] L. Fidkowski and A. Kitaev, Effects of interactions on the topological classification of free fermion systems, Phys. Rev. B 81, 134509 (2010).
  * Potter  _et al._ [2016] A. C. Potter, T. Morimoto, and A. Vishwanath, Classification of interacting topological floquet phases in one dimension, Phys. Rev. X 6, 041001 (2016).
  * Else and Nayak [2016] D. V. Else and C. Nayak, Classification of topological phases in periodically driven interacting systems, Phys. Rev. B 93, 201103 (2016).