Multifractality without fine-tuning in a Floquet quasiperiodic chain

Periodically driven, or Floquet, disordered quantum systems have generated many unexpected discoveries of late, such as the anomalous Floquet Anderson insulator and the discrete time crystal. Here, we report the emergence of an entire band of multifractal wavefunctions in a periodically driven chain of non-interacting particles subject to spatially quasiperiodic disorder. Remarkably, this multifractality is robust in that it does not require any fine-tuning of the model parameters, which sets it apart from the known multifractality of $critical$ wavefunctions. The multifractality arises as the periodic drive hybridises the localised and delocalised sectors of the undriven spectrum. We account for this phenomenon in a simple random matrix based theory. Finally, we discuss dynamical signatures of the multifractal states, which should betray their presence in cold atom experiments. Such a simple yet robust realisation of multifractality could advance this so far elusive phenomenon towards applications, such as the proposed disorder-induced enhancement of a superfluid transition.


Introduction
Multifractal wavefunctions are beautifully complex states, extended yet non-ergodic, comprising both rare high peaks and long polynomial tails of wavefunction amplitudes. The physics of multifractality is commonly associated with critical wavefunctions at Anderson localisationdelocalisation and quantum Hall plateau transitions [1][2][3]. Multifractality also appears in hierarchical and infinite-dimensional systems like random regular graphs, Bethe lattices, and more generally in fully connected random matrix ensembles and network models [4][5][6][7][8][9][10][11][12]. The presence of long-ranged physics diagnosed via correlation and localisation lengths unifies these two contexts. Hence, realising multifractality in an inherently short-ranged system, specifically systems with short-ranged hoppings and interactions unlike those for example, represented by power law banded or infinite ranged random matrices [5,9], without fine-tuning to criticality poses not only an interesting and important theoretical challenge but is also desirable for a robust experimental realisation of multifractality and consequent applications.
We find that multifractality in a short-ranged system requires only relatively simple ingredients, namely a time-periodic modulation of a spatially-quasiperiodic system possessing a single particle mobility edge. Periodically driven systems, also known as, Floquet systems [13] have witnessed much interest recently with significant advances [14] in the understanding of their statistical mechanics [15][16][17] and phase structures [18] and in their experimental realisations with cold atoms [19]. Technically, the eigenfunctions of the Floquet unitary timeevolution operator over one period, U , encode the full information about the stroboscopic dynamics of the system, much like the eigenfunctions of the Hamiltonian of a static system [14]. At the same time, it has also been realised that single particle mobility edges occur naturally in simple incommensurate bichromatic potentials [20,21]. At the level of onedimensional lattice systems, this is related to the single particle mobility edges that generally exist in deformations of the Aubry-André model [22][23][24][25][26][27][28].
The central finding of this work is that, when the periodic drive hybridises the localised and delocalised states on either side of a mobility edge in a one-dimensional system with quasiperiodic potential, it gives rise to a band of multifractal eigenstates of the corresponding Floquet operator U . Remarkably, this multifractality exists in a finite range of parameters and thence requires no fine-tuning, while the states nonetheless show anomalous algebraic multifractal correlations similar-in some but not all-respects to the critical ones [29][30][31][32][33]. We present an effective random matrix Hamiltonian, which captures the numerically obtained multifractality remarkably well, bearing a family resemblance to the Rosenzweig-Porter random matrix ensemble [34], generalisations of which are known to host multifractal eigenstates [9,11,[35][36][37][38][39][40].  (1), where the colour shows the scaling of the inverse participation ratio with system size, with green corresponding to the delocalised states (∼ L −1 ) and blue, localised states (∼ L 0 ). The red line denotes the mobility edge. (b) The energy levels corresponding to V 0 (denoted by the black dashed line in (a)) where a periodic drive with frequency Ω chosen to be slightly smaller than the bandwidth couples the delocalised and localised states approximately within the gray shaded windows.

Model and numerical results
Our starting point is a variant of the one-dimensional Aubry-André Hamiltonian It comprises a simple nearest-neighbour hopping term alongside a potential quasiperiodic on account of its incommensurate wavevector, which we set to the golden mean, κ = ( √ 5 + 1)/2. The model exhibits a mobility edge [28] at an energy, ε ME = 2sgn(V )(|J| − |V |/2)/µ as shown in Fig. 1(a). We set J = 1 and µ = −0.6 throughout. Eigenstates with energies above and below ε ME are completely delocalised and exponentially localised, respectively. In all numerical analysis, we average the data over various values of the θ which is analogous to disorder averaging.
The system is driven by a time-periodic modulation of the amplitude of the quasiperiodic potential in the form of a square wave with frequency Ω, mean V 0 , and amplitude ∆V . Such a protocol allows for the exact computation of the Floquet eigenstates, denoted henceforth as |φ via numerical diagonalisation of the U which can be calculated relatively straightforwardly as U = e −iH + π/Ω e −iH − π/Ω where H ± denotes the Hamiltonian in the two steps of the square wave.
A common diagnostic for localisation properties of wavefunctions is their inverse participation ratio, IPR = I 2 = x |φ(x)| 4 which scales with system size L as L −1 (L 0 ) for delocalised(localised) states in one dimension. The first signs of Floquet multifractality appear in the scalings of IPRs of the Floquet eigenstates. As Ω is chosen to be slightly smaller the bandwidth of the spectrum of the static Hamiltonian (1), with V = V 0 (see Fig. 1), the drive primarily couples states close to the top and bottom of the undriven spectrum, leaving largely unaffected all the localised and delocalised states in between. These latter two, together with The inverse participation ratio (IPR), shown for the Floquet eigenstates sorted in increasing order of I 2 , for different system sizes L. The collapse of different segments of the data for different L, when the IPR is scaled with (a)L 0 , (b)L 1/2 , and (c)L reflects the presence of localised, multifractal, delocalised states, respectively. The data in red in (a) shows the IPRs of static eigenstates (for L = 8192) for reference. (d) Averaging the moments over Floquet eigenstates in different windows highlighted by the vertical shaded regions in (b), τ q is plotted as a function of q, where the colour corresponds to the window. While the localised (blue) and delocalised (yellow) states show the expected standard behaviour, the multifractal states (shades of green) have τ q ≈ D(q − 1) at q 1 with D 1/2 (red dashed line). When averaged over all the multifractal states, the IPRs scale as L −τ 2 , where τ 2 = 0.55 ± 0.04 close to D 1/2. The τ q s are extracted as the slope of a linear fit of log I q versus log L. Representative fits are shown in Fig. 3. (e) The corresponding spectrum of fractal dimensions, f (α) as function of α clearly shows the multifractal states distinct from both localised and delocalised cases. The system parameters are V 0 = 2J, ∆V = J/2, and Ω = 2.74πJ, where the bandwidth of the undriven spectrum is ≈ 2.76πJ (see Appendix B for effects of lower frequencies).
our newly discovered multifractal states, are evident in Fig. 2(a)-(c), which now shows three distinct scalings of the IPR. In disordered systems, since the energy spectrum varies across disorder realisations, labelling the Floquet eigenstates in increasing order of their IPRs as in Fig. 2 turns out to be rather convenient. However, we also study the quasienergy resolved IPRs by appropriately binning the data (see Appendix A).
A more complete characterisation of multifractality is via a generalised IPR and its scaling exponent τ q , where D q = τ q /(q − 1) is known as the fractal dimension. For delocalised and localised states, D q = 1 and 0, respectively (for q > 0), whereas any other behaviour of D q implies multifractality. Multifractality is thus evidenced in τ q shown in Fig. 2(d), where the localised and delocalised states show their standard behaviour. The multifractal states on the other hand seem to have show a good agreement with τ q ≈ D(q − 1) with the q-independent D q = D ≈ 1/2, although one must notice that there is a spread in the behaviour of τ q across the window of all multifractal states. When averaged over all the multifractal states, D q turns out to be 0.55 ± 0.04 for q 1.
To extract τ q shown in Fig. 2(d), we do a linear fit of log I q versus log L for all values of q by selecting states in the shaded windows shown in Fig. 2(b). In Fig. 3, we show some examples of such fits for the delocalised and multifractal states for some representative values of q.
An equally fundamental measure of multifractality is the spectrum of fractal dimensions, f (α), which is defined via: the number of sites in a lattice system with total L sites where the wavefunction intensity |φ(x)| 2 ∼ L −α scales as L f (α) [1]. f (α) is a rather powerful measure as it formally contains the information of all the τ q s via a Legendre transform, f (α) = q * α − τ q * where q * is the solution of α = dτ q /dq. f (α) for the Floquet multifractal states is shown in Fig. 2(e) which is strikingly distinct from that of a localised (f (α) = lim αmax→∞ α/α max ) and delocalised (f (α) = 1 for α = 1 and −∞ otherwise) state. Figure 4: Power law decay of multifractal correlations. Multifractal spatial correlations C(r; p, q) in space shown as a function of the distance r for different values of (p, q), where the collapse of the data suggests an algebraic scaling form C(r; p, q) ∼ L −τ p+q −1 r −bp,q . The evident algebraic decay of the correlations is faster than for the previously studied critical multifractal states (red dashed lines). Turning to spatial correlations, multifractal wavefunctions exhibit algebraic behaviour concomitant with usual notions of criticality. To tease out multifractal behaviour, one again takes variable powers of the wavefunctions, to define the correlators C(r; p, q) = |φ(x)| 2p |φ(x+ r)| 2q x,disorder , which are averaged both spatially and over disorder realisations. These then have a scaling form ∼ L −ap,q r −bp,q , where a p,q = τ p+q + 1 as shown in Fig. 4. This is similar to known critical multifractal wavefunctions but in our case b p,q is larger than the reported value of τ p + τ q − τ p+q + 1 [1,30]. These states are thus genuinely fractal, not just mimicking fractality in their moments as in certain random-energy models [41].

Spectral decomposition of Floquet multifractal states
The underlying mechanism of Floquet generation of multifractality is the hybridisation of the localised and delocalised eigenstates of the undriven Hamiltonian close to the bottom and top of the spectrum, respectively. In this section, we provide evidence in form of the spectral decomposition of the Floquet eigenstates in terms of the those of the undriven Hamiltonian. To this order, we define the spectral density of states at energy ε for a Floquet eigenstate |φ as where ε ψ is an eigenvalue of the undriven Hamiltonian corresponding to the eigenstate |ψ , N is a normalization factor to ensure ρ( ) = 1 and η is a small broadening factor. As expected, ρ( ) for typical Floquet multifractal states chosen at random has overwhelming contributions Figure 6: Dependence of multifractality on driving frequency. Absence and presence of multifractal states in the Floquet spectrum when the frequency of the driving is larger and smaller than the bandwidth of the undriven spectrum, respectively. The two rows correspond to the two values of Ω, also shown schematically by the arrows next to the undriven spectrum, whereas the columns correspond to different scalings of IPR with L. The presence (absence) of multifractal states can be identified from the presence (absence) of finite fraction of states with IPR∼ L −1/2 in the second column.
from the localized and delocalized states near the bottom and top of the undriven spectrum as shown in Fig. 5.
To further corroborate this, we also show that the multifractal states appear only when the frequency of the driving, Ω is smaller than the bandwidth of the undriven spectrum. This is shown in Fig. 6 where the IPRs of the Floquet eigenstates show only two scalings ∼ L 0 and ∼ L −1 when Ω is larger than the bandwidth. On the contrary as soon as Ω is tuned below the bandwidth, the multifractal states show their presence which can be identified by the IPR of the finite fraction of Floquet eigenstates scaling approximately as L −1/2 .

Effective random matrix model
We now turn towards understanding the Floquet multifractality within a random matrix framework. A central ingredient is the coupling between localised states {|l }, mediated by the delocalised states, {|d }, to which the localised ones couple through the driving. That the Floquet drive strongly couples eigenstates of the undriven Hamiltonian which are resonant (separated in energy by the frequency of the drive), is elegantly represented in the so called Shirley picture [42], where the time-periodic problem is mapped onto a static problem of hopping on a ladder, whose legs are copies of the chain and whose 'transverse' coupling is provided by the time-periodic drive. With our Ω just below the bandwidth, resonant coupling occurs between states close to the edges of the undriven spectrum on the localised and delocalised sides in the undriven spectrum. This can be modelled by a two-leg truncation of the Shirley ladder as couplings to higher legs come with an energy denominator of the order of the bandwidth and are therefore parametrically suppressed.

Derivation of effective random matrix Hamiltonian
We start with deriving the offdiagonal matrix elements of the effective random matrix Hamiltonian.
According to the Bloch-Floquet theorem, the eigenstates of a time periodic Hamiltonian where ω is called the quasienergy and |φ(t) is itself periodic in time with frequency Ω. Expressing |φ(t) in terms of its Fourier components |φ(t) = n |φ n e inΩt , the Schrödinger equation for |φ n is given by where the H m denote Fourier components of H(t) = m H m e imΩt . One can choose to work in the eigenbasis of H 0 denoted by {|ψ } such that in which case, the Schrödinger equation can be recast as Since our driving frequency is only slightly smaller than the bandwidth, resonances can occur only between states near the top and the bottom of the static spectrum which are separated in energy by Ω. Hence, we employ a simple two-leg Shirley ladder to analyze the system, as further legs correspond to processes involving multiples of Ω, to which there are no corresponding resonances. Thus, we only keep H ±1 ∼ ∆V x v xnx and only the n = −1 and n = 0 sectors. Also, in our notation |ψ = x ψ(x)c † x |0 . Hence the matrix element ψ |H ±1 |ψ = ∆V x ψ * (x)ψ(x)v(x). These results can be put back in the equations for c n,ψ (t) as We know that the multifractal states come from the hybridisation of the localised ({|l }) and delocalised states ({|d }) near the bottom and the top of the spectrum, respectively. Hence, in the two-leg Shirley ladder [42], the only states with relevant contributions are the delocalised ones from the n = −1 sector and the localized ones from the n = 0 sector. This is schematically shown in Fig. 7 where the participating undriven states are marked with a gray shaded window.
Hence the coefficients of interest are c −1,d and c 0,l . The equation for c −1,d reads The gray shaded regions highlights the window of states which hybridize between the two legs when the driving is turned on.
Assuming that a localized state |l i is δ-function localized at x i , Eq. (9) can be plugged in Eq. (8) to obtain the equation for the coefficients of c 0,l i as where M l i l j denotes the off-diagonal matrix elements of the effective Hamiltonian, which is the leading effective matrix element determining a resulting Floquet eigenstate at quasienergy ω.
Here, a localised eigenstate of the undriven Hamiltonian |l i = |x i is assumed to be δfunction localised at x = x i , and the primed sum denotes a sum over resonant delocalised states, highlighted by the gray shaded window in Fig. 1(b). This leads to a fully connected random matrix Hamiltonian within the localised states, with the undriven eigenenergies on the diagonal, and the M l i l j as the off-diagonal matrix elements. As an aside, we note that this model formally resembles the Rosenzweig-Porter random matrix ensemble, unlike which, however, it has a probability distribution of M (Eq. (11)), denoted by P (M ) which is not Gaussian as we show in the following section 4.3.
The effective random matrix model allows us to connect the multifractality of the wavefunctions to the statistical properties of the Floquet-generated matrix elements M and their scalings with system size.

Perturbative calculation of wavefunction intensity distributions from effective random matrix
Similar to the analysis in Ref. [9], we treat the corrections to the δ-function localised eigenstates perturbatively φ l i (x j ) = M l i l j /(ε l i − ε l j ). From the statistics of the perturbed wavefunctions, f (α) can be extracted as we will show now. We derive the leading behavior of the distribution of wavefunction intensities P L|φ| 2 . Since, the offdiagonal terms of the effective random matrix Hamiltonian (11) typically decay with system size, they can be treated perturbatively as the thermodynamic limit is approached, in an analysis similar to Ref. [9] The localized eigenstates of the unperturbed effective Hamiltonian are approximated to be δ-function localized in space. As mentioned before, the localized state denoted by |l i is assumed to be localized at x i . To leading order in M , the wavefunction intensity at any site can then be written as Since we are interested in the probability distribution, P L|φ| 2 , we consider its generating function where · denotes the average over sites and disorder realizations. The regular part of the generating function G(s) coming from the perturbative couplings is given by For simplicity, we consider the energy differences (ε l −ε l ) 2 = ∆ 2 l,l belonging to a Gaussian distribution P ∆ = e −∆ 2 /2σ 2 ∆ / 2πσ 2 ∆ , the width of which is assumed to have no scaling with L. With these assumptions, the regular part of G(s) can be expressed as Our quantity of interest, P L|φ| 2 , is the inverse Fourier transform of G reg (s), In the thermodynamic limit, in the integral over s in Eq. (16), the range of s which contributes is such that sL|ψ| 2 ∼ 1. This implies that, if we are interested in the probability of the wavefunction intensity scaling as L −α , we must consider s ∼ L α−1 . Assuming a singleparameter scaling of the distribution P M = P (m = M/M 0 )×M 0 with the scaling M 0 ∼ L −γ/2 , and the convergence of the integral in Eq. (15) as one can expand the latter in a series and truncate at the first order for α < γ with an L-independent constant c. We will discuss the question of the convergence of the above mentioned series and scaling of the distribution function P M for the next section. The leading behavior in P L|φ| 2 is which in terms of α can be expressed as where C 1 depends at most logarithmically with L. The normalization of the probability distribution and the wavefunction intensities put further bounds on α.

Scaling of distribution of off-diagonal elements of random matrix Hamiltonian
In order to justify the assumptions in Sec. 4.2, we analyse the distribution P M in detail. We focus on the distribution of the absolute value |M | as the probability distribution of M l i l j is symmetric with respect to the sign of the matrix element. Constructing the probability distribution P (M ) numerically, Fig. 8(a), shows a strongly non-Gaussian distribution with polynomial tails consistent with the Levy random matrices [37,43].
Rescaling the distribution as with M 0 = L −ν gives a reasonable collapse for ν = 0.9 ± 0.2 in the considered range of system sizes L = 2 9 − 2 13 , (see Fig. 8 gives b 0.5, m 0 4, c = a − b/2, and a 1.1 for ν = 0.9 (see Fig. 8(b)). This confirms the first assumption of Sec. 4.2. However, the accuracy of the extracted parameter ν which assumed to be equal to γ/2 is of order of 20 %.
The integral (15) of the form with Re[A] > 0 converges both at m = 0 (as b < 1) and at m → ∞ (due to exponential decay of the integrand). The first-order expansion (18) in AM 0 is valid for the parameters a > 1 corresponding to the converging first moment |M | in the limit AM 0 → 0. As the best fit gives a 1.1 this justifies the second and the final assumption of Sec. 4.2.
For a more accurate estimate of ν, we calculate the mean |M | and the typical |M | typ ≡ exp ln |M | values of the distribution as both of them should be governed by M 0 . The numerical calculations show the algebraic decay of both mean and typical with the exponents ν mean = 0.78 and ν typ = 0.69 rather close to the expected value ν = γ/2 = 1 − D/2 0.725 (see Fig. 8(c)). The difference between ν mean and ν typ should be considered as an error bar estimate as the points for different L scattered within this interval. As mentioned previously, since τ q and f (α) are related via Legendre transformation one obtains from this analysis and Eq. (22), τ q = (2 − γ)(q − 1) with 2 − γ = 0.53 ± 0.11. This is in rather close agreement with the result numerically obtained in Fig. 2(d) thus validating the random matrix model. Thus, the scaling analysis of the distribution P M confirms both assumptions of Sec. 4.2.
The underlying origin of multifractality is hence the non-trivial mixing of localised states mediated by the delocalised states as an effect of the Floquet drive, thus linking our fully short-range model to an effective long-ranged random matrix ensemble.

Wavepacket dynamics
How can this physics be probed experimentally? An auspicious setting is provided in cold atom experiments, where incommensurate potentials in one-dimension have already been realised, for instance by superposing optical lattices of wavelengths 532nm and 738nm [44] and periodic drives have recently been prominently investigated by sinusoidally modulating laser intensities [19]. As this naturally permits dynamical measurements, we address a conceptually simple process, the spreading of an initially localised wavepacket, |ψ 0 , focusing on signatures of multifractality. Spreading is conveniently quantified by σ(n) via averaged over disorder, where |ψ n = U n |ψ 0 is the wavepacket after n driving periods. The presence of an extensive number of delocalised states in the eigenbasis of U leads to a ballistic leading behaviour σ 2 (n) ∼ n 2 .
In the presence of multifractal states, subleading behaviour emerges as σ 2 (n) ∼ λ 1 n 2 + λ 2 n β . This becomes increasingly visible with a growing ∆V which, as we observe from our numerics, increases the fraction of multifractal states. This is visible in the plot of σ 2 (n)/n 2 , Fig. 9(a), where the amplitude of the ballistic growth is continuously suppressed with increasing ∆V . As a matter of principle, to accurately capture β, it is desirable to remove the dominant contribution of the delocalised states, which can be achieved in theory by removing the projection of the initial wavefunction onto them.
The (normalised) projected initial state |ψ 0 , now has the leading contribution to the spreading from the multifractal states. It shows a much slower growth of σ 2 as shown in Fig. 9(b). The dynamics is in fact subdiffusive with β ≈ 0.72, Fig. 9(c). A collapse of the data suggests a scaling form σ 2 (n) ∼ L 2 F(n/L 2/β ) where F(x) ∼ x β in the scaling regime and F(x) ∼ 1 as x → ∞. This is not unlike the results obtained on hierarchical lattices [45].
The fact that the subleading behavior in the wavepacket spreading due to the multifractal states becomes stronger with increasing ∆V is a consequence of the fact that the fraction of multifractal states in the spectrum of U increases with increasing ∆V . We confirm this by providing the scaling of the Floquet eigenstates for different ∆V , see Fig. 10.

Unequal time density correlators and wavefunction moments
Alternatively, we note that the time-averaged unequal time density correlators can reproduce all moments of the eigenstate wavefunctions and hence the full multifractal spectrum as we  show in this section. Since, the initial state is taken to be δ-function localized at x 0 , |ψ 0 = |x 0 , the quantity of interest is the n-time density correlator measured at the site x 0 , where we use φ j |n x 0 |φ i = φ * j (x 0 )φ i (x 0 ). The infinite time average of R is related to the moments of the eigenstates averaged over the spectrum as follows, In the next step, we take an average over the initial conditions which essentially gives So Eq. (29) implies that R (n) ∞ is the 2(n + 1) th moment of eigenstates averaged over all the eigenstates. For instance n = 0 gives just the normalization, n = 1 gives the IPR, and so on. ∞ obtained by subtracting out the leading behavior by extrapolating the data to L → ∞.
The moments of the eigenstates, averaged over the eigenstates, can thus be useful because they can carry non-trivial information about the presence of multifractal states in the spectrum. For example, consider that the Floquet spectrum has f 1 L localized states, f 2 L delocalized states, and f 3 L multifractal states where f i s denote the respective fractions and f 1 +f 2 +f 3 =1. In such a scenario, let us consider R (1) ∞ which consists of the information of the IPRs of the Floquet eigenstates. The numerical results presented in Fig. 2 suggest that the IPRs of the multifractal states approximately scale as L −τ 2 , where τ 2 = 0.55 ± 0.04. Hence, R ∞ is expected to have an approximate form So by extrapolating the data as function of 1/L to zero, one can obtain the L → ∞ value which can be subtracted, and the leading behavior with L can be obtained. As shown in Fig. 11 the procedure yields a slope of −0.54, which is indeed within the error bars of τ 2 . A similar analysis can be done for the higher moments.

Conclusion
In conclusion, periodically driving a system with a single particle mobility edge can yield robust multifractal states. Note that a single particle mobility edge exists rather generically for two mutually incommensurate potentials in the continuum, with the Aubry-André lattice only a limiting case [20,21]. Thus, with both incommensurate potentials and periodic drives available in state of the art cold atom experiments, our work opens up a new avenue towards realisation of wavefunction multifractality. No fine-tuning is required, with multifractality exhibited by a finite fraction of the Floquet eigenstates for a range of driving strengths. Avenues for future research are evident. For instance, what are the precise properties of the sub-diffusion exhibited by the Floquet multifractal states, and are the observed nontrivial exponents universal? Quite broadly, a key question for Floquet multifractality is the role of dimensionality, known to be a central ingredient for the physics of Anderson localisation. More narrowly, localisation in 'infinite-dimensional' hierarchical systems, where multifractality is ubiquitous, is often considered as toy model for many-body localisation [46] owing to the hierarchical nature of the Fock space. Hence, at a conceptual level, one may ask to what extent a Floquet system hosting multifractal states can be likened to models of many-body localisation, where interestingly the quasiperiodic nature of disorder can have a nontrivial influence on the localisation transition [47].
From a practical perspective, perhaps the most tantalising prospect is to ask how one can use robust multifractal states as basis for the realisation of other interesting phenomena, the possibilities of which are already hinted at by their potential role in enhancing the superconducting transition temperature in a quasi-1D superconductor via algebraic spatial correlations of multifractal states [48].
Note: During the consideration of the manuscript we have become aware of the other work [43] where the similar effective random matrix model (called preferred basis Levy matrix ensemble) was obtained for a completely different system. Our results are consistent with the ones presented in [43] at γ = 1.45 ± 0.04.  Figure 12: IPRs as a function of quasienergies. The top row panels show the IPRs unscaled, scaled with √ L and L, respectively as a function of quasienergy ω. The gray shaded windows denotes the multifractal states, IPRs of which show reasonable collapse for various L when scaled with √ L. To clarify that there are a macroscopic number of multifractal states in the narrow quasienergy band, in the bottom panel we plot the same data but as function of quasienergy index. Figure 13: IPRs at lower driving frequency. The (scaled) IPRs similar to Fig. 2 for a much lower value of driving frequency Ω = 2.5πJ. The persistence of the multifractal states shows that no fine-tuning is needed in the driving parameters.

B Robustness of Floquet multifractality to driving frequency
In order to show that the Floquet multifractality is robust to the driving frequency (Ω) and Ω does not need to be fine tuned to close to but less than the bandwidth of the undriven Hamiltonian's spectrum, in this section, we show numerical evidence for the persistence of multifractality at lower frequencies as well. Recalling that the bandwidth of the undriven spectrum was approximately 2.76πJ, here we choose Ω = 2.5πJ and in fact see an enhancement in the fraction of multifractal states in the Floquet spectrum. The results are shown in Fig. 13 in a similar fashion as Fig. 2. Note that there are almost no states whose IPRs scale as 1/L. This is due to the fact that the lower frequency of the driving forces all the delocalised states in the narrow band at the top of the undriven spectrum ( Fig. 1) to participate in hybridisations. On the other hand, the collapse of the IPRs when scale with √ L is worse. This might be attributed to higher order resonances for which the perturbative treatment based on the two-leg Shirley ladder is insufficient and its detailed analysis constitutes the topic for a future work.