Non-equilibrium steady state of a driven levitated particle with feedback cooling

Laser trapped nanoparticles have been recently used as model systems to study fundamental relations holding far from equilibrium. Here we study, both experimentally and theoretically, a nanoscale silica sphere levitated by a laser in a low density gas. The center of mass motion of the particle is subjected, at the same time, to feedback cooling and a parametric modulation driving the system into a non-equilibrium steady state. Based on the Langevin equation of motion of the particle, we derive an analytical expression for the energy distribution of this steady state showing that the average and variance of the energy distribution can be controlled separately by appropriate choice of the friction, cooling and modulation parameters. Energy distributions determined in computer simulations and measured in a laboratory experiment agree well with the analytical predictions. We analyse the particle motion also in terms of the quadratures and find thermal squeezing depending on the degree of detuning.

Abstract. Laser trapped nanoparticles have been recently used as model systems to study fundamental relations holding far from equilibrium. Here we study, both experimentally and theoretically, a nanoscale silica sphere levitated by a laser in a low density gas. The center of mass motion of the particle is subjected, at the same time, to feedback cooling and a parametric modulation driving the system into a nonequilibrium steady state. Based on the Langevin equation of motion of the particle, we derive an analytical expression for the energy distribution of this steady state showing that the average and variance of the energy distribution can be controlled separately by appropriate choice of the friction, cooling and modulation parameters. Energy distributions determined in computer simulations and measured in a laboratory experiment agree well with the analytical predictions. We analyse the particle motion also in terms of the quadratures and find thermal squeezing depending on the degree of detuning.

Introduction
In a macroscopic system, thermodynamic quantities such as the work carried out during a thermodynamic transformation or the heat exchanged with a heat bath have well defined values due to the statistics of large numbers. For instance, if we repeatedly carry out a certain thermodynamic transformation always starting from the same initial state and following the same protocol, the work performed on the system will always be the same. In small systems, on the other hand, thermodynamic quantities typically fluctuate. Then the work and heat of a thermodynamic transformation, carried out, for instance, by stretching a single biomolecule in solution, need to be characterised with a statistical distribution rather than a single value. Even small systems, however, are subject to the basic laws of thermodynamics and, on the average, obey the second law usually formulated in terms of inequalities. As realised by Jarzynski, more specific results can be derived for the fluctuations of work and other quantities that transform the inequalities of thermodynamics into equalities [1,2,3,4], which remain valid arbitrarily far from equilibrium. Such so-called fluctuation theorems have now been derived for several quantities, such as heat, work and entropy [5,6], shedding new light on the significance of irreversibility and the second law at the nanoscale [7,8]. Besides their fundamental importance, fluctuation theorems also provide the basis for the interpretation of single-molecule experiments [9,10,11] as well as for the development of novel non-equilibrium computer simulation methods [12].
Experimentally, fluctuation relations have been studied in a variety of systems mainly in the over-damped regime, such as a particle dragged through a liquid [13] or a biomolecule in solution [10], where the system is strongly coupled to a thermalising environment. Recently, several experimental setups for the investigation of nonequilibrium fluctuations under low-damping conditions were proposed [14,15,16,17]. Due to their weak coupling to the heat bath, such systems hold the promise to enable investigation of the statistics of non-equilibrium fluctuations in the quantum regime. Also, the precise control over the dynamics that can be achieved in such systems permits to construct situations in which microscopic reversibility does not hold.
Here, we study, using theory, simulation and experiment a levitated nanoparticle in the low-friction regime [18]. In particular, we derive analytical expressions for the energy and phase-space distribution of the system in non-equilibrium steady states. Based on these distributions one can relate heat, entropy and energy to each other, thereby providing additional insight into the physics underlying the fluctuation theorems. The particle, consisting of a dielectric material, oscillates in a laser trap and is surrounded by a low-density gas, which exerts frictional and random thermal forces on the particle. The amount of friction can be controlled by changing the pressure of the gas. In addition, the particle is subjected to a nonlinear feedback cooling mechanism and a parametric modulation. Together, these effects allow to bring the oscillating particle into a variety of non-equilibrium steady states with tuneable parameters, turning such nano-mechanical oscillators into ideal test-systems for studies of stochastic thermodynamics. Based on a Langevin equation written for the oscillating particle, we derive analytical expressions for the energy distribution in the stationary states and find that, under appropriate circumstances, our theoretical predictions agree very well with the energy distributions observed in the simulations. In addition, we find that in our experiments parameter fluctuations dominate the noise contribution from Brownian motion, which leads to and additional broadening of the experimental distributions.
In addition to the levitated nanoparticle considered here, our model applies to other nonlinear oscillators, including ultra high-Q nano-mechanical oscillators fabricated from silicon nitride [19] and carbon nanotubes and graphene resonators [20]. The latter naturally exhibit nonlinear damping that is formally identical to our feedback mechanism. Thus, in addition to providing insights into thermodynamics on the nanoscale, the work presented here provides insight into the interaction of noise with inherent nonlinearities of nano-mechanical oscillators and the resulting amplitude and phase noise. Most notably phase noise, despite being an active topic of research for many decades, is still a pertinent topic today [21,22,23], since it plays a prominent role for the application of such systems as sensors and in timing and frequency control.
The remainder of the article is organised as follows. In Sec. 2 we lay out the theory for the energy distribution of a nano-mechanical oscillator subject to friction, nonlinear feedback cooling and parametric modulation. Computer simulations are then used, in Sec. 3, to verify the theoretical predictions and probe the limits of the theory. In Sec. 4 we first describe our experimental setup and explain how we determine the relevant system parameters. We then present energy and phase distributions and discuss how they compare with theory and simulations. Some conclusions and an outlook are provided in Sec. 5.

Equation of motion
We consider a particle of mass m oscillating in a trap with a Duffing potential where q specifies the position of the particle, k is the trap stiffness, and ξ is the Duffing parameter, which quantifies how strongly the trap deviates from a purely harmonic potential. Using the frequency Ω 0 = k/m of the harmonic case, the total energy of the oscillator is given by where p = mq is the momentum of the particle. The force due to the trap is hence given by Since the particle is immersed in a low density gas of temperature T , it experiences also a frictional force F friction = −Γp (4) and the related fluctuating random force where Γ is the friction constant, k B is the Boltzmann constant and w(t) is white noise. A feedback of strength η, acting on the particle with force used to control the effective temperature of the center of mass motion of the particle and cool it far below the gas temperature T [18]. In addition, the particle is driven parametrically by periodically modulating the trap stiffness with frequency Ω m leading to the force where the modulation depth ζ determines the intensity of the parametric driving. Taken together, these forces yield the following stochastic equations of motion for the motion of the particle in the trap, Here, W (t) is the Wiener process with Note that W 2 (t) = t for any time t ≥ 0 and, thus, for an infinitesimal time interval dt one has (dW) 2 = dt. The white noise w(t) appearing in the random force can be viewed as the time derivative of the Wiener process, w(t) = dW (t)/dt. In order to determine the energy distribution of the oscillator in the steady state, we now examine the time evolution of the energy generated by the stochastic equations of motion. To avoid multiplicative noise, i.e., a noise term with an amplitude depending on the current value of the energy, we consider the square root of the energy rather than the energy itself, Applying Ito's formula [24] for the change of variables to ǫ(q, p) we find that the change dǫ during a short time interval is given by where is the sum of the non-conservative forces consisting of the frictional force F friction , the feedback force F feedback and the driving force F drive . Note that the conservative forces, including the force due to the non-linear Duffing term in the energy, do not contribute to the energy change.
The stochastic equation of motion for ǫ, Equ. (13), explicitly depends on the position and momentum of the particle. To eliminate this dependence and obtain a closed equation depending only on ǫ, we observe that the particle settles into a periodic motion with a frequency Ω that is not necessarily equal to the frequency Ω 0 of the unperturbed oscillator. Integrating Equ. (13) over one oscillation period τ = 2π/Ω we we obtain the change ∆ǫ = τ 0 dǫ of ǫ during the time τ , To compute the integrals, we assume that during this time, which at low friction is short compared to the time for energy relaxation, the particle performs an undisturbed harmonic oscillation evolving according to where the amplitude R of the oscillation is related to ǫ by R = 2/m(ǫ/Ω). The phase φ accounts for a possible phase shift with respect to the driving force, which is proportional to cos(Ω m t). Note that the oscillation frequency Ω is not necessarily the same as the frequency Ω 0 of the unperturbed harmonic oscillator. The central assumption, which allows to treat the motion of the system as that of an undisturbed oscillator during one oscillation period and eliminate the dependence on the rate of energy change on the phase space variables q and p by integration, is that the system evolves at nearly constant energy during one oscillation period. This condition is met if there is a separation of time scales between the time scale of the oscillation and the time scale for energy loss/gain. In other words, the relative change in energy ∆E/E occurring during one oscillation period should be much smaller than unity. The stochastic differential equation derived below for the time evolution of the energy, Equ. (27), provides a way to estimate for which ranges of the parameters Γ, η and ζ this condition holds. Analyzing each term on the right hand side of Equ. (27) individually, we find that the separation of time scale requires that Γ/Ω ≪ 1, ζ ≪ 1 and ηk B T eff /mΩ 2 ≪ 1, where T eff is the effective temperature of the oscillator.
Carrying out the integrals over t, the first three terms in Equ. (15) yield The change in ǫ resulting from the driving (fourth term in Equ. (15)) is given by ∆ǫ ′′ = − ǫζΩ 2 0 sin(π Ωm Ω ) cos(2φ) sin(π Ωm Ω ) − ( Ωm 2Ω ) sin(2φ) cos(π Ωm Ω ) Ωπ 4 − Ω 2 m Ω 2 τ. (18) This expression is independent of time only if after one oscillation period the relative phase of the oscillation with respect to the periodic driving force is the same as at the beginning of the period. For the parameters studied here and a modulation frequency of Ω m ≈ 2Ω 0 , the oscillator locks to the modulation and oscillates with Ω = Ω m /2. We limit our considerations to this case in the following. Carrying out the limit Ω m → 2Ω in the above equation, one finds Finally, the last term in Equ. (15), is a stochastic integral due to the noise term in the equations of motion. As a weighted sum of Gaussian random variables, ∆ǫ ′′′ is also a Gaussian random variable with mean and variance Thus, the random variable ∆ǫ ′′′ can be written in terms of the Wiener process as Putting things together, one obtains Since the oscillation period τ is short compared to all the time scale on which the energy changes, one can finally write the following stochastic differential equation for the square root of the energy ε The corresponding Fokker-Planck equation [25] governing the time evolution of the probability density function P ǫ (ǫ, t) is given by In writing these two equation we have implicitly assumed that the phase φ between the modulation and the particle oscillation is fixed (or at least that it changes only very slowly in time). As we will show below, this condition is met very well particularly at low friction. Equation (25) implies that the time evolution of ε can be viewed as a Brownian motion in the high friction limit under the influence of an external force. Note that due to the integration over one oscillation period, this equation has ǫ as its only time dependent variable while the dependence on other variables has been removed. In the following section we will use this equation to determine the energy distribution as well as the phase space distribution of the steady state generated by the parametric modulation and the feedback mechanism. Changing variables from ǫ to E = ǫ 2 and applying Ito's formula [24] yields the corresponding stochastic differential equation for the energy, In contrast do stochastic equation of motion for ǫ, here the noise is multiplicative, i.e., its amplitude is energy dependent. The corresponding Fokker-Planck equation for the probability density function P E (E, t) is given by

Energy distribution
The stochastic differential equation (25) has the form of the equation of motion describing the time evolution of a one-dimensional Brownian particle under the external force f (x) with large friction ν at temperature T , where x is the position of the Brownian particle. The motion resulting from this equation of motion is known to sample the Boltzmann-Gibbs distribution where β = 1/k B T is the reciprocal temperature and U(x) is the potential corresponding to the external force, f (x) = −dU/dx. By virtue of this isomorphism with over-damped Brownian motion, established by setting ν = 4/Γ and identifying ǫ with x, the determination of the energy in the nonequilibrium steady state of the driven oscillator turns into an equilibrium problem. One can then immediately infer that Equ. (25) samples the distribution where the potential generates the force acting on the variable ǫ. As a result, the systems samples the ǫ-distribution Note that a small friction Γ corresponds to large friction ν determining the time evolution of ǫ and, thus, the energy E of the oscillator. By a change of variables from ǫ to E, we finally obtain the probability density function of the energy E, The normalisation factor Z = P E (E)dE is given by where the function h(x) is defined as and erfc(x) is the complementary error function. Thus, the energy distribution is that of an equilibrium system with effective energy and configurational partition function Z. While the term proportional to E 2 is caused by the feedback cooling, the term proportional to E is affected only by the parametric modulation. According to Equ. (35), the energy distribution is Gaussian with a cutoff at E = 0. The maximum of the Gaussian is located at while its variance (neglecting the cutoff) is given by Hence, the width of the Gaussian does neither depend on the driving parameters nor on the phase φ.

Phase space distribution
Since for low friction the energy of the oscillator changes slowly, one can also obtain the full phase space density P qp (q, p) from the energy density P E (E). To determine the phase space density P qp (q, p), we consider the micro-canonical phase space distribution P mc (q, p;Ẽ) of the oscillator evolving at a given constant total energyẼ, where δ(x) is the Dirac delta function and we have denoted the fixed value of the energy withẼ to distinguish it from the energy function E(q, p), which depends on the position q and the momentum p. The normalising factor g(Ẽ) is the micro-canonical density of states, The phase space distribution of Equ. (41) would be observed for an oscillator evolving freely in the absence of feedback and without coupling to a heat bath. Since for the parameter ranges studied here the energy is essentially constant over many oscillation periods, the total phase space density P qp (q, p) can be written by averaging the microcanonical distribution over the energy distribution, This linear superposition of micro canonical distributions is valid as long as the energy changes slowly on the time scale of the oscillation period. For the low friction constants and the small feedback strength studied here this assumption is met even under nonequilibrium conditions. Carrying out the integral yields As further approximation, we now use the density of states g(E) = 2π/Ω 0 for the harmonic oscillator, thus neglecting the Duffing term of the potential in this part of the calculation, and obtain Inserting the energy distribution from Equ. (35) into this equation we finally find the phase distribution function Note, however, that while we have neglected the Duffing term in the expression for the density of states, it is included in the energy appearing in the argument of the exponential on the right hand side of the above equation.
From the phase space density P qp (q, p) one can obtain the distribution P q (q) of the position by integration over the momenta, In the absence of parametric modulation (ζ = 0), one finds by carrying out the integral where K 1/4 is a generalised Bessel function of the second kind. For simplicity, we have considered the case Ω 0 = Ω here. A similar expression can also be derived for the momentum distribution.

Relative entropy change
As shown recently, a fluctuation theorem holds for the relative entropy change ∆S for a system relaxing towards equilibrium starting from the non-equilibrium steady state prepared by feedback cooling and parametric driving [15]. In this process, the feedback and the driving are turned off during the relaxation such that the system evolves freely and the dynamics is microscopically reversible. The relative entropy change ∆S is defined as the logarithmic ratio of the probability P [u(t)] to observe a certain trajectory u(t) and the probability P [u * (t)] of the time reversed trajectory u * (t), Here, u(t) denotes an entire trajectory of length t including position and momentum of the oscillator and u * (t) denotes the trajectory that consist of the same states visited in reverse order with inverted momenta. Since during the relaxation detailed balance is obeyed, for the quantity ∆S a detailed fluctuation can be proven, where P t (∆S) is the probability density to observe the value ∆S at time t as determined over many repetitions of the relaxation experiment. For the relaxation process considered in Ref. [15] the relative entropy change is given by where ] is the energy absorbed by the bath during the relaxation, and E 0 and E t are the energy of the oscillator at time 0 and t, respectively. The quantity φ(q, p) is defined as as the logarithm of the stationary phase space distribution and ∆φ is the difference of φ at the beginning and the end of the trajectory, Hence, the relative entropy change ∆S depends on the state of the system at the beginning and the end of the trajectory. In general, the steady distribution P qp (q, p) necessary to compute ∆φ is unknown. However, from the distribution derived for our model, Equ. (46), we find that for the relaxation from a non-equilibrium steady state generated by nonlinear feedback and parametric modulation, the relative entropy change is given by Thus, our stochastic model allows us to express the relative entropy change during a relaxation trajectory in terms of the energy at the beginning and the end of that trajectory. Note that since no work is performed on the system, the heat Q h exchanged along a trajectory equals the energy lost by the system. Thus, in the absence of nonlinear feedback cooling, the relative entropy change is proportional to the heat and resembles relaxation form a thermal bath with effective temperature T eff = T /(1 − ζΩ 2 0 sin(2φ) 2ΓΩ ). By choosing parameters, one can therefore switch from a purely thermal situation with the phase space distribution of a harmonic oscillator (but with changed temperature) to a truly non-equilibrium steady-state with non-linear effects controlled by the feedback parameter η.

Quadratures
Parametrically driven nano-mechanical oscillators have been shown to support classical squeezed states in which the amplitude of the vibration in one phase is reduced with respect to the thermal equilibrium amplitude. To probe our oscillator for squeezed states we analyse its motion in terms of the so-called quadratures. For the oscillator driven by the parametric modulation F drive = ζmΩ 2 0 cos(Ω m t)q, we write the time evolution of the oscillator position as where Ω is the frequency of the particle oscillating at half the frequency of the driving, Ω = Ω m /2. Here, R(t) and φ(t) are the amplitude and the phase of the particle, respectively, and the phase is measured with respect to the driving signal. Using the addition theorem for the sine-function, Equ. (55) can be written as the sum of two contributions, one in-phase with the driving signal and one out-of-phase, where the second line defines the in-phase component X(t) = R(t) cos φ(t) and the quadrature Y (t) = R(t) sin φ(t). Together, X and Y are referred to as the quadratures. The quadratures can be computed from the time evolution of the position q(t) and the momentum p(t). The momentum of the particle is given by: where we neglected the time derivatives of the amplitude and phase, since for an oscillator at low friction both the amplitude and the phase vary slowly in time.
Combining this equation with Equ. (56) yields corresponding to transformation to a coordinate system that rotates clockwise with frequency Ω with respect to the (q, p/mΩ)-plane [26,27]. In this coordinate system, a sinusoidal oscillation of frequency Ω is represented by a static point. Note that the amplitude and phase can be expressed in terms of the quadratures, and that is the energy of a harmonic oscillator with frequency Ω.

Simulations
In this section we verify the analytical expressions for the distributions of energy and positions by comparing them with simulation results. The simulations were performed for parameter values close to those of the experiments, which we will present and discuss subsequently.

Simulation methods
In our simulations, we integrated the Langevin equation of motion with the OVRVO algorithm of Sivak, Chodera and Crooks [28], which can be viewed as a stochastic generalisation of the velocity Verlet algorithm for deterministic dynamics [29]. This discrete time integration scheme uses a time step rescaling in the deterministic update step for positions and momenta to satisfy a number of desiderata proposed in the literature for stochastic integrators [30]. In all simulations we used a time step of ∆t = 0.01 in reduced units. This time step is about 1/628 of the oscillation period. Test runs carried out with smaller time steps (∆t = 0.001) yielded identical results up to statistical errors. In most cases, the total simulation time was t = 10 7 corresponding to about 3 ×10 6 modulation cycles. For some parameters we carried out longer simulations of up to 3 × 10 10 steps corresponding to a total simulation time of t = 3 × 10 8 . All simulations were carried out for k B T = 1, m = 1, and k = 1.
To facilitate comparison of the results of theory/simulation and experiments, in the following we use the thermal energy E = k B T , the inverse frequency T = 1/Ω 0 and the particle mass M = m as our basic units of energy, time and mass, respectively. Accordingly, distances are measured in units of L = (1/Ω 0 ) k B T /m and velocities in units of V = k B T /m. Hence, the unit of length is given by the variance of the position of the harmonic oscillator, q 2 = k B T /mΩ 2 0 = L 2 and the unit of energy is the average energy of the harmonic oscillator E = k B T = E. The friction constant is given in units of Ω 0 such that it equals the inverse of the quality factor, Q = Ω 0 /Γ = 1/ΓT . The feedback strength η and the Duffing coefficient ξ have the dimension of 1/area and are measured in units of 1/L 2 . The modulation depth ζ is dimensionless. In the following, we use reduced units in which E = T = M = 1.  We first consider the oscillator without parametric modulation (ζ = 0.0) but subjected to feedback cooling. Without driving, the phase φ is not a relevant parameter and the expression for the energy distribution simplifies considerably,

Oscillator with feedback cooling but without parametric modulation
where we have assumed that the particle oscillates with Ω = Ω 0 . The first term in the exponential is the same as that of the uncooled oscillator, but the second term proportional to E 2 is due to the feedback loop and strongly penalises high energy states thereby cooling the system. The cooling effect is stronger for weak friction Γ and small frequencies Ω 0 . Several energy distributions obtained from simulations together with the corresponding predictions of Equ. (61) are shown in the left panel of Fig. 1. The simulations were carried out for a friction of Γ = 0.0001 and a Duffing parameter of ξ = −0.022. Without feedback, η = 0, the energy distribution is exponential, but for η > 0 the E 2 term caused by the feedback suppresses high energies leading to a parabolic shape of the distribution in the logarithmic representation. In all cases, the theoretical predictions agree very well with the simulation results. Positions distributions for the same set of parameters are shown in the right panel of Fig. 1. While without feedback the position distribution is Gaussian, the feedback quenches large deviations leading to a narrowing of the distributions. Also in the case of the position distributions the agreement between theory and simulation is excellent.  We next turn to the oscillator with parametric driving but without feedback cooling. In this case, the energy deposited in the system by the modulation is removed only by the coupling to the gas as quantified by the friction constant Γ. If the particle oscillation is locked to the driving with a fixed phase φ, the resulting energy distribution following from Equ. (35) is expected to be exponential,

Oscillator with parametric modulation but without feedback cooling
where we have assumed that the modulation frequency is Ω m = 2Ω 0 . For a vanishing Duffing parameter ξ = 0.0, i.e., for a perfectly harmonic trap, the phase is expected to be φ = −π/4 in the absence of thermal fluctuations [31]. If this is the case, the decay constant of the exponential is β(1 − ζΩ 0 /2Γ). Hence, the decay constant is positive only for ζ < 2Γ/Ω 0 . If the modulation depth ζ exceeds this limit, the friction cannot remove the energy pumped into the oscillator by the modulation such that the oscillator energy keeps growing preventing the system from settling in a steady state. We indeed find in our simulations that for ζ > 2Γ/Ω 0 the energy continuously increases. For weak driving, on the other hand, the energy distribution is expected to be exponential with the decay constant predicted by Equ. (35). Several energy distributions for this case are shown in Fig. 2. Note that we performed these calculations for a relatively large friction constant of Γ = 0.01, because for lower friction it takes exceedingly long to sample all relevant energies. For weak driving, ξ = 0.001 (red symbols), the energy distribution is exponential as predicted by the theory. The negative slope of this distribution in the logarithmic representation is, however, slightly too large. The reason for this discrepancy is that the oscillation does not lock to the parametric driving as can bee seen in the distribution of the phase φ shown in the right panel of Fig. 2. The theory developed above, on the other hand, assumes a fixed phase of φ = −π/4 (for ξ = 0). For ξ = 0.001, the phase distribution is essentially flat implying that there is no preferred phase. As a consequence, essentially no heating occurs and the energy distribution is indistinguishable from the equilibrium distribution (black symbols). As the strength of the parametric driving is increased, a pronounced phase relation between driving and oscillation develops and two distinct peaks appear in the phase distribution at equivalent positions, one at φ = −π/4 and one at φ = −π/4 + π. Since the phase relation is more pronounced at high energies, in this regime the energy distributions shown in the left panel of Fig. 2 converge to the form predicted by theory. In the figure, the theoretical distributions are indicated by lines with logarithmic slope of −β(1 − ζΩ 0 /2Γ). For low energies, the phase relation is lost and the energy distributions have the logarithmic slope of the equilibrium distribution. Thus, the energy injected into the system by the parametric driving results in a longer tail in the energy distribution where it has the right phase relationship with the oscillation. In contrast at low energies, the form of the distribution is essentially unchanged with respect to the equilibrium distribution.  Next, we consider the oscillator with parametric driving and feedback cooling. To understand the energy distributions for this case, we first take a closer look at the statistics of the phase φ. In the derivation of the analytical energy distribution, Equ. (35), we have assumed a fixed phase φ between the modulation and the particle  oscillation. In practice, however, the phase φ follows a statistical distribution with a position and width that depend on the parameters, particularly on the Duffing parameter ξ and the friction constant Γ. Several distributions of the phase obtained from our simulations for Γ and ξ are shown in Fig. 3. These simulations were carried out for a modulation depth of ζ = 0.03 and and a feedback strength of η = 0.022, because these values can be realised in experiments. For all parameters considered here, the phase distributions are strongly peaked at a particular phase. The peaks are narrow for small friction and small Duffing parameters and broaden for increasing friction and non-linearity. Note that the Duffing parameters considered here are negative because the non-linearity is due to the shape of the focal intensity distribution, which is approximately Gaussian [32]. Without non-linearity, ξ = 0.0, the peak is located at φ = −π/4 for all values of the friction constant. As one turns on the non-linearity by making the Duffing parameter more negative, the peaks become broader and shift towards more negative values.

Oscillator with feedback cooling and parametric driving
A closer analysis of how the phase depends on the Duffing parameter is shown in Fig. 4. The left panel of the figure shows the positions of the maximum of the phase distribution. i.e., the most likely phase φ max , as a function of the Duffing parameter ξ for different friction constants Γ. As can be inferred from the figure, the most likely phase φ max determined from the simulations (symbols) follows exactly the form predicted by secular perturbation theory [31] (solid lines). While this theory neglects thermal fluctuations and cannot predict the entire phase distribution, it yields an accurate location of the maximum.
Due to the thermal fluctuations, which lead to a broadening of the phase distribution, the oscillator might entirely loose the lock with the driving modulation and regain it only after falling behind by one entire turn of 2π. For the lowest friction studied here this never happens during a simulation of total time t = 10 7 , but for higher frictions, and in particular for large Duffing parameters, the oscillation may fall behind the parametric modulation several times. The number of times this occurs in the course of the simulations is shown in the right panel of Fig. 4 for different friction constants as a function of ξ.
We now compare the energy distribution determined in our simulations for the oscillator with parametric driving and feedback cooling with the theoretical prediction of Equ. (35). To do that, we identify the phase φ occurring in the theoretical expression with the most likely phase φ max determined in the simulations. Energy distributions obtained for friction constants ranging from Γ = 10 −5 to Γ = 10 −3 are shown in the left panel of Fig. 5. In all cases, the system was driven at Ω m = 2Ω 0 and the Duffing parameter, the feedback strength and the modulation depth were ξ = −0.022, η = 0.022, ζ = 0.03, respectively. While for high friction the theoretical predictions deviate considerably from the energy distributions determined in the simulations, most likely due to the lack of a stable phase relation, very good agreement is obtained for low friction, where phase distributions are strongly peaked. This excellent correspondence is confirmed by the energy distributions shown along with theoretical predictions in the right panel of Fig. 5 for different Duffing parameters at low friction. Thus, the position and the width of the energy distribution in the non-equilibrium steady state generated by driving and cooling at the same time can indeed be controlled independently by an appropriate choice of parameters.

Quadratures
Finally, we take a look at the distribution of the quadratures X and Y for different driving frequencies. Scatter plots of the quadratures obtained at different driving frequencies and for different values of the friction constant are shown in Fig. 6. From left to right, the driving frequency Ω m is slightly below 2Ω 0 , equal to 2Ω 0 and slightly above 2Ω 0 . As in previous simulations, the parameters were ξ = −0.022, η = 0.022, Figure 6. Scatter plot of the quadratures X and Y for the friction constants Γ = 0.00001 (black), 0.0001 (red) and 0.01 (blue) for Ω m = 1.98 Ω 0 (left), Ω m = 2 Ω 0 (center) and Ω m = 2.02 Ω 0 (right). The simulations were carried out for η = 0.022, ξ = −0.022, and ζ = 0.03. and ζ = 0.03. At Ω m = 2Ω and low friction the quadratures of the driven system are Gaussian with equal width along the two quadrature axes. Thus, they resemble a thermal state, albeit, displaced from the origin. In contrast, for driving frequencies off 2Ω 0 , the distributions are deformed, indicating the occurrence of classical squeezing.

Experiments
In this section, we discuss how to retrieve the energy and phase of a trapped nanoparticle from discrete measurements of the particle positions. From the retrieved energies and phases we reconstruct the energy and phase distributions and compare them to the theory and simulation results presented in the previous sections. This allows us to extract the experimental parameters, which are detailed in Table 1. While the maxima of the distributions are in good agreement with our theory and simulations, the width of the experimental distributions is significantly broader due to parameter fluctuations not taken into account in the theoretical considerations.

Experimental configuration
In our experiments we use a silica nanoparticle trapped at the focus of a single beam optical tweezers. The optical tweezer is formed by a 1064 nm laser beam (∼ 35mW) focused by a NA = 0.9 objective, which is mounted inside a vacuum chamber. The particle motion is recorded with an additional colinear laser (780nm) and three balanced photodetectors. A home-built electronic circuit is used to generate the feedback signal (η), while a frequency generator serves as the parametric modulation signal (ζ). The approximately Gaussian shape of the optical potential is responsible for the trap anharmonicity (ξ) [32]. The detectors and the size of the nanoparticle are calibrated from measurements of the power spectral density of the particle motion at 5.1mBar. At this pressure the Q-factor is high enough to resolve the three spatial modes, while broadening effects due to nonlinear mode coupling are negligible [32]. For further details of the experimental configuration and calibration procedure see Refs. [33,34] measurements are carried out at 1.2 × 10 −5 mBar. While our theoretical model is one-dimensional, the particle in the experiment moves in three dimensions along three main axes. The three axes are determined by the symmetry of the laser focus. However, there is no direct coupling between the three spatial modes. In addition, feedback cooling reduces the amplitude such that also the nonlinear coupling becomes very weak. Therefore, our one-dimensional model is a very good approximation for the particle motion along one of the three main axes.

Amplitude and phase estimation
The particle oscillation frequencies along the three main axes are well separated and dont overlap. Therefore, we can apply the maximum likelihood estimation for a single tone signal, that is a signal containing only one frequency component. The maximum likelihood estimation of the oscillation amplitude and phase of a single tone signal q(t) is given by [35] R ML = |A q (Ω)| (63) where Ω/2π is the estimated frequency of the signal, t 0 is the time origin and is the discrete Fourier transform of q evaluated at ω. Here, q n = q(t n ) is the measurement sample of the time trace at time t n = t 0 + n∆t, N is the number of samples entering the estimation and ∆t is the sampling interval. The estimation of the amplitude and the phase relies on precise estimation of the frequency Ω. We estimate Ω by maximising (65) with respect to ω, i.e. A(Ω) = max(A(ω)). The width of the function A(ω), and thereby our ability to localise the maximum, depends on the length of the time trace q(t). Therefore, we use a long time trace measured over T meas. = 0.1 s and sampled at 625 kilosamples/second to estimate Ω. Subsequently, we use that value of Ω and Equs.
(63) and (64) to estimate the instantaneous amplitude and phase from short parts of that same time trace. The short parts of the time trace contain N = 160 samples, corresponding to an integration over 32 particle oscillations. This constitutes a good compromise between sufficient data points for an accurate estimation of R and φ, and fast time resolution to resolve the dynamics of the energy and phase fluctuations. Note that maximising (65) allows us to estimate the frequency with much better accuracy than 1/T meas. . The absolute phase of a harmonic oscillator is a time delay with respect to some time reference. Without such a time reference the absolute phase is arbitrary and has no meaning. However, the relative phase between two oscillators is meaningful, because one oscillator serves as a time reference to determine the phase of the other oscillator with respect the first oscillator. Formally, this is expressed as where A p and A m are the Fourier transforms of the two signals, respectively (c.f. (65)), and k is an integer which takes into account that the phase is only determined up to modulo 2π. Note that the exponent Ω p /Ω m takes care that (66) does not depend on t 0 . Without loss of generality, we set φ m = 0, i.e. we choose our time origin such that it coincides with a maximum of the signal with frequency Ω m . For the special case of a parametrically driven particle, which oscillates at half the frequency of the parametric modulation (Ω m = 2Ω p ), we get ∆φ = φ p − kπ. Therefore, the above method allows to estimate the relative phase between the particle oscillation and the parametric modulation up to a multiple of π.

Parameter estimation
We measure the distribution of the energy and phase for modulation at Ω m /2π = 247, 248, 249, and 250 kHz. Each distribution is obtained from 100 time traces of 0.1s duration. Fig. 7 shows the maximum values of the energy and phase distributions shown in Fig. 8 and a fit to secular perturbation theory [34,31]. While independent fits to the energy and phase, shown in blue and red, respectively, yield excellent agreement with the theoretical model, we cannot fit a set of parameters that would agree with both the energy and the phase. Note that the phase fit includes a constant phase offset φ 0 = 50 • to account for the finite response time of the intensity modulator and delays in the electronics. Averaging the results from the independent fits to energy and phase yields ξ = −5.4 ± 1.1 µm −2 , η = 3.9 ± 1.3 µm −2 and ζ = 16.1 ± 5.7 × 10 −3 . The theoretical curve for the parameters obtained by the energy and phase is shown in green together with numerical simulations using the parameters summarised in Table 1.  Figure 7.
Left: Most likely energy. Right: Most likely phase. The black and green circles are the experimental data points and simulation results, respectively. The blue and red solid lines are the theoretical predictions for parameters obtained from independent fits to the energy and phase, respectively and the green solid line is the theoretical prediction for the averaged parameters. The main uncertainty in the determination of the experimental parameters arises from the estimation of the particle mass and the resulting uncertainty in the voltage calibration and from parameter fluctuations, which we discuss in the next section. As an independent measurement, we also measure the energy distribution without parametric modulation (ζ = 0). A fit of the energy distribution to Eq. (61) yields η = 4.5±0.9 µm −2 , in good agreement with the previously determined value. Fig. 8 shows the experimental energy and phase distributions fitted with a Gaussian. As predicted by our theory and simulations, the distributions are Gaussian and their widths depend only weakly on the modulation frequency. Fig. 9 shows the widths of the distributions obtained from the Gaussian fits in Fig. 8 and from numerical simulations. For comparison, we also show the theoretical prediction according to Equ. (40). The broadening of the distributions has two contributions, thermal motion and parameter fluctuations.  Thermal motion of the resonator, caused by residual air molecules, enters directly as a random white noise, which we considered in our theoretical model. In addition, it enters indirectly through amplitude-phase conversion [36]. The latter contribution has not been considered in our theoretical model but is naturally present in the numerical simulations. Amplitude-phase conversion refers to the interdependence of energy and phase (c.f. (39)). Therefore, fluctuations in the phase cause fluctuations in the energy and vice versa. This leads to a broadening of the distributions near the instability boundaries, where the deviation of the numerical simulation from our model is largest. Within this range, on the other hand, this interplay manifests itself as sidebands in the power spectral density of the particle position [34].

Distributions
In addition to Brownian motion, parameter fluctuations broaden the experimental distributions [37]. The experimental parameters fluctuate due to laser intensity and polarisation fluctuations and also due to the nonlinear coupling with the other two degrees of freedom, which were not considered in our model [32,34]. Noise in the feedback electronics and modulator gives rise to further broadening. In general, broadening due to fluctuating parameters dominates broadening due to Brownian motion. As a consequence, the measured width of the energy and phase distributions σ E = 78 ± 3 × 10 −3 k B T and σ φ = 1.7 ± 0.1 × 10 −3 π, respectively, are approximately one order of magnitude larger than the theoretical values 5.1 × 10 −3 k B T and 0.15 × 10 −3 π, averaged over the range of detunings of the experimental data. To identify the noise sources responsible for the deviation from theory one can deliberately introduce noise and systematically study its effect on the measurement outcome.

Conclusion
We have developed a stochastic model for the dynamics of the energy of a nonlinear nanomechanical oscillator subject to parametric modulation and nonlinear damping.
Under these conditions the oscillator attains a non-equilibrium steady state. Our model allows us to predict the energy distribution of the steady state. The steady state distribution is intimately related to fluctuation theorems, which describe the statistical properties of the system for transitions between different states [15]. Consequently, our work opens the door to test these fluctuation theorems in different scenarios.
We confirmed the validity of the model by extensive numerical simulations and found excellent agreement with our theory. In addition, we performed experiments with a levitated nanoparticle. While the measured mean energy and phase are in close agreement with the numerical simulations, their distributions are broadened due to parameter fluctuations that are not accounted for in the theory and are subject to further investigation. Besides quantifying additional noise sources experimentally, future work includes the development of a more generalised model including a stochastic model for the phase and incorporating other white and non-white noise sources, resulting from fluctuating parameters [23].