
\section{Floquet topological phase of the 1D ladder}

The 1D ladder model introduced in the main text is an example of an anomalous Floquet topological phase, one in which topologically protected end modes appear even though the bulk bands have vanishing topological invariants. To see this, we write the Hamiltonian and the Floquet operator in momentum space by ``straightening out'' the ladder, as shown in Fig.~\ref{fig:straighten}.

\begin{figure}[h]
 \includegraphics[width=0.4\textwidth]{straighten.pdf}
 \caption{For an infinite 1D ladder (top), we can determine the momentum-space Hamiltonian by rearranging the site doublets denoted by $j$ such that they become the unit cells of an infinite chain (bottom). The hoppings on the rungs of the ladder (blue) then connect neighboring unit cells. Half of the hoppings on the legs of the ladder (green) connect sites within a unit cell, while the other half (orange) become next-nearest-neighbor hoppings. \label{fig:straighten}}
\end{figure}

The resulting $2\times2$ momentum-space Hamitonian reads
\begin{equation}\label{eq:bulkH}
 H(k, t)=\begin{pmatrix}
    0 & q(k, t) \\
    q^*(k, t) & 0 \\
   \end{pmatrix},
\end{equation}
where $q(k, t)=J^4(t)+J^{1,3}(t) e^{i k}+J^2(t) e^{2ik}$. As in the main text, we denote by $J^{1,3}$ the hoppings on the rungs of the ladder, by $J^4$ the hopping within doublets, and by $J^2$ the one across doublets. 

Before imposing a periodic modulation of the $J^i$, we examine the symmetries of the Hamiltonian Eq.~\eqref{eq:bulkH} in the time-independent case. Since we consider real hopping amplitudes, the Hamiltonian Eq.~\eqref{eq:bulkH} exhibits time-reversal ${\cal T}$, particle-hole ${\cal P}$, as well as sublattice $\Gamma$ symmetries. Introducing Pauli matrices $\sigma_{x, y, z}$ acting on the sublattice degree of
freedom, the symmetry operators are $\Gamma = \sigma_z$, ${\cal T} = {\cal K}$, and ${\cal P} = \sigma_z {\cal K}$, with ${\cal K}$ complex conjugation, such that
\begin{equation}\label{eq:bulkHsym}
\begin{split}
 H(k) & = H^*(-k) \\
 H(k) & = - \sigma_z H^*(-k) \sigma_z \\
 H(k) & =  -\sigma_z H(k) \sigma_z \\
\end{split}
\end{equation}
 Notice
that ${\cal T}^2={\cal P}^2=+1$, such that $H(k)$ belongs to class BDI in the Altland-Zirnbauer classification \cite{Altland1997}.

Next, we consider the case of time-periodic hopping amplitudes. As discussed in the main text, during each of the four steps of the driving protocol only
one of the hoppings is active, such that the hopping $J^i=J$ for $(i-1)T/4 \leq t < iT/4$ and all other hoppings vanish.
Owing to the simple form of the Hamiltonian and of the driving protocol, the Floquet operator takes the form
\begin{equation}\label{eq:Floquetbulk}
U_{0\rightarrow T}(k)=U_F(k) = e^{-i \frac{T}{4} H_4(k)}e^{-i \frac{T}{4} H_3(k)}e^{-i \frac{T}{4} H_2(k)}e^{-i \frac{T}{4} H_1(k)},
\end{equation}
where $U_{0\rightarrow T}(k)$ denotes the time-evolution operator from time $t=0$ to $t=T$. The Floquet operator
can be computed exactly, although its form is involved for generic parameter values. At the resonant driving point however, when the strength of the active hopping $J$ is such that $JT/4=\pi/2$, the bulk Floquet operator becomes the identity operator $U_F(k)=\mathds{1}$. As such, the bulk Floquet bands are flat and positioned at zero quasienergy. In addition, this form of $U_F$ explains the anomalous nature of the topological phase: since Floquet eigenstates become momentum-independent at the resonant driving point, they must be topologically trivial. Away from the fine-tuned value $JT/4=\pi/2$, bulk bands are dispersive and show a parabolic band touching point at $k=\pi$, as shown in Fig.~\ref{fig:bulkbands}.
\begin{figure}[h]
 \includegraphics[width=0.3\textwidth]{bulkbands.pdf}
 \caption{Bandstructure of the Floquet operator Eq.~\eqref{eq:Floquetbulk} for $JT/4=\pi/2$ (dashed) and $JT/4=\pi/3$ (solid). At the resonant driving point the bulk bands are flat, whereas they disperse otherwise and form a parabolic touching point at $k=\pi$ and $\varepsilon=0$.\label{fig:bulkbands}}
\end{figure}

As discussed in the main text, the 1D ladder hosts a topological Floquet phase, such that end states appear in a finite system despite the topologically trivial bands. The presence of end states at quasienergies $\varepsilon=\pm\omega/2$ can be deduced directly from the real-space Floquet system at the resonant driving point, but remain pinned to these quasienergies also away from it. This behavior is linked to the symmetries of the Floquet operator. To determine them, it is convenient to choose a different origin of time when computing the Floquet operator and go to a so called \emph{symmetric time frame} by defining:
\begin{equation}\label{eq:Floquetsym}
U_{\frac{3T}{8}\rightarrow T+\frac{3T}{8}}(k)=\widetilde{U}_F(k) = e^{-i \frac{T}{8} H_2(k)}e^{-i \frac{T}{4} H_1(k)}e^{-i \frac{T}{4} H_4(k)}e^{-i \frac{T}{4} H_3(k)}e^{-i \frac{T}{8} H_2(k)}.
\end{equation}

For the time-shifted Floquet operator, time-reversal, particle-hole, and sublattice symmetry take the respective forms:
\begin{gather}\label{eq:FTRS}
 {\widetilde{U}}^{\phantom\dag}_F(k) = {\widetilde{U}}^T_F(-k), \\
\label{eq:FPHS}
 {\widetilde{U}}^{\phantom\dag}_F(k) = \sigma_z {\widetilde{U}}^*_F(-k) \sigma_z, \\
\label{eq:FSLS}
 {\widetilde{U}}^{\phantom\dag}_F(k) = \sigma_z {\widetilde{U}}^\dag_F(k) \sigma_z,
\end{gather}
where the superscript $T$ denotes transposition.
Note that out of these three constraints only particle-hole symmetry does not reverse the order of the product in Eq.~\eqref{eq:Floquetsym}, since it relates the Floquet operator to its complex conjugate. This means that the same constraint Eq.~\eqref{eq:FPHS} holds for all choices of initial time, and in particular also for the Floquet operator $U_F(k)$ of Eq.~\eqref{eq:Floquetbulk}, while time-reversal and sublattice symmetries acquire complicated, $k$-dependent forms.
As such, we will focus in the following on the un-shifted Floquet operator $U_F(k)$ and on the role played by particle-hole symmetry in protecting its end states.
Due to this symmetry, momentum eigenstates must come in pairs, related by $\varepsilon\to-\varepsilon$ and $k\to-k$, as seen in Fig.~\ref{fig:bulkbands}. In real space, particle-hole symmetry relates states at the same position but opposite values of quasienergy, explaining the robustness of the end modes. Even if the system is perturbed away from resonant driving, the end modes cannot couple due to the bulk gap so they cannot shift in quasienergy away from the particle-hole symmetric $\varepsilon=\omega/2=-\omega/2$.

Since the driven system shows particle-hole, time-reversal, as well as sublattice symmetries, it belongs to symmetry class BDI \cite{Altland1997}, which was shown to have a $\mathbb{Z}$ topological classification in one dimension. To determine the topological invariant responsible for the presence of end states, we use the method of Ref.~\cite{Jiang2011}, which is based on evaluating the time-evolution operator at all times during the driving cycle. Writing
\begin{equation}
 U_{0\rightarrow t}(k) = \mathsf{T} \exp \left( {-i \int_0^t dt' \, H(k,t')} \right),
\end{equation}
the topological invariant can be determined as the total number of times the eigenphases of both $U_{0\rightarrow t}(k=0)$ and $U_{0\rightarrow t}(k=\pi)$ cross $\pi$ in the interval $t\in[0,T]$. Due to the simple form of the piecewise constant Hamiltonian Eq.~\eqref{eq:bulkH}, these eigenphases can be evaluated exactly, since $H(k=0)=(J^4+J^{1,3}+J^2)\sigma_x$ and $H(k=\pi)=(J^4-J^{1,3}+J^2)\sigma_x$. Therefore, the two eigenphases of $U_{0\rightarrow t}(k=0,\pi)$ are opposite and increase or decrease linearly in each driving step, depending on whether $J^i$ enters the Hamiltonian with a positive or negative sign. At resonant driving, for $k=0$ the eigenphases are monotonic as a function of $t$ and show a single crossing at $\pi$, when $t=T/2$. For $k=\pi$, the eigenphases reverse direction in steps 1 and 3 as compared to steps 2 and 4, such that they never cross $\pi$. As such, the total number of $\pi$-crossings is nonzero, as shown in Fig.~\ref{fig:topcrossing}, implying a topologically non-trivial phase.
\begin{figure}[h]
 \includegraphics[width=0.6\textwidth]{topcrossing.pdf}
 \caption{Eigenphases of $U_{0\rightarrow t}(k=0)$ (left) and $U_{0\rightarrow t}(k=\pi)$ (right) as $t$ is advanced from $0$ to $T$, using $JT/4=\pi/2$. Pairs of opposite eigenphases are denoted by solid and dashed lines. On the left, eigenphases are monotonic, since all $J^i$ enter the Hamiltonian with the same sign for $k=0$. On the right, the eigenphases reverse direction in each driving step, since $J^{1,3}$ and $J^{2,4}$ have opposite signs for $k=\pi$. The total number of $\pi$-crossings is odd, signaling a topologically non-trivial phase.\label{fig:topcrossing}}
\end{figure}

\section{Timing noise in Floquet systems with piecewise constant driving}
 \label{sec:timing-noise-floquet-pc-drive}

Here we derive the discrete-time Floquet-Lindblad equation (FLE) that describes the
time evolution of periodically driven systems with timing noise. First we
consider the simplest case of binary driving (or, equivalently, a system
subject to periodic kicks), and then generalize our results to driving cycles
comprising an arbitrary number of steps. We discuss the conditions under
which the FLE can be reduced to a classical master equation. Then, we apply
the general formalism to the model system considered in the main text.
 
\subsection{Binary piecewise constant driving}
\label{sec:binary-pc-driv}

Let us begin with the simplest case, which is a binary and piecewise constant
driving protocol. We first discuss perfectly periodic driving, and then
modifications due to timing noise. Without noise, a driving cycle consists of
the following two steps: first, the Hamiltonian $H_1$ is applied for a time
$T_1$, and then another Hamiltonian $H_2$ is applied for a time $T_2$. The full
duration of the driving cycle is thus $T = T_1 + T_2$, and the Floquet operator
$U_F$, which describes the evolution of the system during one period, is given
by
\begin{equation}
  \label{eq:UF_binary}
  U_F = \mathsf{T} e^{-i \int_0^T dt \, H(t)} = U_2 U_1 = e^{-i T_2 H_2} e^{-i
    T_1 H_1}.
\end{equation}
Here, $\mathsf{T}$ denotes time ordering, and $H(t)$ is the time-dependent
Hamiltonian 
\begin{equation}
  H(t) =
  \begin{cases}
    H_1, & n T \leq t < n T + T_1, \\
    H_2, & n T + T_1 \leq t < \left( n + 1 \right) T.
  \end{cases}
\end{equation}
The integer $n$ counts the number of driving cycles. A Floquet operator with the
same form as in Eq.~\eqref{eq:UF_binary} arises in periodically ``kicked''
systems, in which $H_2$ is applied as an instantaneous pulse at multiples of the
period $T$. Such a scenario is described by the following time-dependent
Hamiltonian:
\begin{equation}
  \label{eq:H_kicked}
  H(t) = H_1 + \lambda \sum_{n \in \mathbbm{N}} \delta(t - n T) H_2.
\end{equation}
For this type of driving, the noiseless evolution is given by the Floquet
operator $U_F = e^{-i \lambda H_2} e^{-i T H_1}$. In both of these driving
schemes, the state of the system at multiples of the driving period,
$\ket{\psi_n} = \ket{\psi(n T)}$, can be obtained by repeated application of
$U_F$, i.e., $\ket{\psi_n} = U_F^n \ket{\psi_0}$, where $\ket{\psi_0}$ is the
initial state of the system.

We now proceed to discuss how the above driving protocols are modified in the
presence of timing noise. For simplicity, we assume that there is noise only in
the first step of the driving cycle. The generalization of our considerations to
include noise in the second step or to a driving protocol that comprises more
than two steps is straightforward and summarized in the next section. If there
is timing noise in the first step, the Hamiltonian $H_1$ is not applied exactly
for a time $T_1$. Instead, in the $n$-th cycle, it is applied for
$T_1 + \tau_n$, where $\tau_n$ is a random number with zero mean,
$\overline{\tau_n} = 0$ and fluctuations $\overline{\tau_n^2} = \tau^2$.
Strictly speaking, causality requires that $\tau_n \geq -T_1$, but formally we
can relax this constraint, and we simply require that $\tau \ll T_1$. (In
particular, this allows us to take the distribution of the $\tau_n$ to be
Gaussian as in the numerical results presented in the main text.) The precise
form of the probability distribution is not important in the
following. Moreover, we assume that the time shifts in different cycles are
uncorrelated, i.e., $\overline{\tau_n \tau_{n'}} = \tau^2 \delta_{nn'}$.  Under
these conditions, after $n$ driving cycles the state of the system is given by
\begin{equation}
  \ket{\psi_n} = e^{-i T_2 H_2} e^{-i \left( T_1 + \tau_n \right) H_1} \dotsb
  e^{-i T_2 H_2} e^{-i \left( T_1 + \tau_1 \right) H_1} \ket{\psi_0} = U_{F, n}
  \dotsb U_{F, 1} \ket{\psi_0}.
\end{equation}
In the last equality, we introduced the noisy Floquet operator $U_{F, n}$, which
describes the evolution of the system during the $n$-th driving cycle. It can be
written as
\begin{equation}
  \label{eq:UF_binary_noisy}
  U_{F, n} = \mathsf{T} e^{-i \int_{t_{n - 1}}^{t_n} dt \, H(t)} =
  e^{-i T_2 H_2} e^{-i \left( T_1 + \tau_n \right) H_1},
\end{equation}
where $t_n = n T + \sum_{n' = 1}^n \tau_n$ depends on all prior time shifts. The
time-dependent Hamiltonian is now given by
\begin{equation}
  \label{eq:noisy_Hamiltonian}
  H(t) =
  \begin{cases}
    H_1, & t_n \leq t < t_n + T_1 + \tau_n, \\
    H_2, & t_n + T_1 + \tau_n \leq t < t_n + T_1 + \tau_n + T_2.
  \end{cases}
\end{equation}
For the kicking protocol defined by Eq.~\eqref{eq:H_kicked}, the noisy
Hamiltonian assumes the form
\begin{equation}
  \label{eq:H_kicked_noise}
  H(t) = H_1 + \lambda \sum_n \delta(t - t_n) H_2,
\end{equation}
where $t_n$ is defined as above, i.e., the duration between two kicks is
$t_{n + 1} - t_n = T + \tau_n$. The corresponding Floquet operator reads as in
Eq.~\eqref{eq:UF_binary_noisy}, only with $U_2 = e^{-i T_2 H_2}$ replaced by
$e^{-i \lambda H_2}$. In both cases, the evolution of the state
$\ket{\psi_n} = \ket{\psi(t_n)}$ during one driving cycle is given by
\begin{equation}
  \label{eq:psi_evolution}
  \ket{\psi_{n + 1}} = U_{F,n + 1} \ket{\psi_n} = U_2 e^{-i \left( T_1 +
      \tau_{n + 1} \right) H_1} \ket{\psi_n}.
\end{equation}
Evidently, the state of the system at time $t_n$ depends on the particular noise
realization, i.e., on the sequence of all prior time shifts
$\tau_1, \dotsc, \tau_n$. Taking the average over noise realizations, the
expectation value of an observable $O$ can be written as
\begin{equation}
  \label{eq:O}
  \overline{\langle O_n \rangle} = \overline{\braket{\psi_n | O | \psi_n}} = \mathop{\mathrm{tr}}
  \! \left( O \overline{\ket{\psi_n} \bra{\psi_n}} \right) = \mathop{\mathrm{tr}}(O \rho_n),
\end{equation}
where $\rho_n = \overline{\ket{\psi_n} \bra{\psi_n}}$ is the density matrix that
describes the noise-averaged state of the system at time $t_n$. Thus, to
evaluate the expectation values of observables, it is sufficient to track the
evolution of $\rho_n$, and in the following we derive an evolution equation for
this quantity in the weak-noise limit. A key point that facilitates this
derivation is that the noise average can be performed for each driving cycle
individually because the $\tau_1, \dotsc, \tau_n$ are statistically
independent. Thus, rewriting Eq.~\eqref{eq:psi_evolution} for the density
matrix, we obtain
\begin{equation}
  \label{eq:rho_evolution}
  \rho_{n + 1} = \overline{U_{F, n + 1} \rho_n U_{F, n + 1}^{\dagger}}.
\end{equation}
For weak noise in the sense that $\tau \norm{H_1} \ll 1$ where
$\norm{\makebox[1ex]{$\cdot$}}$ is a suitably chosen operator norm, we can
expand the noisy Floquet operator~\eqref{eq:UF_binary_noisy} in the time
shift. Keeping terms up to second order in $\tau_{n + 1}$, we find
\begin{equation}
  \label{eq:UF_expansion}
  U_{F, n + 1} = U_F \left( 1 - i \tau_{n + 1} H_1 - \frac{\tau_{n + 1}^2}{2} H_1^2 \right),
\end{equation}
where $U_F = U_2 U_1$ is the Floquet operator for perfectly periodic
driving. Inserting this form in Eq.~\eqref{eq:rho_evolution}, the average over
noise becomes straightforward, and we obtain (dropping terms of cubic and higher
order in $\tau_{n + 1}$)
\begin{equation}
  \label{eq:discrete_Lindblad_binary}
  \begin{split}
    \rho_{n + 1} & = U_F \overline{\left( 1 - i \tau_{n + 1} H_1 - \frac{\tau_{n
            + 1}^2}{2} H_1^2 \right) \rho_n \left( 1 + i \tau_{n + 1} H_1 -
        \frac{\tau_{n + 1}^2}{2} H_1^2 \right)} U_F^{\dagger} \\ & = U_F
    \overline{\left[ \rho_n - i \tau_{n + 1} \left[ H_1, \rho_n \right] +
        \tau_{n + 1}^2 \left( H_1 \rho_n H_1 - \frac{1}{2} \left\{ H_1^2, \rho_n
          \right\} \right) \right]} U_F^{\dagger} \\ & = U_F \left[ \rho_n - i
      \overline{\tau_{n + 1}} \left[ H_1, \rho_n \right] + \overline{\tau_{n +
          1}^2} \left( H_1 \rho_n H_1 - \frac{1}{2} \left\{ H_1^2, \rho_n
        \right\} \right) \right] U_F^{\dagger} \\ & = U_F \left( \rho_n + \tau^2
      \mathcal{D}[L_1] \rho_n \right) U_F^{\dagger}.
  \end{split}
\end{equation}
This discrete-time evolution equation combines the usual noise-free coherent
stroboscopic Floquet evolution, described by the Floquet operator $U_F$, with
dissipative dynamics as familiar from quantum master equations in Lindblad
form. In particular, in analogy to the usual continuous-time Lindblad equation,
we identify the quantum jump operator $L_1 = H_1$ and the dissipator
\begin{equation}
  \label{eq:dissipator}
  \mathcal{D}[L] \rho = L \rho L - \frac{1}{2} \left\{ L^2, \rho \right\} =
  \frac{1}{2} \left[ \left[ L, \rho \right], L \right].
\end{equation}
In the Floquet-Lindblad equation (FLE)~\eqref{eq:discrete_Lindblad_binary},
these two elements of the evolution --- coherent evolution and dissipation ---
are applied in a staggered fashion, i.e., the map $\rho_n \mapsto \rho_{n + 1}$
is a composition of $\rho \mapsto \rho + \tau^2 \mathcal{D}[L_1] \rho$ and
$\rho \mapsto U_F \rho U_F^{\dagger}$. For comparison, the continuous-time form
of the Lindblad equation reads
\begin{equation}
  \frac{d \rho}{d t} = -i [H, \rho] + \gamma \mathcal{D}[L] \rho.  
\end{equation}
We note that while in the present context of noise in Floquet systems the jump
operator $L_1 = H_1$ is always Hermitian, this is not the case in general. Then,
the dissipator should be modified to
$\mathcal{D}[L] \rho = L \rho L^{\dagger} - \frac{1}{2} \left\{ L^{\dagger} L,
  \rho \right\}$.
In the usual (continuous-time) quantum master equation in Lindblad form, the
dissipator describes the effect of an environment or bath on the system
dynamics, and assumes the above time-local form if (i) the system-bath coupling
is weak and (ii) the bath correlation time is much shorter than the time scales
of the system dynamics (see, e.g.,~\cite{Gardiner2000}). Conditions (i) and (ii)
justify the Born and Markov approximations, respectively, which are made in
derivations of the Lindblad equation. In the derivation of the FLE, the
expansion in the noise strength $\tau$ is analogous to the Born approximation,
and the noise we consider is Markovian (i.e., uncorrelated on the intrinsic time
scale $T$ of the evolution of system) by assumption. Just like the usual master
equation, the FLE can immediately be seen to be trace-preserving and completely
positive.

\subsection{Multi-step piecewise constant driving}
\label{sec:multi-phase-pc}

The above derivation can be generalized straightforwardly to extended driving
protocols, defined in terms of a sequence of Hamiltonians $H_i$ with
$i = 1, 2, \dotsc, M$ which are applied for times $T_i$ so that the duration of
a full driving cycle is $T = \sum_{i = 1}^M T_i$. Assuming --- as in the main
text --- that there is timing noise in each step of the driving cycle, the FLE
that generalizes Eq.~\eqref{eq:discrete_Lindblad_binary} takes the form
\begin{equation}  
  \label{eq:discrete_Lindblad}
  \rho_{n + 1} = U_F \left( \rho_n + \tau^2 \sum_{i = 1}^M \mathcal{D}[L_i]
    \rho_n \right) U_F^{\dagger},
\end{equation}
where $U_F = U_M \dotsb U_1$ with $U_i = e^{-i T_i H_i}$. Here we take the
shifts $\tau_{ni}$ in different steps and driving cycles to be uncorrelated and
identically distributed,
$\overline{\tau_{ni} \tau_{n' i'}} = \tau^2 \delta_{nn'} \delta_{ii'}$. The jump
operators are given by
\begin{equation}
  \label{eq:L}
  L_i = U_F^{\dagger} U_M \dotsb U_i H_i U_{i - 1} \dotsb U_1 = U_1^{\dagger}
  \dotsb U_{i - 1}^{\dagger} H_i U_{i - 1} \dotsb U_1.
\end{equation}
As above, the jump operators are Hermitian, $L_i^{\dagger} = L_i$.

\subsection{Reduction to classical master equation}
\label{sec:reduct-class-master}

In many cases of practical interest, the discrete Lindblad equation for the
density matrix can be reduced to a classical master equation for the diagonal
elements of the density matrix written in the basis of Floquet eigenstates
$\ket{\alpha}$. The latter are the right eigenvectors of the Floquet operator,
$U_F \ket{\alpha} = e^{-i T \epsilon_{\alpha}} \ket{\alpha}$, and
$\epsilon_{\alpha}$ is the respective quasienergy. In this basis, the density
matrix can be written as
\begin{equation}  
  \rho = \sum_{\alpha, \beta} \rho^{\alpha \beta} \ket{\alpha}
  \bra{\beta}. 
\end{equation}
The condition for obtaining a classical master equation is that upon repeated
application of $U_F$ in Eq.~\eqref{eq:discrete_Lindblad}, the off-diagonal
elements of the density matrix dephase more rapidly than they are repopulated
from the diagonal. Then, for the diagonal elements
$\rho^{\alpha} = \rho^{\alpha \alpha}$ we obtain the discrete master equation
\begin{equation}
  \label{eq:master_equation}
  \rho^{\alpha}_{n + 1} =  \sum_{\beta} W_{\beta \to \alpha} \rho^{\beta}_n,
\end{equation}
where the transition probabilities are given by
\begin{equation}
  \label{eq:transition_probabilities}
  W_{\beta \to \alpha} = \delta_{\alpha \beta} + \tau^2 \sum_i \left( \abs{\braket{\alpha | L_i |
        \beta}}^2 - \braket{\alpha | L_i^2 | \alpha} \delta_{\alpha \beta}
  \right).
\end{equation}
The transition probabilities are symmetric,
$W_{\alpha \to \beta} = W_{\beta \to \alpha}$, and conservation of probability
--- in other words, conservation of the trace of the density matrix --- results
in $\sum_{\beta} W_{\beta \to \alpha} = 1$. Thus, the probability to remain in
state $\alpha$ can be written as $W_{\alpha \to \alpha} = 1 - \sum_{\beta \neq
  \alpha} W_{\beta \to \alpha}$, and Eq.~\eqref{eq:master_equation} can be
recast as
\begin{equation}  
  \rho^{\alpha}_{n + 1} = \rho^{\alpha}_n + \sum_{\beta \neq \alpha} W_{\beta
    \to \alpha} \left( \rho^{\beta}_n - \rho^{\alpha}_n \right),
\end{equation}
which is the form quoted in the main text.

\section{Decay of an End-State in a Noisy Floquet Topological Chain}

Here, we apply the formalism derived above to derive the time evolution of an
imperfectly driven ladder as introduced in the main text. We show that timing
noise leads to diffusion exactly on resonance and to exponential decay away from
resonance. The microscopic parameter entering the discussion is the hopping
parameter $J = 2\pi/T + \delta J$, with $\delta J$ being a measure for the
bulk's bandwidth. We assume that $\delta J$ is small, enabling us to make
analytical progress. In the following, we will distinguish the point of resonant
driving $\delta J=0$ from the limit $\delta J \to 0$ and show how they capture
the diffusive and exponential decay, respectively. Working in the eigenbasis of
the Floquet operator, we find that transitions between Floquet states are
induced by the noise at order $J^2\tau^2$. For small deviations from resonance,
in the calculation of the transition probabilities we can set $\delta J = 0$
(thus, ignoring terms $O(J^2\tau^2 \delta J^2)$). The crucial point is that we
have to be careful in taking limit $\delta J \to 0$: starting from finite
$\delta J$ and sending it to zero, the bulk states remain delocalized with an
arbitrarily small but finite quasienergy and bandwidth; on the other hand,
working \emph{exactly} on resonance, we should take the bulk states to be
localized with an exact degeneracy of all states at quasienergy $0$.

Employing the formalism of the Floquet-Lindblad equation derived above, we note
that the periodic drive considered in the main text comprises $M = 4$ individual
steps. Thus, the Floquet operator is $U_F = U_4 U_3 U_2 U_1$, with the
individual steps $U_i = e^{-iT H_i/4}$ where
$H_i = -\sum_{\mu \nu}J_{\mu \nu}^i \left( c^\dagger_{\mu} c_{\nu} + \mathrm{H.c.}
\right)$,
as introduced in the main text. The time evolution during one step can thus be
evaluated assuming single particle states and using the fact that the sum over
bonds in $H_i$ is over mutually disconnected pairs of neighboring lattice
sites. Using for $n \geq 1$ the identities
\begin{align}
  H_i^{2n} = J^{2n} \sum_{\braket{\mu \nu}_i} \left( n_\mu + n_\nu \right)
  \qquad \text{and} \qquad
  H_i^{2 n - 1} = J^{2 \left( n - 1 \right)} H_i,
\end{align}
with the sum running over all bonds affected by the hopping in step $i$, we find
\begin{align}
	U_i = 1 + \left( \cos(\phi) - 1 \right) \sum_{\braket{\mu \nu}_i} (n_\mu
  + n_\nu) - i \sin(\phi) \, H_i/J,
\end{align}
with $\phi = JT/4$. Equivalently, we can represent the $U_i$ in the single
particle basis with the projector on all involved lattice sites
$P_i = \sum_{\braket{\mu \nu}_i} \left(\ket{\mu}\bra{\mu} +
  \ket{\nu}\bra{\nu}\right)$ and its complement $Q_i = 1 - P_i$:
\begin{align}
	U_i = P_i \left(\cos(\phi) - i \sin(\phi) \, H_i /J\right) P_i + Q_i ,
\end{align}
as used in the main text. In steps $1$ and $3$ all lattice sites are involved
and $Q_1=Q_3=0$, while in steps $2$ and $4$ one site at each end drops out of
the respective Hamiltonian and thus
$Q_2 = \ket{1, + }\bra{1, + }+\ket{L-1,- }\bra{L-1,-}$ and
$Q_4 = \ket{0,- }\bra{0,- }+\ket{L,+ }\bra{L,+}$.  The jump operators $L_i$
introduced above are $L_1 =H_1$, $L_2 = U_1^{\dagger}H_2 U_1$,
$L_3 =U_1^{\dagger} U_2^{\dagger}H_3 U_2 U_1$, and
$L_4 =U_1^{\dagger} U_2^{\dagger} U_3^{\dagger}H_4 U_3 U_2 U_1$.

The form of the jump operators simplifies further exactly on resonance, when
$\phi = \pi/2$ and during each step of the driving protocol a particle is fully
transferred from one lattice site to another. Then, we find $L_1 =H_1$,
$L_2=H_1H_2H_1$,
$L_3=H_1 \left( H_2-iQ_2 \right) H_3 \left( H_2+iQ_2 \right)H_1$, and
$L_4=H_1 \left( H_2-iQ_2 \right) H_3H_4H_3 \left( H_2+iQ_2 \right)H_1$.
Moreover, on resonance the Floquet eigenstates take a particularly simple form:
The end states are localized to the outermost lattice sites
$\ket{e_l} = \ket{0,-}$ and $\ket{e_r}=\ket{L,+}$, while for the bulk states we
can choose the lattice site basis $ \ket{b} = \ket{j,s}$ with $j=1,\dots,L-1$
and $s=\pm$. In the lattice basis, the jump operators are represented as
$L_1 = J \sum_{j} \left( \ket{j,+}\bra{j-1,-} + \mathrm{H.c.} \right)$,
$L_2 = J \sum_{j} \left( \ket{j,+}\bra{j,-} + \mathrm{H.c.} \right)$,
$L_3 = J \sum_{j} \left( \ket{j,+}\bra{j+1,-} + \mathrm{H.c.} \right)$, and
$L_4 = J \sum_{j} \left( \ket{j,+}\bra{j,-} + \mathrm{H.c.} \right)$. Focusing on the left
end state $\ket{e} \equiv \ket{e_l}$, we immediately see that the diagonal
matrix elements of the jump operators vanish, $\braket{e|L_i|e} = 0$, and that
the left end state is an eigenstate of the square of the jump operators:
\begin{equation}  
  \label{eq:left_edge_L2}
  L_i^2 \ket{e} =
  \begin{cases}
    J^2 \ket{e}, & i=1,3, \\
    0, & i = 2,4.
  \end{cases}
\end{equation}
This, in turn, implies that $\braket{e|L_i^2|e}=J^2$ for $i=1,3$ and
$\braket{e|L_i^2|e}=0$ otherwise, while $\braket{b | L_i^2 | e} = 0$. Below, we
investigate the stability of the chain's end state very close to as well as
exactly on resonance. In both cases, we can use the above expressions for matrix
elements of the jump operators which we derived for resonant driving. The
crucial difference comes from the matrix elements $\braket{b | L_i | e}$, which
we discuss in detail below.

In the following, we restrict our attention to the left end of the chain which
we denote by $\ket{e}$, and we disregard the other one, assuming a semi-infinite
system. Then, the system's density matrix can be represented in the Floquet
eigenbasis as
\begin{align}
  \rho_n = \rho_n^e \ket{e}\bra{e} 
  + \sum_b \left(\rho_n^{eb} \ket{e}\bra{b} + \rho_n^{be} \ket{b}\bra{e}  \right) 
  +\sum_{bb'}\rho_n^{bb'}\ket{b}\bra{b'} .
\end{align}
Evolving the density matrix by one driving period with the
FLE~\eqref{eq:discrete_Lindblad}, we find for the occupation of the chain's left
end state:
\begin{multline}
  \rho^e_{n+1} = \left(1 + \tau^2 \sum_i \left( \left| \braket{e|L_i|e}
      \right|^2 - \braket{e|L_i^2|e}\right) \right) \rho^{e}_n \\ + \tau^2
  \sum_i \sum_b \left[ \left( \braket{e | L_i | e} \braket{b | L_i | e} -
      \frac{1}{2} \braket{b | L_i^2 | e} \right) \rho_n^{eb} + \left( \braket{e
        | L_i | b} \braket{e | L_i | e} - \frac{1}{2} \braket{e | L_i^2 | b}
    \right) \rho_n^{be} \right] \\ + \tau^2 \sum_i \sum_{bb'} \braket{e|L_i|b}
  \braket{b'|L_i|e} \rho^{bb'}_n.
\end{multline}
This can be simplified using the matrix elements evaluated above, and we obtain
(cf.\ Eq.~(5) of the main text)
\begin{align}
  \label{eq:evo_rho_exp}
  \rho^{e}_{n+1} 
  = \left(1- 2 \tau^2 J^2 \right) \rho^{e}_n
  + \tau^2 \sum_i \sum_{bb'} \braket{e|L_i|b} \braket{b'|L_i|e} \rho^{bb'}_n,
\end{align}
where, in particular, coherences between the end state and bulk states drop out.
    
\subsection{Exponential Decay for a Dispersive Bulk} 
\label{sub:exponential_decay}

We first analyze the noisy stroboscopic time evolution of the end states in
Eq.~\eqref{eq:evo_rho_exp} for a dispersive bulk, i.e., $\delta \phi \neq 0$,
for which the bulk states acquire a finite band width $\propto \delta \phi$ and
are delocalized over the full length of the system. An analytical treatment is
possible only close to resonant driving, and hence we focus on the limit
$\delta \phi \to 0$. This allows us to approximate the end state $\ket{e}$ and
the jump operators with their forms at resonant driving. Importantly, the bulk
states have finite albeit arbitrarily small quasienergies.

For the matrix elements connecting the end to the bulk states, $\braket{b|L_i|e}$,
we need to keep in mind that the jump operators transfer electrons only between
lattice sites within a short range of each other. Assuming that the bulk states
are evenly spread out over the entire chain, they carry a normalization factor
of $\sim 1/\sqrt{L}$, leading to the following scaling of the matrix elements
with system size:
\begin{align}
  \braket{b|L_i|e}\sim {1\over \sqrt{L}} .
\end{align}

Moreover, it is important to note that in the delocalized basis the off-diagonal
part of the bulk density matrix $\rho_n^{bb'}$ for $b\neq b'$ is non-zero. These
matrix elements are getting populated by excitations from the end state at a
rate $\tau^2 \braket{b|L_i|e}\braket{e|L_i|b'}\sim \tau^2/L$ and, to leading
order, evolve coherently with $U_F$ per period. The off-diagonal elements
therefore pick up phase factors of $e^{i T \left( \varepsilon_b-\varepsilon_b' \right)}$
for each full period. No matter how small the quasienergies, these phases
accumulate and average to zero for long enough times. Only the diagonal elements
are static and contribute to the sum over bulk states in
Eq.~\eqref{eq:evo_rho_exp}, which then decays with the system size as can be
seen from a simple dimensional analysis
\begin{align}
	\sum_{b} \rho^{bb}_n \left|\braket{e|L_i|b}\right|^2 \propto {1\over L} .
\end{align}
For large enough systems and long enough times, the contribution from this sum
can be neglected and leaves us with an exponential decay of the end state (cf.\
Eq.~(6) of the main text)
\begin{align}
  \rho^{e}_{n+1} 
  = \left(1- \tau^2 J^2 \right) \rho^{e}_n + O(1/L),
\end{align}
confirming the intuitive picture that a dispersive bulk would carry away any
excitation out of the edge. This contrasts the case of the localized bulk at
resonant driving in which excitations get stuck close to the end state and have
a finite return probability, as we discuss in the next section.


\subsection{Diffusive Decay at Resonant Driving} 
\label{sub:diffusive_decay_at_resonant_driving}

We now focus on the dynamical evolution of the density matrix for a flat bulk
band at resonant driving. As explained above, in this case we can work in the
basis of lattice sites, in which the end states are
$\ket{e_l} = \ket{0,-} \equiv \ket{e}$ and $\ket{e_r}=\ket{L,+}$, while the bulk
states are give by $ \ket{b} = \ket{j,s}$ with $j=1,\dots,L-1$ and $s=\pm$. Note
that the lattice site basis for the bulk is an arbitrary choice since the bulk
band is fully degenerate.

In the equation for the end state population~\eqref{eq:evo_rho_exp}, only $L_1$
and $L_3$ connect the end state to the bulk states, such that
$\braket{e|L_i|b}=J$ for $i=1,3$ and $\ket{b_1}=\ket{1,+}$ and
$\ket{b_3}=\ket{1,-}$. All other elements $\braket{e|L_i|b}$ are zero. Then, the
dynamical equation reduces to
\begin{align}
  \label{eq:edge_diffusive}
	\rho^{0,-}_{n+1} =
	\left(1 - 2 \tau^2 J^2\right) \rho^{0,-}_{n} 
	+ \tau^2 J^2 \left( \rho^{1,+}_{n} + \rho^{0,-}_{n}\right) ,
\end{align}
where we introduced a notation for the density matrix diagonal labeled by
lattice sites $\rho_n^{bb}=\rho_n^{j,s}$. We note that all off-diagonal elements
of the density matrix drop out of Eq.~\eqref{eq:edge_diffusive} --- a property
that is also true for the evolution equation of the population of a generic bulk
state. Thus, starting from an initially diagonal density matrix,
$\rho_n = \ket{e} \bra{e}$, the density matrix remains diagonal at all
times. Importantly, the time evolution for the end state is strongly coupled to
its neighboring bulk states. To understand its behavior in the long-time limit
we can alternatively look at the dynamical equation for a generic bulk state,
\begin{equation}
  \begin{split}
  \rho^{j,s}_{n+1} &= \left(1 - \tau^2 \sum_i \sum_{l,s'}
                     \braket{l,s'|L_i|j,s}^2 \right) \rho^{j,s}_n + \tau^2 \sum_i \sum_{l,s'}
                     \braket{l,s'|L_i|j,s}^2 \rho^{l,s'}_n
                      \\
                   &=  \left(1 - 4\tau^2 J^3\right) \rho^{j,s}_n
                     + \tau^2 J^2 \left( 2 \rho_n^{j,-s} + \rho_n^{j+s,-s} +
                       \rho_n^{j-s,-s} \right).
\end{split}
\end{equation}
The last equation is Eq.~(7), and can straightforwardly be reduced to the
diffusion equation~(8) for the total occupation of doublets
$\rho^j_n = \rho^{j, +}_n + \rho^{j,-}_n$. As bulk and end state can be treated
in the same way at resonant driving, being localized to one lattice site each,
the behavior of the end state for long times can be inferred immediately. It
obeys a diffusion equation against a hard wall boundary for a particle
initialized right next to this boundary.

Note that also in the resonant driving case we can choose to describe the bulk
in a delocalized basis, rendering the off-diagonal elements $\rho_n^{bb'}$
non-zero. However, in contrast to the dispersive case, these elements are static
in the stroboscopic time-evolution due to the flat bulk band and can thus not be
neglected.


\section{Numerical simulations and disorder effects}

Adding a chemical potential disorder to the Floquet system breaks particle-hole symmetry, such that the end modes are no longer topologically protected. However, since in the clean case the end modes are well localized at the ladder boundaries and separated in quasienergy from bulk states, one can expect these features to persist for small enough disorder strength. We show here that this is indeed the case, by determining the average localization lengths and positions of Floquet eigenstates in the presence of on-site disorder. On a disordered ladder consisting of $L=200$ rungs, we compute the real space probability distribution of each Floquet eigenstate, $|\psi(x)|^2$ with $x=1,\ldots, L$ indexing the rungs of the ladder, and find the associated localization length by fitting with an exponential decay of the form
\begin{equation}\label{eq:loclen}
 f(x) = A \exp \left( -\frac{|x-x_0|}{\xi} \right).
\end{equation}
The localization length $\xi$ is extracted from the fit, while the amplitude $A$ and the wavefunction center $x_0$ are taken to be the value and position of the maximum of each $|\psi(x)|^2$. From the values of $x_0$ we also determine the average displacement of states from the center of the chain, $d=|x_0 - L/2|$. Since $x_0 \in [1,L]$, the displacement takes values $d\in [0, L/2]$. States located on the boundaries of the system are therefore expected to show a large displacement $d\to L/2$, while states that are uniformly distributed throughout the bulk should have an average displacement $d\to L/4$.

The results shown in Fig.~\ref{fig:loclen} confirm our expectations. States with quasi energies close to $\varepsilon=0$ are uniformly distributed throughout the chain. They become more localized and spread out in quasienergy with increasing disorder strength.
However, states close to $\varepsilon=\pm\omega/2$ are localized exclusively on the chain boundaries and have a much smaller localization length as compared to bulk states. As such, even with a small on-site disorder strength, edge modes remain well localized on the system boundaries and separated in quasi-energy from the bulk states.

\begin{figure}[tb]
 \includegraphics[width=0.9\textwidth]{statestat.pdf}
 \caption{Average displacement $d$ (left) and localization length $\xi$ (right) of Floquet eigenstates on ladder of $L=200$ rungs, using $\phi=1.45$ as in the main text Fig. 2. The on-site disorder strength is set to $VT/4=0.2$ (solid orange squares) and $VT/4=0.5$ (empty blue squares). Each point represents an average over states in a range of quasienergies centered at $\varepsilon$ and having a width of $\omega/101$. A total of $10^4$ independent disorder realizations were used for each disorder strength. As disorder strength is increased, bulk states around $\varepsilon=0$ spread in quasienergy and become more localized, but remain uniformly distributed throughout the system, showing an average displacement $d\simeq L/4=50$. Edge modes close to $\varepsilon=\pm\omega/2$ show a much smaller localization length as compared to bulk states, and remain localized at the ends of the system, with an average displacement $d\simeq L/2=100$. \label{fig:loclen}}
\end{figure}


As a final point, we show numerically that hopping
disorder also slows the end mode decay down to a diffusive process. This is
indicated in Fig.~\ref{fig:1dhopping}, and is a consequence of the fact that
hopping disorder in this model leads to a localization of all bulk states. While
localization is expected for on-site disorder, which breaks the particle-hole
symmetry, this is not immediately obvious in the case of random hoppings, since
the particle-hole symmetry of the Floquet operator, Eq.~\eqref{eq:FPHS}, is
preserved for every disorder realization. One-dimensional particle-hole
symmetric chains may enter a so-called critical phase, characterized by the
presence of delocalized states \cite{Brouwer2000, Motrunich2001, Brouwer2003,
  Gruzberg2005}. In the 1D ladder model, however, Floquet bulk bands form a
parabolic (as opposed to linear) band touching point at the particle-hole
symmetric quasienergy $\varepsilon=0$, meaning that the bulk state velocity
vanishes already in the clean limit. In this way, the 1D critical phase is
avoided, all bulk states become localized, and the decay is diffusive as in the
case of on-site disorder.

\begin{figure}[tb]
 \includegraphics[width=0.4\textwidth]{1d_disorder_hopping.pdf}
 \caption{As in the main text, we plot the survival probability of an initial
   end mode, $s=\overline{\braket{e | \psi_n}^2}$ as a function of the
   number of noisy driving cycles, $n$. All system parameters, as well as the
   blue and yellow curves are the same as in the main text Fig.~2. The orange
   curve however is obtained by introducing disorder in the hopping amplitudes,
   using $V_{\rm hop}T/4=0.75$, and averaging over 4000 independent realizations
   of disorder and noise. Since hopping disorder also localizes all bulk states,
   the decay is again diffusive, as indicated by the dashed
   line.\label{fig:1dhopping}}
\end{figure}

