Floquet theory for temporal correlations and spectra in time-periodic open quantum systems. Application to squeezed parametric oscillation beyond the rotating-wave approximation

Open quantum systems can display periodic dynamics at the classical level either due to external periodic modulations or to self-pulsing phenomena typically following a Hopf bifurcation. In both cases, the quantum fluctuations around classical solutions do not reach a quantum-statistical stationary state, which prevents adopting the simple and reliable methods used for stationary quantum systems. Here we put forward a general and efficient method to compute two-time correlations and corresponding spectral densities of time-periodic open quantum systems within the usual linearized (Gaussian) approximation for their dynamics. Using Floquet theory we show how the quantum Langevin equations for the fluctuations can be efficiently integrated by partitioning the time domain into one-period duration intervals, and relating the properties of each period to the first one. Spectral densities, like squeezing spectra, are computed similarly, now in a two-dimensional temporal domain that is treated as a chessboard with one-period x one-period cells. This technique avoids cumulative numerical errors as well as efficiently saves computational time. As an illustration of the method, we analyze the quantum fluctuations of a damped parametrically-driven oscillator (degenerate parametric oscillator) below threshold and far away from rotating-wave approximation conditions, which is a relevant scenario for modern low-frequency quantum oscillators. Our method reveals that the squeezing properties of such devices are quite robust against the amplitude of the modulation or the low quality of the oscillator, although optimal squeezing can appear for parameters that are far from the ones predicted within the rotating-wave approximation.

Open quantum systems can display periodic dynamics at the classical level either due to external periodic modulations or to self-pulsing phenomena typically following a Hopf bifurcation. In both cases, the quantum fluctuations around classical solutions do not reach a quantum-statistical stationary state, which prevents adopting the simple and reliable methods used for stationary quantum systems. Here we put forward a general and efficient method to compute two-time correlations and corresponding spectral densities of time-periodic open quantum systems within the usual linearized (Gaussian) approximation for their dynamics. Using Floquet theory we show how the quantum Langevin equations for the fluctuations can be efficiently integrated by partitioning the time domain into one-period duration intervals, and relating the properties of each period to the first one. Spectral densities, like squeezing spectra, are computed similarly, now in a two-dimensional temporal domain that is treated as a chessboard with one-period × one-period cells. This technique avoids cumulative numerical errors as well as efficiently saves computational time. As an illustration of the method, we analyze the quantum fluctuations of a damped parametrically-driven oscillator (degenerate parametric oscillator) below threshold and far away from rotating-wave approximation conditions, which is a relevant scenario for modern low-frequency quantum oscillators. Our method reveals that the squeezing properties of such devices are quite robust against the amplitude of the modulation or the low quality of the oscillator, although optimal squeezing can appear for parameters that are far from the ones predicted within the rotating-wave approximation.

I. INTRODUCTION
In recent years, with the development of new quantum technologies, more complex protocols to control and manipulate quantum devices have been proposed. Usually, those devices are made up of nearly isolated quantum systems (atoms, solid-state defects, superconducting circuits, mechanical elements, etc) that interact coherently with the electromagnetic field (at optical or microwave frequencies) via an input port, where a driving is applied, and an output port, where the detection is performed. Also, the considered quantum system can interact with its environment, usually leading to an incoherent exchange of excitations which manifests as noise (thermal, electronic, etc). A proper engineering of all these processes is key to design quantum technologies for applications in quantum computation, simulation, communication, and metrology.
Taking all these things into consideration, it is clear that periodic modulations play a major role in many different fields of contemporary quantum physics. In this work, we will concentrate on periodically-driven open quantum-optical systems. There are two standard mathematical descriptions of such quantum systems: (i) via a set of coupled quantum Langevin equations, which are Heisenberg (differential) equations for the operators supplemented by dissipation terms and input quantum noises, or (ii) via a master equation for the density operator, which consists on the von Neumann equation for the state, to which Lindblad terms accounting for irreversible quantum jumps are added. Let us remark that master equations can be mapped to a set of stochastic Langevin equations by resorting to phasespace representations of the density operator (like the Wigner function or, more robustly, the positive P distribution [85][86][87]). Hence, in both approaches a set of Langevin equations can be ultimately obtained, which provide a route towards the numerical analysis of dynamical features.
Due to the generally nonlinear nature of such equations, exact solutions are hard or impossible to find, except for very specific cases. Analytic or semi-analytic insight is usually gained by using the so-called standard linearization technique, which typically provides sensible results, except close to phase transitions or in the presence of spontaneous symmetry breaking, which nevertheless can be treated with suitable generalizations of such technique [50,[88][89][90][91][92][93][94]. Within this approach, one considers small quantum fluctuations around a reference classical state, leading to a linear system of Langevin equations for the fluctuations, which is easily handled only if the classical reference is time independent. This can occur in problems involving a constant pump, like the laser, or even in the presence of a monochromatic drive, as in optical parametric oscillators [88] and optomechanical cavities [95] (but only if a rotatingwave approximation can be invoked). However, even in such cases, the stationary classical solutions can become unstable (e.g. via a Hopf bifurcation), and spontaneous oscillations can emerge in the classical dynamics, leading to nontrivial linearized Langevin equations for the fluctuations, which in particular will contain now time-periodic coefficients that make them hard to treat. The same happens if the drive contains more than one frequency (or if the rotating-wave approximation cannot be used), in which case it is in general impossible to obtain time-independent classical states.
For this type of linear Langevin equations with time-periodic coefficients, common strategies are based on Fourier expansions [4,96,97]. Recently, however, we put forward a more compact approach based on the Floquet theorem [50], which transforms linear homogeneous differential equations with periodic coefficients into equivalent equations with constant coefficients. There, however, we focused on the determination of the asymptotic covariance matrix for long times, and its connection with the steady state of the master equation after diffusion around the limit cycle has taken over. In this work we go deeper into the general Floquet method for linearized systems, in particular using it to develop an efficient method for the computation of experimentally-relevant quantities such as two-time correlation functions and the corresponding spectral densities. Specifically, we show how these quantities can be evaluated just from knowledge of the behavior of the system during a single period, which is crucial in order to avoid significant errors in the computation of such observables: a large numerical effort can be employed at a low cost to perform highly precise integrations along one period, and then propagate that information algebraically over the long term.
As a practical example, we use the theory to analyze the squeezing properties of degenerate parametric oscillators beyond the rotating-wave approximation, which has become a timely issue, since such a model can be implemented nowadays in low-frequency superconducting oscillators well within the quantum regime [29]. Our results support the robustness of squeezing against the modulation amplitude or the bad quality of the oscillator. Moreover, we show that once counter-rotating terms are incorporated, optimal squeezing is achieved for modulation amplitudes below the oscillation instability, contrary to the rotating-wave predictions, for which optimal squeezing always occurs at the instability.
The manuscript is organized as follows. In Sec. II we briefly review the description of open quantum systems via linearized Langevin equations, and introduce the Floquet-based method for the determination of their solutions. In Secs. III-V we use the solutions to manipulate two-time correlations and the corresponding spectral densities, producing compact expressions solely based on dynamics over a single period. Finally, in Sec. VI we apply the theory to degenerate parametric oscillation beyond the rotating-wave approximation.

II. LINEARIZATION IN TIME-PERIODIC OPEN QUANTUM SYSTEMS: FLOQUET THEORY
Consider an open quantum system furnished with a set of D operatorsr = (r 1 , ...,r D ) T , where the symbol T denotes transposition. The system evolves according to its own dynamics as well as to interactions with its environment, which in general is composed of several reservoirs with which the system exchanges energy. In the Heisenberg picture, which we adopt, and assuming standard Markovian conditions, the system operators evolve generically according to some quantum Langevin equations dr dt = A(λ;r) + B(λ;r)ξ(t).
Here A = (A 1 , ..., A D ) T accounts for the deterministic, Hamiltonian or not, part of the dynamics and depends on the system operators and also on a set of control parameters generically denoted by λ (e.g. the amplitude and frequency of a driving field). These can be time dependent thereby inducing periodic dynamics in the classical limit. The fluctuations fed by the reservoirs into the system, responsible for irreversible quantum jumps, enter the dynamics through a noise term. With full generality, we write it as a (D × N ) matrix B (that might depend as well on the control parameters and the system operators) acting on a vectorξ(t) = (ξ 1 (t), ...,ξ N (t)) T composed of Gaussian white noises with ξ n (t) = 0 and two-time correlators ξ m (t)ξ n (t ) = G mn δ(t − t ), which define a noise-correlation matrix G. Note that the number N of independent noises needs not equal D (e.g. in an optical cavity there are several input vacua per mode, even if in many instances one can ignore all but one). Let us remark that the stochastic Langevin equations naturally obtained from the Schrödinger picture through phase-space representations such as the positive P [85,86] have the same form as Eq. (1), but replacing operators by suitable stochastic variables. Hence, the theory that we are going to put forward applies also to such an alternative, but common approach to open quantum systems.
Note that throughout this work we use bold fonts for vectors, e.g. r, which by default correspond to columns with components denoted by r m , so that r T corresponds to a row vector; also, a dagger will denote the conjugate-transpose as usual, e.g. r † := r * T . On the other hand, we use calligraphic fonts for matrices, e.g. G, whose components we denote by G mn .
We follow the standard linearization procedure that starts by splitting each operatorr m as its mean field r m plus a fluctuationx m , i.e.r = r +x. In the semiclassical approximation that is commonly adopted, the mean field r , to be denoted as r, is ruled by the dynamical system of equations dr/dt = A (λ; r), obtained from Eq. (1) by substituting operators by their mean values and ignoring noises. These correspond to the classical limit, and we are here interested in the case where such classical dynamics is periodic, i.e. r = R(t), with R (t + T ) = R(t) some periodic function with period T . This can happen either when the control parameter λ is periodically modulated in time, or following a dynamical (typically Hopf) bifurcation occurring at some critical value λ = λ osc which marks the onset of self-sustained oscillations (see [50] for a detailed example).
The dynamics of the fluctuations is governed by the original quantum Langevin equations (1), which after linearization with respect to fluctuations and noises are written as where we denote B[λ; R(t)] simply by B(t), and the (D × D)-matrix L is the Jacobian of the classical dynamical equations, with elements L mn (t) = ∂A m (λ; r) /∂r n | r=R(t) . The Jacobian L depends on the parameters and on the classical solution, and thus it is explicitly T -periodic, L(t + T ) = L(t), as is the matrix B(t). Hence (2) is a non-autonomous dynamical system of linear equations, which prevents its analytical solving. However application of Floquet theory allows us to transform Eq. (2) into a system with a time-independent Jacobian, which is more amenable to analytical or semi-analytical treatments. Let us review here the procedure, which we will exploit throughout the rest of the work. We start by defining the principal fundamental matrix F(t) through the initial-value problem where I D×D is the (D × D) identity matrix. Note that the choice of 0 as the initial time is arbitrary, and any other choice connects with it by a trivial time shift. Next, we construct a constant matrix M through which serves to decompose the fundamental matrix in its so-called Floquet normal form, where P(t) is a T -periodic invertible matrix. Defining a transformed fluctuation vector the non-autonomous Eq. (2) turns into the autonomous one This constitutes an example of Floquet's theorem.
The system of equations (7) can be formally solved in terms of the eigensystem of matrix M. Let us denote by S the (D × D) matrix that diagonalizes M through the similarity transformation The eigenvalues {µ α } D α=1 are known as Floquet (or characteristic) exponents. Note that in previous works [50] we have used a slightly less compact notation, where we defined the set of right and left eigenvectors of M, satisfying Mv α = µ α v α , w † α M = µ α w † α , and orthonormality relations w † α v β = δ αβ . These two notations are connected by It proves convenient to define the auxiliary matrix Upon multiplying (7) by S −1 from the left, and defining the projectionŝ we get which are a set of decoupled linear equations for the components of c, whose formal solution can be put aŝ Here we assumed that all the eigenvalues µ α have negative real part (i.e., the analyzed semiclassical state is linearly stable), hence the integral (13) is bounded. Expressions (3), (4), (5), (8), and (13) constitute the basis of our analysis, as they allow computing the fluctuation vectorx in terms of the noise integrals that depend only on the auxiliary matrix K(t) and the Floquet exponents {µ α } D α=1 . To conclude this section we note that computing matrix M is not required at any step. Instead, we can use the so-called monodromy matrix F(T ), which is diagonalized by the same similarity transformation (8), and possesses eigenvalues {φ α } D α=1 related to the Floquet exponents by µ α T = ln φ α .

III. COMPUTATION OF TWO-TIME CORRELATIONS
Our goal is the computation of physical quantities related to the quantum fluctuations of the system x around the stable, periodic semiclassical solution r = R(t). Within the linearized approximation that we are using, which is equivalent to assuming the state to be Gaussian [93,101], the most general quantities that one can consider are two-time correlations, since for Gaussian distributions any higher-order correlation can be reduced to products of two-time ones. Hence, the most general correlators we want to compute are where we remind thatx is a column vector, so X is a matrix. As a first result of this work, we provide here a simple expression for this two-time correlation matrix that exploits the periodic nature of the problem. We start by using (14) to rewrite it as where are elementary correlations that we work out in Appendix A. We relegate the technical derivations to that appendix, and summarize here only the final compact expressions. Note first that the projected noises (11b) are delta correlated as with a projected-noise correlation matrix With this definition at hand, we show in Appendix A that the correlation matrix (17) can be worked out to yield the components where with Equations (19)- (21) are the first main result of this work as they allow us to compute any two-time correlation in terms of integrals of functions evaluated just in the interval t ∈ [0, T ]. The importance of this result emerges when long measurement times are involved, as those required for the computation of spectral densities (see next section), because, apart from being numerically demanding, the errors accumulated in finding the Floquet matrix P(t) at long times can be large enough to invalidate the results.
Note that, for numerical purposes, it is typically more efficient to evaluate ν αβ (τ ) from the equivalent initial-value problemν rather than from the integral (24).

IV. COMPUTATION OF SPECTRAL DENSITIES
Another important tool for characterizing quantum fluctuations are the spectral densities associated to two-time correlations. A relevant example is the light squeezing spectrum, which is the spectral variance of the (quantum) noise carried by a light beam, and can be measured experimentally via balanced homodyne detection [88,98] or alternative correlation measurements [99]. In the usual stationary case, i.e. when two-time correlations are a function only of the two-time difference, these densities are just plain Fourier transforms. However, when such correlations are not stationary one has to use a different definition in order to match the experimentally detected spectral density [88,98], namely where O(t, t ) is the considered two-time correlation and T d is the detection time. In general the measurable densities will be linear combinations of S(ω) and S(−ω), as we will see later through a practical example. We then consider spectral densities of the form (23), being P α (t) and P β (t ) generic T -periodic functions whose meaning is as follows. When such functions are chosen as K mα (t) and K nβ (t ), respectively, and summing over α and β, one can compute spectral densities corresponding to the correlations X mn (t, t ), see (16). The choice P α (t) = P β (t ) = 1 is also interesting as it provides the spectral densities corresponding to the elementary correlations C αβ (t, t ), which in some cases are proportional to measurable quadratures [88][89][90][91][92]. Finally, whenr is formed of annihilation and creation operators, with the choice P α (t) = Λ α (t)K mα (t) and P β (t ) = Λ β (t )K nβ (t ), (24) allows the computation of spectral densities corresponding to homodyne-detection experiments when the local oscillator is a T -periodic function [100], in which case Λ α (t) and Λ β (t ) are proportional to the amplitude (or its complex conjugate) of that local oscillator. At first sight it seems easy to solve the problem once the correlation functions C αβ (t, t ) have been expressed in Eqs. (19) in terms of the first period. However, for long measurement times, as realistically needed, the integrals are still numerically demanding and can carry important numerical errors. In order to avoid this, we have worked out Eq. (24) by exploiting the properties of the integral's kernel, and managed to simplify it into a few integrals defined only over a single period. As we did in the previous section, we relegate the technical derivations to Appendix B, presenting here the final result. Moreover, we focus on the common situation of a long detection time that contains very many periods, that is, T d T . In this limit, as proven in Appendix B, the general spectral density (24) is simplified as where Hence, again we have been able to write spectral densities in terms of first-period objects only, which comes with all the numerical benefits that we highlighted above. Hence, this is the second main result of our work, which provides a compact way of evaluating arbitrary spectral densities in periodic systems from knowledge of the Floquet eigensystem over a single period. Similarly to what we did in the previous section with ν αβ (t) in Eq. (22), it is useful for numerical efficiency to find the integrals defined above from their equivalent differential equations. In the case of I αβ (ω) and I αβ (ω), both are of the integral form I = T 0 dtf (t) T 0 dt h(t ). Hence, defining two independent initial-value problemṡ we get I = F (T )H(T ). On the other hand, I αβ (ω) and I αβ (ω) are of the nested type I = T 0 dtf (t) T t dt h(t ), which makes their differential form a bit more intricate, but equally efficient from a numerical standpoint. In this case, we first solve the initial-value problemḢ backwards in time in the domain t ∈ [0, T ], and next the initial-value probleṁ so that I = F (T ).

V. CROSS-CORRELATIONS AND CROSS-SPECTRAL DENSITIES WITH THE NOISE
In the previous sections we focused on the two-time correlations and spectral densities of the variablesx (or, equivalently, the projectionsĉ). However, in many situations one also needs objects related to the cross-correlations between the variables and the noisesξ. A most prominent case is related to the evaluation of quantities related to the field leaking out of the open system by using input-output relations. We will showcase this in the practical example that we consider in the next section. This section is then devoted to provide compact expressions for these type of cross-correlations and spectral densities. Again we make all technical derivations in Appendix C, and offer here just the final results.
We start by providing the two-time cross-correlators between the projectionsĉ and the noisesξ, which are easily worked out as and C (ξc) where we have defined the matrices χ (cξ) (t) = K −1 (t)B(t)G and χ (ξc) (t) = GB T (t)K −1T (t). The corresponding spectral densities, defined, respectively, by replacing C αβ (t, t ) in (24) by C with As an application of the method developed above, we consider now the degenerate parametric oscillator as an example. In essence, it consists of a lossy quantum-mechanical harmonic oscillator whose frequency is modulated periodically at twice its natural frequency. This model serves as the canonical one for the study of quantum squeezing, and has been traditionally explored experimentally with nonlinear optical cavities [88]. Since in this context accessible modulation amplitudes are much smaller than optical frequencies, one can perform a rotating-wave approximation that maps the problem to an effective time-independent one. In contrast, modern implementations based on lowfrequency oscillators (e.g., in superconducting circuits [29] or optomechanical devices [30]) allow to explore the regime where the modulation amplitudes are a significant fraction of the oscillation frequencies. Under such conditions, the predictions derived within the rotating-wave approximation require corrections, and it is our purpose to study these here.

A. The model
Consider an oscillator of mass m and intrinsic frequency Ω, with positionq and momentump, such that [q,p] = i . We can describe the modulated case by the Hamiltonian with (normalized) modulation amplitude ε. Let us write the position and momentum in terms of annihilation and creation operators asq with [â,â † ] = 1, and move to a picture rotating at frequency Ω, where the system evolves according to the Hamiltoniañ where we have removed terms that are not operators, and therefore play no role in the dynamics. In the limit ε 1, one can invoke the rotating-wave approximation, which allows neglecting the rapidly-oscillating terms, leading to the time-independent HamiltonianH ≈ i Ωε(â †2 −â 2 )/8. This is the usual Hamiltonian employed to analyze degenerate parametric oscillators. Here, in contrast, we use the theory developed in the previous sections to study the full Hamiltonian (35).
In order to include losses, we consider the interaction between the oscillator and a bosonic environment at zero temperature. Assuming that the standard Born-Markov approximation holds, one can integrate out the environment obtaining the quantum Langevin equation where we have defined the normalized modulation amplitude σ = εΩ/4γ, andâ in (t) is the so-called input operator, which is Gaussian and characterized by the following statistical properties: Let us also remark that the operators in (36) are slowly-varying or rotating-picture operators, but we keep the same notation as before for simplicity, and because these are actually the operators that homodyne detection is sensitive to, so they are the ones we will use to compute the relevant spectral densities. It is convenient to introduce dimensionless timet = γt, which we adopt in the following but removing the tilde for notational simplicity. Let us further define the vectorâ = (â,â † ) T andâ in = (â in ,â † in ) T , from which we build the quadrature vectorsx = Tâ andx in = Tâ in / √ γ, with T = 1 1 In terms of quadratures, Eq. (36) is then written as the linear system < l a t e x i t s h a 1 _ b a s e 6 4 = " 4 q l 6 M t Q m O q g y 1 I Q 2 e Y D 6 d j B 7 P 6 k = " t 7 + Q f n w q K X j V B H a J D G P V S f E m n I m a d M w w 2 k n U R S L k N N 2 O L 6 d + e 0 n q j S L 5 Y O Z J D Q Q e C h Z x A g 2 V n r s a T Y U + M Z z a / 1 y x X O 9 O d A q 8 X N S g R y N f v m r N 4 h J K q g 0 h G O t u 7 6 X m C D D y j D C 6 b T U S z V N M B n j I e 1 a K r G g O s j m F 0 / R m V U G K I q V L W n Q X P 0 9 k W G h 9 U S E t l N g M 9 L L 3 k z 8 z + u m J r o K M i a T 1 F B J F o u i l C M T o 9 n 7 a M A U J Y Z P L M F E M X s r I i O s M D E 2 p J I N w V 9 + e Z W 0 L l y / 6 t b u q 5 X 6 d R 5 H E U 7 g F M 7 B h 0 u o w x 0 0 o A k E J D z D K 7 w 5 2 n l x 3 p 2 P R W v B y W e O 4 Q + c z x 9 2 q p A b < / l a t e x i t > = 0.7 < l a t e x i t s h a 1 _ b a s e 6 4 = " w 1 S 4 + v f j N s U n u W + n H Z 1 e 5 I V j 1 t c = " t 7 + Q f n w q K X j V B H a J D G P V S f E m n I m a d M w w 2 k n U R S L k N N 2 O L 6 d + e 0 n q j S L 5 Y O Z J D Q Q e C h Z x A g 2 V n r s a T Y U + M Z z a / 1 y x X O 9 O d A q 8 X N S g R y N f v m r N 4 h J K q g 0 h G O t u 7 6 X m C D D y j D C 6 b T U S z V N M B n j I e 1 a K r G g O s j m F 0 / R m V U G K I q V L W n Q X P 0 9 k W G h 9 U S E t l N g M 9 L L 3 k z 8 z + u m J r o K M i a T 1 F B J F o u i l C M T o 9 n 7 a M A U J Y Z P L M F E M X s r I i O s M D E 2 p J I N w V 9 + e Z W 0 L l y / 6 l 7 e V y v 1 6 z y O I p z A K Z y D D z W o w x 0 0 o A k E J D z D K 7 w 5 2 n l x 3 p 2 P R W v B y W e O 4 Q + c z x 9 5 s p A d < / l a t e x i t > = 1.04 < l a t e x i t s h a 1 _ b a s e 6 4 = " 3 q d O n H 8 u y F m X p Q K b s D z v p 3 S 2 v r G 5 l Z 5 u 7 K z u 7 d / U D 0 8 a p s k 0 5 S 1 a C I S 3 Y 2 I Y Y I r 1 g I O g n V T z Y i M B O t E 4 7 u Z 3 3 l i 2 v B E P c I k Z a E k Q 8 V j T g l Y K e g Z P p T k 1 n e 9 e r 9 a 8 1

0.8
< l a t e x i t s h a 1 _ b a s e 6 4 = " z X n g e K L W z L 9 t 9 4 l / W / w u 9 6      , and Jacobian

< l a t e x i t s h a 1 _ b a s e 6 4 = " H 8 M i h C P P 7 8 e J e k Q q u 9 N N 6 s v J 7 Q k = " > A A A B 6 n i c b V A 9 S w N B E J 2 L X z F + R S 1 t F o M Q m 3 A X F C 0 s A j a W E c 0 H J E f Y 2 8 w l S / b 2 j t 0 9 I Y T 8 B B s L R W z 9 R X b + G z f J F Z r 4 Y O D x 3 g w z 8 4 J E c G 1 c 9 9 v J r a 1 v b G 7 l t w s 7 u 3 v 7 B 8 X D o 6 a O U 8 W w w W I R q 3 Z A N Q o u s W G 4 E d h O F N I o E N g K R r c z v / W E S v N Y P p p x g n 5 E B 5 K H n F F j p Y c y O + 8 V S 2 7 F n Y O s E i 8 j J c h Q 7 x W / u v 2 Y p R F K w w T V u u O 5 i f E n V B n O B E 4 L 3 V R j Q t m I D r B j q a Q R a n 8 y P 3 V K z q z S J 2 G s b E l D 5 u r v i Q m N t B 5 H g e 2 M q B n q Z W 8 m / u d 1 U h N e + x M u k 9 S g Z I t F Y S q I i c n s b 9 L n C p k R Y 0 s o U 9 z e S t i Q K s q M T a d g Q / C W X 1 4 l z W r F u 6 h c 3 l d L t Z s s j j y c w C m U w Y M r q M E d 1 K E B D A b w D K / w 5 g j n x X l 3 P h a t O S e b O Y Y / c D 5 / A I v f j U w = < / l a t e x i t > (b) < l a t e x i t s h a 1 _ b a s e 6 4 = " g d T + T O C j q J a y C 2 1 i N u B Q L y K Y + S s = " > A A A B 6 n i c b V A 9 S w N B E J 2 L X z F + R S 1 t F o M Q m 3 A X F C 0 s A j a W E c 0 H J E f Y 2 + w l S / b 2 j t 0 5 I Y T 8 B B s L R W z 9 R X b + G z f J F Z r 4 Y O D x 3 g w z 8 4 J E C o O u + + 3 k 1 t Y 3 N r f y 2 4 W d 3 b 3 9 g + L h U d P E q W a 8 w W I Z 6 3 Z A D Z d C 8 Q Y K l L y d a E 6 j Q P J W M L q d + a 0 n r o 2 I 1 S O O E + 5 H d K B E K B h F K z 2 U g / N e s e R W 3 D n I K v E y U o I M 9 V 7 x q 9 u P W R p x h U x S Y z q e m 6 A / o R o F k 3 x a 6 K a G J 5 S N 6 I B 3 L F U 0 4 s a f z E + d k j O r 9 E k Y a 1 s K y V z 9 P T G h k T H j K L C d E c W h W f Z m 4 n 9 e J 8 X w 2 p 8 I l a T I F V s s C l N J M C a z v 0 l f a M 5 Q j i 2 h T A t 7 K 2 F D q i l D m 0 7 B h u A t v 7 x K m t W K d 1 G 5 v K + W a j d Z H H k 4 g V M o g w d X U I M 7 q E M D G A z g G V 7 h z Z H O i / P u f C x a c 0 4 2 c w x / 4 H z + A I p a j U s = < / l a t e x i t > (a) < l a t e x i t s h a 1 _ b a s e 6 4 = " 5 B h L E N 4 A H C Q 2 r p G A n S 9 c 9 f T Y f a A = " > A A A B 6 n i c b V A 9 S w N B E J 2 L X z F + R S 1 t F o M Q m 3 A X F C 0 s A j a W E c 0 H J E e Y 2 + w l S / b 2 j t 0 9 I Y T 8 B B s L R W z 9 R X b + G z f J F Z r 4 Y O D x 3 g w z 8 4 J E c G 1 c 9 9 v J r a 1 v b G 7 l t w s 7 u 3 v 7 B 8 X D o 6 a O U 0 V Z g 8 Y i V u 0 A N R N c s o b h R r B 2 o h h G g W C t Y H Q 7 8 1 t P T G k e y 0 c z T p g f 4 U D y k F M 0 V n o o 4 3 m v W H I r 7 h x k l X g Z K U G G e q / 4 1 e 3 H N I 2 Y N F S g 1 h 3 P T Y w / Q W U 4 F W x a 6 K a a J U h H O G A d S y V G T P u T + a l T c m a V P g l j Z U s a M l d / T 0 w w 0 n o c B b Y z Q j P U y 9 5 M / M / r p C a 8 9 i d c J q l h k i 4 W h a k g J i a z v 0 m f K 0 a N G F u C V H F 7 K 6 F D V E i N T a d g Q / C W X 1 4 l z W r F u 6 h c 3 l d L t Z s s j j y c w C m U w Y M r q M E d 1 K E B F A b w D K / w 5 g j n x X l 3 P h a t O S e b O Y Y / c D 5 / A I j V j U o = < / l a t e x i t >
with periodicity T = π/Q in terms of the normalized frequency Q := Ω/γ, which coincides with the resonator quality factor. Note that we have obtained a linear system of equations directly, because our initial Hamiltonian (33) was quadratic. This is however an idealization that works only in a limited range of parameters, whose breakdown is signaled by the equations becoming unstable. For example, within the common rotating-wave approximation valid when ε = 4σ/Q 1 as mentioned above, and obtained from (38) by neglecting the oscillatory terms L non-RWA (t) in (39), the Jacobian takes the diagonal form L RWA with eigenvalues −(1 ± σ). Hence, this idealized linear picture is valid only for σ < 1. Beyond such point, the modulation cannot be treated as a given ε sin(2Ωt) term anymore, and needs a dynamical treatment of its own, for example as a dynamical variable that feels some backaction from the oscillator (known as pump depletion in optical implementations). Similar behavior is to be expected beyond the rotating-wave approximation, but this time signaled by the real part of some Floquet exponent µ α becoming positive.

B. Spectral covariance matrix
In order to understand the squeezing properties of this system, we will consider the spectral covariance matrix, which is the standard object recovered via homodyne detection of the excitations that leak out of the oscillator (e.g., the light exiting the cavity through a partially transmissive mirror in a degenerate parametric oscillator). Introducing the output operatorx

l R F f S Z 1 k k V 3 Q d Z G U p u Q = " > A A A B 6 n i c d V D J S g N B E K 1 x j X G L e v T S G A R P Q 8 8 4 c T k I A S 8 e E z Q L J E P o 6 f Q k T X o W u n u E E P I J X j w o 4 t U v 8 u b f 2 E l G U N E H B Y / 3 q q i q F 6 S C K 4 3 x h 7 W 0 v L K 6 t l 7 Y K G 5 u b e / s l v b 2 m y r J J G U N m o h E t g O i m O A x a 2 i u B W u n k p E o E K w V j K 5 n f u u e S c W T + E 6 P U + Z H Z B D z k F O i j X R b v 3 J 7 p T K 2 X e x g f I m w j e c w 5 N S t n H k e c n K l D D l q v d J 7 t 5 / Q L G K x p o I o 1 X F w q v 0 J k Z p T w a b F b q Z Y S u i I D F j H 0 J h E T P m T + a l T d G y U P g o T a S r W a K 5 + n 5 i Q S K l x F J j O i O i h + u 3 N x L + 8 T q b D C 3 / C 4 z T T L K a L R W E m k E 7 Q 7 G / U 5 5 J R L c a G E C q 5 u R X R I Z G E a p N O 0 Y T w 9 S n 6 n z R d 2 / H s S t 0 r V 7 0 8 j g I c w h G c g A P n U I U b q E E D K A z g A Z 7 g 2 R L W o / V i v S 5 a l 6 x 8 5 g B + w H r 7 B N A K j X Y = < / l a t e x i t > Q = 3 < l a t e x i t s h a 1 _ b a s e 6 4 = " G I O v C e / t W N X R r P m r 8 a E 6 M h C E B H Q = " > A A A B 6 n i c d V D L S g N B E O y N r x h f U Y 9 e B o P g a Z n N R h M P Q s C L x w R N D C R L m J 3 M J k N m H 8 z M C i H k E 7 x 4 U M S r X + T N v 3 G S r K C i B Q 1 F V T f d X X 4 i u N I Y f 1 i 5 l d W 1 9 Y 3 8 Z m F r e 2 d 3 r 7 h / 0 F Z x K i l r 0 V j E s u M T x Q S P W E t z L V g n k Y y E v m B 3 / v h q 7 t / d M 6 l 4 H N 3 q S c K 8 k A w j H n B K t J F u m p d u v 1 j C t l v G N b e K s I 0 X M K T s u v j 8 A j m Z U o I M j X 7 x v T e I a R q y S F N B l O o 6 O N H e l E j N q W C z Q i 9 V L C F 0 T I a s a 2 h E Q q a 8 6 e L U G T o x y g A F s T Q V a b R Q v 0 9 M S a j U J P R N Z 0 j 0 S P 3 2 5 u J f X j f V Q c 2 b 8 i h J N Y v o c l G Q C q R j N P 8 b D b h k V I u J I Y R K b m 5 F d E Q k o d q k U z A h f H 2 K / i f t s u 1 U 7 L N m p V R 3 s z j y c A T H c A o O V K E O 1 9 C A F l A Y w g M 8 w b M l r E f r x X p d t u a s b O Y Q f s B 6 + w T g G Y 2 A < / l a t e x i t > Q = 4 < l a t e x i t s h a 1 _ b a s e 6 4 = " Z F S M Y q S g I T Z / J f L f w O C e p S 6 w F G Y = " > A A A B 6 n i c d V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 h q o l 6 E g h e P L d o P a E P Z b D f t 0 s 0 m 7 G 6 E U v o T v H h Q x K u / y J v / x m 0 b Q U U f D D z e m 2 F m X p h y p r T j f F i F l d W 1 9 Y 3 i Z m l r e 2 d 3 r 7 x / 0 F J J J g l t k o Q n s h N i R T k T t K m Z 5 r S T S o r j k N N 2 O L 6 e + + 1 7 K h V L x J 2 e p D S I 8 V C w i B G s j X T b u P L 6 5 Y p j e 7 7 n + u f I s Z 0 F D H F 9 t 3 p W R W 6 u V C B H v V 9 + 7 w 0 S k s V U a M K x U l 3 X S X U w x V I z w u m s 1 M s U T T E Z 4 y H t G i p w T F U w X Z w 6 Q y d G G a A o k a a E R g v 1 + 8 Q U x 0 p N 4 t B 0 x l i P 1 G 9 v L v 7 l d T M d X Q Z T J t J M U 0 G W i 6 K M I 5 2 g + d 9 o w C Q l m k 8 M w U Q y c y s i I y w x 0 S a d k g n h 6 1 P 0 P 2 l V b d e z / Y Z X q b l 5 H E U 4 g m M 4 B R c u o A Y 3 U I c m E B j C A z z B s 8 W t R + v F e l 2 2 F q x 8 5 h B + w H r 7 B N d l j X g = < / l a t e x i t > Q = 6
< l a t e x i t s h a 1 _ b a s e 6 4 = " g p k N 5 6 H Y Y t H u 8 w W H + N U l E L h I 7 z k = " > A A A B 6 n i c d V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B Z B P I S k 9 E M P Q s G L x x Z t L b S h b L a b d u l m E 3 Y 3 Q g n 9 C V 4 8 K O L V X + T N f + O 2 j a C i D w Y e 7 8 0 w M 8 + P O V P a c T 6 s 3 M r q 2 v p G f r O w t b 2 z u 1 f c P + i o K J G E t k n E I 9 n 1 s a K c C d r W T H P a j S X F o c / p n T + 5 m v t 3 9 1 Q q F o l b P Y 2 p F + K R Y A E j W B v p p n V Z G x R L j l 2 v 1 M v l G n J s Z 4 E 5 c a s X T h W 5 m V K C D M 1 B 8 b 0 / j E g S U q E J x 0 r 1 X C f W X o q l Z o T T W a G f K B p j M s E j 2 j N U 4 J A q L 1 2 c O k M n R h m i I J K m h E Y L 9 f t E i k O l p q F v O k O s x + q 3 N x f / 8 n q J D s 6 9 l I k 4 0 V S Q 5 a I g 4 U h H a P 4 3 G j J J i e Z T Q z C R z N y K y B h L T L R J p 2 B C + P o U / U 8 6 Z d u t 2 N V W p d Q 4 y + L I w x E c w y m 4 U I c G X E M T 2 k B g B A / w B M 8 W t x 6 t F + t 1 2 Z q z s p l D + A H r 7 R P l P I 1 8 < / l a t e x i t > Q = 5

10
< l a t e x i t s h a 1 _ b a s e 6 4 = " e P l Y j 3 a y n X f R 9 B 4 M 5 9 B y M 0 O 1 Y j w = " > A A A B 6 X i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m k Y o 8 F L x 6 r 2 A 9 o Q 9 l s J + 3 S z S b s b o Q S + g + 8 e F D E q / / I m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g o b m 1 v b O 8 X d 0 t 7 + w e F R + f i k r e N U M W y x W M S q G 1 C N g k t s G W 4 E d h O F N A o E d o L J 7 d z v P K H S P J a P Z p q g H 9 G R 5 C F n 1 F j p w X M H 5 Y p b d R c g 6 8 T L S Q V y N A f l r / 4 w Z m m E 0 j B B t e 5 5 b m L 8 j C r D m c B Z q Z 9 q T C i b 0 B H 2 L J U 0 Q u 1 n i 0 t n 5 M I q Q x L G y p Y 0 Z K H + n s h o p P U 0 C m x n R M 1 Y r 3 p z 8 T + v l 5 q w 7 m d c J q l B y Z a L w l Q Q E 5 P 5 2 2 T I F T I j p p Z Q p r i 9 l b A x V Z Q Z G 0 7 J h u C t v r x O 2 l d V r 1 a 9 v q 9 V G v U 8 j i K c w T l c g g c 3 0 I A 7 a E I L G I T w D K / w 5 k y c F + f d + V i 2 F p x 8 5 h T + w P n 8 A e g I j O 0 = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " J e q G y x v l Y z C / y / Z I h 9 P T U J E 1 6 F r p r x D A E L / 6 K F w + K e P U r v P k 3 d p a D J j 4 o e L x X R V U 9 N x Z c o W V 9 G 5 m l 5 Z X V t e x 6 b m N z a 3 s n v 7 v X U F E i G d R Z J C L Z c q k C w U O o I 0 c B r V g C D V w B T X d w P f a b 9 y A V j 8 I 7 H M b g B L Q X c p 8 z i l r q 5 g 8 a 3 X K x E w X Q o 5 f W S f u 0 g / C A q X c 1 c r r 5 g l W y J j A X i T 0 j B T J D r Z v / 6 n g R S w I I k Q m q V N u 2 Y n R S K p E z A a N c J 1 E Q U z a g P W h r G t I A l J N O X h i Z x 1 r x T D + S u k I 0 J + r v i Z Q G S g 0 D V 3 c G F P t q 3 h u L / 3 n t B P 1 z J + V h n C C E b L r I T 4 S J k T n O w / S 4 B I Z i q A l l k u t b T d a n k j L U q e V 0 C P b 8 y 4 u k U S 7 Z l d L F b a V Q r c 7 i y J J D c k S K x C Z n p E p u S I 3 U C S O P 5 J m 8 k j f j y X g x 3 o 2 P a W v G m M 3 s k z 8 w P n 8 A r w q W X Q = = < / l a t e x i t > Figure 2. Zero-frequency spectrum of the squeezed quadrature, V2(ω = 0), as a function of the normalized modulation amplitude σ. Different solid lines correspond to different values of the quality factor Q, with the dashed yellow line shows the rotating-wave approximation (Q → ∞) limit. Note that V2 is provided in −dB (i.e. we plot −10 log 10 V2), and hence larger values correspond to better squeezing. Note also that for finite Q the squeezing is maximum far from the instability (further the smaller Q is), which for each value of Q corresponds to the value of σ where the curve halts. We show the optimal squeezing and the corresponding modulation amplitude as functions of Q in Fig. 3.
the spectral covariance matrix is defined as with We remind that we are working with a dimensionless time, and therefore, the detection frequency ω in this equation is also dimensionless, with the real detection frequency given by γω. The spectral covariance matrix (41) is subject, for all ω, to the usual constrains of the standard covariance matrix of Gaussian states [101]. For example, it is real, symmetric, and must posses positive eigenvalues (corresponding to the spectral density of the normal quadratures of the problem), and it satisfies the condition det{V(ω)} ≥ 1 linked to Heisenberg's uncertainty relations. Note that A(ω) has the same form as the generic spectral density that we defined in (23), just replacing the generic correlation function O(t, t ) by the output correlation matrix C out (t, t ) = x out (t)x T out (t ) . Hence, we now proceed to rewrite it in terms of the spectral densities that we have defined in the previous sections. First, note that C out (t, t ) can be written in terms of the previously defined correlations (17) and (30) as where we have used (40), (14), and the two-time correlators of the noises as defined after Eq. (1). Using now the definitions for the spectral densities that we introduced in (24) and (31), the components of A(ω) are rewritten as an expression that is readily evaluated using the results of the previous sections. Specifically, we first solve the Floquet problem (38) numerically, that is, we determine the Floquet exponents {µ α } α=1,2 and K(t) over one period, and then use the simplified expressions of the spectral densities as given in (25) and (31).

C. Squeezing properties
Let us start discussing the results within the rotating-wave approximation. As mentioned above, in this limit the Jacobian in Eq. (39) is time-independent and has the diagonal form L RWA . The particularization of the expressions above to such case easily leads to the following well-known expression for the spectral covariance matrix of Eq. (41):

k t s G W 4 E d h O F N A o E d o L J 7 d z v P K H S P J a P Z p q g H 9 G R 5 C F n 1 F j p w a s N y h W 3 6 i 5 A 1 o m X k w r k a A 7 K X / 1 h z N I I p W G C a t 3 z 3 M T 4 G V W G M 4 G z U j / V m F A 2 o S P s W S p p h N r P F p f O y I V V h i S M l S 1 p y E L 9 P Z H R S O t p F N j O i J q x X v X m 4 n 9 e L z V h 3 c + 4 T F K D k i 0 X h a k g J i b z t 8 m Q K 2 R G T C 2 h T H F 7 K 2 F j q i g z N p y S D c F b f X m d t K + q 3 n W 1 d n 9 d a d T z O I p w B u d w C R 7 c Q A P u o
For σ = 0 this is just the covariance matrix of vacuum for all ω as expected, as in the absence of modulation, the oscillator simply relaxes to its ground state. As σ increases, V 1 (ω) gets larger and larger, while V 2 (ω) gets smaller and smaller, corresponding to quantum squeezing in the momentum quadrature. Eventually, at σ = 1 (the socalled "threshold"), we get V 2 (ω = 0) = 0 and V 1 (ω = 0) = ∞, signaling perfect momentum squeezing, and the breakdown of our ideal linear model. Note that the system remains in a minimum uncertainty state for all σ, since det{V(ω)} = V 1 (ω)V 2 (ω) = 1.
In this work we have studied the deviations of the full V(ω) with respect to this rotating-wave picture. In particular, we summarize our main results through Figs. 1 to 3. Following the notation introduced above within the rotating-wave approximation, let us denote by {V j (ω)} j=1,2 the eigenvalues of the spectral covariance matrix V(ω) with V 1 > V 2 for definiteness. In Fig. 1 we plot these as a function of the dimensionless detection frequency ω, for different values of σ (as indicated in the figure) and Q = 3 (similar behavior is found for any other value of Q ≥ 1; we don't consider extremely-bad oscillators with Q < 1, for which the Born-Markov quantum Langevin description adopted here ceases to make sense). The first thing that we can appreciate from Figs. 1a-c is that even for a finite Q, the optimal squeezing is still found at ω = 0, and is degraded with respect to its rotating-wave value, that is, V 2 (ω) > V RWA 2 (ω). In addition, the spectra show sidebands at ω = ±2nQ, with n ∈ N, as expected for an output field carrying a modulation of period T = π/Q. The sidebands are relatively broad, and have a shape that departs more and more from Lorentzian as σ approaches the instability at which a Floquet eigenvalue becomes zero. We denote such value of σ by σ ins , which we show in Fig. 3b as a function of Q. Remarkably, once very close to the unstable point, the sidebands of V 2 develop a secondary sharper peak that diverges at σ = σ ins (see Fig. 1c). Let us remark that the sidebands do not show squeezing for any value of the parameters; on the contrary, they simply add noise. Moreover, we have also found that the oscillator is not in a minimum uncertainty state anymore, that is, V 1 (ω)V 2 (ω) > 1 for any finite Q. Of course, for any value of the rest of parameters, the product V 1 V 2 approaches 1 as Q increases.
Knowing that maximum squeezing occurs at ω = 0, in Fig. 2 we plot V 2 (ω = 0) as a function of σ for different values of Q. Note that we plot it in −dB units, defined as −10log 10 V 2 , such that higher values correspond to larger squeezing, with 10dB equivalent to 90% of quantum noise reduction or V 2 = 0.1. Contrary to the rotating-wave case, squeezing is not maximized at σ = σ ins , but at an optimal value σ opt that can be rather small for small Q. This is appreciated in Fig. 3b, where we plot σ opt as a function of Q, which of course tends to 1 (the rotating-wave instability) as Q → ∞. Note also that even for moderate values of Q the optimal squeezing is quite large, e.g., ∼10dB at Q = 2, as shown in Fig. 3a, so our theory shows that squeezing in parametric oscillation is quite robust against the quality of the oscillator and the modulation amplitude.

VII. CONCLUSIONS
In this work we have provided an efficient tool for the evaluation of two-time correlation functions and related spectral densities of time-periodic open quantum systems. In particular, using an approach based on the Floquet theorem, we have shown that these quantities can all be related to simple integrals over a single period, which can be efficiently evaluated. Among other applications, this provides a compact and robust tool for the systematic analysis of the corrections that may arise when generating effective dynamics via periodic modulations. In addition, it is a tool that will find applications in the determination of the quantum properties of systems undergoing limit-cycle motion with the corresponding spontaneous breaking of time-translational invariance.
As a testbed for the method, we have studied the quantum properties of a damped parametrically-driven oscillator under conditions where the rotating-wave approximation cannot be invoked. This regime is easily attainable nowadays in low-frequency superconducting or mechanical oscillators that work in the quantum regime. Our results show that even for relatively large modulation amplitudes or low-quality oscillators, large levels of squeezing prevail. However, the optimal squeezing levels occur for a modulation amplitude far below the oscillation instability, which is where rotating-wave results predict optimal squeezing.
Expression (A4) can be rewritten as C αβ (t, t ) = Υ(µ α + µ β )C αβ (t, t ), where which has the same form as Eq. (19) in the main text, except for the fact that the integral (A5b) still requires knowledge of N αβ outside the first period because of its augmented argument. In order to keep the evaluation restricted to the first period, the integral can be worked out as we explain next. First, we perform the variable change t 3 = t 2 + τ , and split the resulting integral as Next we perform the change of variable t 4 = t 3 − T in the second integral, which is the one that extends beyond the first period. Noting that N αβ (t 4 + T ) = N (t 4 ), we obtain N αβ (τ ) = e (µα+µ β )τ which coincides with Eq. (21) in the main text.