Theory for Cavity Cooling of Levitated Nanoparticles via Coherent Scattering: Master Equation Approach

TL;DR

This paper presents a comprehensive theoretical framework for cavity cooling of levitated nanoparticles via coherent scattering, deriving a master equation and demonstrating the potential for ground-state cooling in three dimensions.

Contribution

We develop an analytical master equation model for cavity cooling of levitated nanoparticles through coherent scattering, explicitly calculating decoherence rates and benchmarking with experiments.

Findings

Reproduces experimental results for 3D cooling

Demonstrates feasibility of ground-state cooling along all axes

Provides explicit decoherence rate calculations

Abstract

We develop a theory for cavity cooling of the center-of-mass motion of a levitated nanoparticle through coherent scattering into an optical cavity. We analytically determine the full coupled Hamiltonian for the nanoparticle, cavity, and free electromagnetic field. By tracing out the latter, we obtain a Master Equation for the cavity and the center of mass motion, where the decoherence rates ascribed to recoil heating, gas pressure, and trap displacement noise are calculated explicitly. Then, we benchmark our model by reproducing published experimental results for three-dimensional cooling. Finally, we use our model to demonstrate the possibility of ground-state cooling along each of the three motional axes. Our work illustrates the potential of cavity-assisted coherent scattering to reach the quantum regime of levitated nanomechanics.

Figures

Tables

Table 1. Table 1: Input parameters for our model. All the values except for the last two are taken from Ref. Windey  et al. . Although the tweezer waist is not directly measured, an estimation of its order of magnitude can be drawn (see main text). Note that, because of the chosen convention, the rate κ 𝜅 \kappa in our theory corresponds to half the value of κ 𝜅 \kappa reported in Ref. Windey  et al. .

Parameter	Value
Tweezer power	$P_t=0.5$ W
Tweezer wavelength	${\lambda }_0=1.55\mu$m
Tweezer waist	$W_t\sim 1\mu$m
Cavity length	$L_c=6.46\mathrm{mm}$
Cavity waist	$W_c=48\mu$m
Cavity linewidth	$2\kappa =2\pi \times 1.06\mathrm{MHz}$
Radius of NP	$R\sim \left(70\pm 20\right)\mathrm{nm}$
${\text{SiO}}_2$ permittivity at $1.55\mu$m  Malitson (1965)	$\epsilon =2.07$
Density of ${\text{SiO}}_2$   Haynes (2011)	$\rho =2200{\text{kg/m}}^3$

Table 2. Table 2: Heating rates contributing to the position localization rate Γ j subscript Γ 𝑗 \Gamma_{j} in Eq.  90 , for the parameters in Table 1 . The second, third, and fourth rows show, respectively, the heating rate due to pressure ( P mbar = subscript 𝑃 mbar absent P_{\rm mbar}= pressure in mbar), the photon recoil heating rate, and the trap displacement heating rate.

Axis	X	Y	Z
${\Gamma }_j^{\left(p\right)}/2\pi \text{ (GHz)}$	$28.6P_{\mathrm{mbar}}$	$24.5P_{\mathrm{mbar}}$	$85.7P_{\mathrm{mbar}}$
${\Gamma }_j^{\left(r\right)}/2\pi \text{ (kHz)}$	$0.09$	$0.15$	$1.89$
${\Gamma }_j^{\left(d\right)}/2\pi \text{ (kHz)}$	$94{\sigma }_x^2$	$110{\sigma }_y^2$	$31{\sigma }_z^2$

Table 3. Table 3: Summary of the most relevant parameters defined in the main text. The corresponding primed symbols denote the same parameters after the re-normalization induced by tracing out the free EM modes.

Symbol	Physical meaning	Definition and related parameters
${\Omega }_j$ $\left({\Omega }_j^{\prime }\right)$	Bare (renormalized) COM motional frequencies	Eq. 25 (Eqs. 169, 158)
$\kappa$ $\left({\kappa }^{\prime }\right)$	Bare (renormalized) COM motional frequencies	Sec. II, main text (Eqs. 172 and 150)
${\delta }^{\prime }$	Renormalized detuning between cavity and tweezer	Eq. 166 (see also Eqs. 52, 56, 132, 141, 147)
${\beta }_j$ $\left({\alpha }_c\right)$	Displacements of the COM operators (of the cavity operators)	Eq. 60 (Eq. 61) Appendix A
${\alpha }_c{g}_{cj}$	Bare OM coupling induced by cavity occupation	Eqs. 34 and 35
$G$	Bare OM coupling induced by photon scattering from the tweezer into the cavity	Eq. 46
$g_j$ $\left(g_j^{\prime }\right)$	Total bare (renormalized) optomechanical coupling	Eq. 65 (Eq. 171, see also Eqs. 167, 170, and 161)
$\gamma$	Gas pressure damping coefficient	Eq. 87
${\Gamma }_j^{\left(r\right)}$ / ${\Gamma }_j^{\left(p\right)}$ / ${\Gamma }_j^{\left(d\right)}$	Heating rates associated to photon recoil / gas pressure / trap displacement noise	Eq. 157 / Eq. 91 / Eq. 83
${\Gamma }_j$	Total heating rate	Eq. 90

Equations

$$\hat{H}=\frac{{\hat{\mathbf{P}}}^2}{2m}+{\hat{H}}_F+{\hat{H}}_{\mathrm{Int}}.$$

$${\hat{H}}_F=\frac{{\epsilon }_0}{2}\int {\text{d}}^3\mathbf{r}\left[{\hat{\mathbf{E}}}^2(\mathbf{r})+c^2{\hat{\mathbf{B}}}^2(\mathbf{r})\right],$$

$${\hat{H}}_{\mathrm{Int}}=-\frac{1}{2}\alpha {\hat{\mathbf{E}}}^2(\hat{\mathbf{R}}),$$

$$\hat{\mathbf{E}}(\mathbf{r})\approx {\mathcal{E}}_{\mathrm{tw}}(\mathbf{r},t)+{\hat{\mathbf{E}}}_{\mathrm{cav}}(\mathbf{r})+{\hat{\mathbf{E}}}_{\mathrm{free}}(\mathbf{r}),$$

$${\mathcal{E}}_{\mathrm{tw}}(\mathbf{r},t)=\frac{1}{2}\left[{\mathcal{E}}_{\mathrm{tw}}(\mathbf{r})e^{\text{i}{\omega }_{0}t}+{\mathcal{E}}_{\mathrm{tw}}^{\ast }(\mathbf{r})e^{-\text{i}{\omega }_{0}t}\right].$$

$${\hat{\mathbf{E}}}_{\mathrm{cav}}(\mathbf{r})=\sum _{\alpha }\sqrt{\frac{\hslash {\omega }_c}{2{\epsilon }_0{V}_c}}\left[{\mathbf{E}}_{c\alpha }(\mathbf{r}){\hat{c}}_{\alpha }+\text{H.c.}\right],$$

$${\hat{\mathbf{E}}}_{\mathrm{free}}(\mathbf{r})=\sum _{\mathbf{k}\epsilon }\sqrt{\frac{\hslash {\omega }_{\mathbf{k}}}{2{\epsilon }_0\mathcal{V}}}\left[{\mathbf{\epsilon }}_{\mathbf{k}}{e}^{\text{i}\mathbf{kr}}{\hat{a}}_{\epsilon \mathbf{k}}+\text{H.c.}\right].$$

$${\hat{H}}_F\approx \hslash {\omega }_c\sum _{\alpha }{\hat{c}}_{\alpha }^{†}{\hat{c}}_{\alpha }+\hslash \sum _{\mathbf{k}\epsilon }{\omega }_{\mathbf{k}}{\hat{a}}_{\epsilon \mathbf{k}}^{†}{\hat{a}}_{\epsilon \mathbf{k}}.$$

$$\begin{split}\hat{H}_{F}&=\hbar\omega_{c}\sum_{\alpha}\hat{c}_{\alpha}^{\dagger}\hat{c}_{\alpha}+\hbar\sum_{\mathbf{k}\varepsilon}\omega_{\mathbf{k}}\hat{a}_{\varepsilon\mathbf{k}}^{\dagger}\hat{a}_{\varepsilon\mathbf{k}}\\ &\!\!\!\!+\!\hbar\int_{0}^{\infty}\!\!\!\text{d}\omega\bigg{[}\omega\hat{a}_{0}^{\dagger}(\omega)\hat{a}_{0}(\omega)+\text{i}\gamma(\omega)\!\left[\hat{a}_{0}(\omega)\hat{c}^{\dagger}\!-\!\text{H.c.}\right]\!\!\bigg{]}.\end{split}$$

$${\mathcal{E}}_{\mathrm{tw}}(\mathbf{r})=E_0\frac{W_t}{W(z)}{e}^{-x^2/W_x^2(z)}{e}^{-y^2/W_y^2(z)}{e}^{\text{i}{k}_{0}z}{e}^{\text{i}{\varphi }_t(\mathbf{r})}.$$

$${\varphi }_t(\mathbf{r})=\arctan \left(\frac{z}{z_R}\right)-\frac{k_{0}z}{2}\frac{x^2+y^2}{z^2+z_R^2}\approx 0,$$

$${\mathcal{H}}_{\mathrm{tw}}(\mathbf{r},t)=\frac{1}{2}{\mathbf{e}}_y\left[{\mathcal{H}}_{\mathrm{tw}}(\mathbf{r})e^{\text{i}{\omega }_{0}t}+{\mathcal{H}}_{\mathrm{tw}}^{\ast }(\mathbf{r})e^{-\text{i}{\omega }_{0}t}\right],$$

$$P_t=\int \text{d}\mathbf{S}{⟨{\mathcal{E}}_{\mathrm{tw}}(\mathbf{r},t)\times {\mathcal{H}}_{\mathrm{tw}}(\mathbf{r},t)⟩}_T,$$

$$E_0=\sqrt{\frac{4P_t}{\pi {\epsilon }_{0}c{W}_t^2{A}_x{A}_y}}.$$

$$\begin{array}{rl}{\mathbf{E}}_{c\alpha }(\mathbf{r})& ={\mathbf{e}}_{\alpha }\frac{W_c}{W(y^{\prime })}{e}^{-(x^{\prime 2}+z^2)/W^2(y^{\prime })}\\ & \times \cos [k_c{y}^{\prime }+{\varphi }_c(\mathbf{r})-\varphi ].\end{array}$$

$$x^{\prime }=x\cos \Theta -y\sin \Theta .$$

$$y^{\prime }=x\sin \Theta +y\cos \Theta .$$

$$W(y^{\prime })=W_c\sqrt{1+{\left(y^{\prime }/y_R\right)}^2},$$

$${\varphi }_c(\mathbf{r})=\arctan \left(\frac{y^{\prime }}{y_R}\right)-\frac{k_c{y}^{\prime }}{2}\frac{x^{\prime 2}+z^2}{y^{\prime 2}+y_R^2},$$

$$\mathbf{e}_{\alpha}=\left\{\begin{array}[]{l}\mathbf{e}_{z}\\ \mathbf{e}_{\perp}=\mathbf{e}_{x}\cos\Theta-\mathbf{e}_{y}\sin\Theta.\end{array}\right.$$

$${\hat{H}}_{\mathrm{Int}}=-\frac{{\epsilon }_0{ϵ}_{c}V}{2}{\left[{\mathcal{E}}_{\mathrm{tw}}(\hat{\mathbf{R}},t)+{\hat{\mathbf{E}}}_{\mathrm{cav}}(\hat{\mathbf{R}})+{\hat{\mathbf{E}}}_{\mathrm{free}}(\hat{\mathbf{R}})\right]}^2.$$

$$\begin{array}{rl}{\hat{H}}_{\mathrm{t}-\mathrm{t}}& =-\frac{{\epsilon }_0{ϵ}_{c}V}{2}{\mathcal{E}}_{\mathrm{tw}}^2(\hat{\mathbf{R}},t).\end{array}$$

$$\begin{array}{rl}{\mathcal{E}}_{\mathrm{tw}}^2(\hat{\mathbf{R}},t)& \approx E_0^2{\cos }^2({\omega }_{0}t)\times \\ & \times \left[1-\left(\frac{2{\hat{X}}^2}{A_x^2{W}_t^2}+\frac{2{\hat{Y}}^2}{A_y^2{W}_t^2}+\frac{{\hat{Z}}^2}{z_R^2}\right)\right]\\ & +E_0^2(k_0\hat{Z})\left[k_0\hat{Z}\cos (2{\omega }_{0}t)-\sin (2{\omega }_{0}t)\right].\end{array}$$

$${\hat{H}}_{\mathrm{t}-\mathrm{t}}=\sum _{j=x,y,z}\frac{1}{2}m{\Omega }_j^2{\hat{R}}_j^2.$$

$$\left[\begin{array}[]{c}\Omega_{x}\\ \Omega_{y}\\ \Omega_{z}\end{array}\right]=\sqrt{\frac{\varepsilon_{0}\epsilon_{c}E_{0}^{2}}{2\rho W_{t}^{2}}}\left[\begin{array}[]{c}\sqrt{2}/A_{x}\\ \sqrt{2}/A_{y}\\ \lambda_{0}/(\pi W_{t})\end{array}\right],$$

$$\frac{{\hat{\mathbf{P}}}^2}{2m}+{\hat{H}}_{\mathrm{t}-\mathrm{t}}=\hslash \sum _{j=x,y,z}{\Omega }_j{\hat{b}}_j^{†}{\hat{b}}_j\equiv {\hat{H}}_{\mathrm{COM}}.$$

$${\hat{R}}_j=r_{j0}\left({\hat{b}}_j^{†}+{\hat{b}}_j\right);{\hat{P}}_j=\text{i}m{\Omega }_j{r}_{j0}\left({\hat{b}}_j^{†}-{\hat{b}}_j\right).$$

$$\begin{array}{rl}{\hat{H}}_{\mathrm{c}-\mathrm{c}}& =-\frac{{\epsilon }_0{ϵ}_{c}V}{2}{\hat{\mathbf{E}}}_{\mathrm{cav}}^2(\hat{\mathbf{R}}).\end{array}$$

$$\begin{array}{r}{\hat{H}}_{\mathrm{c}-\mathrm{c}}=-\hslash \frac{{\omega }_c{ϵ}_{c}V}{2V_c}|E_c(\hat{\mathbf{R}}){|}^2\left(\sum _{\alpha }{\hat{c}}_{\alpha }^{†}{\hat{c}}_{\alpha }+\frac{1}{2}\right).\end{array}$$

$$\begin{array}{rl}{\hat{H}}_{\mathrm{c}-\mathrm{c}}& =-\hslash \frac{{\omega }_c{ϵ}_{c}V}{2V_c}{\cos }^2(\varphi )\sum _{\alpha }{\hat{c}}_{\alpha }^{†}{\hat{c}}_{\alpha }\\ & -\hslash \frac{{\omega }_c{ϵ}_{c}V}{2V_c}{k}_c\sin (2\varphi ){\hat{Y}}^{\prime }\left(\sum _{\alpha }{\hat{c}}_{\alpha }^{†}{\hat{c}}_{\alpha }+\frac{1}{2}\right).\end{array}$$

Full text

I Introduction

Initially conceived as a way to minimize clamping losses in mechanical resonators Romero-Isart et al. (2010); Chang et al. (2010), the study of levitated nanoparticles (NPs), or levitodynamics, has branched into a wide research field in the last years. On the one hand, this is due to the wide range of particles available for levitation, such as dielectrics Li et al. (2011); Gieseler et al. (2012); Kiesel et al. (2013); Millen et al. (2015); Fonseca et al. (2016); Jain et al. (2016), nanocrystals containing quantum emitters Kuhlicke et al. (2014); Rahman and Barker (2017); Delord et al. (2018); Conangla et al. (2018); Alda et al. (2016); Frimmer et al. (2017), nanomagnets Rusconi et al. (2017) and superconducting spheres Romero-Isart et al. (2012); Cirio et al. (2012); Pino et al. (2018), as well as the large variety in NP shapes Asenbaum et al. (2013); Hoang et al. (2016); Ahn et al. (2018); Kuhn et al. (2017a). On the other hand, many recent experiments have demonstrated a very precise control of both the center of mass (COM) motion Jain et al. (2016); Diehl et al. (2018); Ahn et al. (2018); Tebbenjohanns et al. and the rotation Reimann et al. (2018); Kuhn et al. (2017b, 2015) of levitated NPs, as well as the integration of levitated NPs with optical emitters Rahman and Barker (2017); Conangla et al. (2018); Delord et al. (2018), and levitation in ultra high vacuum Jain et al. (2016). Such achievements pave the way toward interesting new possibilities, such as using NPs as inertial or force sensors Ranjit et al. (2016); Monteiro et al. (2017); Hebestreit et al. (2018), studying microscopic thermodynamics of the COM motion of a NP Gieseler et al. (2014); Millen et al. (2014); Gieseler and Millen (2018), or the potential to use levitated NPs for optomechanics Romero-Isart et al. (2010); Chang et al. (2010); Romero-Isart et al. (2011a) and for the preparation of large quantum superposition states Romero-Isart et al. (2011b); Romero-Isart (2011); Bateman et al. (2014).

The promising prospects of levitodynamic-based quantum applications rely on the ability of cooling the COM motion down to the ground state. Although this milestone has not yet been attained, significant advances toward this goal have been achieved, such as parametric feedback cooling Li et al. (2011); Gieseler et al. (2012); Diehl et al. (2018) or cavity-assisted optomechanical cooling Asenbaum et al. (2013); Kiesel et al. (2013); Millen et al. (2015). Among these, a particularly interesting option, initially proposed for cold atoms and molecules Vuletić et al. (2001); Domokos et al. (2002); Ritsch et al. (2013), is to cool the COM motion via coherent light scattering into an optical cavity, as recently demonstrated in two experiments Windey et al. ; Delić et al. . In this configuration, a NP is optically trapped by an optical tweezer Ashkin (1970); Ashkin and Dziedzic (1976), whose photons are scattered into a blue-detuned cavity, reducing the mechanical energy of the COM in the process. This cooling method presents several advantages such as reduced technological complexity, a high level of control and the possibility of cooling the COM motion along the three motional axes. It therefore appears as a strong candidate for realizing the COM ground state in levitated nanomechanical systems.

Previous theoretical works have addressed the method of cavity-assisted COM cooling by coherent scattering: either semiclassically for cold atoms Hechenblaikner et al. (1998); Domokos et al. (2001); Vuletić et al. (2001), or using quantum optomechanical theory to describe the cooling of the COM of atoms along one or two of the trapping axes Lev et al. (2008); Vukics et al. (2004), and of optically levitated dielectric NPs along one trapping axis Romero-Isart et al. (2011a). However, a full quantum theory of three-dimensional cooling of optically levitated dielectric NPs via coherent scattering, as well as a detailed study of the relevant heating and decoherence mechanisms, is lacking so far. On the one hand, developing such a theory will contribute to understand the limitations of present experiments Windey et al. ; Delić et al. and the way to overcome them in order to achieve ground state cooling. On the other hand, once the ground state is reached, a full quantum theory including a detailed quantum description of the decoherence sources will be essential not only to describe this state, but to implement further quantum protocols and applications to prepare, for instance, non-Gaussian quantum states.

In this paper, we develop a quantum theory of three-dimensional COM cooling via cavity-assisted coherent scattering. First, in Sec. II we derive the full optomechanical Hamiltonian in the long wavelength approximation. In Sec. III, we obtain an effective equation of motion for the reduced subsystem formed by cavity and COM motion, and introduce three relevant heating rates, namely the recoil heating, the heating related to the background gas pressure, and the trap displacement noise. In Sec. IV, we focus on a case study and characterize the system dynamics for a realistic experimental setup Windey et al. . We continue in Sec. V by analyzing the possibility of ground state cooling along each motional axes in state-of-the-art experiments. Finally, our conclusions are presented in Sec. VI.

II Hamiltonian of the total system

The system under study, schematically depicted in Fig. 1, consists of a single dielectric NP with radius $R$, and homogeneous and isotropic relative permittivity $\varepsilon$. The NP is trapped at the focus of a high-intensity optical tweezer propagating along the $z$ axis whose frequency is $\omega_{0}=2\pi c/\lambda_{0}$, where $c$ is the vacuum speed of light and $\lambda_{0}$ the tweezer wavelength. Additionally, the NP is placed inside an optical cavity and is coupled to a mode with frequency $\omega_{c}$. The cavity axis is orthogonal to the propagation direction of the tweezer. We include the two degenerate polarization states of the non-birefringent cavity with frequency $\omega_{c}$, although, as we will see below, only one will significantly contribute to COM cooling. In this section we describe the system and the fundamental light-matter Hamiltonian governing its dynamics. We then derive, for a NP confined close to the tweezer focus, a second quantization form for such Hamiltonian. Finally, we reduce this Hamiltonian to a quadratic form by transforming into a frame where non-quadratic terms become negligible.

The system detailed above, and depicted in Fig. 1, can be described by a Hamiltonian containing the kinetic energy of the free NP, the energy of the electromagnetic (EM) field, and the NP-field interaction, i.e.,

[TABLE]

Here, $\hat{\mathbf{P}}$ denotes the COM momentum operator, and $m$ denotes the mass of the NP. The free Hamiltonian of the EM field is given by

[TABLE]

with the transverse electric and magnetic field operators $\hat{\mathbf{E}}(\mathbf{r})$ and $\hat{\mathbf{B}}(\mathbf{r})$ at position $\mathbf{r}$. The NP-EM field interaction term in Eq. 1 reads Romero-Isart et al. (2011a)

[TABLE]

where we have introduced the NP polarizability $\alpha=\varepsilon_{0}\epsilon_{c}V$, with the NP volume $V$, the vacuum permittivity $\varepsilon_{0}$, and $\epsilon_{c}\equiv 3(\varepsilon-1)/(\varepsilon+2)$. The above expression is valid in the long wavelength approximation, i.e., when $R\ll\lambda$ for all relevant wavelengths $\lambda$ Novotny and Hecht (2006). The motion-light interaction arises from evaluating the electric field operator at the COM position of the NP, $\hat{\mathbf{R}}$. Note that matter is treated classically in this work, as the optical response of the NP is described exclusively by its polarizability.

In order to write the Hamiltonian in second quantization, we need to determine the structure of the EM field. Both in Eqs. 2 and 3, the electric and magnetic field operators should be written in terms of the eigenmodes of the corresponding EM structure. Since solving Maxwell’s equations in all space in the presence of a cavity is rather involved, it is common in photonics to approximate the EM field as independent free-space modes plus some extra modes representing the structure, i.e. the cavity. Here we use the same approximation, thus writing the electric field as

[TABLE]

whose three components correspond to the tweezer field, the cavity field, and the free EM field, respectively. Since we assume the tweezer field to be in a highly populated coherent state, we approximate it as a classical quantity by its mean value in the rotating frame,

[TABLE]

On the other hand, the electric fields of the cavity and the free EM modes are described quantum mechanically, i.e.,

[TABLE]

[TABLE]

Here, $V_{c}$, $\mathbf{E}_{c\alpha}(\mathbf{r})$ and $\hat{c}_{\alpha}$ represent the cavity mode volume, the normalized classical electric field of the cavity mode, and the annihilation operator for a cavity mode with polarization $\alpha$, respectively. Analogously, $\mathcal{V}$ and $\hat{a}_{\varepsilon\mathbf{k}}$ represent the quantization volume for the free EM field and the annihilation operator of a free EM mode with wavevector $\mathbf{k}$ and polarization $\boldsymbol{\varepsilon}_{\mathbf{k}}$. Note that, in our approximation, the Hamiltonian of the EM field in Eq. 2 can be written as

[TABLE]

However, in order to account for the cavity losses due to transmission through the mirrors, which play a relevant role in the cooling, it is convenient to refine the above Hamiltonian by explicitly including an additional set of EM modes to which the cavity is coupled Romero-Isart et al. (2011a); Gardiner and Zoller (2004), i.e., we take

[TABLE]

Here, the continuum of modes described by the photonic operators $\hat{a}_{0}(\omega)$ would in principle represent a fraction of the free-space modes. Thus, by considering them a separate degree of freedom, we are double-counting some of the states. This overcounting, very common in quantum optics Dutra (2005), has nevertheless a negligible effect on the correctness of the solution, since the ensemble of extra modes $\hat{a}_{0}(\omega)$ has zero measure Romero-Isart et al. (2011a). The coupling constant $\gamma(\omega)$ is directly related to the cavity linewidth $\kappa$, here defined as the cavity field decay rate, and can be considered constant across a wide frequency range centered at $\omega_{c}$ with a value $\gamma(\omega)=\sqrt{\kappa/\pi}$ Gardiner and Zoller (2004). This equality can be readily certified by obtaining the Born-Markov master equation of the cavity in the presence of the environment composed by the modes $\hat{a}_{0}(\omega)$ (see e.g. Appendix B).

In order to calculate the second-quantized interaction Hamiltonian $\hat{H}$ using the electric fields described above, it is convenient to determine the explicit form of the tweezer and the cavity field profiles. Let us start by the tweezer field, namely $\boldsymbol{\mathcal{E}}_{\rm tw}(\mathbf{r})$ in Eq. 5. Our first step is to fix the $x$ and $y$ axis in our coordinate system. Those axes are aligned with the symmetry axes of the tweezer-formed trapping potential. As pointed out in studies with atomic clouds Vuletić et al. (2001), the possibility of 3D cooling requires that none of these axes is parallel to the cavity axis. In the present case, this is achieved for tweezer polarizations being neither purely parallel nor purely orthogonal to the cavity axis, as schematically depicted in Fig. 2. Since the tweezer polarization axis determines the trapping directions of the levitated particle, it is convenient to choose a coordinate system centered on the tweezer intensity maximum, and aligned with the tweezer polarization. We therefore choose the tweezer to be polarized along the $x$ axis (see Fig. 2), i.e., $\boldsymbol{\mathcal{E}}_{\rm tw}(\mathbf{r})=\mathbf{e}_{x}\mathcal{E}_{\rm tw}(\mathbf{r})$, and modelled in the paraxial approximation by a classical, zeroth order, propagating Hermite-Gaussian beam Erikson and Singh (1994),

[TABLE]

In the above equation, $W(z)=W_{t}\sqrt{1+(z/z_{R})^{2}}$, where $W_{t}$ is the tweezer waist at the focus and $z_{R}=\pi W_{t}^{2}/\lambda_{0}$ is the Rayleigh range. The asymmetry of the transverse beam profile, which becomes more relevant for strongly focused tweezers, is encoded in the different extensions of the beam along the transverse axis, $W_{x,y}(z)=A_{x,y}W(z)$ with $A_{x,y}$ adimensional. The phase factor $\phi_{t}(\mathbf{r})$ is given by

[TABLE]

in the vicinity of the origin, i.e., for $x,y,z\ll z_{R}\approx 1\,\mu\rm m$ in typical experiments Windey et al. . Finally, the field amplitude $E_{0}$ can be related to the tweezer power $P_{t}$ by first obtaining the paraxial expression for the associated magnetic field Erikson and Singh (1994),

[TABLE]

where $\mathcal{H}_{\rm tw}(\mathbf{r})=\sqrt{\varepsilon_{0}/\mu_{0}}\mathcal{E}_{\rm tw}(\mathbf{r})$. One then defines the tweezer power in terms of the Poynting vector as

[TABLE]

where $d\mathbf{S}$ is the surface element perpendicular to the Poynting vector (in this case, $d\mathbf{S}=dxdy\mathbf{e}_{z}$), and $\langle\rangle_{T}$ denotes the time average. Using this expression we find

[TABLE]

In this way we can relate all the beam parameters to experimental input parameters Windey et al. , such as power, wavelength and beam waist (or, equivalently, numerical aperture of the focusing lens).

Let us now shift our attention to the cavity field $\mathbf{E}_{c\alpha}(\mathbf{r})$ in Eq. 6. We will also approximate its profile by a zeroth-order Hermite-Gauss beam, this time in a standing wave configuration. Because of our choice of coordinate system (Fig. 2), the axis of this beam will not be parallel to $x$ or $y$, but instead lie along along $\mathbf{e}_{\rm cav}=\sin\Theta\mathbf{e}_{x}+\cos\Theta\mathbf{e}_{y}$. The expression for the two degenerate polarization modes is thus given by

[TABLE]

Here, $k_{c}=\omega_{c}/c$, and $W_{c}$ is the cavity waist which, together with the cavity length $L_{c}=\pi c/(2\kappa\mathcal{F})$ ($\mathcal{F}$ being the Finesse), determines the mode volume as $V_{c}=\pi W_{c}^{2}L_{c}/4$. The above beam contains the rotated coordinates $x^{\prime}$ and $y^{\prime}$, which are given by

[TABLE]

[TABLE]

Importantly, the phase factor $\phi$ in Eq. 15 determines the field intensity at the position of the NP, i.e. at the origin of coordinates. For $\phi=0(\pi/2)$ the particle is at a maximum(minimum) of the cavity intensity profile. Finally, the functions $W(y^{\prime})$ and $\phi_{c}(\mathbf{r})$ in Eq. 15 are given by

[TABLE]

[TABLE]

with $y_{R}=\pi W_{c}^{2}/\lambda_{c}$ being the Rayleigh range of the cavity beam. Since we will be interested in the cavity field close to the origin, i.e., for $x,y\ll y_{R}\sim\,$mm in typical experiments Windey et al. , we will assume $\phi_{c}(\mathbf{r})\approx 0$. Finally, the polarization vectors in Eq. 6 are the two unit vectors orthogonal to the cavity axis, and will be labeled as

[TABLE]

Note that, by taking $\Theta\to 0$, we recover the usual Hermite-Gauss standing wave along the $y$ axis.

We are now in a position to derive the second quantized form of the total Hamiltonian, assuming that the NP always remains close to the origin of coordinates, i.e., close to the tweezer intensity maximum. We will follow the treatment done in Ref. Romero-Isart et al. (2011a), although in this paper our main focus will be a different optomechanical coupling, namely that arising from the coherent scattering from the tweezer into the cavity. We start by introducing the expression for the total electric field operator, Eq. 4, into the interaction Hamiltonian, Eq. 3, to obtain

[TABLE]

The expression above results in six different terms, which give rise to the different physical interactions in this system. Let us analyze them separately.

First, the tweezer-tweezer term $\propto\boldsymbol{\mathcal{E}}^{2}_{\rm tw}(\hat{\mathbf{R}},t)$ gives rise to the trapping of the NP. The interaction energy arising from this term is

[TABLE]

Since the NP will be confined close to the origin of coordinates, we approximate the electric field of the tweezer by its expansion close to $\hat{\mathbf{R}}=0$. Using Eqs. 5 and 10, this expansion is carried out to yield

[TABLE]

As the tweezer frequency $\omega_{0}$ is very large compared to the typical coupling rates between the different degrees of freedom, we can perform a Rotating Wave Approximation by neglecting all the rapidly oscillating terms $\exp(\pm 2i\omega_{0}t)$. Within this approximation, we introduce the above electric field into Eq. 22 and, after neglecting constant energy contributions, we obtain the usual harmonic potential,

[TABLE]

This potential represents the trapping of the NP by the tweezer, with trapping frequencies

[TABLE]

where $\rho=m/V$ is the mass density of the NP. Note that the presence of this trapping potential allows for the definition of the COM phonons. Specifically, we can combine the above interaction term with the kinetic energy of the COM, namely the first term in Eq. 1, to obtain

[TABLE]

Here, the bosonic operators $\hat{b}_{j}$ are defined in terms of the position and momentum operators through

[TABLE]

where we have defined the zero-point motion $r_{j0}=(\hbar/2m\Omega_{j})^{1/2}$.

The second contribution we consider is given by the square of the cavity field,

[TABLE]

Introducing Eq. 6 together with the field $\mathbf{E}_{c\alpha}(\mathbf{r})=\mathbf{e}_{\alpha}E_{c}(\mathbf{r})$ from Eq. 15, and keeping only the non-rotating terms 111This Rotating Wave Approximation (RWA) is analogous to the one taken for the tweezer field, and equally valid since $\omega_{c}\approx\omega_{0}$. Note that we take it at this point for simplicity, but one could keep the counter-rotating terms until Eq. 52, where they would vanish under the general RWA undertaken immediately afterwards., we find

[TABLE]

We now expand the square modulus in powers of the COM position, $\hat{\mathbf{R}}$, up to first order, obtaining

[TABLE]

where we have defined for convenience the rotated $y$ operator

[TABLE]

The interaction Hamiltonian $\hat{H}_{\rm c-c}$ contains three physically different effects: first, a shift of the cavity frequency that depends on the NP position inside the cavity,

[TABLE]

Second, an optomechanical coupling between COM and cavity modes,

[TABLE]

with

[TABLE]

[TABLE]

Lastly, a force on the COM along the $X-Y$ plane,

[TABLE]

For common experimental values Windey et al. , we find that $g_{\rm cj}\approx 2\pi\times 1$Hz. This allows us to neglect this third term, as it corresponds to a negligible shift in the equilibrium position of the COM motion ($\sim g_{\rm cj}/\Omega_{j}\sim 10^{-5}$ times the zero point motion). Note that the optomechanical coupling in Eq. 33 could still be relevant if the cavity occupation is large, so this term must be retained. We thus write the contribution of the term $\hat{\mathbf{E}}^{2}_{\rm cav}(\hat{\mathbf{R}})$ as

[TABLE]

The third contribution to Eq. 21 stems from the square of the free EM field, and is given by

[TABLE]

where we have used Eq. 7. This term has been proven to become negligible for particles smaller than the relevant wavelengths Pflanzer et al. (2012), and thus will be ignored hereafter.

Finally, we address the contributions arising from the three cross terms, namely

[TABLE]

[TABLE]

and

[TABLE]

which result in interactions between the three system components. They are all constructed in the same way, namely by expanding the corresponding electric fields close to the origin,

[TABLE]

[TABLE]

[TABLE]

and keeping terms of up to first order in the COM position. Among the three cross terms, the most important is the tweezer-cavity contribution Eq. 39, which will responsible for the cooling via coherent scattering. It is given by

[TABLE]

where

[TABLE]

Note that $\hat{H}_{\rm t-c}$ contains both a displacement of the cavity modes induced by the trapping field, and an interaction between cavity and COM degrees of freedom, the latter of which will ultimately be responsible for the cooling. Note that it follows from Eq. 20 that the coupling $G_{z}$ vanishes, i.e. the $z-$polarized cavity mode plays a negligible role in the dynamics since it will not be populated by the $x-$polarized tweezer. Because of this argument, we will neglect this mode and reduce the cavity to a single mode polarized along $\mathbf{e}_{\perp}$.

The second cross term, namely the one coming from the tweezer-free field product, reads

[TABLE]

with the coupling rate

[TABLE]

The term $\hat{H}_{\rm t-f}$ contains a displacement of the free space EM modes, and a tweezer-mediated interaction between free space modes and COM motion. As we will see below, the latter interaction will result in the recoil heating of the NP Romero-Isart et al. (2011a); Jain et al. (2016).

Finally, the product between cavity and free EM field adds a contribution

[TABLE]

with a coupling factor

[TABLE]

This interaction term is more involved than the previous two. On the one hand, it contains a quadratic interaction between cavity modes and free field, responsible for extra losses of the cavity into free space. Since these losses are mediated by the NP, they depend on its position through $\cos\phi$. Note that this term is small for sub-wavelength particles, as otherwise the cavity linewidth would be modified by the presence of the NP, an effect not observed in experiments Windey et al. . On the other hand, $\hat{H}_{\rm c-f}$ contains a non-quadratic, three body interaction between free space photons, cavity photons, and COM phonons. This contribution, however, can be safely neglected since it is of order $\sim kr_{j0},k_{c}r_{j0}\sim 10^{-6},10^{-5}$ times smaller than the quadratic contribution at optical frequencies. This argument might fail if the cavity is largely populated by the tweezer, where this three-body interaction results in an effective free field-COM coupling which contributes to the recoil heating. However, this extra, cavity-induced recoil will be negligible compared to the recoil induced by the much highly occupied tweezer (Eq. 47), which is known to dominate Romero-Isart et al. (2011a); Jain et al. (2016). Therefore, we can safely approximate

[TABLE]

The full Hamiltonian of the system is obtained by adding up all the above contributions, namely Eqs. 9, 26, 37, 45, 47, and 49. As discussed above, from now on we will ignore the $z-$polarized cavity mode and refer to the remaining creation operators simply as $\hat{c}_{\perp}\equiv\hat{c}$. In terms of the COM creation and annihilation operators, the full Hamiltonian reads

[TABLE]

where $G\equiv G_{\perp}$, $G(\mathbf{k})\equiv G_{\perp}(\mathbf{k})$, and $\tilde{\omega}_{c}=\omega_{c}-\Delta_{c}$.

The above Hamiltonian, though simpler than the original, remains very challenging to solve, but it can be further simplified. First, we elliminate the time dependence by virtue of the same Rotating Wave Approximation undertaken after Eq. 23. In order to do this, we first perform a unitary transformation into a frame rotating with the tweezer frequency, i.e.,

[TABLE]

with

[TABLE]

After applying this transformation, the Hamiltonian will contain two different kinds of terms, namely the non-rotating terms which do not depend on time, and contributions rotating at frequencies $\pm 2\omega_{0}$. We then take the Rotating Wave Approximation by neglecting the latter. After such approximation, the Hamiltonian is reduced to

[TABLE]

where we have defined the detunings

[TABLE]

[TABLE]

[TABLE]

Note that the Hamiltonian contains a displacement of the cavity and the free EM modes (second line of Eq. 55). As usual in quantum optics, it is convenient to remove such displacements by means of a second unitary transformation, which displaces all the system modes at once:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

After substitution of the above operators into Eq. 55, we neglect the constant energy shifts and set the terms proportional to $\hat{a}_{\varepsilon\mathbf{k}}$, $\hat{b}_{j}$, $\hat{c}_{\alpha}$, and $\hat{a}_{0}(\omega)$ to zero, so that any displacement will vanish from the transformed Hamiltonian. In this way we obtain a system of equations relating the coefficients $\alpha_{\mathbf{k}}$, $\beta_{j}$, $\gamma_{\alpha}$, and $\alpha_{0}(\omega)$, whose solution is given in Appendix A together with the general form of the transformed Hamiltonian.

After the above displacement, the operators appearing in the Hamiltonian represent the fluctuations of the corresponding degree of freedom above a classical value. Since these fluctuations are small, we can neglect the remaining non-quadratic terms in the transformed Hamiltonian. Then, in the transformed frame, the Hamiltonian can be written in a very compact form as

[TABLE]

Here, we are describing the system in the usual notation in open quantum systems Breuer et al. (2002), where we divide our degrees of freedom into a system S, composed by the cavity and the COM modes, and two reservoirs, RA and RB, composed by the free EM modes and the output modes of the cavity. Such distinction facilitates the procedure outlined in the next section, namely the tracing out of the reservoir modes to get a reduced equation of motion only for the cavity and COM subsystem. In this system plus reservoir picture, the Hamiltonian of the system S is defined as

[TABLE]

where $\hat{V}_{0}=\sum_{j}g_{j}\hat{c}^{\dagger}(\hat{b}^{\dagger}_{j}+\hat{b}_{j})+\text{H.c.}$ is the interaction between the system degrees of freedom, with coupling rates

[TABLE]

On the other hand, the two reservoirs are governed by the Hamiltonians

[TABLE]

[TABLE]

Finally, the interaction between system and reservoirs is given by three independent terms, namely

[TABLE]

[TABLE]

and

[TABLE]

where the coupling rates are given, respectively, by

[TABLE]

[TABLE]

By using the above definitions, the quadratic Hamiltonian describing both the system and the reservoirs, Eq. 63, takes a suitable form for tracing the latter out.

III Reduced dynamics of the cavity and the center-of-mass modes

The Hamiltonian in Eq. 63 describes an infinite system of coupled harmonic oscillators. Since we are only interested in the reduced dynamics of the system formed by COM and cavity, in this section we trace out the two sets of continuum EM modes to obtain an effective equation of evolution for such subsystem. The procedure for tracing out the reservoir degrees of freedom starts by transforming into the Interaction Picture with respect to the free evolution of all degrees of freedom, i.e. with respect to $\hat{H}_{RA}+\hat{H}_{RB}+\hat{H}_{S}-\hat{V}_{0}$. In this picture, the evolution of the total density matrix $\rho$ is given by the von Neumann Equation Breuer et al. (2002),

[TABLE]

where $\hat{V}_{0}(t)$ and $\hat{V}(t)=\hat{V}_{A}(t)+\hat{V}_{B1}(t)+\hat{V}_{B2}(t)$ represent the two interaction potentials in the Interaction Picture. In order to solve the above equation, we undertake the weak coupling or Born approximation, which consists on approximating the full density matrix as a product state, i.e., $\hat{\rho}(t)=\hat{\mu}(t)\otimes\hat{\rho}_{R}$. Here, $\hat{\mu}(t)=\text{Tr}_{R}\hat{\rho}(t)$ is the reduced density matrix of the system, $\text{Tr}_{R}$ denoting the partial trace over the reservoir modes. On the other hand, $\hat{\rho}_{R}$ is the reduced density matrix of the reservoir, which is assumed constant on the basis of the environment being composed by an infinitely large amount of modes, whose state is therefore only negligibly modified by the presence of the system. We will assume that the density matrix describing the free EM reservoirs is a thermal state at room temperature,

[TABLE]

Note that, in principle, the assumption of a thermal state for the reservoirs is only legitimate in the original frame, and such thermal state should be transformed accordingly into the displaced frame given by the transformation $\{\hat{a}_{\mathbf{\varepsilon k}},\hat{a}_{0}(\omega)\}\to\{\hat{a}_{\mathbf{\varepsilon k}}+\alpha_{\mathbf{k}},\hat{a}_{0}(\omega)+\alpha_{0}(\omega)\}$. However, this state remains a good approximation in the transformed frame since, as evidenced by our results in Appendix A, the contribution of the reservoir modes is only relevant at frequencies close to $\omega_{0}$. At these frequencies, the presence of a highly occupied classical tweezer field makes any other dynamics of the EM field negligible. Because of this reason, the thermal state in Eq. 74 remains a good approximation for reproducing experimental results, as we will see below.

Under the Born approximation and with the reservoirs in a thermal state, we can formally solve Eq. 73, reinsert it into itself, and take the trace over the reservoir modes, obtaining

[TABLE]

The second argument inside the commutator in the first line above can be identified as $\hat{\mu}(t)$, as can be readily checked by taking the partial trace and formally solving for $\hat{\mu}(t)$ in Eq. 73. Moreover, in the second line of Eq. 75, it is customary to undertake the Markov approximation, which assumes the reservoir correlation functions decay much faster than the system-reservoir interaction rate. Formally, this is equivalent to approximating $\hat{\mu}(s)\approx\hat{\mu}(t)$ and extending the upper integration limit to infinity Breuer et al. (2002). The final Master equation reads

[TABLE]

The Master Equation 76 now has the desired form, namely an equation of motion involving only the system degrees of freedom, i.e. the cavity and the COM modes. This equation contains the free evolution of the system plus some extra terms, describing the effect of the reservoirs (second line of Eq. 76). The explicit calculation of such terms, carried out in Appendix B, leads to the following Master Equation in the Schrödinger Picture,

[TABLE]

Here, the Hamiltonian $\hat{H}_{S}^{\prime}$ has the same form as the Hamiltonian $H_{S}$, but with all the involved frequencies re-normalized (shifted) by the effect of the reservoirs,

[TABLE]

where we define $\hat{q}_{j}\equiv\hat{b}_{j}^{\dagger}+\hat{b}_{j}$, and the expressions for $\delta^{\prime}$, $\Omega_{j}^{\prime}$, and $g_{j}^{\prime}$ are given in Appendix B. Note that the shifts in the frequencies and coupling rates represent the conservative part of the reservoir-induced system dynamics. On the other hand, the term $\mathcal{D}[\hat{\mu}]$ contains the incoherent (or non-conservative) dissipators, which include cavity losses at a rate $\kappa^{\prime}$, recoil heating of the COM modes at a rate $\Gamma_{j}^{(r)}$, and incoherent cavity-COM interaction at a rate $\Upsilon$:

[TABLE]

where the curly brackets denote the anti-commutator. Here, the cavity linewidth $\kappa^{\prime}$ contains in principle also a small correction due to the presence of the NP. The explicit expressions for $\kappa^{\prime}$, $\Gamma_{j}^{(r)}$, and $\Upsilon$ are given in Appendix B.

III.1 Other noise sources

In a typical levitodynamics experiment, the noise sources associated to thermal free photons are not the only ones. Indeed, some noise sources, stemming from different degrees of freedom not accounted for in our original Hamiltonian (Eq. 1), can play a relevant role in the system dynamics. Therefore, we must include such sources in our equation of motion for the cavity and COM degrees of freedom. In this work we focus on two particular decoherence channels, namely displacement noise in the trap and residual gas pressure.

We first focus our attention on the displacement noise associated with the “shaking” of the center of the trap. This mechanism has been discussed in detail in the literature Pino et al. (2018); Joos et al. (2003); Breuer et al. (2002); Henkel et al. (1999) and results in an extra dissipator in the Master Equation of the position localization or Brownian motion form,

[TABLE]

where $\hat{R}_{j}=r_{j0}\hat{q}_{j}$ is the $j-$component of the COM position operator. The magnitude $\Lambda_{j}r_{j0}^{2}$ is the corresponding dissipation rate, and can be related to observables in the following way: let us describe the trap displacements along each direction $j$ by means of independent classical fluctuating variables $\xi_{j}(t)$, such that $\hat{R}_{j}\to\hat{R}_{j}+\xi_{j}(t)$. We assume such variables to have zero mean and nonzero fluctuations, i.e. $\langle\xi_{j}(t)\rangle_{T}=0$ and $\langle\xi_{j}(t)\xi_{j}(t^{\prime})\rangle_{T}\neq 0$ Pino et al. (2018). In terms of these variables, we can define the following two-sided power spectral density (PSD)

[TABLE]

which has units of m${}^{2}/$Hz. For the particular case $S_{jj}^{(d)}(\Omega_{j})=S_{jj}^{(d)}(-\Omega_{j})$ one can demonstrate, by averaging over the stochastic force generated by $\xi_{j}$, that the displacements result in the dissipator Eq. 80, with

[TABLE]

The dissipation rates associated to the trap displacement, $\Gamma_{j}^{(d)}\equiv\Lambda_{j}r_{0j}^{2}$, can therefore be written as

[TABLE]

i.e., the ratio $\Gamma_{j}^{(d)}/\Omega_{j}$ is proportional to the PSD in units of $r_{0j}^{2}/\Omega_{j}$, whose square root we define as $\sigma_{j}$. Note that the above rate scales linearly with the NP mass, therefore taking much higher values for a NP than for a trapped atom. Indeed, as we will see below, in recent levitodynamic cooling experiments Windey et al. , this mechanism could represent a relevant source of heating.

Let us now focus on the effect of gas pressure which, by inducing an extra heating of the COM motion, limits the cooling power of the experiment. The pressure of the gas is modelled through a combination of two dissipators Joos et al. (2003); Rodenburg et al. (2016),

[TABLE]

Here, the first term also takes the form of a position localization dissipator,

[TABLE]

with $T$ the temperature of the gas, whereas the second term describes viscous friction, and is given by

[TABLE]

where we define the conjugate variable $\hat{p}_{j}\equiv\hat{b}^{\dagger}_{j}-\hat{b}_{j}$. The rate $\gamma$ can be obtained from the kinetic theory of gases Beresnev et al. (1990); Li et al. (2011), and reads

[TABLE]

where $\eta_{g}=\bar{l}P(2m_{0}/\pi K_{B}T)^{1/2}$ is the viscosity of the gas, $\bar{l}$ the mean free path of the air molecules, $P$ the pressure and $m_{0}$ the molecular mass of the gas.

Both the first dissipator associated with the gas pressure, $\mathcal{D}_{R}[\mu]$, and the displacements of the trap center, $\mathcal{D}_{\text{d}}[\mu]$, have the same form as the recoil heating term in Eq. 79, and can thus be grouped together. The final Master Equation then reads

[TABLE]

where $\hat{H}_{S}^{\prime}$ is given by Eq. 78, and the final dissipator reads

[TABLE]

with

[TABLE]

and

[TABLE]

Equation 88 is the final equation of evolution for the system formed by cavity and COM motion. As a reading guide, a compilation of the most relevant parameters governing the system evolution is shown in Table 3 (Appendix D). Note that, although we are neglecting any parametric noise induced by fluctuations in the trapping frequencies, such noise could be incorporated in our model by means of a Brownian motion dissipator for the operators $\hat{R}_{j}^{2}$ Schneider and Milburn (1999). This, however, lies beyond the scope of the present work since, as will be shown below, the three heating sources introduced here, namely gas pressure, recoil heating, and displacement noise, are both necessary and sufficient for our model to be compatible with experimental observations.

IV Results: 3D cavity cooling in recent experimental setup

In this section we study the behavior of the system for realistic experimental parameters. In order to illustrate our model, we take as a case study the recent experiment with a SiO2 NP by Windey et al. Windey et al. . The values of the parameters measured in this experiment are shown in Table 1, together with the permittivity and density of Silica extracted from the literature. Note that not all the parameters in our model are measured directly, leaving some of them free to fit the experimental results. For instance, the radius of the NP is only known within a relatively wide range, and in the following we will take $R=50\,\rm nm$ for the sake of definiteness.

Our first task is to evaluate all the parameters appearing in our effective equation of motion, Eq. 88, starting by the COM trapping frequencies $\Omega_{j}^{\prime}$. To obtain them, we need to fix the values of the tweezer waist, $W_{t}$, and the asymmetry factors $A_{x}$ and $A_{y}$. Reasonable values for the former lie on the order of $W_{t}\sim\lambda_{0}/(\pi\text{NA})\sim 1\,\mu$m, where $\text{NA}\approx 0.8$ is the numerical aperture of the lens Novotny and Hecht (2006). On the other hand, since the tweezer cross section is not expected to deviate significantly from a cylindrically symmetric spot, both $A_{x}$ and $A_{y}$ should be of the order of $1$. Within these bounds, we choose $W_{t}=1.08\,\mu$m, $A_{x}=1.03$, and $A_{y}=0.89$, which result in the following mechanical frequencies,

[TABLE]

These values agree with the measurements in Ref. Windey et al. .

We now focus on the cavity parameters in Eq. 88, namely the re-normalized detunings $\delta^{\prime}$ and linewidth $\kappa^{\prime}$. Note that, due to the large cavity occupation induced by the high tweezer power, the cavity frequency can be significantly modified by the presence of the NP, thus allowing for the definition of two different detunings. Indeed, the detuning measured without the NP is defined as $\delta_{\rm bare}=\omega_{c}+\Delta_{A}-\omega_{0}$, whereas the detuning measured with the NP inside the cavity is given by $\delta^{\prime}=\delta_{\rm bare}-\Delta_{c}+\Delta_{B1}-2g_{\rm cx}\beta_{x}-2g_{\rm cy}\beta_{y}$, where the expressions for all undefined shifts are given in Appendix B. For the parameters in Table 1, the difference between these two quantities is dominated by $\Delta_{c}$, and can reach a maximum value of $|\delta^{\prime}-\delta_{\rm bare}|\approx 2\pi\times 3.4$ kHz when the NP is placed at an intensity maximum. Since this value is $\sim 10\%$ the mechanical frequency $\Omega_{z}$, a proper identification of which detuning is measured experimentally could be relevant for resolved sideband cooling Aspelmeyer et al. (2014). Following the measurements in Ref. Windey et al. , where the detuning is measured with the NP inside the cavity, in the following we will refer to $\delta^{\prime}$ as the cavity detuning and set it to $\delta^{\prime}=2\pi\times 400$kHz unless stated otherwise 222For the chosen parameters, a smaller detuning would result in the system becoming dynamically unstable, and the NP would abandon the trap. This phenomenon has also been experimentally observed in Ref. Windey et al. .. Finally, regarding the cavity linewidth, we find that $\kappa^{\prime}=\kappa+\kappa_{B1}\approx\kappa$, since the correction to the bare cavity linewidth is negligible, $\kappa_{B1}=2\pi\times(5\cos^{2}\phi)$Hz.

The determination of the mechanical frequencies and the cavity parameters allows us to compute the displacements of the cavity, $\alpha_{c}$, and of the COM motion, $\beta_{j}$, which arise after introducing the trapped NP into the optical cavity. The former can be expressed in terms of the tweezer-induced photonic occupation of the cavity mode, $|\alpha_{c}|^{2}$, whereas the latter quantifies the cavity-induced displacement of the COM equilibrium positions, $2\beta_{j}r_{0j}$. These displacements are shown in Fig. 3 versus the NP position inside the cavity (panel a) and the tweezer polarization angle (panel b). All the displacements vanish when the NP sits at a cavity node ($\phi=\pi/2$), since the cavity field is zero at such position, and when the tweezer polarization is parallel to the cavity axis ($\Theta=\pi/2$), as no photons are scattered into the cavity direction by the dipolar NP. In general, the displacements of the COM equilibrium positions are orders of magnitude above the zero-point amplitudes $r_{0j}\sim 10^{-12}\,\rm m$. Moreover, the shift in the COM position along the tweezer axis ($Z$) has a larger magnitude and negative sign. Such traits can be ascribed, respectively, to the larger mechanical frequencies and to the reduction in the tweezer scattering force along $Z$, caused by the coherent scattering of tweezer photons into the cavity. Finally, note that the cavity might be largely occupied, containing up to $\sim 10^{6}$ photons. This might result in a significant modification of the optomechanical coupling (see Eq. 65).

Using the displacements calculated above, we determine the remaining parameters in the equation of motion. First, we show in Fig. 4 the coherent optomechanical coupling rates between the cavity mode and the COM motion, namely $g_{j}^{\prime}$ in the Hamiltonian Eq. 78. As shown in panel 4a, for the NP sitting at the node of the cavity ($\phi=\pi/2$), the coupling vanishes for the $Z-$COM coordinate and reaches its maximum value for the transverse ($X-Y$) coordinates, while the opposite trend is observed at a cavity anti-node ($\phi=0$). Such behavior evidences the different origin of the couplings for the $Z$ and for the $X-Y$ directions, see e.g. Eq. 45. Indeed, the interaction along the $Z-$ direction is proportional to the electric field intensity of the cavity, whereas the interaction along the $X-Y$ direction is proportional to its derivative. In Fig. 4b, we show the dependence of the optomechanical couplings with the polarization angle of the tweezer. As expected, when the tweezer is polarized along the cavity axis ($\Theta=\pi/2$), no tweezer photons can be scattered into the cavity, and the interaction rates along all three directions vanish. Moreover, as mentioned above, a tweezer polarization purely perpendicular to the cavity axis ($\Theta=0,\pi$) results in the decoupling of the $X-$COM motion, and thus prevents the possibility of 3D cooling. Note that the coupling rates can reach values comparable or even higher than the mechanical frequencies, which may result in the system reaching dynamically unstable regimes Kustura et al. .

Finally, we calculate the remaining rates appearing in the dissipator Eq. 89, for the parameters in Table 1. First, the incoherent coupling rate between the $Z-$COM mode and the cavity is given by $\Upsilon=2\pi\times 82\cos\phi\cos\Theta\text{ Hz}$, and depends on the tweezer power as $P_{t}^{1/4}$. This rate can safely be neglected since it is much smaller than the corresponding coherent coupling rate, i.e., $\Upsilon\ll g_{z}^{\prime}$. Regarding the dissipation associated to the viscous friction term, $\mathcal{D}_{p}$, we find that

[TABLE]

where $P_{\rm mbar}$ is the gas pressure in mbar. From here, we can immediately determine the rates of position localization $\Gamma_{j}$ in Eq. 90, which are independent on the angles $\phi$ or $\Theta$. The three contributions to these rates are shown in Table 2. Note that the displacement noise depends critically on the adimensional PSD, $\sigma_{j}$. In order to reproduce the experimental observations, these parameters are chosen as $\{\sigma_{x},\sigma_{y},\sigma_{z}\}=\{0.67,0.26,18.6\}$, which correspond to displacement PSD values of \{\big{[}S_{xx}^{(d)}\big{]}^{1/2},\big{[}S_{yy}^{(d)}\big{]}^{1/2},\big{[}S_{zz}^{(d)}\big{]}^{1/2}\}\approx\{6,2,500\}\times 10^{-15}\,\text{m}\,\text{Hz}^{-1/2}. These values are of the same order as those estimated in Ref. Windey et al. .

Once the parameters of the equation of motion Eq. 88 have been calculated, we are in a position to study the system dynamics. The expected value of any system observable $\hat{O}$ can be obtained from the density matrix, $\hat{\mu}$, through the usual trace expression, $\langle\hat{O}\rangle(t)=\text{Tr}[\hat{O}\hat{\mu}(t)]$ Breuer et al. (2002). In the present case, there is no need to compute the whole density matrix, since our Master Equation is quadratic in the system creation and annihilation operators. Thus, any information about the state of the system is encoded in the first- and second-order momenta, namely the expected value of single system operators ($\langle\hat{c}\rangle,\langle\hat{c}^{\dagger}\rangle,\langle\hat{b}_{x}\rangle,\langle\hat{b}_{x}^{\dagger}\rangle,...$) and the expected value of all their quadratic combinations ($\langle\hat{c}^{2}\rangle,\langle\hat{c}^{\dagger}\hat{c}\rangle,\langle\hat{c}^{\dagger}{}^{2}\rangle,\langle\hat{c}\hat{b}_{x}\rangle,\langle\hat{c}\hat{b}_{x}^{\dagger}\rangle...$). As detailed in Appendix C, from Eq. 88 we derive the equations of motion for all the above expected values, which form a closed system of $44$ differential equations. After numerically solving such system of equations we can compute any observable as a function of their solutions, for instance the position operator \langle\hat{R}_{j}\rangle=r_{0j}\big{(}\langle\hat{b}_{j}\rangle+\langle\hat{b}^{\dagger}_{j}\rangle\big{)}. Specifically, the temperature of the motion along a given axis is defined as $T_{j}\equiv\hbar\Omega_{j}^{\prime}\langle\hat{b}^{\dagger}_{j}\hat{b}_{j}\rangle/k_{B}$.

Our results for the steady-state temperature along each of the trapping axes as a function of the pressure of the surrounding gas are shown in Fig. 5. As expected, the dependence of the temperatures on the position inside the cavity, $\phi$, mimics that of the optomechanical couplings in Fig. 4. Specifically, the cooling along the $Z-$axis is maximally efficient at a cavity anti-node and vanishes at the cavity node, whereas the cooling along the $X-$ and $Y-$ directions obeys an opposite trend. For the parameters chosen above, the temperatures reached at $P=10^{-5}$mbar and optimum cooling conditions are $T_{x}\approx 100\,{\rm mK}$, $T_{y}\approx 3\,{\rm mK}$, and $T_{z}\approx 60\,{\rm mK}$, which are consistent with the experimental observations Windey et al. . According to our model, the saturation of $T_{z}$ at low pressures is limited by the trap displacement heating $\Gamma_{z}^{(d)}$, whereas at large pressures the gas heating rates $\Gamma_{j}^{(p)}$ dominate for all axes and prevent the motional cooling. Interestingly, we have confirmed that for the chosen parameters the system is very close to a dynamical instability Kustura et al. due to the high tweezer powers, which can result in strong sensitivity to the system parameters. For instance, a slight increase in the size of the NP results in large heating rates at low pressures due to the critical ratio between the coupling rates and the optomechanical frequencies. Regarding an experimental implementation, such dynamical instabilities must be taken into account as they limit the parameter values for which efficient cooling is attainable.

It is common in optomechanics to define the cooling power by means of a cooling rate, which determines both the minimal achievable temperature and the time required to reach it. This quantity can be defined by calculating the power spectral density (PSD) of the COM motion, defined as

[TABLE]

where the sub-index ss denotes the steady state. As detailed in Appendix C, the motional PSD is computed by using the Quantum Regression Theorem Carmichael (1999) which, in its simplest form, relates the two-time correlation function $\langle\hat{q}_{j}(t+\tau)\hat{q}_{j}(t)\rangle_{\rm ss}$ to its value at zero delay, namely $\langle\hat{q}_{j}^{2}\rangle_{\rm ss}$. The calculation of the PSD is straightforward from here, since the zero-delay expected value has been already determined by numerically solving the system equations of motion.

The motional PSD, Eq. 94, is shown in Fig. 6a for $\phi=\pi/4$. As evidenced here, the PSD along a given axis, $j$, shows a series of Lorentzian peaks, the most relevant of which is centered in the vicinity of $\Omega_{j}^{\prime}$. The full width at half maximum (FWHM) of these peaks, shown in Fig. 6b, determines the total rate at which the temperature evolves in time Windey et al. . Note that, as demonstrated in Appendix C, such FWHM is independent of the position localization rates $\Gamma_{j}^{(d)},\Gamma_{j}^{(p)}$, and $\Gamma_{j}^{(r)}$, which only affect the steady state temperature. Furthermore, the remaining source of heating provided by the rate $\gamma$ (Eq. 93) becomes negligible at low pressures. Therefore, we can identify the cooling rates as the low-pressure limits of the curves in Fig. 6, that is, $\gamma_{x}\sim 2\pi\times 30\,\rm Hz$ and $\gamma_{y,z}\sim 2\pi\times 1\,\rm kHz$ for optimal cooling angles $\phi$. These values are consistent with experimental observations Windey et al. and, remarkably, they are of the same order as the rates predicted by adiabatic elimination of the cavity Wilson-Rae et al. (2008) ($\gamma_{j}\approx|g_{j}^{\prime}|^{2}\kappa/(\kappa^{2}+(\delta^{\prime}-\Omega_{j}^{\prime})^{2})$), despite the fact that the assumptions for such elimination, namely $\kappa\gg|g_{j}^{\prime}|$, are not fulfilled in this case.

From the above argumentation and the curves in Fig. 6, we can predict that, at pressures $P\lesssim 5\cdot 10^{-3}\,\rm mbar$ and optimal cooling positions $\phi$, the temperatures will decay to $1/e$ times their initial value in around $\sim 1/\gamma_{x}\sim 5\,\rm ms$ along the $X-$ axis, and $\sim 1/\gamma_{y,z}\sim 0.2\,\rm ms$ along the $Y-$ and $Z-$ axes. This trend is indeed observed in the system dynamics, displayed in Fig. 7. Clearly, the dependence of the optomechanical couplings on the position of the NP inside the cavity not only influences the steady-state temperature, but the whole time evolution. Moreover, the efficiency of the cooling depends critically on $\phi$, with no cooling at all occurring in the less optimal configuration. In order to certify that the coupling to the cavity is responsible for the cooling, we display in Fig. 8 the reheating of the COM motion along the $Z-$axis (similar curves are obtained for every axis). Here, we have initialized the system in its cooled steady state and effectively turned off the cavity cooling by setting the cavity detuning to \delta^{\prime}=2\pi\times 20\,\rm MHz$$\gg\kappa^{\prime},\Omega_{j}^{\prime}. As expected, in this situation the dynamics is cavity-independent and governed only by gas pressure reheating, thus occurring on a timescale $\gamma^{-1}\sim 50\,\rm ms$, in agreement with the experiment Windey et al. . Also in agreement with the above discussion, the reheating dynamics is independent on the position localization heating rates $\Gamma_{j}^{(d)},\Gamma_{j}^{(p)}$, and $\Gamma_{j}^{(r)}$. Finally, note that in this limit, namely $g_{j}^{\prime}\to 0$, the equations of motion can be analytically solved in terms of $\gamma$, the steady-state temperature $T_{\rm ss}$, and the initial temperature $T_{0}$. The resulting analytical expression is represented by the dashed line in Fig. 8, and shows excellent agreement with the full numerical solution. Our results in Figs. 7 and 8 quantitatively agree with experimental measurements Windey et al. .

Because the above defined cooling rates are a function of the optomechanical couplings $g_{j}^{\prime}$, they depend strongly on the cavity detuning $\delta^{\prime}$, a dependence that stems from the tweezer-induced displacement $\alpha_{c}$ (see Eq. 114). In order to characterize this dependence, we show in Fig. 9(a-c) the steady-state temperature along the three motional axes as a function of $\delta^{\prime}$. As expected, the cooling becomes significantly less efficient when the detuning becomes larger than the cavity linewidth, $\delta^{\prime}\gtrsim\kappa^{\prime}\approx 2\pi\times 0.5$ MHz, and completely disappears at around $\delta^{\prime}\approx 10\kappa^{\prime}$. Moreover, the steady-state temperature saturates at small detunings due to the large optomechanical couplings $|g_{j}^{\prime}|\approx\Omega_{j}^{\prime}$. Indeed, as mentioned above, for such large couplings the system is near a dynamical instability regime, where further reduction of the detuning toward resonance with the COM modes pushes the system closer to the unstable regime. This results in additional heating (not shown in Fig. 9) and, eventually, in the loss of the NP from the trap as the system becomes unstable. For the current parameters, such instability might hinder the realization of ground state cooling, since reaching the ground state usually requires the cavity detuning to be in resonance with the mechanical frequency.

The dependence of the cooling rates on the optomechanical coupling rates $g_{j}^{\prime}$ also allows for their tuning by means of the tweezer power. In Fig. 9d, we show the dependence of the steady-state temperature with the tweezer power, at $\delta^{\prime}=2\pi\times 400\,\rm kHz$ and at optimal cooling positions of the NP. Here, we observe an increase of the cooling efficiency with higher powers, as both the COM frequencies and the optomechanical couplings become more relevant with respect to the heating rates. As the power is increased, however, the decrease of the steady-state temperature becomes less pronounced, as the enhancement in the cooling rate starts to be compensated by the displacement heating rates $\Gamma_{j}^{(d)}\propto P_{t}^{2}$. Above a certain power ($P_{t}\approx 2\,\rm W$ for the parameters of Fig. 9), such heating rates become dominant, and the steady-state temperatures become a quadratically increasing function of the tweezer power.

As shown by the results in this section, our model is able to reproduce all the experimental observations reported in Ref. Windey et al. , both qualitatively and quantitatively. This shows the validity of such model to describe cavity-assisted cooling via coherent scattering.

V Requirements for ground state cooling

The parameter values chosen for the above section do not allow for ground state cooling, defined as the reduction of the COM phonon occupation to $\langle\hat{n}_{j}\rangle\equiv\langle\hat{b}^{\dagger}_{j}\hat{b}_{j}\rangle<1$. The first reason is the high values of the gas pressure and the tweezer power, which result in large heating rates $\Gamma_{j}^{(d)}$ and $\Gamma_{j}^{(p)}$. Moreover, the optomechanical couplings $|g_{j}^{\prime}|$, especially along the $Z-$axis, are comparable to the COM frequencies. As a result, a desirable condition for cooling, namely the cavity being in resonance with the COM modes ($\delta^{\prime}=\Omega_{j}^{\prime}$), cannot be reached without the system becoming dynamically unstable Kustura et al. . Another reason hindering ground state cooling is the unresolved sideband regime, i.e. the large linewidth of the cavity with respect to the mechanical frequencies. In this section we use our model, which has been tested by reproducing experimental observations, to determine the conditions for reaching the mechanical ground state along any of the motional axes.

Let us first focus on cooling the motion only along one of the axes, regardless of the dynamics along the remaining two. In order to simplify the physical interpretation of the results, for cooling along the $Y-$($Z-$) axis we take $\{\Theta,\phi\}=\{0,\pi/2\}$ ($\{\Theta,\phi\}=\{0,0\}$), such that the motion along the remaining two directions is uncoupled from the cavity. On the other hand, for cooling along the $X-$ axis, we can only uncouple the $Z-$ motion by choosing $\{\Theta,\phi\}=\{\pi/4,\pi/2\}$. In this section, we aim at modifying the system parameters to achieve ground state cooling in the resolved sideband regime Aspelmeyer et al. (2014), for which three requisites must be fulfilled: first, the total heating rates $\Gamma_{j}=\Gamma_{j}^{(r)}+\Gamma_{j}^{(p)}+\Gamma_{j}^{(d)}$ have to be minimized. Second, the cavity linewidth must be able to resolve the motional sidebands, i.e., $\kappa^{\prime}\ll\Omega_{j}^{\prime}$. Third, the optomechanical coupling must remain small, $|g_{j}^{\prime}|\ll\Omega_{j}^{\prime}$, so that the motional phonons and cavity photons do not hybridize appreciably. Let us focus on these three requisites independently.

We begin by the minimization of the heating rates. We will assume, on the one hand, a reduced gas pressure of $P=10^{-9}$mbar, a value within reach of current ultra-high vacuum technology. On the other hand, since the calculated COM temperatures in the previous section were limited by trap displacement heating, it is necessary to incorporate some isolation in our system, in order to reduce the corresponding noise PSDs, $S_{jj}^{(d)}$ (Eq. 81). In Fig. 10, we display the partial heating rates $\Gamma_{j}^{(r)}$, $\Gamma_{j}^{(p)}$, $\Gamma_{j}^{(d)}$, the total heating rate $\Gamma_{j}$, the frequencies $\Omega_{j}^{\prime}$, and the optomechanical couplings $|g_{j}^{\prime}|$ as a function of tweezer power and for three different values of the displacement PSD. Note that both the recoil heating rate and the displacement heating rate increase with the tweezer power, whereas $\Gamma_{j}^{(p)}$ obeys the opposite behavior, resulting in an optimal tweezer power for which the heating rates are minimized. At high values of the displacement PSD, the heating rate is completely dominated by the displacement heating rate, which become comparable or even larger than the mechanical frequencies. In this limit, ground state cooling becomes impossible. On the other hand, for \big{[}S_{jj}^{(d)}\big{]}^{1/2}\sim 10^{-16}\,m$/\text{Hz}^{-1/2}$, the trap displacement noise is effectively eliminated and the heating is dominated by the interplay between gas pressure and recoil heating. In this regime (lower row in Fig. 10), the total heating rates reach minimum values of $\Gamma_{j}\sim 2\pi\times 100\,\rm Hz$, orders of magnitude smaller than the optomechanical couplings. Our results in Fig. 10 suggest that choosing an adequate tweezer power is crucial for ground state cooling. We remark that all the results in Fig. 10 are practically independent on the cavity linewidth along a wide interval, at least $100\,\text{Hz}\lesssim\kappa^{\prime}/2\pi\lesssim 1\,\rm MHz$.

We now focus on the tuning of the cavity linewidth into the resolved sideband regime. Using our results in Fig. 10, we calculate the steady state occupation along the $X$ and $Y$ axes (the $Z$ axis is more involved, and will be analyzed later) as a function of the cavity linewidth, choosing the values of the tweezer power that minimize the total heating rate. Such occupations are displayed in Fig. 11(a-b) for three different values of the displacement PSD. For high values of the PSD, the system becomes dynamically unstable at small $\kappa^{\prime}$, as the dissipation rate cannot compensate for both the large heating rate and the high ratio $|g_{j}^{\prime}|/\Omega_{j}^{\prime}\approx 1$. On the other hand, for low values of the displacement PSD, ground state cooling is observed along both of the axes, with occupation numbers $\langle\hat{n}_{x}\rangle\approx 0.09$ and $\langle\hat{n}_{y}\rangle\approx 0.14$ respectively. The optimal cooling is not attained for arbitrarily small linewidths, but reaches a minimum at $\kappa^{\prime}\sim|g_{j}^{\prime}|$. As we will demonstrate below, this can be ascribed to the balance between the photon scattering rate into the cavity and the rate of photon loss through the mirrors. We emphasize that the minimal occupations observed in Fig. 11(a-b) are not necessarily the lowest occupations achievable in this system, as we will see below.

The last requirement for cooling is the reduction of the ratio $|g_{j}^{\prime}|/\Omega_{j}^{\prime}$. For the $X$ and $Y$ axes this becomes manifest in the fact that the tweezer powers chosen in Fig. 11(a-b) are not the optimal ones for cooling despite minimizing the heating rate. Indeed, one can check that for higher values of $P_{t}$, the advantageous reduction of $|g_{j}^{\prime}|/\Omega_{j}^{\prime}\sim P_{t}^{1/4}$ overcomes the detrimental increase in the heating rate, and gives rise to lower occupation numbers than the ones shown in Fig. 11(a-b). However, the critical importance of reducing the ratio $|g_{j}^{\prime}|/\Omega_{j}^{\prime}$ is most evidenced when cooling along the $Z-$axis, where this ratio takes very large values (See Fig. 10). Such large values, combined with the recoil heating rate being $\sim 10\times$ larger along the $Z-$axis, make it impossible to cool the motion for tweezer powers below $P_{t}\sim 1-2$W, as the system becomes unstable for any value of $\kappa^{\prime}$ within the resolved sideband regime. However, in order to efficiently cool the motion along the $Z$ axis, the ratio $|g_{j}^{\prime}|/\Omega_{j}^{\prime}$ can be reduced by either increasing the tweezer power or changing the position of the NP inside the cavity, $\phi$. A combination of both methods allows for ground state cooling along the $Z$ axis as well, as shown in Fig. 11c. Due to the larger recoil heating rates, the minimum occupations reached along $Z$, $\langle\hat{n}_{z}\rangle\approx 0.66$, are roughly one order of magnitude larger than along the remaining two axes. Additionally, since recoil dominates the heating for displacement PSDs below $\sim 10^{-15}\,$m$/\sqrt{\text{Hz}}$ (see Fig. 10), reducing such PSD below that value has little effect on the phonon occupation. Remarkably, for the parameters of Fig. 11c we find a regime of uncertain stability (shaded area), where some solutions of the polynomial equation for $\beta_{y}$ (see Appendix A) yield stable equations of motion and some do not. Physically, this is a regime where the dynamics is extremely sensitive to small perturbations, as evidenced by the fact that the dependence of the equations of motion on $\beta_{y}$ is weak. Discerning whether the system is stable or not in such regime would likely require to refine the treatment of the problem, including previously disregarded terms in the Hamiltonian. This, however, lies beyond the scope of the present work, as such uncertain stability regimes are seldom encountered in the parameter space.

Regardless of the cooling axis, all the steady-state occupation curves in Fig. 11 follow a general trend, reaching a minimum at $\kappa^{\prime}\approx|g_{j}^{\prime}|$. This behavior can be explained from our equations of motion in the low-pressure limit $\gamma\ll\Omega_{j}^{\prime},\kappa^{\prime},|g_{j}^{\prime}|$, where, assuming only one motional axis is coupled to the cavity, its steady-state occupation along the chosen direction can be analytically calculated as

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

and we use the simplified notation $g\equiv|g_{j}|$, $\delta\equiv\delta^{\prime}$, $\Gamma\equiv\Gamma_{j}$, $\kappa\equiv\kappa^{\prime}$, and $\Omega\equiv\Omega_{j}^{\prime}$. Assuming the cavity is tuned at $\delta=\Omega$, we can find three simple limits for the above expression:

[TABLE]

The above limits are easily understood as the ratio between the heating rate and a cooling rate, the latter of which depends on two factors: on the one hand, the ability of the NP to scatter photons inside the cavity, $|g_{j}^{\prime}|$, and on the other hand, the ability of the cavity to dissipate such photons, $\kappa^{\prime}$. For $\kappa^{\prime}\ll|g_{j}^{\prime}|$, the cooling power is limited by the small cavity loss, which is not enough to efficiently reduce the cavity occupation. This results in the increase in steady-state occupation at low $\kappa^{\prime}$ observed in Fig. 11. In the opposite limit, $\kappa^{\prime}\gg|g_{j}^{\prime}|$, the cooling rate is given by $|g_{j}^{\prime}|^{2}/\kappa^{\prime}$, which coincides with the rate obtained through adiabatic elimination Wilson-Rae et al. (2008). In this case, the cooling efficiency is limited by the low value of the optomechanical coupling and, for the system under study, ground state cooling becomes ultimately inefficient as $\kappa^{\prime}$ increases. Finally, in the intermediate regime $\kappa^{\prime}\approx|g_{j}^{\prime}|$ the cooling is optimized, as the rate at which photons are scattered into the cavity is equal to the rate at which they are lost through the cavity mirrors. This explains the optimal values of $\kappa^{\prime}$ found in Fig. 11, and gives an order-of-magnitude theoretical estimation for the lowest achievable occupations, namely $\langle\hat{n}_{x},\hat{n}_{y},\hat{n}_{z}\rangle\approx 2\Gamma_{j}/|g_{j}^{\prime}|=\{0.02,0.02,0.2\}$.

So far, we have demonstrated ground state cooling along each of the three motional axes, in order to have a deeper understanding of the underlying physics and relevant parameters. However, the results obtained above (e.g. in Fig. 11) are ill suited from an experimental point of view, as the motion along the remaining two axes is either in an unfavourable regime for cooling or even completely uncoupled from the cavity. As a result, the occupations along these two axes usually remain at their room temperature values, and in some cases drastically increase as the system is close to a dynamically unstable regime. From a practical perspective, it is convenient to find a configuration in which one of the motional axis is cooled to the ground state while the remaining two are cooled at least below $\sim 1$K, such that their motion does not probe the non-linearity of the trapping potential. This can be achieved by taking our results in Fig. 11 as a starting point in parameter space, and modifying the angles $\{\Theta,\phi\}$ in order to couple the remaining axes to the cavity. As shown in Fig. 12, this method leads to regimes in which near-to-ground state cooling is attained along the three axes. Note that the $X-$ and $Y-$ motional directions can be simultaneously cooled to the ground state (Fig. 12a). In both panels of Fig. 12, the axes that do not reach the ground state are cooled below $10\,\mu$K. These results show that ground state cooling along each of the motional directions is experimentally feasible.

All the results in this section highlight the relevance of isolation from trap displacement noise for reaching the ground state, and provide strict bounds for the displacement noise level above which ground-state cooling becomes impossible. However, note that reducing the trap displacement noise might not be sufficient for reaching the ground state, as the system might be subject to some additional heating sources that, similar to the trap displacement, have no effect on the motional PSD. Thus, a proper identification and suppression of the relevant heating sources in a given experimental setup is crucial for ground state cooling.

VI Conclusion

In this work, we have developed a theoretical model for the COM cooling of a levitated NP via coherent scattering into a cavity. We have done so by, first, obtaining a quadratic Hamiltonian for the NP and the electromagnetic field and, second, by tracing out the free field to obtain a Master Equation for the cavity and COM degrees of freedom. Heating due to gas pressure and trap displacement has been included. Our model has been shown to reproduce recent experimental observations, and is thus revealed as a useful tool for levitodynamics.

We have also used our model to study the possibility of ground state cooling of a levitated NP. The relevance of the heating due to trap displacement, which is often ignored, has been emphasized and evidences the necessity of isolation in order to reach the ground state. Our model demonstrates the potential of current experimental setups to reach the mechanical ground state of the COM motion of a NP. The requirements for ground state cooling given in this paper could bring about the quantum regime of levitodynamics, where the full quantum theory we have developed represents a necessary and powerful theoretical tool.

Acknowledgements.

C. G. B. acknowledges funding from the EU Horizon 2020 program under the Marie Skłodowska-Curie grant agreement no. 796725. This research was partly supported by the Swiss National Science Foundation (no. 200021L$\_$169319) and ERC-QMES (no. 338763). R. R. acknowledges funding from the EU Horizon 2020 program under the Marie Skłodowska-Curie grant agreement no. 702172.

We thank Mathieu L. Juan for insightful discussions.

Appendix A Calculation of the displacements

In this Appendix we find the equations for the coefficients $\alpha_{\mathbf{k}},\beta_{j},\gamma_{\alpha}$ and solve them. We start by substituting the displaced operators 59, 60, 61, and 62 into the Hamiltonian Eq. 55. The transformed Hamiltonian reads $\hat{H}_{0}+\Delta\hat{H}$, where $\hat{H}_{0}$ contains no displacements:

[TABLE]

The non-quadratic term in the first line of the above equation, $\sim\hat{c}^{\dagger}\hat{c}(\hat{b}_{j}+\text{H.c.})$, can be safely neglected, since in the displaced frame the cavity is in the vacuum state and thus the contribution of this term is of order $g_{cj}\sim 2\pi\times 1$Hz. On the other hand, $\Delta\hat{H}$ contains all the single-operator terms:

[TABLE]

We can now make $\Delta\hat{H}=0$ by setting the coefficients of each creation operator to zero, obtaining the following system of equations:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where we have used the fact that the coefficients $\beta_{j}$ can be chosen real without loss of generality. From Eqs. 106 and 107, we immediately obtain Romero-Isart et al. (2011a)

[TABLE]

[TABLE]

where P.V. stands for the Cauchy principal value. We now introduce the above results into Eq. 105 to solve for $\alpha_{c}$. In the process, we must calculate two integrals. The first one is given by

[TABLE]

The principal value integral above can be carried out by using the fact that $\gamma(\omega)\approx\sqrt{\kappa/\pi}$ on a broad range of frequencies. After changing variables from $\omega$ to $\omega-\omega_{0}$, we can approximate the second term of the above equation by

[TABLE]

The second integral to calculate is given by

[TABLE]

where we have transformed the sum into an integral, $\sum_{\varepsilon\mathbf{k}}\to\sum_{\varepsilon}\mathcal{V}/(2\pi)^{3}\int\text{d}\Omega_{\mathbf{k}}\int dkk^{2}$, $\Omega_{\mathbf{k}}$ being the solid angle in wavevector space. Using the closure relation of the polarization vectors,

[TABLE]

we can perform explicitly all the angular integrals above, and solve for the coefficient $\alpha_{c}$ as

[TABLE]

where we have defined for compactness

[TABLE]

[TABLE]

[TABLE]

In Eq. 116, we have introduced a high-energy cutoff $K_{c}$ to prevent the integral from diverging. Physically, this cutoff restricts the integration domain to values of $k$ for which the long-wavelength approximation (and thus our model) remains valid. In our calculations we choose $K_{c}=0.1/R$.

The last step is to introduce Eqs. 108, 109, and 114 into the three equations for the coefficients $\beta_{j}$. In the same fashion as above, we first calculate the sums

[TABLE]

where we define

[TABLE]

The resulting system of equations for $\beta_{j}$ can be further simplified. Indeed, since the expression for $\alpha_{c}$ is linear in $\beta_{z}$, we can express this coefficient in terms of $\beta_{x}$ and $\beta_{y}$ as

[TABLE]

where we define

[TABLE]

Introducing all the above identities into the equations for $\beta_{x}$ and $\beta_{y}$ we obtain the following system of equations,

[TABLE]

[TABLE]

Both the above equations are polynomial equations of degree $5$ in the variables $\beta_{x}$ and $\beta_{y}$. However, they can be reduced by introducing the explicit expressions for $g_{\rm cx}$ and $g_{\rm cy}$ (Eqs. 34 and 35). This allows us to take the ratio of both equations and show that, for $\Theta\neq\pi/2$, the two variables are easily related through

[TABLE]

After introducing this solution into either Eq. 122 or 123, we find a fifth degree polynomial equation for the remaining coefficient, $\beta_{y}$, which has to be solved numerically. Note that imaginary solutions for such polynomial equation are disregarded as $\beta_{y}\in\mathbb{R}$ by assumption. This univocally determines $\beta_{y}$ since in the stable regions the polynomial usually has a single real solution. After numerically obtaining $\beta_{y}$, the remaining coefficients can be easily computed using Eqs. 108, 109, and 114, 120, and 124.

Appendix B Tracing out of the continuum reservoirs

Our starting point in this section is Eq. 76 and, more specifically, the second term

[TABLE]

As mentioned in the main text, the interaction potential contains three terms, resulting in $9$ independent contributions for the above commutator. However, these interaction terms correspond to two reservoirs considered independent, and therefore, in a thermal state, the cross terms between them will vanish after tracing,

[TABLE]

We thus reduce the problem to calculating four independent dissipators,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Let us illustrate the calculation for the cavity dissipation given by $D_{A}$. The potential in the Interaction Picture is given by

[TABLE]

where for simplicity we have defined

[TABLE]

We now introduce this potential into Eq. 127 and take the trace over the bath modes explicitly. In order to do so, we note that in a thermal state,

[TABLE]

[TABLE]

where $\bar{n}_{\omega}$ is the mean occupation of the bath at frequency $\omega$, given by the Bose distribution. After tracing out we obtain

[TABLE]

with

[TABLE]

We now perform the integral in $s$ by using the identity

[TABLE]

obtaining

[TABLE]

in terms of the Lindblad dissipator Breuer et al. (2002); Gardiner and Zoller (2004)

[TABLE]

where the curly brackets denote the anti-commutator. To obtain the result Eq. 138, we have used the same approximations used for Eqs. 110 and 111. In order to estimate the Lamb shift of the cavity in Eq. 138, we assume that, as most physical spectral densities, the coupling $\gamma$ satisfies $\gamma(\omega)\approx 0$ at low frequencies. Then, in the integration, the only relevant contribution will arise from frequencies close to $\omega=\delta_{0}+\omega_{0}\approx\omega_{0}$. In this interval, we can approximate $\bar{n}_{\omega}\approx\exp(-\hbar\omega_{0}/k_{B}T)$ since $\exp(\hbar\omega_{0}/k_{B}T)\gg 1$. Making also use of the identity $\gamma(\omega_{0})=\sqrt{\kappa/\pi}$, we find after changing variables

[TABLE]

with $u_{0}=\hbar(\omega_{0}+\delta_{0})/k_{B}T\gg 1$. The integral above is named the Exponential Integral Wong (2001) $Ei(u_{0})$, and its asymptotic expansion for $u_{0}\gg 1$ is $E_{i}(u_{0})\approx e^{u_{0}}[u_{0}^{-1}+O(u_{0}^{-2})]$. This identity allows us to approximate in this limit

[TABLE]

where we have assumed $\omega_{0}\gg\delta_{0}$, as is usual in cavity cooling experiments. Using this result, and noting again that, since the occupation of the environment at the frequency of the tweezer is negligible, i.e., $\exp(\hbar\omega_{0}/k_{B}T)\gg 1$, we can finally write the dissipator $D_{A}$ in Eq. 138 as

[TABLE]

The above dissipator contains a frequency shift of the cavity and a Lindblad-type dissipation, which adds an exponential decay of the cavity photon occupation at rate $2\kappa$. Note that this frequency shift is not measurable independently or, in other words, a measurement of the frequency of the cavity in the absence of the NP would yield $\omega_{c}+\Delta_{A}$, this being the true frequency of the bare cavity.

The remaining contributions are calculated in the same fashion. First, the contribution $D_{B1}$ can be shown to have a similar form as Eq. 142, as the operator $\hat{V}_{B1}$ also couples the cavity modes to the free EM field. This results again in both a frequency shift of the cavity and a Lindblad dissipation, this time induced by the presence of the NP:

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

By substituting above the expression for the couplings, Eq. 71, and transforming the sum into an integral, we can determine the above coefficients as

[TABLE]

[TABLE]

[TABLE]

For usual parameter values Windey et al. , $\bar{n}_{\omega_{0}+\delta_{0}}\sim 10^{-13}$ and $\kappa_{B1}\approx 2\pi\times 100\cos^{2}\phi\text{Hz}$, which allows us to neglect $\kappa_{B1}^{\prime}$ and express

[TABLE]

The final contribution $D_{B1}$ thus reads

[TABLE]

Regarding the contribution $D_{B2}$, its calculation involves the approximation that both the coupling rates $g_{j\varepsilon\mathbf{k}}$ and the occupation numbers $\bar{n}_{\omega}$ are smooth enough functions to approximate their values at $\omega_{0}\pm\Omega_{j}$ by their values at $\omega_{0}$. This approximation is extremely good in the regime $\Omega_{j}\ll\omega_{0}$, and has been checked numerically to be accurate up to $8$ digits in the final dissipation rates, for the parameters in Ref. Windey et al. . By further neglecting the terms $\bar{n}_{\omega_{0}}\approx 0$ as done above, the resulting contribution reads

[TABLE]

with $\hat{q}_{j}\equiv\hat{b}_{j}^{\dagger}+\hat{b}_{j}$, dissipation rates

[TABLE]

and frequency shifts

[TABLE]

Although Eq. 152 has an apparent off-diagonal structure, by explicitly performing the angular integrals in the above rates we find that

[TABLE]

Thus, we might write the contribution $D_{B2}$ in the more compact way

[TABLE]

which includes a frequency shift of the COM oscillations, and a dissipation of the Brownian Motion form Gardiner and Zoller (2004) which we identify with the photon recoil heating and thus label with the superscript $(r)$. This kind of noise, different from the Lindblad dissipation appearing in $D_{A}$ and $D_{B1}$, induces heating without modifying the rate at which the COM occupations decay. The recoil heating rates are given by

[TABLE]

and as expected are larger for the motion along the tweezer axis. The frequency shifts, on the other hand, read

[TABLE]

Finally, let us focus on the contribution $D_{B12}$. Under the same approximations undertaken above, this term can be shown to be

[TABLE]

i.e., it contains a dissipative interaction between COM and cavity, with rate

[TABLE]

and a coherent interaction at rate

[TABLE]

where we have defined the adimensional constant

[TABLE]

The contribution $D_{B12}$ is therefore simplified to

[TABLE]

For typical experimental values Windey et al. , $C_{0}\approx 10^{-13}$ so $\Upsilon\sim 2\pi\times 100$Hz.

After combining all the terms of this section into the von Neumann equation Eq. 76 and transforming back into the Schrödinger Picture, we find that we can regroup all the terms into the equation

[TABLE]

The coherent part $\Delta\hat{H}$ contains all the coherent interactions and frequency re-normalizations,

[TABLE]

The first two terms in this equation can be reabsorbed into the frequencies and the couplings of $\hat{H}_{S}$, i.e., we define

[TABLE]

[TABLE]

For the third term, we can do the same by redefining the phonon operators, i.e., we extract from $\hat{H}_{S}$ the Hamiltonian of the oscillators and add it up to this term, writing

[TABLE]

From this equation, we can see that this term results in a re-normalization of the COM oscillator frequencies, which now read

[TABLE]

We have also defined new bosonic operators $\hat{\tilde{b}}_{j}$ relative to this frequency, so that the physical meaning of the position and momentum operators remains unchanged. In terms of these operators, every time a term $\hat{b}^{\dagger}_{j}+\hat{b}_{j}$ appears in the Hamiltonian, an extra multiplying factor

[TABLE]

must be added. This substitution allows us to write the full coherent part as

[TABLE]

where $g_{j}^{\prime}\equiv g_{tj}\chi_{j}$, and we have explicitly written the expression for $\hat{V}_{0}$. In the main text, we remove the tildes from the operators in order to simplify the notation.

Regarding the dissipators in Eq. 164, we can write them as

[TABLE]

where we have defined $\kappa^{\prime}\equiv\kappa+\kappa_{B1}$, $\Upsilon=\Upsilon_{z}\chi_{z}$, and the recoil heating $\Gamma_{j}^{(r)}\equiv\Gamma_{j}\chi_{j}^{2}$. Again, in the main text we remove the tildes from the operators to obtain the Master Equation 77.

Appendix C Equations of motion and motional PSD

In this appendix we detail how to solve the equation of motion 88 for the compound cavity+COM system,

[TABLE]

where $\hat{H}_{S}^{\prime}$ is given by Eq. 78, and the dissipator reads

[TABLE]

Here, we have neglected the incoherent interaction $\propto\Upsilon$ since, in the cases we consider in the main text, it is negligible. Otherwise, including this term is straightforward. In principle, one could aim at solving the above equation for the full density matrix $\mu$, but this is not necessary in the present case. Indeed, we note that the above Master Equation is quadratic, i.e. it only contains quadratic combinations of creation and annihilation operators. This property implies that the Gaussian character of any initial state is preserved, i.e., if the initial state is Gaussian, any system observable will be determined exclusively by the first- and second- order momenta ($\langle\hat{c}\rangle,\langle\hat{c}^{\dagger}\rangle,...$ and $\langle\hat{c}^{\dagger}\hat{c}\rangle,\langle\hat{c}^{2}\rangle,...$, respectively) at all times Gardiner and Zoller (2004). We thus can derive a closed system of differential equations containing only the expected values of single creation/annihilation operators, and their quadratic combinations.

In order to obtain the desired system of equations, we start by determining the time evolution of the expected value for a general operator $\hat{O}$ in the Schrödinger Picture, given by

[TABLE]

Introducing Eq. 173 and using the cyclic invariance of the trace, we can write

[TABLE]

This relation allows us to substitute any operator $\hat{O}$ and get its equation of motion, by using the fundamental commutation relations $[c,c^{\dagger}]=1$ and $[b_{i},b^{\dagger}_{j}]=\delta_{ij}$. In this way, we obtain two independent systems of equations, namely one for the single-operator expected values, and a second for the quadratic combinations. The first one is given by the equations

[TABLE]

[TABLE]

and their Hermitian conjugates, which together form an $8\times 8$ system of equations. Note that all the expected values decay with time and, since there is no independent term to these equations, the steady-state value of the single-operator expected values will be zero. On the other hand, the second system of equations is given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and their Hermitian conjugates (here, $C.c.$ stands for complex conjugate). This second system contains $36$ equations, and has to be solved numerically.

The above equations of motion allow us to calculate the power spectral density of the COM motion, defined as

[TABLE]

where “ss” refers to the steady state. In order to calculate the above PSD, we first split the integral into two parts, and use the time-translation invariance of the steady state, i.e., $\langle\hat{A}(t)\hat{B}(t+\tau)\rangle_{\rm ss}=\langle\hat{A}(t-\tau)\hat{B}(t)\rangle_{\rm ss}$, to write in a more convenient notation

[TABLE]

with $S_{jj}^{-}(\omega)=[S_{jj}^{+}(\omega)]^{*}$, and

[TABLE]

Expressed in this form, the PSD can be calculated directly by applying the Quantum Regression Theorem. We start by defining the vector containing the single operators of our equation of motion,

[TABLE]

whose expected value is governed by the matrix equation

[TABLE]

where the matrix $M_{0}$ can be directly extracted from Eqs. 177 and 178. According to the Quantum Regression Theorem Carmichael (1999), any two-time correlation function involving a general operator $\hat{A}$ obeys the equation

[TABLE]

which can be easily solved. For the steady-state values appearing in Eq. 187, we find

[TABLE]

where we have defined the eigenvalues and eigenvectors of the matrix $M_{0}$,

[TABLE]

and the vector of coefficients $\mathbf{c}_{0}$ obeys the equation

[TABLE]

Here, $\Lambda$ is the matrix whose columns are the eigenvectors $\boldsymbol{\lambda}_{j}$. Although the PSD has to be calculated numerically, two important properties can be noted right away. First, from Eq. 191 we see that any two-time correlation function can be written as a sum of exponential functions which, after integration in $\tau$, yield a sum of Lorentzian profiles. Indeed, if we define the position of the operators $\hat{b}_{j}$ and $\hat{b}^{\dagger}_{j}$ inside the vector $\hat{\mathbf{v}}$ as $p_{j}$ and $\bar{p}_{j}$ respectively, we can easily obtain

[TABLE]

where the vector $\mathbf{e}_{l}$ represents the length $8$ vector with components $\mathbf{e}_{l}^{(k)}=\delta_{lk}$. From the above equation, we can conclude that, in general, the motional PSD will be composed by a set of Lorentzian peaks centered at $\omega\propto\pm\text{Im}\lambda_{l}$ and with a linewidth $\Delta\omega\propto\text{Re}\lambda_{l}$. Furthermore, note that the only system parameters determining the positions and widths of the peaks are the entries of the matrix $M_{0}$ or, equivalently, the coefficients of Eqs. 177 and 178. In these equations, the position diffusion coefficient $\Gamma_{j}$ of the dissipator does not appear and, henceforth, these two parameters do not depend on $\Gamma_{j}$. However, the steady-state values of the two-operator products do depend on $\Gamma_{j}$, which thus influences the relative weight of the PSD peaks.

Appendix D Summary of system parameters