Quantum noise properties of multiphoton transitions in driven nonlinear resonators

We investigate the quantum noise properties of a weakly nonlinear Duffing resonator in the deep quantum regime, where only few quanta are excited. This regime is dominated by the appearance of coherent multiphoton resonances in the nonlinear response of the resonator to the modulation. We determine simple expressions for the photon noise spectrum and find that the multiphoton resonances also induces a multiple peak structure in that noise. When the corresponding multiphoton Rabi oscillations are underdamped, zero temperature quantum fluctuations determine comparable populations of all quasienergy states which belong to a resonant multiphoton doublet. Most interestingly, the quantum fluctuations probe the multiphoton transitions by inducing several peaks in the noise spectrum of the resonator observables. In particular, the noise of the photon number contains complete information about the multiphoton states and their stationary populations via pairs of nearly symmetric peaks at opposite frequencies. Their widths are determined by the damping of the Rabi oscillations and their heights are proportional to the stationary nonequilbrium populations. A finite detuning from a multiphoton resonance generates a quasielastic noise peak at zero frequency. In addition, we relate the stationary populations of the quasienergy states with an effective quantum temperature and discuss the role of a finite temperature.


Introduction
Coupling a driven quantum mechanical oscillator to environmental fluctuations allows the oscillator dynamics to reach a stationary state. In the stationary state, energy is coherently absorbed from the pump and leaks into the environment via random dissipative transitions, which inevitably induce noise in the resonator. This occurs even at zero temperature where only environmental zero-point fluctuations (quantum noise) exist. The noise properties of a nonlinear oscillator determine many fundamental nonequilbrium phenomena such as quantum heating [1,2,3] and quantum activation [4].
Nonlinear oscillators are used as basic elements for quantum state detection or amplification. Examples of those are the Josephson bifurcation amplifier [5,6,7,8,9,10] and the cavity bifurcation amplifier [8,11]. In this context, the noise properties of the resonator, which is used as detector or amplifier, determine the backaction of the measurement or amplification on the system itself [12,13,14]. Clearly, it is desirable to keep the backaction as weak as possible, while on the other hand, a significant coupling of the amplification or measurement device to the system is useful in order to achieve a sufficiently strong detection or amplification efficiency. A fundamental lower limit for the introduced disturbance, however, will be set by the quantum noise. Hence, in order to design useful concepts for quantum state detection based on nonlinear resonators in the deep quantum regime, their quantum noise properties have to be addressed.
The Josephson bifurcation amplifier takes advantage of the dynamically induced bistability due to the nonlinearity of the resonator. The eigenstates of the qubit are mapped onto the coexisting stable vibrational states of the resonator, which have different amplitudes and phase relations relative to the phase set by the external drive. Hence, they allow for a large discrimination power. Up to present, these amplifying devices mostly operate in a regime where many quanta in the resonator are excited. This implies that pure quantum fluctuations are typically small on average. Nevertheless, some experiments have been realized at low temperature where the relevant fluctuations are quantum mechanical in nature [8,10]. The regime of weak fluctuations has been the subject of intense theoretical investigation [1,2,3,14,15,16,17,18,19]. It has been shown, that the stationary distribution over the quasienergy states of the driven oscillator at zero temperature has the form of an effective Boltzman distribution, thereby allowing to introduce the concept of an effective quantum temperature implying quantum heating even at T = 0. Signatures of the onset of quantum fluctuations can be seen in the relative intensities of the lines of the resonator noise spectrum [3,15,16,17,18] and in the appeareance of a fine structure in the spectral lines of resonators with comparatively large nonlinearities and large quality factors [3,17]. The spectral fine structure yields detailed information on the quasienergy distribution [3,17].
Recently it has been proposed that nonlinear quantum detectors which operate in the regime of few quanta (deep quantum regime) would bring different advantages, such as a small back action, a large discrimination power with an enhanced readout fidelity, and a sufficiently large measurement efficiency [19]. In the deep quantum regime, the frequency-resolved nonlinear response of the oscillator to the external driving with frequency ω ex shows a rich fine structure [1,2,20,21,22] which is mainly generated by few-photon transitions in the resonator. The splitting of the typical Lorentzian resonance of a harmonic oscillator into a series of non-Lorentzian resonances and antiresonances reflects the intrinsic nonequidistance of the energy levels E n of a nonlinear oscillator. N-photon transitions with the resonance condition E N − E 0 = N ω ex , N = 1, 2, . . ., and the subsequent drift down along the ladder of the few-photon Fock states generate a pronounced nonequilbrium quasienergy distribution which is strongly different from the Boltzman-type [1,2,21]. Peaks or dips in the nonlinear response are a direct consequence of the nonequilbrium distribution over states with different oscillation amplitude and phase [1,2,21]. The signatures of such a characteristic non-Lorentzian lineshape of the response has been observed experimentally in a circuit-cavity QED set-up [23].
In this work, we investigate the noise properties of modulated nonlinear oscillators in the deep quantum regime. We consider the simplest example of a monostable anharmonic oscillator which has a quartic nonlinearity (Duffing oscillator). Such a weakly nonlinear Duffing oscillator has a remarkable symmetry: its energy levels E n with n ≤ N are pairwise resonant for the same driving frequency ω ex , E N −n − E n = (N − 2n) ω ex . An example of the energy spectrum for the case N = 3 is sketched in Fig. 1 (a). After preparing the oscillator in its n-th excited state n ≤ N, it displays periodic quantum oscillations between the n-th and the N − n-th excited states. During these oscillations, |N − 2n| photons are being exchanged between the oscillator and the modulation field. The oscillations of the photon numbern are usually referred to as multiphoton Rabi oscillations. Their characteristic frequency, the Rabi frequency Ω nN , depends on the intensity of the driving field and on the number of photons exchanged. The Rabi frequency Ω 0N for the N-photon oscillations is the smallest Rabi frequency. The multiphoton Rabi oscillations with N −n photons involved are underdamped if their Rabi frequency Ω nN exceeds the dissipative rate of photon leaking into the environment. The latter is the oscillator relaxation rate γ. For γ ≪ Ω 0N all the Rabi oscillations are in general underdamped. The periodically driven resonator reaches its stationary state on the timescale γ −1 .
In the stationary state, quantum noise induces -even at zero temperaturefluctuations in the photon numbern. The dynamics of these fluctuations is characterized by multiphoton oscillations which manifest themselves as peaks in the noise spectrum S(ω) ofn, located at plus/minus the Rabi frequencies Ω nN . In the underdamped regime, the dissipative dynamics of the driven oscillator is most appropriately described in terms of random transitions between the oscillator quasienergy states. When the driving is resonant, the pairs of oscillator Fock states with n-and N − n-photons are resonantly superposed. The corresponding oscillator quasienergy states are a symmetric and an antisymmetric superposition of the two Fock states. Their splitting in quasienergy is given by the Rabi frequency Ω nN . The corresponding peak in the noise spectrum at (−)Ω nN is due to random transitions from the state with (highest) lowest to that with the (lowest) highest quasienergy of the doublet. The peak intensity is proportional to the stationary occupation probability of the initial quasienergy state. Therefore, the noise spectrum offers a convenient way to directly probe the stationary distribution over all the quasienergy states. Moreover, for weak driving and exactly zero detuning from the multiphoton resonance, the noise spectrum of then-photon transition is symmetric, i.e., S(ω) = S(−ω) and two inelastic peaks are signatures of an oscillatory decay of the fluctuations towards the stationary state. States belonging to a multiphoton doublet then have the same stationary occupation probabilities. For a weakly detuned modulation or a stronger driving, the spectrum becomes asymmetric. Besides, an additional quasielastic peak appears at zero frequency which represents incoherent relaxation of the fluctuations towards the stationary state. These features have some analogy in the spectral correlation function of a (static) quantum mechanical two-level system weakly coupled to a dissipative harmonic bath [24]. There, the spin correlation function is a sum of three Lorentzians. The two inelastic peaks are symmetrically located at finite frequencies and their width determines the inverse of the dephasing time. In addition, the quasielastic peak at zero frequency represents incoherent relaxation with the inverse relaxation time given by its width. In the driven system, the appeareance of a quasielastic peak depends on the intriguing interplay between the nonlinearity, the driving strength and the dissipation strength.

Multiphoton Rabi oscillations of the Duffing oscillator
We consider a periodically modulated quantum oscillator with mass m, eigenfrequency ω 0 and a quartic (Kerr) nonlinearity described by the Hamiltonian The modulation amplitude F is assumed to be so small that it induces only weakly nonlinear vibrations. This is guaranteed by the condition αA 2 ≪ mω 2 0 , with A(F ) being the typical amplitude of the nonlinear vibrations. The modulation frequency ω ex is chosen to be close to the oscillator eigenfrequency ω 0 such that the detuning δω is small, i.e., Our theory applies to hard as well as to soft nonlinearities α ≶ 0, but for concreteness we will focus on the case of a hard nonlinearity, α > 0. The quantum dynamics of the weakly detuned and weakly nonlinear driven oscillator is most conveniently described in terms of the oscillator ladder operators a and a † , in a rotating frame determined by the unitary transformation In the rotating frame, the typical time scale of the resonator dynamics is given by δω −1 , so that terms oscillating with frequencies ±2ω ex and ±4ω ex average out and can be We emphasize that relaxational transitions at zero temperature typically occur in both directions, i.e., downwards and upwards along the quasienergy surface, which is in striking contrast to dissipative transitions in static potential surfaces, where only "downward relaxation" is possible. An escape due to "upward relaxation" is known as quantum activation [4].
neglected in the transformed Hamiltonian R(t) H(t) R † (t) − i R(t)Ṙ † (t). Thereby, we obtain the RWA Hamiltoniañ wheren ≡ a † a is the photon number operator, ν and f are the frequencies associated with the Kerr nonlinearity and the external field amplitude at the quantum scale x ZPF = /mω 0 , i.e., ν = 3αx 4 ZPF /4 and f = F x ZPF / √ 2 . In order to keep the notation compact we have set = 1 in Eq. (4) and in the remainder of the paper. The oscillator quasienergies ε n and quasienergy states |ψ n are the eigenvalues and eigenvectors of the rotating wave Hamiltonian,H|ψ n = ε n |ψ n . For vanishing driving, the quasienergy spectrum is given by ε n = δω n + νn(n + 1)/2 for f → 0 .
We are primarily interested in studying the noise spectrum in presence of multiple multiphoton resonances E N −n − E n = (N − 2n)ω ex for n ≤ N, or equivalently ε N −n = ε n for f → 0. From Eq. (5) we find the resonant condition Up to leading order in the driving, the quasienergy eigenstates |ψ n for n ≤ N = N/2 are a resonant superposition of the pair of oscillator Fock states |n and |N − n , i.e., We choose the signs − and + for n < N/2 and N/2 < n ≤ N, respectively. In the following, we refer to the resonant superposition of a pair of Fock states as resonant states or multiphoton states. The states |ψ n which are not involved in a multiphoton transition (n > N and n = N/2 for N even) can be approximated as the corresponding Fock states |ψ n ≈ |n . The Rabi frequency Ω nN −n of the multiphoton oscillations within the pair of Fock states |n and |N − n is given by the splitting of the corresponding levels ε n and ε N −n [25,21] as The resonance condition in Eq. (6) is not renormalized by a finite driving within the RWA. Only for a comparatively larger driving f ∼ ν ≪ ω 0 , the multiphoton transitions have to be reinterpreted as tunneling transitions between semiclassical states [26,27].
As we shall detail in Section 4, the multiphoton Rabi oscillations induce peaks in the spectral densities of oscillator observables only when the Rabi frequency Ω 0N for the multiphoton transition from the zero-photon ground state is larger than the noise-indcued level broadening of the relevant quasienergy levels ε 0 and ε N . In the next section, we will pave the wave for the calculation of the noise spectrum in this regime, by formulating the master equation for a weakly nonlinear oscillator and by evaluating the stationary occupation populations over the quasienergy states.

Stationary dynamics in the deep quantum regime
In the presence of a weak bilinear coupling to the fluctuations of a bosonic bath, the assumptions of small detuning and weak nonlinearity that underly the RWA naturally lead to a Liouville-von Neuman quantum master equation in the Lindblad form for the density matrixρ of the weakly damped oscillator in the rotating frame [21,22], Here, L and D[O] are the Liouville and the Lindblad superoperators, respectively. Moreover, γ is the oscillator damping rate for which we assume that γ ≪ ω 0 . It results from a standard Ohmic bath spectral density J(ω) = γω. In addition,n is the bosonic occupation number at frequency ω 0 and temperature T and is given bȳ

The stationary distribution
For long times, the density matrix in the rotating frameρ relaxes to a stationary statê ρ ∞ , satisfying When the oscillator decay rate γ is larger than the driving, γ ≫ f , the width of the resonant quasienergy levels ε n induced by the bath fluctuations are larger than the corresponding Rabi frequency Ω nN −n of the multiphoton transitions. Then, the multiphoton resonances are smeared out and the coherent effects associated with multiphoton oscillations are strongly suppressed. Hence, dissipation sets a lower limit for the driving strength, f ≫ γ, which has to be overcome in order to observe multiphoton transitions. On the other hand, for comparatively larger driving f ∼ ν, quantum fluctuations are significantly reduced and the oscillator is latched to a classical attractor at asymptotic times which are much larger than the typical relaxational time scale, γ −1 .
Here, we restrict our analysis to the deep quantum regime where the driving is larger than the damping but smaller than the nonlinearity, γ ≪ f ≪ ν. Thereby, we have implicitly assumed a comparatively large nonlinearity ν ≫ γ.

Underdamped regime:
We start our discussion assuming that all Rabi oscillations are underdamped. Put differently, we assume that the smallest Rabi frequency Ω 0N is larger than the relevant level broadening. We refer to this regime as the fully underdamped regime. Then, the off-diagonal matrix elements ofρ ∞ projected onto the quasienergy basis |ψ l are negligible and we can set them to zero, i.e., we perform a secular approximation, Then, a balance equation for the stationary occupation probabilities ρ ∞ ll follows from Eqs. (9) and (10) according to Here, W l,k is the transition rate from state |ψ k to state |ψ l , and γ l is the width of quasienergy level ε l given by γ l ≡ k =l W l,k . We can now formulate more precisely the condition for underdamped Rabi oscillations to occur within the pair forming the narrowest resonance, which is Ω 0N ≫ γ 0 . The solution for stationary occupation probabilities up to leading order in the small parameters f /ν andn is given in Ref. [21]: The pair of multiphoton states |ψ n and |ψ N −n in Eq. (7) have equal stationary population, i.e., ρ ∞ nn = ρ ∞ N −n N −n . The pair with the narrowest resonance has the occupation probabilities ρ ∞ 00 = ρ ∞ N N . The occupation probability grows algebraically with n < N/2 as The states |ψ l with l > N have vanishing occupation probability, ρ ∞ ll = 0. As follows from the discussion above, the degeneracy ρ ∞ 00 = ρ ∞ N N is approximate and is lifted for higher order in f /ν.

Quasienergy distribution close to a multiphoton resonance:
One can easily generalize the above expressions to the case where the detuning δω does not exactly match the resonant condition, δω = δω N . Since the Rabi frequencies for the different pairs of resonant transitions in Eq. (8) are exponentially different, we can choose |δω−δω N | ≪ Ω 1N −1 , so that all the pairs of Fock states |n and |N −n with 1 < n < N/2 are still resonantly superposed, except for . The corresponding solution for the stationary density matrix close to resonance is [21] ρ

3.1.3.
Partially underdamped regime: Next we consider a comparatively large relaxation rate γ, so that the narrowest Rabi resonance is overdamped but the remaining resonances are still underdamped, Ω 0N ≪ Nγ ≪ Ω 1N −1 . We refer to this regime as the partially underdamped regime. Then, incoherent multiphoton transitions from the ground state |0 to state |N with a small rate Ω 2 0N /(Nγ) and the subsequent emission of excitations into the bath determines a small but finite occupation of the resonant states ρ ∞ nn , n ≥ 1. Formally, the stationary distributionρ ∞ can be obtained by setting all the off-diagonal elements of ρ ∞ lk to zero except for ρ ∞ N 0 and ρ ∞ 0N and solving Eq. (10). Thereby, we find The crossover between this solution and the fully underdamped solution Eq. (14) is given in Ref. [21]. Both stationary nonequilbrium distributions are determined by quantum fluctuations and are very different from the equilibrium Boltzmann-type distribution when a driven resonator is latched to a classical attractor.

The nonlinear response of the oscillator
In the steady state regime, t ≫ γ −1 , the oscillator state is described by the timeindependent density matrixρ ∞ in the rotating frame and the oscillator dynamics is embedded in the time-dependent reference frame R(t). The mean value of an observables O is Therefore, the stationary oscillations of the position expectation value x(t) ∞ are sinusoidal, It has been shown that the nonlinear response x(t) ∞ of the oscillator as a function of ω ex shows resonances and antiresonances in the deep quantum regime [20,21,22]. The response is proportional to the transmitted amplitude in a heterodyne measurement scheme and it has already been measured for a weakly nonlinear oscillator [23]. Clearly, such a measurement scheme, or more general, any measurment scheme which probes stationary mean values as opposed to correlations does not allow to resolve the different degenerate resonances separately. Neither, they allow us to access the stationary distribution ρ ∞ ll directly. This becomes possible only when correlations, e.g., via noise spectra are measured. In the next section, we show that this can indeed be achieved by measuring the spectrum of the photon number noise.

Definition of the noise spectrum
The Lindblad master equation (9) in general also allows to investigate transient phenomena and correlation functions. Its formal solution for a given initial stateρ 0 is given byρ(t) = e Ltρ 0 . Moreover, a general correlator O ′ (t ′ )O(t) can be evaluated as the mean value of the operator O ′ at time t ′ with the virtual operator R † (t)OR(t)ρ(t) at time t. This view has been established several decades ago by the Lax formula [28,29] according to For long times t ≫ γ −1 , we find that In general, such correlators are periodic functions of the preparation time t. The noise spectrum is defined as a double average over quantum fluctuations and the time t.
Since this correlator does not depend on the initial time t as a consequence of the RWA, we can define the noise spectrum in terms of a single average over quantum fluctuations according to It is useful to separate the contributions to S(ω) into those coming from the expectation value ofn, and those from its fluctuations, i.e., Here, δn is the operator for the photon number fluctuations, i.e., δn =n − n ∞ . Our path to compute the noise spectrum consists in three steps: i) We express the virtual preparationnρ ∞ in terms of right eigenvectors of the superoperator L. ii) We plug the resulting decomposition into Eq. (22). Then, each term decays exponentially with a different exponent which is given by the corresponding eigenalue of L. iii) We compute the Fourier integral in Eq. (23), which thereby yields a sum over (overlapping) Lorentzians.
The general expression, which is useful for a concrete numerical evaluation, for the noise spectrum given in terms of the eigenvectors and the eigenvalues of L is derived in Appendix A. In the next section, we consider the special case of underdamped multiphoton Rabi oscillations.

Noise spectrum in the underdamped regime
When all the multiphoton Rabi oscillations are underdamped, Ω N 0 ≫ Γ N , the coherences |ψ N −n ψ n | and |ψ n ψ N −n | are approximate eigenvectors of the Liouvillian L. Then, with the level widths being given as Γ n = γ n = γ(n + 1/2)N + γn for n < (N − 1)/2. For N odd, Γ (N −1)/2 = γ(1 + 2n)(5N + 1)/8 + γn. Up to leading order in f /ν, the decomposition of the virtual preparationnρ ∞ in terms of right eigenvectors of L has the simple expression Clearly, each term of the above decomposition yields a Lorentzian peak in the noise spectrum S(ω). The first term yields the contribution to S(ω) from the expectation value ofn, (N/2) 2 δ(ω). The remaining terms yield inelastic peaks associated to random transitions between quasienergy states belonging to the same multiphoton doublet. Since the populations ρ ∞ nn and ρ ∞ N −nN −n are approximately equal, peaks at opposite frequency have approximately equal intensity. By putting together Eqs. (22), (23), (25), and (26), we find S(ω) = (N/2) 2 δ(ω) + δS(ω) with Hence, the Lorentzians are centered at the multiphoton Rabi frequencies Ω nN −n and have a resonance width of Γ n . The factor (N − 2n) 2 /4 is the leading order expression for the squared matrix element | ψ n |n|ψ N −n | 2 . Remarkably, the line intensities depend only weakly on the driving f and on the temperature through the stationary distribution ρ ∞ nn . Up to leading order, the driving f enters only in the splitting of the lines through the Rabi frequencies. Notice that Eq. (27) is valid only in the vicinity of a multiphoton peak since terms of order γ are not taken into account. In order to evaluate the tails of the peaks more precisely, one has to take into account the contribution stemming from all eigenvectors of L, see Appendix A.
In the left and right panels of Fig. 2, we show the noise spectrum S(ω) for the cases N = 2 and N = 3, respectively. The noise spectrum for N = 2 shows a pair of symmetric peaks which correspond to the transitions |ψ 0 ↔ |ψ 2 . Likewise, the noise spectrum for N = 3 displays two pairs of symmetric peaks corresponding to the transitions |ψ 0 ↔ |ψ 3 and |ψ 1 ↔ |ψ 2 . The green dashed lines mark the results from our approximate analytical formula in Eq. (27) while the yellow solid lines show the data obtained by numerically evaluating the expression in Eq. (A.2). An excellent agreement is found.
In Fig. 2a), additional smaller side peaks of the order of f /ν are also visible, see the gray lines representing a ten-fold zoom. They are not associated to any resonant transition between multiphoton states and are thus not captured by the leading order expression given in Eq. (27). The particular subleading peaks in Fig. 2a) belong to the transitions |ψ 0 ↔ |1 .
These features have a direct analogy in the spectral correlation function of a static quantum mechanical two-level system which is weakly coupled to a dissipative harmonic bath [24]. For a general biased two-state system with anticrossing energy levels, the pair correlation function is a sum of three Lorentzians. The two inelastic peaks are symmetrically located at finite frequencies and their width determines the inverse of the dephasing time. For a biased static two-level system away from resonance, an additional quasielastic peak at zero frequency appears which represents incoherent relaxation with the inverse relaxation time given by its width. Since we consider here the case strictly at resonance (in the RWA), no zero-frequency peak is present.

Photon anti-bunching
In general, the photon emission characteristics of a quantum mechanical resonator can show peculiar nonclassical features. For instance, counterintuitive correlation phenomena such has photon antibunching can occur, where the photon number correlation function for short delay times is smaller than the one for classical, uncorrelated photons. This implies that the probability for photons to arrive in pairs is suppressed [30]. Our approach provides a natural framework to investigate a possible non-Poissonian statistics of the multiphoton events in the nonlinear resonator. Therefore, we consider the normalized photon number correlation function or secondorder coherence function defined as For long delay times τ , the counts of two photons with a delay time τ are statistically independent events, g (2) (τ → ∞) = 1. For vanishing delay times, we have Photon antibunching corresponds to the case g (2) (τ = 0) < 1. For the fully underdamped case, we find the expression which represents the known result of the second-order correlation function of the electromagnetic field [30]. Hence, the oscillator displays photon antibunching close to a multiphoton transition. The second-order coherence of the stationary state of the quantum Duffing oscillator at the N-th multiphoton resonance has the same value as the second order coherence for an oscillator prepared in the single Fock state |N , in spite of its fluctuations over the quasienergy states.

Lineshape of the noise spectrum close to a multiphoton resonance
In presence of a small detuning from the multiphoton resonance, δω − δω N ∼ Ω 0N , the states |ψ 0 and |ψ N are no longer a resonant superposition of the Fock states |0 and |N . Hence, the corresponding stationary occupation probabilities ρ 00 and ρ N N , given in Eq. (16), become significantly different. In turn, the pair of peaks S 0 (ω) and S N (ω), which are associated to the transitions |ψ 0 ↔ |ψ N , become asymmetric such that S 0 (ω) = S N (−ω). This behavior is shown in shown in Fig. 4a) for the case around the 3-photon resonance. The peak lineshapes can readily been evaluated and we find Their distance increases with the quasienergy splitting, ε N − ε 0 = sgn(δω − δω N )(Ω 2 0N + N 2 |δω − δω N | 2 ) 1/2 , whereas the peak width does not change close to the multiphoton resonance, δω −δω N ∼ Ω 0N . The asymmetry is determined by the stationary occupation probabilities ρ ∞ 00 and ρ ∞ N N . From Eq. (16), we find (33) The above expression is valid for ω close to the center of the largest peak, ω ∼ ε N − ε 0 , and |δω − δω N | not too large such that S(±ω) ≫ γ.
In addition to the peaks at finite frequencies (which induce decaying coherent multiphoton Rabi oscillations), also a zero frequency peak appears. This quasielastic peak is associated to incoherent relaxational decay of the multiphoton Rabi oscillations and is also known for the noise correlation function of a static biased quantum two level system [24].
In Fig. 3b), we show the logarithm of the asymmetry ratio given in Eq. (33). The asymmetry shows a clear maximum at approximately ε 3 − ε 0 .
To further illustrate the asymmetry in the peak heights, we show in Fig. 3c) the peak maxima associated to the transitions |ψ 0 → |ψ 3 and |ψ 3 → |ψ 0 . At the 3-photon resonance (black dashed vertical line), both peaks are equal in height (symmetric noise spectrum). Away from the resonance, the low (high) frequency branch aquires more spectral weight for negative (positive) detuning.

Photon noise at zero frequency
Fluctuations of an oscillator (quasi)energy induce a broad (with width ∝ γ) zero frequency peak in the noise spectrum of an observable whose mean value depends on the (quasi)energy [31]. For weak driving f ≪ ν and at a resonance |δω − δω N | ≪ Ω 0N , the quasienergy states of the Duffing oscillator have large fluctuations as several quasienergy Figure 3. (a) Asymmetric structure of the photon noise spectrum at frequency δω = δω 3 + δ, i.e., out of resonance for a detuning δ = 1.6 × 10 −4 ν for the same parameters used in Fig. 2b) (orange solid line). In addition, we show in the background the symmetric photon noise at the resonant frequency δω 3 (grey shadowed area). Moreover, we depict the inverted case δ → −δ, which shows a symmetric behavior under the reflection ω → −ω (green solid line). (b) Noise asymmetry via the logarithm of Eq. (33) for the same parameters as in (a). (c) Height of the photon noise peak for the transition |ψ 0 → |ψ 3 (orange solid line), and |ψ 3 → |ψ 0 (green solid line) as a function of the external frequency. The peak maximum is located at δω 3 ± δ.
states have comparable occupation probabilities even at T = 0. However, the mean value ofn becomes independent from the quasienergy, ψ n |n|ψ n ≈ N/2 for n ≤ N. As a consequence, the contribution to the noise spectrum ofn coming from fluctuations δS(ω) does not have a peak at zero frequency since δS(0) ∝ γ. Close to resonance, when |δω − δω N | ∼ Ω 0N , two dynamical effects compete: on one hand, the quasienergy fluctuations quickly decrease for increasing detuning, i.e., moving away from resonance as the occupation probability of the state |ψ 0 approaches one. On the other hand, the mean value ofn becomes strongly dependent on the quasienergy. As a result of this competition, the intensity of the zero frequency noise plotted as a function of δω has two maxima at the two opposite sides of the resonant value δω N . In Fig. 4, we show the zero frequency noise for the special case N = 2. The yellow solid line represents the intensity at zero frequency computed numerically, while the green dashed line is the leading order contribution (in f /ν)

Noise spectrum towards the semiclassical regime
Next, we investigate the noise spectrum for larger driving strengths, f ν. In order to illustrate how the noise spectrum changes for increasing driving, we show the intensities of the brightest peaks as a function of the driving strength for the N = 5-photon resonance in ; see Fig. 5a). In Fig. 5b), we also show the quasienergy spectrum, and the  noise spectrum for a comparatively large value of the driving amplitude f = ν is shown in Fig. 5c). A peak in the noise spectrum at frequency ω = ε l − ε k is associated to a single transition |ψ k → |ψ l and is given by Hence, the relative intensities of a pair of peaks at opposite frequencies is still related to the occupation probability of the corresponding initial states through S(ε l − ε k )/S(ε k − ε l ) = ρ kk /ρ ll . For weak driving, we have three pairs of approximately symmetric peaks as described by Eq. (27). Each peak corresponds to a transition between two states belonging to a multiphoton doublet of quasidegenerate states: |ψ 0 ↔ |ψ 5 , |ψ 1 ↔ |ψ 4 , and |ψ 2 ↔ |ψ 3 . For increasing driving, the spectrum becomes increasingly asymmetric. For moderate values of the driving, the noise spectrum undergoes two major qualitative changes: i) the peak at zero frequency becomes clearly visible; ii) a pair of peaks corresponding to the transitions |ψ 1 ↔ |ψ 3 acquires a significant intensity. For f = ν, the peak associated with the transition |ψ 3 → |ψ 1 is even the second brightest peak.
These qualitative changes can be explained in terms of a semiclassical description valid beyond the weak driving limit. The RWA Hamiltonian in Eq. (4) can be rewritten in terms of rotating quadratures, and interpreted as a quasienergy surface in phase space [26,27]. It has the shape of a tilted mexican hat and is sketched in Fig. 1 (c) for two values of f . The larger f is, the stronger is the induced tilt. The local maximum and the minimum of the quasienergy surface are the classical attractors. In the static frame, they describe stationary oscillations with a small and a large amplitude, respectively. In the vicinity of the attractors the vibration amplitude and the slow part of the oscillation phase display slow vibrations with frequency ∝ δω. In absence of resonant transitions, each quasienergy state can be associated to a quantized quasiclassical orbit which lies on the internal surface around the local maximum, on the external surface, or along the quasienergy well around the minimum. For very weak driving, f ≪ ν/ 2(N + 1), the quantum mechanical Fock states |n with n < N/2 are associated to quasiclassical trajectories on the internal surface around the local maximum, whereas the Fock states with photon number n larger than N/2 are associated to semiclassical orbits on the external surface. Within this representation, the multiphoton transitions can then be reinterpreted as tunneling transitions between the internal and the external parts of the surface [26,27]. For comparatively larger driving, the zero-point quasienergy associated to the slow vibrations around the minimum (∝ δω) becomes smaller than the dynamical barrier height. Then, quasienergy states appear which are localized in the quasienergy well. In turn, the noise spectrum becomes qualitatively different from the one for weak driving. The small quantum fluctuations around the minimum of the quasienergy surface can be described in terms of an effective auxiliary oscillator with ladder operators b and b † and are given by Here, a h is the amplitude of the stationary oscillations rescaled by √ 2x ZPF [17,18]. They can be mimicked by a local effective quantum temperature T e = (2k B ln coth r * h ) −1 which depends on the squeezing factor r * h [1,2,17,18]. For f = ν, the states |ψ 2 , |ψ 3 , and |ψ 1 can be identified with the groundstate and first two excited states of the auxiliary oscillator (but in the remainder of this discussion we keep the same labels for the states as in the weak driving limit). The level spacing ε 3 − ε 2 is of the order of the frequency of the slow classical oscillations of the amplitude and slow part of the phase.
Such oscillations appear in the noise spectral density of a classical oscillator as a pair of peaks. In a nonlinear quantum oscillator whose quasienergy levels are not equidistant and their distance exceeds the damping strength, the classical peaks have a "quantum" fine structure [3]. In the present case of the Duffing oscillator, the classical noise peak is splitted into two peaks associated to the nearest neighbor transitions between the ground state and the first excited state, and the first and the second excited state, |ψ 2 ↔ |ψ 3 and |ψ 3 ↔ |ψ 1 , respectively. Their peak height is proportional to the square of the rescaled vibration amplitude a h and to the occupation of the initial state ρ ∞ nn . The latter, in particular, is governed by the quantum temperature T e . For the ratio of the peak heights, we find [3] Next nearest neighbor transitions can also yield peaks in the noise spectra of a Duffing oscillator [18]. In the present case, the transitions |ψ 2 ↔ |ψ 1 yield a pair of dimmer peaks, however, located at frequencies outside the frequency range shown in Fig. 5.
In the weak damping, weak driving regime discussed so far, the quasienergy well around the minimum is still very shallow, and the oscillator can escape from the small amplitude attractor via tunneling. Therefore, the oscillator is not latched to any of the attractors and the noise spectral density has also peaks which are associated to intrawell transitions. In particular, the pair of peaks with the smallest splitting describes coherent tunneling oscillations between the internal and the external part of the quasienergy surface (coherent dynamical tunneling or multiphoton Rabi oscillations).
Before closing this section, we mention that for the stronger driving f = ν, also a zero frequency peak appears in the noise spectrum, see Fig. 5c), although the frequency detuning has been fixed to the 5-photon resonance δω = δω 5 . However, as discussed above, this resonance condition is only valid for small f ≪ ν, which is obviously not fullfilled. So the larger driving induces an effective small detuning away from the exact avoided quasienergy level crossing and generates an effective bias. Then, a relaxation pole appears in the relevant self energy [24] which corresponds to a quasielastic relaxation peak at zero frequency.

Dependence of the noise spectrum on damping and temperature
So far, we have analyzed the case of zero temperature and small damping,n ≪ 1 and γ ≪ Ω 0N . In this section, we briefly address how the noise spectrum is modified for larger damping and finite temperature by presenting numerical results of the spectrum in a broad parameter range.
In Fig. 6a), we show S(ω) for different values of the damping for the 3−photon resonance where δω = δω 3 . As expected, the peaks in the noise spectrum get broader for increasing damping. Outside the fully underdamped regime, the two peaks of the pair associated with the transitions |ψ 0 ↔ |ψ 3 start to overlap and eventually merge into a single peak at zero frequency. Thereby, the zero frequency noise is no longer suppressed S(ω ≈ 0) ∝ γ −1 , since incoherent relaxation prevails over coherent decay for large damping. The peaks associated with the underdamped transitions |ψ 1 ↔ |ψ 2 are still described by Eq. (27), even when the spectrum has a peak at zero frequency. The decrease in the peak intensities reflects the decrease of the populations ρ ∞ 11 and ρ ∞

22
in the partially underdamped regime. The dependence of the noise spectrum on temperature is shown in Fig. 6b) and behaves qualitatively similarly. For small temperaturesn ≪ 1, the spectrum is described by Eq. (27). The temperature dependence enters in the line widths of the quasienergy levels as well as in the stationary distribution ρ ∞ nn . For larger temperatures, the two lowfrequency peaks merge into a single peak at zero frequency and the side peaks becomes increasing broader as expected.

Conclusions
In recent years, the rich phenomenology of driven and damped nonlinear quantum oscillators has been impressively consolidated, including their nonlinear response behavior in form of resonant and antiresonant amplification, quantum coherent multiphoton Rabi oscillations, quantum activation and quantum heating. Gradually, the nontrivial effects visible in noise correlation functions have also moved to the focus of interest. Those become relevant whenever a nonlinear quantum oscillator is used as a central element in an amplifier or quantum measurement device. In this work, we have analyzed the noise properties of the quantum Duffing oscillator in the regine when only few quanta are excited. Then, the nonlinear response shows pronounced multiphoton peaks which are associated to resonant multiphoton Rabi oscillations. The noise properties of these multiphoton transitions show a rich phenomenology. To obtain the noise spectrum by analytical means, we invoke the Lax formula for the autocorrelation function of the photon number at different times and calculate its Fourier transform. Exactly at a multiphoton resonance, the noise spectrum consists in a collection of pairs of related resonances which are located at opposite frequencies and which are equal in height. Each pair is associated to a multiphoton doublet. In spite of large fluctuations over the oscillator quasienergy, no quasielastic peak occurs at zero frequency. This is a consequence of a special symmetry of the quantum Duffing oscillator: all quasienergy states which are associated to a multiphoton doublet have the same mean value of the photon numbern.
Slightly away from a multiphoton resonance, the noise spectrum becomes asymmetric and the two resonances are no longer equal in height. In addition, as the mean values ofn become different for quasienergy states with comparable occupations, the quasielastic peak emerges. Since the quasienergy fluctuations are suppressed away from a multiphoton resonance, the intensity of the quasielastic peak as a function of the detuning displays a maximum at the two opposite sides of the resonant value δω N .
Our approach also allows us to evaluate the transition to the semiclassical regime by increasing the photon number by a larger driving amplitude. Then, a quasiclassical quasipotential landscape in phase space is a convenient tool to understand the stationary nonequilibrium dynamics. This view directly leads to quantum mechanical squeezed states which exist close to the local minimum of the quasienergy landscape. A harmonic expansion allows us to characterize the quantum fluctuations via an effective quantum temperature. At larger (real) temperature and damping strengths, all these quantum coherent features are washed out.
Although the time-resolved detection of noise properties of quantum observables of driven resonators requires considerably more experimental effort, we are confident that future experiments will soon elucidate the importance of quantum noise in these systems.