Using the transient trajectories of an optically levitated nanoparticle to characterize a stochastic Duffing oscillator

We propose a novel methodology to estimate parameters characterizing a weakly nonlinear Duffing oscillator represented by an optically levitating nanoparticle. The method is based on averaging recorded trajectories with defined initial positions in the phase space of nanoparticle position and momentum and allows us to study the transient dynamics of the nonlinear system. This technique provides us with the parameters of a levitated nanoparticle such as eigenfrequency, damping, coefficient of nonlinearity and effective temperature directly from the recorded transient particle motion without any need for external driving or modification of an experimental system. Comparison of this innovative approach with a commonly used method based on fitting the power spectrum density profile shows that the proposed complementary method is applicable even at lower pressures where the nonlinearity starts to play a significant role and thus the power spectrum density method predicts steady state parameters. The technique is applicable also at low temperatures and extendable to recent quantum experiments. The proposed method is applied on experimental data and its validity for one-dimensional and three-dimensional motion of a levitated nanoparticle is verified by extensive numerical simulations.

Scientific RepoRtS | (2020) 10:14436 | https://doi.org/10.1038/s41598-020-70908-z www.nature.com/scientificreports/ and manipulation in the field of nanotechnology and chemistry. It is therefore desirable that the methodology characterizing parameters of a nonlinear oscillator is applicable also in quantum mechanics. Moreover, the methodology should apply to transient effects, in contrast to the widespread approach based on PSD. The levitating nanoparticle system is specific because the nanoparticle position is typically the only directly observable quantity. The velocity can be in principle determined using an independent Doppler velocimetry 38 or by homodyne detection 39 but this technique requires a high speed detection and data acquisition system. An object velocity or momentum is most frequently estimated by a numerical differentiation of measured positions of a levitating particle. Therefore, a method which uses information mainly from the particle position is generally more broadly applicable.
In this paper, we present a suitable approach estimating all desired parameters of a weakly nonlinear Duffing oscillator assuming that the particle radius is known. Our estimation method is based on a post-processing of acquired time records of particle positions and averaging trajectories in the statistical ensemble. We compare our results with the commonly used PSD method and, unlike Refs. 18,40 , we obtain the Duffing coefficient of nonlinearity without any external driving force which could affect the system parameters and would be undesirable for experimental studies of quantum effects. Moreover, the determined damping coefficient follows well the theoretical prediction down into low pressures. Our method is capable of the determination of parameters on the time scale that is faster than the heating rate which is crucial for experimental testing of transient stochastic phenomena with levitating particles.

experimental set-up and data processing
A nanoparticle (NP) is trapped in a focused laser beam in all three orthogonal directions. A scheme of an optical trap in a transversal direction is shown in Fig. 1a. Since the non-conservative scattering force is proportional to a sphere radius a as a 6 41, 42 , it can be neglected with respect to the conservative gradient forces acting upon an NP. Thus the spatial profile of the trapping potential U x, y, z follows an inverted shape of the optical intensity I x, y, z as U x, y, z ≃ −I x, y, z . Since the real lateral or longitudinal intensity profile is close to the Gaussian or Lorentzian shape 43 , respectively, the potential energy of the NP displaced further from the equilibrium position is lower in the real optical trap comparing to the ideal parabolic potential (compare solid red and dashed blue curves in Fig. 1a). Such slight deviations from an ideal shape are frequently neglected for an NP with cooled translational motion 28 because it moves in the dominantly parabolic bottom of the potential well, and the NP trajectory is analyzed following the theory for an ideal harmonic oscillator 29 . However, if NP's motion is cooled but coherently excited, the nonlinearities in the potential arise and drag the NP motion out of ideal harmonic oscillations. It is the same nonlinearity affecting the NP at higher temperatures. Assuming just the first lowereven-order terms in Taylor series expansion of the real profile of the potential near its minimum 18,44 , see also Supplementary Information, it gives the Duffing type nonlinearity. The nonlinearity in the system can be characterized by the stiffness of the Duffing spring that depends on the position x of the object as κ D = κ 1 − ξ x 2 , Figure 1. (a) A scheme of a nanoparticle (NP) levitating in one-dimensional optical potential along the lateral axis x . Dashed blue and solid red curves denote the parabolic shape of the ideal harmonic oscillator and nonlinear Duffing potential, respectively. Initial conditions of the particle motion are denoted as (x 0 , v x0 ) . (b) Experimental set-up. A low-noise y-polarized laser beam (wavelength 1064 nm, Mephisto 2000NE) is 3× expanded using lenses L1 and L2 (Thorlabs AC127-025-C, AC254-075-C). The trapping power is controlled by a rotating half-wave plate / 2 (Thorlabs WPH10M-1064) placed in front of a polarizing beamsplitter PBS (Thorlabs PBS253). The beam is focused by an aspheric lens L3 (Lightpath 355330, NA = 0.77) placed inside a vacuum chamber and maximal beam power 100 mW can be obtained here. Silica NPs (170 nm in diameter, Bangs Laboratories, Inc.) are launched from a silicon substrate towards the beam focus by a focused laser pulse (wavelength 532 nm, energy 4 mJ in a single pulse, Continuum Minilite MLII) 60 under a chamber pressure of 7 mBar. The unscattered and scattered light by the NP trapped near the beam focus is collected by an aspheric lens L4 (Thorlabs C240-TME) in forward direction, demagnified by a telescope L5 and L6 (Thorlabs AC254-250-C, AC254-030-C), and detected by a quadrant photodetector QPD (Hamamatsu G6849). The QPD signal is processed by a home-made electronics and the NP positions in all three axes are acquired by National Instruments FPGA card (NI FPGA NI 5783 adapter module, FlexRIO FPGA) with the sampling frequency 1.78 MHz.
Scientific RepoRtS | (2020) 10:14436 | https://doi.org/10.1038/s41598-020-70908-z www.nature.com/scientificreports/ where κ is the stiffness of the ideal harmonic oscillator and ξ represents the coefficient of Duffing nonlinearity. In the real optical traps the coefficient ξ is positive and thus the stiffness of the optical spring gets weaker with an increasing deviation from the equilibrium position x = 0 . Such an oscillator is called a softening Duffing oscillator 14,18,45,46 . Technical details of the whole experimental set-up are described in Fig. 1b. An example of a 1D trajectory of an optically levitated NP in a potential well formed by a single strongly focused Gaussian laser beam is plotted in Fig. 2a. It shows that the NP oscillates with an amplitude that varies in time. Figure 2b shows two magnified regions of Fig. 2a (blue) and corresponding NP instantaneous velocity (orange) determined as a central difference v x (t) = [x(t + �t) − x(t − �t)]/2�t , where t is the time step in data acquisition given by the sampling frequency �t = 1/f sample = 1/1.78 μs which is sufficiently high comparing to the examined processes (~ 90 kHz). In the first example, the NP "wiggles" with a low amplitude ( |x| < σ x , where σ x is the standard deviation of the particle position) close to the bottom of the potential (center of the optical trap) while the second part of the trajectory shows the sustained oscillations with a higher amplitude ( |x| > σ x ). At the same time scale, small-amplitude motion is strongly affected by noise whereas the largeamplitude oscillations are rather modified by the nonlinear potential softening. Furthermore, Fig. 2c shows a two dimensional histogram of the Gaussian shape obtained from a long-time data set (30 s) showing probability density function of an NP being in a given volume of a phase space-i.e. having given position and velocity in the x axis.
A standard steady-state analysis of NP motion in a harmonic potential is based on the power spectral density (PSD) of NP positions described by the following function 10,24,26,28,47,48 where k B , T SP , m , � 0 = 2πf 0 , ω , Ŵ , and P s denote the Boltzmann constant, effective spectral temperature of the thermal bath driving an NP with a white noise distribution, particle mass, NP eigenfrequency (in rad/s), frequency, medium damping, and technical measurement noise 29 , respectively. PSDs calculated from the long data set for two different ambient pressures are shown in Fig. 2d. For ambient pressure of 10 mBar and higher we obtained almost perfect fit of Eq. (1) to experimental data because the peak broadening is caused predominantly by the damping Ŵ and the resonant peak is not visibly influenced by the nonlinear behaviour or crosstalks with other axes. In contrast in the lower pressures noticeable differences appear: (1) the measured peak is wider and lower than the theoretically predicted one, (2) there is a clear asymmetry of this peak and it exhibits a visible negative skewness. Such a deviation of the experimental PSD profile from the ideal one can be caused by the softening of the experimental potential profile from the assumed ideal harmonic trap 18,29,44,49 . If we consider the lowest order of nonlinearity characteristic for a Duffing oscillator (see Fig. 1a), the oscillation frequency depends on the amplitude of oscillations 18,44,45 as experimental data in Fig. 2f demonstrate. This behaviour leads to a continuum of new frequency components that contribute to the PSD at frequencies below the PSD resonant frequency (main peak) and cause asymmetry of the PSD peak observed in Fig. 2d for the lower pressure 49 . The secondary peaks observed in Fig. 2d are caused by crosstalks between the measured x axis and other coordinate axes y and z and their higher harmonics. Figure 2e reveals experimental evolution of the main peak for more values of the ambient pressure with two highlighted levels of pressure discussed above.
Obviously, fitting the steady-state PSD of a real nonlinear oscillator with Eq. (1), which is valid only for the ideal harmonic oscillator, has to lead to an imperfect determination of parameters of the real oscillator 29 , as we show later in Fig. 5. Moreover, this PSD approach is not suitable for current investigation and use of transient out-of-equilibrium coherent effects faster than any heating of the motion 40 . After the cooling of levitating systems to the ground state 37 , such estimation of the Duffing oscillator from fast transient effects are crucial for upcoming studies of quantum effects [50][51][52][53][54] .
Therefore, in this work, we propose a novel method for determination of parameters of an optically levitated NP. Unlike the methods based on equipartition theorems in equilibrium steady-state or PSD methods that are independent of initial conditions, we post-process measured stochastic trajectories for selected suitable initial conditions and follow an NP trajectory during a short transient dynamics to determine required parameters. However, we study a stochastic system where a random noise influences the individual trajectories of the NP probing a nonlinear potential. Thus each transient process starting at the same point in the phase space (x 0 , v x0 ) leads to a different trajectory. Therefore, we opted for an analysis of the moments (e.g. mean and variance) of such ensemble of NP trajectories that start in the same volume of the phase space. Statistical moment dynamics can be further used up to the ground state of mechanical motion.
The consequent post-processing algorithm for one axis (shown for the x axis) proceeds in the following steps: 1. Record of QPD signals in voltage is transformed to the NP positions x(t) by the calibration factor C V→m described in Methods. 2. The NP positions acquired at low pressures ( p < 1 mBar) are filtered by a frequency domain bandpass filter (passband for the x axis: 87-96 kHz, passband for the y axis: 72-90 kHz) in order to keep only the leading oscillation in the selected axis because at low pressures crosstalks between different axes become visible in PSD and influence the analyses of particle motion in a selected axis. 3. Using the filtered NP positions x(t) determined in equidistant time steps t , we calculate NP velocities v x (t) , where t is given by the sampling frequency during the data acquisition �t = 1/f sample = 1/1.78 μs. Furthermore, we verified by the computer simulations that the experimental sampling frequency is sufficient to calculate the instantaneous particle velocity even at atmospheric pressure.
(1)  www.nature.com/scientificreports/ 4. We look for an event when both position and velocity pass through a small region of the phase space where both x 0 and v x0 are small but slightly bigger than the experimental uncertainty of measured positions and calculated velocities. If we find several consecutive points in this region, the selected one gives the minimal separation from x 0 . The uncertainties were estimated by the high frequency tail of the PSD at frequencies above 660 kHz (see Calibration of quadrant photodetector in Methods) and they correspond to the white noise having standard deviation ∼ 2 nm, which is about 3% of the particle standard deviation from the optical trap center. 5. When such an event is detected, we assume it is the beginning of a new trajectory corresponding to t = 0 and we take a snippet of an on-going trajectory of length δt (typically δt ≃ 100 μs). This snippet is added into a statistical ensemble of trajectories starting at the vicinity of the point (x 0 , v x0 ) in the phase space. 6. Steps (4) and (5) are applied on the remaining part of the acquired NP record while over-lapping parts of the trajectories are not included. 7. Averaged trajectory x(t) , velocity v x (t) and their variances Several examples showing such averaged trajectories starting at points (x 0 , 0) are shown in Fig. 2f where the damping of oscillator amplitudes can be clearly seen as well as the frequency dependence on the oscillator amplitude. Furthermore, Fig. 2g shows the increase in the position variance Var(x) determined for the particular initial condition (x 0 , v x0 ) = (0, 0) ± 2 nm, 1 mm s −1 . After a short transient process the variance achieves a saturated value which corresponds to the thermalization of the particle motion to the surrounding bath.

Methods for determination of parameters of the Duffing oscillator
Below we present a description of four methods we use for estimation of the parameters of the nonlinear Duffingtype oscillator. At the end of each subsection describing one method the procedure for the parameter determination is summarized.
The Duffing oscillator approximation (DOA). In many types of nonlinear oscillators with a harmonic potential minimum, a nonlinearity only weakly modifies harmonic oscillations. In this case it is sufficient to add only the first nonlinear term to the harmonic potential profile description given by a Taylor series expansion 45,55 . On time scales shorter than any thermalization, the motion of the so-called Duffing oscillator can be described by the deterministic Duffing equation (DDE) in the following form where x denotes the particle position along one axis, Ŵ = γ /m denotes the medium damping with m and γ being the mass of the oscillator (e.g. a levitating NP) and drag coefficient of the medium (e.g. air at low pressure), respectively. � 0 = √ κ/m denotes the eigenfrequency of the ideal harmonic oscillator with the oscillator stiffness κ (e.g. the stiffness of the optical trap). ξ represents the coefficient of Duffing nonlinearity and is related to the third order coefficient of the Taylor expansion of the optical force near the potential minimum at x = 0 . The softening Duffing oscillator described in this way will be subject of our following analyses.
Local analysis near x = 0 of such Duffing oscillator for weak damping Ŵ leads to a solution of the lowest perturbation order in the form of a damped oscillator 56 where x 0 and θ 0 are given by the initial conditions at t = 0 . The frequency D can be expressed as 56 Considering only a weakly damped nonlinear oscillator we can define the eigenfrequency of the damped harmonic oscillator H depending on its initial amplitude as In reality, the micro-or nano-oscillators are driven by thermal noise which is generally assumed as the white noise and the DDE modifies to the Langevin type of stochastic Duffing equation (SDE) in the following form 18,45 (2) denotes correlation of functions within brackets and δ(x) denotes the delta function. This equation is later used for numerical simulations of the dynamics, results are presented in Supplementary Information. Even though Eq. (7) assumes the driving force with a white noise spectrum, the simulations and experimental activities can be extended to a driving force with a coloured noise spectrum, e.g. using an external force with well-controlled frequency spectrum.
Idealized harmonic oscillations. For initial conditions close to the potential minimum, the weakly damped nonlinear oscillator behaves predominantly as harmonic 18 and a transient motion can be well described by a harmonic oscillator with an eigenfrequency 0 continuously damped with rate Ŵ and excited by thermal surrounding characterized by the effective transient temperature T TR . It is sufficient to select initial conditions at the potential minimum with zero speed and observe heating and random build-up of many linearized small oscillations. We use the post-selection process described in the previous section, in order to virtually cool down the NP and afterwards the NP motion under the influence of the heat transfer from surroundings is analyzed. This is equivalent to a direct observation of the particle heating in experimental systems 57 . If we select the initial conditions at (x 0 , v x0 ) = (0, 0) ± 2 nm, 1 mm s −1 and determine Var(x) from the experiment for long enough acquisition time to cover the thermalization, we can fit the variance profile very well using the analytical solution of the linearized harmonic oscillator 24 From the fit to the experimental data (see Fig. 3a) we can determine all three parameters 0 , Ŵ and T TR . These parameters remain constant even if the NP enters the nonlinear region. As a cross-check, we can use the variance of velocities, which provides us with the same results. In a weakly damped system with Ŵ ≪ � D , we can use a simplified formula without oscillations in order to quantify the heating of the system, i.e. medium damping as well as the effective temperature using Such fit to the experimental data for a weakly damped oscillator is shown in Fig. 3b.
A few nonlinear oscillations. A nonlinearity in the potential profile influences the NP motion if an NP initial position is placed far from the potential minimum. If the damping is weak, i.e. Ŵ ≪ � 0 , and initial oscillator amplitude is sufficiently large, the NP amplitude does not change significantly over a few periods neither due to the damping nor due to the thermal heating (see Fig. 2f). In this regime, the averaged trajectories or phase space portraits follow very well deterministic nonlinear damped oscillations with negligible stochastic driving term, as described by Eqs. (2)(3)(4)(5)(6).  Figure 4b confirms the dependence of H on the initial conditions (x 0 , v x0 ) and shows the expected parabolic profile described by Eq. (6).
Determination of Ŵ : Fitting Eq. (9) to the variance Var(x) of trajectory ensembles starting at (x 0 , v x0 ) = (0, 0) ± 2 nm, 1 mm s −1 we obtained the damping coefficient Ŵ and the saturated-level pre-factor k B T/ m� 2 0 (see Fig. 3). Alternatively for higher ambient pressures we obtained the damping coefficient Ŵ , eigenfrequency 0 and bath temperature T by fitting Eq. (8) to the same dataset.
Determination of D : Fitting Eq. (3) to the mean position of trajectory ensembles x(t) over a few periods ( ∼ 10 ) we determined the oscillation frequency D (see Fig. 2f). The initial amplitudes x 0 may reach up to ∼ 2× the standard deviation of particle positions in the trap. For even bigger initial displacements the size of trajectory ensemble becomes very small and the results of analysis are unreliable.
Determination of 0 and ξ : Employing Eq. (5) with D and Ŵ determined above gives the eigenfrequency of a damped Duffing oscillator H . Fitting parabolic dependence of Eq. (6) to H for different initial conditions x 0 gives the eigenfrequency 0 and the coefficient of Duffing nonlinearity ξ (see Fig. 4b).
Determination of T : Determination of Ŵ in the first step above gave the pre-factor k B T/ m� 2 0 or directly the effective temperature T . Since 0 is known from the previous step, the pre-factor gives the effective temperature T (knowing the NP diameter 170 nm and its density 2000 kg/m 3 ).

Numerical solution of deterministic Duffing equation (DDE).
To check the accuracy of the approximative solutions given by Eqs. (3-6), we solve DDE given by Eq. (2) using the direct numeric integration (Matlab, ODE45 function) under the selected initial conditions. We use this procedure to search for parameters ( Ŵ , 0 , and ξ ) giving the best coincidence between the numerical solution of DDE and experimental data corresponding to a few nonlinear oscillations presented in Sect. 1 (see also an example in Fig. 2f).
Determination of Ŵ, 0 , and ξ : We fitted numerical solution of Eq.
(2) to all averaged trajectories x(t) (for various initial conditions (x 0 , 0) ) at once. The contribution of trajectories was weighted as x −2 0 so that all trajectories contributed to the residual sum with the same weight.

Power spectral density (PSD).
For comparison with other methods we also used this classical method.
PSD does not depend on initial conditions, therefore, it is not proper to characterize the transient dynamics. It can be used only to observe similarities and difference between the oscillator parameters at the short-time and www.nature.com/scientificreports/ long-time scale. Furthermore, this method does not allow to determine the coefficient of nonlinearity because only the harmonic potential is assumed to get Eq. (1). Determination of Ŵ and 0 : Fitting Eq. (1) to experimental P xx gives 0 and Ŵ. Determination of T : In the case that the calibration constant of the position detector is determined by other means we can use the velocity PSD P vv , see Eq. (11) and Methods, in order to determine the spectral effective temperature. In such a case the spectral temperature can be expressed as 29 For this approach the limiting factor is the determination of high frequency shot noise, see Methods, or the actual limits of the integral achievable from the experiment.

Numerical solution of stochastic Duffing equation (SDE).
In order to verify the described methods used for data processing and determination of parameters of the Duffing oscillator, we simulate the random particle motion in 1D described by the SDE given by Eq. (7) based on the Verlet scheme 58 . This approach is also generalized to 3D and corrected for influences of the other axes, as we described in more details in Supplementary Information.
We simulated 200 trajectories of an NP with random initial conditions with a time step corresponding to the experimental sampling frequency (total duration of a single trajectory was 1 s). The data were processed in the same way as experimental data and the obtained parameters were compared with the input ones.

Results and discussion
We have acquired data at several levels of ambient pressure, processed in a way described in the previous section for different initial conditions and determined oscillator's parameters using methods DOA, DDE and PSD. Further, we obtained trajectories using SDE by means of computer simulations for parameters determined from the experiment and processed in the same way as the measured trajectories, see also Supplementary Information. Eigenfrequency 0 . 0 obtained by DOA and DDE methods are about 1 kHz higher comparing to the PSD method. This is caused by the fact that 0 for the Duffing oscillator corresponds to the extrapolated oscillation frequency with zero amplitude x 0 which is the highest frequency, see Eq. (6). In contrast, the PSD method averages oscillation frequencies of all amplitudes.
Only for the highest pressure (10 mBar) values of 0 obtained by all methods are in agreement. In this case the dominant effect in the PSD method is the peak broadening by the damping and the resonant peak is not influenced by the nonlinear behaviour at all. The frequency drift through different pressures (same in both axes) is caused by the fluctuation of the total optical power in the trapping laser beam.
The analysis of the 1D simulated data follows the same trends as in the case of experimental data. In case of DDE method a weak bias ∼ 0.3 % towards higher frequencies appears. On the contrary the analysis of 3D simulations shows that values of eigenfrequency 0 obtained by all methods are shifted by a similar value towards lower frequencies, involving corrections described in Supplementary Information. Similar results were also confirmed by the analysis of SDE simulated data with various input conditions (see Supplementary Information). Therefore, we are persuaded the DOA and DDE methods provide valid 0 for weakly nonlinear systems of optically levitated particles.
The coefficient of Duffing nonlinearity ξ. DDE method gives slightly higher values of coefficient of nonlinearity than DOA with the mean value and standard error in the x and y axes equal to ξ x = (2.7 ± 0.2) µm −2 and ξ y = (2.1 ± 0.2) µm −2 , respectively. Discrepancy between both methods increases for ambient pressure above 2 mBar where values from DOA method drop down significantly. At these pressures the oscillations are strongly damped and the eigenfrequency H can not be determined with sufficient precision (see the yellow curve in Fig. 4b) where the fit by Eq. (6) fails to give reliable results.
If the DOA and DDE methods were applied on the simulated data (both 1D and 3D), the DDE method returns values corresponding to the input values of the simulations. The DOA method gives the same values as DDE but deviates in the same way as experimental data for higher pressures.
The values obtained with DOA and DDE can be compared with a theoretical estimate of ξ for the Gaussian beam with the beam waist w 0 under Rayleigh approximation 18, 44 using ξ = 2/w 2 0 that gives ξ = 2.76 µm −2 . The beam waist in the focal plane w 0 = 0.85 µm was determined from Zemax software package in accordance to the aberrations of the focusing lens at the trapping wavelength but ignoring polarization effects in such nonparaxial beam. Thus the obtained beam waist is the same for the x and y axes and approximately two times wider comparing to the direct utilization of numerical aperture of ideal focusing optics. Moreover, the particle diameter of 170 nm is slightly bigger than Rayleigh approximation allows, therefore the force profiles deviate from the expected simple gradient of optical intensity dependence which may also slightly influence the value of the coefficient of nonlinearity. Consequently the value of ξ predicted from the Gaussian beam waist corresponds quite well with the values determined from the experimental trajectories. Only the DOA method gives a very good coincidence with the theoretical model in the whole range of investigated pressures. DDE and PSD methods follow well the theory for pressures above 0.5 mBar however for lower pressures Ŵ values obtained by those two methods differ from the theoretical predictions by one order or even more. Since the peak in experimental PSD broadens and gets asymmetric at lower pressures due to the nonlinearity in the system (see Fig. 2d,e), its fit by a function derived for an ideal harmonic oscillator results in a larger Ŵ . In the case of 1D SDE simulated data the DDE method follows well the theoretical predictions (see red circles in row 3 of Fig. 5) while the analysis of 3D simulated trajectories by the DDE method (red crosses) reveals the same trend as the analysis of experimental data. On the contrary the DOA analysis of variance increase in the experimental data as well as in 1D or 3D SDE simulated trajectories leads to the almost perfect coincidence with the values predicted by the theoretical model. www.nature.com/scientificreports/ Thus we conclude that the increased damping obtained by DDE method is caused by the crosstalks between coordinate axes (especially between lateral and longitudinal ones).
Effective temperature T. The last parameter describing dynamics of levitated NP in a Duffing type potential is the effective temperature T. We compared the temperature obtained by the position variance in DOA method (referred to as transient temperature T TR ) and by the integration of the area under the oscillation peak in velocity PSD (PSDvv), referred to as the spectral temperature T SP

29
. We note we fixed the value of the calibration constant obtained at pressures higher than 0.4 mBar, as we explain in Methods. Both methods provide us with comparable results due to the fact that the energy of particle motion contained in PSDvv is independent of the nonlinear effects. We see that T is constant for pressures above 0.5 mBar and corresponds to the ambient temperature. For lower pressures the effective temperature firstly slowly increases which is followed by the steep increase for even lower pressures up to ∼ 500 K. This behaviour corresponds to the temperature increase observed in 57 and was explained by a decrease of the heat dissipation by conduction at low pressures. Unfortunately, our system does not allow us to reach even lower pressures ( < 10 −2 mBar) where one would expect the saturation of temperature since all absorbed heat would be dissipated only by the radiation independent of ambient pressure. The increase in the temperature is evident from the saturation level of the position variance shown in Fig. 3. For a higher pressure of 10 mBar (Fig. 3a) the saturation level is lower compared to the variance at pressure of 0.1 mBar shown in Fig. 3b. The temperature obtained by the analysis of the simulated trajectories corresponds well to the experimental values that were used as simulation inputs.

conclusions
Parameters of a nonlinear oscillator are usually estimated from steady states but the oscillator parameters found during the transient dynamics can be significantly different. The ability to thoroughly characterize the nonlinear oscillator on the short time scale is crucial for experimental testing of transient dynamics of levitating particles used in nanosensing and thermodynamical engines. The parameters are predicted even before the heating rate considerably affects the experiment which is desirable e.g. for quantum experiments with a prepared initial state. We introduced new transient methods to characterize parameters of an optically levitated NP behaving as a weakly nonlinear Duffing oscillator. The novel approach is based on post-processing of the acquired experimental data in such a way that for each selected initial state in the phase space, an ensemble of short evolutions is collected. From mean position and momentum and their variances, a Duffing oscillator can be fully characterized through its eigenfrequency, damping, coefficient of nonlinearity or temperature. The described methodology can be used also for very low temperatures down to quantum mechanical motion.
We developed two methods, i.e. Duffing oscillator approximation (DOA) and deterministic Duffing equation (DDE), and we applied them on experimental data and on datasets obtained from stochastic Duffing equation (SDE) for a wider range of experimentally accessible parameters. The comparisons between the methods and with the widely used steady-state method of power spectral density (PSD) assuming only an ideal harmonic oscillator are revealed in Fig. 5 and discussed in detail in the previous section. The comparison with SDE verified the reliability of parameters extracted by the methods for one-dimensional and three-dimensional motion of the NP. Focusing only on the parametric region corresponding to the experimental parameters we conclude that both DOA and DDE determine eigenfrequency 0 with a precision better than 1 % and nonlinearity ξ with a precision better than 14 % for pressures lower than 1 mBar. DOA gives damping factor Ŵ with a precision better than 1 % for all considered pressures, nonlinearities and temperatures. In contrast velocity PSD method gives temperature estimate within 1 % while DOA within 14 % for pressures between 0.01 and 10 mBar.
Moreover, the trajectories simulated by the SDE could be in principle used directly to obtain the parameters of the experimental system together with the concept of machine learning 59 . Here the artificial neural network (ANN) would be trained using the simulated trajectories or some features acquired of such trajectories and later the ANN would predict parameters of the experimental system based on the measured data. Yet this approach is extremely time consuming and computationally demanding and contradicts our approach in this paper to develop a simple and fast method to characterize oscillator parameters and optical trap properties. In conclusion, the presented transient methods can reliably characterize all important parameters describing the Duffing oscillator especially for lower pressures (under 1 mBar), where the nonlinearity plays a significant role and a peak profile in PSD starts to deviate from the ideal Lorentzian shape.

Methods
Calibration of quadrant photodetector. For the position calibration of the quadrant detector we used a method based on an integrating of the velocity PSD P vv 29 which was calculated directly from the position PSD P xx as P ∞ xx is the average level of the white (shot) noise at high frequencies in P xx . The calibration constant is then calculated as where T = 295 K is the room temperature. To determine the noise level P ∞ xx we analyzed P xx for frequencies ω/2π > 0.75f Nyq , where f Nyq = f sample /2 is the Nyquist frequency and the sampling frequency (11) P vv = P xx − P ∞ xx ω 2 . where A , B , and C are fitted parameters. If B > 0 we set the noise level as P ∞ xx ≡ C , otherwise we set the noise level to the mean value of the high frequency part ( ω/2π > 800 kHz) of P xx .
The recovered calibration factor C V→m forms a constant level plateau for pressures above 0.4 mBar (see Fig. 6) which is consistent with the previously published results 29 . For lower pressures below 0.4 mBar the calibration constant C V→m starts to drop. This can be attributed to the change of the NP temperature leading to a higher temperature of the NP center of mass motion than the room temperature used for calculation of the red curves in Fig. 6. Therefore, we applied the value of calibration constant equal to mean value of C V→m in the pressure range from 0.4 mBar up to 10 mBar to all studied pressures.
List of all possible methods for determination of parameters of the Duffing oscillator. In the following section we present an overview of several alternative ways that can be used for determination of parameters describing a Duffing oscillator. The mentioned methods are explained in more details in the main text.
1. Fitting Eq. (9) to the variance of trajectory ensembles Var(x) the damping coefficient Ŵ and the saturatedlevel pre-factor k B T/ m� 2 0 are obtained. 2. Fitting Eq. (8) to the variance of trajectory ensembles Var(x) for higher ambient pressure, the damping coefficient Ŵ , eigenfrequency 0 and effective temperature T are obtained. 3. Fitting Eq. (3) to the mean position of trajectory ensembles x(t) over a few periods, medium damping Ŵ and the oscillation frequency D are determined. 4. Fitting Eq. (1) to the position PSD gives Ŵ assuming an ideal harmonic oscillator. 5. Fitting Eq. (11) to the velocity PSDvv gives Ŵ assuming an ideal harmonic oscillator. 6. Fitting numerical solution of Eq. (2) to all averaged trajectories x(t) (for various initial conditions (x 0 , 0) ) at once gives values of Ŵ , 0 and ξ.
Determination of D .
8. Employing Eq. (5) with D and Ŵ determined above gives the eigenfrequency of a damped Duffing oscillator H . Fitting parabolic dependence of Eq. (6) to H for different initial conditions x 0 gives the eigenfrequency 0 and the coefficient of Duffing nonlinearity ξ . 9.
Step (4) gives also 0 .   Figure 6. Pressure dependence of calibration constant C V→m calculated from the PSD P vv (red) and a value of calibration constant that was actually used for the data processing (blue) for the x (left) and y (right) axes.