Following Floquet states in high-dimensional Hilbert spaces

An iterative algorithm is established which enables one to compute individual Floquet states even for many-body systems with high-dimensional Hilbert spaces that are not accessible to commonly employed conventional methods. A strategy is proposed for following a Floquet state in response to small changes of a given system's Hamiltonian. The scheme is applied to a periodically driven Bose-Hubbard chain, verifying the possibility of pseudoadiabatic Floquet state following. In particular, it is demonstrated that a driving-induced Mott insulatorlike target Floquet state can be populated with high efficiency if the driving amplitude is turned on smoothly but not too slowly. We conclude that the algorithm constitutes a powerful tool for the future investigation of many-body Floquet systems.


I. INTRODUCTION
The experimental and theoretical investigation of periodically time-dependent many-body quantum systems, nowadays often dubbed Floquet systems, has turned into a remarkably fruitful area of physics in recent years, comprising, among many others, the dynamics of cold atomic quantum gases in periodically driven optical lattices [1], the principles underlying Floquet time crystals [2,3], and fundamental aspects of nonequilibrium statistical physics [4,5].
Advances on the theoretical side of this field appear to be impeded by the fact that, while it is still feasible to solve the time-dependent Schrödinger equation numerically for many such systems of interest, there is a notable lack of powerful methods for computing the systems' Floquet states. Since these states constitute a natural basis which fully incorporates the periodic time dependence, thus providing a key for understanding both the shorttime and the long-time behavior of periodically driven quantum systems, lack of knowledge of these states tends to limit one to a mere description of the results of numerical calculations, precluding deeper conceptual insight.
So far, numerical methods for treating large Floquet systems either make heavy use of specific properties of the respective system [6] or require supercomputing facilities [7], typically enabling one to perform exact calculations for systems with a Hilbert space with a dimension on the order of 10 4 . In the present paper we establish a general computational strategy for calculating the Floquet states of even higher-dimensional systems with only modest numerical effort. This progress is achievable if one does not require all Floquet states and their quasienergies, but wishes to obtain information on specific individual states. In contrast to a variational principle which has been suggested recently for the same purpose [8], here we propose an iterative algorithm that does not require a variational ansatz and can be executed whenever a sufficiently efficient scheme for propagating states in time is available. As an experimentrelated application, we employ the exact Floquet states obtained in this manner for a periodically driven finite Bose-Hubbard chain in order to investigate the peculiarities of pseudoadiabatic Floquet-state preparation in the system's high-frequency regime [9]. This is a somewhat subtle topic, since a proper adiabatic limit, referring to a turn-on of the driving amplitude that proceeds "infinitely slowly", cannot be expected to exist when the Bose-Hubbard system becomes large, owing to the fact that the gap condition required by the standard adiabatic theorem [10] cannot be satisfied; moreover, for an infinite system the quasienergy eigenvalues probably are nowhere differentiable with respect to the adiabatically changing parameter [11,12]. We argue that, nonetheless, there may be a window of opportunity involving driving amplitudes which vary so fast that these pathologies remain almost unresolved, but still sufficiently slow to enable effectively adiabatic following at least to a high degree; the detailed experimental verification of this scenario with cold atoms in shaken optical lattices might constitute a rewarding challenge in the near future.
We proceed as follows: In Sec. II we review some basic elements of the Floquet approach to periodically timedependent quantum systems, thereby establishing our notation in the form it will be required later. We then introduce our iterative algorithm for computing Floquet states of high-dimensional systems in Sec. III, and discuss a method for "following" Floquet states in response to small changes of the system's Hamiltonian in Sec. IV. This type of computational following is still not directly related to adiabatic following but based on the likeness of Floquet states. In Sec. V we apply these concepts to a periodically driven one-dimensional Bose-Hubbard chain, documenting both the viability of the iterative algorithm and the usefulness of the Floquet state-following scheme. In Sec. VI we then tie open ends together and scrutinize what we call pseudoadiabatic following, demonstrating that a periodically driven many-body system may respond to a slowly changing driving amplitude by actually taking the Floquet path sorted out by the computational following procedure. Some conclusions regarding possible experimental implications are drawn in the final Sec. VII.

II. THE FLOQUET PICTURE
Consider a quantum system with a periodically timedependent HamiltonianĤ(t) =Ĥ(t + T ) acting on a Hilbert space H, giving rise to the time-dependent Schrödinger equation Then the unitary time-evolution operator of that system possesses the Floquet product form [13][14][15][16] where the unitary operatorP (t) =P (t + T ) inherits the periodic time dependence of the Hamiltonian, witĥ P (0) =1, and the operatorĜ is self-adjoint. Thus, the transformed states obey the Schrödinger equation in whichĜ appears as an effective time-independent Hamiltonian. The quasienergies ε n constitute the spectrum ofĜ [15]. Requiring H to be of finite dimension, the "stroboscopic" approach to periodically time-dependent quantum systems rests on the eigenvalue problem U (T, 0)|n = exp(−iε n T / )|n (5) which is posed by the unitary one-cycle evolution opera-torÛ (T, 0) = exp(−iĜT / ) on H. Since all eigenvalues exp(−iε n T / ) fall on the unit circle of the complex plane, the quasienergies are thus determined up to an integer multiple of 2π/T ≡ ω.
If H is not of finite dimension, the eigenvalue problem may become quite difficult. For example, the quasienergy spectrum of a linearly driven harmonic oscillator [17] is pure point if the driving is nonresonant but absolutely continuous if the driving frequency matches the oscillation frequency of the undriven oscillator. More generally, the question under which conditions the quasienergy spectrum of a periodically driven quantum system may become continuous has been termed the "quantum stability problem" in the literature [18]; some sophisticated theorems have been developed which allow one to exclude the presence of a continuous spectrum in particular cases [19][20][21]. In order to avoid such complications we require H to be of large, but finite dimension from here on.
The "extended" viewpoint emerges when introducing the Floquet functions |u n (t) =P (t)|n (6) which are T -periodic by construction, |u n (t) = |u n (t + T ) .
By virtue of the representation (2), the Floquet states |ψ n (t) = |u n (t) exp(−iε n t/ ) (8) then constitute a set of particular solutions to the Schrödinger equation (1) which is orthogonal and complete in H at each instant t. Every solution to Eq. (1) can be expanded with respect to these states with timeindependent coefficients a n , |ψ(t) = n a n |u n (t) exp(−iε n t/ ) , showing that the Floquet functions adopt a role which is conceptually similar to that of the energy eigenfunctions of time-independent systems. Moreover, inserting a Floquet state (8) into the Schrödinger equation, one readily derives This is an eigenvalue equation for the quasienergies akin to the stationary Schrödinger equation, posed by the quasienergy operator on an extended Hilbert space L 2 [0, T ] ⊗ H of timeperiodic functions [22]. Denoting the scalar product on H by · | · , the scalar product on L 2 [0, T ] ⊗ H is naturally given by since the time t is an additional coordinate in this extended space. When regarding a Floquet function |u n (t) no longer as a time-dependent function on H but rather as an element of L 2 [0, T ] ⊗ H it is written as |u n , so that the quasienergy eigenvalue equation (10) takes the formK |u n = ε n |u n .
There is a seemingly simple but important implication that marks a crucial difference between this quasienergy eigenvalue problem and the more familiar energy eigenvalue problems encountered with timeindependent Hamiltonian operators: If |u n (t) is a Floquet function with quasienergy ε n , and if m is any positive or negative integer, then |u n (t)e imωt is a further T -periodic Floquet function with quasienergy ε n + m ω, where, again, ω = 2π/T ; all these different solutions are required for the completeness relation pertaining to the eigenfunctions ofK in L 2 [0, T ] ⊗ H. On the other hand, solutions differing by m only give rise to the same Floquet state in H, since Expressed pictorially, the spectrum of the quasienergy operator (11) is obtained by "unrolling" the eigenvalues exp(−iε n T / ) of the one-cycle evolution operator from the unit circle to the infinite real energy axis. More technically, a "quasienergy" should not be regarded as a single eigenvalue, but rather as an infinite set of equivalent representatives spaced by ω, implying that quasienergies cannot be ordered with respect to their magnitude without additional conventions. The spectrum ofK thus consists of infinitely many identical Brillouin zones of width ω, with each zone containing one quasienergy representative of each Floquet state. When the dimension of H becomes high, the distance between the quasienergy representatives within one Brillouin zone necessarily becomes small. On the other hand, being the eigenvalues of a Hermitian operator, the quasienergies are subject to the von Neumann-Wigner noncrossing rule and hence generally do not cross when only one parameter of the quasienergy operator is varied, unless such a crossing is permitted by a symmetry [23]. Hence, excepting integrable systems, which actually do possess smooth quasienergies such as the nonresonantly driven harmonic oscillator [17], the quasienergies of a high-dimensional periodically time-dependent quantum system form an intricate, dense net of avoided crossings when viewed as functions of one of its parameters, such as the amplitude of a driving force; the full complexity of this net may remain unresolvable when limited to the accuracy attainable by numerical computations [11,12]. The question how to "follow" an individual Floquet state on parameter changes in such a net is, therefore, far from trivial.

III. ITERATIVE COMPUTATION OF FLOQUET STATES
Our present approach to computing individual Floquet states of periodically time-dependent quantum systems is based on the following theorem, resorting to the stroboscopic viewpoint: LetÛ 1 =Û (T, 0) be the one-cycle evolution operator of a quantum system defined on a finite-dimensional Hilbert space H, which possesses a time-periodic Hamiltonian H(t) =Ĥ(t + T ), and consider the functional Γ γ on H which is given by Then one has, for any γ ∈ R and any normalized |ψ ∈ H, This is easily shown: After expanding |ψ with respect to the eigenvectors |n = |u n (0) = |u n (T ) ofÛ 1 , obtaining |ψ = n a n |n , one finds = n a n 2 + 2 cos(ε n T / − γ) |n (18) and, hence, From this representation the theorem follows immediately, since 1 ≥ cos(ε n T / − γ) ≥ −1, and |ψ is assumed to be normalized, n |a n | 2 = 1.
In particular, the maximum Γ γ [ |ψ ] = 4 of the functional Γ γ is attained if |ψ = |n equals one of the eigenvectors ofÛ 1 , and γ = ε n T / equals the corresponding phase. This observation enables one to invoke iterative methods for computing eigenvectors possessing the largest eigenvalue, updating the phase γ at each step. Here we employ a power method based on the scheme combined with subsequent normalization of |ψ new . The real parameter α can be adjusted empirically in order to speed up the convergence; evidently, one requires α > −2 for filtering out the desired eigenvector of the kernel of the functional Γ γ which belongs to the largest eigenvalue. To be precise, we propose the following three-step algorithm: Step 1: Choose a convenient initial state |ψ 0 and propagate this state in time over one period T , obtaininĝ U 1 |ψ 0 . If then where δ > 0 is a predefined small error tolerance, |ψ 0 already is an eigenvector |n ofÛ 1 , to the accuracy thus specified. Hence, the algorithm terminates, and the corresponding quasienergy is obtained from the relation Otherwise, that is, if the condition (21) is not satisfied, compute yielding the vector Step 2: Now perform a propagation backward in time over one period T to obtainÛ † 1 |ψ 1 and compute Step 3: Compute and normalize, obtaining With this new |ψ 0 go back to Step 1 and repeat until the algorithm terminates.
Observe that this scheme for computing |n = |u n (0) and the quasienergy eigenvalue ε n can be executed already if a sufficiently efficient routine for propagating states in time is available. Thus, it is well applicable even if the dimension of H is so large that the computation and subsequent diagonalization of the one-cycle evolution operator, or alternatively the diagonalization ofK in the extended Hilbert space, are rendered impracticable.

IV. COMPUTATIONAL FOLLOWING OF FLOQUET STATES
Now suppose that a Floquet function |u 1 of some quasienergy operatorK 1 with eigenvalue ε 1 is already known,K Next, suppose that the quasienergy operator is modified by adding a piece λV , where λ is a dimensionless parameter, givingK 2 =K 1 + λV . It then appears natural to seed the algorithm devised in Sec. III, searching for a Floquet function of the new operatorK 2 , with the old Floquet function |u 1 pertaining toK 1 . This is based on the continuity assumption that for sufficiently small λ there should be an eigenfunction ofK 2 which closely resembles |u 1 . Would it then be possible to make a useful a priori statement concerning the performance of the algorithm? To this end, consider the expression whereŶ is an operator acting on the extended Hilbert space, |u is an element of that space, z is a scalar, and double angular brackets indicate the scalar product (12) on L 2 [0, T ] ⊗ H. In view of Eq. (28), one evidently has this identity implies a variational principle for Floquet states [8]. Inserting the old, known Floquet function |u 1 and its quasienergy ε 1 , but the new quasienergy operator K 2 into this expression (29), it jumps to the nonzero value More elaborately, computing and inserting this ε instead of ε 1 , one derives We emphasize that this strategy for following an individual Floquet state in parameter space does not necessarily result in following by continuity with respect to λ. If, for instance, the relevant quasienergy functions ε n (λ) are broken by a large number of tiny avoided crossings, indicating weak multiphoton-like resonances, and if the increment ∆λ is larger than the typical width of these avoided crossings, the initial state will be followed diabatically, that is, ignoring the avoided crossings as if they did not exist. If, on the other hand, ∆λ is comparable to the size of the avoided crossings, they will be resolved, and the numerically computed Floquet state will follow its quasienergy continuously, that is, adiabatically. It remains to be explored whether this potential sensitivity of the computational following procedure to the stepsize ∆λ can be exploited for extracting useful information about the physics of the respective system under investigation.

V. APPLICATION: THE PERIODICALLY DRIVEN BOSE-HUBBARD MODEL
The Bose-Hubbard model constitutes an idealized lattice system of theoretical many-body physics, embodying the competition between delocalization due to kinetic energy and localization due to repulsive potential energy, thus giving rise to a superfluid-Mott insulator quantum phase transition [24]. Here we consider a one-dimensional Bose-Hubbard chain, as specified by the Hamiltonian with nearest-neighbor tunnelinĝ and on-site interaction where j labels the chain's sites in consecutive order,â j (â † j ) annihilates (creates) a Bose particle at site j, implying the canonical commutator [â j ,â † k ] = δ jk , whilê n j =â † jâj denotes the number operator at that site. Moreover, J is the nearest-neighbor hopping matrix element, and U is the repulsion energy contributed by each pair of particles occupying a common site. In the limit of infinite chain length, the phase transition occurs at (J/U ) c = 0.297 ± 0.01 for unit filling, that is, when the chain is occupied by one particle per site [25,26].
In order to equip this model with a periodic time dependence, thus admitting an additional wealth of Floquet physics, we introduce a monochromatic homogeneous driving force described bŷ and investigate the periodically driven Bose-Hubbard this driving scheme can be implemented experimentally with ultracold atoms in periodically shaken optical lattices [27][28][29][30]. Some salient features of this system (38) can already be deduced by calculating the matrix elements of its quasienergy operatorK for the Floquet-Fock functions [27] |{n j }, m spanning its extended Hilbert space, where |{n j } denotes a Fock state with n j particles occupying the site labeled by j, and m is the integer having shown up before in Eq. (14): While the on-site contributions remain diagonal, the matrix elements of the tunneling operator are found to read sinceĤ tun transfers one particle by one site along the chain, with the plus (minus) sign referring to a tunneling process to the left (right). Invoking the Jacobi-Anger identity for the Bessel functions J of integer order , this becomes In the high-frequency regime, in which the width ω of the quasienergy Brillouin zone becomes the dominant energy scale and "small" avoided crossings of quasienergy representatives belonging to different m can be ignored, one may neglect all Bessel function factors J m −m (R/ ω), except for m − m = 0. Thus, one arrives at a time-independent effective Hamiltonian which differs from the original, undriven Bose-Hubbard model (34) only through the replacement of the hopping matrix element J by the "renormalized" hopping strength [27] J eff = J J 0 (R/ ω) ; this effective Hamiltonian is a simple, but often apparently sufficient approximation to the exact operatorĜ introduced on general grounds in Eqs. (2) and (4). Since the iterative scheme established in Sec. III does not require any such approximation, we are now in a position to subject this high-frequency approximation to a "hard" numerical test.
The interest in such a test stems from the following deliberation: Assume that J/U > (J/U ) c , so that the undriven chain (34) is in its superfluid ground state. If then the scaled driving amplitude R/( ω) is increased from zero toward the first zero j 0,1 ≈ 2.405 of the Bessel function J 0 , the effective hopping strength (45) decreases monotonically to zero. This implies that the effective time-independent model enters the Mott insulator regime at a certain "critical" driving strength, at which J eff /U = (J/U ) c ; for even stronger driving, the energy ground state of the effective model is separated from the excited states by a finite gap. Thus, the effective model predicts the emergence of a driving-induced Mott insulator state [27]. The full Floquet system (38), however, does not possess such a gapped ground state. Considering an infinitely long chain for the sake of the argument, the quasienergy of a Floquet state originating on activation of the drive from the undriven chain's Mott insulator ground state would be embedded in the continuum of excited states due to the Brillouin zone structure of the quasienergy spectrum. Therefore, this state should turn into a resonance, that is, into a decaying state characterized by a Lorentzian peak in the spectral density with a certain width determining its life time, akin to the familiar Floquet resonances of atomic states in strong laser fields [31], so that the decay width may be interpreted as an inverse characteristic heating time.
Returning to large but finite driven Bose-Hubbard chains which are more realistic from an experimental point of view, the quasienergy continuum resulting from the excited states gives way to a huge number of discrete quasienergy eigenvalues filling each Brillouin zone. Thus, on variation of the driving amplitude, a Floquet state associated with a Mott insulator state would not be protected by a gap but would have to avoid a plethora of other states instead, to the effect that its quasienergy cannot depend smoothly on the driving amplitude, but is "broken" by countless tiny avoided crossings. Indeed this "roughness" of the quasienergy functions should enable one to estimate the corresponding decay or heating times, by a procedure mimicking the so-called L 2 stabilization method [32,33]. Here we do not attempt to determine these heating times, but focus on the first and foremost question: If one goes beyond the convenient, but essentially uncontrolled high-frequency approximation (45), would it actually be possible to identify drivinginduced Mott insulatorlike Floquet states despite the above caveats?
To this end we study a finite Bose-Hubbard chain (34) possessing M = 11 sites ranging from j min = −5 to j max = 5, and specify unit filling, so that the chain is occupied by N = 11 Bose particles. The dimension d of this system's Hilbert space H then is determined by the binomial coefficient posing a serious challenge to more traditional methods commonly used for investigating Floquet systems. Adopting the above reasoning, we fix the parameter J/U = 1/3 > (J/U ) c , placing the undriven system in the superfluid regime. The energy spectrum of this chain ranges from E min /U = −3.64 to E max /U = 55.25. Moreover, we select the driving frequency ω/U = 14/3, thus aiming for the high-frequency regime in which the approximation (45) should be viable [27]. This gives E min /( ω) = −0.78 and E max /( ω) = 11.84, implying that the chain's energy spectrum covers more than 12 quasienergy Brillouin zones. With this choice of parameters one finds J eff /U = (J/U ) c for R/( ω) ≈ 0.67. Hence the driven system should enter a Mott insulatorlike regime at about this value of the the driving amplitude, and stay therein for stronger driving. The quasienergy operatorK of the driven chain remains invariant under a generalized parity transformation in its extended Hilbert space, corresponding to the spatial reflection j → −j combined with a shift in time by half a period, t → t + T /2. Since the Floquet functions are odd or even under this generalized parity, and eigenvalues belonging to the same parity class are subject to the von Neumann-Wigner noncrossing rule, each of the d quasienergy respresentatives falling into one Brillouin zone inevitably is perforated by a large number of avoided crossings when, for example, the driving amplitude is varied.
As a preliminary application of the strategy put forward in Sec. IV, we try to follow the Floquet state which develops from the superfluid energy ground state of the undriven chain in response to a time-periodic drive (37) with increasing driving amplitude, so that the scaled amplitude R/( ω) now plays the role of the parameter λ. The calculations were performed on a standard laptop computer equipped with a GPU (GeForce GTX 1060 Mobile). Figure 1 depicts the performance of the algorithm when scanning the interval 0 ≤ R/( ω) ≤ 7 with stepsize ∆R/( ω) = 7/120 and tolerance δ = 10 −4 ; here, the parameter α employed in Eq. (26) has been set to α = −1.9. These results are fairly encouraging: Indeed, the procedure always converges to a Floquet state after less than 30 iterations. Most significantly, the expected entrance into a Mott insulatorlike regime is reflected by the number of iterations; in the initial superfluidlike regime, the state is increasingly harder to follow when R/( ω) is enhanced. This behavior changes after entering the potential Mott regime; for R/( ω) > 2, very few iterations are required for settling down to the new Floquet state. Although there can be no sharp transition here, related features are observed when monitoring the functional Var 1 (Ĥ drive (t)/R), as defined by Eq. (33). In particular, the low values of this functional for large driving amplitudes appear to match the incompressibility which is characteristic of a genuine Mott insulator [24], since the system's compressibility is determined by the variances of the particle-number operators showing up in H drive (t).
In Fig. 2 we compare the quasienergies of the exact, numerically computed Floquet states obtained by this following procedure to the corresponding ground state energy of the effective Hamiltonian with the appropriately renormalized hopping strength (45). The agreement is quite impressive except for scaled driving amplitudes R/( ω) close to the zeros j 0,1 ≈ 2.405 and j 0,2 ≈ 5.520 of the Bessel function J 0 : At these zeros, the ground-state energy of the effective time-independent model vanishes, while the exact quasienergies rise to slightly higher values. Still, it needs to be kept in mind that the exact Floquet state considered here is coupled to a background of many others by large numbers of avoided crossings; yet, these avoided crossings are too narrow to be individually resolved on the scale imposed in Fig. 2. In contrast, the energy of the effective model's ground state varies in a perfectly smooth manner with the driving amplitude. Under the present conditions the Floquet state following procedure thus provides a coarse-grained image of the exact eigenvalue, selecting the followers diabatically, that is, according to maximum likeness of the states. The fact that no "roughness" of the quasienergy can be detected in Fig. 2 also implies that the corresponding Floquet resonance is fairly long-lived.
On the level of the states, the fair agreement between the effective and the full, periodically time-dependent system in the high-frequency regime is emphasized fur- ther by the following comparison. Figure 3 shows contour plots of the momentum distribution functions of the effective model, where k is a wave number and d indicates the lattice constant; |ϕ 0 is the energy ground state of the effective model at the respective value of R/( ω). Such momentum distributions are accessible to measurement with the help of time-of-flight interference experiments [34][35][36]. One observes marked peaks of this distribution function at k = 0 mod 2π/d in the superfluid regime at low driving amplitudes, since the system tends to condense into the single-particle state at the bottom of the energy band of the noninteracting system with U = 0. For scaled driving amplitudes R/( ω) in the vicinity of the zeros j 0,1 and j 0,2 of J 0 this distribution function is almost flat, as it should, since the effective hopping matrix element (45) is small there. Interestingly, between these zeros one finds a weak precursor of a new superfluid state conforming to a condensate in the single-particle state k = π/d mod 2π/d, as corresponding to the upper edge of the band (48); this is naturally explained by the fact that the sign of the effective hopping strength is inverted between the zeros. When the energy eigenstates |ϕ 0 entering the functions (47) are replaced by the exact Floquet functions |u(0) computed by following the ground state of the undriven chain with increasing driving amplitude, providing the exact distribution functions C exact (k), one obtains a plot that is almost indistinguishable from Fig. 3. There-  (47) predicted by the effective model and that obtained for the exact Floquet states computed by the following scheme, averaged over k. Note how the scale here differs from that in Fig. 3.
fore, we depict in Fig. 4 the squared difference averaged over k, revealing that the deviation between the two functions figures on the subpercent level. Thus, under the conditions considered in this section, the simple time-independent, effective model performs fairly well when compared with the full, periodically driven Bose-Hubbard chain (38), and the question posed above can be answered in the affirmative: The driving-induced Mott insulatorlike Floquet state predicted by the highfrequency approximation (45) is no artifact.

VI. PSEUDOADIABATIC FOLLOWING OF FLOQUET STATES
Going beyond the computational Floquet statefollowing strategy suggested in Sec. IV and practiced for one particular example in Sec. V, there also is the more experiment-related question of whether the actual state of a system would be able to follow one of its Floquet states in an adiabatic manner when its parameters are changing in real time. Indeed there is an adiabatic principle for Floquet states [37,38], which formally resembles the celebrated adiabatic theorem of quantum mechanics [10]. This principle allows one, for instance, to define generalized π-pulses together with a generalized area theorem, which govern transitions in periodically driven multilevel ladder systems in close analogy to their two-level analogs [39]. However, for many-body Floquet systems with a high-dimensional Hilbert space, as considered here, the inevitable multitude of avoided crossings renders adiabatic following, if understood in the formal mathematical sense, impossible. Yet, a closely related option emerges: If the relevant quasienergies are punctured only by narrow avoided crossings below a certain scale, and if the parameter variation, while sufficiently slow, still proceeds so fast that these narrow avoided crossings are traversed almost diabatically by Landau-Zener-like Floquet-state transitions [40], highly diabatic quantum motion actually appears as effectively adiabatic motion on coarse-grained quasienergy eigenvalue surfaces [16,37]. This is termed pseudoadiabatic following here, enabling one to design the parameter variation such that desired Floquet target states are populated with high probability.
For demonstrating the feasibility of this concept we now perform a series of numerical experiments: Initially, at time t = 0, the state |ψ(t = 0) of the system is given by the energy ground state of the undriven Bose-Hubbard chain (34). Then the driving amplitude is smoothly ramped up according to the protocol reaching its final value R max at time t = zT , that is, after z driving cycles T = 2π/ω. The resulting state |ψ(t = zT ) , as computed by numerical integration of the Schrödinger equation, is then projected on the Floquet function |u(0) calculated for R = R max by means of the procedure outlined in Sec. V, yielding the overlap Similar studies have been reported by Poletti and Kollath in Ref. [9], the difference being that here we target an exact Floquet state instead of the ground state of the effective Hamiltonian. Figure 5 depicts numerical results obtained in this manner for R max /( ω) = 4, again with J/U = 1/3 and ω/U = 14/3: While the absolute value of the overlap (51) is only small if the turn-on of the driving force proceeds within a few cycles, because the initial superfluid state is given insufficient time to adjust to the Mott insulatorlike target state, it becomes appreciably higher than 0.95 when the final amplitude is reached after some 10 cycles. Yet, the inset of Fig. 5 reveals that this trend cannot be extended to much longer turn-on times: For z > 70 the population of the target Floquet state starts to decrease. This is the expected consequence of the many tiny avoided crossings that have remained unresolved in Fig. 2 but nonetheless are effective: Even if each individual avoided crossing is still traversed almost diabatically when z is increased, they conspire by their sheer number to divert an increasing fraction of the evolving state into the anticrossing Floquet states; this fraction cannot reach the target state. This is, in a nutshell, a visualization of the window of opportunity referred to in the Introduction: Pseudoadiabatic following of Floquet states is possible if (i) their quasienergies do not exhibit large avoided crossings, and (ii) the parameter variation proceeds reasonably slowly, but not too slowly.
A variation of this theme is shown in Fig. 6, which depicts the corresponding results for R max /( ω) = 1: Although the final amplitude is four times smaller here than that employed in the preceding Fig. 5, so that the target state resembles the initial state more closely, the decrease of the target-state population now becomes notable for z ≈ 20 already. The reason for this unexpected behavior can be spotted in Fig. 1: With R max /( ω) = 1 the final amplitude falls into the regime where Var 1 (Ĥ drive (t)/R) is large, signaling a relatively strong sensitivity of the Floquet states to the driving amplitude; this is felt by the state if it spends too much time in this regime.
Finally, Fig. 7 shows data obtained for R max /( ω) = 7, beyond the second zero j 0,2 of J 0 . Here the overlap first decreases with increasing length of the turn-on before it increases again. This behavior again may allow one to draw deductions concerning the instantaneous Floquet states involved: Now the initial state already exhibits a certain likeness to the target Floquet state, as is reflected by Fig. 3; this likeness is advantageous in case of an almost sudden turn-on. However, the time-evolving state is forced to undergo a substantial structural change for j 0,1 < R/( ω) < j 0,2 , where the sign of the effective hopping matrix element (45) is inverted, requiring a sufficiently moderate growth rate of the driving amplitude.

VII. CONCLUSIONS
The iterative algorithm proposed in this work for calculating Floquet states of periodically time-dependent quantum systems does not yield the systems' full spectrum, but provides selected individual states and their quasienergy eigenvalues. Nonetheless, its attractive feature stems from the fact that it requires neither diagonalization of the quasienergy operator in the extended Hilbert space [6] nor computation of the one-cycle evolution operator [7]. Thus, it provides a powerful tool for scrutinizing many-body Floquet systems that are not accessible to these older standard techniques due to the high dimensions of their Hilbert spaces. Since a Floquet system does not possess a proper ground state, and since it is neither possible nor even desirable to compute all Floquet states of a truly highdimensional system, one necessarily faces the question of how to select suitable, hopefully typical Floquet states for inspection. Here we have suggested to "follow" a Floquet state in parameter space by slightly modifying the quasienergy operator, seeding the iterative algorithm with the old Floquet state, and letting it relax to the new one; this strategy amounts to following according to maximum likeness of the Floquet states.
Our explorative application of this strategy to a Bose-Hubbard chain subjected to a high-frequency drive was facilitated by the observation that here the quasienergy which emerges from the system's ground-state energy un-dergoes only tiny avoided crossings, unresolved on the scale of Fig. 2 when regarded as function of the driving amplitude. In such favorable cases one encounters pseudoadiabatic Floquet state following in response to a driving amplitude that changes smoothly but still so fast that all these tiny avoided crossings are traversed diabatically. The attempt to reach a formal adiabatic limit, as corresponding to a change of the driving amplitude that proceeds "infinitely slowly", would be fruitless [11] because countless Landau-Zener-like transitions would then divert the evolving state into undesired channels.
It should also be pointed out that the comparatively simple dynamics encountered in the present study within the high-frequency regime will give way to more complicated dynamics at lower frequencies, when resonances make themselves felt [9]. With the tools provided here, which do not rely on any high-frequency approximation, the Floquet analysis also of such resonant dynamics can now be performed routinely.
From an experimental viewpoint, the paradigmatic results depicted in Figs. 5 -7 suggest what may be termed "turn-on spectroscopy" with cold atoms in optical lattices: Prepare a superfluid or a Mott insulator state, then start shaking the lattice according to a deliberately designed turn-on protocol, and perform time-of-flight imaging after the shaking has become strictly periodic in time. Results obtained in this manner will depend on details of the respective turn-on protocol, and model calculations of the type discussed in this paper may be of profound help for relating the different signatures observed for different protocols to the underlying quantum nonequilibrium many-body dynamics.