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

We consider a particle of mass m oscillating in a trap with a Duffing potential

$$\begin{eqnarray}&&V(q)=\frac{1}{2}k{{q}^{2}}+\frac{1}{4}\xi k{{q}^{4}},\end{eqnarray} \tag{ 1 }$$

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 ${{\Omega }_{0}}=\sqrt{k/m}$ of the harmonic case, the total energy of the oscillator is given by

$$\begin{eqnarray}&&E(q,p)=\frac{1}{2}m\Omega _{0}^{2}{{q}^{2}}+\frac{1}{4}\xi m\Omega _{0}^{2}{{q}^{4}}+\frac{{{p}^{2}}}{2\;m},\end{eqnarray} \tag{ 2 }$$

where $p=m\dot{q}$ is the momentum of the particle. The force due to the trap is hence given by

$$\begin{eqnarray}&&{{F}_{{\rm trap}}}=-m\Omega _{0}^{2}q-\xi m\Omega _{0}^{2}{{q}^{3}}.\end{eqnarray} \tag{ 3 }$$

Since the particle is immersed in a low density gas of temperature T, it experiences also a frictional force

$$\begin{eqnarray}&&{{F}_{{\rm friction}}}=-\Gamma p\end{eqnarray} \tag{ 4 }$$

and the related fluctuating random force

$$\begin{eqnarray}&&{{F}_{{\rm random}}}=\sqrt{2\;m\Gamma {{k}_{{\rm B}}}T}w(t),\end{eqnarray} \tag{ 5 }$$

where Γ is the friction constant, ${{k}_{{\rm B}}}$ is the Boltzmann constant and w(t) is white noise. A feedback of strength η, acting on the particle with force

$$\begin{eqnarray}&&{{F}_{{\rm feedback}}}=-{{\Omega }_{0}}\eta {{q}^{2}}p,\end{eqnarray} \tag{ 6 }$$

is 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 ${{\Omega }_{m}}$ leading to the force

$$\begin{eqnarray}&&{{F}_{{\rm drive}}}=\zeta m\Omega _{0}^{2}{\rm cos} ({{\Omega }_{m}}t)q,\end{eqnarray} \tag{ 7 }$$

where the modulation depth ζ determines the intensity of the parametric driving. Taken together, these forces yield the following stochastic equations of motion for the particle in the trap,

$$\begin{eqnarray}&&{\rm d}q=\frac{p}{m}{\rm d}t,\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}\begin{array}{rcl} {\rm d}p & = & \left[ -m\Omega _{0}^{2}q-\xi m\Omega _{0}^{2}{{q}^{3}}-\Gamma p-{{\Omega }_{0}}\eta {{q}^{2}}p+\zeta m\Omega _{0}^{2}{\rm cos} ({{\Omega }_{m}}t)q \right]{\rm d}t \\ {} & {} & +\sqrt{2\;m\Gamma {{k}_{{\rm B}}}T}\;{\rm d}W. \\ \end{array}\end{eqnarray} \tag{ 9 }$$

Here, W(t) is the Wiener process with

$$\begin{eqnarray}&&\langle W(t)\rangle =0,\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\langle W(t)W(t^{\prime} )\rangle ={\rm min} (t,t^{\prime} ).\end{eqnarray} \tag{ 11 }$$

Note that $\langle {{W}^{2}}(t)\rangle =t$ for any time $t\geqslant 0$ and, thus, for an infinitesimal time interval ${\rm d}t$ one has $\langle {{({\rm dW})}^{2}}\rangle ={\rm d}t$. The white noise w(t) appearing in the random force can be viewed as the time derivative of the Wiener process, $w(t)={\rm d}W(t)/{\rm d}t$.

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,

$$\begin{eqnarray}&&\epsilon (q,p)=\sqrt{E(q,p)}.\end{eqnarray} \tag{ 12 }$$

Applying Ito's formula [24] for the change of variables to $\epsilon (q,p)$ we find that the change ${\rm d}\epsilon$ during a short time interval is given by

$$\begin{eqnarray}&&{\rm d}\epsilon =\left[ \frac{pF(q,p,t)}{2\;m\epsilon }+\frac{\Gamma {{k}_{{\rm B}}}T}{2\epsilon }\left( 1-\frac{{{p}^{2}}}{2\;m{{\epsilon }^{2}}} \right) \right]{\rm d}t+\sqrt{2\;m\Gamma {{k}_{{\rm B}}}T}\frac{p}{2\;m\epsilon }{\rm d}W,\end{eqnarray} \tag{ 13 }$$

where

$$\begin{eqnarray}&&F(q,p,t)=-\Gamma p-{{\Omega }_{0}}\eta {{q}^{2}}p+\zeta m\Omega _{0}^{2}{\rm cos} ({{\Omega }_{m}}t)q\end{eqnarray} \tag{ 14 }$$

is the sum of the non-conservative forces consisting of the frictional force ${{F}_{{\rm friction}}}$, the feedback force ${{F}_{{\rm feedback}}}$ and the driving force ${{F}_{{\rm 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 , equation (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 ${{\Omega }_{0}}$ of the unperturbed oscillator. Integrating equation (13) over one oscillation period $\tau =2\pi /\Omega$ we obtain the change $\Delta \epsilon =\int _{0}^{\tau }{\rm d}\epsilon$ of during the time τ,

$$\begin{eqnarray}\begin{array}{rcl} \Delta \epsilon & = & -\Gamma \int _{0}^{\tau }\frac{{{p}^{2}}}{2\;m\epsilon }{\rm d}t-{{\Omega }_{0}}\eta \int _{0}^{\tau }\frac{{{q}^{2}}{{p}^{2}}}{2\;m\epsilon }{\rm d}t+\Gamma {{k}_{{\rm B}}}T\int _{0}^{\tau }\frac{1}{2\epsilon }\left( 1-\frac{{{p}^{2}}}{2\;m{{\epsilon }^{2}}} \right){\rm d}t \\ {} & {} & +\zeta m\Omega _{0}^{2}\int _{0}^{\tau }\frac{{\rm cos} ({{\Omega }_{m}}t)qp}{2\;m\epsilon }{\rm d}t+\sqrt{2\;m\Gamma {{k}_{{\rm B}}}T}\int _{0}^{\tau }\frac{p}{2\;m\epsilon }{\rm d}W. \\ \end{array}\end{eqnarray} \tag{ 15 }$$

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

$$\begin{eqnarray}&&q(t)=R{\rm cos} (\Omega t+\phi )\qquad \qquad p(t)=-m\Omega R{\rm sin} (\Omega t+\phi ),\end{eqnarray} \tag{ 16 }$$

where the amplitude R of the oscillation is related to by $R=\sqrt{2/m}(\epsilon /\Omega )$. The phase ϕ accounts for a possible phase shift with respect to the driving force, which is proportional to ${\rm cos} ({{\Omega }_{m}}t)$.

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 $\Delta 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, equation (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 equation (27) individually, we find that the separation of time scale requires that $\Gamma /\Omega \ll 1$, $\zeta \ll 1$ and $\eta {{k}_{{\rm B}}}{{T}_{{\rm eff}}}/m{{\Omega }^{2}}\ll 1$, where T_eff is the effective temperature of the oscillator given by the average energy, ${{k}_{{\rm B}}}{{T}_{{\rm eff}}}=\langle E\rangle$.

Carrying out the integrals over t, the first three terms in equation (15) yield

$$\begin{eqnarray}&&\Delta {{\epsilon }^{\prime }}=-\frac{\Gamma \epsilon }{2}\tau -\frac{\eta {{\epsilon }^{3}}{{\Omega }_{0}}}{4\;m{{\Omega }^{2}}}\tau +\frac{\Gamma {{k}_{{\rm B}}}T}{4\epsilon }\tau .\end{eqnarray} \tag{ 17 }$$

The change in resulting from the driving (fourth term in equation (15)) is given by

$$\begin{eqnarray}&&\Delta {{\epsilon }^{\prime\prime }}=-\frac{\epsilon \zeta \Omega _{0}^{2}{\rm sin} (\pi \frac{{{\Omega }_{m}}}{\Omega })\left[ {\rm cos} (2\phi ){\rm sin} (\pi \frac{{{\Omega }_{m}}}{\Omega })-(\frac{{{\Omega }_{m}}}{2\Omega }){\rm sin} (2\phi ){\rm cos} (\pi \frac{{{\Omega }_{m}}}{\Omega }) \right]}{\Omega \pi \left( 4-\frac{\Omega _{m}^{2}}{{{\Omega }^{2}}} \right)}\tau .\end{eqnarray} \tag{ 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 ${{\Omega }_{m}}\approx 2{{\Omega }_{0}}$, the oscillator locks to the modulation and oscillates with $\Omega ={{\Omega }_{m}}/2$. We limit our considerations to this case in the following. Carrying out the limit ${{\Omega }_{m}}\to 2\Omega$ in the above equation, one finds

$$\begin{eqnarray}&&\Delta {{\epsilon }^{\prime\prime }}=-\frac{\epsilon \zeta \Omega _{0}^{2}{\rm sin} (2\phi )}{4\Omega }\tau .\end{eqnarray} \tag{ 19 }$$

Finally, the last term in equation (15),

$$\begin{eqnarray}&&\Delta {{\epsilon }^{}}=\sqrt{2\;m\Gamma {{k}_{{\rm B}}}{{T}_{0}}}\int _{0}^{\tau }\frac{p}{2\;m\epsilon }{\rm d}W,\end{eqnarray} \tag{ 20 }$$

is a stochastic integral due to the noise term in the equations of motion. As a weighted sum of Gaussian random variables, $\Delta {{\epsilon }^{}}$ is also a Gaussian random variable with mean

$$\begin{eqnarray}&&\langle \Delta {{\epsilon }^{}}\rangle =\sqrt{2\;m\Gamma {{k}_{{\rm B}}}T}\int _{0}^{\tau }\frac{p}{2\;m\epsilon }\langle {\rm d}W\rangle =0\end{eqnarray} \tag{ 21 }$$

and variance

$$\begin{eqnarray}\begin{array}{rcl} \langle {{(\Delta {{\epsilon }^{}})}^{2}}\rangle & = & 2\;m\Gamma {{k}_{{\rm B}}}T\int _{0}^{\tau }\int _{0}^{\tau }\frac{p(t)p(t^{\prime} )}{4\;{{m}^{2}}{{\epsilon }^{2}}}\langle {\rm d}W{\rm d}W^{\prime} \rangle \\ {} & = & 2\;m\Gamma {{k}_{{\rm B}}}T\int _{0}^{\tau }\frac{{{p}^{2}}}{4\;{{m}^{2}}{{\epsilon }^{2}}}{\rm d}t=\frac{\Gamma {{k}_{{\rm B}}}T}{2}\tau . \\ \end{array}\end{eqnarray} \tag{ 22 }$$

Thus, the random variable $\Delta {{\epsilon }^{}}$ can be written in terms of the Wiener process as

$$\begin{eqnarray}&&\Delta {{\epsilon }^{}}=\sqrt{\frac{\Gamma {{k}_{{\rm B}}}T}{2}}W(\tau ).\end{eqnarray} \tag{ 23 }$$

Putting things together, one obtains

$$\begin{eqnarray}&&\Delta \epsilon =\left[ -\frac{\Gamma \epsilon }{2}-\frac{\eta {{\epsilon }^{3}}{{\Omega }_{0}}}{4\;m{{\Omega }^{2}}}+\frac{\Gamma {{k}_{{\rm B}}}T}{4\epsilon }-\frac{\epsilon \zeta \Omega _{0}^{2}{\rm sin} (2\phi )}{4\Omega } \right]\tau +\sqrt{\frac{\Gamma {{k}_{{\rm B}}}T}{2}}W(\tau ).\end{eqnarray} \tag{ 24 }$$

Since the oscillation period τ is short compared to the time scale on which the energy changes, one can finally write the following stochastic differential equation for the square root of the energy

$$\begin{eqnarray}&&{\rm d}\epsilon =\left[ -\frac{\Gamma \epsilon }{2}-\frac{\eta {{\Omega }_{0}}{{\epsilon }^{3}}}{4\;m{{\Omega }^{2}}}+\frac{\Gamma {{k}_{{\rm B}}}T}{4\epsilon }-\frac{\epsilon \zeta \Omega _{0}^{2}{\rm sin} (2\phi )}{4\Omega } \right]{\rm d}t+\sqrt{\frac{\Gamma {{k}_{{\rm B}}}T}{2}}{\rm d}W.\end{eqnarray} \tag{ 25 }$$

The corresponding Fokker–Planck equation [25] governing the time evolution of the probability density function ${{P}_{\epsilon }}(\epsilon ,t)$ is given by

$$\begin{eqnarray}\begin{array}{rcl} \frac{\partial {{P}_{\epsilon }}(\epsilon ,t)}{\partial t} & = & \frac{\partial }{\partial \epsilon }\left[ \frac{\Gamma \epsilon }{2}+\frac{\eta {{\Omega }_{0}}{{\epsilon }^{3}}}{4\;m{{\Omega }^{2}}}-\frac{\Gamma {{k}_{{\rm B}}}T}{4\epsilon }+\frac{\epsilon \zeta \Omega _{0}^{2}{\rm sin} (2\phi )}{4\Omega } \right]{{P}_{\epsilon }}(\epsilon ,t) \\ {} & {} & +\frac{\Gamma {{k}_{{\rm B}}}T}{4}\frac{{{\partial }^{2}}}{\partial {{\epsilon }^{2}}}{{P}_{\epsilon }}(\epsilon ,t). \\ \end{array}\end{eqnarray} \tag{ 26 }$$

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 holds accurately 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={{\epsilon }^{2}}$ and applying Ito's formula [24] yields the corresponding stochastic differential equation for the energy

$$\begin{eqnarray}\begin{array}{rcl} {\rm d}E & = & \left[ -\Gamma (E-{{k}_{{\rm B}}}T)-\frac{\eta {{\Omega }_{0}}{{E}^{2}}}{2\;m{{\Omega }^{2}}}-\frac{E\zeta \Omega _{0}^{2}{\rm sin} (2\phi )}{2\Omega } \right]{\rm d}t \\ {} & {} & +\sqrt{2E\Gamma {{k}_{{\rm B}}}T}{\rm d}W. \\ \end{array}\end{eqnarray} \tag{ 27 }$$

In contrast to the 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

$$\begin{eqnarray}\begin{array}{rcl} \frac{\partial {{P}_{E}}(E,t)}{\partial t} & = & \frac{\partial }{\partial E}\left[ \Gamma (E-{{k}_{{\rm B}}}T)+\frac{\eta {{\Omega }_{0}}{{E}^{2}}}{2\;m{{\Omega }^{2}}}+\frac{E\zeta \Omega _{0}^{2}{\rm sin} (2\phi )}{2\Omega } \right]{{P}_{E}}(E,t) \\ {} & {} & +\Gamma {{k}_{{\rm B}}}T\frac{{{\partial }^{2}}}{\partial {{E}^{2}}}E{{P}_{E}}(E,t). \\ \end{array}\end{eqnarray} \tag{ 28 }$$

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,

$$\begin{eqnarray}&&{\rm d}x=\frac{1}{\nu }f(x){\rm d}t+\sqrt{\frac{2{{k}_{{\rm B}}}T}{\nu }}{\rm dW},\end{eqnarray} \tag{ 29 }$$

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

$$\begin{eqnarray}&&{{P}_{x}}(x)\propto {\rm exp} \{-\beta U(x)\},\end{eqnarray} \tag{ 30 }$$

where $\beta =1/{{k}_{{\rm B}}}T$ is the reciprocal temperature and U(x) is the potential corresponding to the external force, $f(x)=-{\rm d}U/{\rm d}x$.

By virtue of this isomorphism with over-damped Brownian motion, established by setting $\nu =4/\Gamma$ and identifying with x, the determination of the energy in the non-equilibrium steady state of the driven oscillator turns into an equilibrium problem. One can then immediately infer that equation (25) samples the distribution

$$\begin{eqnarray}&&{{P}_{\epsilon }}(\epsilon )\propto {\rm exp} \{-\beta U(\epsilon )\},\end{eqnarray} \tag{ 31 }$$

where the potential

$$\begin{eqnarray}&&U(\epsilon )={{\epsilon }^{2}}+\frac{\eta {{\Omega }_{0}}{{\epsilon }^{4}}}{4\;m\Gamma {{\Omega }^{2}}}-{{k}_{{\rm B}}}T{\rm ln} \epsilon +\frac{{{\epsilon }^{2}}\zeta \Omega _{0}^{2}{\rm sin} (2\phi )}{2\Gamma \Omega }\end{eqnarray} \tag{ 32 }$$

generates the force

$$\begin{eqnarray}&&f(\epsilon )=-\frac{{\rm d}U(\epsilon )}{{\rm d}\epsilon }=-2\epsilon -\frac{\eta {{\Omega }_{0}}{{\epsilon }^{3}}}{m\Gamma {{\Omega }^{2}}}+\frac{{{k}_{{\rm B}}}T}{\epsilon }-\frac{\epsilon \zeta \Omega _{0}^{2}{\rm sin} (2\phi )}{\Gamma \Omega }\end{eqnarray} \tag{ 33 }$$

acting on the variable . As a result, the system samples the -distribution

$$\begin{eqnarray}&&{{P}_{\epsilon }}(\epsilon )\propto \epsilon {\rm exp} \left\{ -\beta \left[ \left( 1+\frac{\zeta \Omega _{0}^{2}{\rm sin} (2\phi )}{2\Gamma \Omega } \right){{\epsilon }^{2}}+\frac{\eta {{\Omega }_{0}}}{4\;m\Gamma {{\Omega }^{2}}}{{\epsilon }^{4}} \right] \right\}.\end{eqnarray} \tag{ 34 }$$

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,

$$\begin{eqnarray}&&{{P}_{E}}(E)=\frac{1}{Z}{\rm exp} \left\{ -\beta \left[ \left( 1+\frac{\zeta \Omega _{0}^{2}{\rm sin} (2\phi )}{2\Gamma \Omega } \right)E+\frac{\eta {{\Omega }_{0}}}{4\;m\Gamma {{\Omega }^{2}}}{{E}^{2}} \right] \right\}.\end{eqnarray} \tag{ 35 }$$

The normalization factor $Z=\int {{P}_{E}}(E){\rm d}E$ is given by

$$\begin{eqnarray}&&Z=\sqrt{\frac{\pi m\Gamma {{\Omega }^{2}}}{\beta \eta {{\Omega }_{0}}}}h\left( \sqrt{\frac{\beta m\Gamma {{\Omega }^{2}}}{\eta {{\Omega }_{0}}}}\left( 1+\frac{\zeta \Omega _{0}^{2}{\rm sin} (2\phi )}{2\Gamma \Omega } \right) \right),\end{eqnarray} \tag{ 36 }$$

where the function h(x) is defined as

$$\begin{eqnarray}&&h(x)={\rm exp} ({{x}^{2}}){\rm erfc}(x)\end{eqnarray} \tag{ 37 }$$

and ${\rm erfc}(x)$ is the complementary error function. Thus, the energy distribution is that of an equilibrium system with effective energy

$$\begin{eqnarray}&&H=\left[ 1+\frac{\zeta \Omega _{0}^{2}{\rm sin} (2\phi )}{2\Gamma \Omega } \right]E+\frac{\eta {{\Omega }_{0}}}{4\;m\Gamma {{\Omega }^{2}}}{{E}^{2}}\end{eqnarray} \tag{ 38 }$$

and configurational partition function Z. While the term proportional to E² is caused by the feedback cooling, the term proportional to E is affected only by the parametric modulation.

According to equation (35), the energy distribution is Gaussian with a cutoff at E = 0. The maximum of the Gaussian is located at

$$\begin{eqnarray}&&\bar{E}=-\frac{2\;m\Gamma {{\Omega }^{2}}}{\eta {{\Omega }_{0}}}\left[ 1+\frac{\zeta \Omega _{0}^{2}{\rm sin} (2\phi )}{2\Gamma \Omega } \right]\end{eqnarray} \tag{ 39 }$$

while its variance (neglecting the cutoff) is given by

$$\begin{eqnarray}&&\sigma _{E}^{2}=\frac{2\;m\Gamma {{\Omega }^{2}}{{k}_{{\rm B}}}T}{\eta {{\Omega }_{0}}}.\end{eqnarray} \tag{ 40 }$$

Hence, the width of the Gaussian does neither depend on the driving parameters nor on the phase ϕ.

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}_{{\rm mc}}}(q,p;\tilde{E})$ of the oscillator evolving at a given constant total energy $\tilde{E}$,

$$\begin{eqnarray}&&{{P}_{{\rm mc}}}(q,p;\tilde{E})=\frac{1}{g(\tilde{E})}\delta [E(q,p)-\tilde{E}],\end{eqnarray} \tag{ 41 }$$

where $\delta (x)$ is the Dirac delta function and we have denoted the fixed value of the energy with $\tilde{E}$ to distinguish it from the energy function $E(q,p)$, which depends on the position q and the momentum p. The normalizing factor $g(\tilde{E})$ is the micro-canonical density of states,

$$\begin{eqnarray}&&g(\tilde{E})=\int {\rm d}q{\rm d}p\ \delta [E(q,p)-\tilde{E}].\end{eqnarray} \tag{ 42 }$$

The phase space distribution of equation (41) is that of 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,

$$\begin{eqnarray}&&{{P}_{qp}}(q,p)=\int {\rm d}\tilde{E}{{P}_{E}}(\tilde{E}){{P}_{{\rm mc}}}(q,p;\tilde{E})=\int {\rm d}\tilde{E}\frac{{{P}_{E}}(\tilde{E})}{g(\tilde{E})}\delta [E(q,p)-\tilde{E}].\end{eqnarray} \tag{ 43 }$$

This linear superposition of micro canonical distributions is valid as long as the energy changes slowly on the time scale of one oscillation period. For the low friction constants and the small feedback strength studied here this assumption is met even under non-equilibrium conditions. Carrying out the integral yields

$$\begin{eqnarray}&&{{P}_{qp}}(q,p)=\frac{{{P}_{E}}[E(q,p)]}{g[E(q,p)]}.\end{eqnarray} \tag{ 44 }$$

As a further approximation, we now use the density of states $g(E)=2\pi /{{\Omega }_{0}}$ for the harmonic oscillator, thereby neglecting the Duffing term of the potential in this part of the calculation, and obtain

$$\begin{eqnarray}&&{{P}_{qp}}(q,p)=\frac{{{\Omega }_{0}}}{2\pi }{{P}_{E}}[E(q,p)].\end{eqnarray} \tag{ 45 }$$

Inserting the energy distribution from equation (35) into this equation we finally find the phase space distribution function

$$\begin{eqnarray}&&{{P}_{qp}}(q,p)=\frac{{{\Omega }_{0}}}{2\pi Z}{\rm exp} \left\{ -\beta \left[ \left( 1+\frac{\zeta \Omega _{0}^{2}{\rm sin} (2\phi )}{2\Gamma \Omega } \right)E(q,p)+\frac{\eta {{\Omega }_{0}}}{4\;m\Gamma {{\Omega }^{2}}}E{{(q,p)}^{2}} \right] \right\}.\end{eqnarray} \tag{ 46 }$$

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,

$$\begin{eqnarray}&&{{P}_{q}}(q)=\int _{-\infty }^{\infty }{\rm d}p\;{{P}_{qp}}(q,p).\end{eqnarray} \tag{ 47 }$$

In the absence of parametric modulation ($\zeta =0$), one finds by carrying out the integral

$$\begin{eqnarray}\begin{array}{rcl} {{P}_{q}}(q) & \propto & \sqrt{2+\eta \frac{{{\Omega }_{0}}}{\Gamma }\left( \frac{{{q}^{2}}}{2}+\xi \frac{{{q}^{4}}}{4} \right)}{\rm exp} \left( -\frac{\beta m{{\Omega }_{0}}}{8\eta }{{\left[ 2+\eta \frac{{{\Omega }_{0}}}{\Gamma }\left( \frac{{{q}^{2}}}{2}+\xi \frac{{{q}^{4}}}{4} \right) \right]}^{2}} \right) \\ {} & {} & \times {{K}_{\frac{1}{4}}}\left( \frac{\beta m{{\Omega }_{0}}}{8\eta }{{\left[ 2+\eta \frac{{{\Omega }_{0}}}{\Gamma }\left( \frac{{{q}^{2}}}{2}+\xi \frac{{{q}^{4}}}{4} \right) \right]}^{2}} \right), \\ \end{array}\end{eqnarray} \tag{ 48 }$$

where ${{K}_{1/4}}$ is a generalized Bessel function of the second kind. For simplicity, we have considered the case ${{\Omega }_{0}}=\Omega$ here. A similar expression can also be derived for the momentum distribution.

As shown recently, a fluctuation theorem holds for the relative entropy change $\Delta \mathcal{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 $\Delta \mathcal{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)$,

$$\begin{eqnarray}&&\Delta \mathcal{S}={\rm ln} \frac{P[u(t)]}{P[{{u}^{*}}(t)]}.\end{eqnarray} \tag{ 49 }$$

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, a detailed fluctuation can be proven for the quantity $\Delta \mathcal{S}$,

$$\begin{eqnarray}&&{{P}_{t}}(-\Delta \mathcal{S})/{{P}_{t}}(\Delta \mathcal{S})={\rm exp} (-\Delta \mathcal{S}),\end{eqnarray} \tag{ 50 }$$

where ${{P}_{t}}(\Delta \mathcal{S})$ is the probability density to observe the value $\Delta \mathcal{S}$ at time t as determined over many repetitions of the relaxation experiment. For the relaxation process considered in [15] the relative entropy change is given by

$$\begin{eqnarray}&&\Delta \mathcal{S}=\beta {{Q}_{h}}+\Delta \phi ,\end{eqnarray} \tag{ 51 }$$

where ${{Q}_{h}}=-[{{E}_{t}}-{{E}_{0}}]$ is the energy absorbed by the bath during the relaxation, and E₀ and E_t are the energy of the oscillator at time 0 and t, respectively. The quantity $\phi (q,p)$ is defined as the logarithm of the stationary phase space distribution

$$\begin{eqnarray}&&\phi (q,p)=-{\rm ln} {{P}_{qp}}(q,p)\end{eqnarray} \tag{ 52 }$$

and $\Delta \phi$ is the difference of ϕ at the beginning and the end of the trajectory,

$$\begin{eqnarray}&&\Delta \phi ={{\phi }_{t}}-{{\phi }_{0}}.\end{eqnarray} \tag{ 53 }$$

Hence, the relative entropy change $\Delta \mathcal{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 $\Delta \phi$ is unknown. However, from the distribution derived for our model, equation (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

$$\begin{eqnarray}&&\Delta \mathcal{S}=-\beta \frac{\zeta \Omega _{0}^{2}{\rm sin} (2\phi )}{2\Gamma \Omega }\left[ {{E}_{t}}-{{E}_{0}} \right]-\beta \frac{\eta {{\Omega }_{0}}}{4\;m\Gamma {{\Omega }^{2}}}\left[ E_{t}^{2}-E_{0}^{2} \right],\end{eqnarray} \tag{ 54 }$$

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 the relaxation resembles that of an oscillator initially coupled to a thermal bath with effective temperature ${{T}_{{\rm eff}}}=T/(1-\frac{\zeta \Omega _{0}^{2}{\rm sin} (2\phi )}{2\Gamma \Omega })$. By choosing parameters appropriately, 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 η.

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 analyze its motion in terms of the so-called quadratures. For the oscillator driven by the parametric modulation ${{F}_{{\rm drive}}}=\zeta m\Omega _{0}^{2}{\rm cos} ({{\Omega }_{m}}t)q$, we write the time evolution of the oscillator position as

$$\begin{eqnarray}&&q(t)=R(t){\rm cos} [\Omega t+\phi (t)],\end{eqnarray} \tag{ 55 }$$

where Ω is the frequency of the particle oscillating at half the frequency of the driving, $\Omega ={{\Omega }_{m}}/2$. Here, R(t) and $\phi (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, equation (55) can be written as the sum of two contributions, one in-phase with the driving signal and one out-of-phase,

$$\begin{eqnarray}\begin{array}{rcl} q(t) & = & R(t){\rm cos} \phi (t){\rm cos} (\Omega t)-R(t){\rm sin} \phi (t){\rm sin} (\Omega t) \\ {} & = & X(t){\rm cos} (\Omega t)-Y(t){\rm sin} (\Omega t), \\ \end{array}\end{eqnarray} \tag{ 56 }$$

where the second line defines the in-phase component $X(t)=R(t){\rm cos} \phi (t)$ and the quadrature $Y(t)=R(t){\rm sin} \phi (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:

$$\begin{eqnarray}\begin{array}{rcl} p(t)/m\Omega & = & -R(t){\rm sin} [\Omega t+\phi (t)] \\ {} & = & -X(t){\rm sin} (\Omega t)-Y(t){\rm cos} (\Omega t), \\ \end{array}\end{eqnarray} \tag{ 57 }$$

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 equation (56) yields

$$\begin{eqnarray}\begin{array}{rcl} X(t) & = & \;\;\;\;q(t){\rm cos} (\Omega t)-\frac{p(t)}{m\Omega }{\rm sin} (\Omega t), \\ Y(t) & = & -q(t){\rm sin} (\Omega t)-\frac{p(t)}{m\Omega }{\rm cos} (\Omega t). \\ \end{array}\end{eqnarray} \tag{ 58 }$$

corresponding to transformation to a coordinate system that rotates clockwise with frequency Ω with respect to the $(q,p/m\Omega )$-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

$$\begin{eqnarray}\begin{array}{rcl} R & = & \sqrt{{{X}^{2}}+{{Y}^{2}}}, \\ \phi & = & {\rm arctan} (Y/X), \\ \end{array}\end{eqnarray} \tag{ 59 }$$

and that

$$\begin{eqnarray}&&\frac{m{{\Omega }^{2}}{{R}^{2}}}{2}=\frac{m{{\Omega }^{2}}}{2}({{X}^{2}}+{{Y}^{2}})\end{eqnarray} \tag{ 60 }$$

is the energy of a harmonic oscillator with frequency Ω.

In our simulations, we integrated the Langevin equation of motion with the OVRVO algorithm of Sivak et al [28], which can be viewed as a stochastic generalization 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 $\Delta 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 ($\Delta 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\times {{10}^{6}}$ modulation cycles. For some parameters we carried out longer simulations of up to $3\times {{10}^{10}}$ steps corresponding to a total simulation time of $t=3\times {{10}^{8}}$. All simulations were carried out for ${{k}_{{\rm 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 $\mathcal{E}={{k}_{{\rm B}}}T$, the inverse frequency $\mathcal{T}=1/{{\Omega }_{0}}$ and the particle mass $\mathcal{M}=m$ as our basic units of energy, time and mass, respectively. Accordingly, distances are measured in units of $\mathcal{L}=(1/{{\Omega }_{0}})\sqrt{{{k}_{{\rm B}}}T/m}$ and velocities in units of $\mathcal{V}=\sqrt{{{k}_{{\rm B}}}T/m}$. Hence, the unit of length is given by the variance of the position of the harmonic oscillator, $\langle {{q}^{2}}\rangle ={{k}_{{\rm B}}}T/m\Omega _{0}^{2}={{\mathcal{L}}^{2}}$ and the unit of energy is the average energy of the harmonic oscillator $\langle E\rangle ={{k}_{{\rm B}}}T=\mathcal{E}$. The friction constant is given in units of ${{\Omega }_{0}}$ such that it equals the inverse of the quality factor, $Q={{\Omega }_{0}}/\Gamma =1/\Gamma \mathcal{T}$. The feedback strength η and the Duffing coefficient ξ have the dimension of $1/$ area and are measured in units of $1/{{\mathcal{L}}^{2}}$. The modulation depth ζ is dimensionless. In the following, we use reduced units in which $\mathcal{E}=\mathcal{T}=\mathcal{M}=1$.

We first consider the oscillator without parametric modulation ($\zeta =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,

$$\begin{eqnarray}&&{{P}_{E}}(E)\propto {\rm exp} \left\{ -\beta \left( E+\frac{\eta }{4\;m\Gamma {{\Omega }_{0}}}{{E}^{2}} \right) \right\},\end{eqnarray} \tag{ 61 }$$

where we have assumed that the particle oscillates with $\Omega ={{\Omega }_{0}}$. The first term in the exponential is the same as that of the uncooled oscillator, but the second term proportional to E² is due to the feedback loop and strongly penalizes high energy states thereby cooling the system. The cooling effect is stronger for weak friction Γ and small frequencies ${{\Omega }_{0}}$. Several energy distributions obtained from simulations together with the corresponding predictions of equation (61) are shown in the left panel of figure 1. The simulations were carried out for a friction of $\Gamma =0.0001$ and a Duffing parameter of $\xi =-0.022$. Without feedback, $\eta =0$, the energy distribution is exponential, but for $\eta \gt 0$ the E² 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 figure 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.

Zoom In Zoom Out Reset image size

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 equation (35) is expected to be exponential

$$\begin{eqnarray}&&{{P}_{E}}(E)\propto {\rm exp} \left\{ -\beta \left( 1+\frac{\zeta {{\Omega }_{0}}{\rm sin} (2\phi )}{2\Gamma } \right)E \right\},\end{eqnarray} \tag{ 62 }$$

where we have assumed that the modulation frequency is ${{\Omega }_{m}}=2{{\Omega }_{0}}$. For a vanishing Duffing parameter $\xi =0.0$, i.e., for a perfectly harmonic trap, the phase is expected to be $\phi =-\pi /4$ in the absence of thermal fluctuations [31]. If this is the case, the decay constant of the exponential is $\beta (1-\zeta {{\Omega }_{0}}/2\Gamma )$. Hence, the decay constant is positive only for $\zeta \lt 2\Gamma /{{\Omega }_{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 $\zeta \gt 2\Gamma /{{\Omega }_{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 equation (35). Several energy distributions for this case are shown in figure 2. Note that we performed these calculations for a relatively large friction constant of $\Gamma =0.01$, because for lower friction it takes exceedingly long to sample all relevant energies. For weak driving, $\xi =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 figure 2. The theory developed above, on the other hand, assumes a fixed phase of $\phi =-\pi /4$ (for $\xi =0$). For $\xi =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 $\phi =-\pi /4$ and one at $\phi =-\pi /4+\pi$. Since the phase relation is more pronounced at high energies, in this regime the energy distributions shown in the left panel of figure 2 converge to the form predicted by theory. In the figure, the theoretical distributions are indicated by lines with logarithmic slope of $-\beta (1-\zeta {{\Omega }_{0}}/2\Gamma )$. 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.

Zoom In Zoom Out Reset image size

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, equation (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 figure 3. These simulations were carried out for a modulation depth of $\zeta =0.03$ and and a feedback strength of $\eta =0.022$, because these values can be realized 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, $\xi =0.0$, the peak is located at $\phi =-\pi /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.

Zoom In Zoom Out Reset image size

A closer analysis of how the phase depends on the Duffing parameter is shown in figure 4. The left panel of the figure shows the positions of the maximum of the phase distribution. i.e., the most likely phase ${{\phi }_{{\rm max} }}$, as a function of the Duffing parameter ξ for different friction constants Γ. As can be inferred from the figure, the most likely phase ${{\phi }_{{\rm 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.

Zoom In Zoom Out Reset image size

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\pi .$ 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 figure 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 equation (35). To do that, we identify the phase ϕ occurring in the theoretical expression with the most likely phase ${{\phi }_{{\rm max} }}$ determined in the simulations. Energy distributions obtained for friction constants ranging from $\Gamma ={{10}^{-5}}$ to $\Gamma ={{10}^{-3}}$ are shown in the left panel of figure 5. In all cases, the system was driven at ${{\Omega }_{m}}=2{{\Omega }_{0}}$ and the Duffing parameter, the feedback strength and the modulation depth were $\xi =-0.022$, $\eta =0.022$, $\zeta =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 figure 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.

Zoom In Zoom Out Reset image size

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 figure 6. From left to right, the driving frequency ${{\Omega }_{m}}$ is slightly below $2{{\Omega }_{0}}$, equal to $2{{\Omega }_{0}}$ and slightly above $2{{\Omega }_{0}}$. As in previous simulations, the parameters were $\xi =-0.022$, $\eta =0.022$, and $\zeta =0.03$. At ${{\Omega }_{m}}=2\Omega$ 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{{\Omega }_{0}}$, the distributions are deformed, indicating the occurrence of classical squeezing.

Zoom In Zoom Out Reset image size

In our experiments we use a silica nanoparticle trapped at the focus of single beam optical tweezers. The optical tweezer is formed by a $1064\;{\rm nm}$ laser beam ($\sim 35\;{\rm mW}$) focused by a ${\rm NA}=0.9$ objective, which is mounted inside a vacuum chamber. The particle motion is recorded with an additional colinear laser ($780\;{\rm nm}$) and three balanced photodetectors. A home-built electronic circuit is used to generate the feedback signal (η), while a frequency generator provides 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.1\;{\rm mBar}$. 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 [33, 34]. Subsequent measurements are carried out at $1.2\times {{10}^{-5}}\;{\rm 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, for small oscillation amplitudes, 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.

The particle oscillation frequencies along the three main axes are well separated and don't 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]

$$\begin{eqnarray}&&{{R}_{{\rm ML}}}=|{{A}_{q}}(\Omega )|,\end{eqnarray} \tag{ 63 }$$

$$\begin{eqnarray}&&{{\phi }_{{\rm ML}}}={\rm arg} \left[ {\rm exp} (-{\rm i}{{\Omega }_{0}}{{t}_{0}}){{A}_{q}}(\Omega ) \right],\end{eqnarray} \tag{ 64 }$$

where $\Omega /2\pi$ is the estimated frequency of the signal, t₀ is the time origin and

$$\begin{eqnarray}&&{{A}_{q}}(\omega )=\frac{1}{N}\mathop{\sum }\limits_{n=0}^{N-1}{{q}_{n}}{\rm exp} (-{\rm i}\omega n\Delta t).\end{eqnarray} \tag{ 65 }$$

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\Delta t$, N is the number of samples entering the estimation and $\Delta t$ is the sampling interval. The estimation of the amplitude and the phase relies on precise estimation of the frequency Ω. We estimate Ω by maximizing (65) with respect to ω, i.e. $A(\Omega )={\rm max} (A(\omega ))$. The width of the function $A(\omega )$, and thereby our ability to localize the maximum, depends on the length of the time trace q(t). Therefore, we use a long time trace measured over ${{T}_{{\rm meas}.}}=0.1\;{\rm s}$ and sampled at 625 kilosamples/second to estimate Ω. Subsequently, we use that value of Ω and equations (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 maximizing (65) allows us to estimate the frequency with much better accuracy than $1/{{T}_{{\rm 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 to the first oscillator. Formally, this is expressed as

$$\begin{eqnarray}&&\Delta \phi ={\rm arg} \left[ {{A}_{p}}\cdot {{\left[ A_{m}^{\frac{{{\Omega }_{p}}}{{{\Omega }_{m}}}} \right]}^{*}} \right]=({{\phi }_{p}}-\frac{{{\Omega }_{p}}}{{{\Omega }_{m}}}[{{\phi }_{m}}+k2\pi ]),\end{eqnarray} \tag{ 66 }$$

where A_p and A_m are the Fourier transforms of the two signals, respectively (cf (65)), and k is an integer which takes into account that the phase is only determined up to modulo $2\pi$. Note that the exponent ${{\Omega }_{p}}/{{\Omega }_{m}}$ takes care that (66) does not depend on t₀. Without loss of generality, we set ${{\phi }_{m}}=0$, i.e. we choose our time origin such that it coincides with a maximum of the signal with frequency ${{\Omega }_{m}}$. For the special case of a parametrically driven particle, which oscillates at half the frequency of the parametric modulation (${{\Omega }_{m}}=2{{\Omega }_{p}}$), we get $\Delta \phi ={{\phi }_{p}}-k\pi$. Therefore, the above method allows to estimate the relative phase between the particle oscillation and the parametric modulation up to a multiple of π.

We measure the distribution of the energy and phase for modulation at ${{\Omega }_{m}}/2\pi =247,248,249$, and $250\;{\rm kHz}$. Each distribution is obtained from 100 time traces of $0.1\;{\rm s}$ duration. Figure 7 shows the maximum values of the energy and phase distributions shown in figure 8 and a fit to secular perturbation theory [31, 34]. 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 ${{\phi }_{0}}=50{}^\circ$ 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 $\xi =-5.4\pm 1.1\;\mu {{{\rm m}}^{-2}}$, $\eta =3.9\pm 1.3\;\mu {{{\rm m}}^{-2}}$ and $\zeta =16.1\pm 5.7\times {{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 summarized in table 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 ($\zeta =0$). A fit of the energy distribution to equation (61) yields $\eta =4.5\pm 0.9\;\mu {{{\rm m}}^{-2}}$, in good agreement with the previously determined value.

Figure 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. Figure 9 shows the widths of the distributions obtained from the Gaussian fits in figure 8 and from numerical simulations. For comparison, we also show the theoretical prediction according to equation (40). The broadening of the distributions has two contributions, thermal motion and parameter fluctuations.

Zoom In Zoom Out Reset image size

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 (cf (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].

In addition to Brownian motion, parameter fluctuations broaden the experimental distributions [37]. The experimental parameters fluctuate due to laser intensity and polarization 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 ${{\sigma }_{E}}=78\pm 3\times {{10}^{-3}}$ k_BT and ${{\sigma }_{\phi }}=1.7\pm 0.1\times {{10}^{-3}}\pi$, respectively, are approximately one order of magnitude larger than the theoretical values $5.1\times {{10}^{-3}}$ k_BT and $0.15\times {{10}^{-3}}\pi$, 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.