Floquet prethermalization and Rabi oscillations in optically excited Hubbard clusters

We study the properties of Floquet prethermal states in two-dimensional Mott-insulating Hubbard clusters under continuous optical excitation. With exact-diagonalization simulations, we show that Floquet prethermal states emerge not only off resonance, but also for resonant excitation, provided a small field amplitude. In the resonant case, the long-lived quasi-stationary Floquet states are characterized by Rabi oscillations of observables such as double occupation and kinetic energy. At stronger fields, thermalization to infinite temperature is observed. We provide explanations to these results by means of time-dependent perturbation theory. The main findings are substantiated by a finite-size analysis.


Floquet prethermalization and Rabi oscillations in optically excited Hubbard clusters
Junichi Okamoto 1,2* & Francesco Peronaci 3 We study the properties of Floquet prethermal states in two-dimensional Mott-insulating Hubbard clusters under continuous optical excitation. With exact-diagonalization simulations, we show that Floquet prethermal states emerge not only off resonance, but also for resonant excitation, provided a small field amplitude. In the resonant case, the long-lived quasi-stationary Floquet states are characterized by Rabi oscillations of observables such as double occupation and kinetic energy. At stronger fields, thermalization to infinite temperature is observed. We provide explanations to these results by means of time-dependent perturbation theory. The main findings are substantiated by a finite-size analysis.
Coherent optical manipulation of matter is a growing field of study due to the development of intense laser sources. Photoinduced phase transitions in many-body systems are now possible with non-thermal processes without loss of quantum nature [1][2][3][4][5][6] . For example, recent mid-infrared pump-probe experiments revealed exciting phenomena such as light-induced superconductivity [7][8][9] , ultrafast structural transitions 10,11 , and metastable charge ordering 12 , which are driven by short optical excitation of phonons. More direct manipulation of electronic states can be realized through the so-called Floquet engineering by continuous optical excitation [13][14][15] . Notable achievements in this direction are the prediction and realization of dynamical localization [16][17][18] and topological band structures [19][20][21][22] .
In the high-frequency limit, Floquet theory provides a good description of low-energy phenomena in terms of effective static Hamiltonians. In this limit, heating is a rather minor effect even in interacting systems, since the drive frequency is away from any characteristic absorption energy of the system [23][24][25] . Thus, long-lived quasistationary states, which appear to be "thermalized", are realized until the much later heating time scales 26,27 ; such states are termed as Floquet prethermal states (FPSs). On the other hand, when the drive frequency is near an energy scale of a generic interacting system, or its submultiples, heating is expected due to possible photon absorption processes. In isolated systems, this leads to thermalization to infinite temperature 28 . Thus, resonant excitation makes the analysis by the Floquet picture more complicated.
As an alternative to the high-frequency limit, Floquet prethermalization is also observed in systems close to integrability, in which quasi-integrals of motion constrain the dynamics for finite but long times [29][30][31][32][33][34][35][36][37][38][39][40][41][42] . In particular, systems in the Mott-insulating phase of the infinite-dimensional Hubbard models were shown to display extremely long-lived prethermal states, even for frequencies close to resonance 36,41 . Further investigation of these long-lived quasi-steady states is of great importance to advance the Floquet engineering protocols for generic frequencies. To this end, we use a realistic finite lattice geometry and investigate the stability and controllability of possible FPSs.
In this work, we study Floquet prethermalization in two-dimensional (2D) Hubbard clusters under continuous optical excitation. The setup can be realized in solid-state systems such as quantum dot arrays under laser fields [43][44][45][46] or in ultracold atoms in shaking optical lattices 26,27,47 . Starting from the Mott-insulating state, we calculate the driven time evolution by exact diagonalization. We find that Floquet prethermal states emerge even at frequencies resonant with absorption-peak energies, provided a small field amplitude. Remarkably, these long-lived quasi-steady states show Rabi oscillations of observables such as double occupation and kinetic energy. The spectral density shows that the system oscillates between the ground state and the one-photon excited state. For stronger excitation, the system goes to the infinite-temperature state, in general, and the spectral density is spread over many excited states. We elucidate the origin of the Rabi oscillations of the FPSs with the aid of

Formalism
We study a half-filled Hubbard cluster in two dimensions under continuous optical excitation, where c iσ and c † iσ are annihilation and creation operators at site i with spin σ , and l ij = R i − R j the bond vector connecting sites i and j at positions R i and R j . Opposite spin densities n iσ = c † iσ c iσ are subject to a local Hubbard repulsion U. We take = c = e = 1 , and use J 0 = 1 and U = 6 , which ensures that the initial ground state is a Mott insulating state. In the following, we focus on a two-dimensional cluster with L = 10 sites with a periodic boundary condition ( Fig. 1) 48 . Results hold qualitatively unchanged for other geometries and sizes as shown below. Optical excitation is induced by electric fields E(t) = −∂A(t)/∂t along the x-axis, where the vector potential A(t) is switched on as The form of the excitation is motivated by experiments with ultracold atoms, where driving along one axis is commonly used 26,27 . Our envelop function f(t) is chosen for computational convenience, while a constant drive amplitude after a linear ramp is used in Refs. 26,27 . We have also checked that fields along the diagonal direction (x +ŷ) bring similar results.
The initial ground state is calculated by the conventional Lanczos method [48][49][50][51] , and the subsequent time-evolution is implemented by the Krylov-space method with time step dt = 0.02 and Krylov dimension M = 20 52-55 . For each time step, we use the midpoint Hamiltonian, while a higher order Magnus expansion is also possible 56 .
To characterize the time-evolved states, we calculate double occupation d(t) = i �n i↑ n i↓ � /L and kinetic energy E kin (t) = �H J (t)� /L . For ultracold atoms in an optical lattice, these quantities can be measured by radiofrequency spectroscopy and time-of-flight measurements, respectively 26,27,57 . Since time evolution is unitary, isolated systems do not thermalize as a whole. For a large enough system size, however, local observables do thermalize, as described by the eigenstate-thermalization hypothesis 28 . In driven isolated systems, this leads towards the infinite-temperature state, which is characterized by d = 0.25 and E kin = 0 . We further introduce two quantities to obtain a detailed picture of the time evolution. One is the spectral density of the wavefunction |ψ(t)�, which quantifies how many of the energy eigenstates are occupied. For example, in the ground state, the spectral density is peaked at zero energy and κ = 1 , whereas if the system oscillates between two levels we have κ ≈ 2 and two peaks in the spectral density. The infinite-temperature limit may be defined as the spreading of the wavefunction over all the states of the Hilbert space H and thus κ = dim H . We note that Loschmidt amplitudes can be used to calculate the spectral density as well 59,60 .

Results
Now we turn to the two-dimensional cluster with L = 10 . As a basis to interpret the time-dependent phenomena, let us first look at the linear absorption spectrum 50,61 where I x is the total current operator along the x-axis, namely along the direction of the optical field [see Eq. (7)]. H 0 is the unperturbed zero-field Hamiltonian in Eq. (1). In Fig. 2, we plot α(ω) for U/J 0 = 6 with small broadening factor η = 1/L . There are three absorption peaks around U at ω = 4.52, 5.67 , and 6.6, which represent transitions to the one-photon excited states.
Double occupation. First, we discuss the average double occupation after the optical fields are turned on.
Here the onset of the excitation is set to be τ 0 = 10 with a ramp of width τ = 0.1 ; we have checked that the steady states are not sensitive to the specific value of τ . This also indicates that the phase of the drive is not decisive for the steady states. The time-averaged values are plotted as a function of drive frequency ω d and amplitude A 0 in Fig. 3. As expected, we find strong absorption at frequencies resonant with the three peaks in the optical absorption spectrum ω d ≈ 4.52 , 5.67, and 6.63. We notice, however, that even at these resonance conditions, heating is still moderate for small field amplitudes. We find additional two-photon resonances at approximately half the  . We also note that there is no resonance at frequencies commensurate to U, i.e., ω d = U/n with n ∈ N , in contrast to 36 , which uses an infinite dimensional lattice and a strong field amplitude. The discrepancy is probably due to the different lattice geometries or the effect of finite temperatures. In Fig. 4a we plot, for various field amplitudes, the time evolution of double occupation at the one-photon resonance ω d = 4.52 . Similar results are obtained at the other one-photon resonances. Despite the resonance condition, for weak fields, we find FPSs where the double occupation oscillates around a constant value up to t = 2000 (not shown in the figure). The oscillations can be considered as Rabi oscillations, since their frequencies, R , increase approximately linearly with drive amplitude A 0 [Fig. 4e]. Below, we explain the linear dependence by time-dependent perturbation theory. Such a FPS with persistent oscillations is one of our main findings. In contrast, for strong fields, the oscillations are on top of an increasing value saturating around the infinite-temperature limit. Kinetic energy approaches to zero in this limit (see the supplementary material). We have confirmed that the same qualitative behavior holds up to L = 14 sites. As the system size grows, one of the one-photon excited states, which has excitonic nature due to the spin-polaron formation, becomes prominent [62][63][64] . Thus, we also expect that larger systems show Rabi oscillations by applying a resonant light to this state. Similar oscillations are found at the two-photon resonance ω d = 2.40 for weak excitation [ Fig. 4b].
In contrast to the one-photon resonance, the Rabi frequency depends nonlinearly on the field amplitude [ Fig. 4f], because of the two-photon absorption process. We find weak dependence of R on A 0 at small field amplitudes, which is possibly due to the slight deviation of ω d from the resonance. In addition to the slow Rabi oscillations, there exist fast oscillations with frequency ω ≈ 4.5 , which presumably comes from the coherent superposition of the ground state and the excited states with absorption-peak energies. At large field amplitudes, heating and thermalization to infinite temperature occur.
Away from the resonances, Rabi oscillations disappear and two distinct behaviors arise instead. At weak fields and low frequencies [ Fig. 4c], the double occupation increases almost linearly in time. We may consider this linear production of double occupation as an analog of the DC response 65,66 . On the other hand, at weak fields and high frequencies [ Fig. 4d], FPSs without Rabi oscillations appear as in 36,41 . Their steady-state values depend on the drive amplitude A 0 .
As shown in the supplemental material, at resonances, we find monotonic growth of the double occupation and the kinetic energy as the drive amplitude increases. Thus, the total energy also increases monotonically until the system reaches the infinite-temperature state. This amplitude dependence contrasts to a driven two-level system, where absorption saturation occurs for intermediate and large drive amplitudes 67 . The difference comes from the lack of spontaneous emission and the available multi-photon excited states in our model. Spectral density. In order to deepen our understanding of the Rabi oscillations and the dependence on the drive frequency and amplitude, we consider the spectral density S(t, ω) of the time-evolved wavefunction, Eq. (3), and the effective dimension κ , Eq. (4). Figure 5 shows the spectral density S(t, ω) at times multiple of the drive period T d = 2π/ω d for drive frequency close to the one-photon resonance ω d = 4.3 . At small amplitude A 0 = 0.03 [ Fig. 5a] the weight oscillates mostly between the ground state and the one-photon excited state. Thus, the effective dimension κ remains below 2. During the oscillations, the system repeatedly absorb a photon and release it via stimulated emission 68 . Since our model treats only the classical light fields, spontaneous emission does not occur. Thus, the spontaneous decay of the population in the excited states is not observed. Increasing to A 0 = 0.06 [ Fig. 5b], weight on the   Fig. 5c], the spectral density spreads over a wide energy range due to multiple-photon absorption and κ > 10 , meaning that the system has reached the infinite-temperature state, i.e., the double occupation and the kinetic energy reaches the infinite-temperature values. We note that due to the limited number of Krylov dimensions (here M = 50 ), the effective dimension calculated by the Lanczos method cannot reproduce the bulk infinite-temperature behavior, κ → dim H. The spectrum density demonstrates the discrete nature of the multi-photon excited states and gradual filling of these high-energy states with increasing drive amplitudes, which resembles the state-filling effect in semiconductor quantum dots. In the latter, the discrete energy levels are formed by strong spatial confinement of electrons and stronger photoexcitation leads to larger populations in the high-energy levels 69,70 .

Time-dependent perturbation theory.
Here we discuss the physical origin of the Rabi oscillations at the resonant excitation and the amplitude dependence using time-dependent perturbation theory. For this purpose, we define the unperturbed zero-field Hamiltonian H 0 and the time-dependent perturbation, where I x and K x are the current and kinetic-energy operators along the x-axis: Assuming weak perturbation |A(t)| ≪ 1 , we can expand V(t) as 50 which represents the paramagnetic and diamagnetic contributions. Applying first-order time-dependent perturbation theory to the ground state |0� , we find the transition amplitude to another eigenstate of the unperturbed Hamiltonian |n� as The first term describes the one-photon absorption process. If the matrix element �n|I x |0� is peaked at one or a few eigenstates with similar energies, the form has the same structure to the two-level system, and gives Rabi oscillations under continuous excitation A(t) = A 0 cos(ω d t) with ω d ≃ ǫ n 68,71 . From the linear absorption spectrum, Fig. 2a, we see that this matrix element has discrete peaks, which then explains the observed Rabi oscillations. Within the rotating wave approximation, the Rabi frequency is given by The second term, in addition to the second-order terms, describes two-photon excitation, www.nature.com/scientificreports/ which gives rise to the nonlinear amplitude dependence ∼ A 2 0 and to the increase of effective dimension by spreading the weight of the wavefunction to more states [see Fig. 5c]. These terms are responsible for the twophoton resonance at ω d ≃ ǫ n /2.
At the lowest order in A(t) (i.e., retaining only the paramagnetic term), the spectral density is found to be Note that the ground state occupation is calculated as | �0|ψ(t)� | 2 = 1 − n>0 | �n|ψ(t)� | 2 to avoid the use of second-order perturbation. In Fig. 5a, the approximate expression, Eq. (10), is compared with the spectral density obtained by the Lanczos method at a weak amplitude, which well reproduces the Rabi oscillations. On the other hand, for resonant excitation, the transition probability becomes too large and makes the perturbative expression invalid. For larger amplitudes, higher orders in perturbation theory are required.
Expressions similar to Eq. (9) can be derived for other types of time-dependent perturbation. An extensively studied example is the driven-interaction protocol 41,58 , For a weak interaction J 0 /U 1 , the matrix element �n|n i↑ n i↓ |0� is not peaked nor the excited states are degenerate, and we do not expect Rabi oscillations. However, in the limit of J 0 /U ≪ 1 , the excited states are nearly degenerate, and the lower and upper Hubbard bands form an effective two-level system leading to Rabi oscillations 58 .

Finite-size analysis
Up to this point, we have elaborated on a two-dimensional cluster with L = 10 sites. The essential question is if the Floquet prethermal states that we find for the L = 10 cluster survive for larger systems. For example, several studies show that the critical drive amplitude to reach the infinite-temperature state vanishes in the thermodynamic limit for different models 38,40 . In order to check the robustness of the results obtained above, here we simulate various sizes of one-dimensional lattices (chains), ladders, and other two-dimensional clusters (see Fig. 1). For each lattice, the drive frequency is taken to be at the lowest one-photon excitation peak in the linear absorption spectrum. In the inset of Fig. 6, we show an example of time evolution of double occupation for the L = 14 ladder. We confirm that the Rabi oscillations appear for weak drive amplitudes, while the system approaches to the infinite-temperature state for larger amplitudes.
We estimate the critical field amplitudes that separate the prethermal regime d(t) < 0.25 and the infinitetemperature regime d(t) ≈ 0.25 , based on the time-averaged double occupation near the end of the simulations t ∼ 1000 . For comparison of different lattices, we convert the drive amplitude A 0 to the field amplitude E 0 = ω d A 0 . The obtained critical field amplitudes for various lattices are plotted in Fig. 6. For L < 10 , due to the limited available absorption processes, the critical field amplitudes are quite large compared to J 0 . This is consistent with the two-site model (see the supplementary material). For larger systems L ≥ 10 , the critical field amplitudes are of the order of J 0 . However, there is no systematic decrease as the system size grows as Refs. 38,40 , and thus we expect that the results are valid for larger systems.

Conclusions
In this work, we have elaborated on the properties of prethermal Floquet states in optically excited Hubbard clusters.We have demonstrated that the prethermal states exist even at resonance, as far as the drive strength is weak. In particular, the Floquet prethermal states at resonance involve Rabi oscillations, whose frequency scales linearly with the field amplitude at one-photon resonances, and with the square amplitude at two-photon resonances. Experimentally, these results can be tested by using ultracold atoms or quantum dot arrays. We have elucidated the origin of the Rabi oscillations by time-dependent perturbation theory. A finite-size analysis has confirmed the robustness of the main results. The observation of Rabi oscillations in this model suggests possible coherent manipulation of quantum manybody states. Rabi-like oscillations are expected to be a general phenomenon for resonantly driven Hubbard-like clusters; for example, similar oscillations are found after short pulses 72,73 . Considering the coherent nature of Rabi oscillations, the mechanism can be combined with different shapes of drive pulses and used for optical control of correlated quantum states. We expect that introducing dissipation or spontaneous emission 68,74 to our model will weaken the Rabi oscillations while stabilizing the prethermal Floquet states, and it is an interesting open question to investigate their effects on the final steady states.