Coherent oscillations of a levitated birefringent microsphere in vacuum driven by nonconservative rotation-translation coupling

Coupling between the rotational and translational motions of an optically levitated microsphere leads to coherent oscillations.

It is polycrystalline, composed of highly oriented nanocrystals (20 − 30 nm in diameter), in such a way that the optical axes are laid in a hyperbolic distribution throughout the whole vaterite crystal [ Fig. S1(b)] [38]. As a result, vaterite is effectively a positive uniaxial birefringent crystal with the birefringence ∆n ≈ 0.1, which allows the exchange of spin angular momentum with circularly polarised light [29]. When vaterite is trapped with a linearly polarised beam, on the other hand, the mean optical axis of the crystal will align to the electric field (x axis), which is perpendicular to the propagation direction of the trapping beam (z axis) [Figs. S1(b,c)]. where the mean optical axis is in the x-direction. With a linearly polarised beam (z axis), the mean optical axis will align to the direction of E along the x axis. (c) Optical microscope image of a vaterite particle (2.2 µm in radius) levitated in vacuum with a linearly polarised beam, where the particle orients its optical axis (x-direction) to the E-field. Due to the birefringence, optical image of the particle may look like an oval shape.

S2. GAS PRESSURE AND VISCOSITY IN A RAREFIED GAS
Gas viscosity around a microsphere in a rarefied gas is not the same as the bulk viscosity but depends on the ratio of r/λ, where r is the radius of the particle and λ is the mean free path of air molecules [19]. To directly measure gas viscosity in a rarefied gas, we have trapped a vaterite microsphere of the same size (r = 2.2 µm) with a circularly polarised beam. Fig. S2 shows the rotation rate of the microsphere depending on the gas pressure from 1000 mBar (atmospheric) to 1 mBar (blue curve). Calibrating air viscosity of 18.1µPa s at atmospheric pressure, gas viscosity around the microsphere (red curve) can be determined by its rotation rate (blue curve), where we obtain the critical viscosity µ X = 0.2 µPa s at the critical pressure of P X = 1.18 mBar.

S3. NUMERICAL SIMULATIONS
Direct simulations are performed in 3d by numerical integration of the following Langevin equation, where q are the coordinates of the centre mass (x, y, z) and orientation of the rotationally symmetric particle. f L (t) is the uncorrelated Langevin force, with amplitude fixed by the fluctuation-dissipation theorem and a mean value of zero.
M the mass (m) and moment of inertia, (I). We use the integration scheme described in [39]. Although this numerical scheme applies only to translational motion, it can be trivially extended to include rotations. Instead of updating positions and velocities alone, we also update orientations and angular velocities. As usual, this is achieved through successive multiplication of rotation matrices (See [40] and references therein). Optical calculations are performed with T-matrix theory [41] and the forces are calculated through integrals of the optical momentum flux, given by the Maxwell stress tensor, through a closed surface surrounding the particle. In particular, optical fields are expanded in vector spherical harmonics labelled by eigenvalue pairs (n, m) with |m| < n. Expansions are truncated at n = N max . The particle itself is modelled as a homogeneous sphere with uniform, birefringent refractive indices [42]. Although this ignores the internal structure of real vaterite particles, it is, nonetheless, a birefringent particle with the same overall symmetry and should be expected to behave in a similar way to a real particle. Numerical simulations are performed with a sphere of nominal parameters corresponding to experimental conditions, i.e. the radius is 2.2 µm, and the ordinary and extraordinary refractive indices are those of bulk vaterite, i.e. (n e = 1.65, n o = 1.55). The modelled beam is a linearly polarised Gaussian beam, rendered by the Richards and Wolf formulation [43], with the same parameters as those used in the experiment. An optical power (i.e. integral of the normal component of the Poynting vector over a plane normal to the beam axis) of 5 mW is used for all simulations. As the simulation runs, a Brownian trail is collected. From this trail we can evaluate power spectral densities and time dependent covariances. Table S1 gives the default parameters of the simulation. Simulations of this sort take ≈ 10 hours on a typical desk top computer. Fig. S3 describes the motion of the centre of mass of the particle. Three conditions are considered, (i) a birefringent particle ∆n = 0.1 at low viscosity µ = 1 × 10 −6 Pa s, (ii) the same birefringent particle at a higher viscosity, µ = 3.25 × 10 −6 Pa s, and (iii) an optically isotropic particle, ∆n = 0 at the lower viscosity, µ = 1 × 10 −6 Pa s. The top left panel shows a scatter plot of the centre of mass in the transverse (xy) plane. The top right shows the power spectral density of the x coordinate of the centre of mass. The lower panel shows the autocorrelation of the x coordinate. As can be seen, cases (ii) and (iii) show similar behaviour; the centre of mass is symmetrically confined by the trap and the power spectral densities show no particularly sharp peaks. The autocorrelation of the x coordinate thermalises quickly. The behaviour for case (i) is dramatically different. As with the experiment, the x coordinate varies much more widely, ± ∼ 1 µm, while the y coordinate remains tightly confined. Very sharp peaks appear in the power spectrum, the strongest (by a factor of ∼ 10 6 ) close to the resonant frequency of the trap. In addition the oscillation of the x coordinate becomes highly coherent, with very slow decay of the autocorrelation function. These findings closely mimic those observed in the experiment, as discussed above and in the main text.
The simulations allow us to understand the rotational motion of the particle. A strong interaction between the position (x) and pitch (θ) of the birefringent sphere is observed, see Fig. S4. The scatter plot of the (x, θ) coordinates, combined with the cross correlation reveal the mechanism responsible for the formation of the high amplitude coherent oscillations. Vibrations in x couple constructively with vibrations in θ and become self-sustaining. In particular, angular oscillations generate force components in the x direction, that drive the x oscillation. Similarly, the oscillation in x generates torques about the y axis, driving the angular oscillation. Thus, angular and translational oscillations feed one another, passing energy from the optical field to the particle motion. At steady state, this accumulated energy is balanced by viscous dissipation into the gas. The resulting, self-sustaining oscillation is very clearly illustrated in the accompanying supplementary movie 2.

S4. THEORY FOR LINEARISED FORCES WITH ROTATION-TRANSLATION COUPLING
For small displacements, all of the degrees of freedom are uncoupled, except x and θ. This follows from basic symmetry arguments [8,42,44,45], and is confirmed by the T-matrix calculations. Uncoupled degrees of freedom are trivial and, for completeness, treated in the final section below. In summary, uncoupled degrees of freedom are effectively one dimensional and, since any 1D force field can be derived from a potential, they show the usual equilibrium behaviour with the thermal and elastic energies being equal, and the variance of the coordinates independent of viscosity. The anomalous behaviour discussed in this paper comes from the coupling between x and θ. Essential features of the system can be captured by restricting attention to these coupled coordinates.
In the following we consider a linearised Langevin equation for the coupled degrees of freedom of the vaterite . The right hand panel shows the cross correlations between x and θ.
particle. As discussed below, uncoupled coordinates are effectively at equilibrium and do not play a part in the effects of interest.
where K is a stiffness matrix given by linearising the external forces (f ext = f opt − mgẑ) about the trapping configuration i.e. K = −[∇f ext (q eq )], Ξ the hydrodynamic friction, M the mass (m) and moment of inertia, (I), and q are the coupled coordinates.
Here, the external (i.e. optical and gravitational) forces and torques have been linearly approximated relative to a position and orientation at which they vanish and are locally restoring. We find the stiffness coefficients by orienting the particle with its optic axis parallel to the x axis, and scanning its centre of mass along the z axis in increments of 0.1 µm. At each stage we calculate the the z component of the force. When it changes sign, we have an interval within which we perform a golden search for the point at which the forces and torques vanish. The stiffness coefficients are then found by numerically calculating the gradient of the forces about the trapping point. For positively birefringent spheres the optic axis aligns with the electric polarisation and the centre of mass lies on the beam axis so that the upward optical scattering forces are equal and opposite to the downward weight of the particle. Optical forces are, in general, non-conservative [6]. As a result K is generally non-symmetric (i.e. K xθ = K θx ) [8]. f L (t) is the uncorrelated Langevin force, with amplitude fixed by the fluctuation-dissipation theorem and a mean value of zero.
In the frequency domain, this becomes, Equivalently, where, (S11) The eigenfrequencies of the matrix A in Eq. (S9) determine the stability of the trap and statistical measures of the stochastic motion of the particle. Below, we describe some general properties of these eigenfrequencies, their dependence on the ambient viscosity and the consequences for the linear stability, power spectral density (PSD) and time correlations of the trapped particle.

A. Eigenfrequencies and stability
The eigenfrequencies of A, ω i are given by the condition |A(ω)| = 0, a quartic in ω, where, Linear stability requires that (ω i ) > 0 for each ω i . Since P (ω) is a quartic, succinct expressions for the ω i are not generally available. However, by examining two extreme conditions, the overdamped limit (where inertia is negligible) and the underdamped limit (where viscosity is negligible), we can anticipate situations in which one of the (ω i ) will become negative as the ambient viscosity is decreased. Below, we derive this condition and give an expression for the viscosity at which (ω i ) = 0, and associated value of the real part of the frequency.

Overdamped limit, inertia negligible
In the low Reynolds number regime, inertia is negligible and Eq. (S12) becomes: which has roots, with, In order for the trap to be unstable in this regime, we require (ω i ) < 0. From Eq. (S15) this means, which reduces to, This is a hard condition to satisfy. The right hand side is positive (if the forces are restoring), which requires comparatively large coupling coefficients of the same sign. In the course of the parametric survey, presented below, this situation never arises. It may be fundamentally impossible for coupling forces and torques to exceed restoring forces and torques in amplitude.
2. Underdamped, low pressure limit, viscosity negligible In low pressure limit, Ξ → 0 and Eq. (S12) becomes, In this case, the roots of P (ω) = 0 are, with, Owing to the inclusion of inertia, the number of eigenfrequencies is twice that in the overdamped limit. Without x/θ coupling, we get: i.e. the spring vibrates at its natural frequencies. In the underdamped limit, (ω i ) < 0 for at least one of the ω i whenever ∆ LP < 0. Compared to the overdamped regime, this condition is relatively easily satisfied. For example, if the resonant frequencies for the x and θ oscillations are close, the bracketed term approaches zero and the stability condition reduces to K xθ K θx < 0. In this case, the size of the coupling coefficients is immaterial, it is only necessary that they are of opposite sign.

Critical viscosity giving (ωi) = 0
When the trap is stable with respect to rotation-translation coupling in the overdamped limit, but unstable in the underdamped limit then the value of at least one of the (ω i ) must change from positive to negative as the viscosity is decreased. The critical viscosity, for which (ω i ) = 0, can be determined, by requiring that ω i = (ω i ) = ω r is purely real in Eq. (S12). Under this condition, the real and imaginary parts of P (ω) must both be zero, Note that, according to Eqs. (S12,S13), when K xθ = K θx = 0, the condition (ω) = 0 requires ξxx m ω r = 0 or ξ θθ I ω r = 0 which when ξ xx = 0 or ξ θθ = 0, which can only be satisfied when the viscosity and, therefore, the pressure are zero. When K xθ = 0 and K θx = 0, Eqs. (S23a,b) must hold simultaneously: Eq. (S23b) gives, where ω X is the real frequency at which the imaginary component vanishes. Since both ξ xx and ξ θθ are directly proportional to the viscosity, i.e. ξ xx = 6πµa ≡ s x µ, ξ θθ = 8πµa 3 ≡ s θ µ, This expression can be substituted into Eq. (S23a) to give a relationship for the corresponding critical viscosity, µ X . rearranging, where µ X is the critical viscosity for which one of the (ω i ) = 0. When K xθ = K θx = 0 this expression gives a negative µ 2 X . The product K xθ K θx must be negative enough to make µ 2 positive. Interestingly, ω X = (ω i ) is completely independent of viscosity and given as a function of the stiffness coefficients, K xx , K θθ , the radius of the particle, its mass and moment of inertia. For a fixed field and particle, optical forces are directly proportional to the optical power. It may seem appear that ω 2 X and µ X are directly proportional to the beam power however this is only approximately correct. The trapping configuration of the particle is determined by a balance between the downward weight and the upward scattering force. As the power in the beam changes, so does the height of the trapping configuration. For this reason the stiffness in the trapping configuration is not directly proportional to the beam power.
Since ω X and µ X can be simply calculated from the stiffness parameters at the equilibrium trapping configuration, we can efficiently assess their variation with respect to the parameters defining the trap. Fig. S5 shows calculated values of ω X and µ X as functions of anisotropy, ∆n = n e − n o for varying beam power, sphere radius, mean refractive index and beam numerical aperture. The mean refractive index, n 2 = (n 2 e + 2n 2 o )/3, is held constant, unless its variation is being considered.
In general, µ X increases rapidly with ∆n, providing that ∆ LP < 0, and varies in the range 0 →≈ 5 × 10 −6 Pa s. ω X can vary from around 400 Hz to 1500 Hz within the range of parameters considered here. The strongest dependence of ω X is with the beam power.

B. Power spectral density
The power spectral density (PSD) for the coupled system can be expressed in terms of the eigenfrequencies, ω i When the preceeding conditions are satisfied, and one of the (ω i ) → 0 as µ → µ X , then the corresponding motional mode (or eigenvector of the motion) begins to dominate. Measures of the stochastic motion of the particle can then be  expressed in terms of this mode alone. Use this definition from now on Power spectra: from above, where, and A is defined in Eq. (S9). The product, CΞC H , is, where P x , P θ are the polynomials given in Eq. (S13). The denominator of Eq. (S28), P (ω)P * (ω), determines the general form of the PSD and can be factorized as follows. First, the quartic, P (ω), can be expressed as, P (ω) = i n a n ω n , with real coefficients a n . It follows immediately that, if ω i is a root of P (ω), then −ω * i is also a root. For example, if P (ω i ) = 0 then, P (−ω * i ) = i n a n (−ω * ) n = (−i) n a n (ω * i ) n ) = i n a n ω n i * = 0. Similarly, if ω i is a root of P (ω), then −ω i is a root of P * (ω). The roots of P (ω)P * (ω) therefore have the following form: Hence, we can write P (ω)P * (ω) as follows, It follows directly that P (ω)P * (ω) = P (−ω * )P * (−ω * ). We note that similar arguments can be made to show that C(ω)ΞC H (ω) * = C(−ω * )ΞC H (−ω * ), as can be confirmed directly from Eqns. (S13,S30).
The form of P (ω)P * (ω) fixes the behaviour of the power spectrum for the coupled system. If (ω i ) → 0 as the viscosity, µ → µ X , then the denominator, P (ω)P * (ω) becomes singular as ω → ω X . In particular, as µ → µ X , ω i ≈ ω X + iδ, for small δ, and (ω − ω 2 i ) ≈ 2ω X δi. As µ → µ X , and ω → ω X , Qualitatively, as the critical viscosity is approached, the power spectrum is increasingly dominated by a sharp peak at ω = ω X , giving an oscillation with a high quality factor at a finite viscosity. As discussed below (section S4 D, Eq. (S41), the behaviour for uncoupled or conservative systems, is similar except that the condition required for the imaginary eigenfrequencies to approach zero is µ → 0. For coupled systems, a similar limit is reached for finite viscosities, µ X . More fundamental differences between coupled (non-conservative systems) and uncoupled (effectively conservative systems) are revealed in the correlations and variances of the particles, discussed immediately below.

C. Time correlations
The time correlations follow from the Weiner-Khinchine theorem: where Q(ω)⊗Q * (ω) given by Eq. (S28). The integral in Eq. (S33) can be evaluated by applying the residue theorem to a contour closed by a semi-circle in the lower half plane for τ > 0 i.e.
where ω j are the poles in the lower half plane. If ω 1 is one such pole, the corresponding residue is, .
Equations (S35,S37), show the key point. The covariance matrix, q(t) ⊗ q(t + τ ) is the sum of two damped oscillations with amplitudes proportional to the reciprocal of the imaginary parts of the eigenfrequencies of the system. As described above, in section S4 A, (ω i ) → 0, under particular conditions. When these conditions pertain, the time dependent covariance of the coordinates becomes increasingly dominated by this mode as µ → µ X . Both the amplitude (the instantaneous covariance) of this oscillation and its decay time become large. In particular, if ω ≈ ω X + iδ as µ → µ X , we have, To clarify the meaning of Eq. (S38), we note that q(t) ⊗ q(t + τ ) is, of course, a symmetric matrix. This matrix consists of the sum of the residues detailed above. Each term in this sum is weighted by a scalar whose denominator contains a term (ω i ), and whose time variation if e −i (ωi)τ e − (ωi)τ , for characteristic frequencies ω i with poles in the lower half plane. As µ → µ X , this sum is increasingly dominated by the terms for which (ω i ) tends to zero. The right hand side of Eq. (S38) shows the limiting form of the scalar function multiplying the well behaved matrix term. The equation shows how each component of the time dependent covariance matrix grows as µ → µ X . This result is fundamentally different from the familiar situation where coordinates are uncoupled (see final section, S4 D) Eq. (S48)). In that case, the amplitude of the time dependent covariance (e.g. x 2 ) is a constant, independent of viscosity, given by the equipartition theorem e.g. 1 2 k B T = 1 2 K xx x 2 [12]. The decay times of the covariance, meanwhile, also increase with decreasing viscosity but they do so much more slowly, and they do not approach zero for any finite viscosity.
We note that the result, x 2 → ∞, with µ → µ X , obtained for the coupled system, does not indicate that the particle amplitude is becomes unbounded in reality. Physically, the particle makes excursions of increasingly large amplitude until the linear range of the trap is exceeded. Beyond this limit, the motion becomes erratic, or even chaotic. In the experiment, the particle may simply be lost.
In the next section we show that the treatment for the coupled system reduces to the familiar results for uncoupled coordinates when K xθ = K θx = 0.

D. Specialisation to uncoupled (1D) systems
Finally, we consider the more familiar case of uncoupled (i.e. one dimensional) systems. In contrast to the coupled systems described above, a one dimensional force field is necessarily conservative. The most conspicuous consequence of this, is that the instantaneous variance of the coordinates, x 2 is independent of viscosity and given by the equipartition theorem, i.e. 1 2 k B T = 1 2 K xx x 2 . Other qualitative differences appear in the PSD and the time dependent covariance. In both the coupled and uncoupled cases, the peak in the PSD is ∝ 1/ (ω i ) and the decay of the time dependent covariance is ∝ exp( (ω i )τ ). However, for the uncoupled case, (ω i ) ∝ ξ xx , which remains finite for finite viscosity, in contrast to the coupled case for which (ω i ) → 0 as the viscosity approaches its critical value, µ X . In the uncoupled case, therefore, the peak in the PSD has a limited height and the time dependent covariance decays more rapidly [Eq. (S48), below], in comparison to the coupled system [Eq. (S38)]. In the following, we show that these familiar results follow from the preceding treatment of the coupled system, when the coupling coefficients are set to zero, i.e. when K xθ = K θx = 0.
First of all, the eigenfrequencies are given by the roots of P (ω), Eq. (S12). When K xθ = K θx = 0, so the roots are given by P x (ω) = 0 and P θ (ω) = 0. Restricting attention to the x coordinate, θ being similar, the roots of P x (ω) = 0 are, Notably, (ω 1 ) = ξ xx /2m, which remains finite for finite viscosity, in contrast to the coupled case. The PSD, given by Eq. (S28), for the x coordinate is, when K xθ = K θx = 0. Factorizing P x P * x gives the PSD, As the viscosity is reduced, ω 1 becomes, where ω 0 is the natural, resonant frequency. Then, the PSD peak height is, which becomes large as µ → 0, but much less rapidly than the coupled case, which becomes large as the viscosity approaches a finite, critical value namely µ X . In particular the peak height of the PSD ∝ 1/µ for the uncoupled case, and ∝ 1/(µ − µ X ) 2 for the vaterite spheres discussed in this article [see Eq. (S32)]. The time dependent covariance, x(t)x(t + τ ) , depends on the residues of the PSD, through Eq. (S34). Here, the poles in the lower half plane are the roots of P * x , i.e. −ω 1 , ω * 1 . The residue at −ω 1 is, = k B T m iω 1 ξ xx |ω 1 | 2 (ω 1 ) (ω 1 ) .
Summing the residues with Eqs. (S36, S37), and combining with the chracteristic frequency, Eq. (S40) the time correlations, [Eq. (S34)] gives, The final expression uses the approximation to the characteristic frequency, Eq. (S44). In terms of algebra, the reason that the instantaneous variance, x 2 , is independent of viscosity is that the residues introduce a factor (ω i ) in the denominator which, in this case is ξ xx [Eq. (S40)]. This factor precisely cancels the term ξ xx that appears in the numerator due to the normalization of the stochastic forces [Eq. (S3)] which itself is a consequence of the fluctuationdissipation theorem. This cancellation does not generally occur when degrees of freedom are coupled; in that case we have (ω i ) → 0 at finite viscosities, which cannot cancel the ξ xx in the numerator which, of course, is finite for finite viscosity. Physically, the qualitative difference between the two forms of behaviour rests in the fact that the coupled system is generally non-conservative, and the uncoupled system is effectively one dimensional and, therefore, effectively conservative.

S5. ADDITIONAL SUPPLEMENTARY MATERIALS
Supplementary movie 1 shows a vaterite microsphere (2.2 µm in radius) oscillating along the polarisation direction with an amplitude ≤ 1 µm at a gas pressure of 1.2 mBar. The frame rate is rendered at 15 fps from 5000 fps.
Supplementary movie 2 shows the simulated stochastic motion of optically trapped spheres. The spheres are shown full size. The double arrow bisecting each sphere shows the direction of the symmetry axis each particle. In the case of the isotropic particle, this axis is arbitrarily fixed in the sphere. The parameters used correspond to those