Higher order correlations in a levitated nanoparticle phonon laser Higher order correlations in a levitated nanoparticle phonon laser

: We present theoretical and experimental investigations of higher order correlations of mechanical motion in the recently demonstrated optical tweezer phonon laser, consisting of a silica nanosphere trapped in vacuum by a tightly focused optical beam [R. M. Pettit et al., Nature Photonics 13 , 402 (2019)]. The nanoparticle phonon number probability distribution is modeled with the master equation formalism in order to study its evolution across the lasing threshold. Up to fourth-order equal-time correlation functions are then derived from the probability distribution. Subsequently, the master equation is transformed into a nonlinear quantum Langevin equation for the trapped particle’s position. This equation yields the non-equal-time correlations, also up to fourth order. Finally, we present experimental measurements of the phononic correlation functions, which are in good agreement with our theoretical predictions. We also compare the experimental data to existing analytical Ginzburg-Landau theory where we ﬁnd only a partial match.

Recently, we have proposed and demonstrated a phonon laser based on a nanosphere levitated in an optical tweezer by applying linear and nonlinear feedback on a center-of-mass mode of oscillation [27]. In the context of optomechanics, phonon lasers have been previously studied in a wide range of different systems as mechanical analogs to optical lasers [28][29][30][31][32][33]. The uses of the phonon laser demonstrated earlier by us are expected to be practical (generating nonclassical mechanical states starting with coherent phonons) as well as conceptual (exploring analogies to the optical laser). Following the latter direction, we recall that in the optical case, correlation functions of order higher than two [34] form the basis for correlation interferometry [35], nonclassical state characterization [36], determining higher order coherence properties [37] and photon-photon interactions [38]. In this article, we consider the higher order correlation functions of the phonon number distribution to deepen the optics-mechanics analogy, to quantify the achievable degree of coherence and the relevant timescales of our phonon laser, and to further test the agreement between theoretical predictions and experiment.
Specifically, we calculate and measure the center-of-mass phonon number probability distribution for an optically trapped silica nanosphere under the influence of both linear feedback amplification and nonlinear feedback cooling. The phonon number probability distribution is modeled and measured in the steady state as the phonon laser transitions across the threshold. The higher order correlation functions of the phonon number operator are then considered up to fourth order to quantify the degree of coherence throughout the transition. Finally, by transforming the master equation model into a quantum Langevin equation, the non-equal-time higher order correlations are obtained up to fourth order and compared with experimental data as well. Good agreement is found between experiment and theory.

Experiment
Realization of our phonon laser is based on the free-space optical levitation of a silica nanoparticle dipole-trapped in vacuum as shown in Fig. 1 (see [27] for details of the experimental apparatus). Light scattering from the trapped particle provides a position measurement that is processed electronically to generate feedback signals that control the particle's dynamics. These signals provide nonlinear parametric cooling and linear amplification of center-of-mass phonons, giving rise to the realization of the phonon laser. The former hinders the particle's motion by increasing the trap stiffness whenever the particle moves away from the trap center and reducing it when the particle falls back toward the trap. This is equivalent to the modulation of the trap power by supplying feedback signals ∝ r i (t) r i (t) in all three positional degrees of freedom where Ω i is the frequency of oscillation, this yields a power modulation at twice the particle's natural oscillation frequency with an appropriate phase shift.
Our feedback cooling electronics are composed of a series of homebuilt analog circuits that take the measured position signals at different directions and generate the required phase-shifted and frequency-doubled signals. These signals are combined in a summing amplifier circuit and fed back to an electro-optic modulator that controls the trap power. The magnitude of this signal is adjusted using a feedback gain controller on the same circuit. The linear feedback amplification of the particle's oscillation is addressed by modulating the trap power about its mean value by a signal ∝ − r i (t) that acts at the frequency of particle oscillation Ω i . This signal is derived from a phase-locked loop (PLL) on a digital lock-in amplifier (Zurich Instruments HF2LI). The PLL tracks the phase and frequency of the particle's motion along a desired axis of the trap, and from the PLL a signal of the same frequency is generated with a fixed phase relationship to the measured particle motion. The phase of the generated signal is set to induce heating of the particle's motion, and the final output amplitude can be set to the desired value. Deriving the amplification signal from a PLL produces a free-running oscillator that is not locked to any fixed external phase reference. The amplification signal only addresses a single spatial degree of freedom of the particle's motion (x-axis here).

Theory
The phonon laser can be modeled by the master equation [27] [27]. The schematic includes an electro-optic modulator (EOM), polarizing beam splitters (PBS1 and PBS2), high numerical aperture focusing lenses (L1 and L2), a detector (D), two switches (S 1 and S 2 ) and a beam dump (BD). A nanoparticle is optically trapped by the single-beam optical tweezer, and its position measured using a probe beam. The particle is addressed with nonlinear feedback cooling (γ c ) and linear feedback amplification (γ a ); associated with this measurement and feedback are also backactions (Γ c and Γ a ) in this system. For the feedback cooling and heating, BP, 2Ω 0 , ∆Φ, G h(c) and PLL denote the electronic bandpass filter, frequency doubler, phase shifter, electronic gain for heating (cooling), and phase-locked loop, respectively [27].
where ρ is the density matrix operator which describes the center-of-mass motion of the nanoparticle along the x -axis, Ω 0 is the natural oscillation frequency of the nanoparticle in the optical harmonic trap, b(b † ) is the annihilation (creation) operator of phonon number, and are the dimensionless position and momentum operators, respectively, along the x axis. The quantities A t , D x , and D p are the rates associated with photon scattering from the trap beam, and position and momentum diffusion, respectively, such that where N 0 = k B T/ Ω 0 is the initial thermal phonon number at temperature T, k B being Boltzmann's constant, m and V are the nanoparticle mass and volume, respectively, c is the dielectric permittivity of the nanoparticle, c is the speed of light, and ω 0 and I t are the optical frequency and intensity, respectively. The parameters γ g , γ a , γ c , Γ a and Γ c denote the rates of gas damping, linear feedback amplification, nonlinear feedback cooling and their respective backaction rates. The experiments are conducted in a regime where feedback backaction is negligible and we therefore do not provide the detailed expressions for Γ a and Γ c . Also, where M a (M c ) is the modulation depth of the amplifying (cooling) feedback signal and N is the mean phonon occupation number. The gas damping is determined by residual air pressure [27]. The Lindblad superoperator D[O] for the operator O acts on ρ according to By utilizing the rotating frame transformation, the dynamical equation for the probability that the nanoparticle center of mass is excited to n phonons, P(n) = n| ρ|n , can be derived from the master equation Eq. (1) as [27] where D t = A t + D p . The various terms in Eq. (4) correspond to one-and two-phonon transitions between the oscillator's levels with n − 1, n and n + 1 phonons, respectively; the contributions due to stimulated emission and spontaneous emission, the steady-state solution and time-dependent phonon number distribution P(n) have been provided in [27].
However, for our present purpose, the modulation-evolution (i.e. the evolution of P(n) with M a ) of the steady state phonon number distribution is also important. In Fig. 2, the numerical calculation, from Eq. (4), of the evolution of the phonon number probability distribution P(n) as a function of M a is shown. Far below threshold, the phonon number distribution satisfies a Maxwell-Boltzmann (M-B) distribution; far above the threshold, the distribution is close to Poissonian. These results expand on those presented previously in [27]. The evolution of the steady state phonon number probability distribution P(n) from below to above threshold as the modulation depth M a of the linear feedback amplification is increased. The points correspond to experimental data and the curves to numerical results. Different colors indicate phonon probability distributions for varying modulation depth. The system parameters are particle radius R = 69.1 nm, Ω 0 = 2π × 128.5 kHz, and γ c /(2π) = 4.4 × 10 −4 Hz; the background pressure is 6 × 10 −5 mbar and T = 300K.

Higher order correlations in phonon laser
We now consider the correlation functions of the optical tweezer phonon laser. Correlation functions in the optical case have been found to be useful for the measurement of quantum coherence, imaging, and other applications [39,40]. The second-, third-, and fourth-order correlation functions can be expressed in terms of the creation (b † ) and annihilation (b) operators as follows [41], where indicates an ensemble average and τ, τ 1 , τ 2 and τ 3 represent time delays. From the definition of higher order correlations, when the time delay τ = 0, we can obtain equal-time correlations as follows [42] Thus, g (n) (t, 0) is the equal-time nth order correlation at time t; for notational convenience we express it below as g (n) 0 (t). For describing experimental and theoretical applications more clearly, the definitions of higher order correlations have been clarified in Appendix A.

Equal-time higher order correlation functions in the steady state
Solving Eqs. (4) and (6) numerically, we obtained the equal-time higher order correlation functions in steady state as well as in transience. The variation of the steady state phonon correlations as a function of the feedback amplification modulation M a is shown in Fig. 3(a). The corresponding experimental results obtained from Eq. (14) with all time delays equal zero are shown in the same figure with M a = δP a /P 0 , where P 0 is the power of the trapping beam and δP a is the power modulation induced by the feedback amplification. Phonon occupation numbers used in these calculations were recorded by monitoring oscillator dynamics in a time window of 20ms and the length of the error bars represents ±1 standard deviation of 100 such measurements. When the modulation is small, g (n) 0 (∞) → n!, as expected for a thermal state, and for large modulation g (n) 0 (∞) → 1, as expected for a coherent state [35]. The theoretical predictions are in good agreement with the experimental data.
In contrast, Fig. 3(b) shows the transient evolution of the system after switching on the amplification modulation M a = 5 × 10 −3 , which takes the laser far above threshold. In this case the equal-time higher order correlations can be recorded in real time. Experimentally, these variations were constructed using the transient mean phonon population calculated from 500 iterations of the switching experiment in which the gain was switched on at time t = 0. After each period of amplification, the gain was switched off and the particle was allowed to cool under the influence of the feedback cooling to re-initialize its thermal state. When t < 0, the phonon number of the particle satisfies the thermal state distribution. After switching on the feedback amplification, the probability distribution of phonon number approaches a coherent state distribution gradually. Thus, when the laser is far below threshold, the higher order correlations will start at n!, but tend to unity at long times. Figures 3(c) and (d) show the corresponding Fano factor F = ( N 2 − N 2 )/ N , and Figs. 3(e) and (f) the Signal-to-Noise Ratio SNR = N /( N 2 − N 2 ) 1/2 . As can be seen, the theoretical predictions agree well with experimental data.

Non-equal-time higher order correlation functions
For obtaining non-equal-time higher order correlations, we transform the master equation of Eq.
(1) into a quantum Langevin equation [43,44], where q = /(2mΩ 0 )Q x is the position of the particle, ξ T is the Brownian noise from the environment, and ξ Fa and ξ Fc are the noises from the feedback amplification and cooling respectively, with zero means and correlations, i.e. ξ i (t)ξ j (t ) = δ i,j δ(t − t ), i, j = T, Fa and Fc. We extract the non-equal-time higher order correlations numerically by solving Eq. (7) for q and using that to find N, as indicated in Appendix A.
We first present results with the feedback cooling and heating turned off. In this case, the phonon statistics satisfy the M-B distribution. The quantities g (2) (τ), g (3) (τ = (τ 2 1 + τ 2 2 ) 1/2 ) and g (4) are obtained by numerically simulating Eq. (7), and as can be seen from Fig. 4, compare well with experiment. When the time delay τ equals zero, the nth order correlation equals n!, which corresponds to the case of equal-time higher order correlation. From the Eq. (7) and Fig. 4, we can see that the noise spectrum is Lorentzian in profile. The corresponding second order correlation is given by [41], where τ c is the coherence time of the thermal state below threshold. From the experimental and numerical results, we find τ c = 40µs. Considering the third order correlation in more detail in the (τ 1 , τ 2 ) plane, the experimental and numerical results are shown in Fig. 5. Mathematically the ridges in the plots can be understood from the symmetries of g (3) (t, τ 1 , τ 2 ) in Eq. (5): they lie on the lines τ 1 = 0, τ 2 = 0, and τ 1 = τ 2 . Physically, the center of the plot corresponds to three-phonon bunching, while the ridges correspond to two-phonon bunching. In the figure, we see quite good agreement between theory and experiment. Now we consider the case for g (2) when feedback cooling and amplification are turned on. In this case the laser is in a coherent state to a good approximation and g (2) (τ) ∼ 1 and the second-order correlation does not yield the coherence time. However, it can be found by inverting the linewidth, which yields the coherence time τ c ∼ 1s above threshold [27].
Moving on to the higher-order correlations, g (3) (τ 1 , τ 2 ) is shown below, near and above threshold in Fig. 6. In Fig. 6, the figures in the upper row are experimental results, and the figures in the lower row are theoretical calculations. When the feedback amplification is turned off, the phonon number distribution satisfies the M-B distribution. When the time delays equal zero, g (3) (τ 1 , τ 2 ) has the maximum value g (3) (τ 1 = 0, τ 2 = 0) = 6. This is in agreement with the results in Fig. 3. When the time delays do not all equal zero, there are three ridges. They intersect the origin along the line τ 1 = 0, τ 2 = 0 and τ 1 = τ 2 . When the modulation increases, there is still one maximum at the origin and three ridges, but g (3) (τ 1 , τ 2 ) decreases until the phonon number distribution approaches Poissonian, i.e. such that g (3) Similarly, we consider the fourth order correlation below, near and above the threshold, and the corresponding properties of g (4) (τ 1 , τ 2 , −τ 1 ) have been shown in Fig. 7. The peak value of g (4) occurs at the point (τ 1 = 0, τ 2 = 0, τ 3 = 0). Different amplification modulation factors lead to different maxima of g (4) with larger M a leading to smaller maxima. When the modulation leads the laser to far above threshold, the fourth order correlation is unity. There are also ridges like in Fig. 6, which can be deduced from the symmetries of g (4) (τ 1 , τ 2 , −τ 1 ) in Eq. (5). Physically the center of the graph corresponds to four-phonon bunching and the ridges to three-phonon bunching.
We explore further the case for g (4) when the phonon laser works below threshold and τ 3 is not always equal to −τ 1 . The corresponding τ 3 evolution of the fourth order correlations is shown in Fig. 8. In this case, there is always a maximum value at the origin, which occurs when τ 1 = τ 2 = τ 3 = 0, g (4) reaches its peak value of 24.

Ginzburg-Landau theory for a phonon laser
In this section we consider the use of higher-order correlation functions to determine how well a standard paradigm of second-order phase transitions, namely the Ginzburg-Landau (G-L) model, can describe our phonon laser [46][47][48]. G-L is an elegant, powerful and useful theory, allowing, for example, for the analytical calculation of correlation functions. It is therefore of interest to see how well these analytic results from a standard theory can model our experimental correlations. A compact method for displaying the range of agreement, independent of detector quantum efficiency, has been suggested recently; this method involves comparing the higher order correlations, rather than the full phonon statistics, and we follow those authors [48].
By utilizing the transformation ρ(t) = ∫ φ|v v|d 2 v, the master equation of Eq. (1) was transformed into the Glauber-Sudarshan P representation and the Fokker-Planck equation for the P-function φ is obtained. For obtaining the standard G-L Fokker-Planck equation, third order terms were neglected, as shown in Appendix B. This gives Utilizing the transformation, v = re iϕ , Eq. (9) is transformed to where β = 6γ c , d = (γ a − γ g )/(6γ c ) and q = (A t + D p + D x + γ a − γ c )/4. From standard Ginzburg-Landau laser theory [48], we can get the Ginzburg-Landau potential, and the phonon correlation function g (l) 0 (∞) as where a = −βd/2q, b = β/4q and D ν (z) is the parabolic cylinder function. Substituting the experimental parameters into the Ginzburg-Landau potential, we have plotted the logarithms of the higher order (l = 3, 4) correlations versus the logarithm of the l = 2 correlation in the Fig. 9(a), both for G-L theory and from our experiment in the manner of [48]. Since the correlation functions all approach unity far above threshold, their logarithms approach zero in that limit. Keeping this in mind, we can see from Fig. 9 that the G-L theory matches experimental results only far above threshold. We find this is due to higher order terms having been neglected to arrive at the G-L form of the Fokker-Planck equation, Eq. (9). The effect can also be seen in Fig. 9(b) where the disagreement between G-L theory and experiment worsens for higher-order correlations at low gain modulations, i.e. below threshold.

Conclusion
In this paper, we have investigated equal-time and non-equal-time higher order correlations for the levitated nanoparticle phonon laser using theoretical as well as experimental methods. Utilizing the master equation, the modulation-evolution of phonon number probability distribution was obtained by numerical methods. The phonon number probability distribution satisfies the thermal state distribution when the optical tweezer phonon laser works far below the threshold. When the phonon laser is operated far above the threshold, the phonon number probability distribution approximates the coherent state distribution. These calculations are verified by experimental results.
On the one hand, the equal-time higher order correlations of a phonon laser are numerically and experimentally studied in steady state as well as transient state. When modulation of feedback amplification is far below the threshold, the phonon distribution obeys the thermal state distribution, and higher order correlations of phonons satisfy g (n) 0 = n!. The higher order correlations gradually decrease with the increment of the modulation of feedback amplification until g (n) 0 = 1 as expected for a Poisson distribution. On the other hand, the non-equal-time higher order correlations of the phonon laser are obtained from a quantum Langevin equation correspondign to our master equation; the higher order correlations thus obtained below, near and above the threshold are compared to experimental data. Overall, all numerical results match well with experimental results. Finally, the Ginzburg-Landau theory is used to calculate analytically the higher order correlations of the optical tweezer phonon laser; we find the G-L theory is suitable for the optical tweezer phonon laser far above threshold only.
Our work extends the validity of the optical analogy for a phonon laser to the regime of higher order correlations. The ability to measure and model these correlations, which we have demonstrated, should be useful in physically interesting situations where mechanical nonlinearities can cause the phonons to self-interact, in establishing the validity of various paradigms of laser action (such as Ginzburg-Landau), and in characterizing conditioned measurements (e.g., g (3) (τ 1 , τ 2 , τ 3 )/g (2) (τ 1 , τ 2 ) gives information about three-phonon coincidence conditioned on two-phonon coincidence).

A. Calculation of higher order correlations for experiment and numerical method
Since the experiments are in the classical regime, all dynamical variables commute, and the higher order correlation functions can be represented as: where τ, τ 1 , τ 2 and τ 3 are the time delays between measurements of N, and · denotes an ensemble average. The phonon number N is related to measurable quantities through the relation N = mΩ 0 q 2 / , with q being the displacement of the particle's center-of-mass with respect to the center of the trap.

B. The derivation of Fokker-Planck equation
The master Eq. (1) can be transferred into P representation, and some terms in Eq. (1) can be transferred into P representation as follows, Therefore, the P representation of master equation is as follows, For getting the standard Fokker-Planck equation, the terms −12γ c |v| 2 ∂ 2 ∂v∂v * φ + 3γ c v ∂ 3 ∂v 2 ∂v * φ + 3γ c v * ∂ 3 ∂v∂v * 2 φ which involve third order derivatives were neglected. The complex amplitude ν can be replaced by ν = x 1 + ix 2 . Therefore, and the Fokker-Planck equation becomes For simplicity, the parameters are redefined as follow: ( 1 4 (A t + D p + D x + (γ a − γ g ))6γ c ) 1/2 t ≡ t (24γ c /(A t + D p + D x + (γ a − γ g ))) 1/4 x ≡ x 2(γ a − γ g ) (6γ c (A t + D p + D x + (γ a − γ g ))) 1/2 ≡ a, (18) and the standard Fokker-Planck equation is obtained as follows, The primes were introduced in Eq. (18) have been dropped again for simplicity. By using polar-coordinates x 1 = r cos ϕ,