# Family-Vicsek dynamical scaling and Kardar-Parisi-Zhang-like superdiffusive
growth of surface roughness in a driven one-dimensional quasiperiodic model

Sreemayee Aditya1, Nilanjan Roy2,3 1Centre for High Energy Physics, Indian
Institute of Science, Bengaluru 560012, India
2Centre for Condensed Matter Theory, Department of Physics, Indian Institute
of Science, Bengaluru 560012, India
3 Division of Physics and Applied Physics, School of Physical and Mathematical
Sciences, Nanyang Technological University, Singapore 637371

###### Abstract

The investigation of the dynamical universality classes of quantum systems is
an important, and rather less explored aspect of non-equilibrium physics. In
this work, considering the out-of-equilibrium dynamics of spinless fermions in
a one-dimensional quasiperiodic model with and without periodic driving, we
report the existence of the dynamical one-parameter Family-Vicsek (FV) scaling
of the “quantum surface-roughness” associated with the particle-number
fluctuations. In the absence of periodic driving, the model is interestingly
shown to host a subdiffusive critical phase with anomalous FV scaling
exponents, separated by two subdiffusive critical lines and a triple point
from other phases. Our analysis on the fate of the critical phase in the
presence of (inter-phase) driving indicates that the critical phase is quite
fragile, and has a tendency to get absorbed into the delocalized or localized
regimes depending on the driving parameters, especially in the slow driving
limit. Interestingly, periodic driving can also conspire to show quantum
Kardar-Parisi-Zhang (KPZ)-like superdiffusive dynamical behavior, which seems
to have no classical counterpart. We further construct an effective Floquet
Hamiltonian, which qualitatively captures this feature occurring in the driven
model.

##  I Introduction

Although the out-of-equilibrium dynamics of closed quantum systems have been
studied quite extensively in the past decade Polkovnikov  _et al._ (2011);
Langen  _et al._ (2015); Oka and Kitamura (2019), the dynamical scaling has
not been addressed enough until very recently. In a series of recent
theoretical works, a dynamical one-parameter scaling, namely, the Family-
Vicsek (FV) scaling is found in a variety of systems, namely, the clean
interacting Bose-Hubbard model Fujimoto  _et al._ (2020), free fermionic
models in presence of random disorder Fujimoto  _et al._ (2021) and
dissipation Fujimoto  _et al._ (2022) by studying the surface roughness
operator (see Eq. 2). In analogy to the surface growth operator in classical
systems, where FV scaling is studied quite extensively Vicsek and Family
(1984); Family and Vicsek (1985); Barabási and Stanley (1995), a similar
“quantum surface-height operator” is introduced in this context Jin  _et al._
(2020), which essentially represents the fluctuations of particle-number over
a given length scale. The measure of the standard deviation of this operator
defines the quantum surface roughness, which plays the key role in giving rise
to the one-parameter FV scaling relation. The surface growth physics in
classical systems primarily demonstrates two dynamical universality classes -
the diffusive or Edwards-Wilkinson (EW) class Edwards and Wilkinson (1982) and
the superdiffusive or Kardar-Parisi-Zhang (KPZ) universality class Kardar  _et
al._ (1986) with dynamical exponent $z$ = 2 and 3/2, respectively. Also an
anomalous diffusion have been recently reported across the Ising and
Kosterlitz-Thouless transitions in classical 2D XXZ model Ruidas and Banerjee
(2021). However, in quantum systems, the most commonly found non-trivial
behavior is either diffusive or ballistic. Nevertheless, KPZ-like dynamical
behavior of different-time correlations has been found recently in interacting
integrable spin chains in presence of the non-abelian symmetries both
theoretically Jin  _et al._ (2020); Ljubotina  _et al._ (2019); Ilievski  _et
al._ (2018); Ye  _et al._ (2022); Keenan  _et al._ (2022) and experimentally
Wei  _et al._ (2022). As opposed to such different-time correlation functions,
quantum surface roughness is a sum of equal time correlations, and can be
defined by a local quantity; further, it directly demonstrates FV scaling
Fujimoto  _et al._ (2020). Moreover, it is an easily accessible quantity using
a microscope in cold atomic experiments Fujimoto  _et al._ (2021), and
provides a direct way of measuring the entanglement entropy in a non-
interacting many-body quantum systems Fujimoto  _et al._ (2021).

Further, out-of-equilibrium dynamics in the presence of quasiperiodic disorder
has been a current subject to intense research Roy and Sharma (2021); Ahmed
_et al._ (2021); Purkayastha  _et al._ (2017); Varma  _et al._ (2017); Iglói
_et al._ (2013). Some of the recent works also look into the interplay between
periodic driving and quasiperiodic disorder Sarkar  _et al._ (2021); Aditya
_et al._ (2023); Roy  _et al._ (2018); Shimasaki  _et al._ (2022); Ray  _et
al._ (2018), which give rise to rich phase diagrams with no equilibrium
counterparts. Moreover, the recent advances in cold atom experiments enable us
to uncover these purely non-equilibrium phases tailored by Floquet engineering
in experimental settings Schreiber  _et al._ (2015); Eckardt (2017); Messer
_et al._ (2018). Considering the intricacies of out-of-equilibrium dynamics in
the quasiperiodic models, one can ask the following questions pertaining to
the surface roughness dynamics: Can we find dynamical universality classes
other than ballistic and diffusive ones in a simple quasiperiodic non-
interacting quantum model? How good does the FV scaling work for the critical
phase having only critical states that lead to critical slowing down of
dynamics Chaikin and Lubensky (1995)? Can we capture the dynamics emerging out
of the interplay between quasiperiodic disorder and periodic driving in a one-
dimensional closed quantum system within the framework of FV scaling? Is it
possible to achieve the emergent KPZ-like dynamics by tuning the parameters
appropriately in the driven system?

In this work, we study surface roughness growth in a one-dimensional model of
spinless fermions with a quasiperiodic potential and a quasiperiodic hopping
in non-equilibrium settings with and without a periodic square pulse protocol.
The single-particle phase diagram of the static model, as shown in Fig. 1,
comprises three phases: extended, critical, and localized along with three
critical lines, which meet at a triple point Han  _et al._ (1994). We
investigate the presence of FV scaling, particularly, on the critical lines
and in the critical phase, where the behaviors are expected to be the most
intriguing. We find that while the vertical critical line ($\mu=2$) containing
the Aubry-André critical point ($\lambda=0$, $\mu=2$) demonstrates (almost)
diffusive growth dynamics, the other two critical lines ($\lambda=1$ and
$\mu=2\lambda$, respectively) and the critical phase exhibit subdiffusive
dynamics with anomalous scaling exponents. A scaling analysis for the inverse
participation ratio (IPR) Aoki (1986); Mirlin  _et al._ (2006); Janssen
(1998); Evers and Mirlin (2008a) of the single-particle Floquet eigenstates of
the driven system reveals that the (inter-phase) drive-induced mixing between
critical and delocalized (localized) phases impels the system toward the
delocalization (localization) in the low driving frequency limit. Furthermore,
the phase-mixing between delocalized and localized phases shows a re-entrant
delocalization transition Aditya  _et al._ (2023); Ray  _et al._ (2018);
Shimasaki  _et al._ (2022) at extremely low driving frequencies. In the
stroboscopic dynamics of surface roughness at intermediate driving
frequencies, we observe an emergent KPZ-like superdiffusive dynamics, which
are primarily governed by the quantum nature Jin  _et al._ (2020); Ljubotina
_et al._ (2019); Ilievski  _et al._ (2018); Ye  _et al._ (2022); Keenan  _et
al._ (2022) of the model with the states evolving according to the quantum
Schrödinger equation as well as the quantum nature of the constituent
particles (fermions in our case). We summarize our results in Fig. 1. Finally,
we construct an effective Hamiltonian using the van-Vleck perturbation theory
Mikami  _et al._ (2016) that describes the KPZ dynamics of the driven model.

Figure 1: Schematic of the main results: Top panel: Phase diagram of our
static system described in Eq. 1 along with the direction of mixing when the
system is subjected to periodic driving. Bottom panel: A table showing the
effect of drive-induced phase mixing in the intermediate to low-frequency
regimes.

The rest of the paper is organized as follows. In Sec. II, we discuss the
static model that we consider throughout our work to study the dynamics in the
absence and presence of periodic driving. In Sec. III, we discuss the
dynamical universality classes of the critical lines and critical phase in the
framework of Family-Vicsek (FV) scaling while studying the exact surface
roughness dynamics in Sec. III.1 and approximate surface roughness dynamics in
relation to dynamics of half-chain entanglement entropy in Sec. III.2 starting
from a many-body initial state. In addition to that, we also discuss the
single-particle wave packet dynamics (charge transport) for the critical lines
and critical phase in Sec. III.3. This is followed by Sec. IV where we examine
the effect of inter-phase driving across the critical lines of the static
model. Next in Sec. V, we discuss the dynamical behaviors of the critical
phases subjected to periodic driving, which involves Kardar-Parisi-Zhang
(KPZ)-like superdiffusive dynamics as detailed in Sec. V.1. We further explain
this emergent behavior in by employing van Vleck expansion and Floquet
perturbation theory in the high-frequency and large driving amplitude regimes
in Sec. V.2 and Sec. V.3, respectively. This is followed by Sec. VI, where we
discuss the re-scaled surface height distribution for different dynamical
phases appearing in the static model as well as in the driven model. Finally,
in Sec. VII, we conclude by discussing the main results and point out possible
experiments that can confirm our theoretical findings.

##  II The static model

We consider a one-dimensional spinless fermionic model with correlated
disorders, which is described by the following Hamiltonian

$\displaystyle H$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{L}\left[J+\lambda\cos\left(2\pi\beta\left(j+1/2\right)+\phi\right)\right]\left(c_{j}^{\dagger}c_{j+1}+\rm{H.c.}\right)$
(1) $\displaystyle+\mu\sum_{j=1}^{L}\cos(2\pi\beta j+\phi)n_{j},$

where $J$ and $\lambda$ denote the uniform and quasiperiodic nearest-neighbor
hopping, respectively, $\mu$ stands for the quasiperiodic potential, $\phi$
denotes a global phase lying between $0$ and $2\pi$, and
$\beta=(\sqrt{5}-1)/2$ is the irrational golden mean. We further set $J$ to be
$1$ for the rest of our analysis. The phase diagram of this static model is
already well-studied in the literature Han  _et al._ (1994), which consists of
three phases (extended, critical and localized) and three critical lines
separating the phases (see Fig. 1). Interestingly, this model supports no
mobility edge. For our future convenience, we use the following terminology
while referring to the three critical lines: vertical ($\mu=2$), horizontal
($\lambda=1$), and tilted ($\mu=2\lambda$), respectively, as shown in Fig. 1.
In addition to that, these three critical lines meet at a triple point as also
indicated in Fig. 1. The other end of the vertical critical line is nothing
but the Aubry-André (AAH) critical point Aubry and André (1980), which is also
a self-dual point as already well-known in the literature. Moreover, the
vertical critical line also falls under the universality class of critical
exponents same as the self-dual point Sokoloff (1985).

##  III The surface roughness and Family-Vicsek scaling

###  III.1 Exact surface roughness dynamics

For dynamical characterization of these phases appearing in the static model,
we employ the Family-Vicsek (FV) scaling of surface roughness while
specifically focusing on the critical lines and critical phase. The main
motivation behind the proposal of “quantum surface roughness” is the
correspondence between $\partial_{x}h(x,t)$ and the fluctuation of $\rho(x,t)$
in classical systems, where $h(x,t)$ and $\rho(x,t)$ denote the surface height
and the local density, respectively. The correlation function of
$\delta\rho(x,t)$ of a one-dimensional nonlinear fluctuating hydrodynamics has
earlier been demonstrated to show the same scaling relation as that of
$\partial_{x}h(x,t)$ in KPZ equation Spohn (2014); Mendl and Spohn (2015);
Kulkarni  _et al._ (2015). The quantum surface-height operator analogous to
the classical systems is defined as
$\hat{h}_{j}=\sum_{i=1}^{j}(\hat{c}_{i}^{\dagger}\hat{c}_{i}-\nu)$, where
$\nu=N/L$ denotes the filling fraction with $N$ and $L$ being the total number
of particles and the length of the system, respectively Fujimoto  _et al._
(2020, 2021, 2022). Then one can calculate the average surface height,
$h_{av}=\frac{1}{L}\sum_{j=1}^{L}Tr(\hat{\rho}(t){\hat{h}_{j}})$, where
$\hat{\rho}(t)$ stands for the density matrix. One can further define the
surface roughness $W(L,t)$, which measures the standard deviations of
$\hat{h}_{j}$, given by

$\displaystyle
W(L,t)=\sqrt{\frac{1}{L}\sum_{j=1}^{L}Tr(\hat{\rho}(t)[\hat{h}_{j}-h_{av}]^{2})}.$
(2)

The proposal of FV scaling says that $W(L,t)$ satisfies the following scaling
relation Vicsek and Family (1984); Family and Vicsek (1985); Fujimoto  _et
al._ (2020, 2021, 2022), that is

$\displaystyle W(L,t)=L^{\alpha}~{}g(t/L^{z}),{}$ (3)

where the function $g(x)\sim x^{\beta}$ and $1$ for $x\ll 1$ and $x\gg 1$,
respectively. Hence, Eq. (3) implies that $W(L,t)\propto~{}t^{\beta}$ for
$t\ll t_{sat}$ and $L^{\alpha}$ for $t\gg t_{sat}$, where the crossover time
$t_{sat}=L^{z}$ with the dynamical exponent $z$ being $\alpha/\beta$. As per
this scaling relation mentioned above, one should see a data collapse onto a
single curve for ordinate $W(L,~{}t)$ and abscissa $t$ normalized by
$L^{\alpha}$ and $L^{z}$, respectively. Thus this scaling relation gives rise
to two famous dynamical universality classes: the EW class Antal  _et al._
(2008); Eisler and Rácz (2013); Hunyadi  _et al._ (2004) with
$(\alpha=1/2,\beta=1/4,z=2)$ and the KPZ class Kardar  _et al._ (1986) with
$(\alpha=1/2,\beta=1/3,z=3/2)$. Moreover, $(\alpha=1/2,\beta=1/2,z=1)$
corresponds to the ballistic dynamics. In our work, $(\beta<1/4,z>2)$ and
$(1/4<\beta<1/2,1<z<2)$ are identified with subdiffusive and superdiffusive
dynamical classes, respectively.

Figure 2: Dynamical scaling of the surface roughness ${\boldmath W(L,t)}$ in
the static model: (a-d): The variation of $W(L,t)$ with time $t$ normalized by
$(L/200)^{\alpha}$ and $(L/200)^{z}$, respectively, for system sizes
$L=400,~{}600$ and $800$. (a): Here we consider the vertical critical line
with $\lambda=0.5$, $\mu=~{}2$. The best fitted scaling exponents for this
case are $(\alpha=0.51,\beta=0.24,z=2.13)$, which closely correspond to
diffusive/EW dynamical class. (b-d): The surface roughness dynamics are
examined for the triple point ($\lambda=1$, $\mu=2$), horizontal ($\lambda=1$,
$\mu=1.5$), and tilted critical lines ($\lambda=1.5$, $\mu=3$), respectively.
The best fitted scaling exponents for these cases are
$(\alpha,\beta,z)=(0.54,0.15,3.60)$, $(0.53,0.17,3.12)$, and
$(0.52,0.13,4.00)$, respectively, which belong to subdiffusive dynamical
class.

To study the dynamical behavior of the static model, we consider the many-body
initial state
$|\Psi_{\text{int}}\rangle=\prod_{j=1}^{L/2}c_{2j}^{\dagger}|0\rangle$ at
half-filling, which is a natural choice Antal and Rácz (1996); Antal  _et al._
(2008); Eisler and Rácz (2013); Hunyadi  _et al._ (2004) due to its easy
accessibility in cold atom experiments. First, we examine the surface
roughness dynamics of the vertical critical line ($\lambda=0.5,~{}\mu=2$) that
demonstrate diffusive dynamical behavior with the best fitted scaling
exponents $(\alpha,\beta,z)=(0.51,0.24,2.13)$, as shown in Fig. 2 (a).
Further, the AAH critical point ($\lambda=0$, $\mu=2$) Aubry and André (1980)
also shows the same exponents as shown in Fig. 14 (a), indicating the EW
universality class continuing to persist along the vertical critical line. In
addition to that, this behavior is consistent with the earlier finding Han
_et al._ (1994), which is that the critical exponent $\nu=1$ at self-dual
critical point of the AAH model Aubry and André (1980) continues to be $\nu=1$
along the vertical critical line Han  _et al._ (1994), which we confirm
through a scaling analysis of entanglement entropy as shown in Figs. 13 (a-b).
We then investigate behavior of the triple point $(\lambda=1$, $\mu=2$), where
we observe subdiffusive dynamical behavior with the scaling exponents
$(\alpha,\beta,z)=(0.54,0.15,3.60)$, as shown in Fig. 2 (b). In Figs. 2 (c-d),
the best fitted scaling exponents for the parameter points on the horizontal
($\lambda=1$ and $\mu=1.5$) and tilted critical lines ($\lambda=1.5,~{}\mu=3$)
are seen to be $(0.53,0.17,3.12)$ and $(0.52,0.13,4.00)$, respectively, which
indicate the presence of subdiffusive dynamical class. However, the different
exponents obtained for the subdiffusive cases are somewhat ‘anomalous’ due to
their rare occurrences in the classical systems and clean non-interacting
quantum systems. The scaling exponents in the critical phase also lie under
the subdiffusive dynamical class, e.g. at ($\lambda=1.2$, $\mu=1.8$) exponents
are $(0.53,0.19,2.79)$ as evinced in Fig. 14 (b).

###  III.2 Approximated expression of the surface roughness and the relation
to the bipartite entanglement entropy

As a further confirmation, we also study the dynamics of the square of the
approximate surface roughness $W_{a}^{2}(L,t)$ and the dynamics of bipartite
entanglement entropy $S_{L/2}(t)$. It can be shown that $W_{a}^{2}(L,t)$ is
closely related to the half-chain entanglement entropy $S_{L/2}(L,t)$ Fujimoto
_et al._ (2021) in the case of non-interacting many-body quantum systems
provided the following conditions are being satisfied: i)
$h_{av}(t)~{}\approx~{}0$, ii)
$W^{2}(L,t)~{}\approx~{}Tr(\hat{\rho}(t)\\{h_{L/2}-h_{av}\\}^{2})$, and, iii)
$\sum_{j=1}^{L/2}Tr(\hat{\rho}(t)n_{j})\approx\frac{\nu L}{2}$, where $h$,
$\nu$ stand for the surface-height operator and the filling fraction of the
quantum many-body initial state, respectively. We also observe that these
quantities obey the previously mentioned conditions for our static model. As
example cases, in Figs. 3 (a-d), we show that the conditions (i) and (iii) are
being satisfied for two different sets of parameter values,
$\lambda=1,~{}\mu=1.5$ (a point on horizontal critical line), and
$\lambda=1.5,~{}\mu=3$ (a point on the tilted critical line) for the static
model, i.e., $h_{\rm{av}}(t)$ and $\sum_{j=1}^{L/2}{\rm Tr}(\rho(t)n_{j})$
remain close to zero and $\frac{\nu L}{2}$, respectively, for all time.

(a) (b) (c) (d)

Figure 3: Time evolution of $\mathbf{h_{\rm av}(t)}$ and
$\mathbf{\sum_{j=1}^{L/2}}\mathbf{{\rm Tr}(\rho(t)n_{j})}$: (a-b) The
variation of $h_{\rm{av}}(t)$ with time for $L=200$. The parameter values
chosen for these plots are as follows: Fig. (a): $\lambda=1,~{}\mu=1.5$ (a
point on the horizontal critical line), Fig. (b): $\lambda=1.5,~{}\mu=3$ (a
point on the tilted critical line). In both cases, $h_{\rm{av}}(t)$ remains
close to zero for all time. (c-d) The variation of $\sum_{j=1}^{L/2}{\rm
Tr}(\rho(t)n_{j})$ with time for the parameter values same as Figs. (a-b). In
both cases, $\sum_{j=1}^{L/2}\rm{Tr}(\rho(t)n_{j})$ remains close to
$\frac{\nu L}{2}$, which is $50$ for our case since we have taken the initial
state with $\nu=1/2$ with the system size being $L=200$.

These assumptions further allow us to define the approximate surface roughness
operator in the following manner Fujimoto  _et al._ (2021)

$\displaystyle
W^{2}(L,t)~{}\simeq~{}W_{a}^{2}(L,t)~{}=~{}Tr\left[\hat{\rho}(t)\left(\sum_{i=1}^{L/2}c_{i}^{\dagger}c_{i}-\frac{L\nu}{2}\right)^{2}\right],$
(4)

which is equivalent to the particle-number fluctuations in the half of the
system. For non-interacting many-body quantum systems, one can also argue that
the half-chain entanglement entropy $S_{L/2}(L,t)~{}$
$\approx~{}3W_{a}^{2}(L,t)$ [See supplementary of Ref. Fujimoto  _et al._
(2021) for a detailed derivation], which implies that the surface roughness
can be a good quantifier of entanglement entropy in such systems.

One can therefore employ FV scaling relations both for $W_{a}^{2}(L,t)$ and
$S(L,t)$ to deduce following scaling relations given below

$\displaystyle S(L,t)~{}=~{}L^{\alpha_{1}}s(t/L^{z_{1}})$
$\displaystyle~{}\propto~{}$
$\displaystyle~{}t^{\beta_{1}}~{}~{}~{}t<<t_{sat},$
$\displaystyle~{}\propto~{}$
$\displaystyle~{}L^{\alpha_{1}}~{}~{}~{}t>>t_{sat}.$ $\displaystyle
W_{a}^{2}(L,t)=L^{2\alpha_{2}}w(t/L^{z_{2}})$ $\displaystyle~{}\propto~{}$
$\displaystyle~{}t^{2\beta_{2}}~{}~{}~{}t<<t_{sat},$
$\displaystyle~{}\propto~{}$
$\displaystyle~{}L^{2\alpha_{2}}~{}~{}~{}t>>t_{sat}.$

Similar to the one-parameter scaling relation for the exact surface roughness
operator $W(L,t)$ shown in Eq. (3), one can further argue that the crossover
time for these cases satisfies $t_{sat}=L^{z_{1}(z_{2})}$, where
$z_{1(2)}=\alpha_{1(2)}/\beta_{1(2)}$. Hence, one should also expect a scaling
collapse of $S(L,t)$ ($W_{a}^{2}(L,t)$) vs $t$, followed by a proper
normalization of ordinate and abscissa by $L^{\alpha_{1}}$ ($L^{2\alpha_{2}}$)
and $L^{z_{1}}$ ($L^{z_{2}}$), respectively. As shown in Fig. (4), we
numerically investigate the dynamics of the bipartite entanglement entropy
$(S_{L/2})$ for the same density-wave (DW) initial state to confirm the
exponents obtained from the scaling analysis of the exact surface roughness
dynamics. In Fig. 4 (a), we consider the vertical critical line ($\lambda=0.5$
and $\mu=2$), where we observe diffusive dynamical behaviors with
$(\alpha_{1},\beta_{1},z_{1})=(1.05,0.45,2.33)$, respectively. In addition to
that, we see a similar behavior for the AAH citical point
($\lambda=0,~{}\mu=2$). In Figs. 4 (b-d), we again consider four parameter
points on the critical phase boundaries, namely, the triple point
($\lambda=1,~{}\mu=2$), horizontal critical line ($\lambda=1,~{}\mu=1.5$) and
tilted critical line ($\lambda=1.5,~{}\mu=3$). In all three cases, the system
demonstrates subdiffusive dynamical behaviors with the best fitted scaling
exponents being $(\alpha_{1},\beta_{1},z_{1})=(1.24,0.33,3.76)$,
$(1.15,0.32,3.59)$, and $(1.07,0.29,3.70)$, respectively. In addition to that,
we observe the critical phase demonstrating similar subdiffusive dynamical
behavior (not shown here). These scaling exponents more or less agree with the
scaling exponents obtained from the exact $W(L,t)$ dynamics. In Figs. 5 (a-d),
we further numerically examine the dynamics of $W_{a}^{2}(L,t)$ as defined in
Eq. (4) for the same many-body initial state, which essentially carries the
information about the particle-number fluctuation in the half of the system.
As shown in Figs. 5, the best fitted exponents obtained from this analysis for
the same sets of parameter values are also found to be close to those obtained
from the previous two dynamical measurements.

Figure 4: FV scaling of entanglement entropy at criticalities: (a-d) The
variation of $S_{L/2}$ with time for $L=400,~{}600$ and $800$. The parameter
values chosen for these plots are as follows: Fig. (a):$\lambda=0.5,~{}\mu=2$
(vertical critical line), Fig. (b): $\lambda=1,~{}\mu=2$ (Triple point), Fig.
(c): $\lambda=1,~{}\mu=1.5$ (horizontal critical line), Fig. (d):
$\lambda=1.5,~{}\mu=3$ (tilted critical line). For almost all cases, we see
good agreement with FV scaling relations with the scaling exponents
$(\alpha_{1},\beta_{1},z_{1})$ being: Fig. Fig. (a): $(1.05,0.45,2.33)$, Fig.
(b): $(1.24,0.33,3.76)$, Fig. (c): $(1.15,0.32,3.59)$, Fig. (d):
$(1.07,0.29,3.70)$. Figure 5:  FV scaling of approximate surface roughness
operator $W^{2}(L,t)$ at criticalities: (a-d) Plots of square of the
approximate surface roughness $W^{2}(L,t)$ vs $t$ for the system sizes being
$L=400,~{}600$ and $800$. The plot parameters chosen for these figures are the
same as Figs. 4 (a-d). The best fitted scaling exponents
$(2\alpha_{2},2\beta_{2},z_{2})$ estimated from this analysis are as follows:
(a) $(1.04,0.45,2.31)$, (b) $(1.18,0.32,3.69)$, (c) $(1.18,0.32,3.69)$, (d)
$(1.08,0.29,3.72)$.

###  III.3 single-particle wave packet dynamics in the static model

In addition to the many-body DW dynamics, we also examine the dynamics of the
single-particle root-mean-squared displacement (RMSD), which is given by

$\displaystyle\sigma(t)$ $\displaystyle=$ $\displaystyle\sqrt{\langle
X^{2}(t)\rangle-\langle X(t)\rangle^{2}},$ $\displaystyle\langle
X^{2}(t)\rangle$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{L}(j-j_{0})^{2}~{}|\psi_{j}(t)|^{2},$
$\displaystyle\langle X(t)\rangle$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{L}(j-j_{0})~{}|\psi_{j}(t)|^{2},$ (6)

where the initial state is a single-particle wave packet localized at site
$j_{0}$, and $\psi_{j}(t)$ is the $j$-th component of the time evolved initial
state with the static Hamiltonian shown in Eq. 1. In Figs. 6 (a-d), we examine
the wave packet dynamics at the critical phase boundaries, namely, vertical
critical point ($\lambda=0.5,~{}\mu=2$), triple point ($\lambda=1,~{}\mu=2$),
horizontal critical line ($\lambda=1,~{}\mu=1.5$), and tilted critical line
($\lambda=1.5,~{}\mu=3$). The scaling exponents
($\alpha_{3},~{}\beta_{3},~{}z_{3}$) found from FV scaling for these cases are
found to be $(0.95,0.49,1.94)$, $(1.10,0.32,3.44)$, $(1.05,0.33,3.19)$ and
$(1.10,0.33,3.33)$, respectively, which also respect the FV scaling relation.
Note that the first dynamical exponents again show diffusive dynamical
behaviors as seen in the previous cases. On the other hand, the last three
dynamical exponents lie in the subdiffusive dynamical regimes, which agree
with the behaviors obtained from different dynamical measurements performed
with the many-body DW initial state. We have checked that the many-particle
RMSD calculated for DW intial state also produces same exponents (not shown
here).

Figure 6:  FV scaling from the single particle wave packet dynamics: (a-d) The
dynamics of the single-particle root-mean-squared displacement
$(\sigma(t)=\sqrt{\langle X^{2}(t)\rangle-\langle X(t)\rangle^{2}})$ vs $t$
for system sizes $L=400,~{}600$ and $800$. The parameter values fo these plots
are the same as Fig. 4 (a-d). The single-particle dynamics also respect the FV
scaling relations as seen in Figs. (a-d) with the best fitted scaling
exponents $(\alpha_{3},\beta_{3},z_{3})$ being : Fig(a): $(0.95,0.49,1.94)$,
Fig(b): $(1.10,0.32,3.44)$, Fig(c): $(1.05,0.33,3.19)$, Fig(d):
$(1.10,0.33,3.33)$.

From our numerical analysis we find that
$\beta_{1}=2\beta_{2}\approx\beta_{3}$ and
$\alpha_{1}=2\alpha_{2}\approx\alpha_{3}$. This shows the equivalence between
the growths of entanglement ($S_{L/2}(t)$) and number fluctuations
($W_{a}^{2}(L,t)$) and also establishes a relationship between surface
roughness growth and transport exponents. A similar relationship has been
reported in a recent study for an interacting quasi-periodic model Bhakuni and
Lev (2023). Although for some instances, it has been observed that the
transport exponents show significant disagreement with the surface roughness
exponents, such as in the random-dimer model (RDM) Fujimoto  _et al._ (2021).
Hence, this relationship as we obtain for our model, may not be a universal
one and can be model-dependent. But this needs more serious investigations in
the future.

##  IV Periodically driven model and its phase diagram

For the driven case, we specifically focus on the inter-phase driving across
the critical lines to systemically explore the phase diagram. To do so, we
drive the static Hamiltonian with a square pulse protocol in the following
manner

$\displaystyle H(t)$ $\displaystyle=$ $\displaystyle
H(\lambda+\delta\lambda~{}f(t),\mu+\delta\mu~{}f(t)),$ (7)

where $f(t)=1~{}(-1)$ for $0\leq t\leq T/2$ $(T/2\leq t\leq T)$;
$\delta\lambda$ and $\delta\mu$ denote the amplitudes of periodically driven
quasiperiodic hopping and onsite potential, respectively. $T=2\pi/\omega$
denotes the time period for the stroboscopic time evolution, and $\omega$ is
the driving frequency. For such a square pulse protocol, the Floquet evolution
operator is given by, $U_{F}(T,0)$ = $e^{-iH_{2}T/2}e^{-iH_{1}T/2}$ =
$\sum_{m}e^{-i\theta_{m}}|m\rangle\langle m|$ Blanes  _et al._ (2009); Sen
_et al._ (2021); Bukov  _et al._ (2015), where $H_{1}$ ($H_{2}$) denotes the
Hamiltonian for $0\leq t\leq T/2$ ($T/2\leq t\leq T$). $e^{-i\theta_{m}}$ is
the $m$-th Floquet eigenvalue, and $|m\rangle=\sum_{j}\psi_{m}(j)|j\rangle$
defines the $m$-th single particle Floquet eigenstate with $j$ denoting the
coordinate of the lattice sites. To study the localization properties of a
normalized Floquet eigenstate, we define the inverse participation ratio,
given by, $I_{m}^{(2)}$ = $\sum_{j}|\psi_{m}(j)|^{4}$ Aoki (1986); Mirlin  _et
al._ (2006); Janssen (1998); Evers and Mirlin (2008a). The quantity
$I_{m}^{(2)}$ scales with system size $L$ as $L^{-\eta_{2}}$, where $\eta_{2}$
is the (multi)fractal dimension. $\eta_{2}=0$ (1) for localized (extended)
states whereas $\eta_{2}$ should lie between $0$ and $1$ for critical states
Evers and Mirlin (2008b); Castellani and Peliti (1986).

As shown in Fig. 7 (a-d), we examine the phase diagram of the periodically
model, driven across the critical lines, by extracting $\eta_{2}$ of the
single-particle spectrum-resolved Floquet eigenstates (sorted in the
increasing order of IPR-values), as a function of the fractional eigenstate
index, $m/L$, and driving frequency, $\omega$. In Fig. 7 (a), we first
consider the vertical critical line ($\lambda=0.5$, $\mu=2$) of the static
model, and allow inter-phase driving between delocalized and localized states
with $\delta\lambda=0$ and $\delta\mu=1.5$. We see the Floquet eigenstates
exhibiting critical behavior in the high-frequency regime just like the static
case, coexisting phases in the intermediate-frequency regime, and finally
delocalized behavior Ray  _et al._ (2018); Aditya  _et al._ (2023); Shimasaki
_et al._ (2022) at extremely small $\omega$-values. In Fig. 7 (b), we show
results for the horizontal critical line ($\lambda=1,~{}\mu=1.5$) for drive-
induced mixing between delocalized and critical regimes with $\delta\lambda$ =
0.8 and $\delta\mu=0$. In this case, we observe that the Floquet eigenstates
first exhibit critical behavior in the high-frequency limit (same as the
static model), which, after showing intriguing features at intermediate
frequencies, gradually moves toward complete delocalization with decreasing
$\omega$-values in the slow driving limit. In Fig. 7 (c), the tilted critical
line ($\lambda=1.5$ and $\mu=3$) is subjected to periodic driving along the
horizontal direction with $(\delta\lambda,~{}\delta\mu)$ = $(0,2.5)$, which
only allows mixing between critical and localized phases. In this case, the
critical regime in the high-frequency regime (like the static case) slowly
flows towards the complete localization in the comparatively slow driving
limit. Finally, we drive the triple point $(\lambda,\mu)=(1,2)$ along the
vertical direction with $(\delta\lambda,\delta\mu)=(0.5,0)$, which gives rise
to only critical Floquet eigenstates, as shown in Fig. 7 (d), for all driving
frequencies due to the drive-induced phase-mixing between critical line
(vertical) and critical phase.

Figure 7: Fractal dimension $\eta_{2}$ of the single-particle Floquet
eigenstates: (a-d) The variation of $\eta_{2}$ as a function of the fractional
single-particle Floquet eigenstate index, $m/L$ (where $m$ denotes the
eigenstate index after arranging total $L$ Floquet eigenstates in the
ascending order of their IPR values), and frequency $\omega$ for
$L=500,~{}1000$, $2000$, $3000$ and $4000$. (a): $\lambda=0.5$, $\mu=2$,
$\delta\lambda=0$, and $\delta\mu=1.5$ (Delocalized and localized phase mixing
across the vertical critical line). (b): $\lambda=1$, $\mu=1.5$,
$\delta\lambda=0.8$, and $\delta\mu=0$ (Critical and delocalized phase mixing
across the horizontal critical line). (c): $\lambda$=1.5, $\mu=3$,
$\delta\lambda$ = 0, and $\delta\mu=2.5$ (Critical and localized phase mixing
across the tilted critical line). (d): $\lambda=1.0$, $\mu=2.0$,
$\delta\lambda=0.5$, and $\delta\mu=0$ (Critical line and critical phase
mixing across the triple point).

We then perform similar analysis for several other choices of parameters (not
shown here) and summarize our general findings about the drive-induced phases
in the table shown in Fig. 1. Our analyses establish that the critical phase
is quite fragile when it competes for its existence with delocalized or
localized phase via periodic driving. Moreover, inter-phase driving between
delocalized and localized phases in general flows towards complete
delocalization in the low-frequency regime Ray  _et al._ (2018); Aditya  _et
al._ (2023); Shimasaki  _et al._ (2022), followed by an intricate
intermediate-frequency regime, which we will discuss later. We thus reach a
hierarchy of phases in terms of the stability in the slow driving limit, i.e,
$\text{delocalized}>\text{localized}>\text{critical}$.

##  V Dynamical scaling in the periodically driven model

###  V.1 Surface roughness dynamics in the periodically driven model

Figure 8: Superdiffusive and subdiffusive dynamics of surface roughness caused
by drive-induced mixing between different phases: (a-d): The stroboscopic
dynamics of $W(L,t)$ with time $t$ normalized by $(L/200)^{\alpha}$ and
$(L/200)^{z}$, respectively, for $L=400,~{}600$ and $800$. (a): The vertical
critical line $(\lambda=0.5,~{}\mu=2)$ is considered to allow drive-induced
mixing between delocalized and localized states with
$(\delta\lambda,~{}\delta\mu,\omega)=(0,1.5,5)$. (b-c): The horizontal
$(\lambda=1,\mu=1.7)$ and tilted critical lines $(\lambda,\mu)=(1.5,3)$ are
subjected to periodic driving with
$(\delta\lambda,\delta\mu,\omega)=(0.8,0,0.7)$ and $(1.2,2.4,0.32)$,
respectively, to intermix critical and delocalized states. (d): The horizontal
critical $(\lambda=1,\mu=1.5)$ is driven with
$(\delta\lambda,\delta\mu,\omega)=(0,1,5)$, which causes the critical and
localized phase to intermix with one another. In Figs. (a-c), we observe KPZ-
like superdiffusive surface roughness dynamics with the best fitted scaling
exponents, $(\alpha,\beta,z)$ = $(0.50,0.33,1.50)$, $(0.49,0.31,1.58)$, and
$(0.50,0.33,1.50)$, respectively. In Fig. (d), we see subdiffusive dynamics
with scaling exponents $(\alpha,\beta,z)$ = $(0.45,0.11,4.10)$.

After exploring the IPR phase diagram, we then examine the stroboscopic
dynamics of $W(L,t)$ Fujimoto  _et al._ (2021) in the driven model for the
same initial many-body state. At first, we analyze the inter-phase mixing
between delocalized and critical/localized states across the critical lines,
which gives rise to coexisting phases in the Floquet spectrum in the
intermediate-frequency regime, as shown in Figs. 7 (a-b). This regime is the
most intriguing from the dynamical perspective since drive-induced
subdiffusive and ballistic surface roughness growth in the high-frequency and
low-frequency regime, respectively, are already present in the static model
due to presence of critical and delocalicalized states. Hence, the most
suitable candidates for drive-induced superdiffusive behavior (not present in
static model) are those where the multiple phases can coexist in the spectrum
with a significant fraction of states being delocalized. This may lead to
suppression of the ballistic growth of delocalized states due to the slow
dynamics of localized/ critical states, and the net result can be
superdiffusive in nature Nair  _et al._ (2023). Noting this fact, in Figs. 8
(a-c), the vertical ($\lambda=0.5$, $\mu=2$), horizontal ($\lambda=1$,
$\mu=1.7$), and tilted $(\lambda,\mu)=(1.5,3)$ critical lines are subjected to
periodic driving with $(\delta\lambda,~{}\delta\mu,~{}\omega)=(0,1.5,5)$,
$(0.8,0,0.7)$ and $(1.2,2.4,0.32)$, respectively, to allow phase mixing
between either delocalized and localized phases (vertical) or delocalized and
critical phases (other two critical lines); further, Floquet spectrum for all
cases embrace the property mentioned earlier. The best fitted
$(\alpha,\beta,z)$ for these cases are found to be $(0.50,0.33,1.50)$,
$(0.49,0.31,1.58)$, and $(0.50,0.33,1.50)$, respectively, which implies an
emergent KPZ-like superdiffusive dynamics. Since the delocalized and critical
phases of the static model exhibit ballistic and subdiffusive behaviors, phase
mixing can, in principle, give rise to superdiffusive dynamics, even with KPZ
scaling exponents if driven appropriately, which we observe from Fig. 8 (b-c).
Similar argument also holds for mixing between delocalized and localized
states, as shown in Fig. 8 (a).

We now discuss the surface roughness growth for phase mixing between critical
and localized phases. As evinced in Fig. 8 (d), the horizontal critical line
($\lambda=1,~{}\mu=1.5$) is subjected to periodic driving with
$(\delta\lambda,\delta\mu,\omega)=(0,1,5)$, which allows the phase mixing
between critical and localized phases in the intermediate-frequency regime. In
this case, we observe subdiffusive behavior Lezama  _et al._ (2019) with
$(\alpha,\beta,z)=(0.45,0.11,4.10)$, where the subdiffusive growth exponent of
the static model ($\beta\sim 0.15$) gets further suppressed due to the
presence of the coexisting localized states in the spectrum. Further, the fate
of such phase mixing in the extremely low-frequency regime is mostly complete
localization with no well-defined FV scaling exponents, as shown in Figs. 17.
whereas in the extremely fast driving limit the scaling exponents agree with
the exponents obtained for the static case as shown in Figs 18.

###  V.2 KPZ-like scaling behavior on the vertical critical line subjected to
periodic driving and the effective van Vleck Hamiltonian

To get an analytical insight into the KPZ-like scaling behavior, we consider a
point on the vertical critical line of the static model with
$\lambda=0.5,~{}\mu=2$, and periodically drive it along the horizontal
direction of the static phase diagram with $\delta\lambda=0$, $\delta\mu=1.5$
and $\omega=5$ (Fig. 4(a) in the main text). This is comparatively fast
driving regime, however not extremely fast so that the signatures of periodic
driving can still be present in the stroboscopic Floquet dynamics. These
parameter values allow us to use the van Vleck perturbative expansion Mikami
_et al._ (2016) to find an effective Floquet Hamiltonian. Before proceeding
further, one should note that the stroboscopic Floquet evolution operator for
the square pulse protocol, shown in the main text, possesses the following
symmetry property, which will be important later in the discussion,

$\displaystyle
U_{F}^{-1}(J,\lambda,\mu,\delta\mu,\delta\lambda)=U_{F}(-J,-\lambda,-\mu,\delta\mu,\delta\lambda).$
(8)

This further implies the following symmetry property for the effective Floquet
Hamiltonian

$\displaystyle
H_{F}(J,\lambda,\mu,\delta\mu,\delta\lambda)=-H_{F}(-J,-\lambda,-\mu,\delta\mu,\delta\lambda).$
(9)

We can now find the Floquet Hamiltonian $H_{F}$ in the relatively high-
frequency regime employing the van Vleck perturbation theory. Note that this
is a perturbative expansion in powers of $1/\omega$, and it has the advantage
over the Floquet-Magnus expansion that it does not depend on the phase of the
driving protocol, i.e., it is invariant under $f(t)\to f(t+t_{0})$. The high-
frequency expansion obtained from the van Vleck perturbation theory to the
third order Mikami  _et al._ (2016) is as follows:

$\displaystyle H_{F}$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}H_{F}^{(n)},$ $\displaystyle H_{F}^{(0)}$
$\displaystyle=$ $\displaystyle H_{0},$ $\displaystyle H_{F}^{(1)}$
$\displaystyle=$ $\displaystyle\sum_{m\neq 0}\frac{[H_{-m},H_{m}]}{2m\omega},$
$\displaystyle H_{F}^{(2)}$ $\displaystyle=$ $\displaystyle\sum_{m\neq
0}\frac{[[H_{-m},H_{0}],H_{m}]}{2m^{2}\omega^{2}}\,+\,\sum_{m\neq
0}\sum_{n\neq 0,m}\frac{[[H_{-m},H_{m-n}],H_{n}]}{3mn\omega^{2}},$
$\displaystyle H_{F}^{(3)}$ $\displaystyle=$ $\displaystyle\sum_{m\neq
0}\frac{[[[H_{-m},H_{0}],H_{0}],H_{m}]}{2m^{3}\omega^{3}}+\sum_{m\neq
0}\sum_{n\neq 0,m}\frac{[[[H_{-m},H_{0}],H_{m-n}],H_{n}]}{3m^{2}n\omega^{3}}$
(10) $\displaystyle+\sum_{m\neq 0}\sum_{n\neq
0,m}\frac{[[[H_{-m},H_{m-n}],H_{0}],H_{n}]}{4mn^{2}\omega^{3}}+\sum_{m,m\neq
0}\frac{[[[H_{-m},H_{m}],H_{-n}],H_{n}]}{12mn^{2}\omega^{3}}]$
$\displaystyle+\sum_{m\neq 0}\sum_{n\neq
0,m}\frac{[[H_{-m},H_{0}],[H_{m-n},H_{n}]]}{12m^{2}n\omega^{3}}+\sum_{m,n\neq
0}\sum_{l\neq 0,m,n}\frac{[[[H_{-m},H_{m-l}],H_{l-n}],H_{n}]}{6lmn\omega^{3}}$
$\displaystyle+\sum_{m,n\neq 0}\sum_{l\neq
0,m-n}\frac{[[[H_{-m},H_{m-n-l}],H_{l}],H_{n}]}{24lmn\omega^{3}}+\sum_{m,n\neq
0}\sum_{l\neq
0,m,n}\frac{[[H_{-m},H_{m-l}],[H_{l-n},H_{n}]]}{24lmn\omega^{3}},$

where $H_{m}$ stands for the $m$-th Fourier component of $H(t)$. The zeroth-
order Floquet effective Hamiltonian given by van Vleck perturbation theory for
our case is

$\displaystyle H_{F}^{(0)}$ $\displaystyle=$
$\displaystyle\frac{1}{T}~{}\int_{0}^{T}dt~{}H(t)=\sum_{j}[(1+\lambda\cos(2\pi\beta(j+1/2)+\phi)~{}c_{j}^{\dagger}c_{j+1}+{\rm
H.c.}]+\mu\sum_{j}\cos(2\pi\beta j+\phi)~{}c_{j}^{\dagger}c_{j}.$ (11)

Now, the $m$-th Fourier component of $H(t)$ turns out to be

$\displaystyle H_{m}$ $\displaystyle=$
$\displaystyle\frac{2i\delta\mu}{m\pi}\sum_{j}\cos(2\pi\beta
j+\phi)c_{j}^{\dagger}c_{j},~{}~{}{\rm for}~{}~{}m~{}~{}\text{odd},$ (12)
$\displaystyle=$ $\displaystyle 0,~{}~{}~{}~{}{\rm
for}~{}~{}~{}~{}~{}~{}~{}~{}~{}m~{}~{}\text{even},$

and $H_{0}=H_{F}^{(0)}$. We then find that the first-order Floquet
Hamiltonian, which is given by $H_{F}^{(1)}=\sum_{m\neq
0}\left[H_{-m},H_{m}\right]/(2m\omega)$, and is zero for all $m$. The second-
order Floquet Hamiltonian consists of two terms

$\displaystyle H_{F}^{(2)}$ $\displaystyle=$ $\displaystyle\sum_{m\neq
0}\frac{[[H_{-m},H_{0}],H_{m}]}{2m^{2}\omega^{2}}$ (13)
$\displaystyle+\sum_{m\neq 0}\sum_{n\neq
0}\frac{[[H_{-m},H_{m-n}],H_{n}]}{3mn\omega^{2}}.$

From our calculation, $H_{F}^{(2)}$ is found to be

$\displaystyle H_{F}^{(2)}$ $\displaystyle~{}=~{}$
$\displaystyle-\frac{\delta\mu^{2}\sin^{2}(\pi\beta)}{4\omega^{2}}\sum_{j}\sin^{2}(\pi\beta(2j+1)+\phi)\times$
$\displaystyle\left(1+\lambda\cos(2\pi\beta(j+1/2)+\phi)\right)c_{j}^{\dagger}c_{j+1}+\rm{H.c}.$

Therefore, the effective Hamiltonian to second order in $(1/\omega)$ will be
$H_{F}=H_{F}^{(0)}+H_{F}^{(2)}$. We further note that the third order
effective Hamiltonian appears to be zero for our model for this particular
choice of the driving protocol, which can be shown by the following argument.
The general expression of the third order van Vleck effective Hamiltonian
contains multiple terms; however, only two terms which are important for our
calculation are as follows (the rest of the terms immediately go to zero due
to the fermionic commutation relations),

$\displaystyle H_{F}^{(3)}$ $\displaystyle=$ $\displaystyle\sum_{m\neq
0}\frac{[[[H_{-m},H_{0}],H_{0}],H_{m}]}{2m^{3}\omega^{3}}$
$\displaystyle+\sum_{m\neq 0,~{}n\neq
0,~{}m}\frac{[[[H_{-m},H_{0}],H_{m-n}],H_{n}]}{3m^{2}n\omega^{3}}.$

Moreover, the first term shown in Eq. (LABEL:HF3) can be shown to be zero due
to the symmetry property discussed in Eq. (9) Aditya  _et al._ (2023).
Furthermore, the second term also reduces to zero using the fact that the
Fourier components, $H_{m}$, $H_{n}$ and $H_{m-n}$ simultaneously can not be
non-zero due to the constraints imposed by the Fourier expansion, as shown in
Eq. (12) Aditya  _et al._ (2023). Therefore, the effective Hamiltonian to
second order in $1/\omega$ is a good approximation to the exact Floquet
effective Hamiltonian. With this effective Hamiltonian, we can now study the
surface roughness dynamics stroboscopically in time.

Figure 9: KPZ-like growth from the effective van Vleck Hamiltonian: (a) The
variation of exact $W(L,t)$ in time for the
$\lambda=0.5,~{}\mu=2,~{}\delta\mu=1.5,~{}\omega=5$, (b) The variation of
approximate $W_{a}^{2}(L,t)$ in time for the same set of parameter values. The
estimated FV scaling exponents are found to be as follows: Fig. (a):
$(\alpha,~{}\beta,~{}z)=(0.45,~{}0.30,~{}1.50)$, Fig. (b):
$(2\alpha_{2},~{}2\beta_{2},~{}z_{2})=(0.91,~{}0.67,~{}1.36)$.

In Fig. 9 (a-b), we analyze the variation of exact surface roughness and the
square of the approximate surface roughness operators, $W(L,t)$ and
$W_{a}^{2}(L,t)$ vs time, for a point on the vertical critical line of the
static model subjected to periodic driving with parameter values,
$\lambda=0.5,~{}\mu=2,~{}\delta\mu=1.5$, and $\omega=5$ with the effective
Floquet Hamiltonian, shown in Eq. (11) and Eq. (LABEL:HF2). Remarkably, the
best fitted scaling exponents obtained from these analysis for these parameter
values are found to be $(\alpha,~{}\beta,~{}z)=(0.45,~{}0.30,~{}1.50)$, and
$(2\alpha_{2},~{}2\beta_{2},~{}z_{2})=(0.91,~{}0.67,~{}1.36)$, respectively,
which show the superdiffusive surface roughness dynamics with close KPZ-like
scaling behavior. Moreover, these scaling exponents are quite close to the
exponents obtained from exact numerical Floquet analysis, revealing that this
effective Floquet Hamiltonian obtained from van Vleck expansion is able to
capture the feature of the exact numerical Floquet analysis.

###  V.3 Effective Floquet Hamiltonian in large driving amplitude limit

In this section, we will discuss both the IPR-based phase diagram and
dynamical behaviors of the periodically driven model in the large driving
amplitude limit. To do so, we will employ the Floquet perturbation theory
(FPT) Soori and Sen (2010); Sen  _et al._ (2021) to find an effective Floquet
Hamiltonian for $\delta\mu\gg J,~{}\mu,~{}\lambda$. We further choose
$\delta\lambda=0$ since this choice of parameter values are easy to handle for
our analytical treatment. For this particular case, the unperturbed
Hamiltonian, $H_{0}$ is as follows

$\displaystyle H_{0}=\delta\mu f(t)\sum_{j}\cos(2\pi\beta j+\phi)n_{j},$ (16)

which is diagonalizable in real space having the instantaneous energy
eigenvalues, $E_{j}=\delta\mu f(t)\cos(2\pi\beta j+\phi)$ in the single-
particle limit, where $j$ stands for real space lattice indices ranging
between 1 and $L$. These instantaneous eigenvalues satisfy the following
condition:

$\displaystyle e^{i\int_{0}^{T}dt\left[E_{j}(t)-E_{j^{\prime}}(t)\right]}=1,$
(17)

where $T=2\pi/\omega$ is the driving period. This condition implies that one
needs to carry out degenerate FPT Sen  _et al._ (2021); Soori and Sen (2010)
to find the effective Floquet Hamiltonian. The eigenfunction corresponding to
a particular $E_{j}$ in the single-particle limit is given by

$\displaystyle|j\rangle=\hat{c}_{j}^{\dagger}|0\rangle.$ (18)

To chart out the main steps of FPT for our periodically driven model, we begin
with the Schrödinger equation

$i\frac{d|\psi\rangle}{dt}~{}=~{}\left(H_{0}+V\right)|\psi\rangle,$ (19)

where $V$ is given by

$\displaystyle V$ $\displaystyle=$
$\displaystyle\sum_{j}\left(1+\lambda\cos(2\pi\beta(j+1/2)+\phi)\right)~{}c_{j}^{\dagger}c_{j+1}+\rm{H.c.}$
(20) $\displaystyle+\mu\sum_{j}\cos(2\pi\beta j+\phi)n_{j}.$

We further assume that $|\psi(t)\rangle$ has the following expansion

$\displaystyle|\psi(t)\rangle~{}=~{}\sum_{j}~{}c_{j}(t)~{}e^{-i\int_{0}^{t}dt^{\prime}E_{j}(t^{\prime})}|j\rangle.$
(21)

Substituting Eq. (21) in Eq. (19), we find that

$\frac{dc_{j}}{dt}~{}=~{}-i~{}\sum_{j^{\prime}}~{}\langle
j|V|j^{\prime}\rangle~{}e^{i\int_{0}^{t}dt^{\prime}(E_{j}(t^{\prime})-E_{j^{\prime}}(t^{\prime}))}c_{j^{\prime}}(t).$
(22)

One can further recast, $c_{j}=\sum_{n}\,c_{j}^{(n)}$ where $c_{j}^{(n)}$
denotes the coefficient proportional to $V^{n}$. One can now solve Eq. (22)
order by order in $V$. Using the expansion of $c_{j}$ in orders of $V$
mentioned earlier and comparing both sides of Eq. (22), one can easily show
that the zeroth order coefficient, $c_{j}^{(0)}$ satisfies

$\displaystyle\frac{dc_{j}^{(0)}(t)}{dt}$ $\displaystyle=$ $\displaystyle 0.$
(23)

In a similar manner, it can be shown that the first order coefficient,
$c_{j}^{(1)}$ satisfies the equation given below

$\displaystyle\frac{dc_{j}^{(1)}(t)}{dt}$ $\displaystyle=$
$\displaystyle-i\sum_{j^{\prime}}\langle
j|V|j^{\prime}\rangle~{}e^{i\int_{0}^{t}dt^{\prime}(E_{j}(t^{\prime})-E_{j^{\prime}}(t^{\prime}))}c_{j^{\prime}}^{(0)}(t).$

One can now integrate both sides of Eq. (23) to show that
$c_{j}^{(0)}(t)=c_{j}^{(0)}(0)$. Moreover, substituting
$c_{j}^{(0)}(t)=c_{j}^{(0)}(0)$ in (LABEL:c1), one obtains the following:

$\displaystyle\frac{dc_{j}^{(1)}(t)}{dt}$ $\displaystyle=$
$\displaystyle-i\sum_{j^{\prime}}\langle
j|V|j^{\prime}\rangle~{}e^{i\int_{0}^{t}dt^{\prime}(E_{j}(t^{\prime})-E_{j^{\prime}}(t^{\prime}))}c_{j^{\prime}}^{(0)}(0).$

Integrating both sides of Eq. (LABEL:c11), $c_{j}^{(1)}(t)$ turns out to be

$\displaystyle c_{j}^{(1)}(t)$ $\displaystyle=$ $\displaystyle
c_{j}^{(1)}(0)-i\sum_{j^{\prime}}\langle
j|V|j^{\prime}\rangle\int_{0}^{t}\,dt^{\prime}\times$ (26) $\displaystyle
e^{i\int_{0}^{t^{\prime}}dt^{\prime}(E_{j}(t^{\prime\prime})-E_{j^{\prime}}(t^{\prime\prime}))}c_{j^{\prime}}^{(0)}(0).$

Without any loss of generality, we further set $c_{j}^{(1)}(0)$ equal to be
zero. Therefore, $c_{j}$ to first order in V after one stroboscopic time
evolution reduces to

$\displaystyle c_{j}(T)$ $\displaystyle=$ $\displaystyle
c_{j}^{(0)}(T)+c_{j}^{(1)}(T),$ $\displaystyle=$ $\displaystyle
c_{j}^{(0)}(0)-i\sum_{j^{\prime}}\langle
j|V|j^{\prime}\rangle\int_{0}^{T}\,dt^{\prime}\,e^{i\int_{0}^{t}dt^{\prime}(E_{j}(t^{\prime})-E_{j^{\prime}}(t^{\prime}))}c_{j^{\prime}}^{(0)}(0),$
$\displaystyle=$
$\displaystyle\left(\delta_{jj^{\prime}}-i\sum_{j^{\prime}}\langle
j|V|j^{\prime}\rangle\int_{0}^{T}\,dt^{\prime}\,e^{i\int_{0}^{t}dt^{\prime}(E_{j}(t^{\prime})-E_{j^{\prime}}(t^{\prime}))}\right)\times$
$\displaystyle c_{j^{\prime}}^{(0)}(0).$

Noting the fact, due to the Floquet degenerate perturbation theory Sen  _et
al._ (2021), $c_{j}(0)=c_{j}^{(0)}(0)$, Eq. (LABEL:cm2) can be further recast
as

$\displaystyle c_{j}(T)$ $\displaystyle~{}=~{}$
$\displaystyle\sum_{j^{\prime}}\left(I-iH_{F}^{(1)}T\right)_{jj^{\prime}}c_{j^{\prime}}(0),$
$\displaystyle\langle j|H_{F}^{(1)}|j^{\prime}\rangle$ $\displaystyle~{}=~{}$
$\displaystyle\frac{\langle
j|V|j^{\prime}\rangle}{T}\,\int_{0}^{T}\,dt\,e^{i\int_{0}^{t}dt^{\prime}(E_{j}(t^{\prime})-E_{j^{\prime}}(t^{\prime}))},$

where $I$ denotes the identity matrix and $H_{F}^{(1)}$ is the Floquet
Hamiltonian to first order in $V$. The non-zero matrix elements of the
effective Hamiltonian $H_{F}$ are as follows:

$\displaystyle\langle j|H^{(1)}_{F}|j+1\rangle$ $\displaystyle=$
$\displaystyle\left(1+\lambda\cos(2\pi\beta(j+1/2)+\phi)\right)P_{j},$
$\displaystyle\text{where}~{}~{}P_{j}$ $\displaystyle=$ $\displaystyle
e^{i\gamma\sin(\pi\beta(2j+1)+\phi)}~{}\frac{\sin(\gamma\sin(\pi\beta(2j+1)+\phi)}{\gamma\sin(\pi\beta(2j+1)+\phi)},$
$\displaystyle\langle j|H^{(1)}_{F}|j\rangle$ $\displaystyle=$
$\displaystyle\mu\cos(2\pi\beta j+\phi),$

where $\gamma=\frac{\delta\mu T}{2}\sin(\pi\beta)$. Hence, the first order
effective Hamiltonian obtained from the first order FPT is given by

$\displaystyle H_{F}^{(1)}$ $\displaystyle=$
$\displaystyle\sum_{j}\left(P_{j}\times\left(J+\lambda\cos(2\pi\beta(j+1/2)+\phi)\right)c_{j}^{\dagger}c_{j+1}+{\rm
H.c.}\right)+\mu\sum_{j}\cos(2\pi\beta j+\phi)n_{j}.$

One can further show that the second order effective FPT Hamiltonian,
$H_{F}^{(2)}$ vanishes since $H_{F}$ can only have terms which are odd powers
in $J,~{}\lambda$, and $\mu$ due to the symmetry property discussed in Eq. (9)
Aditya  _et al._ (2023). Consequently, the correction next to first order will
be third order, and hence, the first-order effective FPT Hamiltonian is a good
approximation to the exact effective Floquet Hamiltonian.

(a) (b)

Figure 10: Fractal dimension of Floquet eigenstates obtained from exact
numerics and effective FPT Hamiltonian: (a-b) The variation of $\eta_{2}$ as a
function of fractional single-particle Floquet eigenstate index, $m/L$ and
driving frequency, $\omega$ for $J=1,~{}\lambda=0.5$, $\mu=2$, and
$\delta\mu=4$. Fig. (a) is obtained from the exact numerical Floquet analysis
while Fig. (b) is found from the first order effective FPT Hamiltonian. The
dependence of $\eta_{2}$ as a function of $m/L$ and $\omega$ for the first
analysis qualitatively agrees with the second one when
$10^{-2}<\omega<10^{2}$. However, the effective FPT Hamiltonian fails to
capture the delocalization phase when $\omega<10^{-2}$ since the effective
Hamiltonian is not uniquely defined at such a low frequency regime due to the
branch cut ambiguities.

In Fig. 10, we compare the phase diagram of single-particle spectrum-resolved
Floquet eigenstates as a function of $\omega$ obtained from the exact
numerical Floquet analysis with the one acquired from the first-order FPT
Hamiltonian for $J=1,~{}\lambda=0.5,~{}\mu=2$ and $\delta\mu=4$. In Figs. 10
(a)-10 (b), we see the variation of $\eta_{2}$ (fractal dimension extracted
from the scaling of IPR, $I_{m}^{(2)}$, which are sorted in the increasing
order of IPR-values for $L=500,~{}1000$ and $2000$) as a function of $\omega$
obtained from the exact numerical Floquet analysis and the first order FPT
Hamiltonian, respectively. Both of the results qualitatively agree with each
other when $\omega>10^{-2}$. Moreover, both of them in the extremely high
frequency regime exhibit critical behavior (similar to the static case as
expected), and then slowly flows towards the localized phase with decreasing
$\omega$-values. However, when $\omega<10^{-2}$, $\eta_{2}$ phase diagram
obtained from the exact numerical Floquet analysis shows almost delocalized
behavior, whereas the first-order FPT Hamiltonian fails to capture this
particular regime. The reason for this discrepancy is due to the extremely
slow driving limit, which forbids to define the effective FPT Hamiltonian,
$H_{F,~{}FPT}$ in an unambiguous manner due to the following reason. The
effective Floquet Hamiltonian, $H_{F}$ in this limit cannot be uniquely
defined since the eigenvalues $e^{-i\theta_{m}}$ of the Floquet operator $U$
will not satisfy $|\theta_{m}|\ll\pi$ for all the Floquet eigenstates $m$
(also see Sec. VI). Hence $H_{F}=(i/T)\ln U$ will suffer from branch cut
ambiguities.

Nevertheless, this re-entrant delocalized behavior in extremely slow driving
regime observed in the exact numerical Floquet phase diagram is somewhat
expected due to the following argument. Note that the driving parameters
chosen for this case ($\delta\mu=4$) allows the phase-mixing between the
positive $\mu$ and negative $\mu$-direction about the vertical critical line
of the static model. To get insight into this phase-mixing, we consider the
maximally expanded phase diagram of the static model along the negative
$\mu$-axis discussed in Sec. A. Taking it into account, one can clearly see
that periodic driving for this particular set of parameter values causes the
delocalized and localized phases to intermix with one another. As discussed in
the main text, this kind of phase mixing initially impels the system toward
the localized behavior before finally reaching drive-induced delocalized
phase, which we exactly see in Fig. 10 (a).

Figure 11: Surface roughness dynamics obtained from exact numerical Floquet
analysis and first order FPT Hamiltonian: (a-b) The stroboscopic time
evolution of the square of the approximate surface roughness operator,
$W^{2}(L,~{}t)$ vs time obtained from the exact numerical Floquet calculation
and first order FPT Hamiltonian, respectively. In both cases, we consider
$J=1,~{}\lambda=0.5,~{}\mu=2.0$, $d\mu=4$, and $L=400,~{}600,~{}800$. In the
first case, we consider $\omega=6.1$, and for the second one, we consider
$\omega=5.7$. The best fitted FV scaling exponents
($2\alpha_{2},~{}2\beta_{2},~{}z_{2}$) for these cases are as follows: Fig.
(a): $(0.92,0.74,1.24)$, Fig. (b): $(0.90,0.75,1.20)$.

(a) (b) (c) (d)(e) (f)

Figure 12: Rescaled height distribution for different dynamical phases: (a-f)
Rescaled surface-height distribution $P(\chi)$ obtained for different
dynamical phases for $L=400$ and for 5000 disorder realizations. We consider
many-body DW state as the initial state to obtain these plots. The parameter
values considered for these cases are (a): $\lambda=0,\mu=3,t=20$ (localized
case), (b): $\lambda=1,\mu=1.5,t=1000$ (subdiffusive dynamics), (c):
$\lambda=0.5,\mu=2,t=500$ (diffusive dynamics), (d):
$\lambda=0.5,\mu=2,\delta\mu=1.5,\omega=5,t=200$ (drive-induced KPZ behavior),
(e): $\lambda=0.5,\mu=2,\delta\mu=4,\omega=5.7$, and $t=500$ (Drive-induced
non-KPZ superdiffusive), (f): $\lambda=0,\mu=1$, $t=100$ (ballistic growth
dynamics). Moreover, in Fig. (d), we fit the left and right tails of $P(\chi)$
to show the agreement of the asymptotic behavior ($\chi\pm\infty$) with the
Tracy-Widom distribution.

In Figs. 11 (a)-11 (b), we also compare the square of the approximate surface
roughness dynamics stroboscopically in time obtained from the exact numerical
Floquet analysis with the one procured from the first order FPT Hamiltonian,
respectively. For both cases, we have taken $J=1,~{}\lambda=0.5,~{}\mu=2$, and
$\delta\mu=4$ (the vertical critical line of the static model subjected to
periodic driving along the horizontal direction with large $\delta\mu$) for
$L=400,~{}600,~{}800$. As shown in Fig. 11, the best fitted FV scaling
exponents obtained from the exact numerical Floquet analysis for $\omega=6.1$,
$(2\alpha_{2},2\beta_{2},~{}z_{2})=(0.92,~{}0.74,~{}1.24)$ are quite close to
those from the first order FPT Hamiltonian for $\omega=5.7$ with
$(2\alpha_{2},~{}2\beta_{2},~{}z_{2})=(0.90,~{}0.75,~{}1.20)$. This
essentially reflects the efficiency of FPT analysis capturing the effective
dynamics in this regime.

##  VI Rescaled height distribution for different dynamical phases

In this section, we study the rescaled surface height distribution for
different dynamical universality classes appearing in our model in the
presence and absence of periodic driving. The rescaled height operator for
this analysis is defined as follows:

$\displaystyle\chi$ $\displaystyle=$ $\displaystyle\frac{h_{L/2}-\langle
h_{L/2}\rangle}{\sigma_{h_{L/2}}},$ (30)

where $h_{L/2}$ stands for the values of quantum surface height for site
$j=L/2$ (defined in the main text) for different disorder realizations.
$\langle h_{L/2}\rangle$ and $\sigma_{h_{L/2}}$ denote the mean and standard
deviation, respectively of the distribution of $h_{L/2}$. In Fig. 12, we
examine the rescaled quantum surface-height distribution $P(\chi)$, with zero
mean and unit standard deviation, for different dynamical phases appeared in
our quantum model. In Fig. 12 (a-f), $P(\chi)$ is shown for the localized,
subdiffusive, diffusive, drive-induced KPZ, drive-induced superdiffusive (non-
KPZ) and ballistic dynamical classes, respectively, with the parameters
suitably chosen as discussed in the caption of Fig. 12. Interestingly,
$P(\chi)$ for each dynamical class appears to be quite distinct from one
another, as can be seen from Fig. 12. In the localized case (Fig. 12 (a)), the
distribution $P(\chi)$ is quite narrow with peak appearing around $\chi\approx
0$ whereas in the subdiffusive class (Fig. 12 (b)) $P(\chi)$ is broadened and
bi-modal with peaks appearing away from $\chi=0$. In the diffusive phase (Fig.
12 (c)) a single-peaked broad distribution is observed. Similar broadened
distributions are shown in Fig. 12 (d) and Fig. 12 (e) for the drive-induced
KPZ and non-KPZ superfiffusive phases, respectively, although they have their
own characteristics which we will discuss later in this section. In the
ballistic class (Fig. 12 (f)), $P(\chi)$ is not very well-behaved and does not
seem to have well-defined asymptotic limits at $\chi\rightarrow\pm\infty$
unlike in the other phases shown here. We note that these distributions, at
the moment, are mere numerical observations and may not be the universal ones.
However, this demands more of similar (and more careful) studies on various
models to come to a conclusion in the future.

Furthermore, in Fig 12 (d), we compare the asymptotic behaviors of $P(\chi)$
for the drive-induced quantum KPZ case with the Tracy-Widom (TW) distribution
function. The fitting analysis shows that the asymptotic behavior for both
$\chi\pm\infty$ follows the TW asymptotes with some fluctuations. Moreover,
the drive-induced KPZ distribution appears to be quite distinct from the
drive-induced non-KPZ superdiffusive case, which seems to have a more
symmetric distribution as shown in Fig. 12 (e). It is already known from the
literature that the asymptotic behaviors of the rescaled height distribution
in the case of classical KPZ universality class follow the TW distribution of
the Gaussian unitary ensembles (GUE) for a particular choice of initial
condition Takeuchi  _et al._ (2011). Recently, this agreement has been also
found in the density fluctuations of single-particle wave packet dynamics in
the case of two-dimensional Anderson localization Mu  _et al._ . This very
similar asymptotic behavior for our drive-induced quantum KPZ case for a many-
body (DW type) initial state, therefore, opens up new possibilities for
further explorations in the future.

##  VII Discussion and conclusion

In summary, we present an intricate dynamical phase diagram of a one-
dimensional fermionic model with correlated disorder, which includes showing
the existence of the diffusive and subdiffusive dynamical universality classes
of surface roughness with anomalous FV scaling exponents, especially in the
critical phase, which is hitherto not reported in the literature due to the
rare occurrences of critical phases in the familiar models. Interestingly, the
static model has been recently realized experimentally on a superconducting
quantum processor with tunable couplers Li  _et al._ (2023) and also using
ultra-cold atoms Xiao  _et al._ (2021), which can be used as test beds for
verifying our results. We also observe that driving this model with periodic
square pulse protocol can induce emergent KPZ-like superdiffusive dynamics for
certain choices of driving parameters. This is the first observation of this
special type of superdiffusive growth dynamics in a periodically driven free
fermionic model. However, this emergent KPZ-like behavior, apparently, can not
be explained by taking the semi-classical limit due to the quadratic nature of
the (driven) model Jin  _et al._ (2020); Ljubotina  _et al._ (2019); Ilievski
_et al._ (2018); Ye  _et al._ (2022); Keenan  _et al._ (2022) and hence opens
up scope for further understanding. One can argue that this behavior is
intrinsically quantum in nature and arises due to the intricate interplay
between periodic driving and the quantum nature of the system (Many-body
fermionic states whose evolution governed by the Schrodinger equation), which
essentially controls the growth in dynamics, and therefore, has no classical
counterpart.

Our study reveals that the interplay between periodic driving and correlated
disorder can give rise to new dynamical universality classes Lezama  _et al._
(2019) of surface roughness, which seem to be different from the known ones,
and therefore, require further understanding. For the future purpose, it would
be worthwhile to examine the surface roughness dynamics in the presence of
interaction, which might open up new possibilities to unravel the surface
growth physics of the many-body localization Basko  _et al._ (2006); Vosk and
Altman (2013); Huse  _et al._ (2013); Oganesyan and Huse (2007); Pal and Huse
(2010); Schreiber  _et al._ (2015) and interacting Floquet quantum matter
Ponte  _et al._ (2015); Bordia  _et al._ (2017); Sierant  _et al._ (2023);
Decker  _et al._ (2020); Khemani  _et al._ (2016); Else  _et al._ (2016);
Zhang  _et al._ (2017).

## Acknowledgments

We would like to thank Diptiman Sen for insightful discussions and for
valuable comments on the manuscript. We thank Sumilan Banerjee for fruitful
discussions. S.A. thanks MHRD, India for financial support through a PMRF. NR
acknowledges support from the Nanyang Technological University and Indian
Institute of Science, where the work has been started.

APPENDIX

##  Appendix A Symmetries of the static model

We begin by discussing the symmetries of the static model (where $\lambda>0$
and $\mu>0$) as considered in the main text, which will allow us to extend the
phase diagram for $\lambda<0$ and $\mu<0$. Moreover, this symmetry analysis
will be useful later in the discussion to understand the drive-induced phase
mixing in the large driving amplitude limit (see Sec. V.3). The symmetries of
the static model are discussed below.

  1. 1.

It should be noted from the phase diagram shown in Fig. 1(a) of the main text
that all the states are delocalized for ($0<\lambda<1$, $0<\mu<2$), localized
for ($\mu>2$, $\mu>2\lambda$), and critical for ($\lambda>1$, $\mu<2\lambda$).
One can further argue that the spectrum will be symmetric about $E=0$ by the
argument given below. Followed by the particle-hole transformation, i.e.,
$c_{j}~{}\rightarrow~{}c_{j}^{\dagger}$ and
$c_{j}^{\dagger}~{}\rightarrow~{}c_{j}$, the static Hamiltonian shown in Eq. 1
of the main text (without loss of generality, we set $\phi$=0 for this
analysis) takes the following form:

$\displaystyle H$ $\displaystyle=$
$\displaystyle\sum_{j}\left(J+\lambda\cos(2\pi\beta(j+1/2))\right)(c_{j}c_{j+1}^{\dagger}+c_{j+1}c_{j}^{\dagger})$
$\displaystyle+$ $\displaystyle\mu\sum_{j}\cos(2\pi\beta
j)~{}c_{j}c_{j}^{\dagger},$ (31) $\displaystyle=$
$\displaystyle-H+\mu\sum_{j}\cos(2\pi\beta j).$

Now, for sufficiently large system sizes, $\sum_{j}\cos(2\pi\beta j)$
necessarily goes to zero, which, therefore, suggests that the energy spectrum
will be symmetric about $E=0$. Moreover, this transformation is equivalent to
the following changes of the parameters of the static model,
$J~{}\rightarrow~{}-J$, $\lambda~{}\rightarrow~{}-\lambda$, and
$\mu~{}\rightarrow~{}-\mu$.

  2. 2.

Next, one can perform the following symmetry operation on the annihilation
operator, namely $c_{j}~{}\rightarrow~{}(-1)^{j}c_{j}$, which is necessarily
equivalent to $J~{}\rightarrow~{}-J$, $\lambda~{}\rightarrow~{}-\lambda$, and
$\mu~{}\rightarrow~{}\mu$. This symmetry transformation along with the
previous one imply $J~{}\rightarrow~{}J$, $\lambda~{}\rightarrow~{}\lambda$,
and $\mu~{}\rightarrow~{}-\mu$.

  3. 3.

There is one more symmetry transformation, which one can perform using the
property of an irrational number. First, one should note that both
quasiperiodic hopping and quasiperiodic potential contain $\cos(2\pi\beta j)$
kind of terms, where $\beta$ is the golden ratio. Under the rational
approximation, $\beta$ is equal to $\frac{P}{L}$, where $L$ is the system size
and $P$ is the integer closest to $(\frac{\sqrt{5}-1}{2})L$. Using this fact,
one can further perform the following transformation, namely,
$j~{}\rightarrow~{}j+\frac{L}{2P}$, followed by which, both quasiperiodic
potential and quasiperiodic hopping change signs. Therefore, this
transformation suggests $J~{}\rightarrow~{}J$, $\mu~{}\rightarrow~{}-\mu$, and
$\lambda~{}\rightarrow~{}-\lambda$.

Taking these three transformations into account, one can now maximally expand
the phase diagram along the negative $\lambda$ and negative $\mu$ direction of
the $\lambda-\mu$ parameter plane. The extension of the phase diagram along
the negative axes indicates the following:

  1. 1.

Taking the second symmetry into account, the phase diagram in the second
quadrant of the $\lambda-\mu$ plane tells us: there will be a delocalized
regime for ($0<\lambda<1$, $-2<\mu<0$), a localized phase for
($\mu<-2,~{}-2\lambda$), and a critical phase for ($2\lambda>2,~{}-\mu$).

  2. 2.

Taking the third symmetry into consideration, the phase diagram can be
extended in the third quadrant of the $\lambda-\mu$ plane: a delocalized
regime for ($-1<\lambda<0$, $-2<\mu<0$), a localized regime for
($\mu<-2,~{}2\lambda$) and a critical regime for ($2\lambda<-2,~{}\mu$).

  3. 3.

Including all three symmetry transformation discussed above implies that phase
diagram will remain unaltered for $J~{}\rightarrow~{}J$,
$\lambda~{}\rightarrow~{}-\lambda$ and $\mu~{}\rightarrow~{}\mu$. With this
transformation, one can further extend the phase diagram in the fourth
quadrant of the $\lambda-\mu$ parameter plane. In this case, the phase diagram
will be the following: a delocalized regime for ($-1<\lambda<0$, $0<\mu<2$), a
localized regime for ($\mu>2,~{}-2\lambda$), and a critical regime for
($2\lambda<-2,~{}-\mu$).

##  Appendix B Finite-size scaling of the steady-state entanglement entropy
and phase transitions in the static model

In this section, we will discuss the finite-size scaling analysis of the
entanglement entropy ($S_{\infty}$) of the nonequilibrium steady-state on the
critical phase boundaries, i.e., the AAH point, the triple point, and three
critical lines. To perform the appropriate scaling analysis, we normalize
$S_{\infty}$ by the thermal page value $[S_{T}=(L/2)\ln 2-0.5]$, which is
expected to be scale invariant at criticality. Therefore, the finite-size
scaling ansatz for the entanglement entropy is as follows Wang  _et al._
(2021)

$\displaystyle\frac{S_{\infty}}{S_{T}}=f(\delta\tau\,L^{1/\nu}).$ (32)

where $S_{\infty}$ is the saturation value of the entanglement entropy after
being approached the long time limit, $L$ stands for the system size. the
correlation length $\xi$ diverges at the critical phase boundaries as
$\xi~{}\sim~{}\delta\tau^{-{\nu}}$, where $\delta\tau=|\tau-\tau_{c}|$ is the
distance from the critical value in the parameter space, and ${\nu}$ denotes
the critical exponent associated with the correlation length.

Figure 13: Finite-size scaling of the saturation value of entanglement entropy
$S_{\infty}/S_{T}$ : Fig. (a): Scaling collapse for transition around the AAH
point $(\lambda=0$, $\mu=2)$, Fig. (b): similar analysis around a point
($\lambda=0.5$, $\mu=2$) on the vertical critical line, Fig. (c): for a point
$(\lambda=1$, $\mu=1.5)$ on the horizontal critical point while approaching
from the delocalized side of the phase diagram, Fig. (d-e): same analysis for
the triple point with $\lambda=1$ and $\mu=2$, while the critical line being
approached from the delocalized and localized sides of the phase diagram,
respectively. Fig. (f): similar analysis for a point $(\lambda=1.5$, $\mu=3)$
on the tilted critical line with with the point being approached from the
localized side of the phase diagram. Both for Figs. (a-b), we see the best
scaling collapse for $\nu=1$. For Fig. (c), we see the best scaling collapse
for ${\nu}=0.64$. For Fig.(d-e), the estimated values of ${\nu}$ are found to
be $\approx~{}0.6$ and $0.44$, respectively, while being approached from the
delocalized and localized sides of the phase diagram. Fig. (f): We observe the
best scaling collapse for ${\nu}=0.45$.

As per Eq. (32), one would expect that $S_{\infty}/S_{T}$ vs $\delta\tau
L^{1/\nu}$ to exhibit a scale collapse to a single curve according to Eq.
(32).

As shown in Fig. 13 (a)-13 (b), we first analyze the finite-size scaling of
the entanglement entropy for the AAH critical point and the vertical critical
line ($\lambda=0.5,~{}\mu=2$) of the static model, separating the extended and
localized phases. For both cases, we see the best scaling collapse for
$\nu=1$, which agrees with the earlier findings Sinha  _et al._ (2019). In
Fig. 13 (c), we repeat the same analysis for the horizontal critical line with
$\lambda=1,~{}\mu=1.5$, which separates the extended phase from the critical
one. For this case, we observe $\nu~{}\approx~{}0.64$ when being approached
from the delocalized side of the critical line. In Fig. 13 (d)- 13 (e), we do
the similar analysis for the triple point with $\lambda=1$ and $\mu=2$ when
being approached from the delocalized and localized sides of phase diagram,
respectively. For these two cases, we see the best scaling collapse for
$\nu\approx 0.6$ and $\nu\approx 0.44$, respectively. In Fig. 13 (f), we
perform the same analysis for the tilted critical line with $\lambda=1.5$ and
$\mu=3$ while approaching from the localized side of the critical line, and
the corresponding value of $\nu\approx~{}0.45$. Interestingly, the triple
point shows two different values of the critical exponent $\nu$ when
approached from the delocalized phase and localized phase, respectively. The
new kind of criticality of the quantum triple point Belitz and Kirkpatrick
(2017) demands a more careful investigation, which goes beyond the scope of
current work. Scaling collapse in the critical phase is not possible due to
strong fluctuations in the data, which presumably originates due to
characteristic huge quantum fluctuations in the critical phase.

##  Appendix C The growth of the exact surface roughness operator at the AAH
critical point and inside the critical phase

Figure 14: The dynamics of the surface roughness operator at the AAH critical
point and inside the critical phase: (a-b) The variation of the exact surface
roughness operator $W(L,t)$ with time for $L=400$, $600$ and $800$. The
parameter values chosen for these plots are as follows: Fig. (a):
$\lambda=0,~{}\mu=2$ (AAH critical point), Fig. (b): $\lambda=1.2,~{}\mu=1.8$
(a point inside the critical phase). In the first cases, the dynamics show
diffusive growth in time with best fitted FV scaling exponents
$(\alpha,\beta,z)$=$(0.51,0.24,2.13)$. In Fig. (b), we see subdiffusive growth
in time with the best fitted FV scaling exponents
$(\alpha,\beta,z)=(0.53,0.19,2.79)$.

Here we show the dynamics of the exact surface roughness operator $W(L,t)$
with time for the AAH critical point and for the critical phase. In Fig. 14
(a-b), we consider the AAH critical point $(\lambda,\mu)=(0,2)$ of the static
model and a point inside the critical phase of the static model with
$\lambda=1.2$ and $\mu=1.8$, respectively. In Fig. 14 (a), we see that the
surface roughness exhibits almost the diffusive dynamical behavior with the
best fitted FV scaling exponents $(\alpha,\beta,z)=(0.51,0.24,2.13)$, which
agrees with the earlier findings Han  _et al._ (1994). On the other hand, the
critical phase demonstrates subdiffusive dynamical behavior with the scaling
exponents $(\alpha,\beta,z)=(0.53,0.19,2.79)$, as shown in Fig. 14(b).

##  Appendix D Fractal dimension of the single-particle Floquet eigenstates of
the periodically driven model

In this section, we will analyze the extraction procedure of the fractal
dimension from the scaling of the inverse participation ratio (IPR) for the
single-particle Floquet eigenstates using the numerical fitting procedure. To
determine if any of the Floquet eigenstates exhibit multifractal behavior, we
calculate IPR, also defined in the main text as

$I_{m}^{(2)}~{}=~{}\sum_{j}~{}|\psi_{m}(j)|^{4},$ (33)

and study its scaling with the system size $L$ Castellani and Peliti (1986);
Evers and Mirlin (2008b). If $I_{m}^{(2)}$ scales as $L^{-\tau_{2}}$, then
$\tau_{2}=1$ for extended states (since $|\psi_{m}(j)|^{2}$ is of the order of
$1/L$ for all $j$ for such states), and $\tau_{2}=0$ for localized states
(since $|\psi_{m}(j)|^{2}$ for such states is of order 1 over a finite region
whose size remains constant as $L\to\infty$). Multifractal or critical states
typically have $0<\tau_{2}<1$. In Figs. 15 (a-c), we consider three single-
particle Floquet eigenstates for three different parameter values, namely, (a)
$J=1$, $\lambda=0.5$, $\mu=2$, $\delta\mu=1.5$, and $\omega=0.005$ (vertical
critical line of the static model subjected to periodic driving in the slow
driving limit), (b-c) $J=1$, $\lambda=1.5$, $\mu=3$, $\delta\mu=2.5$,
$\omega=30$, and $0.5$, respectively (tilted critical line of the static model
subjected to periodic driving in the fast and slow driving regime,
respectively). We then extract the scaling exponent of IPR ($\tau_{2}$) by
performing a linear fit of $log(I_{m}^{(2)})$ vs $log(L)$ plot. As shown in
Fig. 15 (a-c), the single-particle Floquet eigenstate show extended,
(multi)fractal and localized behavior with $\tau_{2}=0.99$, $0.62$, and $0$,
respectively, which have been extracted from the linear fitting analysis. This
is how the dynamical phases, as shown in Figure 3 of the main text, has been
determined.

Figure 15:  The scaling exponents of the inverse participation ratios of the
single-particle spectrum-resolved Floquet eigenstates: (a-c) The inverse
participation ratio ($I_{m}^{(2)}$) of the single-particle spectrum-resolved
Floquet eigenstates are plotted for three different system sizes,
$L=500,~{}1000$ and $2000$, $3000$ and $4000$ in log-log scale. The scaling
exponent ($\tau_{2}$) can be easily extracted from the slope of
$log(I_{m}^{2})$ vs $log(L)$ plot by a linear fitting analysis. Fig. (a): In
this case, we consider a single-particle Floquet eigenstate with
$m/L=0.5$$J=1$, $\lambda=0.5$, $\mu=2$, $\delta\mu=1.5$, and $\omega=0.005$.
Fig. (b): Here we consider $J=1$, $\lambda=1.5$, $\mu=3$, $\delta\mu=2.5$, and
$\omega=30$. The same analysis as Fig. (a) is performed for the single-
particle Floquet eigenstate with $m/L=0.7$. Fig. (c): For this particular
case, we consider $J=1$, $\lambda=1.5$, $\mu=3$, $\delta\mu=2.5$,
$\omega=0.5$, and the single-particle Floquet eigenstate with $m/L=0.5$. The
values of $\tau_{2}$ extracted for these cases obtained from the numerical
fitting analysis are $1~{}(\text{extended}),~{}0.62~{}(\text{critical})$, and
$0.00~{}(\text{localized})$, respectively.

##  Appendix E  More details on the dynamics of the driven model

###  E.1 Equivalence between the approximate surface-roughness operator and
half-chain entanglement entropy in the case of periodically driven model

Just like the static case, the equivalence between approximate surface
roughness operator, $W_{a}^{2}(t)$ and von Neumann entanglement entropy,
$S_{L/2}(t)$ Fujimoto  _et al._ (2021) holds even for the periodically driven
model in the case of stroboscopic time evolution. To show that, in Figs. 16
(a-d), we examine the stroboscopic time evolution of $h_{\rm{av}}(t)$ and
$\sum_{j=1}^{L/2}{\rm Tr}(\rho(t)n_{j})$ for the parameter values as follows:
Figs. 16 (a) - 16 (c): $\lambda=1,~{}\mu=1.5,~{}\delta\lambda=0.8$ (driving
the horizontal critical line of the static model along the vertical direction
with finite $\delta\lambda$), and Figs. 16 (b) - 16 (d):
$\lambda=1.5,~{}\mu=3,~{}\delta\mu=2.5$ (driving the tilted critical line of
the static model along the horizontal direction with finite $\delta\mu$). In
both cases, we consider periodic driving ranging from slow to fast driving
limit with $\omega=1,~{}5$ and $10$ as for examples. During the stroboscopic
time evolution, $h_{\rm{av}}(t)$, as shown in Figs. 16 (a-b), and
$\sum_{j=1}^{L/2}{\rm Tr}(\rho(t)n_{j})$, as shown in Figs. 16 (c-d), remain
close to zero and $\frac{\nu L}{2}$ ($\nu$ being the filling fraction of the
many-body initial state) similar to the static case, revealing that
stroboscopic Floquet dynamics satisfy the conditions (i) and (iii) as
mentioned earlier in the static case. These conditions lead to the equivalence
between the square of the approximate surface roughness and von Neumann
entanglement entropy even during the stroboscopic dynamics.

(a) (b) (c) (d)

Figure 16: Stroboscopic time evolution of $\mathbf{h_{\rm{av}}(t)}$ and
$\mathbf{\sum_{j=1}^{L/2}{\rm Tr}}\mathbf{(\rho(t)n_{j})}$ of the periodically
driven model: (a-b) The stroboscopic evolution of $h_{\rm{av}}(t)$ vs $t$ for
the parameter values, Fig. (a): $\lambda=1,~{}\mu=1.5,~{}\delta\lambda=0.8$,
Fig. (b): $\lambda=1.5,~{}\mu=3,~{}\delta\mu=2.5$. For both cases, we consider
slow to fast driving limits with $\omega=1,~{}5,~{}10$, and for $L=200$. In
both cases, $h_{\rm{av}}(t)$ remains close to zero for all time. (c-d) The
stroboscopic time evolution of $\sum_{j=1}^{L/2}{\rm Tr}(\rho(t)n_{j})$ are
shown with time for the parameter values same as Figs. (a-b). In both cases ,
$\sum_{j=1}^{L/2}{\rm Tr}(\rho(t)n_{j})=\frac{\nu L}{2}\simeq 50$ during the
stroboscopic dynamics just as the static model, which is expected since
$\nu=1/2$ and $L=200$ for these particular cases.

###  E.2 Localized states with no well-defined FV scaling exponents

Here we will discuss the FV scaling in the case of drive-induced localized
states. As shown in Figs. 17 (a-b), the tilted critical line of the static
model ($\lambda=1.5,~{}\mu=3$) is subjected to inter-phase driving (critical
and localized phases) with $(\delta\lambda,\delta\mu)=(0,2.5)$ in the
comparatively slow driving regime ($\omega=1$ and $\omega=5$). In Fig. 17 (c),
we drive the horizontal critical line in the identical manner with
$(\delta\lambda,\delta\mu,\omega)=(0,1.2,1)$. In all three cases, we see no
well-defined scaling exponents since the surface roughness saturates to a
value that does not scale up with the system size properly. Hence, here FV
scaling is absent for localized phases.

Figure 17: Absence of FV scaling in localized phase of the periodically driven
model: (a-b) The stroboscopic variation of approximate $W^{2}(L,t)$ with time
for the tilted critical line ($\lambda=1.5,~{}\mu=3.0$) driven with
$\delta\lambda=0$ and $\delta\mu=2.5$ in comparatively slow driving regime
($\omega=1$ and $\omega=5$, respectively). (c) The same variation for the
horizontal critical line ($\lambda=1$, $\mu=1.5$) subjected to periodic
driving with $(\delta\lambda,\delta\mu,\omega)=(0,1.2,1)$. In all three cases,
we see localized behavior with no well-defined FV scaling behavior.

Figure 18: Dynamics of driven systems in extremely fast driving regime: (a-c)
The stroboscopic time evolution of approximate $W^{2}(L,t)$ vs time where the
horizontal, tilted and vertical critical lines are subjected to periodic
driving with Fig. (a): $\lambda=1$, $\mu=1.5$, $\delta\lambda=0.8$,
$\delta\mu=0$, and $\omega=20$, Fig. (b): $\lambda=1.5$, $\mu=3$,
$\delta\lambda=0$, $\delta\mu=2.5$, and $\omega=20$, Fig. (c): $\lambda=0.5$,
$\mu=2$, $\delta\lambda=0$, $\delta\mu=1.5$, and $\omega=20$, respectively.
The best fitted FV scaling exponents, $(2\alpha_{2},2\beta_{2},z_{2})$ for
these cases are found to be: Fig. (a): $(1.08,0.32,3.34)$, Fig. (b):
$(0.92,~{}0.31,~{}3.00)$, Fig. (c) $(1.07,~{}0.45,~{}2.38)$, which are quite
close to the static case.

###  E.3 The surface-roughness dynamics in the extremely fast driving regime

The dynamics of the periodically driven model in the extremely high frequency
regime converge to the static model, which can be seen in Fig. 18. As shown in
Figs. 18 (a-c), the horizontal ($\lambda=1,~{}\mu=1.5$), tilted
($\lambda=1.5,~{}\mu=3$) and vertical critical lines of the static model are
subjected to periodic driving in the high-frequency regime with the parameter
values, $(\delta\lambda,\delta\mu,\omega)=(0.8,0,20)$, $(0,2.5,20)$, and
$(0,1.5,20)$, respectively. In first two cases, the square of the approximate
surface roughness operator show subdiffusive dynamics with the best fitted FV
scaling exponents
$(2\alpha_{2},~{}2\beta_{2},~{}z_{2})$=$(1.08,~{}0.32,~{}3.34)$ and
$(0.92,0.31,3.00)$, respectively, while the last one demonstrates almost
diffusive dynamics with FV scaling exponents
$(2\alpha_{2},~{}2\beta_{2},~{}z)=(1.07,0.45,2.38)$. These behaviors are found
to be the same as observed in the static model, which one can anticipate due
to reduced amount of drive-induced phase mixing in the high-frequency limit.

##  Appendix F Entanglement entropy of the driven model using the correlation
matrix method

In this section, we will briefly discuss the formulation of the correlation
matrix method to analyze the stroboscopic dynamics of a generic non-
interacting Floquet system. Before proceeding further, we first briefly
discuss how to compute the dynamics of the entanglement entropy of a static
model with a non-interacting Hamiltonian $H$ starting from an many-body
initial state $|\psi_{\rm{in}}\rangle$, which is not the many-body eigenstate
of $H$. To do so, let us first consider a generic non-interacting Hamiltonian
given by

$\displaystyle H=\sum_{i\neq j}t_{ij}c_{i}^{\dagger}c_{j}+\rm{H.c.},$ (34)

where $c_{i}$ and $c_{i}^{\dagger}$ are the fermionic annihilation and
creation operator, respectively. Followed by the diagonalization, such a
Hamiltonian takes the following form
$H=\sum_{k=1}^{N}\epsilon_{k}\tilde{c}^{\dagger}_{k}\tilde{c}_{k}$, where
$\tilde{c}_{k}=\sum_{j=1}^{N}\psi_{k}(j)c_{j}$. Setting $\hbar=1$, the time
evolution of the annihilation operator can be found using the Heisenberg
equation of motion
$\frac{d\tilde{c}_{k}(t)}{dt}=\frac{1}{i}[\tilde{c}_{k},H]=\frac{1}{i}\epsilon_{k}\tilde{c}_{k}$,
which implies $\tilde{c}_{k}=e^{i\epsilon_{k}t}\tilde{c}_{k}(0)$. To study the
dynamics, we consider a many-body density-wave (DW) type of initial state
given below

$\displaystyle|\psi_{\rm{in}}\rangle=c_{2}^{\dagger}c_{4}^{\dagger}......c_{N}^{\dagger}|0\rangle,$
(35)

where $N$ being the system size with an even number of lattice sites, and we
further assume that the system is at half-filling, i.e., $N_{p}=N/2$. Next, we
will construct the correlation matrix Peschel (2003, 2012); Peschel and Eisler
(2009); Roy and Sharma (2018) in a sub-system (A) of size $L$ for this initial
many-body quantum state, which is,
$C_{ij}(t)=\langle\psi_{\rm{in}}|c_{i}^{\dagger}c_{j}|\psi_{\rm{in}}\rangle$,
where $i,~{}j\in A$. Following the usual procedure Peschel (2003, 2012);
Peschel and Eisler (2009); Roy and Sharma (2018), one can easily deduce that
the elements of the correlation matrix of the full system takes the following
form

$\displaystyle\langle\psi_{\rm{in}}|c_{i}^{\dagger}c_{j}|\psi_{\rm{in}}\rangle=\sum_{k,k^{\prime}=1}^{N}\sum_{i^{\prime}=2,4,....}$
$\displaystyle\psi_{i}(k)\psi_{i^{\prime}}^{\ast}(k^{\prime})\psi_{i^{\prime}}(k)\psi_{i^{\prime}}^{\ast}(k^{\prime})$
(36) $\displaystyle~{}\times e^{-i(\epsilon_{k}-\epsilon_{k^{\prime}})t}$

The entanglement entropy $S_{A}(t)$ followed by the diagonalization of a
$L\times L$ block of the correlation matrix associated with the sub-system A
is given by

$\displaystyle
S_{A}(t)=-\sum_{m=1}^{L}\left[\lambda_{m}\ln\lambda_{m}+(1-\lambda_{m})\ln(1-\lambda_{m})\right],$

where $\lambda_{m}$ stands for the $m$-the eigenvalue of the $L\times L$ block
of the correlation matrix.

Now, we will generalize this procedure for the stroboscopic closed Floquet
dynamics. For a generic periodically driven system, one can define an
effective time-independent Floquet Hamiltonian,
$H_{F}=\frac{i}{T}\ln(U_{F}(T))$ as long as the Floquet eigenvalues,
$e^{-i\theta_{m}}$ satisfy the condition, $|\theta_{m}|<<\pi$ for all the
Floquet eigenstates, where $T$ and $U_{F}$ being the driving period and the
Floquet evolution operator, respectively. Therefore, in this case, the static
formalism will hold with a slight modification, i.e., the Hamiltonian of a
static model, $H$ will get replaced by the effective Floquet Hamiltonian,
$H_{F}$ as far as the stroboscopic Floquet dynamics is concerned.
Interestingly, the correlation matrix method holds even if the condition,
$|\theta_{m}|<<\pi$ is violated (mostly in extremely slow driving regime,
which is closely connected to the adiabetic limit), and the reason for this is
as follows. Note that in the last line of Eq. 36, the eigenvalues of a generic
static Hamiltonian appear only in the exponential factors, i.e.,
$e^{-i\epsilon_{k}t}$ factors, which turns out to be crucial for its
generalization to the stroboscopic time evolution of a generic closed non-
interacting Floquet quantum system. In this scenario, although $H_{F}$ is not
uniquely defined, the Floquet eigenvalues, $e^{-i\epsilon_{m}T}$ (also well-
defined for all driving frequencies) still can be defined by diagonalizing the
Floquet evolution operator, $U_{F}$ (setting $\hbar=1$),

$\displaystyle U_{F}|m\rangle=e^{-i\epsilon_{m}T}|m\rangle,$ (38)

where $|m\rangle$ denotes the $m$-th single-particle Floquet eigenfunction.
Substituting the Floquet eigenvalues and single-particle Floquet eigenstates
obtained from Eq. (38) in Eq. (36), one can easily generalize the correlation
matrix formalism to analyze the stroboscopic Floquet dynamics even if the
effective Floquet Hamiltonian is ambiguous.

## References

  * Polkovnikov  _et al._ (2011) A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Reviews of Modern Physics 83, 863 (2011).
  * Langen  _et al._ (2015) T. Langen, R. Geiger, and J. Schmiedmayer, Annual Review of Condensed Matter Physics 6, 201 (2015).
  * Oka and Kitamura (2019) T. Oka and S. Kitamura, Annual Review of Condensed Matter Physics 10, 387 (2019).
  * Fujimoto  _et al._ (2020) K. Fujimoto, R. Hamazaki, and Y. Kawaguchi, Phys. Rev. Lett. 124, 210604 (2020).
  * Fujimoto  _et al._ (2021) K. Fujimoto, R. Hamazaki, and Y. Kawaguchi, Physical Review Letters 127, 090601 (2021).
  * Fujimoto  _et al._ (2022) K. Fujimoto, R. Hamazaki, and Y. Kawaguchi, Phys. Rev. Lett. 129, 110403 (2022).
  * Vicsek and Family (1984) T. Vicsek and F. Family, Phys. Rev. Lett. 52, 1669 (1984).
  * Family and Vicsek (1985) F. Family and T. Vicsek, Journal of Physics A: Mathematical and General 18, L75 (1985).
  * Barabási and Stanley (1995) A.-L. Barabási and H. E. Stanley,  _Fractal Concepts in Surface Growth_ (Cambridge University Press, 1995).
  * Jin  _et al._ (2020) T. Jin, A. Krajenbrink, and D. Bernard, Physical Review Letters 125, 040603 (2020).
  * Edwards and Wilkinson (1982) S. F. Edwards and D. R. Wilkinson, Proc. R. Soc. Lond. A 381 (1982).
  * Kardar  _et al._ (1986) M. Kardar, G. Parisi, and Y.-C. Zhang, Phys. Rev. Lett. 56, 889 (1986).
  * Ruidas and Banerjee (2021) S. Ruidas and S. Banerjee, SciPost Phys. 11, 087 (2021).
  * Ljubotina  _et al._ (2019) M. Ljubotina, M. Žnidarič, and T. c. v. Prosen, Phys. Rev. Lett. 122, 210602 (2019).
  * Ilievski  _et al._ (2018) E. Ilievski, J. De Nardis, M. Medenjak, and T. Prosen, Phys. Rev. Lett. 121, 230602 (2018).
  * Ye  _et al._ (2022) B. Ye, F. Machado, J. Kemp, R. B. Hutson, and N. Y. Yao, Phys. Rev. Lett. 129, 230602 (2022).
  * Keenan  _et al._ (2022) N. Keenan, N. Robertson, T. Murphy, S. Zhuk, and J. Goold, arXiv:2208.12243 (2022).
  * Wei  _et al._ (2022) D. Wei, A. Rubio-Abadal, B. Ye, F. Machado, J. Kemp, K. Srakaew, S. Hollerith, J. Rui, S. Gopalakrishnan, N. Y. Yao, I. Bloch, and J. Zeiher, Science 376, 716 (2022).
  * Roy and Sharma (2021) N. Roy and A. Sharma, Journal of Physics: Condensed Matter 33, 334001 (2021).
  * Ahmed  _et al._ (2021) A. Ahmed, N. Roy, and A. Sharma, Phys. Rev. B 104, 155137 (2021).
  * Purkayastha  _et al._ (2017) A. Purkayastha, A. Dhar, and M. Kulkarni, Phys. Rev. B 96, 180204 (2017).
  * Varma  _et al._ (2017) V. K. Varma, C. de Mulatier, and M. Žnidarič, Phys. Rev. E 96, 032130 (2017).
  * Iglói  _et al._ (2013) F. Iglói, G. Roósz, and Y.-C. Lin, New Journal of Physics 15, 023036 (2013).
  * Sarkar  _et al._ (2021) M. Sarkar, R. Ghosh, A. Sen, and K. Sengupta, Phys. Rev. B 103, 184309 (2021).
  * Aditya  _et al._ (2023) S. Aditya, K. Sengupta, and D. Sen, Physical Review B 107, 035402 (2023).
  * Roy  _et al._ (2018) S. Roy, I. M. Khaymovich, A. Das, and R. Moessner, SciPost Phys. 4, 025 (2018).
  * Shimasaki  _et al._ (2022) T. Shimasaki, M. Prichard, H. E. Kondakci, J. Pagett, Y. Bai, P. Dotti, A. Cao, T.-C. Lu, T. Grover, and D. M. Weld, arXiv: 2203.09442 (2022).
  * Ray  _et al._ (2018) S. Ray, A. Ghosh, and S. Sinha, Phys. Rev. E 97, 010101 (2018).
  * Schreiber  _et al._ (2015) M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lüschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch, Science 349, 842 (2015).
  * Eckardt (2017) A. Eckardt, Rev. Mod. Phys. 89, 011004 (2017).
  * Messer  _et al._ (2018) M. Messer, K. Sandholzer, F. Görg, J. Minguzzi, R. Desbuquois, and T. Esslinger, Phys. Rev. Lett. 121, 233603 (2018).
  * Chaikin and Lubensky (1995) P. M. Chaikin and T. C. Lubensky,  _Principles of Condensed Matter Physics_ (Cambridge University Press, 1995).
  * Han  _et al._ (1994) J. H. Han, D. J. Thouless, H. Hiramoto, and M. Kohmoto, Phys. Rev. B 50, 11365 (1994).
  * Aoki (1986) H. Aoki, Phys. Rev. B 33, 7310 (1986).
  * Mirlin  _et al._ (2006) A. D. Mirlin, Y. V. Fyodorov, A. Mildenberger, and F. Evers, Phys. Rev. Lett. 97, 046803 (2006).
  * Janssen (1998) M. Janssen, Physics Reports 295, 1 (1998).
  * Evers and Mirlin (2008a) F. Evers and A. D. Mirlin, Rev. Mod. Phys. 80, 1355 (2008a).
  * Mikami  _et al._ (2016) T. Mikami, S. Kitamura, K. Yasuda, N. Tsuji, T. Oka, and H. Aoki, Phys. Rev. B 93, 144307 (2016).
  * Aubry and André (1980) S. Aubry and G. André, Ann. Israel Phys. Soc 3, 18 (1980).
  * Sokoloff (1985) J. Sokoloff, Physics Reports 126, 189 (1985).
  * Spohn (2014) H. Spohn, Journal of Statistical Physics 154, 1191 (2014).
  * Mendl and Spohn (2015) C. B. Mendl and H. Spohn, Journal of Statistical Mechanics: Theory and Experiment 2015, P08028 (2015).
  * Kulkarni  _et al._ (2015) M. Kulkarni, D. A. Huse, and H. Spohn, Phys. Rev. A 92, 043612 (2015).
  * Antal  _et al._ (2008) T. Antal, P. L. Krapivsky, and A. Rákos, Phys. Rev. E 78, 061115 (2008).
  * Eisler and Rácz (2013) V. Eisler and Z. Rácz, Phys. Rev. Lett. 110, 060602 (2013).
  * Hunyadi  _et al._ (2004) V. Hunyadi, Z. Rácz, and L. Sasvári, Phys. Rev. E 69, 066103 (2004).
  * Antal and Rácz (1996) T. Antal and Z. Rácz, Phys. Rev. E 54, 2256 (1996).
  * Bhakuni and Lev (2023) D. S. Bhakuni and Y. B. Lev, arXiv: 2309.03273 (2023).
  * Blanes  _et al._ (2009) S. Blanes, F. Casas, J. Oteo, and J. Ros, Physics Reports 470, 151 (2009).
  * Sen  _et al._ (2021) A. Sen, D. Sen, and K. Sengupta, Journal of Physics: Condensed Matter 33, 443003 (2021).
  * Bukov  _et al._ (2015) M. Bukov, L. D’Alessio, and A. Polkovnikov, Advances in Physics 64, 139 (2015).
  * Evers and Mirlin (2008b) F. Evers and A. D. Mirlin, Rev. Mod. Phys. 80, 1355 (2008b).
  * Castellani and Peliti (1986) C. Castellani and L. Peliti, Journal of Physics A: Mathematical and General 19, L429 (1986).
  * Nair  _et al._ (2023) P. S. Nair, D. Joy, and S. Sanyal, arXiv: 2302.14053 (2023).
  * Lezama  _et al._ (2019) T. L. M. Lezama, S. Bera, and J. H. Bardarson, Phys. Rev. B 99, 161106 (2019).
  * Soori and Sen (2010) A. Soori and D. Sen, Phys. Rev. B 82, 115432 (2010).
  * Takeuchi  _et al._ (2011) A. K. Takeuchi, M. Sano, T. Sasamoto, and H. Spohn, Scientific Reports 1, 34 (2011).
  * (58) S. Mu, J. Gong, and G. Lemarié, arXiv:2306.08272 (2023) .
  * Li  _et al._ (2023) H. Li, Y.-Y. Wang, Y.-H. Shi, K. Huang, X. Song, G.-H. Liang, Z.-Y. Mei, B. Zhou, H. Zhang, J.-C. Zhang, _et al._ , npj Quantum Information 9, 40 (2023).
  * Xiao  _et al._ (2021) T. Xiao, D. Xie, Z. Dong, T. Chen, W. Yi, and B. Yan, Science Bulletin 66, 2175 (2021).
  * Basko  _et al._ (2006) D. Basko, I. Aleiner, and B. Altshuler, Annals of Physics 321, 1126 (2006).
  * Vosk and Altman (2013) R. Vosk and E. Altman, Phys. Rev. Lett. 110, 067204 (2013).
  * Huse  _et al._ (2013) D. A. Huse, R. Nandkishore, V. Oganesyan, A. Pal, and S. L. Sondhi, Phys. Rev. B 88, 014206 (2013).
  * Oganesyan and Huse (2007) V. Oganesyan and D. A. Huse, Phys. Rev. B 75, 155111 (2007).
  * Pal and Huse (2010) A. Pal and D. A. Huse, Phys. Rev. B 82, 174411 (2010).
  * Ponte  _et al._ (2015) P. Ponte, Z. Papić, F. m. c. Huveneers, and D. A. Abanin, Physical Review letters 114, 140401 (2015).
  * Bordia  _et al._ (2017) P. Bordia, H. Lüschen, U. Schneider, M. Knap, and I. Bloch, Nature physics 13, 460 (2017).
  * Sierant  _et al._ (2023) P. Sierant, M. Lewenstein, A. Scardicchio, and J. Zakrzewski, Phys. Rev. B 107, 115132 (2023).
  * Decker  _et al._ (2020) K. S. C. Decker, C. Karrasch, J. Eisert, and D. M. Kennes, Phys. Rev. Lett. 124, 190601 (2020).
  * Khemani  _et al._ (2016) V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. 116, 250401 (2016).
  * Else  _et al._ (2016) D. V. Else, B. Bauer, and C. Nayak, Phys. Rev. Lett. 117, 090402 (2016).
  * Zhang  _et al._ (2017) J. Zhang, P. W. Hess, A. Kyprianidis, P. Becker, A. Lee, J. Smith, G. Pagano, I.-D. Potirniche, A. C. Potter, A. Vishwanath, N. Y. Yao, and C. Monroe, Nature 543, 217 (2017).
  * Wang  _et al._ (2021) Y. Wang, C. Cheng, X.-J. Liu, and D. Yu, Phys. Rev. Lett. 126, 080602 (2021).
  * Sinha  _et al._ (2019) A. Sinha, M. M. Rams, and J. Dziarmaga, Phys. Rev. B 99, 094203 (2019).
  * Belitz and Kirkpatrick (2017) D. Belitz and T. Kirkpatrick, Physical review letters 119, 267202 (2017).
  * Peschel (2003) I. Peschel, Journal of Physics A: Mathematical and General 36, L205 (2003).
  * Peschel (2012) I. Peschel, Journal of Physics A: Mathematical and General 42, 267 (2012).
  * Peschel and Eisler (2009) I. Peschel and V. Eisler, Journal of Physics A: Mathematical and Theoretical 42, 504003 (2009).
  * Roy and Sharma (2018) N. Roy and A. Sharma, Phys. Rev. B 97, 125116 (2018).