3D sympathetic cooling and detection of levitated nanoparticles

Cooling the center-of-mass motion of levitated nanoparticles provides a route to quantum experiments at mesoscopic scales. Here we demonstrate three-dimensional sympathetic cooling and detection of the center-of-mass motion of a levitated silica nanoparticle. The nanoparticle is electrostatically coupled to a feedback-cooled particle while both particles are trapped in the same Paul trap. We identify two regimes, based on the strength of the cooling: in the first regime, the sympathetically cooled particle thermalizes with the directly cooled one, while in the second regime, the sympathetically cooled particle reaches a minimum temperature. This result provides a route to efficiently cool and detect particles that cannot be illuminated with strong laser light, such as absorptive particles, and paves the way for controlling the motion of arrays of several trapped nanoparticles.

Levitated particles are a promising experimental platform for testing fundamental physics and detecting weak forces [1][2][3]. Having recently entered the quantum regime [4][5][6][7], levitated particles offer a means to study quantum mechanics at mesoscopic scales and to build sensors with quantum-enhanced sensitivity. In particular, it has been proposed to use levitated particles with internal degrees of freedom, such as magnetic particles or nanodiamonds with color centers, to create macroscopic superposition states [8,9] and to study spin-mechanical coupling [10]. However, optical trapping of such absorptive particles in high and ultra-high vacuum is challenging due to the increase in their bulk temperature, which can lead to burning and graphitization [11][12][13]. Moreover, in the quantum regime, elevated internal temperature is associated with the decoherence of the centerof-mass (CoM) motion [14]. Electrodynamic levitation removes the trapping laser field as a heating mechanism and has already been employed to trap absorptive particles [12,[15][16][17]. Nonetheless, state-of-the-art experiments with electrodynamic traps still rely on laser fields to cool [18,19] and interrogate the particle [15,20].
Sympathetic methods provide a solution: if an absorptive particle is electrostatically coupled to one or more transparent particles, cooling and detection of the absorptive particle can be achieved when only the transparent particles are exposed to laser light; the absorptive particle remains in the dark. Sympathetic cooling has been realized in experiments with neutral atoms [21], atomic ions [22], molecular ions [23], positrons [24] and proteins [25], while sympathetic detection has made possible indirect observation of ions in a Paul trap [23,26]. Recently, one-dimensional sympathetic cooling has been demonstrated with two levitated micro-and nanoparticles coupled via electrostatic [27,28] or light forces [29]. However, stabilization of the motion in the two remaining dimensions is still necessary to counteract trapping instabilities [4][5][6][7]. Moreover, in order to understand whether sympathetic cooling will be sufficient to reach the quantum regime with absorptive particles, the limitations of the technique need to be studied. Finally, sympathetic detection has not previously been demonstrated for mesoscopic particles.
Here we report on sympathetic cooling and detection along all three spatial axes in a system of two levitated mesoscopic particles. We investigate different regimes of the feedback cooling strength and identify the optimum strategy to achieve the minimum temperature of the sympathetically cooled particle. Finally, we reconstruct the motion of a levitated particle via interaction with an auxiliary particle present in the trap. In this proof-ofprinciple demonstration, both particles are transparent to laser radiation, but it is possible to extend the same techniques to particles with arbitrary optical properties. Optimal detection and cooling along all three axes of motion thus paves the way for experiments with absorptive particles in high and ultra-high vacuum and provides a tool to control the motion of arrays of several particles, which could be used in the search for dark matter [30][31][32] and in investigations of quantum gravity [33].
A schematic of the experimental setup is shown in Fig. 1a. We levitate a pair of silica spheres 300 nm in diameter (Bangs Laboratories, Inc.) in a linear Paul trap at pressures between 1 × 10 −3 mbar and 1 × 10 −2 mbar. Each particle carries 950(50) elementary charges [34]. The distance between opposing radiofrequency (RF) electrodes is 2r 0 = 3 mm, while the distance between the endcap electrodes is 2z 0 = 7.9 mm. We drive the trap with a peak-to-peak voltage of 470 V at 10.1 kHz and a DC offset on the RF electrodes of 11.5 V. The DC voltage on the endcaps is varied between 40 V and 100 V, depending on the experiment. The endcap voltage sets the distance between the two particles and thus the strength of their electrostatic interaction.
In order to verify sympathetic cooling and detection, we measure the position of each particle independently via optical detection. The endcap electrodes are hollow cylinders, through which a weakly focused 100 mW laser beam at 780 nm illuminates the particles (Fig. 1a) IG. 1. Trapping and detection of two particles. (a) Schematic overview of the experimental setup. The particles are illuminated with a weakly focused laser beam at 780 nm (detection laser). Light from this laser scattered by the particles is collected with aspheric lenses, coupled to optical fibers, and used to detect the individual motion of each nanoparticle at photodiodes APD1 and APD2. A laser beam at 1064 nm (feedback laser) applies an optical feedback force to particle P1. Inset: trap axes and the alignment of the detection system with respect to the particles. (b) An example of the power spectral density (PSD) of the center-of-mass motion of particle P2. Each peak corresponds to one of six collective oscillation modes, illustrated above.
lenses and guided to fiber collimators. One lens collects light from the top particle, P1, while the second lens collects light from the bottom particle, P2; each lens and collimator image their particle on the end facet of a fiber. The power coupled into each fiber depends on the particle position [35,36] and is measured with avalanche photodiodes (Thorlabs APD440A). A linear relationship between the position of each particle and the photodiode signal is achieved by displacing the image to one side of the fiber core. Cross-talk between the detection channels is below Weak cooling regime. Temperature of the sympathetically cooled particle P2 as a function of the temperature of the directly cooled particle P1 for (a) radial modes xout, xin, yout, yin and (b) axial modes zin, zout. The line Tsymp = T direct is plotted as a guide to the eye. Error bars correspond to one standard deviation and are determined from the fit. (c) Power spectral densities (PSDs) of the motion of P2 at room temperature (blue curve) and with feedback cooling applied simultaneously to all six collective oscillation modes (orange curve). When feedback cooling is applied, the mean temperature of the six modes is 22 K. the detection noise. We choose this detection scheme for its simplicity; we find that it has roughly the same sensitivity (2 × 10 −10 m/ √ Hz) as that of the standard scheme based on the interference of forward-scattered light and a Gaussian reference beam, although we have not made a detailed comparison of the two schemes' efficiencies. From the equations of motion of two coupled harmonic oscillators, one finds two eigenmodes of oscillation along each spatial axis {x, y, z}, corresponding to the particles' collective in-phase and out-of-phase motion [34]. Figure 1b shows an example of the power spectral density (PSD) of particle P2's centerof-mass motion, in which each peak corresponds to an eigenmode.
A second laser beam at 1064 nm is coupled through a dichroic mirror and focused onto particle P1. The laser power is modulated with an acousto-optic modulator (not shown) controlled via a field-programmable gate array (FPGA). This configuration allows us to apply optical cold damping to P1 [18].
Having reviewed the setup, we now present experimental results on sympathetic cooling of particle P2. All data presented in the current manuscript were ob-tained with the same pair of nanoparticles; the nanoparticles have the same size and the same charge [34]. We start with an analysis of the weak cooling regime. In this case, the rate γ fb,q of feedback cooling applied to particle P1 along axis q is much smaller than the frequency splitting between the collective oscillation modes: Here q ∈ {x, y, z}, j q is a coupling constant and ω q,0 is the resonance frequency of a single particle in the Paul trap along axis q [34]. We study the cooling of each mode individually. The motional temperature of P1 is set via the strength of the cold damping, which is controlled by the feedback gain on the FPGA. For each setting of the feedback gain, we determine the temperatures of both P1 (directly cooled) and P2 (sympathetically cooled) by fitting the PSD with a Lorentzian function [34]. Figures 2a and 2b show the dependence of the temperature of P2 on the temperature of P1 for each collective mode. In this regime, we find that the temperatures of the effective CoM motion of the two particles are equal, that is, the temperature of P1 is fully transferred to P2. This result is consistent with the transfer of energy between the particles being significantly faster than the cooling rate of the directly cooled particle. The minimum temperatures shown in Fig. 2 are determined by the feedback gain set via the FPGA and by the background gas pressure. We choose to work at 2.7 × 10 −3 mbar for convenience: this pressure makes the trapping robust against accidental excitation of the motion (and subsequent loss of the particles) while still allowing us to work in the weak cooling regime. There are no intrinsic limitations to extending sympathetic cooling to lower pressures [28].
It is possible to sympathetically cool not only individual modes but also all six degrees of freedom simultaneously. To this end, we apply cold damping via six FPGAs, each configured to address a single mode. The six output signals of the FPGAs are mixed and used to modulate the power of the 1064 nm laser. Figure 2c compares the PSDs of the sympathetically cooled particle at room temperature and under the influence of feedback cooling. The room-temperature PSD is measured at 2.7 × 10 −2 mbar, while the PSD of the cooled particle is obtained at 2.3 × 10 −3 mbar. (We intended to take both data sets at the same pressure, but due to a citywide power failure, the particles were lost from the trap before we could do so.) The mean CoM temperature of the particle under feedback cooling is 22 K; just as for the individual modes, it is determined by the feedback gain and the background gas pressure.
Next, we examine the regime of strong cooling. In this regime, the cooling rate is comparable to or larger than the normal-mode frequency difference: γ fb,q /2π |ν in q − ν out q | ≈ j q /2πω 0 . We consider motion along the x axis. By reducing the endcap voltage and thus increasing the particle separation, we reduce the mode splitting to ν in x − ν out x ≈ 70 Hz; we then vary the cooling rate from 2 Hz to Error bars correspond to one standard deviation and are determined from the uncertainty in the calculation of the area under the velocity PSD curves [34]. The thickness of the theory lines expresses the temperature uncertainty due to the estimation of the gas damping parameter. 220 Hz and record PSDs of both particles' motion. The cooling rate is determined as follows: In the weak cooling regime, fits of the mode linewidths allow us to extract cooling rates for a range of electronic gains [37]. We then have a known relationship between the electronic gain and the cooling rate that also applies in the strong cooling regime.
The PSD of the sympathetically cooled particle P2 as a function of the cooling rate is shown in Figs. 3a and 3b. For weak cooling, we observe two distinct peaks of the resonant motion. As the cooling rate is increased into the strong cooling regime, the peaks broaden and merge into a single peak, the height of which increases. A detailed theoretical analysis based on the model of coupled harmonic oscillators [34] shows that the two modes are still present. However, only one of them is not overdamped and thus is still visible in the PSD. We extend the model by coupling the oscillators to a thermal bath [34] and calculate PSDs for different feedback cooling rates. The resulting plot is shown in Fig. 3c. We find good agreement between the measured and the calculated evolution of the PSD, suggesting that the model captures the essential features of the cooling process.
With this understanding of the feedback cooling pro- cess in hand, we analyze the particles' temperatures as a function of the cooling rate. Since the two individual modes cannot be resolved in the strong cooling regime, we calculate the effective temperature of motion along the x axis from the area under the velocity PSD [34]. Figure 3d shows the extracted temperatures as a function of the cooling rate applied to the directly cooled particle P1. For cooling rates below the frequency of the normal-mode splitting, the two particles have equal temperatures, similar to the data shown in Fig. 1. As the cooling rate increases, the temperature of P1 continues to decrease monotonically, while the temperature of P2 reaches a minimum. We also plot the temperatures extracted from the calculated PSDs of Fig. 3c, which are in agreement with the experimental data.
To gain intuition about the origin of the temperature minimum under sympathetic cooling, one can consider the efficiency of the energy exchange between the two particles as coupled harmonic oscillators. First, this efficiency is limited by the coupling strength. In the strong cooling regime, we extract energy from P1 faster than it flows from P2 to P1. At this point, increasing the cooling rate further does not lower the temperature of P2. Second, two coupled oscillators exchange energy most efficiently when their bare resonance frequencies, i.e., in the absence of coupling, are equal. Damping applied to P1 shifts the particle's bare resonance to a lower frequency than that of P2, making the energy exchange inefficient. As a result, in the strong cooling regime, the flow of energy from P2 to P1 is reduced.
Not only can an auxiliary particle in an ion trap be exploited for sympathetic cooling, but also it can be used to detect the motion of the particle of interest, which is especially useful when that motion cannot be read out directly. We implement sympathetic detection, for which the roles of particles P1 and P2 are reversed: P2 is used to detect the motion of P1. We excite the motion of P1 by applying a sinusoidal optical force to it via the 1064 nm laser beam. No external driving force is applied to P2. We scan the frequency of the force and measure the amplitude and phase response of both particles, which are shown in Fig. 4. The amplitude response of P2 to the external drive of P1 allows us to calculate a transfer function Θ (ω) that connects the motion of P1 and P2 [34]. From the transfer function, we reconstruct the amplitude and the phase response of P1 from the measurements of P2 via where {A P1 , A P2 } and {φ P1 , φ P2 } are the amplitude and phase responses of the particles. The reconstructed response functions A P1 and φ P1 , plotted in Fig. 4, are in excellent agreement with the independently measured values of P1, establishing that sympathetic detection can be used to reconstruct the motion of P1. The dip between the in-phase and out-of-phase modes in the amplitude response of P1 is due to destructive interference of the motion in those modes when one particle is driven, an effect analogous to electromagnetically induced transparency [38].
In conclusion, we have demonstrated both sympathetic cooling and detection of nanoparticles levitated in a Paul trap and coupled via electrostatic interaction. If the cooling rate is smaller than the normal-mode frequency difference, the temperature of the sympathetically cooled particles equals the temperature of the directly cooled particle. For higher cooling rates, the temperature of the sympathetically cooled particle reaches a minimum value due to reduced efficiency of the energy exchange between two coupled harmonic oscillators. Both methods pave the way for experiments with absorptive particles in high and ultra-high vacuum in the absence of light. For example, a nanodiamond with nitrogen vacancies could be co-trapped with a silica nanoparticle; both particles' charge and mass could be calibrated at low vacuum, where the diamond is not significantly heated by laser light; subsequently, at ultra-high vacuum, the silica nanoparticle could be used to measure and control the motion of the nanodiamond. Furthermore, the technique provides a tool to detect and control the 3N modes of collective oscillations of a string of N particles by acting on just a single particle in the chain. Building on the spinmechanical coupling recently demonstrated with a single levitated diamond hosting nitrogen-vacancy defects [39], our work opens up new experimental possibilities for real-izing a quantum spin transducer based on levitated particles. Such particles could be used as long-lived quantum memories for quantum information applications [40,41] or for force sensing below the standard quantum limit in the search for new physics [30][31][32]. In this section, the motion of nanoparticles is analyzed in two cases: free evolution and driven oscillations. We model two levitated and electrostatically interacting nanoparticles as coupled harmonic oscillators. We assume that feedback cooling is applied to the first particle, referred to in the main text as P1 and here denoted with 1 in subscript. For displacements that are small compared to the distance between the particles, the equations of motion can be written as where q ∈ {x, y, z} is the axis of motion, and the particles at rest lie along the z axis; m 1 and m 2 are the masses of the nanoparticles; F q,1 and F q,2 are the external forces acting on the first and second particle, respectively; ω 0q,1 and ω 0q,2 are the resonance frequencies of each particle's center-of-mass (CoM) motion; γ 01 and γ 02 are the rates at which the particles are damped by background gas, which depend on their radii and mass; γ fb,q is the cooling rate; µ = m 1 /m 2 is the mass ratio; and j q is a coupling constant. For axes x and y, the coupling constant is j x = j y = Q 1 Q 2 / 4π 0 d 3 m 1 ; for axis z, the coupling constant is j z = −2Q 1 Q 2 /4π 0 d 3 m 1 , where Q 1 and Q 2 are the particles' charges and d is the distance between the particles. The resonance frequencies ω 0q of a charged particle in a linear Paul trap can be calculated as [1]  with the a and q parameters defined as where Ω d is the frequency of the AC voltage applied to the trap radiofrequency (RF) electrodes, V pp is the peak-topeak amplitude of the AC voltage, V off is the DC offset applied to the RF electrodes, V end is the DC voltage applied to the trap endcap electrodes, κ RF and κ end are geometrical factors, r 0 is half the distance between opposing RF electrodes, and z 0 is half the distance between the endcap electrodes. Our parameters are V pp = 470 V, V off = 11.6 V, κ RF = 0.77, κ end = 0.25, z 0 = 4 mm, r 0 = 1.5 mm, m = 2 × 10 −17 kg, and Q = 950e. In our experimental setup, the AC voltage is applied to all four RF electrodes but with a phase shift of π between pairs of electrodes, often referred to as an RF-RF configuration, as depicted in Fig. 1. The DC offset V off is applied to one pair of electrodes.

Free oscillations
In the absence of external forces, that is, for F q,1 = F q,2 = 0, the collective motion consists of damped oscillations [q 1 , q 2 ,q 1 ,q 2 ] T = e exp (λ q t) with four complex eigenfrequencies λ q ∈ λ ± q,in , λ ± q,out , defined as the eigenvalues of matrix M in Eq. A2, and corresponding eigenvectors e ∈ a q,in , b q,in , λ ± q,in a q,in , λ ± q,in b q,in T , a q,out , b q,out , λ ± q,out a q,out , λ ± q,out b q,out T . The index "in" corresponds to the mode of collective oscillations with both oscillators moving in phase, while the index "out" corresponds to the out-of-phase mode. We can write the complex eigenfrequencies in terms of their real and imaginary parts: The real parts set the damping rate, and the imaginary parts determine the frequency of oscillation. Without loss of generality, we consider only eigenvalues with a positive imaginary part. Our system thus has six modes of collective oscillation, two along each spatial axis. The eigenvalues of matrix M are solutions of a quartic algebraic equation and have cumbersome expressions. We calculate them with a script in Wolfram Mathematica (version 12.1). For identical particles, meaning ω 0x,1 = ω 0x,2 = ω 0x , γ 01 = γ 02 = γ gas and µ = 1, and under the conditions γ fb,x + γ gas ω 0x and j/2ω 0x ω 0x , the frequencies can be approximated by [2] Instead of analyzing the general expression, we focus on the modal structure of the system of two nanoparticles used in the experiments described in the main text. For comparison with the data presented in Fig. 3 of the main text, we restrict ourselves to motion along the x axis. The system parameters are given in Table I. The evolution of the modal frequencies ω x,in /2π and ω x,out /2π is plotted in Fig. 2a for cooling rates up to γ fb,x = 300 Hz. We plot both the approximate solutions A13 and the exact eigenvalues of M in Eq. A2. As γ fb,x is increased, the frequency of the out-of-phase mode grows, while the frequency of the in-phase mode drops. When γ fb,x reaches the value 2j x /ω 0x = 135 Hz, the approximate solutions yield equal frequencies for the in-phase and out-of-phase modes. The exact solutions show that the frequency of the in-phase mode stabilizes at ω 0x − j x /2ω 0x = 2π × 2020 Hz, while the frequency of the out-of-phase mode approaches this value but then starts to decrease and reaches 0 Hz in the limit of an infinite feedback cooling rate. The evolution of the modal damping rates γ x,in and γ x,out is shown in Fig. 2b. For values of γ fb,x below 2j x /ω 0x , the modal damping rates are equal and increase monotonically. For values of γ fb,x above 2j x /ω 0x , the damping rate of the out-of-phase mode continues to increase while the damping rate of the in-phase mode approaches γ gas .
Next, we analyze the energies of the individual particles. Each particle contains some fraction of the energies of the in-phase and out-of-phase modes. These fractions can be calculated as [3] D in = a 2 x,in where D is the energy fraction in the directly cooled particle; I is the energy fraction in the sympathetically cooled particle. Figures 2c and 2d show the energy fractions in the out-of-phase and in-phase modes. As the cooling rate is increased, the energy of the directly cooled particle becomes concentrated in the out-of-phase mode, while the energy of the sympathetically cooled particle becomes concentrated in the in-phase mode.

Appendix B: Temperature estimation
This section describes how to extract the temperature of the CoM motion of the particles along axis q using two different methods [4]. The first method is based on fitting power spectral densities (PSDs) S qq,1 (ω) and S qq,2 (ω) with Lorentzian functions, while the second method relies on calculating the area under the curves Sqq ,1 (ω) = ω 2 S qq,1 (ω) and Sqq ,2 (ω) = ω 2 S qq,2 (ω), where the subscripts 1 and 2 indicate the particle number. We use the first method to calculate the temperatures of the individual modes and the second method to obtain the total kinetic energy of a single particle along a given spatial axis x, y, or z. In the latter case, the motion comprises two modes: the in-phase mode and the out-of-phase mode.

Temperature estimation of a mode of collective oscillation
In the weak cooling regime, defined as γ fb,q j q /ω 0q,1 and γ fb,q j q /ω 0q,2 , we can consider each mode as an independent oscillator. Each PSD of particle motion can be written as the sum of two Lorentzian functions: where S (th) q,in = 4 are the PSDs of the thermal force acting on each mode, and are the effective masses of the modes. We extract the modal temperature of each particle by fitting the PSD for each mode to a Lorentzian [4,5]. While the measured PSDs can be expressed in units of m 2 Hz −1 , it is not necessary for temperature estimation, so we use arbitrary units for power spectral density plots.

Temperature estimation of the center of mass of an individual particle
If the individual modes of collective oscillations cannot be resolved, for example, when the damping of the motion is larger than the normal mode splitting, the particle's temperature along axis q can be calculated based on the average kinetic energy [4] The variance q 2 can be calculated as where Sqq(ω) is the velocity PSD and S qq (ω) is the displacement PSD. However, since the experimentally determined S qq contains a detection-noise component, integration over all frequencies will lead to an overestimate of the velocity variance. Therefore, we limit the integration range to frequencies at which the detection signal is dominated by the particle motion.
After one or more nanoparticles have been loaded into the trap, for the purposes of charging, we set the laser's pulse energy below the threshold at which nanoparticles are desorbed from the foil. We observe that for particles 300 nm in diameter, the charging process saturates at ∼ 900 e. Thus, for two particles, if both are charged to the saturation level, we achieve an equal number of charges on each one. indicates the mean value of the power ratio, and the shaded region indicates the uncertainty. Inset: camera images of the particles at different illumination powers. The power ratio allows the ratio r1/r2 of the particle radii to be determined. (b) Amplitude response of the first particle as a function of the frequency of an external electric field (in addition to the Paul trap potential). (c) Amplitude response of the second particle as a function of the external field frequency.
In this section, we first determine the particles' relative masses, for which we measure the power P 1 and P 2 of the light scattered from each particle. The power is proportional to the sixth power of the particle radius; the ratio of the radii r 1 and r 2 is thus r 1 /r 2 = (P 1 /P 2 ) 1/6 . In Fig. 4, the ratio P 1 /P 2 is plotted for a range of illumination powers. The scattered power is determined from a camera image: the values of the pixels' brightness are summed in a region of interest around each particle. From these ratios, we find r 1 /r 2 = 1.00 ± 0.01. The inset in Fig. 4 shows examples of the images of two particles for different illumination powers.
Next, we determine the particles' relative charge, for which we measure the particles' responses to an external sinusoidal electric force, plotted in Figs. 4b and c. Only resonances corresponding to the particles' in-phase motion are observed. The particles' motion is coupled to the applied electric field through the particles' charge. The driving electric field acts on both particles in phase. For equal charge on the particles, the force will only couple to the in-phase mode, which is consistent with our observation. If the charge on the two particles were different, the outof-phase modes at 2.0 and 2.5 kHz would also be excited. A calculation of the response function shows that for the experimental parameters used to obtain the data in Figs. 4b and c, a 5% difference in charge number would result in an excitation of the out-of-phase mode with a peak height of 3% of the in-phase mode, or −15 dB for the data shown in Figs. 4b and c. We conclude that the charge on the two particles is equal within 5%.

Appendix F: Transfer function for sympathetic detection
This section provides details on the reconstruction of particle P1's response to an external sinusoidal force, applied only to P1, from a measurement of particle P2's motion. We focus on motion along the radial axes x and y; the case of motion along z is analogous. In the final experiment described in the main text, the force is applied to P1 in the directionx +ŷ, and the particles' motion is detected in the directiond =x −ŷ, as shown in Figs. 5a   An intensity-modulated 1064 nm laser beam exerts a force on the top particle, P1. Each particle's motion is measured with an independent detection system, as described in the main text. (c) Measured (solid blue line) and fitted (dashed orange line) response of the particle P2 to the drive of the particle P1. (d) Calculated amplitude and phase transfer functions that connect the motion of the particle P1 and P2.
From Eqs. A19, A21, A23 and A25 and the fit parameters, we can fully determine the transfer function Θ. The amplitude and phase components of this function are plotted in Fig. 5d. Note, however, that in order to compare the detected motion of the particles, as in Fig. 4 of the main text, we need to normalize the detection sensitivity for each particle. To that end, we compare the thermally driven oscillations of the particles at half of the frequency of the slowest mode. We then include this normalization factor in the transfer function Θ.