\section{introduction}
Kondo-lattice systems are prototypical strongly correlated
materials, in which thermal and quantum fluctuations can drive the electronic
states in spin, charge, or orbital channels and induce
various phases. In the simplest model, the ground states of  heavy-fermion
compounds are determined by the competition between the Ruderman-Kittel-Kasuya-Yosida (RKKY)
interaction and the Kondo effect, which evolve differently with
temperature, magnetic field, and pressure \cite{Stewart,Hilbert,Gegenwart,Weng2016}. Adjusting  non-thermal control parameters can tune the strengths of both the RKKY coupling, which is the intersite exchange interaction between localized moments leading to magnetic order, and the Kondo coupling corresponding to onsite antiferromagnetic exchange  between localized moments and itinerant conduction electrons, giving rise to a nonmagnetic singlet ground state. To date, many novel phenomena, such as non-Fermi-liquid behavior and unconventional
superconductivity, have been found in the vicinity of quantum
critical points (QCPs), where the Kondo coupling prevails and the magnetic ordering temperature is smoothly suppressed to zero. However, whether a universal description of quantum critical behavior can be found still needs to be determined, and in particular it is necessary to investigate  how the emergent phenomena evolve upon tuning Kondo lattice systems from localized to itinerant states.

The Kondo lattice compound CeRhIn$_5$  provides a good opportunity to
study the change from localized to delocalized $4f$ states. In zero-field at ambient pressure, the Ce 4$f$
electrons in CeRhIn$_5$ are localized and order antiferromagnetically at around
3.8~K \cite{Hegger,Shishido2002,Harrison2004}. In high magnetic fields, evidence for a field-induced
spin-density-wave (SDW) type QCP has been found near the critical field
$B_{c0}\approx50$~T \cite{Jiao}.
Theoretically, this kind of conventional SDW-type QCP involves only three-dimensional fluctuations of the
antiferromagnetic (AFM) order parameter around zero temperature
\cite{Hertz,Moriya,Millis}, which has been identified in some heavy-fermion compounds,
such as CeCu$_2$Si$_2$ \cite{Arndt}, and La-doped CeRu$_2$Si$_2$ \cite{Kambe}.
In this scenario, the $f$ electrons already are delocalized in the magnetically ordered
region, and only a few ``hot spots/lines" on the Fermi surface (FS) become critical
at the QCP \cite{Gegenwart}. For CeRhIn$_5$, a field-induced change of the FS was
observed inside the AFM state at around $B^{\ast}$$\approx$~30~T, as concluded
from sharp changes in the dHvA frequencies, Hall effect, and magnetoresistance
\cite{Jiao,Moll,Jiao2}. The expansion of the FS volume (or charge carrier concentration) is attributed
to the $4f$ electrons becoming delocalized at $B>B^*$ \cite{Jiao}, which is consistent with the presence of a  conventional SDW-type QCP at  $B_{c0}$. Such a field-induced delocalization of the $4f$ electrons is different to that generally anticipated for Kondo lattice systems, where the $4f$ electrons are often expected to be more localized in high fields \cite{Aoki1993,Harrison1993,Harrison1998,Ce142}.

To date, no clear anomalies have been detected around $B^*$ in the specific heat \cite{Jiao}, torque
magnetometery \cite{Jiao,Moll}, and magnetic susceptibility \cite{Takeuchi}. Therefore, the evolution of the
Ce 4$f$ electrons in high magnetic fields is still not determined, and the nature of
the quantum phase transition at $B^*$ is not well understood \cite{Jiao}. Moreover, evidence was recently provided
for a spin nematic state in CeRhIn$_5$ at high fields, from resistivity measurements of microstructured devices for certain field orientations \cite{Ronning2017}. A pronounced in-plane
anisotropy was found in the resistivity when a field was applied 20$^\circ$ away from the $c$~axis, where
the anisotropy showed a sudden increase above $B^*$. Obviously, the relationship between this enigmatic
phase transition and the local to itinerant transition of the $4f$ electrons needs to be resolved.
Here, we report a detailed study of the magnetostriction, specific heat, and de Haas-van
Alphen (dHvA) effect of CeRhIn$_5$ in high magnetic fields. These results enable us
to trace the evolution of the correlated electrons upon approaching $B^*$.

\section{Experimental details}

Single crystals of CeRhIn$_5$ were grown using In flux, as described elsewhere \cite{Jiao}.
Measurements of the absolute values of the specific heat are extremely challenging in pulsed
magnetic fields. These were performed in a long-pulsed magnet at the
International Megagauss Science Laboratory of ISSP in Kashiwa up to 43.5~T. To
minimize the influence associated with the field instability of the pulsed field,
we generated highly stabilized ($\pm$100 Oe) magnetic fields with a 100 ms timescale
using a field-feedback controller \cite{Yoshi1} and measured the absolute value of
the specific heat in the stabilized field by applying the heat-pulse method \cite{Yoshi2}.
The resolution of the specific heat data were successfully improved \cite{sup}, and a
detailed analysis of the data is possible up to 43.5~T. The optical fibre Bragg
grating (FBG) method \cite{Daou,Rotter} was used to measure the magnetostriction in
pulsed magnetic fields up to 60 T with $B$~$\parallel$~$a$ and $B$~$\parallel$~$c$.
Measurements of the dHvA effect of CeRhIn$_5$ were conducted utilizing a
torque technique between 5~K and 300~mK in dc fields up to 45~T \cite{Jiao}.

\section{Results}
\subsection{High-field magnetostriction}

The magnetostriction of a metal is proportional to the
conduction-electron density of states at the Fermi level \cite{Ekreem}. Therefore,  so as to gain insights into the interactions between the
conduction electrons and the lattice, we performed measurements of the  longitudinal magnetostriction of
CeRhIn$_5$ for $B$~$\parallel$~$a$ and $B$~$\parallel$~$c$ at 0.7~K,  where the relative length changes $\Delta L/L$ and the corresponding field derivatives as a function of applied field are displayed in Fig.~\ref{Fig1}. It can be seen that at lower fields up to around 15~T, the in-plane $\Delta L/L$  is nearly field independent, before beginning to decrease more rapidly at higher fields. Meanwhile at very high fields in the vicinity of the critical field $B_c$, $\Delta L/L$ shows a change in slope at 50 and 52~T for fields along the $a$- and $c$-axes respectively. While no clear anomaly is observed near $B^*\approx30$~T along the $c$~axis, there is a kink in $\partial(\Delta L/L)/\partial B$ along the $a$-axis [Fig.\ref{Fig1}(b)].  These results give further evidence for the existence of $B^*$ from in-plane measurements, but also suggest that the
coupling to the lattice of the Fermi surface transition $B^*$ is rather weak. While the critical fluctuations associated with  phenomena such as CDW transitions and valence
changes generally couple strongly to the lattice, those corresponding to SDW order commonly couple more weakly \cite{Stockert2012}. We note that weak anomalies have also been reported in the magnetostriction at $B^*$ when fields are applied at 11$^{\circ}$ and 20$^{\circ}$ to the $c$-axis \cite{Ronning2017,Rosa2018}.  It is possible that the magnetostriction change at $B^*$ is too weak to be resolved in our measurements with fields along the $c$-axis.

\begin{figure}[tb]\centering
 \includegraphics[width=0.85\columnwidth]{Fig1Arx.eps}
\caption{Magnetostriction  of CeRhIn$_5$ in high magnetic fields, where (a) the fractional length change $\Delta$$L$/$L$, and (b) the derivative are displayed for $B$~$\parallel$~$a$ and $B$~$\parallel$~$c$ at 0.7~K, and  the former data have been scaled. The arrows mark the onset of the mass enhancement $B_{\rm onset}$, $B^*$ and the critical field $B_c$.}
\label{Fig1}
\end{figure}

\subsection{Specific heat  in high magnetic fields}

Figure~\ref{Fig2}(a) shows the temperature dependence of the specific heat of
CeRhIn$_5$ [$C(T)$] at different magnetic fields with $B$~$\parallel$~$c$, reaching a higher maximum field, and lower temperatures in the high-field region than previous studies of $C(T)$ \cite{Kim}; here the level of precision allows for a detailed study of the evolution of electronic correlations.
\begin{figure}[tb]\centering
 \includegraphics[width=0.99\columnwidth]{Fig2.eps}
\caption{(a) Temperature dependence of the specific heat $C$/$T$ of CeRhIn$_5$ in
various magnetic fields with $B$~$\parallel$~$c$. Data between 0--15 T were measured in a PPMS while
the high-field data was measured in pulsed magnetic fields. The orange curve at low $C/T$ is the
specific heat of LaRhIn$_5$ at $B$ = 0. (b) $C$/$T$ as a function of $T/T_N$ in zero field and 41.3~T, one
 field above and one below $B^*$. $C$/$T$ is displayed as a function of $T/T_N$  for applied fields (c) up to 16.6~T, and (d) at fields higher than 16.6~T.}
\label{Fig2}
\end{figure}
The AFM transition is clearly observed at all fields, which is
manifested by a pronounced peak in $C(T)/T$. With increasing magnetic
field, $T_N$ smoothly shifts to lower temperatures, consistent with previous
reports \cite{Jiao,Kim}. The accuracy of the current data obtained in pulsed
magnetic fields was checked by comparing to those measured in dc fields \cite{Kim,sup}. $C/T$ in the
displayed temperature range is much larger than in the La homologue, and therefore the phonon contribution is negligible.

Figure~\ref{Fig2}(b) displays $C/T$ as a function of $T/T_N$ for two representative curves [bold in Fig.~\ref{Fig2}(a)], at $B=0$ and 41.3~T, which show distinctly different behavior. At high fields, the transition is sharper, with a larger jump, while the values of $C/T$ just above $T_N$ are reduced. This suggests that compared to the low-field results, at high fields there is a reduction
in the contribution from short-range fluctuations which appear above the magnetic ordering temperature. In Fig.~\ref{Fig2}(c), $C/T$ is displayed for fields up
to 16.6~T as a function of $T/T_N$  \cite{sup}. One can see that the data below $T_N$ overlap fairly well, demonstrating that $T_N$
is the only energy scale which determines the specific heat in this field range. On the other hand, the higher field data do not
scale so well as a function of $T/T_N$, suggesting the role of additional parameters. This may also be taken as an indication of changes in the underlying correlated state, such as a variation of the effective mass, as described below.

We further analyzed the data well below $T_N$ using $C$/$T$ = $\gamma$ + $C_m$/$T$, where $\gamma$ is the Sommerfeld coefficient, and $C_m$ is
the contribution of the AFM magnons to the specific heat \cite{Continentino,sup}. This equation can be fitted well to
the low-temperature data [$T$ $<$ 0.7$T_N$]  below 35.3~T. Above 35.3~T,
the fit is less reliable due to the limited accessible data  for $T$ $<$ 0.7$T_N$. A
detailed discussion of the specific heat fits is presented in the Supplemental Material \cite{sup}, where with increasing magnetic field, $\gamma$($B$) changes little at lower fields,
but shows an enhancement onsetting at fields above  $B_{\rm onset}\approx17$~T.

\begin{figure}[tb]\centering
 \includegraphics[width=6cm]{Fig3.eps}
\caption{(a) and (b) FFT spectra of dHvA effect oscillations of CeRhIn$_5$ for the field range
35--45 T at various temperatures with for $B \parallel c$. (c) Temperature dependence of dHvA oscillation
amplitudes of the $\alpha_3$ orbit at various magnetic fields. The field values correspond to the midpoints of the 10~T wide field windows used in the FFT. The dashed lines
display fits to the spin-dependent LK formula. The inset shows the field dependence
of the cyclotron masses of the two spin orientations (open blue and green squares).
Experimental data obtained at low magnetic fields from
Refs.~\onlinecite{Cornelius00,Hall} are also shown by the red solid squares.}
\label{Fig3}
\end{figure}

\subsection{Analysis of the effective cyclotron masses from the dHvA effect}

To directly probe the effective mass of the charge carriers, we also measured the dHvA effect.
It has been established by band-structure calculations and quantum-oscillation
measurements that the FS of CeRhIn$_5$ primarily arises from three bands ($\alpha$,$\beta$,$\gamma$), where the dHvA effect is dominated by the extremal orbits denoted $\alpha_1$, $\alpha_2$, $\alpha_3$, $\beta_1$, and
$\beta_2$ \cite{Hall,Shishido2002}. In our previous work
(Refs.~\onlinecite{Jiao,Jiao2}), we studied the magnetic field dependence of the dHvA
frequencies at 0.33 K, which yielded clear evidence for an abrupt expansion of the FS volume at $B^*$. Here, we further investigate the temperature and magnetic-field
dependence of the dHvA amplitudes, from which the effective cyclotron masses of the charge carriers are determined, focussing on the $\alpha_3$ orbit, which corresponds to the strongest oscillation.

In a conventional metal, the dHvA oscillation amplitudes are given by the
Lifshitz–Kosevich (LK) formula \cite{Shoenberg}, where in the nonspin-dependent case the temperature and field dependence are determined by two adjustable parameters, the effective cyclotron mass ($m_c$) and the Dingle temperature ($T_D$),  the analysis of the latter being presented in the Supplemental Material \cite{sup}. Figures~\ref{Fig3}(a) and (b) display  the fast Fourier transform (FFT) spectra of the dHvA oscillations of CeRhIn$_5$ at various
temperatures for 35~T~$\leq B \leq 45$~T with $B\parallel c$.
\begin{figure}[tb]\centering
 \includegraphics[width=7cm]{Fig4.eps}
\caption{Temperature-field phase diagram of CeRhIn$_5$ where the field dependence of $\gamma$ and the effective mass, $m_c$, of the heavy spin orientation of the FS sheet for $B \parallel c$ [see inset of Fig.~\ref{Fig3}(c)] are  shown by the olive circles and red squares respectively. The scales for these quantities with corresponding colors are on the right side of the panel.  The positions of the FS reconstruction ($B^*$) and the field at which the effective mass enhancement onsets ($B_{\rm onset}$) are also displayed. The triangular and diamond shaped symbols represent $T_N$ as obtained here and in Ref. \onlinecite{Jiao}, respectively.}
\label{Fig4}
\end{figure}
If the data were well described by the LK  formula, the amplitude of the oscillations
would monotonically increase with decreasing temperature. As shown in Fig.~\ref{Fig3},
this is the case for the $\alpha_3$ orbit above 1~K only, but below this temperature
the dHvA amplitude reaches a maximum before it decreases with decreasing temperature. This decrease cannot arise due to an increased scattering rate, since it is found that $T_D$ decreases strongly in this temperature and magnetic field range \cite{sup}. Such a phenomenon was previously observed in a
pulsed-field experiment in the range 31--50 T and was attributed to
the formation of SDW order \cite{Cornelius00}. However, this explanation cannot account for the similar  observations of the $\alpha_3$ orbit of CeCoIn$_5$, where it was ascribed to spin-dependent
mass enhancements of the FS \cite{McCollam}. This phenomenon may arise in heavy fermion systems when a sufficiently high magnetic field lifts the spin degeneracy of the renormalized bands by the Zeeman effect. Due to the asymmetry of the electronic bands near the Fermi level in the heavy fermion state, this can lead to  different effective masses for the two spin species \cite{McCollam}. Therefore, we studied the
spin-dependent effective masses of the electrons on the $\alpha_3$ orbit by performing
FFT on the dHvA oscillations at various temperatures in 10~T
intervals.

Following Refs.~\onlinecite{Takashita,McCollam}, a spin-dependent LK formula was
used to fit the derived amplitudes given by $\tilde{M}~=~\sqrt{(\tilde{M}_{\uparrow}+\tilde{M}_{\downarrow})^2\textrm{cos}^2
\theta+(\tilde{M}_{\uparrow}-\tilde{M}_{\downarrow})^2\textrm{sin}^2\theta}$, where $\theta$ is the phase term and $\tilde{M}_{\uparrow}$ ($\tilde{M}_{\downarrow}$)
is the dHvA amplitude of the spin up (down) FS described by the LK formula. This model fits the experimental results well [dashed-dotted
lines in Fig.~\ref{Fig3}(c)]. As shown in the inset of Fig.~\ref{Fig3}(c), the fit reveals a relatively large field-dependent effective mass
for one of the spin orientations but a small, nearly field-independent mass for the
other one. It should be noted that we cannot determine which spin orientation corresponds to the carriers with the larger effective
mass. The mass of the lighter band is around 1.6$m_e$ in the measured field
range, which is slightly larger than for the corresponding orbit in
LaRhIn$_5$~\cite{Shishido2002}, while the mass of the heavier band increases
with increasing magnetic field.  Note that this mass reaches a value of about 20$m_e$ when the midpoint field is 35~T, at which the dHvA amplitude displays a significantly steeper temperature dependence. Such a heavy mass at a field above $B^*$ is similar to that of the corresponding orbit of CeCoIn$_5$ \cite{McCollam}, where the $4f$ electrons also contribute to the Fermi surface \cite{Harrison2004}. A similar suppression of the dHvA amplitude at low temperatures is
observed for the $\beta_2$ orbit \cite{sup}. Since dHvA signals of the $\alpha_2$, $\alpha_1$, and $\beta_1$
orbits only appear below 2.5~K and above 30~T, the field dependences of these cyclotron masses cannot be estimated.

\section{Discussion}

Our main results are  summarized in Fig.~\ref{Fig4}, which show that although $T_N$ is smoothly suppressed as a function of applied field, there is an enhancement of the $\gamma$ coefficient and effective
cyclotron mass $m_c$ which onsets above  $B_{\rm onset}\approx17$~T, well below $B^*\approx30$~T. Combined with our previous results \cite{Jiao,Moll,Jiao2}, we are able to identify a number of distinct behaviors across different field ranges at low temperature. Between zero field and $B^*$, the $4f$ electrons are well localized and the Fermi surface is small. At lower fields up to $B_{\rm onset}$,  $\gamma$ barely changes with increasing field, and there is also little change in the in-plane magnetostriction. However, upon further increasing the field, the specific heat and dHvA measurements are consistent in revealing a clear enhancement of the effective quasiparticle mass, which likely signals an increasing fraction of the large FS in the fluctuation spectrum of the FS crossover \cite{JMMMXX}. We note that an anomaly at $B_M\approx18$~T was  previously found in Hall-resistivity measurements \cite{Jiao}, but whether this is related to the onset of the electronic mass enhancement requires further determination. At the subsequent crossover at $B^*$ $\approx$ 30~T, the Fermi surface becomes reconstructed from small to large.
Here, a maximum of the effective mass at $\approx$ 30~T may be anticipated, corresponding to the finite-temperature signature of a zero-temperature divergence  \cite{Gegenwart}. However, it is difficult to confirm whether such a maximum occurs, due to  the extremely large error bar of the cyclotron mass for the field-window centered around 35~T, in comparison to the preceding and following data points at 30~T and 40~T, respectively. It can be seen  in Fig.~\ref{Fig3}(c) that the 35~T data shows the steepest increase with decreasing temperature and is less well described by the  LK-formula for spin-dependent cyclotron masses, consistent with a rapid field-induced change of the electronic structure in this field-window. On the other hand, the absence of such a maximum  might be related to there being  a sizeable mass enhancement only for \textit{one} spin orientation. The high-field mass enhancement  significantly affects $C(T)$, where the analysis of the low-temperature data suggests a pronounced weakening of the RKKY coupling strength above $B_{\rm onset}$ \cite{sup}.  In addition, the short-ranged correlations above $T_N$ appear to be substantially reduced, although the magnetic structure likely remains unchanged to fields above $B^*$ \cite{UrbanoPoster}.

The field-induced increase of  the effective mass, together with  the observed changes to $C(T)$ in the magnetic state and the drastic changes of some dHvA frequencies upon crossing  $B^*$ with others remaining at a nearly constant value \cite{Jiao},  implies that the Fermi surface reconstruction is driven by the delocalization of the $4f$ electrons. Whether this delocalization transition occurs as an orbitally  selective Mott transition or other  forms of Kondo destruction is an interesting open question. Addressing it may require relating the observed changes of certain dHvA frequencies to details of the Fermi surface, and thus disentangling the effects of the $f$-electron degrees of freedom, momentum dependence of the hybridization and band structure of the conduction electrons. This may shed new light on how orbital and spin degrees of freedom are entangled in strong magnetic fields. Recent results indicating the emergence of an anisotropy of the electrical resistivity due to spin nematicity \cite{Ronning2017} at  $B$ $\approx$ $B^*$ in CeRhIn$_5$ also show a close relationship to the SDW phase with itinerant $4f$ electrons. Indeed, in the low-field regime no resistivity anisotropy is observed but at intermediate fields \textit{below} 30~T,  there is an onset of this resistivity anisotropy  near $B_{\rm onset}$ \cite{sup}.  This in turn suggests that the presence of itinerant, rather than localized, $4f$-electron states is related to the emergence of this new electronic (nematic) order.

\section{summary}

To conclude, we have probed the nature of the Fermi surface reconstruction associated with the  abrupt change of the Fermi surface volume of CeRhIn$_5$ at $B^*\approx30$ T. We find a field-induced enhancement of the effective mass which onsets  at  $B_{\rm onset}\approx17$~T and significantly increases with increasing field below $B^*$, as well as distinct changes to the temperature dependence of the specific heat.  These new insights into the Fermi surface  reconstruction which takes place inside the magnetically ordered states of the Kondo-lattice system CeRhIn$_5$ should  be relevant to other heavy-fermion compounds, as well as a broader range of correlated metals in proximity to Mott insulating states, such as high-$T_c$ cuprates or organic charge-transfer salts.

\begin{acknowledgments}
We acknowledge valuable discussions with S.~Arsenijevic, P.~J.~W.~Moll, Q.~Si, F.~Ronning, S.~Paschen, O.~Stockert and C.-L.~Huang.
Work at Zhejiang University was supported by the National Natural Science Foundation of China (Grants No. U1632275 and No. 11474250), National Key R\&D Program of China (Grants No. 2017YFA0303100 and No. 2016YFA0300202), and the Science Challenge Project of China (No.~TZ2016004). Work at Los
Alamos National Laboratory was performed under the auspices of the US
Department of Energy (DoE). A portion of this work was performed at the National High Magnetic Field Laboratory, which is funded by the National Science Foundation through
DMR-1157490 and the US Department of Energy and the State of Florida. J. S. appreciates funding from Basic Energy Sciences, U. S. Department of Energy FWP {\it Science in 100 T} program. We acknowledge the support of the HLD at HZDR, member
of the European Magnetic Field Laboratory (EMFL), and support by the ANR-DFG grant Fermi-NESt.
L.J. acknowledges support by the Alexander-von-Humboldt foundation.
\end{acknowledgments}


\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}


\section{Introduction}

The problem of a particle confined in a one-dimensional infinite square well
potential lies at the heart of non-relativistic quantum mechanics, being the simplest problem that
illustrates how the wave nature of bound particles implies that their energy is
quantized. So it is no surprise that introductory quantum mechanics courses present and discuss it as
a powerful pedagogical tool to introduce the students the particular features of the spectra of
quantum systems.

The generalization of this problem to three dimensions is used for the statistical description
of a gas of quantum particles (see, for example, ref. \cite{Pathria}) and thus is presented also in
statistical physics courses.

In relativistic quantum mechanics, where quantum mechanics and special relativity are combined, new features of the energy spectra
appear coming from the relativistic energy-momenta relations and spin, as well as the difficulties which arise when one tries to describe one-particle quantum systems similar to ones used in non-relativistic
quantum physics courses. The problem of a relativistic particle in a box is again a direct way to address those questions and thus is an useful
tool to discuss, in advanced quantum mechanics courses, the issues arising when one extends quantum mechanics to incorporate special relativity
(see, for instance, \cite{Li_Blinder}).

The relativistic formulation and solution of the problem of a
spin-1/2 fermion with mass $m$ confined in a one-dimensional
square well potential was done by Alberto, Fiolhais and Gil
\cite{afg} and then extended to a three-dimensional square well potential by Alberto, Das and
Vagenas \cite{adv}.

The problem of a relativistic
particle confined in an infinite square well potential is traditionally not
considered in the textbooks of Relativistic Quantum Mechanics, even in the
most comprehensive ones such as Greiner's \cite{Greiner}, which considers only the case of a finite square
well with a vector (energy coupled) potential. However, as discussed in \cite{afg,adv}, in relativistic
problems one cannot simply let the well height go to infinity, because this would lead to problems
related to the Klein paradox, i.e., the emergence of negative energy solutions and thus compromise
the one-particle picture that one wishes to retain. We address this problem again below.
Other problems particular to the relativistic problem concern the boundary conditions, which cannot be
the same as in the non-relativistic problems, as discussed in refs. \cite{VAlonso} and
\cite{menon}.

Here we extend further the problem of relativistic particle in a box to
spin-0 particles, i.e., solutions of the Klein-Gordon equation in an infinite well,
both in one- and three-dimensional space.
To this end, we use the Feshbach-Villars formalism \cite{fv}, by which one can have
a Schr\"odinger-type equation of motion, linear in the time derivative and thus having
a Hamiltonian defined for relativistic spin-0 particles.
This allows us to follow a similar procedure to obtain an infinite well solution
as was done for the Dirac equation, in particular achieving a clear separation between particle and anti-particle
solutions. Although not common, the Feshbach-Villars formalism has been applied in other contexts,
as in computing creation rates for particle-antiparticle pairs produced
by a supercritical force fields \cite{Cheng_pra}.
\par\noindent
Finally, we compare the solutions obtained with those for the Dirac equation and draw the conclusions.

\section{Klein-Gordon equation in Feshbach-Villars formalism}
The Feshbach-Villars formalism (FVF) consists essentially in obtaining a first-order differential equation in time from
the original second-order Klein-Gordon equation.
This is achieved by defining two wave functions $\varphi$ and $\chi$.
Following the original Feshbach-Villars paper \cite{fv}
but with a slight different notation and definitions, and also inspired by the textbook of Greiner \cite{Greiner}, one
may write
\bea
\psi &=& \varphi + \chi  \la{psia} \\
i\hbar\frac{\partial\psi}{\partial t}&=&mc^2(\varphi - \chi) \la{psib}
\eea
where $\psi$ is the conventional Klein-Gordon wavefunction.
If now one defines the two-component spinor
\beq
\Psi = \left( \begin{array}{c}
\varphi \\[2mm]
\chi
 \end{array} \right)\ ,
\eeq
the free Klein-Gordon equation can be written as a Schr\"odinger-type equation of motion
\bea
 i \hbar \frac{\partial \Psi}{\partial t}&=&H\Psi \la{kgeq1}\\
H &=& \le( \tau_3 + i\tau_2 \ri) \frac{\hat{\vec p}^{\;2}}{2m} + mc^2 \tau_3 ~, \la{kgham1}
\eea
where $\tau_k$ ($k=1,2,3$) are Pauli matrices and $\hat{\vec p}$ is the momentum operator,
given as $\hat{\vec p}=-i\hbar\nabla$.
Note that the Hamiltonian (\ref{kgham1}) is similar to the one of a non-relativistic particle
plus its rest energy.
Equation (\ref{kgeq1}) is equivalent to two coupled equations for $\varphi$ and $\chi$
\bea
 i \hbar \frac{\partial \varphi}{\partial t} &=&\frac{\hat{\vec p}^{\;2}}{2m}(\varphi + \chi)+mc^2\varphi   \\
 i \hbar \frac{\partial \chi}{\partial t}&=&-\frac{\hat{\vec p}^{\;2}}{2m}(\varphi + \chi)-mc^2\chi \ .
\eea
One can check, from these last equations and equations (\ref{psia}) and (\ref{psib}), that one gets the free Klein-Gordon equation for $\psi$
\beq
\frac1{c^2}\frac{\partial^2 \psi}{\partial t^2}-\nabla^2\psi+\frac{m^2c^2}{\hbar^2}\psi=0 \ .
\eeq
\par\noindent
The conserved charge and current assume the following forms
\bea
\rho &=&\varphi^\star \varphi - \chi^\star\chi =\Psi^\dagger \tau_3 \Psi \\
\notag
\vec J &=&
\frac{\hbar}{2im}\le[
(\varphi^* + \chi^*)\nabla(\varphi + \chi)
- (\varphi + \chi)\nabla(\varphi^* + \chi^*)\ri]\\
&=& \frac{\hbar}{2im}\le[
\Psi^\dagger\tau_3 \le( \tau_3 + i\tau_2 \ri) \nabla\Psi
- \le( \nabla \Psi^\dagger \ri) \tau_3 \le( \tau_3 + i \tau_2 \ri) \Psi
\ri]
\ .\la{kgcurrent}
\eea
\par\noindent
Similarly to spin-1/2 particles, there is a charge conjugated spinor given by
\beq
\Psi_c=\left( \begin{array}{c}
\chi^* \\[2mm]
\varphi^*
 \end{array} \right)=\tau_1\Psi^*
\eeq
which changes particle into anti-particle solutions and vice-versa. One can check
that the charge conjugated density and current are
\bea
\rho_c &=&-\rho  \\
\vec J_c &=&- \vec J\ ,\la{kgcurrent_c}
\eea
meaning that these should be interpreted as a charge density and charge current, as is well-known
to be the case for the Klein-Gordon equation. One may also note that here $\tau_3$ is related to the particle charge, so
that the physical significance of Pauli matrices $\tau_i$  is similar to that of the Pauli matrices used, for instance,
for the isospin degree of freedom of nucleons.
\section{Solution of the Klein-Gordon Equation for squared well potentials}
The solution of the free Klein-Gordon equation (\ref{kgeq1}) in the FVF which has good linear momentum, i.e.,
which is an eigenstate of the linear momentum operator, is given by
\cite{fv}
\beq
\Psi^{\pm}_{\vec p}(\vec x,t )=A_{\pm} \left( \begin{array}{c}
\varphi_0^{\pm} (\vec p ) \\[2mm]
\chi_0^{\pm} (\vec p)
 \end{array} \right) e^{ i(\mp E_p t-\vec p \cdot \vec x)/\hbar} \equiv
  A_{\pm}\Psi_0^{\pm}(\vec p)e^{ i(\mp E_p t-\vec p \cdot \vec x)/\hbar}
\label{psi_p}
\eeq
where $\pm$ denotes the positive (negative) energy solutions, i.e.,
$E=\pm E_p$, with $E_p = \sqrt{p^2c^2 + m^2 c^4}$, $A_{\pm}$ are normalization constants, and functions
$\varphi_0^{\pm} (\vec p )$ and $\chi_0^{\pm} (\vec p )$ are written as
\bea
\label{phi0}
\varphi_0^{\pm} (\vec p) &=& \frac{\pm E_p+mc^2}{2\sqrt{mc^2 E_p}} \\
\label{chi0}
\chi_0^{\pm} (\vec p) &=& \frac{mc^2\mp E_p}{2\sqrt{mc^2 E_p}}~.
\eea
Note that $\varphi_0^\pm(\vec p )$ and $\chi_0^\pm(\vec p )$ depend on $\vec p$ only through its magnitude
$|\vec p|$ and $\big[\varphi_0^{\pm}(\vec p )\big]^2-\big[\chi_0^{\pm}(\vec p )\big]^2 =\pm 1$. We also have
$\Psi^{\pm}_{\vec p}(\vec x,t )_c=\Psi^{\mp}_{-\vec p}(\vec x,t )$.
\subsection{Solution for one-dimensional squared well potential}
In order to obtain the solution of the Klein-Gordon equation for an infinite squared well potential, one must tackle the
problem that arises with the one-particle description in Relativistic Quantum Mechanics when particles are subject to
very strong external potentials. This shows up clearly when one computes the transmission of a plane wave through a
vector (energy-coupling) potential barrier, known as the Klein paradox. While it takes a different form for spin-1/2 particles than for spin-0 particles, one common feature is that in both cases a
strong vector potential leads to a non-zero transmission to a classically forbidden region, which is related to particle/anti-particle pair
production. A good historical review can be found in
ref. \cite{cd_K_P}. In order to avoid this problem and retain the one-particle quantum description of spin-0 particles with wave functions which are
solutions of the Klein-Gordon equation, one must consider the one-dimensional potential defined in a one-dimensional box
${\mathbb L} \equiv \{ 0 \leq x \leq L \}$ with length $L$ as a potential invariant under Lorentz transformations, i.e.,
a position dependent mass term, as follows
\beq
m(x)=\left\{
\begin{array}{lc}
  m &,\  x \in {\mathbb L}\\[2mm]
  M \to \infty &,\  x \notin {\mathbb L}
\end{array}\right.\ ,
\label{1d_pot}
\eeq
as was done in the Dirac equation case \cite{afg}. In this case, one can consider the time-independent solutions of the Klein-Gordon equation in the FVF, that is, one has
\beq
H\tilde\Psi(\vec x)=E\tilde\Psi(\vec x) \quad {\rm where}\quad \Psi(\vec x,t)=e^{-i E t/\hbar}\tilde\Psi(\vec x) \  .
\eeq
\par\noindent
From equation (\ref{psi_p}), the general solution for a time-independent positive-energy solution in this one-dimensional box is given by,
with $p>0$,
\beq
\Psi_B(x)=\left\{
\begin{array}{cc}
 A \Psi_0^{+}(p)e^{ ip x/\hbar}+B\Psi_0^{+}(p)e^{- i p x/\hbar} &,\  x \in {\mathbb L}\\[2mm]
0 &\ ,  x \notin {\mathbb L}
\end{array}\right.\ .
\label{1d_wf}
\eeq
The continuity of the wave function at the box borders yields the box boundary conditions
$\Psi_B (0) = 0 = \Psi_B(L)$. Together with the normalization
$\int_{\mathbb L} \rho \,dx =1$, one has for the wave function inside the box
\beq
\Psi_B (x)= \sqrt{\frac{2}{L}}\Psi_0^+(p) \sin(p\,x/\hbar) \quad {\rm with}\quad p\,L=n\pi\hbar\quad n=1,\ldots \ .
\eeq
This gives rise to the quantized energy
\bea
E_B(n) &=&\sqrt{\frac{n^2\pi^2\hbar^2 c^2}{L^2} + m^2c^4 } \\
&\approx&
mc^2 + \frac{n^2 \pi^2 \hbar^2}{2m L^2} + \dots
\eea
where a small $p/(mc)=n\pi\hbar/(L mc)$ expansion has been made in the last line, with the second term being that
which arises by solving the Schr\"odinger equation, $mc^2$ is the rest energy,
and $\dots$ represent second and higher order terms.

The usual Klein-Gordon wavefunction $\psi$ and its time derivative, given by (\ref{psia})
and (\ref{psib}) respectively, can be recovered noting that in those equations
\bea
\varphi&=&\sqrt{\frac{2}{L}}\,\varphi_0^+(p_n)\, e^{-i E_B(n) t/\hbar}\sin(p_n x/\hbar)\\
\chi&=&\sqrt{\frac{2}{L}}\,\chi_0^+(p_n)\, e^{-i E_B(n) t/\hbar}\sin(p_n x/\hbar)
\eea
in which $p_n=n\pi\hbar/L$ and $\varphi_0^+(p_n)$ and $\chi_0^+(p_n)$ are given
respectively  by
(\ref{phi0}) and (\ref{chi0}) with $E_p$ replaced by $E_B(n)$.

Finally, the current as given by (\ref{kgcurrent}) vanishes in this case, since it is composed
of a left and a right moving wave with symmetric amplitudes,
satisfying at the same time current conservation throughout all space.
\subsection{Solution for the three-dimensional squared well potential}
For a three dimensional box defined by ${\mathbb V} \equiv \{ 0 \leq x_i \leq L_i \},~i=1,2,3$,
where $L_i$ are the lengths of the box in each dimension,
the scalar infinite well potential would be given by
\begin{equation}
m(\vec x)=\left\{
\begin{array}{cc}
  m &,\  \vec x\in {\mathbb V}\\[2mm]
  M \to \infty &,\   \vec x\not\in {\mathbb V}
\end{array}\right.\ .
\label{3d_pot}
\end{equation}
Since the potential inside the box is constant, we expect that the wave function in the box would be a linear combination of
wave functions with good momenta of the type (\ref{psi_p}). One may note also that the potential (\ref{3d_pot}) can be written
as separable potential in terms of the one-dimensional well potentials: $m(\vec x)=m+(M-m)\prod_{i=1}^3[\theta(x_i-L_i)+\theta(-x_i)]$, where $\theta(x)$ is the Heaviside step function.
This means that the space-dependent part of the wave function may be written as a product of a linear combination of
one-dimensional plane waves travelling both ways in each direction, or a linear combination of the eight plane waves
\beq
e^{i \sum_{i=1}^3 \epsilon_i p_i x_i/\hbar}\quad {\rm where}\quad \epsilon_i=\pm 1\ .
\eeq
Since the magnitude of the momenta of all of these plane waves $(\pm p_1, \pm p_2, \pm p_3)$
is the same, one can write the general time-independent positive-energy wave function as
\beq
\Psi_B(\vec x)=\left\{
\begin{array}{cc}
 \displaystyle
 \Psi_0^{+}(\vec p)\sum_{\vec\epsilon}A_{\vec\epsilon}\,
 e^{i \sum_{i=1}^3 \epsilon_i p_i x_i/\hbar} &,\   \vec x \in {\mathbb V}\\[2mm]
0 &,\   \vec x \notin {\mathbb V}
\end{array}\right.\ ,
\label{3d_wf}
\eeq
where $\vec\epsilon$ represents all the 8 sets of values of $(\epsilon_1,\epsilon_2,\epsilon_3)$,
that is, $(\pm 1,\pm 1,\pm 1)$ and $A_{\vec\epsilon}$ is
the amplitude for each plane wave.
Note the formal similarity with the Dirac wave function for a relativistic spin-$1/2$
particle in a three dimensional box, equation (33) of \cite{dva} and equation (24) of \cite{adv}.

The boundary conditions and normalization are now
\bea
\Psi_B (0,0,0) &=& 0 = \Psi_B (L_1,x_2,x_3) = \Psi_B (x_1,L_2,x_3) = \Psi_B (x_1,x_2,L_3) \\
\int_{\mathbb V} \rho \,d^3\vec x &=& 1
\eea
giving rise to the interior wave function
\beq
\Psi_B(\vec x)=\sqrt{\frac{2^3}{V}}~\Psi_0^+ (\vec p)\prod_{i=1}^3 \sin(p_i x_i/\hbar)\ ,
\eeq
where $V=L_1 L_2 L_3$ is the volume of the box and one has $p_i L_i/\hbar = n_i \pi$, with $n_i=1,\ldots$, and $i=1,2,3$.

\bea
\label{KG_box3D}
E_B(\vec n) &=& \sqrt{\sum_{i=1}^3 \frac{ n_i^2\pi^2\hbar^2 c^2}{L_i^2} + m^2 c^4 } \\
&\approx& mc^2 + \sum_{i=1}^3 \frac{n_i^2 \pi^2\hbar^2}{2mL_i^2}+ \dots
\eea
where $\vec n=(n_1,n_2,n_3)$, and the wave functions $\varphi$ and $\chi$ are now
\bea
\varphi&=&\sqrt{\frac{2^3}{V}}\,\varphi_0^+(\vec p(\vec n))\, e^{-i E_B(\vec n) t/\hbar}
\prod_{i=1}^3 \sin(p_i(\vec n)) x_i/\hbar)\\
\chi&=&\sqrt{\frac{2^3}{V}}\,\chi_0^+(\vec p(\vec n))\, e^{-i E_B(\vec n) t/\hbar}
\prod_{i=1}^3 \sin(p_i(\vec n)) x_i/\hbar)
\eea
in which $p_i(\vec n)=n_i\pi\hbar/L_i$ while  $\varphi_0^+(\vec p(\vec n))$ and $\chi_0^+(\vec p(\vec n))$ are
given by (\ref{phi0}) and (\ref{chi0}) respectively, with $E_p$ replaced
by $E_B(\vec n)$. Once again the charge current $\vec J$ is zero.
\section{Comparison between the energy spectra for spin-0 and spin-1/2 relativistic particles in a box}
\subsection{One-dimensional box}
The energy solutions of the Klein-Gordon and Dirac equation for particles in a one-dimensional box of size $L$ are given,
in terms of the allowed wave numbers $k_n$, by (see \cite{afg})
\bea
\label{KG_BC_box1D}\hbox{Klein-Gordon}\ &&E_B(n)= \sqrt{\hbar^2 c^2 k_n^2 + m^2c^4 }\  ,\  \ k_n=n\pi/L\ ,\\
\label{Dirac_BC_box1D}\hbox{Dirac}\ &&E_B(n)= \sqrt{\hbar^2 c^2 k_n^2 + m^2c^4 }\  ,\ \ \tan(k_n L)=-\frac{\hbar k_n}{mc} \ .
\eea
In the case of the Dirac equation, $k_n$ denotes the $n$th solution (excluding the trivial one $k=0$) of the
transcendental equation. It is convenient to write these quantities in terms of dimensionless quantities and consider the
scaled kinetic energy ${\cal T}=T/(mc^2)
=E/(mc^2)-1$
instead of the total energy. Defining $x_n=k_n\lambda_C$ and $L_C=L/\lambda_C$, where $\lambda_C=\hbar/(mc)$ is the Compton wavelength
of the relativistic particle of mass $m$, one gets
\bea
\hbox{Klein-Gordon}\ &&{\cal T}^{KG}_B(n)= \sqrt{x_n^2 + 1}-1\ ,\quad x_n=n\pi/L_C\ ,\\
\label{Dirac_box1D}{\rm Dirac}\ &&{\cal T}^D_B(n)= \sqrt{x_n^2 + 1}-1\ ,\quad\ \tan(L_C x_n)=-x_n \ .
\eea
One sees that the main difference between Klein-Gordon and Dirac spectra is the value of the quantity $x_n$.
Writing the transcendental equation for $x_n$ in (\ref{Dirac_box1D}) as
\beq
\tan(y_n)=-\frac{y_n}{L_C}\,, \qquad y_n=x_n L_C\ ,
\eeq
it is clear that, for lengths of the box much bigger than the Compton wavelength, i.e., $L_C\gg 1$, one has $y_n\sim n\pi$ and thus $x_n\sim n\pi/L_C$. Since, on the other hand, $x_n$, at least for values
of $n$ not too high, would be small in this case, we would be in the non-relativistic limit. Therefore, one can state that the effect of spin on the
energy of a relativistic particle in a box is itself coming from relativity, or the covariance of the Dirac equation, since this is how the Dirac plane wave spinor is constructed, and thus the origin of the transcendental equation above. This is similar to what happens to a spin-1/2 particle in a spherical box, for which the spin-orbit coupling disappears for big (non-relativistic) boxes \cite{afo}.

The difference between  the two equations (\ref{KG_BC_box1D}) and (\ref{Dirac_BC_box1D}) can be traced back to
the boundary conditions obeyed by the wave function for a particle-in-a-box problem in Klein-Gordon and Dirac equations,
respectively. In the case of Klein-Gordon equation, being a second-order differential equation in the coordinates, one
requires, as in the Schr\"odinger equation, that the wave function vanishes at the box borders because of the continuity of the wave function.
However, this cannot be done in the Dirac equation, as mentioned before and discussed in refs. \cite{afg} and \cite{adv}.
Indeed, being a first-order differential equation in the coordinates, the former requirement cannot be made, and it is replaced
by a boundary condition which is equivalent to demand that the current flux $\vec j\cdot\hat n$ vanishes at the
 box borders ($\hat n$ is the outward unit vector at the borders). Of course, this happens in the Klein-Gordon equation
as well -- the current vanishes at the borders -- but in the Dirac equation the current $\vec j$ is fundamentally
different from the Klein-Gordon one because of the spinor structure of the wave function.
Therefore, this structure, related to the spin 1/2 nature of the Dirac particles, affects the energy spectrum
for a spin-1/2 relativistic particle in a box problem. It turns out that both boundary condition are equivalent in the non-relativistic limit, because in this limit one recovers the non-relativistic solutions of the same problem (see \cite{afg, adv}). This spin effect in relativistic problems does not happen necessarily with other interactions, namely with the interaction with a pure vector electromagnetic field (except for the Stern-Gerlach term) or for certain combination of scalar and vector potentials \cite{prc_75_047303}.

For the electron, $\lambda_C\sim 3.86\times 10^{-13}\ {\rm m}=3.86\times 10^{-3}$ {\AA},
 a relatively small value at atomic scale, but one could envisage an electron confined in smaller boxes
than an atom-sized box, given a sufficient high electric field. For spin-0 mesons like the charged pions $\pi^\pm$,
 one has $\lambda_C\sim 1.41\ {\rm fm}=1.41\times 10^{-5}$ {\AA}.

In the following figure we plot the values of ${\cal T}^{KG}_B(n)$ and  ${\cal T}^{D}_B(n)$ for $n=1,\ldots,4$,  in logarithmic scale,
for several values of $L_C$.
The non-relativistic value ${\cal T}^{nr}_B(n)=x_n^2/2$ is presented for the highest value of $L_C=300$, which is about $1.15$ {\AA} for an electron.
\begin{figure}[!ht]
\begin{center}
\includegraphics[width=12cm]{fig1_KG.eps}
\end{center}
\caption{Scaled kinetic energy spectrum of the first four solutions of Klein-Gordon and Dirac equations
for four values of the ratio $L_C$ for a one-dimensional box.}
\label{Fig1}
\end{figure}
\subsection{Three-dimensional box}
For the three-dimensional infinite well we have now for the Klein-Gordon and Dirac energies \cite{adv}
\bea
\hbox{Klein-Gordon}\ && E_B(\vec n)= \sqrt{\hbar^2 c^2 |\vec k(\vec n)|^2 + m^2c^4 }\ ,\ k_i(\vec n)=\frac{n_i\pi}{L_i}\ ,\ i=1,2,3\ ,\\
{\rm Dirac}\ && E_B(\vec n)= \sqrt{\hbar^2 c^2 |\vec k(\vec n)|^2 + m^2c^4 }\ ,\nonumber\\
&&\hspace*{0.0cm}\tan(k_i(\vec n) L_i)=-\frac{2(E_B(\vec n)+mc^2)\hbar k_i(\vec n)\,c}
{\hbar^2c^{2} k_i^2(\vec n)-(E_B(\vec n)+mc^2)^2}\ ,\ i=1,2,3 \ .
\eea
As in the previous section, $\vec n$ denotes the set of quantum numbers $(n_1,n_2,n_3)$, all non-zero positive integers.
In the case of the Dirac equation, they are used to
label the different solutions of the three coupled transcendental equations, as discussed in \cite{adv}.
Considering now a cubic box, i.e., $L_1=L_2=L_3=L$, and, as in the one-dimensional case, using the scaled kinetic energy and defining
$x_i(\vec n)=k_i(\vec n)\lambda_C$, $L_C=L/\lambda_C$, we get
\bea
\hspace*{-.7cm}\hbox{Klein-Gordon}&&\hspace*{-.3cm}{\cal T}^{KG}_B(\vec n)= \sqrt{|\vec x(\vec n)|^2 + 1}-1\,,
\ x_i(\vec n)=\frac{n_i\pi}{L_C}\,,\  i=1,2,3\ ,\\
\hbox{Dirac}&&\hspace*{-.3cm}{\cal T}^D_B(\vec n)= \sqrt{|\vec x(\vec n)|^2 + 1}-1\ ,\nonumber \\
&&\tan(x_i(\vec n) L_C)=\frac{2({\cal T}^D_B(\vec n)+2)x_i(\vec n)}
{x_i^2(\vec n)-({\cal T}^D_B(\vec n)+2)^2}\ ,\ i=1,2,3 \ .
\eea
Again, defining $y_i(\vec n)=x_i(\vec n)L_C$ one has for the transcendental equations
\beq
\tan(y_i(\vec n))=L_C\frac{2({\cal T}^D_B(\vec n)+2)y_i(\vec n)}
{y_i^2(\vec n)-L_C^2({\cal T}^D_B(\vec n)+2)^2} \ .
\eeq
When $L_C\gg 1$, one has $\tan(y_i(\vec n))\sim 0$ and thus $x_i\sim n_i\pi/L_C$, and thus, as before, we get the Klein-Gordon
(and in this case, non-relativistic) result. It is interesting to remark that, if one $x_i(\vec n)$ is much bigger than the others, then
$x_i(\vec n)\sim|\vec x(\vec n)|$ and the respective equation becomes
\bea
\tan(x_i(\vec n) L_C)&=&\frac{2({\cal T}^D_B(\vec n)+2)x_i(\vec n)}
{x_i^2(\vec n)-({\cal T}^D_B(\vec n)+2)^2}\nonumber\\
&=& \frac{2({\cal T}^D_B(\vec n)+2)x_i(\vec n)}
{x_i^2(\vec n)-|\vec x(\vec n)|^2-2({\cal T}^D_B(\vec n)+2)}\sim -x_i(\vec n)
\eea
exactly as the one-dimensional case.

For the 3-dimensional box, there is an energy degeneracy regarding the levels with quantum numbers $(n_1,n_2,n_3)$: any permutation
of the $n_i$'s gives rise to the same energy, as it is evident from the Klein-Gordon energy expression (\ref{KG_box3D}) for a
cubic box and was already noted in \cite{adv} for a Dirac particle in a cubic box as well. Therefore, one has 3 types of degeneracies:
1-fold for $n_1=n_2=n_3$; 3-fold when two $n_i$'s are the same but the third is different;
$3!=6\,$-fold degeneracy when the $n_i$'s are all different.

The following figure is similar to Fig.~\ref{Fig1} for the three-dimensional box. In this case the first four distinct energy levels correspond to $\vec n$ being equal to
$(1,1,1)$, $(1,1,2)$, $(1,2,2)$ and $(1,1,3)$ in increasing energy order. Because of the degeneracies mentioned above, one is actually depicting 10 energy levels. Actually, in the Dirac case, one still has an additional 2-fold degeneracy for each level, corresponding to the two independent spin polarizations.

\begin{figure}[!ht]
\begin{center}
\includegraphics[width=12cm]{fig2_KG.eps}
\end{center}
\caption{Scaled kinetic energy spectrum of the first four solutions of Klein-Gordon and Dirac equations
for four values of the ratio $L_C$ for a three-dimensional box.}
\label{Fig2}
\end{figure}

Several aspects are common to the Klein-Gordon and Dirac on
three-dimensional box spectra: the energy levels for spin-1/2 are in general lower than the ones for spin-0 particles and,
for larger boxes, measured in units of the Compton wavelength, they approach the non-relativistic values, as could be expected, since larger boxes mean smaller momenta for the particle confined.
In the three-dimensional case, there are more levels up to a certain energy and they are closer to each other than in the one-dimensional case.
This can be understood by the fact that with increasing $n$, the number of combinations of $(n_1,n_2,n_3)$ such that
$n_1^2+n_2^2+n_3^2\leq n^2$ increases more than linearly with $n$, actually as $n^3$ for sufficiently high $n$.

As discussed in \cite{adv}, the energy levels computed for a three-dimensional box allows us to compare the density of states of a
relativistic fermion gas to the density of a non-relativistic fermion gas. Since
the separation between $k_i$ values is smaller than $\pi/L_i$ for volumes $V\sim \lambda_C^3$, the relativistic density is higher
than the non-relativistic one. This can also be seen in Fig.~\ref{Fig2}. For a spin-0 gas, in spite of the fact the separation between allowed
$k_i$ value is identical with the non-relativistic case, one must take into account the Bose-Einstein statistics, by which one has,
at zero temperature, as many states as one wishes in the same energy state.

\section{Conclusions}
In this paper, we have solved the problem of a relativistic spin-0 particle in a one- and three-dimensional
box using the Klein-Gordon equation in the Feshbach-Villars formalism. By using position-dependent mass for the
infinite square well potential we were able to avoid the Klein paradox problem and thus retain a one-particle description
of spin-0 particles.
We compare these solutions with those previously found for the Dirac equation with the same potentials.
We find that the spin has indeed an effect on the energy spectra for relativistic spin-1/2 particles as compared with
relativistic spin-0 particles with the same mass, decreasing the energy of the former compared in the latter.
As expected, both kinds of particles tend to the same non-relativistic spectrum, if the size of the box is sufficiently
large in units of Compton wavelength.

\vskip1cm
\noindent\textbf{Acknowledgments}
\vskip.1cm
\par\noindent
This work was supported in part by the Natural
Sciences and Engineering Research Council of Canada.