Spectral analysis and parameter estimation in levitated optomechanics

TL;DR

This paper develops a spectral analysis framework for parameter estimation in levitated optomechanics, enhancing the precision of force measurements and thermodynamic quantity inference using maximum likelihood methods.

Contribution

It introduces a theoretical spectrum model and an optimized apparatus for sensitive measurements, accounting for sinusoidal interferometric signals and enabling improved parameter inference.

Findings

Accurate parameter estimation from spectra using maximum likelihood.

Enhanced sensitivity along one motional degree of freedom.

Successful extraction of thermodynamic quantities in complex regimes.

Abstract

Optical levitation of nanoscale particles has emerged as a platform for precision measurement. Extremely low damping, together with optical interferometric position detection, makes possible exquisite force measurement and promises low-energy tests of fundamental physics. Essential to such measurement is an understanding of the confidence with which parameters can be inferred from spectra estimated from the indirect measurement provided by interferometry. We present an apparatus optimized for sensitivity along one motional degree of freedom, a theoretical model of the spectrum, and maximum likelihood estimation. The treatment accounts for the sinusoidal dependence of interferometric signal on particle position, and we use the technique to extract thermodynamic quantities in a regime where simpler treatments are confounded.

Equations

$$\cos \left[\theta -{\varphi }_0\sin \left(\Omega t\right)\right]=\cos \theta J_0+$$

$$2\cos \theta \sum _{n\ge 2}^{\text{even}}{J}_n\cos (n\Omega t)+2\sin \theta \sum _{n\ge 1}^{\text{odd}}{J}_n\sin (n\Omega t)$$

$$\ddot{z}+\Gamma \dot{z}+{\Omega }^{2}z=w(t)$$

$$S_{zz}(\omega )=\frac{2k_{B}T}{M}\frac{\Gamma }{{\left({\omega }^2-{\Omega }^2\right)}^2+{\Gamma }^2{\omega }^2}.$$

$$R_{vv}(t)=v_0^2\exp \left[R_{\varphi \varphi }(t)-{\Phi }^2\right]$$

$${\sigma }_{vv}(f)=e^{-{\Phi }^2}\sum _{n=0}^{\infty }\frac{{\Phi }^{2n}}{n!}{\sigma }_{zz}(f)\stackrel{n}{⊛}{\sigma }_{zz}(f)$$

$$R_{vv}^{\text{(narrow)}}(t)=v_0^2{e}^{-{\Phi }^2}\left[I_0+2\sum _{n=1}^{\infty }{I}_n\cos n\Omega t\right]$$

$$\mathcal{L}(\mathbf{\alpha })=\sum _i\log S(\mathbf{\alpha })+\hat{S}/S(\mathbf{\alpha })$$

$$S^{\prime }(\Phi ,\Omega ,\Gamma )=\frac{1}{\sqrt{2\pi }R}\int (1+r)e^{-r^2/(2R^2)}\times$$

$$S(\Phi /\sqrt{1+r},\Omega \sqrt{1+r},\Gamma )\mathrm{d}r.$$

$$R_{vv}(t)=v_0^2\exp \left[R_{\varphi \varphi }(t)-{\Phi }^2\right]$$

$$R_{vv}(t)=v_0^2{e}^{-{\Phi }^2}{e}^{-{\Phi }^2\cos (\Omega t)}.$$

$$e^{iz\cos \theta }=J_0(z)+2\sum _{n=1}^{\infty }{i}^n{J}_n(z)\cos (n\theta )$$

$$e^{{\Phi }^2\cos (\Omega t)}$$

Full text

I Introduction

Levitated optomechanics uses a nanoparticle, trapped in vacuum by the optical dipole force, as a harmonic oscillator in a thermal bath Bhattacharya et al. (2017). Despite the apparent simplicity, interferometric position readout Gieseler et al. (2012) combined with optical control, including intensity modulation Gieseler et al. (2012); Vovrosh et al. (2017) and shaping of the focus Ricci et al. (2017), yields a rich platform for exploring nanoscale dynamics Gieseler and Millen (2018); Jain et al. (2016), and promises ultra-precise measurements Ranjit et al. (2016); Hebestreit et al. (2018a) and low energy tests of fundamental physics Geraci et al. (2010); Bateman et al. (2014); Wan et al. (2016); Kaltenbaek et al. (2016).

Properties of the oscillator are not measured directly, but rather are inferred from optical measurements, often with the intermediate step of estimating spectral density from time-series data. This optical measurement is interferometric and therefore not linear with position: the linear approximation is appropriate for geometries where forward scattered light is used Gieseler et al. (2012), and care must be taken when extracting parameters in backwards scattered arrangements Mestres et al. (2015); Rashid et al. (2016); Vovrosh et al. (2017), where phase to position sensitivity can be much greater and there may be a non-zero phase offset. Backscatter is desirable because of increased sensitivity and separation from the laser light, but even when used at low centre of mass temperatures, where the phase excursion is small, the offset can render measurement non-linear.

Furthermore, spectral overlap of oscillator modes, cross-coupling, and effects not accounted for in this simplistic description, notably rotation Monteiro et al. (2018); Rashid et al. (2018); Reimann et al. (2018), can pollute the spectrum and affect parameters extracted by, for example, integrating a truncated region of the spectrum. For example, temperature, or a quantity proportional to it, can be found by integrating the spectral density associated with a given mechanical mode Clerk et al. (2010), but this is only possible when the mode is spectrally resolved, precluding use of this simple technique at high pressure ($\gtrsim 10\,\textrm{mbar}$).

We describe an experimental system which is optimized for detection of one motional degree of freedom (longitudinal) while rejecting, to first order, signal from the other two (transverse). We use backscatter, which offers a larger phase-shift for a given displacement than does forward-scatter, but which in our experiment introduces an uncontrolled and slowly drifting offset phase-shift; and we use an optical fibre system which collects a significant fraction of the scatter and which guides this, without diffractive loss, to detection electronics. This collection and guiding offers alignment stability and contrasts strongly with free-space detection where small photodiode area and the comparatively large free-space beams means that much of the collected light is unused. Using scatter more efficiently is essential if we are to approach the standard quantum limit in these systems Chang et al. (2010).

In this manuscript, we describe our experiment, present a theoretical spectrum, including narrow-band limit and comparison with numerics, and then apply a maximum likelihood approach to parameter extraction from simulated and experimental data. Through this, we observe a centre of mass heating effect at intermediate pressure, where spectral overlap confounds simpler methods, and finally we describe the experimental limits which we encounter in our specific implementation.

II Experimental apparatus

Our apparatus, illustrated in figure 1, consists of a single glass nanoparticle held, by the optical dipole force, in the focus of hemispherical parabolic mirror. Gaussian beam illumination of this mirror is provided by a triplet collimator attached to a single-mode optical fibre. Laser light is provided by a narrow linewidth telecommunications wavelength ($\lambda=1550\,\textrm{nm}$) all-fibre laser, which is amplified by an erbium-doped fibre amplifier and sent, via an optical circulator, to the output collimator. The trapped particle explores $\lesssim\lambda/10$ around the focus, and the potential is therefore well-approximated, for our purposes, as harmonic, with a different natural frequency along each of the Cartesian axes.

The particle Rayleigh scatters a small fraction of the trapping light, half of which is backscattered towards the parabolic mirror. This light is collimated by the mirror and directed towards the fibre optic collimator, some fraction of which is then coupled into the optical fibre. Unscattered light, and forward scattered Rayleigh light, diverges strongly as it travels towards the collimator, and only a small fraction of this is coupled back into the fibre. Imperfections in components of the fibre optic network mean that other stray reflections, notably from the fibre output facet ($\sim 10^{-6}$) and the optical circulator ($\sim 10^{-5}$), make their way, with varying amplitudes and phases, to the photodiode. Overall, light at this frequency, which in other schemes would provide a reference for interferometry, drifts in phase and amplitude because of the macroscopic path difference. This motivates our introduction of phase-coherent light at a shifted frequency, implemented by the acousto-optical modulator (AOM) and known as heterodyne detection, as discussed in detail below.

This arrangement differs in two important ways from the forward-scatter, free-space detection scheme, in which laser light and forward-scattered Rayleigh light falls onto a quadrant detector, with overall signal giving longitudinal information and the difference between left/right or top/bottom pairs giving transverse position information. Firstly, in this backscatter scheme the optical phase of scattered light changes more rapidly with particle position along the optical axis; secondly, the geometry and aggressive spatial filtering by the optical fibre mean that, when correctly aligned, this scheme is first-order insensitive to particle motion transverse to the optical axis. For some applications in precision measurement, this larger sensitivity to one axis with strong rejection of others could be a strong advantage. The fibre optic approach also means that, once coupled into the fibre, there are no diffractive losses and a large fraction of guided light can be directed onto a fast photodiode.

II.1 Interferometric detection

The photodiode signal $V$ arises from the interference of a reference field $E_{\textrm{ref}}$ (phase $\theta$) and scatter from the particle $E_{\textrm{sca}}$ (phase $\phi$). The signal is proportional to the modulus squared of the total field, and contains offset (proportional to the sum of the squares of the individual fields) and a term sinusoidal with phase difference $\theta-\phi$ i.e. $V\propto E_{\textrm{ref}}^{2}+E_{\textrm{sca}}^{2}+2E_{\textrm{ref}}E_{\textrm{sca}}\cos{\left(\theta-\phi\right)}$.

We model the phase shift from particle displacement as linear $\phi=\kappa z$ where the sensitivity $\kappa$ depends on the geometry, discussed below, and $z$ is particle position, along the optical axis, relative to the focus. In the forward-scattered case, $\theta=-\pi/2$ via the Gouy shift, expansion to first-order in $\phi$ is justified because the overall phase excursion is typically small, and hence $V$ is approximately linear with position. For backwards scatter, $\theta$ may drift, and phase excursion is not necessarily small.

II.2 Position to phase sensitivity

Through the focus of a Gaussian beam, the optical phase evolves more slowly than would a comparable non-focussed beam. This Gouy shift, near the focus, modifies the rate of change of phase from $k=2\pi/\lambda$ to $k^{\prime}=k-1/z_{R}$ where $z_{R}$ is the Rayleigh range. Hence, for forward scatter, the rate of change (with particle position) of phase difference between Rayleigh scatter and laser field is $\kappa\approx k-k^{\prime}=1/z_{R}$. For back scatter, laser light accrues phase delay while travelling to the particle, and accrues additional phase delay on its return, giving the larger sensitivity $\kappa\approx k+k^{\prime}=2k+1/z_{R}$.

When focussing with numerical aperture $\gtrsim 0.5$, the paraxial approximation underlying Gaussian beam treatments is not valid Varga and Török (2000). For uniform illumination of a circular aperture Pang et al. (2011), with numerical aperture $\textrm{NA}=1$ as in our parabolic mirror, we find $\kappa=2k-1.41\pi/\lambda$. However, illumination would better be approximated by a truncated Gaussian, and we have not modelled the reflection and collimation of scattered light.

Experimental studies of similar parabolic mirrors have revealed strong sensitivity to alignment and manufacturing tolerances Bahk et al. (2005); Alber et al. (2017). For lower aperture optical systems, or those employing compound lens microscope objectives, we could be more confident of the calculated value, and thereby convert the extracted phase modulation depth into a measurement of the temperature to mass ratio. The proportionality of phase modulation to temperature remains, and we use this result later in this work.

II.3 Heterodyne detection

Interferometric detection in levitated optomechanics is usually by homodyne detection, where scatter interferes with light at the same frequency, i.e. $\theta$ is constant. Here, implementation of this scheme is complicated by the presence of several contributions to light at this frequency, and a consequent drift in phase and amplitude because of the macroscopic (metre-scale) optical paths.

Instead, we employ heterodyne detection, where scattered light is interfered with phase-coherent, frequency-shifted light. In our experiment, this light is derived from the same laser and frequency shifted using an acousto-optical modulator at $80\,\textrm{MHz}$. It is guided to the photodiode and mixed with the scattered light using a 1:1000 coupler, which allows us to retain $\sim 99.9\%$ of the scattered light at this component. The electrical signal from the photodiode is shifted down in frequency to $3\,\textrm{MHz}$ using a radio frequency components and $77\,\textrm{MHz}$ local oscillator. Radio frequency mixing and filtering in any realistic device degrades the signal to noise; implementations without this stage are possible using, for example, a electro-optical phase shifter at $3\,\textrm{MHz}$.

The principle advantage of heterodyne detection in our experiment is that spectrum remains stable regardless of slow drifts in the offset of the phase $\theta$, which now evolves at a rate fast compared with the dynamics of the system. It also allows us to use results from radio-frequency communication to understand the spectrum.

III Theoretical spectrum

For illustration, although we will soon improve this description, we consider purely harmonic particle motion $z=z_{0}\sin(\Omega t)$, which gives phase modulation $\phi=\phi_{0}\sin(\Omega t)$ where $\phi_{0}=\kappa z_{0}$, with a phase offset $\theta$ treated as constant over the timescale of observation. We find

[TABLE]

where $J_{n}$ is the n${}^{\textrm{th}}$ order Bessel function evaluated at $\phi_{0}$. This description has been used to understand optomechanical spectra, including a case where $\theta$ was varied systematically Vovrosh et al. (2017). The (co)sinusoidal dependence of (even)odd harmonics is a feature of this expansion being around zero frequency, and a useful interpretation is of negative orders being reflected about the origin to overlap and interfere with their positive frequency counterparts.

We now improve the description of particle motion. To obtain a spectrum, we necessarily observe for a time long compared with the relaxation time. Therefore, it is not accurate to treat motion as purely harmonic. Hereafter, we treat the particle as a stochastic harmonic oscillator, which leads to an important difference. Moreover, to avoid interference of negative and positive orders, we expand about a frequency large compared with the width of the spectrum, by using $\theta=\omega_{0}t+\theta_{0}$.

The differential equation describing a stochastic harmonic oscillator, with natural frequency $\Omega$ and damping $\Gamma$, is

[TABLE]

where $z$ is the particle position, over-dot means differential with time, and $w$ is a Wiener process with volatility $k_{B}T\Gamma/M$ chosen to agree with equipartition; $k_{B}$ is Boltzmann’s constant, $T$ is centre of mass temperature, and $M$ is the particle mass. From the fluctuation–dissipation theorem and linear response theory, the spectrum of fluctuations is

[TABLE]

For forward scatter at low temperature, where $\theta=-\pi/2$ and $\phi_{0}\ll 1$, detection is approximately linear and it is often sufficient to approximate the signal as a scaled version of equation 3. For backward scatter, where these conditions are not satisfied, previous work has combined equations III and 3 in a heuristic way, by considering a scaled delta function at each frequency, broadened by the position spectrum Vovrosh et al. (2017). Here, instead, we use a result from radio communications which gives an exact result.

The correlation function $R_{vv}(t)=\langle v(t)v(t-\tau)\rangle_{\tau}$ of a signal $v(t)=v_{0}\sin\left[\omega_{0}t+\phi(t)\right]$ phase-modulated by a Gaussian random process $\phi(t)$ with correlation function $R_{\phi\phi}(t)$ and variance $\Phi^{2}=\kappa^{2}\langle z^{2}\rangle$ is Godone et al. (2008)

[TABLE]

the Fourier transform of which is the spectrum and can be represented as Middleton (1951); Brochard et al. (2017)

[TABLE]

where $\sigma_{zz}(f)=2\pi\,S_{zz}(2\pi f)/\langle z^{2}\rangle$ is the normalized position spectrum (such that $\int_{-\infty}^{+\infty}\sigma_{zz}(f)df=1$), and $\overset{n}{\circledast}$ is the $n^{\textrm{th}}$ order convolution infix operator defined recursively via $a\overset{n}{\circledast}a=a\circledast a\overset{n-1}{\circledast}a$, with $a\overset{2}{\circledast}a=a\circledast a$ the standard convolution, $a\overset{1}{\circledast}a=a$, and $a\overset{0}{\circledast}a=\delta$ the Dirac delta.

The overall appearance is of peaks near integer multiples of natural frequency $\Omega$, as in equation III. Higher order terms give contributions to lower order harmonics, and hence the energy (and information) contained within a given harmonic is a result of a series summation, not just the contribution from a given order. This expression valid for under- and over-damped oscillators.

III.1 Illustrative, typical spectra

An example of a typical spectrum is shown in figure 2. Also shown here is a spectrum estimated from a simulated realisation of a random process; spectral estimation is discussed further in Section IV.

The spectrum is centred on the chosen modulation frequency, $\omega_{0}=2\pi\,\times\,3\,\textrm{MHz}$, with peaks spaced at integer multiples of $\Omega$ symmetrically either side. The amplitude of higher order peaks reduces monotonically and series convergence is assured because of the factorial in the denominator.

The overall, broad pedestal arises from the increasingly broad contribution of higher order convolutions to the overall spectrum. The presence of this pedestal makes clear that the integrated area under the first order peak is not necessarily an accurate representation of the variance in $\phi(t)$, especially because the shape of the spectrum will change as this variance is e.g. reduced by cooling.

III.2 Narrow-band limit

Often, with this kind of system, one has a narrow-band process. Modelling this as purely sinusoidal gives Equation III; instead, we observe that the limiting case for low damping is that the process has a sinusoidal correlation function $R_{\phi\phi}=\Phi^{2}\cos\Omega t$ for which, using Equation 4 and the Jacobi–Anger identity with imaginary amplitude, we find

[TABLE]

where $I_{n}$ is the modified Bessel function of the first kind evaluated at $\Phi^{2}$. Details of the algebra are given in Appendix A. Previous workMestres et al. (2015); Vovrosh et al. (2017) used the (non-modified) Bessel functions to estimate physical parameters from observed spectra, which agrees with Eq. 6 for small $\Phi$.

The ratio of predicted peak amplitudes, relative to the first order, is illustrated in figure 3, and the treatments agree for small $\Phi$. Notably, treating motion as purely harmonic predicts that the first order will vanish at $\Phi\approx 3.8$, the first zero of $J_{1}$; the modified Bessels have no such zero crossings.

IV Parameter estimation from time series

Our goal is to estimate model parameters describing the process from the time-series measurements of a phase-modulated signal arising from a realisation of this process. Ideally, we would compute the likelihood of the time-series data for given model parameters, and thereby infer the probability density for these parameters Sivia and Skilling (2006). However, the experimental spectrum, while relatively clean, contains features not described by this model, and these are easily filtered in the spectral domain. Therefore, we first make a non-parametric estimate of the spectrum, and then compute the likelihood of this spectrum, over the relevant regions, for given model parameters. Future work may address this parametric estimation problem without the intermediate non-parametric spectral estimation step. There is interest in this approach, but to our knowledge the existing treatments (e.g. autoregressive maximum likelihood; other Bayesian methods Singh et al. (2018)) are not applicable when the measurement of the process is non-linear.

We base our spectral-domain approach on Whittle’s approximate log-likelihood Taniguchi and Kakizawa (2000):

[TABLE]

where $\bm{\alpha}$ are the model parameters, $S$ is the theoretical spectrum, and $\hat{S}$ is an estimator of the spectrum from the time-series. Summation is over the discrete frequencies at which the spectrum is estimated. The probability density for the parameters given the data Sivia and Skilling (2006) is $\textrm{prob}(\bm{\alpha}|\hat{S})\propto\exp\left[-\mathcal{L}(\bm{\alpha})\right]$.

The Whittle likelihood is an approximation of the true likelihood for a stationary Gaussian time-series model. While the underlying harmonic oscillator can be described by such a model, our measurement of it cannot. Therefore, it is unclear whether the Whittle likelihood will give an accurate estimate in this case. We make a slight adaption (described below) and verify the effectiveness numerically in our use case.

The Whittle likelihood was originally formulated for the periodogram, the modulus squared of the discrete Fourier transform. This estimator is asymptotically unbiased, but it is not consistent: the variance does not decrease for a large number of points. We trade points for consistency by averaging the periodogram using Bartlett’s method: split the time-series data of length $N$ into many $M$-length segments, each of which is windowed using the Tukey–Hanning window, and then compute the average periodogram over these segments. The distribution of this estimator relative to the true value is Priestley (1981) $\nu\hat{S}/S\sim\chi^{2}_{\nu}$ where this is the $\chi$-squared distribution, and the degrees of freedom $\nu$ is twice the number of segments, $\nu=2N/M$; this is illustrated in figure 4. For large $\nu$, this tends to a Gaussian distribution $\hat{S}/S-1\sim\mathcal{N}(0,\sqrt{2/\nu})$

IV.1 Parameter estimation with simulated trajectories

We compute the likelihood function and consequent probability density for a simulated trajectory with known parameters $\bm{\alpha}_{0}$, as we vary $\bm{\alpha}$ about the true value. The physical parameters describing our spectrum are $\bm{\alpha}=(\Phi,\Omega,\Gamma)$. For illustration, since $\Omega$ is well constrained, we vary $(\Phi,\Gamma)$ while holding $\Omega=\Omega_{0}$. Details of the simulation are described in Appendix B, and the calculated probability density is illustrated in figure 5. The central spectral peak, the delta function resulting from $n=0$ in equation 5, contains no information about the random process, is relatively large, and depends critically on the spectral windowing function; we therefore exclude a small region near to this peak from our likelihood calculations.

To assess whether this Whittle probability density is an accurate representation of the information which can be extracted through this process, we we create an ensemble of simulations with known $\bm{\alpha}$, find the maximum likelihood estimate in each case $\bm{\alpha}_{\textrm{MLE}}$, and estimate properties of the assumed Gaussian probability distribution from which these estimates are picked. This ensemble estimate of the probability distribution is compared with that obtained by computing $\exp\left[-\mathcal{L}(\bm{\alpha})\right]$ directly.

Rather than raster the parameter space, we compute the profile likelihood, by constraining one element of $\bm{\alpha}$ and fitting all others. This gives access to the marginal probability, effectively integrating over the other (nuisance) parameters. The difference between computing marginal probability density $\textrm{prob}(\Phi)=\int\textrm{prob}(\Phi,\Gamma)\,\mathop{}\!\mathrm{d}\Gamma$ and evaluating $\textrm{prob}(\Phi,\Gamma=\Gamma_{0})$ is apparent in figure 5, where the extremal $\Phi$ accessible for a given probability density is larger if we integrate over $\Gamma$ (marginal distribution) or allow $\Gamma$ to be adjusted (profile likelihood), compared with constraining $\Gamma=\Gamma_{\textrm{MLE}}$.

V Application to experimental data

Experimental measurements are subject to additional complications not captured by the simulations. The collected optical power and the responsivity of the photodetector give some scaling to the recorded voltage signal, and an advantage of the technique is that information is encoded in the spectral shape, not the absolute scale. In addition, the detection system introduces measurement noise, which we treat as white, and which therefore manifests spectrally as a constant offset. We treat this scaling and offset as nuisance parameters by finding their maximum likelihood values for each calculation of $S(\bm{\alpha})$, thereby computing the profile likelihood as described in Section IV.1.

An example of a typical maximum likelihood fit at a reasonably high pressure, where overlap between peaks is significant, is shown in figure 6. Normalised residuals are shown and, because we are averaging over a large number of periodograms, their distribution approaches Gaussian.

V.1 Probability densities and parameter uncertainties

Once parameters $\bm{\alpha}_{\textrm{MLE}}$ have been found which maximize the Whittle likelihood, we compute the probability density for a given parameter by constraining this parameter (to some value near to the maximum likelihood value) and then maximizing probability by adjusting all other parameters. Moreover, when the fitted parameters depend monotonically on the constrained parameter, then this calculation also reveals the probability density for the fitted parameters. An example of probability densities computed in this way is shown in figure 7.

For spectra with relatively high damping ($\Gamma\gtrsim 1\,\textrm{kHz}$) our model describes the experimental spectrum well, and we are able to extract $\Phi$ and the associated uncertainty in a regime where spectral overlap with harmonics would confound the naive approach of integrating under peaks. For spectra with lower damping, the intensity stability of our experiment affects the peak shape, principally through the square-root dependence of $\Omega$ on power; this limitation is particular to our apparatus, does not limit the technique in general, and is discussed further in Section V.3. For $\Gamma\lesssim 150\,\textrm{Hz}$ we expect an effect from windowing of our finite-length time series, and this could be included in a description of the theoretical spectrum.

V.2 Observing heating at intermediate pressure

We apply our estimation technique to spectra obtained for a nominally $100\,\textrm{nm}$ silica particle in a dipole trap, as we reduce the gas pressure in the vacuum chamber. Probability density for $\Phi$ is found by computing the Whittle likelihood as this parameter is constrained and all others are adjusted to maximize this likelihood. The best-estimate and uncertainty are calculated from this probability distribution in the standard way, and results are plotted as a function of pressure in figure 8. The three points at highest pressure (around $100\,\textrm{mbar}$) are derived from probability distributions which are not Gaussian, but are sufficiently broad so as to not be misleading when described by symmetric errorbars.

As discussed in Section II.2, the phase to position sensitivity is not known accurately in this system. The phase modulation depth is proportional to temperature $\Phi^{2}\propto T/(M\Omega^{2})$ and hence, because $\Omega$ is well-constrained and we assume particle mass $M$ remains unchanged, we interpret the slope near $10\,\textrm{mbar}$ as an increase in temperature. Comparing the relative plateaus near $20\,\textrm{mbar}$ and $2\,\textrm{mbar}$, we estimate an increase in temperature of approximately $90\,\textrm{K}$.

Decrease in mass would also manifest as an increase in $\Phi$. Recent experiments on similar systems have revealed that micron-sized particles are porus and can contain significant water Blakemore et al. (2019). However, the strong absorption of $1550\,\textrm{nm}$ light by water suggests that there would be low content even at atmospheric pressure. We cannot rule out the possibility that mass reduces, and future experiments, perhaps cycling pressure with different background gases, may be informative.

Some experiments in levitated optomechanics have suffered from an unexpected increase in particle loss probability at these intermediate pressures, and there has been work to understand temperature in this settingMillen et al. (2014); Hebestreit et al. (2018b). Experiments using telecommunications wavelength $1550\,\textrm{nm}$ (rather than $1064\,\textrm{nm}$) appear to suffer less from this unexplained loss, and this is assumed to be because of the lower material absorption of silica at this wavelength. This technique may be a useful tool to estimate temperature changes in this regime, with applications including diagnosing material properties of fabricated nanoparticles designed to minimize heating caused by laser absorptionFrangeskou et al. (2018).

V.3 Intensity noise and non-linear broadening

An additional complication, relevant at low $\Gamma$, is relative intensity noise: the laser intensity at the focus has some small, low-frequency drift, and this affects properties of the spectrum. For a fibre laser, the noise spectrum is extremely quiet at or above the particle oscillation frequency $\Omega$, making this a good choice to minimize parametric heatingSavard et al. (1997); it is more significant at low frequencies, corresponding to slow drifts during data collection.

From direct measurement of the intensity at points in our fibre network we constrain $R\lesssim 1\%$ and therefore this is significant only for $\Gamma\lesssim\Omega/100\sim 10\,\textrm{kHz}$.

To model this effect, we assume that drifts in intensity are slow compared with the relaxation time of the oscillator $\Gamma^{-1}$ and that the distribution of intensity is Gaussian with some width $R$. This model is unlikely to be sufficient for precision measurements, and future experimental work must be undertaken to minimize this drift. Under these assumptions, we can describe the observed spectrum as an average over the intensity distribution.

The intensity affects several aspects of spectrum: the overall scale (linearly); the natural frequency $\Omega$ (square root); and the modulation depth $\Phi$ (square root) via the natural frequency because the spatial extent of the thermal state depends inversely on $\Omega$. Therefore, the overall broadened spectrum is

[TABLE]

An example spectral peak, with experimental, unbroadened theoretical, and broadened theoretical according to equation V.3 is shown in figure 9. The model parameters are found by maximizing the likelihood over the entire spectrum, and a zoomed in region, where the effect is most visible, is shown for illustration. While this Gaussian broadening captures the behaviour well, and quite often fits our experimental results, care must be taken if this is used for parameter extraction because the duration of the measurement is not sufficiently long that the intensity distribution can be reliably approximated as Gaussian; the dynamics of the drift are too slow. Sometimes, for example, intensity undergoes a linear drift, which results in an asymmetric peak; this might explain examples in the literature, e.g. figure 3 in Ref. Vovrosh et al., 2017. Estimates based on the integrated area are unaffected, but the use of this whole spectrum method must account more carefully for any such drift.

An additional source of broadening is the thermal average of the Duffing non-linear frequency shift Gieseler et al. (2013). This effect, and the consequent distinctive asymmetric peak shapes, is masked in our system by the slightly larger intensity noise. It is straightforward to include this non-linear broadening by averaging the heterodyne spectrum over the Boltzmann distribution, with frequency shift $\Delta\Omega\propto E$ and oscillation amplitude variance $\Phi\propto\sqrt{E}$, similarly to equation V.3. In constrast with slow intensity drifts, it is reasonable to sample sufficiently long that the Boltzmann average is a good approximation, and therefore we expect, although cannot currently verify, that our parameter estimation approach remains valid.

VI Conclusions

We have presented a technique for extracting, with confidence intervals, thermodynamic quantities from interferometric position measurements of a levitated nanoparticle by careful treatment of the estimated spectral density. The techniques relies on the shape of the spectrum, and is indifferent to calibration of the photodiode responsivity or changes in the signal amplitude. We have demonstrated this technique with experimental apparatus which is long-term stable and optimized for sensitivity along one direction, with strong rejection of others, giving a spectrum well described by the model. The technique allows extraction when spectral features are not well resolved, and permits spectral windowing, which means it can be used when experimental spectra are cluttered by unmodeled features.

We have used this technique to observe centre of mass heating at intermediate pressure, where the simpler technique of integrating under a peak is not appropriate. This technique may find use in diagnosing temperature dependence in this pressure range, for example when characterising unwanted heating in carefully fabricated extremely pure nanodiamonds. Alternatively, it may be used if the temperature is known and the mass changing, such as by deliberate evaporation of a nanoparticle to obtain small trapped nanoparticles, with thermalisation via buffer gas.

This implementation was limited by intensity noise, but this is not a fundamental limitation of the technique. Further, one could in principle record for sufficiently long that the histogram of intensity fluctuations is well described by a Gaussian, for which the model would then be expected to fit, but a reduction in the noise is a more efficient approach.

The high aperture optical trap, with possible manufacturing and experimental imperfections and an incomplete model of the focussing and collection optics, means that we cannot, with confidence, calculate the phase to position sensitivity in our experiment; if this were better known, by either direct measurement, calculation, or using different optics, the technique would allow for direct calculation of the temperature to mass ratio from the extracted phase modulation depth.

This work focussed on heterodyne detection, which can be implemented fully optically with little loss in signal quality. However, many existing experiments use homodyne detection, and future work is to extend the formalism to cover this case.

Acknowledgements

This work was supported by startup funding from the College of Science at Swansea University. The authors would like to thank George Winstone for helpful discussions.

Appendix A Narrow band spectrum

Starting with equation 4 for the correlation function $R_{vv}$ of the phase-modulated signal $v$,

[TABLE]

we use a narrow-band process $\phi$ which has correlation function $R_{\phi\phi}=\Phi^{2}\cos{\Omega t}$ to find

[TABLE]

The modulated exponential term can be expressed using the Jacobi–Anger identity with the replacements $iz=\Phi^{2}$ and $\theta=\Omega t$:

[TABLE]

where $J_{n}$ is the $n^{\textrm{th}}$ order Bessel function. Hence,

[TABLE]

where $I_{n}$ is the $n^{\textrm{th}}$ order modified Bessel function, and we have used the identities $J_{n}(ix)=i^{n}I_{n}(x)$ and $I_{n}(-x)=(-1)^{n}I_{n}(x)$.

Appendix B Ensemble simulation comparison

We simulate the process using the Euler–Maruyama method with $\delta t=1\,\textrm{ns}$, and sampled at $10\,\textrm{MS}/\textrm{s}$ for $1\,\textrm{s}$ to generate $10\,\textrm{M}$ points. (The crudeness of the numerical method makes necessary the short time-step.) For an ensemble of trajectories ($\sim 40$), we estimate the spectrum with Bartlett’s method as described in the main text using $M=2^{16}$, and find the maximum likelihood estimate of the parameters, $\bm{\alpha}_{\textrm{MLE}}$.

We use $\Omega=2\pi\times 100\,\textrm{kHz}$, $\Gamma=(10,20,50,100)\,\textrm{kHz}$, and use $\Phi=0.1,0.2,0.5,1.0$ for calculating the phase-modulated spectrum. We use a temperature of $300\,\textrm{K}$ and a mass of $10^{-18}~{}\textrm{kg}$; these are typical experimental values, and only their ratio enters the simulation. In keeping with the experimental system, we use a centre frequency $f_{0}=3\,\textrm{MHz}$. Since the physical parameters are encoded in the shape of the spectrum, and are not reliant on any scaling, we use a unity amplitude sinusoid.

For this simulation we expect full agreement across the spectrum. For the experimental system, regions between the peaks are slightly polluted with either second-order transverse motional peaks, or cross-coupling terms. Therefore, for low $\Gamma$ where these additional unmodeled peaks become visible, it becomes necessary to window the spectrum. This has little effect on the probability density, since these regions, at low $\Gamma$, do not depend strongly on the model parameters, and so have little influence (other than an inconsequential constant offset) on the Whittle log-likelihood. We explored this windowing for low $\Gamma$, but found our experiment limited by intensity noise, as described in Section V.3, and so have not used any selective windowing in the work presented.