
\subsection{Details on the BCH expansion}
Given the two-step form of our drive, we can derive an effective Hamiltonian in the high-frequency regime with a Baker-Campbell-Hausdorff (BCH) formula. If we choose $t_0=0$ as our initial time,  the effective (Floquet) Hamiltonian $\operator{H}_F$ is defined via
\begin{equation}
\operator{F} =  \hat{U}\left(\tau , 0\right) \equiv e^{-i \operator{H}_F \tau} = e^{-i \operator{H}_1 \frac{\tau}{2}} e^{-i \operator{H}_2 \frac{\tau}{2}} \;
\end{equation}
where $H_{1/2}$ are 
\begin{equation}
\begin{aligned}
H_{1/2} = \sum_i (\frac{m}{\sqrt{L}} (a_i^{\dagger} a_1^{\phantom{\dagger}} + a_1^{\dagger} a_i^{\phantom{\dagger}})+\epsilon_i^{\phantom{\dagger}} a_i^{\dagger} a_i^{\phantom{\dagger}}) \pm A a_1^{\dagger} a_1^{\phantom{\dagger}} \;.
\end{aligned}
\end{equation}
Using the standard BCH formula we obtain, up to second order in $1/\omega$,%
\begin{equation}
\begin{split}
H_F= \frac{H_1 + H_2}{2} +& i2\pi\frac{\lbrack H_1, H_2 \rbrack }{8\omega}  + \\ 
+ \frac{(2\pi)^2}{96\omega^2} &(- \lbrack H_1 \lbrack H_1, H_2\rbrack \rbrack + \lbrack H_2 \lbrack H_1, H_2 \rbrack \rbrack ) \;.
\end{split}
\end{equation}
The commutators are easily computed as
\begin{equation*}
\begin{aligned}
&\lbrack H_1 , H_2 \rbrack 
= - 2 A \frac{m}{\sqrt{L}} \sum_i ( a_i ^{\dagger} a_1 ^{\phantom{\dagger}} - a_1 ^{\dagger} a_i    ) \;,\\
&\lbrack H_1 \lbrack H_1, H_2\rbrack \rbrack - \lbrack H_2 \lbrack H_1, H_2 \rbrack \rbrack =
4 \frac{ A^2 m}{\sqrt{L}} (a_1^{\dagger} a_i^{\phantom{\dagger}} + \mathrm{h.c.} ) \;,
\end{aligned}
\end{equation*}
yielding the Floquet Hamiltonian
\begin{equation}
\begin{aligned}
H_F =&   \sum_i \left(\frac{m}{\sqrt{L}} \left(a_i^{\dagger} a_1^{\phantom{\dagger}} +\mathrm{h.c.}\right)+\epsilon_i^{\phantom{\dagger}} a_i^{\dagger} a_i^{\phantom{\dagger}} \right)+\\ & -\frac{2\pi}{4\omega} i A  \frac{m}{\sqrt{L}} \sum_i \left( a_i ^{\dagger} a_1 ^{\phantom{\dagger}} - \mathrm{h.c.}   \right) + \\&- \frac{ A^2 (2\pi)^2 }{24 \omega^2}\frac{m}{\sqrt{L}} \left(a_1^{\dagger} a_i^{\phantom{\dagger}} + \mathrm{h.c.}\right)  + \mathcal{O}\left(\frac{1}{\omega^3}\right).
\end{aligned}
\end{equation}
As we mention in the main text, the first order correction gives an imaginary hopping term, breaking time-reversal symmetry (TRS). However, this is only apparent, because this term depends on the initial time, while the long-time dynamics should not. Changing the initial time is equivalent to a gauge transformation~\cite{Eckardt_NJP15} on $\operator{H}_F$. Indeed, if one, for example, moves the initial time at $\tau/2$, this term changes its sign. Many quantities must be gauge independent, e.g. should not depend on the initial time of the driving. One example are the quasienergies, which we have used as a measure of localization through their gap statistics. If one is instead interested in an expansion which is gauge invariant, a tripartition of the Floquet operator must be employed~\cite{Goldman_PRX14,Eckardt_NJP15}: $\operator{F}_{t_0} = \operator{S}^{\dagger}\exp \left(-i\operator{\tilde{H}}_{F} \tau \right)\operator{S}^{\phantomsection{\dagger}}(t_0)  $. $\operator{S}$ are called kick operators, while the new Hamiltonian $\operator{\tilde{H}}_{F}$ does not depend on the initial time $t_0$ anymore. For our model, this procedure essentially consists of keeping in $\operator{\tilde{H}}_F$ only the second order real term. On the other hand, the kick operators will contain, among others, the first order term. In this way, one can separate between long time effects and intraperiod ones. For simplicity, we display in the main text the BCH result, as the localization behaviour does not change qualitatively.



