Non-equilibrium Floquet steady states of time-periodic driven Luttinger liquids

Time-periodic driving facilitates a wealth of novel quantum states and quantum engineering. The interplay of Floquet states and strong interactions is particularly intriguing, which we study using time-periodic fields in a one-dimensional quantum gas, modeled by a Luttinger liquid with periodically changing interactions. By developing a time-periodic operator algebra, we are able to solve and analyze the complete set of non-equilibrium steady states in terms of a Floquet-Bogoliubov ansatz and known analytic functions. Complex valued Floquet eigenenergies occur when multiples of driving frequency approximately match twice the dispersion energy, which correspond to resonant states. In experimental systems of Lieb-Liniger bosons we predict a change from powerlaw correlations to dominant collective density wave excitations at the corresponding wave numbers as the frequency is lowered below a characteristic cut-off.

In this Letter we now provide the many-body eigenstate solution and report resonance phenomena in onedimensional (1D) interacting quantum systems with time-periodically driven parameters. Our analysis applies to time-periodic driving of generic Tomonaga-Luttinger liquids (TLL), which describe a large class of effectively 1D many-body systems [48] and can also be realized using ultra-cold gases [49][50][51][52]. The timeevolution of an initially prepared state in a TLL has been calculated before [53][54][55][56][57][58][59][60], but much less is known about the nature of possible non-equilibrium steady states under periodic driving. It is therefore desirable to obtain the full eigenbasis of the Floquet eigenvalue problem, which gives systematic information about all stable steady states and their corresponding dominant correlations. We now obtain the explicit steady state solutions of the time-dependent Schrödinger equation of a quantum many-body system in terms of a time-periodic operator algebra, which not only allow a full analysis using closed analytic functions, but also show regions of instabilities in frequency and momentum space. We therefore predict large-amplitude density waves at the characteristic wave vectors in trapped ultracold boson systems.
Model. We will develop a Floquet-Bogoliubov ansatz for general driven TLL models. To make concrete predictions for decoupled 1D tubes of interacting bosons in optical lattices [49][50][51][52] we choose the Lieb-Liniger Hamiltonian [61] as a starting point where = 1 and g = 2a0 ma ⊥ (a ⊥ −1.03a0) is the 1D onsite interaction strength, which is tunable via the 3D scattering length a 0 and the perpendicular confinement length a ⊥ [62][63][64][65]. The static system is integrable and correlations are known to be well described by a TLL model in the long-wavelength limit q < q c [48,63,64] are SU(1,1) generators [66][67][68] in terms of bosonic operators b † L/R,q , which create left-and right-moving density waves at wave-vector q [48]. For the Lieb-Lininger model the Luttinger parameter K is known exactly [63][64][65]. It depends only on the ratio mg/n, where n is the 1D particle density. The cutoff q c above which the TLL description fails has also been determined [64]. The scattering amplitudes g 2 and g 4 are rescaled from the traditional "g-ology" scheme [69] and related to v and K by [48]  where v F = πn/m. Values of g 2 = g 4 = (1/K 2 − 1)/2 for the Lieb-Liniger model are shown in Fig. 1.
We now turn to systems with time-periodically changing control parameters a 0 and a ⊥ , which will result in time-periodic couplings g, g 2 , and g 4 , all of which can be determined exactly. Any desired time-periodic couplings can be created by suitable fields given by the inverted relations, including a pure sinusoidal behavior [70] 2g 2 (t) = 2g 4 (t) =ρ + ρ cos ωt (4) with constant parametersρ and ρ. We will later consider more general behavior. For the Lieb-Liniger model it is known that vK = v F and K > 1 [63,64], so that from Eq. (3) −1/2 < g 2 = g 4 < 0 as shown in Fig. 1. Floquet ansatz. We now seek to solve the timedependent Schrödinger equation ı∂ t |Ψ(t) = H q (t)|Ψ(t) for each momentum q separately, which remains a good quantum number and can be omitted in the following. According to Floquet theory [27,46,47,71] there exists a complete set of quantum numbers n for steady state solutions |Ψ n (t) = e −ıǫnt |u n (t) . Here |u n (t) = |u n (t + T ) with T = 2π/ω obey the Floquet equation where ǫ n are the Floquet quasienergies. We now wish to map this problem onto a static eigenvalue problem [72] Floquet theory has been reviewed extensively [27,46,47,71], but the ansatz (6) goes beyond the usual timeevolution approach since it makes the problem static, diagonalizes it in the original Hilbert space, and determines all steady states for all times in one single unitary transformation Q(t), which is an ambitious goal. The relation of Q = n |n u n (t)| to Floquet concepts is discussed in the Appendix: While the time-evolution operator W (t) is not the topic here, it can be simply obtained W (t) = Q † (t)e −iHt Q(0). However, it is not possible to construct Q using W . Likewise, the so-called Floquet Hamiltonian [27,46,47,71] H F = Q † (0)HQ(0) can be found using Q. We now proceed to find an explicit expression for Q(t) for the model in Eq. (2). Floquet Bogoliubov solution. The goal is to find a static eigenbasis in the rotating frame, which can be achieved ifH becomes diagonal and time-independent. The interacting model H q (t) in Eq. (2) is defined in left and right oscillator Hilbert spaces χ = L, R, so a static solution must be of the formH where we have used a general Floquet-Bogoliubov ansatz for Q in Eq. (8) with the canonical constraint |γ 1 (t)| 2 − |γ 2 (t)| 2 = 1. The defining relation in Eq. (7) provides differential equations for the time-periodic coefficients γ 1,2 where λ 1 = v F q(1+g 4 ) and λ 2 = v F qg 2 . The relation (9) applies to general TLL, but for the Lieb-Liniger model it simplifies since λ 1 − λ 2 = qv F is constant due to Galilean invariance. Using f ± (t) = e ı∆t (γ 1 (t) ± γ 2 (t)) and Eq. (4) we obtain a Mathieu equation and f + = −ıḟ − /qv F . The solution can be expressed as where a = 4 and C(a, p, τ ), S(a, p, τ ) are even and odd Mathieu functions normalized with C(a, p, 0) = S(a, p, π) = 1. The coefficients c 1,2 are determined by the time-periodicity of steady states |u n (t) and operators β(t), which also fixes the quantization condition for ∆: We use Floquet's theorem to write the solution of Eq. (10) f − (t) = e ıντ P ν (τ ) with P ν (τ ) = P ν (τ + π) [73]. Since γ 1/2 are periodic, we find that the Mathieu characteristic exponent is ν = 2∆/ω, which must be real for stable steady states, just like for the Mathieu stability chart [74] of Paul traps [75]. From the normalization above follows cos(πν) = C[a, p, π], which gives (see Appendix) ∆ = arccos[C(a, p, π)]/T, and c 1 is fixed by |γ 1 | 2 − |γ 2 | 2 = 1. Last but not least, we can use the solutions of γ 1,2 to uniquely define three real time periodic functions θ, φ, r, which parametrize an explicit expression of Q(t) in Eq. In the Appendix it is shown that Q(t) in Eq. (13) and the form of the transformed ground state |u 0 (t) =Q † |0 is provided (see Appendix), which obeys β L,R (t)|u 0 (t) = 0 ∀t. Therefore, from Eqs. (7) and (8) all Floquet modes |u n (t) with ǫ n = (n L + n R + 1)∆ are found by application of (β † L (t)) nL (β † R (t)) nR on |u 0 (t) . Instability regions. Before calculating physical observables, we need to analyse the stability of the differential equations, which may not always have a solution due to the periodicity constraint. In Fig. 2(top) we plot the value of ν = 2∆/ω as a function of rescaled momenta qv F /ω usingρ = −0.6 and amplitude ρ = 0.25. We observe that for certain regions of momenta there are no real solutions for ∆. These "instability regions" will have interesting physical implications as discussed below. The stable regions are shown as a function of amplitude ρ in Fig. 2(bottom) forρ = −0.6. For small ρ the instability regions are equally spaced at integer values ℓ ∈ N corresponding to a = ℓ 2 or ℓ = 2 ity regions therefore correspond to integer multiples of frequency which match twice the interacting dispersion relation ℓω = 2vq ℓ , so the physical cause can be traced to resonant excitations on the linear branches of left movers from −vq ℓ to right movers atvq ℓ and vice versa. As shown in Fig. 3, the region of instabilities also occur for more general TLL models where the restriction g 2 (t) = g 4 (t) in Eq. (4) is lifted [60] and/or contain higher harmonics. A general analytic solution remains elusive, but the corresponding differential equation (9) is still valid, which we have solved numerically by Fourier decomposition for several parameters. Instability regions are always expected since the problem is analogous to forbidden energy regions in a band structure of a periodic potential [71], which is of course generic. In Fig. 3 we show the behavior of ν as a function of qv F /ω for the case that only the g 2 scattering process is periodically modulated in time withρ = −0.6 and ρ = 0.25 in Eq. (4) while g 4 = −0.4. While quantitative changes compared to Fig. 2(top) can be identified, the regions of instabilities are again found at resonant wave vectors. In Figs. 2(top) and 3 we see that ν → 1 near the unstable regions and the ratio c 2 /c 1 in Eq. (12) becomes singular.
To understand the physical significance of the instability regions, it is essential to consider damping. Intrinsic damping is always present in the TLL description due to higher order boson-boson interaction terms [48], which lead to a finite quasiparticle life-time. A corresponding broadening of spectral peaks is seen numerically for finite energies and in finite systems [76,77]. The size of damping is not universal since it depends on microscopic details including the system size, but it can be assumed to be smaller than all other energy scales. In Ref. [74] it was shown that solutions of damped Mathieu equations become always stable for amplitudes below a given threshold. We also find that a finite life-time τ 0 in form of an imaginary energy correction Im λ 1 = −1/τ 0 leads to convergence of instabilities as discussed below.
Results. We are now in the position to calculate physical observables. The main effect of the time periodic driving is the excitation of density waves in the steady state. The number of density excitations b † χq b χq (χ = L or R) in the transformed ground state |u 0 (t) is given by In Fig. 4 we plot the time averageη q . For small q we find thatη q approaches the static limit, but a strong divergence is observed as the instability region around q ℓ is approached. In the inset of Fig. 4 we exemplarily show that a finite life-time τ 0 = 10 4 /v F q turns the divergences ofη q into large maxima around q ℓ . The height of the maxima can be tuned by the product ρτ 0 . A universal physical picture emerges analogous to a resonance catastrophe: A finite life-time has little effect away from resonance, but the resonance response is overwhelmingly large and proportional to τ 0 . If q ℓ = ℓω/2v < q c is in the TLL regime, such maxima will therefore dominate the correlations. We find that q c ∼vm/2 is a good estimate for the cutoff.
It is well known how TLL correlation are calculated [48], which is reviewed in the Appendix for the example of density-density correlations. An overwhelming maximum ofη q will dominate the correlations and lead to longrange density order (see Appendix) For large driving amplitudes ρ the magnitude of the induced density waves can become larger than the background density, which may lead to fragementation into irregular density grains.
Discussion. The three energy scales ω,vq c , and v F q determine the behavior of the system, which undergoes three different regimes as the frequency is changed: 1.) High frequencies: For ω vq c the instability regions are outside the TLL regime, so the physical relevant region is free of resonances. The transformation Q results in a systematic change ofη q shown in Fig. 4, which approaches the static limit as q → 0. The famous power-law correlations [48] are corrected for intermediate distances, but the asymptotic static limit is recovered. 2.) Intermediate frequencies: As the frequency is lowered, the resonant wave-numbers q ℓ = ℓω/2v drop below the cutoff q c into the TLL regime. The number of density wavesη q ℓ becomes very large, dominating the correlations in Eq. (16). Instead of powerlaw correlations, standing density waves at wave-numbers q ℓ become stable throughout the system. 3.) Very low frequencies: For ω ≪vq c extended regions of instability will lead to a large number of excitations and heating, destroying the correlations.
Using typical experimental parameters for a 1D 87 Rb gas from Ref. [51] of n = 6.2×10 6 /m and mg/n = 0.6, we arrive atK ≈ 4 and a cutoff of ω c =vq c ≈ 2π × 1.4kHz in the middle of the trap. Driving the perpendicular confinement with a frequency of ω = 2π × 500Hz results in a resonance at q 1 = ω/2v = 444×10 3 /m. We therefore predict a standing density wave with wavelength λ = 2π/q 1 in the µm range, which is observable in real space with optical methods or an electron beam [51,52].
The confining trapping potential leads to lower local densities n near edges [51,52] and reduced velocities v F = πn/m. Everywhere n agrees with the local density approximation (LDA) of TLL correlations for the local trap potential [51,52]. The trapping potential is therefore turning into an advantage: Instead of changing the frequency ω, different regimes can be reached using the changing density n. As a function of n we knowv = ng/m [63][64][65], which in turn determines the resonant wave-vectors q ℓ = ℓω m/ng/2 and the cutoff q c = √ ngm/2. Therefore, we move into the high frequency regime q ℓ /q c ∝ ω/ng as the density is lowered. Note, that the density wavenumbers q ℓ increase near edges in contrast to Fermionic Friedel density wavenumbers, which decrease with lower densities in a trap [78].
In the proposed experiment, we therefore predict standing density waves at λ ∼ 14µm in the middle of the trap, which become shorter λ ∝ √ n and weaker near the edge. It is an interesting open problem if significant corrections would be observed when going beyond the present LDA analysis for a typical trap size of 120 µm in Ref. [51].
Interesting many-body density excitations have been experimentally observed in driven 1D and 2D systems [79,80]. For 1D elongated bosonic 7 Li gases µm-size density grains emerge at 2π×80Hz driving, which were identified as stable many-body effects [79]. Experimental images show grains that appear smaller and weaker near edges which resemble features predicted above, but in a random pattern [79]. All correlations disappear for very low frequencies ω. A future grain size analysis as a function of ω and n may clarify if there is a relation to TLL density waves in Eq. (16).
Conclusion. We have considered time-periodically driven interacting systems in the steady state, corresponding to generic TLL models in general or the Lieb-Liniger model in particular, which e.g. applies to 1D confined atoms in ultra-cold gas experiments with tunable parameters. As we have shown, this setup is one of the very rare cases where the combination of non-equilibrium steady states with many-body physics can be analyzed in great detail. In particular, we have developed a Floquet-Bogoliubov approach by constructing time-periodic creation and annihilation operators, which solve the eigenvalue equation for the steady state by acting on the entire Floquet space. We also identify regions in frequencymomentum space where damped resonant behavior leads to a large number of density excitations. The known static powerlaw correlations [48] are recovered for large distances ≫v/ω, but for frequencies below the cutof vq c characteristic density waves at integer-spaced resonant wave numbers q ℓ = ℓω/2v will become dominant. We emphasize that the proposed Floquet-Bogoliubov algebra is completely general and can be used to solve any time-periodically driven model with Bogoliubov-type interactions exactly. The explicitly known transformation Q maps all steady states onto a diagonal static oscillator basis for all times, which paves the way for a complete analysis of time-dependent effects in strongly interacting systems using a combination of powerful experimental, analytic, and numerical techniques.

APPENDIX
Here we give details on the Floquet Bogoliubov transformation, its relation to Floquet theory, the explicit form of the transformed ground state state, densitydensity correlations, and the application of Floquet's theorem to Mathieu functions.

Relation of the time-dependent transformation to Floquet theory
The goal is to find all possible steady state solutions |u n (t) = |u n (t + T ) under time-periodic driving at each time t, which are defined by the Floquet eigenvalue equation where ǫ n are real quasi energies. It should be noted that it is not always possible to find steady state solutions, but if they exist they form a complete basis in the original Hilbert space. The underlying Floquet theory has been discussed in a number of review articles [27,46,47,71], where different approaches are presented: By Fourier transforming into frequency space, the eigenvalue problem becomes static in an extended Hilbert space. Different frequency components can be perturbatively decoupled using a Magnus expansion, which is helpful in defining a so-called Floquet Hamiltonian H F . The Floquet Hamiltonian is useful since it determines the quasi-energies and the stroboscopic time evolution. The eigenstates of H F are the steady states |u n (0) at one instant in time only, so for the full time evolution it is necessary to additionally know the micromotion operator U (t) = n |u n (t) u n (0)|, which is in general more difficult.
Our novel approach is now to solve the Floquet eigenvalue problem in one single step by mapping it to a static problem in the original Hilbert spacẽ If solutions to the original problem in Eq. (17) exist the unitary transformation Q can formally always be written as which transforms the entire basis of steady state solutions at each time into a diagonal static basis. This new transformation Q therefore does three things at once: It maps the system to a static problem in the original Hilbert space, it diagonalizes the eigenvalue problem, and it provides the time-dependent steady states for all times. All this is done without using a Fourier transform into an extended Hilbert space. Needless to say, each of the above steps is normally highly non-trivial, so finding such a transformation Q into a diagonal rotating frame is very ambitious indeed. Note, that Q(t) = Q(t + T ) is time periodic, but we need not assume that Q(t) becomes the identity at the initial time or any other time.
The operator Q must therefore not be confused with the time-evolution operator W (20) which can be used to study the time-dependence of a given initial state. In particular, knowing the time evolution cannot be used to construct Q, but the time evolution can always be expressed as Moreover, the Floquet Hamiltonian can be obtained by H F = Q † (0)HQ(0), but again just knowing H F cannot be used to extract the steady states for all times unless Q is known. Finally, also the micromotion operator U (t) = Q † (t)Q(0) and all steady states |u n (t) = Q † (t)|n can be obtained with Q, so such a transformation truely contains a complete solution of the many-body driven system.

Explicit form of the Floquet Bogoliubov transformation
The model of interest can conveniently be expressed in terms of SU(1,1) generators where and λ 1 = v F q(1 + g 4 ) and λ 2 = v F qg 2 are the timeperiodic coupling parameters. For the static case it is known that the transformation U 1 = e r(J+−J−) can be used for diagonalization, using the following relations for (28) For the time-dependent transformation, we need a more general ansatz parametrized in terms of three real time-periodic parameters θ, φ, r Using relations Eqs. (24)-(28) together with gauge transformations, we find that the general time-dependent Bogoliubov transformation can be written as with χ = L, R and With this parametrization the transformed operators Λ = QΛQ † can again be straightforwardly derived from Eqs.
Note, that the three real parameters θ, φ, r give a general one-to-one parametrization of the complex functions γ 1 and γ 2 which obey |γ 1 | 2 − |γ 2 | 2 = 1. The functions γ 1 and γ 2 have been extensively discussed in the paper so the transformation Q is already explicitly known, but what is left to show in the following is that the Hamiltonian in Eq. (22) indeed becomes static and diagonal when using those functions. The defining differential equation is given in Eq. (9) of the paper in terms of γ 1 and γ 2 where ∆ is a real constant which is fixed by the constraint that both γ 1 and γ 2 are periodic as discussed in the paper. In terms of the parametrization θ, φ, r, the differential equations become after multiplying by exp(−i θ±φ 2 ) respectively The imaginary parts of both equations give the same relationṙ The real parts give 0 = (∆ +θ/2 − λ 1 −φ/2) cosh r + λ 2 cos φ sinh r(43) 0 = (∆ +θ/2 + λ 1 +φ/2) sinh r − λ 2 cos φ cosh r(44) For later use we take (43)× sinh r−(44)× cosh r, which gives 0 = −(λ 1 +φ/2) sinh 2r + λ 2 cos φ cosh 2r (45) Likewise (44)× sinh r−(43)× cosh r gives ∆ = −θ/2 + (λ 1 +φ/2) cosh 2r − λ 2 cos φ sinh 2r (46) We now turn to identify the different parts in the transformed Hamiltoniañ Collecting all the terms ofH from Eqs. (35)-(39) we find that the prefactor of the diagonal part 2J 0 reads (λ 1 +φ 2 ) cosh 2r − λ 2 cos φ sinh 2r −θ 2 which is exactly ∆ according to Eq. (46) and therefore time-independent. The prefactor of the off-diagonal part e iθ J + is given by Using Eq. (42) for the imaginary part and Eq. (45) for the real part, we see that this expression is indeed zero, so that we have shown that the model in Eq. (22) transforms toH where the constant ∆ is determined by the constraint of periodicity and Floquet's theorem as described in the text.
(54) With e iφ |γ 2 |/|γ 1 | = γ 2 /γ 1 we finally find an explicit expression for the transformed ground state It is important to note that while the form of state (55) is similar to the results of a static Bogoliubov transformation [68] here all parameters are time-dependent. Using the transformation Q the state (55) solves the Floquet Eq. (5) in the main article with ǫ 0 = ∆. Moreover, we can show explicitly that the transformed ground state |u 0 (t) obeys condition Eq. (51) by applying β L (t) = γ 1 (t)b L + γ 2 (t)b † R to Eq. (55), which reads As the first bracket in (56) vanishes trivially, the state (55) is indeed the ground state of the β L (t) operator obeying Eq. (51) and analogously also for β R (t). This is an important result, as |u 0 (t) serves as base case for generating the entire set of steady states |u n (t) by application of (β † L (t)) nL (β † R (t)) nR using Eq. (7) in the main article.

Correlation functions
It is well known how to calculate correlation functions of physical operators in terms of the diagonal boson modelH [48,63,64]. Of particular interest for ultra-cold gases is the density-density correlation, which we will consider here to exemplify the calculation. The fluctuating density is given in terms of the bosonic field n(x) = ∂ x φ(x)/π, which has the mode expansion [48,63,64] For the density-density correlation function we find in the transformed ground state |u 0 (t) = Q † (t)|0 where we have used Eq. (31). If the parameters γ 1,2 are constant we recover the known asymptotic powerlaw behavior 1 2π 2 |γ 1 + γ * 2 | 2 /|x − y| 2 [48,63,64]. However, if a resonance q ℓ = ℓω/2v is part of the linear TLL regime, the parameters γ 1,2 will become very large as discussed in the main article. Therefore, the sum in Eq. (58) will be dominated by the corresponding instability region, leading to a long-range density order of the form u 0 |n(x)n(y)|u 0 ∝ cos q ℓ (x − y).

Floquet solution in terms of Mathieu functions
The solution of the Mathieu equation is usually discussed in terms of even and odd solutions, known respectively as Mathieu cosine C and Mathieu sine S functions. A general solution can be therefore written as with τ = ωt/2. Floquet's theorem states that the solutions of a time-periodic differential equation can always be written in the form with P ν (τ ) = P ν (τ ± π). We want to use the quantum number ν, which is commonly referred to as Mathieu characteristic exponent. Therefore, in this section we clarify the relation between the latter and the Mathieu functions. Comparing Eqs. (61) and (62) and employing the periodicity of P ν (τ ), we get the following relation c 1 C(a, p, τ ) + c 2 S(a, p, τ ) = e ∓ıνπ (c 1 C(a, p, τ ± π) + c 2 S(a, p, τ ± π)) .

Static Bogoliubov transformation
For the time-independent case, Hamiltonian (2) in the main article can be expressed as For the sake of simplicity, in the following we will drop the index q. We observe that J i (i = 0, ±) form a su(1, 1) algebra; therefore the Tomonaga-Luttinger Hamiltonian can be diagonalized by the Schrieffer-Wolff transformation obtained through the unitary operator The transformed HamiltonianH is diagonal in the old bosonic basis with ∆ = v F q (1 + g 4 ) 2 − g 2 2 if tanh (2ϑ) = tanh (2θ) ≡ g 2 /(1 + g 4 ), implying cosh(2θ) = v F q(1 + g 4 )/∆, sinh(2θ) = v F qg 2 /∆. Notice that we have used a passive transformation, where the operators are rotated, while the states defined by the original creation and annihilation operators stay the same.
The classical equations of motions also give insight to the role of damping, which generally reduces the regions of instability of the Mathieu equation. In fact, in [74] it is shown that the presence of linear damping pushes the instability zone upwards in Fig. 2, so that below a critical value of amplitude ρ stable solution are always possible and no instabilities occur. In [74] it is also shown that nonlinear effects can generate subharmonic stable motions. These cubic terms correspond to band curvature in the original band structure, so in a real system the nonlinearity of the band further stabilizes the system. While a large number of density waves is still expected to occur at the critical q−values, the sum over all momenta becomes well defined leading to the predicted density order of scenario 2) in the main paper.

Floquet theory
Analogously to the Bloch theorem, the Floquet theorem asserts that the Schrödinger equation for a timeperiodic Hamiltonian admits steady-state solutions of the form where the modes |u(t) = |u(t + T ) inherit the periodicity from the Hamiltonian, and the quantity ǫ is the so-called Floquet quasienergy. Indeed the steadystate Schrödinger equation can be recasted in the form of an eigenvalue equation for the quasienergy operator H = H(t) − ı∂ t in the extended Hilbert space generated by the product of the state space of the quantum system and the space of square-integrable T -periodic functions: H |u(t) = ǫ|u(t) .