Universal fluctuations of Floquet topological invariants at low frequencies

TL;DR

This paper reveals that Floquet topological invariants in low-frequency driven systems exhibit universal Gaussian fluctuations with a width scaling as 1/√Ω, indicating structured universal behavior even near the adiabatic limit.

Contribution

It demonstrates analytically and numerically that topological invariants in low-frequency Floquet systems fluctuate universally, with a Gaussian distribution and a specific scaling law.

Findings

Topological invariants follow a Gaussian distribution at low frequencies.

The width of fluctuations scales as 1/√Ω.

The quasienergy gap remains finite and scales as Ω².

Abstract

We study the low-frequency dynamics of periodically driven one-dimensional systems hosting Floquet topological phases. We show, both analytically and numerically, in the low-frequency limit $\Omega\to0$, the topological invariants of a chirally-symmetric driven system exhibit universal fluctuations. While the topological invariants in this limit nearly vanish on average over a small range of frequencies, we find that they follow a universal Gaussian distribution with a width that scales as $1/\sqrt{\Omega}$. We explain this scaling based on a diffusive structure of the winding numbers of the Floquet-Bloch evolution operator at low frequency. We also find that the maximum quasienergy gap remains finite and scales as $\Omega^2$. Thus, we argue that the adiabatic limit of a Floquet topological insulator is highly structured, with universal fluctuations persisting down to very low frequencies.

Figures

Equations

$$\nu =\mathcal{W}[d]\equiv \frac{1}{2\pi i}{\int }_{-\pi }^{\pi }\frac{\partial }{\partial k}\ln (d_k)dk.$$

$${\nu }^{♮}=\mathcal{W}[B]\text{and}{\nu }^{♭}=\mathcal{W}[D].$$

$$\cos \frac{2\pi {ϵ}_k}{\Omega }=\cos \frac{2\pi {\bar{E}}_k}{\Omega }{\cos }^2\frac{{\theta }_k}{2}+\cos \frac{2\pi {\breve{E}}_k}{\Omega }{\sin }^2\frac{{\theta }_k}{2},$$

$$D_k=e^{-i{\theta }_k/2}\left(\cos \frac{\pi {\bar{E}}_k}{\Omega }\cos \frac{{\theta }_k}{2}+i\cos \frac{\pi {\breve{E}}_k}{\Omega }\sin \frac{{\theta }_k}{2}\right)$$

$${\tilde{B}}_k=e^{-i{\theta }_k/2}\left(\sin \frac{\pi {\bar{E}}_k}{\Omega }\cos \frac{{\theta }_k}{2}+i\sin \frac{\pi {\breve{E}}_k}{\Omega }\sin \frac{{\theta }_k}{2}\right).$$

$$\Gamma U(\tau ,\tau +T)\Gamma =U^{-1}(\tau ,\tau +T)$$

$$H_{\text{eff}}(\tau )\equiv \frac{i}{T}\ln U(\tau ,\tau +T)=\left[\begin{array}{cc}0& h(\tau )\\ h^{†}(\tau )& 0\end{array}\right],$$

$${\nu }_{\tau }=\mathcal{W}[h(\tau )]=\frac{1}{2\pi i}{\int }_{-\pi }^{\pi }dk\frac{\partial }{\partial k}\ln [\text{mbox}det{h}_k(\tau )].$$

$${\nu }^{♮}=\frac{1}{2}\left(\nu +{\nu }^{\prime }\right),{\nu }^{♭}=\frac{1}{2}\left(\nu -{\nu }^{\prime }\right),$$

$$\hat{H}(t)=\sum _r[w+(-1{)}^r\delta (t)]{\hat{c}}_{r+1}^{†}{\hat{c}}_r^{}+\text{h.c.},$$

$$U_k(\pi /\Omega )={\Pi }_{j=0}^n{e}^{-i{h}_k(j\pi /(\Omega n))\pi /(\Omega n)}.$$

$$\hat{H}=\sum _{r=1}^{N-1}\left(w{\hat{c}}_{r+1}^{†}{\hat{c}}_r^{}+\Delta {\hat{c}}_{r+1}^{†}{\hat{c}}_r^{†}+h.c.\right)-\mu \sum _{r=1}^N{\hat{c}}_r^{†}{\hat{c}}_r^{},$$

$$h_k\to R{h}_k{R}^{†}=\left(2w\cos k-\mu \right){\sigma }_x+|\Delta |\sin k{\sigma }_y,$$

$$H=-\sum _x$$

$$-\Delta \left[(1+\delta /w){\hat{c}}_{B,x}^{†}{\hat{c}}_{A,x}^{†}+(1-\delta /w){\hat{c}}_{A,x+1}^{†}{\hat{c}}_{B,x}^{†}\right]-\frac{\mu }{2}\left({\hat{c}}_{A,x}^{†}{\hat{c}}_{A,x}^{}+{\hat{c}}_{B,x}^{†}{\hat{c}}_{B,x}^{}\right)\}+\text{h.c.},$$

$$h_k=\left[\begin{array}{cccc}-\mu & z_k& 0& {\Delta }_k\\ z_k^{\ast }& -\mu & -{\Delta }_k^{\ast }& 0\\ 0& -{\Delta }_k& \mu & -z_k\\ {\Delta }_k^{\ast }& 0& -z_k^{\ast }& \mu \end{array}\right],$$

Full text

The behavior of a periodically driven system can be qualitatively different from its equilibrium behavior. Manifestations of such behavior in classical physics include resonances, dynamical stabilization of new steady states, and the period-doubling approach to chaos kapitza1951 ; Fei78a ; ref1d . In quantum systems, the effective Floquet dynamics of a driven systems has been employed as a powerful way to engineer designer Hamiltonians, e.g. by using laser sequences in cold atomic gases. In this way, novel phases of matter have been proposed and realized JotMesDes14a ; AidLohSch15a ; zhang2017x ; choi2017x ; Eck17a .

More recently, it has been understood that a driven system can also exhibit essentially non-equilibrium topological phases, dubbed Floquet topological phases oka2009 ; kitagawa2011 ; lindner2011 ; jiang2011 . Drive parameters, such as the frequency $\Omega$ or the shape of the drive (“drive protocol”) have been proposed zhenghao2011 ; dora2012 ; KunSer13a ; KunFerSer14a ; perez2014 ; ref1a ; ref1c ; titum2015 ; KunFerSer16a and used in the lab rechtsman2013 ; WanSteJar13a ; SieLuiLee17a ; ChePanWan18a to engineer a rich array of topological phases not possible in equilibrium systems. The non-equilibrium dynamics at large frequencies is relatively well understood, e.g. within rotating-wave approximation, as a renormalization of the equilibrium parameters of the system blanes2009 ; rahav2003 ; rahav2003b ; bukov2015 ; EckAni15a . The low-frequency regime, on the other hand, remains largely unexplored GomPla13a . This is the relevant regime in solid-state systems driven by ac potentials lindner2011 . It is also important as a way to reduce unwanted heating in the system LazDasMoe14a ; DAlRig14a ; WeiKna17a ; StrEck16a . At a more basic level, it relates to the adiabatic limit as $\Omega\to 0$. Numerical studies have reported nonzero Floquet topological invariants as frequency is lowered KunSer13a ; LiuLevBar13a ; TonAnGon13a ; GomDelPla14a ; MikKitYas15a ; ref1b . This raises questions on the nature of adiabatic limit in Floquet topological phases.

In this paper, we study the low-frequency limit of one-dimensional model driven systems that exhibit a rich Floquet topological phase diagram jiang2011 ; KunSer13a . Assuming the driven systems are chirally symmetric Asbóth et al. (2014); ref1b ; dallago2015 , we derive analytical expressions for the Floquet topological invariants and evaluate them numerically over several decades of the drive frequency. We find that these topological invariants not only remain nonzero at low frequencies, but increasingly fluctuate. While at any fixed frequency the invariants are deterministic, over a range of frequencies $\delta\Omega\ll\Omega$, the invariants distribute pseudorandomly. We argue that this distribution is universal and in our models is given by a Gaussian, whose width is $\sigma(\Omega)\sim 1/\sqrt{\Omega}$. We explain this universal behavior by revealing a diffusive process in the evaluation of the invariants and confirm our results numerically.

Specifically, we study one-dimensional driven systems with periodic boundary conditions, with a Hamiltonian of the form $\hat{H}=\int\hat{c}^{\dagger}_{k}h^{\vphantom{\dagger}}_{k}\hat{c}^{\vphantom{\dagger}}_{k}\frac{dk}{2\pi}$, where $k\in[-\pi,\pi]$ is the crystal momentum, $\hat{c}_{k}$ is a two-component spinor field, and $h_{k}=d_{kx}\sigma_{x}+d_{ky}\sigma_{y}$ with $d_{kx}+id_{ky}\equiv d_{k}$ a model-dependent function. For example, in the Su-Schrieffer-Heeger (SSH) model su1979 ; supp $d_{k}=2e^{ik/2}\left(w\cos\frac{k}{2}+i\delta\sin\frac{k}{2}\right),$ where $w$ ($\delta$) is the hopping (modulation) amplitude. In the Kitaev model Kit01a ; supp , after a suitable rotation in the Nambu space, one finds $d_{k}=2w\cos k-\mu+i\Delta\sin k$, where $\mu$ is the chemical potential and $\Delta$ is the nearest-neighbor pairing amplitude.

These Hamiltonians are particle-hole symmetric, $\sigma_{z}h_{-k}^{*}\sigma_{z}=-h_{k}$, with eigenvalues $\pm|d_{k}|$. In equilibrium, there are two topologically distinct phases: a topological phase, for $\delta/w>0$ in the SSH and $|\mu|<2|w|$ for Kitaev model, and a trivial phase otherwise. These two phases are distinguished on the lattice with open boundary conditions by the presence of zero-energy bound states in the topological phase. With periodic boundary conditions, the phases are distinguished by an integer topological invariant $\nu=0$ or $1$, equal to the winding number

[TABLE]

For a multi-band system, e.g. the SSH-Kitaev WatEzaTan14a ; supp , $d_{k}$ is matrix-valued and the topological invariant is found by $\mathcal{W}[\det d]$.

When the system is periodically driven, the full dynamics is obtained by solving the Floquet-Schrödinger equation $[h_{k}(t)-i\partial_{t}]\phi^{\pm}_{k}(t)=\pm\epsilon_{k}\phi^{\pm}_{k}(t)$ (we are setting $\hbar=1$) for the periodic steady states $\phi^{s}_{k}(t)=\phi^{s}_{k}(t+2\pi/\Omega)$, $s=\pm$, with the quasienergy $s\epsilon_{k}$, which we take to be in the Floquet zone $[-\Omega/2,\Omega/2]$. The Bloch evolution operator can then be written as $U_{k}(t)=\sum_{s=\pm}e^{-is\epsilon_{k}t}\phi_{k}^{s}(t)\phi_{k}^{s\dagger}(0)$. The full-period evolution operator $U_{k}(2\pi/\Omega)$ has eigenstates $\phi^{s}_{k}(0)$ with eigenvalues $e^{-2s\pi i\epsilon_{k}/\Omega}$. Since the quasienergy is a modular quantity, even a two-band model is characterized by two gaps at Floquet zone center (${\epsilon}^{\natural}=0$) and Floquet zone edge (${\epsilon}^{\flat}=\Omega/2$) Note3 . Thus, for periodic boundary conditions there are two independent topological invariants defined for the quasienergy gaps at Floquet zone center, ${\nu}^{\natural}$, and edge, ${\nu}^{\flat}$. For open boundary conditions, the corresponding invariants are the number of midgap steady bound states at Floquet zone center and edge supp .

To simplify our discussion, we take the drive protocol to satisfy the chiral reflection symmetry, $\delta(t+\pi/\Omega)=\delta(-t+\pi/\Omega)$; then, the two topological invariants are found Asbóth et al. (2014); supp from the half-period evolution operator U_{k}(\pi/\Omega)\equiv\left(\begin{array}[]{cc}A_{k}&B_{k}\\ C_{k}&D_{k}\end{array}\right), as

[TABLE]

In the static case, $D_{k}$ is constant and $B_{k}\propto d_{k}$, thus one finds ${\nu}^{\flat}=0$ and ${\nu}^{\natural}=\nu$ as expected. For concreteness, we present our results for the SSH model in the following and for other models in the Supplemental Material supp .

At symmetry points $k_{s}=0,\pm\pi$, $h_{k_{s}}\propto\sigma_{x}$ and the half-period evolution operator takes simple forms, $U_{0}(\pi/\Omega)=e^{-i\frac{2\pi}{\Omega}w\sigma_{x}}$ and $U_{\pi}(\pi/\Omega)=e^{i\frac{2\pi}{\Omega}\bar{\delta}\sigma_{x}},$ where $\bar{\delta}=(\Omega/2\pi)\int_{0}^{2\pi/\Omega}\delta(s)ds$ is the average hopping modulation through one drive cycle. The values $D_{\pm\pi}=\cos(2\pi\bar{\delta}/\Omega)$, $D_{0}=\cos(2\pi w/\Omega)$ and $B_{\pm\pi}=i\sin(2\pi\bar{\delta}/\Omega)$, $B_{0}=-i\sin(2\pi w/\Omega)$ can be used to anchor their winding.

To understand the changes in the winding number ${\nu}^{\flat}$ (${\nu}^{\natural}$) we analyze the contour of $D_{k}$ ($B_{k}$) in the complex plane as frequency varies (see Fig. 1). At high enough frequency the contour of $D_{k}$ ($B_{k}$) is a loop with two crossing points on the real (imaginary) axis at $D_{\pm\pi}$ and $D_{0}$ ($B_{\pm\pi}$ and $B_{0}$); as frequency is lowered the loop twists and untwists, thus changing the number of crossing points on the real (imaginary) axis via two processes: a pair of crossings are “emitted” from $D_{0}$ ($B_{0}$) whenever $D^{\prime}_{0}=0$ ($B^{\prime}_{0}=0$), where the prime denotes $\partial/\partial k$; on the other hand, a pair of crossings are “absorbed” into $D_{\pm\pi}$ ($B_{\pm\pi}$) when $D^{\prime}_{\pm\pi}=0$ ($B^{\prime}_{\pm\pi}=0$). While the rates of these processes depend on the drive protocol, they all scale with $1/\Omega$; thus the number of crossings generically grows as $1/\Omega$. As $\Omega$ is lowered, all crossings move back and forth within the unit disk along the real (imaginary) axis at a speed that scales with $1/\Omega$. When a crossing point of $D_{k}$ ($B_{k}$) passes through the origin, the winding ${\nu}^{\flat}$ (${\nu}^{\natural}$) changes. The inversion symmetry of the SSH model ensures that except $D_{\pm\pi}$ and $D_{0}$ ($B_{\pm\pi}$ and $B_{0}$), all other crossings are doubled.

Denoting the momenta at crossing points with ${k}^{\circ}_{c}(\Omega)$, where $\circ={\natural,\flat}$, the total number of crossings is ${N}^{\circ}={N}^{\circ}_{+}+{N}^{\circ}_{-}$, where ${N}^{\flat}_{\pm}=\sum_{c}\Theta(\pm D_{{k}^{\flat}_{c}}\,\mathrm{Im}\,D^{\prime}_{{k}^{\flat}_{c}})$ and ${N}^{\natural}_{\pm}=\sum_{c}\Theta(\pm iB_{{k}^{\natural}_{c}}\,\mathrm{Re}\,B^{\prime}_{{k}^{\flat}_{c}})$. The winding numbers, on the other hand, are given by ${\nu}^{\circ}=\frac{1}{2}({N}^{\circ}_{+}-{N}^{\circ}_{-})$. At any given frequency, $\Omega$, the values of ${N}^{\circ}_{\pm}$ may be computed deterministically from the number of crossings emitted, absorbed, and moved on the corresponding real or imaginary axis. However, as $\Omega\to 0$, these numbers grow in an increasingly complex way; thus, over a frequency interval $\delta\Omega\ll\Omega$ the distribution of crossing points appears random. We posit that this distribution can be modeled by a universal stochastic process of emission, absorption, and motion of crossing points of $D_{k}$ ($B_{k}$) Note0 . In the low-frequency limit, our numerics show generically that ${N}^{\circ}_{\pm}$ are equally distributed. Taking this to be true, we may think of ${N}^{\circ}_{\pm}$ as the number of steps taken by a one-dimensional random walker in opposite directions, with $2{\nu}^{\circ}$ the distance from the starting point. Thus, winding numbers are diffusive variables with a protocol-dependent diffusion constant ${\mathcal{D}}^{\circ}=2\Omega\sqrt{\langle{{N}^{\circ}_{+}}^{2}\rangle}$. Here, $\langle\cdots\rangle$ stands for the average in the stochastic model or, equivalently, the average over the interval $\delta\Omega$. The winding numbers acquire a Gaussian distribution with a width ${\sigma}^{\circ}(\Omega)=\sqrt{{\mathcal{D}}^{\circ}/\Omega}$. This is our main result.

Changes in the winding number are concomitant with quasienergy gap closings. This is easy to see at symmetry points $k_{s}$, where, for our chirally symmetric protocols, the full-period evolution operator is the square of the half-period evolution operator. At these points, ${\nu}^{\flat}$ (${\nu}^{\natural}$) change by one when $D_{0}$ and $D_{\pm\pi}$ ($B_{0}$ and $B_{\pm\pi}$) vanish, respectively, at ${\Omega}^{\flat}_{0}=\frac{4w}{2m-1}$ and ${\Omega}^{\flat}_{\pi}=\frac{4\bar{\delta}}{2m-1}$ (${\Omega}^{\natural}_{0}=\frac{2w}{m}$ and ${\Omega}^{\natural}_{\pi}=\frac{2\bar{\delta}}{m}$) for integer $m$. Noting that quasienergies at symmetry points are given by $\epsilon_{0}\equiv\pm 2w\>\mathrm{mod}\>\Omega$ and $\epsilon_{\pi}\equiv\pm 2\bar{\delta}\>\mathrm{mod}\>\Omega$, it is easy to see they are equal to ${\epsilon}^{\flat}$ (${\epsilon}^{\natural}$) exactly at frequencies where ${\nu}^{\flat}$ (${\nu}^{\natural}$) changes. Of course, changes in ${\nu}^{\flat}$ (${\nu}^{\natural}$) are also caused at any frequency ${\Omega}^{\flat}_{*}$ (${\Omega}^{\natural}_{*}$) and non-symmetry momenta ${k}^{\flat}_{*}\equiv{k}^{\flat}_{c}({\Omega}^{\flat}_{*})$ [${k}^{\natural}_{*}\equiv{k}^{\natural}_{c}({\Omega}^{\natural}_{*})$], where $D_{{k}^{\flat}_{*}}$ ($B_{{k}^{\natural}_{*}}$) vanishes and the gap at ${\epsilon}^{\flat}$ (${\epsilon}^{\natural}$) closes. Due to inversion symmetry, the winding numbers at these gap closings change by two. We note that the frequencies ${\Omega}^{\flat}_{*}$ and ${\Omega}^{\natural}_{*}$ depend on the drive protocol.

To proceed quantitatively, we choose a periodic two-step drive protocol in the SSH model given by $\delta(t)=\delta_{1}$ for $0<t<2\pi p_{1}/\Omega$ and $\delta_{2}$ for $2\pi p_{1}/\Omega<t<2\pi/\Omega$. Here, $0<p_{1}<1$ is the dimensionless fraction of the period for the first step of the drive. This family of protocols simplifies the numerical calculations, and allows us to obtain both analytically exact and numerically reliable results over a wide range of frequencies. Note that the modulation is chiral symmetric. This is explicit if we take the origin of time to be at $\pi p_{1}/\Omega$. Calculating the full-period evolution operator, the quasienergies are given by

[TABLE]

where the average and difference bands $\bar{E}_{k}=p_{1}|d_{1k}|+p_{2}|d_{2k}|$, $\breve{E}_{k}=p_{1}|d_{1k}|-p_{2}|d_{2k}|$ with index $a=1,2$, indicating $\delta=\delta_{a}$. Here, $\theta_{k}$ is the angle between the complex variables $d_{1k}$ and $d_{2k}$. Without loss of generality, we assume $|\delta_{1}|>|\delta_{2}|$. Gap closings at ${\epsilon}^{\circ}$ for ${k}^{\circ}_{*}\neq 0,\pi$ are obtained when $\cos(2\pi\bar{E}_{{k}^{\circ}_{*}}/\Omega)=\cos(2\pi\breve{E}_{{k}^{\circ}_{*}}/\Omega)={(-1)}^{\circ}$, where ${(-1)}^{\flat}=-1$ and ${(-1)}^{\natural}=1$. This is a resonant condition leading to an implicit equation for $\Omega$, which we solve numerically. Furthermore, for $|\delta_{2}/\delta_{1}|p_{2}<p_{1}<p_{2}$, there exist ${k}^{\natural}_{*}$ where $\breve{E}_{{k}^{\natural}_{*}}=0$; at these points, the gap at ${\epsilon}^{\natural}$ closes for $\bar{E}_{{k}^{\natural}_{*}}/{\Omega}^{\natural}_{*}\equiv 0\>\mathrm{mod}\>1$.

The winding numbers are found from

[TABLE]

and $B_{k}=(d_{1k}/|d_{1k}|)\tilde{B}_{k}$,

[TABLE]

Since $\theta_{\pi}=0$ or $\pi$ and $\theta_{0}=0$, $D_{\pi}$ ($\tilde{B}_{\pi}$) and $D_{0}$ ($\tilde{B}_{0}$) are real. We note $\mathcal{W}[B]=\mathcal{W}[d_{1}]+\mathcal{W}[\tilde{B}]$, i.e. ${\nu}^{\natural}=\nu_{1}+\mathcal{W}[\tilde{B}]$.

Focusing on ${\nu}^{\flat}$ for concreteness and using the explicit forms of $D_{k}$, we find that the crossing points that contribute to either $N_{+}$ or $N_{-}$ are emitted when $2p_{1}w/\Omega$ or $2p_{2}w/\Omega$ is an integer, and they are absorbed when $p_{1}|\delta_{1}|/\Omega$ or $2p_{2}|\delta_{2}|/\Omega$ is an integer. For $p_{1}\approx p_{2}$ and for small enough frequency, we may assume the motion of crossing points yields a nearly uniform distribution along the real axis. Since the winding number varies by 2 only when the crossing two points are on different halves of the real axis, the diffusion constant may be obtained by ${\mathcal{D}}^{\flat}\approx 2|p_{1}(w-|\delta_{1}|)-p_{2}(w-|\delta_{2}|)|$. In the following, we assume $\delta_{1}<0$ for simplicity.

Our analytical expressions for the two-step drive allow for the exact determination of gap closings; however, in general, quasienergy gaps and topological invariants can only be obtained by numerical solutions. In the special case of a symmetric drive, $\delta_{2}=-\delta_{1}$ and $p_{1}=p_{2}$, we can calculate the topological invariants exactly: ${\nu}^{\natural}=0$, and ${\nu}^{\flat}=-1$ when $\lfloor\frac{2w}{\Omega}-\frac{1}{2}\rfloor$ is even and 0 otherwise.

For the numerical solutions, we consider three distinct protocols: the asymmetric protocol, $\delta_{2}=-\delta_{1}-\varrho>0$, which is a periodic switch between the equilibrium trivial and the topological phases with the asymmetry parameter $\varrho$; the critical protocol, $\delta_{2}=0$, i.e. a periodic switch between the equilibrium trivial and the critical point of the system; and finally, the trivial protocol, $\delta_{1}<\delta_{2}<0$, such that the systems is in the equilibrium trivial phase at all times. Our numerical results for ${\nu}^{\flat}$ are summarized in Fig. 2; our results for ${\nu}^{\natural}$ are similar supp .

The quasienergy gaps ${\Delta}^{\circ}=\min_{k}|\epsilon_{k}-{\epsilon}^{\circ}|$ exhibit self-similar patterns, with peaks that scale as $\sim\Omega^{2}$. We have benchmarked our numerical calculation with the exact analytical expressions for gap closings, shown on the same plot; the agreement is extremely good. The $\Omega^{2}$ scaling can be understood within adiabatic perturbation theory RigOrtPon08a ; RodLenSer18a , where the frequency is used as a perturbation parameter. The first-order correction to the quasienergy is the Berry phase of the steady states, which vanishes for our chirally symmetric protocols RodLenSer18a ; MarMoi15a .

The winding number ${\nu}^{\flat}$ fluctuates as frequency is lowered with an increasing relative amplitude. For an asymmetric protocol, as in Fig. 2(c) center panel, when $\Omega>2\varrho$, we observe the regular step-wise behavior as in the symmetric protocol. However, when $\Omega<2\varrho$, the same fluctuating pattern sets in. We have carried out a detailed analysis of ${\nu}^{\flat}$ over a wide range of low frequencies. For each frequency $\Omega$, we show ${N}^{\flat}$, ${N}^{\flat}_{+}$, ${\alpha}^{\flat}_{+}={N}^{\flat}_{+}/{N}^{\flat}$, and ${\nu}^{\flat}$. Both ${N}^{\flat}$, and ${N}^{\flat}_{+}$ scale linearly with $1/\Omega$, confirming our general arguments. The ratio ${\alpha}^{\flat}_{+}$ approaches $1/2$ as the frequency decreases, indicating the diffusive behavior familiar from a random-walk process. Moreover, the range of $|{\nu}^{\flat}|$ scales as $\sim 1/\sqrt{\Omega}$. For a range of frequencies much lower than other energy scales in the system, we have determine the probability $P({\nu}^{\flat})$ of finding ${\nu}^{\flat}$ in our numerical histogram. As shown, $P({\nu}^{\flat})$ follows a Gaussian distribution with a width that is given by $\sigma({\nu}^{\flat})=\sqrt{\mathcal{D}^{\flat}/\Omega}$ over several decades of frequency, confirming our general result.

To test our general arguments for the universality of the fluctuations, we have studied other drive protocols and other models, including a model with multiple bands. These numerical studies support our results in all cases. Details are found in the Supplemental Material supp ; here, we present our results for a multi-step protocol in the SSH model approximating $\delta(t)=\bar{\delta}+A\cos(\Omega t)$. The analytical calculations become increasingly difficult as the number of steps $n$ in the drive increases; however, we can still calculate the topological invariants numerically. A typical sampling of our results are shown in Fig. 7. While ${\nu}^{\flat}$ fluctuates both in magnitude and sign as $n$ is varied, the ratio ${\alpha}^{\flat}_{+}\approx 1/2$, again indicating a diffusive process. Collecting good statistics over a wide frequency range quickly becomes too expensive. However, since the fluctuations in the winding number result from the twisting and untwisting of the contour $D_{k}$, we expect that varying the number of steps $n$ should have a similar effect. Indeed, as shown in Fig. 7(b), after averaging over $n$, the root-mean-square $\delta{\nu}^{\flat}\sim 1/\sqrt{\Omega}$.

In conclusion, we have found universal fluctuations in the topological invariants characterizing a Floquet topological phase. We explained these fluctuations by positing a pseudorandom distribution of crossing points of the complex function whose winding number gives the topological invariant. This distribution follows from the diffusive process of emission, absorption, and motion of crossing points as frequency is lowered. Our results show that the limit $\Omega\to 0$ has a rich structure that is distinct from the simple adiabatic limit: while the topological invariant vanishes Note2 on average, consistent with the adiabatic limit, its fluctuations diverge. These fluctuation may be observed in the noise spectra of relevant quantities such as voltage noise RodFerSer18a , or by spectroscopic measures of the number of Floquet edge modes as recently observed in a photonic crystal emulator ChePanWan18a .

Universal fluctuations in Chern numbers have been studied in quantized classically-chaotic and random matrix theories LebKurFei90a ; WalWil95a . By contrast, we study periodically driven systems, where topology is characterized not just by Chern numbers of a static Hamiltonian, but by independent winding numbers through a drive cycle. In this context, it would be interesting to study if driven systems with different symmetries (say, other than chiral symmetry) can support other universality classes of fluctuations of Floquet topological invariants.

I Supplemental Material

In this Supplemental Material, we recap, for completeness, the derivation of the Floquet topological invariants in one-dimensional driven systems with a chiral reflection symmetry, and present the details of our numerical studies for driven Su-Schrieffer-Heeger (SSH), Kitaev, and SSH-Kitaev models, as well as the multi-step drive protocol in driven SSH model.

I.1 Topological invariant for driven one-dimensional chiral-symmetric systems

In this section, we briefly outline the derivation of the topological invariant for one-dimensional chiral symmetric systems, following Ref. Asbóth et al. (2014). A driven system has chiral symmetry if there exist a unitary, hermitian, local operator $\Gamma=\Gamma^{\dagger}=\Gamma^{-1}$, such that

[TABLE]

where the time-ordered exponential $U(\tau,\tau+T)=\mathbb{T}\exp\left[-i\hbar\int_{\tau}^{\tau+T}H(s)ds\right]$ is the full-period evolution operator with the initial time $\tau$ and period $T=2\pi/\Omega$ for the Hamiltonian $H(t)=H(t+T)$. A topological invariant characterizing a chiral-symmetric one-dimensional system with perioidic boundary conditions can be written in the diagonal basis of $\Gamma$, where

[TABLE]

as

[TABLE]

where $k$ is the lattice momentum.

Now, we say the drive protocol in $H(t)$ has chiral reflection symmetry if for some $\tau$ and $\tau^{\prime}\in(0,T)$, one has $U(\tau+\tau^{\prime},\tau+T)=\Gamma F^{\dagger}\Gamma$, where $F=U(\tau,\tau+\tau^{\prime})$. That is, the drive protocol starting at $\tau$ has a reflection symmetry about $\tau+\tau^{\prime}$. Then, $U(\tau,\tau+T)=\Gamma F^{\dagger}\Gamma F$ and $U(\tau+\tau^{\prime},\tau+\tau^{\prime}+T)=F\Gamma F^{\dagger}\Gamma$ both satisfy satisfy Eq. (6).

Thus, we can define two topological invariants $\nu=\nu_{\tau}$ and $\nu^{\prime}=\nu_{\tau+\tau^{\prime}}$. Physically, they can be interpreted as a bulk “sublattice” polarization at times $\tau$, and $\tau+\tau^{\prime}$. Given that states with quasienergy ${\epsilon}^{\flat}=\Omega/2$ switch sublattice when they evolve from $\tau$ to $\tau+\tau^{\prime}$, neither $\nu$ nor $\nu^{\prime}$ alone are related to the number of edge states in an open system, instead

[TABLE]

where $\nu^{\natural}$ and $\nu^{\flat}$ are the number of edge states with quasienergy ${\epsilon}^{\natural}=0$, and ${\epsilon}^{\flat}=\Omega/2$, respectively. Using algebraic properties of winding numbers in the diagonal basis of $\Gamma$, the authors of Ref. Asbóth et al. (2014) showed that $\nu^{\natural}$ and $\nu^{\flat}$ can be derived from the diagonal and off-diagonal blocks of $F=\begin{bmatrix}A&B\\ C&D\end{bmatrix}$ as $\nu^{\natural}=\mathcal{W}[B]$, $\nu^{\flat}=\mathcal{W}[D]$.

I.2 Driven SSH model

The tight-binding Hamiltonian for driven the Su-Schrieffer-Heeger (SSH) model is given by the

[TABLE]

where $\hat{c}^{\dagger}_{r}$ creates a fermion at lattice site $r$, $w$ is the unmodulated hopping amplitude and $\delta(t)=\delta(t+2\pi/\Omega)$ is the hopping modulation, periodic in time with frequency $\Omega$. For periodic boundary conditions, the crystal momentum $k$ is a good quantum number; defining the spinor $\hat{c}^{\dagger}_{k}=\sum_{x}e^{ikx}(\hat{c}^{\dagger}_{2x},\hat{c}^{\dagger}_{2x+1})$, where $x=\lfloor{r/2}\rfloor$ indexes the unit cells, we have $\hat{H}=\int\hat{c}^{\dagger}_{k}h^{\vphantom{\dagger}}_{k}\hat{c}^{\vphantom{\dagger}}_{k}\frac{dk}{2\pi}$, with $k$ in the Brillouin zone $[-\pi,\pi]$, $h_{k}=d_{kx}\sigma_{x}+d_{ky}\sigma_{y}$ and $d_{kx}+id_{ky}\equiv d_{k}=2e^{ik/2}\left(w\cos\frac{k}{2}+i\delta\sin\frac{k}{2}\right).$ Discrete symmetries of this Hamiltonian include inversion $\sigma_{x}h_{k}\sigma_{x}=h_{-k}$, sublattice $\sigma_{z}h_{k}\sigma_{z}=-h_{k}$, and particle-hole symmetry $\sigma_{y}h_{k}^{*}\sigma_{y}=-h_{k}$, which place the static SSH model in the BDI class Chiu et al. (2016). The instantaneous eigenvalues are $\pm|d_{k}|$ with $|d_{k}|=2\sqrt{w^{2}\cos^{2}\frac{k}{2}+\delta^{2}\sin^{2}\frac{k}{2}}$.

I.2.1 Quasienergy spectrum

In this section, we consider a driven SSH model with open boundary conditions, and calculate the number of edge states present in the system at given frequency. Fig. 4 shows the quasienergies as a function of the inverse frequency in the range $\Omega/w\in[0.1,1]$. The system has $N=1000$ sites and is driven by a two-step protocol switching between $\delta_{1}/w=-0.75$ and $\delta_{2}/w=-0.25$. States with quasienergies ${\epsilon}^{\natural}=0$, and ${\epsilon}^{\flat}=\Omega/2$ are localized at the edge of the system.

Top panel of Fig. 5(a) [(b)] shows $n^{\flat}$ ($n^{\natural}$), the number of edge mode pairs with quasienergies ${\epsilon}^{\flat}$ (${\epsilon}^{\natural}$), and are compared with the corresponding topological invariant $\nu^{\flat}$ ($\nu^{\natural}$). We use the same parameters as in Fig. 4. In order to correctly resolve the quasienergies, we consider a systems with $N=1000$ sites for the range $w/\Omega\in[0.3,1]$, $N=2000$ for $w/\Omega\in[1,3]$, and $N=3000$ for $w/\Omega\in[3,3.38]$.

Both quantities are in agreement, as expected from the bulk-boundary correspondence. Deviations are observed at frequencies in proximity to gap closings, where a larger system is required to correctly resolve the spectrum. Bottom panel of Fig. 5(a) [(b)] shows a close-up in the quasienergy spectrum around ${\epsilon}^{\flat}$ (${\epsilon}^{\natural}$). Changes in $\nu^{\flat}$ ($\nu^{\natural}$) occur only when the gap around ${\epsilon}^{\flat}$ (${\epsilon}^{\natural}$) closes.

I.2.2 Topological invariant $\nu^{\natural}$

In this section, we discuss the low-frequency limit of $\nu^{\natural}$, the topological invariant related to the number of edge states with quasienergy ${\epsilon}^{\natural}$. As discussed in the main text, $\nu^{\natural}$ also presents fluctuations in the low-frequency regime, with exceptions in special cases. For example, in the case of a symmetric drive, with $\delta_{2}=-\delta_{1}$ and $p_{1}=p_{2}$ we obtain ${\nu}^{\natural}=0$.

In general, we have to evaluate the invariant ${\nu}^{\natural}$ numerically. Fig. 6 summarizes our results. The winding number ${\nu}^{\natural}$ fluctuates as frequency is lowered with a relative amplitude that grows as $\Omega\to 0$. Both ${N}^{\natural}$ and ${N}^{\natural}_{+}$ scale linearly with $1/\Omega$. The ratio ${\alpha}^{\natural}_{+}={N}^{\natural}_{+}/{N}^{\natural}$ approaches $1/2$ as the frequency decreases. Moreover, the range of $|{\nu}^{\natural}|$ scales as $\sim 1/\sqrt{\Omega}$. All of these properties are similar to the low-frequency behavior of ${\nu}^{\flat}$. For a range of frequencies much lower than other energy scales in the system, we have determine the probability $P({\nu}^{\natural})$ of finding ${\nu}^{\natural}$ in our numerical histogram. As shown, $P({\nu}^{\natural})$ follows a Gaussian distribution with a width that is given by $\sigma({\nu}^{\natural})=\sqrt{\mathcal{D}^{\natural}/\Omega}$ over a decade of frequencies. When the drive protocol is entirely in the static topological (trivial) phase, the probability $P({\nu}^{\natural})$ is centered at ${\nu}^{\natural}=1$ (${\nu}^{\natural}=0$).

I.2.3 Multi-step drive protocol

In this section, we consider a multi-step driving protocol as an approximation to the smooth protocol $\delta(t)=\bar{\delta}+A\cos(\Omega t)$. For a protocol with $n$ steps, the half-propagator is given by

[TABLE]

In the limit $n\rightarrow\infty$ we recover the smooth protocol. In Fig. 7(a) we plot the invariant $|{\nu}^{\flat}|$ and the ratio ${\alpha}^{\flat}_{+}$ as a function of the frequency $\Omega$ for a protocol with $n=20$ steps. As in the case of a two-step drive, we find that the ration ${\alpha}^{\flat}_{+}$ approaches $1/2$ as the frequency decreases, suggesting a diffusive process. Also, $|{\nu}^{\flat}|$ presents fluctuations that grow as $\sim 1/\sqrt{\Omega}$. For large $n$, it becomes difficult to accurately calculate the invariant, and collecting enough statistics as a function of frequency to determine the scaling in the low-frequency limit becomes more challenging. However, varying $n$ changes the structure of complex function $D_{k}$, creating new twists. Therefore, we compute ${\nu}^{\flat}$ for different frequencies $\Omega$ as a function of $n$ in the range $[2,50]$ and consider the average over $n$. In Fig. 7(b), we plot the root-mean-square $\delta{\nu}^{\flat}$ as a function of frequency and find that it also scales as $1/\sqrt{\Omega}$.

I.3 Driven Kitaev Model

In this section we consider the driven Kitaev model, and study the behaviour of the topological invariants in the low-frequency regime. The Hamiltonian is given by Kitaev (2001)

[TABLE]

where $w$ and $\Delta$ are the nearest-neighbor hopping and pairing amplitudes respectively, and $\mu$ is the chemical potential. Imposing periodic boundary conditions, and introducing the Nambu spinors $\hat{\psi}_{k}^{\raisebox{-1.0pt}{{\scriptsize{T}}}}=\left(\hat{c}^{\vphantom{\dagger}}_{k},\hat{c}^{\dagger}_{-k}\right)$, the Hamiltonian can be written as $\hat{H}=\frac{1}{2}\sum_{k}\hat{\psi}^{\dagger}_{k}h^{\vphantom{\dagger}}_{k}\hat{\psi}^{\vphantom{\dagger}}_{k}-\mu N$, where $h_{k}=\left(2w\cos k-\mu\right)\sigma_{z}+2|\Delta|\sin k\left(\cos\phi\>\sigma_{y}-\sin\phi\>\sigma_{x}\right).$ This has the particle-hole symmetry $h_{k}=-\sigma_{x}h_{-k}^{*}\sigma_{x}$. Upon the rotation $R=\exp(-i\sigma_{y}\pi/4)\exp(i\phi\sigma_{z}/2)$, the Hamiltonian is transformed to

[TABLE]

as shown in the main text. In this representation, the Kitaev model has the discrete symmetries $h_{k}=-\sigma_{z}h^{*}_{-k}\sigma_{z}$, and $h_{k}=h^{*}_{-k}$. The instantaneous eigenvalues are $\pm\sqrt{\left(2w\cos k-\mu\right)^{2}+|\Delta|^{2}\sin^{2}k}$.

In a system with open boundary conditions we find one pair of Majorana states at zero energy for $|\mu|<2|w|$, and zero for $|\mu|>2|w|$. At $|\mu|=2w$ the gap closes. The static topological invariant gives $\nu=1$ for $|\mu|<2w$ and zero otherwise.

Choosing a suitable driving protocol, the driven Kitaev model also has chiral symmetry, and the topological invariants are given by ${\nu}^{\flat}$ and ${\nu}^{\natural}$, as for the SSH model. In this work, we choose to drive the chemical potential $\mu(t)$ with the two-step protocol $\mu(t)=\mu_{1}$ for $0<t<2\pi p_{1}/\Omega$ and $\mu_{2}$ for $2\pi p_{1}/\Omega<t<2\pi/\Omega$, where $0<p_{1}<1$. We fix $p_{1}=0.5$. The bottom panels in figure 8 show the quasienergy spectrum for a chain with $N=1000$ sites in the frequency range $\Omega/w\in[0.5,1.8]$ around the Floquet zone edge and center. Top panel (a) [(b)] shows the topological invariant ${\nu}^{\flat}$ (${\nu}^{\natural}$) underneath the number of states with quasienergy ${\epsilon}^{\flat}$ (${\epsilon}^{\natural}$), shown in red. These quantities are in agreement. Differences can be found at frequencies where the corresponding quasienergy gap is closed.

We now evaluate the topological invariants numerically in the low-frequency regime. For concreteness, we focus on ${\nu}^{\flat}$. Our results are summarized in Fig. 9. As for the SSH model, we find that the range of $|{\nu}^{\flat}|$ scales as $\sim 1/\sqrt{\Omega}$, both ${N}^{\flat}$ and ${N}^{\flat}_{+}$ scale linearly with $1/\Omega$, and the ratio ${\alpha}^{\natural}_{+}$ approaches $1/2$ as the frequency decreases. Furthermore, the probability $P({\nu}^{\flat})$ follows a Gaussian distribution with width $\sigma({\nu}^{\flat})=\sqrt{\mathcal{D}^{\flat}/\Omega}$.

I.4 Four-band driven Kitaev-SSH Model

In this section we consider the Kitaev-SSH model Wakatsuki et al. (2014). The Hamiltonian is given by

[TABLE]

where $x$ labels the unit cells, and $A$, $B$ label the sublattices in each unit cell, $w+\delta$ ($w-\delta$) and $\Delta$ are the intra-cell (inter-cell) hopping and pairing amplitudes respectively and $\mu$ is the chemical potential. If we consider periodic boundary conditions, and define the spinor $\hat{\psi}^{\dagger}_{k}=\left(\hat{c}^{\dagger}_{kA},\hat{c}^{\dagger}_{kB},\hat{c}^{\vphantom{\dagger}}_{-kA},\hat{c}^{\vphantom{\dagger}}_{-kB}\right)$, the momentum-dependent Hamiltonian takes the form Wakatsuki et al. (2014)

[TABLE]

where $z_{k}=-\left[(w+\delta)+(w-\delta)e^{-ik}\right]$, and $\Delta_{k}=-\Delta\left[(1+\delta/w)-(1-\delta/w)e^{-ik}\right]$. For our purpose, the most important feature of this model is that it has four energy bands, given by $\pm\sqrt{\mu^{2}+|z_{k}|^{2}+|\Delta_{k}|^{2}\pm 2\sqrt{\mu^{2}|z_{k}|^{2}+(4\Delta\delta)^{2}}}$

For $\delta/w<-0.2$ we have two pairs of zero-energy states localized at the edge, while for $-0.2<\delta/w<0.2$ we have only one pair. For $\delta/w>0.2$ there are no zero-energy states. For simplicity, now on we consider the case $\mu_{1}/w=\mu_{2}/w=0$, although is possible to also consider finite chemical potentials. The chiral basis for this model is defined by the transformation $\Gamma=\frac{1}{2}(\tau_{0}\otimes\sigma_{0}+\tau_{\alpha}\otimes\sigma_{\alpha}),$ where $\sigma_{\alpha}$ and $\tau_{\alpha}$ for $\alpha=x,y,z$ are the Pauli matrices in particle-hole and sublattice space, and $\tau_{0}=\sigma_{0}=\mathbb{1}$.

We choose the driving protocol $\delta(t)=\delta_{1}$ for $0<t<2\pi p_{1}/\Omega$ and $\delta_{2}$ for $2\pi p_{1}/\Omega<t<2\pi/\Omega$, where $0<p_{1}<1$ to ensure that the driven system posses chiral symmetry. As before, we will fix $p_{1}=0.5$. In Fig. 10, the bottom panels show the quasienergy spectrum for a driven chain with $N=1000$ sites around the Floquet zone edge (a) and center (b). The rest of the parameters used are $\delta_{1}/w=3.0$, $\delta_{2}/w=-2.0$, and $|\Delta|/w=0.2$. Top panel (a) [(b)] shows the topological invariant ${\nu}^{\flat}$ (${\nu}^{\natural}$) underneath of ${n}^{\flat}$ (${n}^{\natural}$), the number of states with quasienergy ${\epsilon}^{\flat}=\Omega/2$ (${\epsilon}^{\natural}=0$), shown in red. These quantities are in agreement, as expected from the bulk-edge correspondence.

The topological invariants are derived from the half-period evolution propagator U_{k}(\pi/\Omega)\equiv\left(\begin{array}[]{cc}A_{k}&B_{k}\\ C_{k}&D_{k}\end{array}\right), where $A_{k},B_{k},C_{k}$ and $D_{k}$ are now $2\times 2$ matrices. The winding numbers are given by $\nu^{\flat}=\mathcal{W}[D]$, and $\nu^{\natural}=\mathcal{W}[B]$ where $\mathcal{W}[h]=\frac{1}{2\pi i}\int_{-\pi}^{\pi}\frac{\partial}{\partial k}\ln\left(\mbox{det}\;h_{k}\right)dk.$ As explained in the main text, for multi-band systems, the invariant is still defined by the winding of a complex functions, $\det D_{k}$ and $\det B_{k}$. The evaluation of $\nu^{\flat}$ is done numerically, and our results are summarized in Fig. 11. As in the other two driven chiral-symmetric systems, the SSH model and the Kitaev model, we find that ${\nu}^{\flat}$ fluctuates as frequency is lowered. The probability $P({\nu}^{\flat})$ also follows a Gaussian distribution, with a width given by $\sigma({\nu}^{\flat})=\sqrt{\mathcal{D}^{\flat}/\Omega}$.