Coherent-scattering two-dimensional cooling in levitated cavity optomechanics

The strong light-matter optomechanical coupling offered by coherent scattering set-ups have allowed the experimental realization of quantum ground-state cavity cooling of the axial motion of a levitated nanoparticle [U. Delić et al., Science 367, 892 (2020)]. An appealing milestone is now quantum two-dimensional (2D) cooling of the full in-plane motion, in any direction in the transverse plane. By a simple adjustment of the trap polarization, one obtains two nearly equivalent modes, with similar frequencies ωx ∼ ωy and optomechanical couplings gx gy—in this experimental configuration we identify an optimal trap ellipticity, nanosphere size, and cavity linewidth which allows for efficient 2D cooling. Moreover, we find that 2D cooling to occupancies nx + ny 1 at moderate vacuum (10−6 mbar) is possible in a “Goldilocks” zone bounded by √ κ /4 gx, gy |ωx − ωy| κ , where one balances the need to suppress dark modes while avoiding far-detuning of either mode or low cooperativities, and κ ( ) is the cavity decay rate (motional heating rate). With strong-coupling regimes gx, gy κ in view one must consider the genuine three-way hybridization between x,y and the cavity light mode resulting in hybridized bright/dark modes. Finally, we show that bright/dark modes in the levitated set-up have a simple geometrical interpretation, related by rotations in the transverse plane, with implications for directional sensing.


I. INTRODUCTION
The coupling between light and matter has led to major milestones in physics, from the Michelson-Morley experiment [1] to the detection of gravitational waves by the LIGO collaboration [2]. The basic scheme relies on the light acting as a probe-offering exceptional sensitivities-which is now routinely done in state-of-the-art optomechanical systems with high-quality mirrors. The latter are themselves interesting systems and have led to the field of cavity optomechanics [3]. A mirror with a motional degree of freedom cooled to its ground state is of particular interest as it becomes a quantum sensor and thus can be used as a detector of weak forces and as a probe of the quantum-to-classical transition [4].
On the other hand, quantum features of an object in all three spatial dimensions-with applications ranging from quantum foundations to directional sensing-can be explored using an optically levitated nanoparticle [5][6][7][8][9]. Initial experimental efforts have been hindered by technical difficulties of stable trapping in high vacuum [10,11] and several implementations have been considered such as hybrid tweezer-cavity traps [12,13], electro-optical traps [14,15], and trapping in the near field of a photonic crystal [16]. Recently, a three-dimensional (3D) coherent scattering (CS) setup was introduced to levitated cavity optomechanics [17,18] using methods adapted from atomic physics [19][20][21][22][23]. In contrast to experiments that consider dispersive coupling, here the cavity is driven solely by the dipole radiation of the optically trapped silica particle. Due to the tight focus of the optical tweezer this scheme yields unprecedentedly high optomechanical coupling rates, which subsequently enabled ground-state cooling of the motion along the cavity axis and thus opened the door to quantum levitated optomechanics [24].
For the purpose of prolonging available free fall experiment times [25], an important future milestone for the coherent scattering setup is the simultaneous ground-state cooling of all three translational degrees of freedom.
There is strong motivation, however, for investigating the cavity cooling of 2D motions in the tweezer transverse plane (x-y plane): The frequencies are similar ω x ≈ ω y and for suitable experimental parameters, g x g y ≡ g = (g x + g y )/2, so both may be strongly coupled to the light. In contrast, g z ∼ 0 at the favorable configuration of trapping at the node of the cavity that minimizes optical heating, a key advantage of CS set-ups. In addition, the frequency of the z motion is typically ω z ω x,y , far from the optomechanical resonance. Finally, since the z motion can be cooled using feedback cooling [26,27] strong 2D x-y cavity cooling lays a path to 3D motional ground-state cooling.
We exploit here recently developed theoretical expressions for the full 3D coherent scattering problem [28]. However, although the numerics presented below are fully 3D, for anal-ysis and understanding we obtain and discuss the 2D x-y problem. We optimize 2D cooling with respect to particle size, trap frequencies, tweezer polarization orientation, as well as detuning between the tweezer frequency and cavity resonance. For readily achieved experimental pressures of p = 10 −6 mbar we identify a "Goldilocks" region √ κ /4 g x , g y |ω x − ω y | κ, where κ ( ) is the cavity decay rate (heating rate). This set of requirements minimizes the formation of decoupled dark modes and optimizes 2D cooling for |ω x − ω y | ∼ κ/2 by using a particle of radius ∼80 nm. While bright/dark modes have been previously investigated in optomechanical systems [29] in the levitated system they have a geometric interpretation in terms of the rotation of the x and y axes of the oscillator, with potential implications for directional sensing. The importance of nondegenerate mechanical frequencies ω x = ω y for successful 2D cooling is a well-known fact in experiments with trapped ions and atoms [22,30].
This work is organized in the following way. We start by reviewing the coherent scattering setup and introducing the relevant experimental parameters (Sec. II). We then illustrate how mechanical modes hybridize with the optical mode, resulting in the formation of bright/dark modes and three-way mixing. In particular, we show how dark modes distort the relation between the displacement and heterodyne spectra, making in general thermometry and sensing nontrivial (Sec. III). In the central part we give a detailed analysis of 2D cooling and discuss how to perform thermometry (Sec. IV) as well as identify the best parameters for 2D cooling by numerically scanning the experimental parameters (Sec. V). We conclude by laying down a path for 3D motional groundstate cooling in the levitated optomechanics-in particular, how 2D cavity-optomechanical cooling can be combined with feedback cooling to achieve the 3D motional ground state of the optically levitated system (Sec. VI).

II. EXPERIMENTAL SETUP
We consider the 3D coherent scattering setup illustrated in Fig. 1(a). The nanoparticle is trapped in an optical tweezer and positioned inside an optical cavity-the cavity is driven entirely by the tweezer light scattered off the nanoparticle, namely coherent scattering with a pattern shown in Fig. 1(a). Such a scheme offers unique versatility with respect to the customary cavity optomechanical system, since the nanoparticle can be placed at any point inside the cavity by displacing the tweezer trap. Here we will consider the case when the nanoparticle is close to a cavity node x (c) 0 ∼ λ/4 (λ: laser wavelength), where the strongest coupling to the nanoparticle x and y motions is achieved. In addition, deleterious effects of cavity photon scattering and recoil heating are minimal.
Linearization of the effective potentials in the CS set-ups [17,18] has shown that where E d = α c tw sin(θ tw ) 2h is the driving amplitude of the cavity, α = 3 0 V s R −1 R +2 ( 0 is the permittivity of free space, V s is the volume of the nanoparticle, R is the relative dielectric permittivity,), c = hω c 2 0 V c (ω c is the cavity frequency and V c is Schematic of coherent-scattering experiments: An adjustment of the tweezer polarization (θ tw = π/4, 3π/4) yields two equivalently coupled x, y mechanical modes, |g x | |g y | g. (b) We find there is a "Goldilocks" region (orange) for 2D ground-state cooling, i.e., n x + n y 1, illustrated for the set-up of Ref. [24] but with θ tw = π/4. The optimal region lies below the blue curve to avoid the formation of decoupled dark modes (i.e., |ω x − ω y |) g) and is bounded from above by the constraint to avoid far-detuning (κ |ω x − ω y |) and from below by the regime of weak quantum cooperativities [i.e., C = 4g 2 /(κγ n B ) 1, where γ (n B ) is the gas damping (mean thermal occupancy)]. Red lines correspond to different particle sizes and indicate R ∼ 80 nm is optimal. the cavity volume), k = ω c /c, and tw = 4P tw w x w y π 0 c (P tw is the tweezer power and w x , w y is the waist of the Gaussian beam along the x or y axis, respectively).
The angle θ tw between the tweezer polarization axis (y axis) and the cavity symmetry axis (x (c) axis) can be arbitrarily set to tune coupling rates g x and g y . Motional 1D ground-state cooling (of a single mechanical degree of freedom) along x (c) has been recently achieved by setting θ tw ∼ π/2. In this case thetweezer-based coordinates (x, y) and the cavity-based coordinates [x (c) , y (c) ] identify the same point in the 2D plane orthogonal to the tweezer symmetry axis, z. However, one obtains g x ≈ g y for θ tw = π/4 and that is the regime we consider for 2D cooling.

III. BRIGHT/DARK MODES
Avoided crossings are ubiquitous in physics. For example, two classical (or quantum) modes, say,x andŷ, approaching an energy degeneracy are universally described by a  Hamiltonian represented in terms of Pauli matrices: where g is the coupling. At the degeneracy, ω x = ω y , the normal modes of the system correspond to the eigenmodes ofσ x and thus have the maximally hybridized formx ±ŷ. The corresponding frequencies are perturbed by ±g.
In regimes of negligible dissipation g κ, γ , the usual optomechanical interaction corresponds to a two-mode avoided crossing. The two-mode crossing was demonstrated experimentally in optomechanics [31] where it was observed that the cavity light mode and mechanical modes were hybridized, with the characteristic 2g splitting, and a more recent study investigated levitated nanoparticles [32].
In the present case, we have a three-mode avoided crossing. These are less frequently encountered but may be discussed in a similar way. In our present case, two mechanical modes,x andŷ, are coupled to the optical mode,Ẑ L , according to the usual position-position form (see Sec. IV for more details): Representing the modes as a vector, [xẐ Lŷ ] , and Eq. (3) in matrix form we write: where we have also included the two mechanical frequencies, ω x , ω y , and the detuning, − . We first consider equal couplings g x = g y ≡ g, and set − = (ω x + ω y )/2. Neglecting the term 1 4 ω x +ω y 2 I, where I is the identity, we write: where nowŜ x ,Ŝ z are spin 1 matrices (divided byh). The associated anticrossing has an enhanced width of 2 √ 2g. For the degenerate case, ω x = ω y , the eigenmodes are simply the textbook eigenvectors of theŜ x matrix: In that case, there are two three-way hybridized "bright" eigenmodes (B ± ) with eigenvalues ± √ 2g, and a two-way hybridized "dark" eigenmode (D) with eigenvalue zero: While the pedagogic Eq. (5) was not used to compute the realistic system eigenmodes, it illustrates the significance of lifting the frequency degeneracy: ω x = ω y introduces aŜ z component that mixes the bright dark modesD,B ± . The case g x = g y (but ω x ∼ ω y ) does not eliminate the bright-dark mode structure: It simply alters the dark eigenvectors at the center of the crossing tô (which is still a dark mode with eigenvalue 0) while the bright modes areB ± = 1 . However, we see that as g y → 0,D →ŷ. This is the quasi-1D dynamics analyzed in Ref. [28] (using θ tw = 85 • ). In this limit, theŷ mode is "dark" simply because it is very weakly coupled.
The true eigenmodes of the coherent scattering set-up for arbitrary were computed numerically from the equations of motion of the system. The Hamiltonian for the reduced 2D case can be put in the following form: where ω x , ω y are the frequencies of the two harmonic motions; x,ŷ (p x ,p y ) are the position (momentum) observables; g x , g y denote the optomechanical couplings [see Eq. (1)]; andẐ L (P L ) denote the amplitude (phase) quadrature of the intracavity field.
The resulting equations of motion, including dissipation and Gaussian noise baths acting on each mode, yield a set of linear coupled equations which are represented in the wellknown form:Ẋ where A is a drift matrix that includes dissipative terms and jth element of the vector (AX . For our discussion we neglect the z motion (but it is included in quantum numerics) so can consider A as a 6 × 6 matrix. X is the vector of the mechanical and optical modes, X = (x,p x ,ŷ,p y ,Ẑ L ,P L ) T , √ represents the diagonal matrix of damping coefficients, while X in (t ) represents the Gaussian noises (gas collisions and optical shot noises).
To obtain the classical normal modes and frequencies of the system, we calculated the eigenvalues and eigenvectors of A as a function of the optical detuning for the 2D g x = g y case. It was sufficient for our classical analysis to consider the case with dissipative terms set to zero ( √ = 0). Such eigenmodes can be represented as unit vectors in the space spanned by the tweezer modesx,ŷ, andẐ L , i.e., as spherical polar angles on a Bloch sphere [see Fig. 2(a)]. As the detuning is varied from = −∞ to ∼ 0 we represented the resulting trajectory traced by each eigenmode using the corresponding spherical polar angles- Fig. 2(b) plots the trajectories as a function of detuning for two realistic scenarios. In the first case (orange lines) we employ the experimental parameters of Ref. [24] and thus a more elliptical tweezer trap w x = 0.6 μm, w y = 0.705 μm where tweezer frequencies differ by about 16%. In the second case (red lines) we set w x = 0.68 μm, so now the frequencies are near-degenerate (near circular trap) and differ by about 3.5%.
at the centre of the crossing = ω x ω y . For the elliptical trap of the experiments [24], however, the "dark" mode (top panel) still rotates to φ = −π/4 at the center of the crossing but mixes appreciably with the optical mode and θ 0.6π . The more general case of g x = g y (and arbitrary detuning) does not eliminate the bright-dark mode structure. In contrast lifting the ω x ∼ ω y degeneracy has a pronounced effect-the bright/dark modes mix and very different trajectories are obtained. Ultimately, this would lead to a decoupling to two independent level crossings, with the associated disadvantage that both modes might no longer be resonant simultaneously-this case depends on κ and is investigated below. Figure 3(a) illustrates the characteristics of a heterodyne power spectral densities (PSD) S het (ω), in regimes of bright dark-modes, where g |ω x − ω y |, κ, γ so dissipation is very low. Specifically, the parameters are set close to experiment of Ref. [24] but with θ tw = π/4, κ → κ/10 and |ω x − ω y | → |ω x − ω y |/4. The figure illustrates the typical structure of a three-level crossing with a coupling in the form of Eq. (4) (two degenerate modesx,ŷ couple indirectly via a third). It also illustrates, however, the difficulty of the usual procedure for thermometry in optomechanics. Both modes are hot, and moreover the usual normalization used to relate the heterodyne measured PSD to the underlying displacement spectra gives very poor results.

IV. TWO-DIMENSIONAL MOTION MODELING
In Ref. [28] the Langevin equations for the full 3D problem was solved, both for the full nonlinear trapping potentials as well as considering linearizations about equilibrium values. For the linearized case general expressions of quantum linear theory (QLT) for the optical modes were obtained as well as the mechanical spectrax 3D ,ŷ 3D , andẑ 3D , including hybridization and encompassing quantum regimes.
From these, PSDs for S xx (ω) = |x 3D (ω)| 2 or S yy (ω) = |ŷ 3D (ω)| 2 were calculated and thence phonon occupancies are related to the area under the PSD curve [33]: for j = x, y, z and compared with optical (heterodyne) PSDs in order to understand the experimental measurements. However, in Ref. [28], the analysis of quantum cooling for particles trapped at the node of the cavity, focused on quasi-1D experimental scenarios of Ref. [24] where only a single mode was strong-coupled to the light. Hybridization with weak-coupled modes, leading to sympathetic cooling or heating was investigated. In the present work we go beyond [28] to investigate the case of two strong-coupled modes, where nontrivial 2D physics arises, including the formation of dark modes.
In addition, we consider also scattering effects not included in Ref. [28]. The Hamiltonian of the system, Eq. (9) is a special case of the Hamiltonians discussed in Refs. [17,28] where the nanoparticle equilibrium position was primarily determined by the gradient force (∼70-nm particles). Here, however, we consider also substantially larger particles (∼100 nm) where the scattering force must be taken into account. In particular, the latter displaces the nanoparticle equilibrium position which leads to new couplings ofx andŷ to the phase quadrature of light,P (in addition, to the coupling to the amplitude quadrature,Ẑ). In Appendix A we show that the Hamiltonian can be transformed back to the form in Eq. (9) by introducing the rotated optical quadraturesẐ L (P L ) with the rotation angle depending on the size of the nanoparticle. In short, all of the results from Refs. [17,28] remain valid even when the scattering force is nonnegligible (but still in the Rayleigh regime) as long as we formally replaceẐ,P with the rotated optical quadraturesẐ L ,P L .
For our 2D analysis, we assume the particle is trapped at the cavity node, i.e., φ tw = π/2 that minimizes deleterious optical heating, and ask the following: What is the optimal angle θ tw between the tweezer polarization axis and the cavity symmetry axis for efficient 2D cooling? The latter controls the couplings g x , g y , and the most natural choice is given by θ tw ∼ π/4 where g x ∼ g y -this maximizes the cooperativities of both the x and y motions. In addition, we have the freedom in choosing the detuning, ( < 0 is red-detuned)-in the first instance this can be set to − = (ω x + ω y )/2. In particular, the perfectly degenerate case, ω x = ω y , where we have the exact relation g x = g y seems the most natural configuration for 2D cooling; however, we will show this is not the case, and nondegenerate frequencies are necessary for efficient 2D cooling.
The 2D equations of motionẊ = AX in Eq. (10), explicitly, are given by where we have for simplicity of presentation omitted the nonconservative terms (damping terms and input noises). We have the optical quadratures, We can express the equations for the modes of the 2D problem in Fourier space: which can be solved in closed form as shown in Ref. [28]. The solutions will be labeled asx 2D (ω),ŷ 2D (ω), andẐ L (ω) [which in general depend on all three input noisesx in (ω),ỹ in (ω) and Z L,in (ω)]. We recall that our calculations fully include the z mechanical motion. However, the 2D analysis below is closely matched by the full 3D theory and thus we can make the identifications,x 2D (ω) ≡x 3D (ω) andŷ 2D (ω) ≡ŷ 2D (ω), for the hybridized modes of the system. The frequency-dependent coupling coefficients are given by where j = x, y, the susceptibilities are given by and κ (γ ) is the cavity decay rate (gas damping). Equations (15)-(17) form a system of coupled equations: The solutionsx(ω),ŷ(ω),â(ω) are function of the input noises x in (ω),ỹ in (ω),ã in (ω). In addition to gas collisions and photon shot-noise we also include recoil heating in the model by adding additional terms tox in (ω) andỹ in (ω) [17,34,35]. The latter can become relevant even at pressure p ∼ 10 −6 mbar when we scan over large values of the couplings g x , g y 100 kHz as we increase the size of the nanoparticle/laser power.
In the following we consider the conditionn (2D) < 1 as a threshold value for 2D motional ground-state cooling. Alternatively we could have required the less restrictive conditionsn x ,n y < 1, i.e., the two motions are separately in the ground state.

A. Optimal frequency difference
The first question we address is what is the optimal frequency difference, δω = |ω x − ω y |, in order to achieve the lowest combined phonon occupancyn (2D) . For concreteness we will consider the parameters from Ref. [24] (but now with θ = π/4) and vary the two waists of the tweezer beam, w x and w y , to scan over the frequencies, ω x , ω y . We find that when |ω x − ω y | > κ simultaneous cooling of thex andŷ modes becomes ineffective-we either cool x motion or y motion but cannot cool both effectively. More surprisingly, we find that when |ω x − ω y | ∼ 0 cooling becomes again ineffective. The optimal frequency difference for efficient 2D cooling is near the midpoint value-when |ω x − ω y | ∼ κ 2 with the detuning set to − ∼ (ω x + ω y )/2 (see Fig. 4).
We can understand qualitatively the reason for the optimal frequency difference |ω x − ω y | ∼ κ 2 by calculating 2D optomechanical cooling formula: opt, j ≡ Im where j = x, k = y or j = y, k = x (from Eqs. (15)-(18) one readily finds the optomechanical cooling formula by calculating the imaginary part of the self-energy [28,36]). Let us consider some special cases. Suppose first that g k ∼ 0 such that we have opt, j ∼ Im[2ig 2 j η(ω j )]-the latter is the usual optomechanical cooling rate which further reduces to opt, j ∼ 4g 2 j /κ. The numerator can be thus associated with the cooling rate from standard 1D cavity optomechanics. On the other hand the denominator depends only on the coupling to the other degree of freedom, ∼g k , and is thus a genuinely 2D effect affecting the j motion.
We are primarily interested in the configuration where both opt,x and opt,y are large. Let us start by considering the perfectly degenerate case, ω x = ω y = − with g = g x ∼ g y . Assuming the regime of strong cooperativity we find that the optomechanical rate reduces to the simple expression opt, j ∼ γ . The gas damping, γ , is tiny at the relevant pressures, and thus we are left only with a negligible optomechanical cooling rate-here γ arises from the denominator in Eq. (23), i.e., from the mechanical susceptibility χ k (ω j ) defined in Eq. (20), and thus the strong suppression of the optomechanical cooling rate can be identified as a 2D effect. Loosely speaking, the energy that is extracted from the x motion (y motion) is immediately fed back to the y motion (x motion) with the optical field mediating this transition. In order to achieve any 2D cooling we thus require some degree of asymmetry, ω x = ω y , in order to disrupt the near-perfect exchange of energy between x andŷ via the optical field, and allow the latter to instead carry the energy away from the system.
We finally note that lowering the finesse does not necessarily improve 2D cooling. This is captured by the optomechanical cooling formula in Eq. (23) through the optical susceptibility η defined in Eq. (21): On the one hand, when we decrease the value of κ we enhance the 1D cooling channel (numerator), but, on the other hand, we also amplify the 2D heating channel (denominator).

B. Optimal particle size
For concreteness we will consider the parameters from Ref. [24] (but now with θ = π/4) which is close to the op- For a given experimental implementation the tweezer power, P tw , is a parameter that can be varied readily. However to minimize deleterious optical heating effects, it is preferable to use the lowest power that meets experimental requirements. We restrict ourselves to P tw 1 W. Nanoparticles of different radii R may also be selected.
We express the relevant parameters for 2D cooling, ω x , ω y , g x , g y as a function of P tw and R. We then find that if the particle size is too small ( 60 nm) the cooperativity remains low and one is limited to values above n x + n y ∼ 1-this is analogous to the requirement for 1D ground-state cooling. However, if the particle size is too large ( 100 nm), then cooling becomes again ineffective when g x , g y |ω x − ω y |. We find that there is a "Goldilocks zone" with the optimal particle size ∼80 nm (see Fig. 5).
We can understand qualitatively the reason for the optimal particle size by looking again at the 2D optomechanical Numerical simulation of 2D phonon occupancy as a function of input power (P tw ) and mean optomechanical coupling [g ≡ (g x + g y )/2]. The values are set as in Ref. [24] but θ = π/4 and variable power and particle radius (which sets the optomechanical couplings). Cooling becomes ineffective if the cooperativity is too low (below the lower green dashed line) as well as if decoupled dark modes are formed (above the upper green dashed line). The red dashed lines plot g for different values of the particle radius, R. We can thus extract the optimal particle size that allows efficient 2D cooling at moderate powers-cooling to the 2D motional ground state, n x + n y < 1, is feasible already for a ∼80-nm particle at P tw ∼ 700 mW.
cooling formula in Eq. (23). We first note that g j ∝ R 3/2 and that g j ∝ P 1/4 tw hence g j ∝ R 3/2 P 1/4 tw . For small (large) values of g j the numerator (denominator) in Eq. (23) is small (large) and cooling becomes inhibited-this illustrates how the "Goldilocks zone" emerges from the competition of the 1D effect in the numerator with the 2D effect in the denominator. In particular, the condition g j |ω x − ω y | emerges from the denominator in Eq. (23)-we have χ k (ω j ) ∼ |ω x − ω y | −1 as well as η(ω j ) ∼ |ω x − ω y | −1 [since − ∼ (ω x + ω y )/2]and hence the the denominator remains suppressed if |g k | |ω x − ω y |, i.e., cooling is not inhibited by the 2D hybridization effect.

C. Reliable thermometry
In the previous sections we have shown that there exists an optimal experimental configuration (Sec. V A) and particle size (Sec. V B) to achieve simultaneous cooling of both x any y motions. However, inferring phonon occupancies from optically detected spectra in the presence of hybridization is not straightforward. Here we show that the same experimental configuration that allows for optimal 2D cooling also allows for reliable thermometry.
Experiments exploiting heterodyne detection have access only to the optical field,â = 1 2 (Ẑ L + iP L ), from which one then extracts the the mechanical displacement spectra. In particular, the heterodyne PSD is given by [33] [24], δω = |ω x − ω y | > g x ∼ g y (upper) with occupancies for a near-circular trap with the same parameters but δω smaller so that |ω x − ω y | < g x ∼ g y (lower). The elliptical trap allows for 2D ground-state cooling and the rescaled heterodyne follows occupancies closely, facilitating thermometry. For the near-circular trap, the modes remain hot and it is difficult to extract occupancies from the optical detection by the usual methods. The particle is positioned at a node (intensity minimum), θ tw = π/4, R = 80 nm, input power P in = 0.8 W, and κ = 193 kHz. (b) Corresponding heterodyne PSDs, with the classical modes overlaid. For the elliptical trap in the upper panels when the detuning is set to − = 400 kHz the modes are cooled to n x + n y < 1.
where LO is the detuning of the local oscillator andâ out = a in − √ κâ is the output field.
In the presence of hybridization and spectral overlaps, extracting displacement spectra S xx (ω) and S yy (ω), from the experimental heterodyne spectra, S het (ω), becomes less straightforward. However, in the optimal case for 2D cooling (see Sec. V A)-when g x g y g, we can write: The corresponding heterodyne PSD [from Eq. (24)] will have contributions not only from independent x, y contributions, but also from interferences. Thermometry is greatly simplified when we neglect interference effects and are able to write: which can be seen as the 2D extension of the familiar textbook relation arising in the 1D case. Clearly more elliptical traps, with less spectral overlap between the x, y terms would be expected to minimize interference contributions. We test the approximation in Eq. (26) in Fig. 6: In Fig. 6(a) we compare the extracted phonon occupancies (using the heterodyne spectra) with the actual ones, and in Fig. 6(b) we show the PSDs for the heterodyne signal with the classical modes overlaid. The upper panels show the result for an elliptical trap lie within the "Goldilocks zone" where frequencies are sufficiently far apart so that both modes are in the quantum regime and moreover, the phonon occupancies inferred from rescaled heterodyne PSD agree reasonably well with those obtained from integrating S xx (ω) + S yy (ω), in contrast with the near-circular trap, which lies outside this zone.
In the latter case (near-circular trap) the modes remain hot and more complicated methods would be required to infer mode occupancies from the optically detected signal.

D. Understanding 2D cooling in terms of geometric bright/dark modes
In our effective 2D system,x 2D andŷ 2D are the mechanical modes. As shown in the previous section, they in general arise from hybridization ofx,ŷ modes and, as we consider regimes of strong-coupling, they also involve hybridization with the optical modeẐ L . Unhybridizedx,ŷ mechanical modes would correspond physically to motions along the tweezer x and y axes, respectively (such as in the case the coupling to the cavity mode would be vanishingly small and we would thus have two completely decoupled mechanical motions).
However, there is another pair of mechanical modes that naturally arises in the coherent scattering setup: These are the bright (dark) modes. We show below that, in our system, the former have an interesting and useful geometric interpretation in terms of modes corresponding to the motion along (orthogonal) to the cavity axis; together with the cavity mode,Ẑ L , they form the cavity-based/geometric modes.
The transformation from the tweezer-based modes to the geometric/cavity-based modes might, in first instance, be understood as a pure 2D rotation in the x-y plane by applying a rotation of angle θ tw . However, care is required when we consider traps with significant ellipticity, where ω x = ω y , where the zero-point motions distort the 2D rotation as we will now show. We start by considering the rotated reference frame where the first (second) axis of the reference frame is parallel (orthogonal) to the cavity axis (see Fig. 1). Specifically, we transform to such "cavity reference frame" by applying a rotation of angle θ tw : where X, Y [X (c) , Y (c) ] are the coordinates in tweezer (cavity) reference frame. The rotation in Eq. (27) in term induces a transformation of the canonical (adimensional) modesx and y. In a nutshell, one first rotates the corresponding physical positions (X andŶ ) to obtain the transformed physical positions (X b andX d ) and thence defines the transformed canonical positions (x b andx d ) by rescaling them with the transformed zero-point motions. Specifically, we perform the following transformations in consecutive order: where x j,zpf = h 2mω j (ω j ) is the zero-point motion (frequency) along the j = b, d axis (for the full details of the derivation see Appendix B 1 and B 2).
The Hamiltonian from Eq. (3) in terms of the new rotated coordinates reduces tô and thus only the modex b is coupled to the light field whilê x d is completely decoupled-we will refer to them as the geometric bright and dark mode, respectively. A similar coupling to the above was obtained in a previous experimental study of hybridization between two mechanical modes [29] leading to bright/dark modes (albeit not in a strong coupling regime and without strong cooling). However, what is new is that in the present case the resulting modes,x b andx d , have a simple geometric interpretation as the motion along and orthogonal to the cavity axis, respectively. Without the geometric interpretation as a guide one can consider also alternative definitions for bright/dark modes-these do not have a geometric interpretation but are otherwise equally valid (see Appendix B 4 for a comparison with the bright/dark mode considered in Ref. [29]).
The couplings in Eq. (31) are given by and the frequencies are given by Importantly, the above derivation of the geometric bright/dark modes is valid for any value of the angle θ tw (see Appendix B 3 for the special case θ tw = π/4 where the expressions simplify further). This is a specific feature of the coherent scattering setup where one can always decompose the motions in the transverse tweezer plane into the motion along/perpendicular to the cavity axis-the corresponding geometric bright/dark modes are by construction coupled/decoupled from the cavity mode. We note that the classical bright/dark eigenmodes of the drift matrix presented in Sec. III correspond closely tox b ,x d in the limit ω x = ω y where g bd = 0 and there is no coupling between them. In such a caseD can be identified withx d whilê B ± corresponds to the hybridization ofx b andẐ L . In general, and in the "Goldilocks zone," g bd = 0, and thus the relation betweenB ± ,D andx b ,x d ,Ẑ L becomes increasingly distorted (as bothx b ,x d start to hybridize withẐ L ).
In Appendix B 5 we also solve for the spectra of the hybridized mechanical modes of the system,x 2D b andx 2D d , analogously to Eqs. (15) and (16), but now given in terms of motions alongx b ,x d (andẐ L ) rather than motions along the tweezer axesx,ŷ (andẐ L ). The exact closed form solutions yields completely equivalent heterodyne spectra, but as we show below decomposing the spectra in terms ofx b ,x d mechanical contributions can be less straightforward.
Although the "Goldilocks zone" is not in the strongly hybridized regime, it is still in an intermediate regime where hybridization nonetheless plays a critical role. It is instructive to reexamine the 2D cooling behavior, but now in terms of thê x 2D b andx 2D d modes. We see from Eq. (31) that the geometric dark mode,x d , is decoupled from the cavity mode-the optomechanical cooling mechanism must rely on hybridization due to the coupling g bd to the bright mode. In other words, we expect to sympathetically coolx 2D d only when it is significantly hybridized and contains contributions fromx b ,x d , andẐ L . We see from Eq. (32) that the coupling g bd depends on |ω x − ω y |: Only when |ω x − ω y | is large can we can expect to cool botĥ In contrast, when ω x = ω y thex 2D d mode cannot be cooled.
The emergence of the "Goldilocks zone" can, in fact also be be analyzed in terms of the 2D optomechanical cooling rates forx 2D b andx 2D d (see Appendix B 5 for the derivation): opt,d ≡ Im where χ j is the mechanical susceptibility defined in Eq. (20) with j = b, d. In particular, we consider the ideal case θ tw = π/4 where we find simple expressions opt,b ∼ 4g 2 b /κ and opt,d ∼ g 2 bd κ/g 2 b . In order to cool effectively in 2D both opt,b and opt,d have to be larger than a certain minimum threshold value, opt,b,d 2 ; these two conditions give rise to the Goldilocks zone. Let us suppose g ≡ g x ∼ g y -we find g b ∼ √ 2g and g bd ∼ (ω y − ω x )/2-which further reduces the 2D optomechanical formulas to opt,b ∼ 8g 2 κ and opt,d ∼ κ 8g 2 (ω y − ω x ) 2 . Combining the two constraints we find the condition for the Goldilocks zone: where can be loosely identified with the total motional heating rate. The motional heating rate has a constant contribution (due to gas collisions) and a power-dependent contribution (from recoil heating). In first instance we can neglect recoil heating for moderate powers at the considered pressure of p ∼ 10 −6 mbar [17], and estimate ∼ γ n B , where γ is the gas damping and n B is the mean thermal occupancy. Using the values in Table I we estimate for /(2π ) ∼ 15 kHz and find a qualitative agreement of Eq. (38) with the Goldilocks boundaries shown in Fig. 5 (the lower and upper dashed green lines). The phonon occupancy of the modesx b ,x d is, however, quantitatively different from the one of the modesx,ŷ-the two sets of modes have different frequencies (ω x , ω y versus ω b , ω d ) which makes a direct comparison of the number of phonons difficult. Even disregarding the modest numerical discrepancies from the frequency differences, estimating the dark mode phonon occupancy from the heterodyne spectra is a nontrivial task. To see this we note that the optical field is proportional to the bright mode (but not to the dark mode): One can then directly extract the bright mode displacement PSD from the normalized heterodyne PSD: where We can thus obtain an occupancy n b by integrating the area under the heterodyne PSD rescaled by the factor g 2 b |η(ω)| 2 [see Eq. (11) with j = x b , x d and define n j ≡ n x j ]. On the other hand, from the measured heterodyne spectra one is not able to measure S x d x d (ω) = |x 2D d (ω)| and thus one cannot directly obtain an estimate for the corresponding phonon occupancy n d . However, by comparing Eqs. (25) and (40) (and using g b ∼ √ 2g at θ tw = π/4) we find S x b x b (ω) ∼ [S xx (ω) + S yy (ω)]/2 and thus n b ∼ (n x + n y )/2. When both the bright/dark mode are cooled with the same optomechanical rate (i.e., opt,b ∼ opt,d ) one has n b ∼ n d , and one can indirectly infer that the area under S x d x d (ω) = |x 2D d (ω)| will give n d ∼ (n x + n y )/2. We thus find n x + n y ∼ n b + n d . This explains why both the tweezerbased and the cavity-based modes approximately agree about the total phonon occupancy and lead to roughly the same Goldilocks zone.

VI. DISCUSSION
A previous theoretical study that investigated hybridization due to optomechanical interaction via coherent scattering by levitated nanoparticles, left an important gap in understanding: that study [28] investigated the quasi-1D dynamics arising for trapping at a node, relevant to recent quantum ground-state cooling experiments. However, the quasi-1D behavior g x g y of one strong-coupled mode, hybridising with one weak coupled mode, obtained for θ tw → π/2, differs fundamentally from the 2D scenario of two strong-coupled modes that arises as θ tw → π/4, where g x g y .
Here we have shown that simultaneous quantum groundstate cooling of both strong coupled modes is no longer achievable simply by the usual 1D strategy of increasing the coupling strengths/cooperativities of the individual modes: One must also factor in the essentially 2D phenomenon of the formation of dark modes that decouple from the optical field, as well as other requirements.
We have investigated the levitated nanoparticle motion in the tweezer transverse (x-y) plane with an optical cavity. We showed that efficient cooling of thex andŷ motion to their quantum ground state must obey certain constraints relating the difference of mechanical frequencies |ω x − ω y | to the coupling rates g x , g y as well as the cavity decay rate κ. We found that cooling and standard thermometry are efficient for a sufficiently elliptical optical trap, while for a more spherical trap the cooling will be hindered by strong three-way mode hybridization with the cavity mode. We found also the optimal particle size that satisfies the conditions, thus allowing for 2D ground-state cooling in the current experimental setup.
The analysis of the 2D levitated nanoparticle problem also found that the dark/bright modes have a geometric interpretation in terms of rotations in the x-y plane. Importantly, the transformation is not a trivial rotation because of the nonequivalence of thex,ŷ phononic modes (we assume trap ellipticity) so the modification for ω x = ω y is discussed. In addition, we also considered in the calculations the effects of scattering force so as to be able to reliably simulate larger particles.
Free-fall experiments that propose recycling of particleswhere particles would be trapped again after a sufficiently long free-fall time-require the nanoparticle energy to be low in all three translational motions. The motion along the optical tweezer axis can be cooled to its ground state via feedback cooling [26,27], thus extending our 2D scheme to fully prepare nanoparticles for free-fall experiments. In addition, the uncoupled three degrees of freedom can be used as a (quantum) sensor of forces acting along a specific direction, such as terrestrial gravity fluctuations [37].
Note added. Recently, we became aware of a new interesting experimental study by Ranfagni et al. [38] which for θ tw 0.4π exhibits a scenario intermediate between the quasi-1D limit and the 2D regimes we investigate here.

ACKNOWLEDGMENTS
This work was presented at the ICTP workshop Frontiers of Nanomechanics (September 2020) and at the conference Quantum Nanophotonics (March 2021). We are extremely grateful to Vladan Vuletić for insightful discussions about bright/dark modes. This work was supported by the Engineering and Physical Sciences Research Council [EP/N031105/1 and EP/L015242/1]. M.T. acknowledges funding by the Leverhulme Trust (RPG-2020-197). U.D. acknowledges support from the European Research Council (ERC CoG QLev4G) and the research platform TURIS at the University of Vienna.

APPENDIX A: NOTES ON SCATTERING FORCE
In this Appendix we look at the modification of the optomechanical interaction due to the shifted equilibrium position of the nanoparticle along the z axis (with respect to the tweezer trap center). In particular, such a shifted equilibrium arises from the scattering radiation pressure force as we increase the particle size. We first find the new equilibrium position (Sec. A 1) and then calculate the new optomechan-ical couplings (Sec. A 2). We finally show that by using appropriately rotated optical quadratures the optomechanical formulas derived for the case of negligible scattering forceappropriate for small nanoparticles-extend also to the case with with a shifted equilibrium position (Sec. A 3).

z-Axis equilibrium position
The competition between the gradient force, F grad , and the scattering force, F scatt , modifies the nanoparticle's equilibrium position, z 0 , along the z axis (the tweezer symmetry axis). In particular, the gradient and scattering force are given by [39]: respectively, where R is the relative dielectic permittivity, c is the speed of light, k = 2π λ , λ is the wavelength, and R is the particle radius. The laser intensity along the tweezer axis is given by where P tw is the laser power at the center of the trap, w x , w y are the beam waists, and z R = π w x w y /λ is the Rayleigh range. We readily find the equilibrium position, z 0 , from the condition F grad + F scatt = 0. Assuming z 0 /z R 1 we find a simple result: but can otherwise numerically solve for the equilibrium position. Importantly, the larger the particle radius, R, the more the equilibrium position, z 0 , is displaced from the Gaussian beam focus [40].

Optomechanical couplings
We start from the coherent scattering interaction potential [17]: where ξ = kz + (z) and (z) = −arctan(z/z R ) is the Gouy phase. The equilibrium position of the nanoparticle with respect to the tweezer trap center will be denoted by (x 0 , y 0 , z 0 ), where we assume x 0 = y 0 = 0, while z 0 is given by Eq. (A4) (for z 0 /z R 1) or obtained by solving numerically for the equilibrium position.
For completeness we list also the z couplings: where g yz = E d k 2 P 0 cos(θ tw )sin(φ)y zpf z zpf .
The couplings g xZ , g yZ , g xy , g zP , and g xz , g yz have been previously obtained by neglecting the scattering force and are valid for small nanoparticles [28].

Rotated optical quadratures
It is instructive to compare the case with negligible scattering force (i.e., z 0 = 0 and ξ = 0) with the case of an arbitrary z axis displacement from the tweezer trap center (i.e., z 0 > 0 and thus ξ > 0). In particular, we introduce the rotated optical quadratures: where the angle of rotation is ξ ∼ kz 0 , andẐ ξ ,P ξ (Ẑ,P) denote the optical quadratures in the case with (without) the z 0 displacement. Let us first consider the mean values. From Eq. (A5), writing the corresponding classical equations of motion, we find that the mean value of the optical quadratures are given by where is the detuning, and κ the cavity decay rate. If we set z 0 = 0 (and hence ξ = 0), then we find the simplified expression for the amplitude and phase quadratures, which we denote by Z 0 and P 0 , respectively. From Eqs. (A20) and (A21) we readily see that Z ξ 0 , P ξ 0 and Z 0 , P 0 are related by the rotation introduced in Eq. (A19).
Using now the rotated quadratures,Ẑ ξ ,P ξ , and the rotated mean values, Z ξ 0 , P ξ 0 , the interaction potential in Eq. (A6) reduces to the expression: where We note that the potential in Eq. (A22) has the same form of the potential previously obtained for the case of small nanoparticles [28] (i.e., where one can neglect the displacement, z 0 , due to the scattering force): One formally replaceŝ Z →Ẑ ξ andP →P ξ . Thus all the formulas obtained in Ref. [28] remain valid also when we consider a significant nonzero displacement along the z axis (displacement from the tweezer trap center), provided we use the rotated mean values, Z ξ 0 and P ξ 0 , given in Eqs. (A20) and (A21), respectively. For the special case considered in the main text the interaction potential remains of the same form as in the case without any z axis displacement: where we have definedẐ L ≡Ẑ ξ , g x ≡ g xZ , and g y ≡ g yZ .

APPENDIX B: BRIGHT/DARK MODES AND GEOMETRIC ROTATIONS
In this section we introduce bright/dark modes and show that they have a simple geometric interpretation as the motion along/perpendicular to the cavity axis (in the transverse tweezer plane). We first recall the usual 2D matrix that rotates the reference frames -we show how it induces a rotation of the physical positions and momenta while it leads to a distorted transformation of the associated canonical position and momenta (Sec. B 1). We then proceed to define the geometric bright/dark modes -we obtain the coupling of the bright mode to the optical field as well as the induced coupling to the dark mode (Sec. B 2) and provide the simplified geometric bright/dark mode couplings at θ tw = π/4 (Sec. B 3). For completeness we also compare the geometric bright/dark modes with the alternative bright/dark modes used in Ref. [29] (Sec. B 4). Finally, we derive the 2D optomechanical cooling formulae for the bright/dark modes (Sec. B 5).

Geometric and canonical rotations
The geometric rotation introduced in Eq. (27) relates the coordinates of the "tweezer reference frame" with the coordinates of the "cavity reference frame" (see Fig. 1). Such a rotation induces a transformation of the position and momenta-we will label physical positions (momenta) with capital letters (X ,Ŷ (P X ,P Y ) for the tweezer frame and X (c) ,Ŷ (c) [P (c) X ,P (c) Y ] for the cavity frame). Specifically, the physical positions transform as and the physical momenta transform as We recall that the physical positions/momenta are related to the canonical positions/momenta (denoted by lowercase letters) through a simple multiplicative rescaling: where the zero-point motions are given by where ω x , ω y , and [ω (c) x , ω (c) y ] are the oscillation frequencies along the axis of the tweezer frame (cavity frame). Hence the rotations in Eqs. (B1) and (B2) can be recast in terms of the canonical positions and momenta: and respectively. We thus see that a geometric rotation of the coordinates introduces a distorted transformation for the canonical variables. Only when the optical trap is perfectly degenerate [i.e., ω x = ω y which implies ω x = ω y = ω (c) x = ω (c) y ] we find that the transformations of the canonical variables in Eqs. (B11) and (B12) reduce to geometric rotations.

Definition of bright/dark modes
We start from the Hamiltonian in Eq. (9) which we recast in terms of the physical positions/momenta (see previous Sec. B 1): We first note that the terms on the first line of Eq. (B13) remain of the same form under the action of a geometric rotation and will thus be omitted in the following analysis [using Eq. (B2) one can readily show thatP 2 We now make the key observation-the coupling between Y (c) andẐ L in the last line of Eq. (B14) vanishes [we recall from Eq. (1) that g x ∼ sinθ tw /x zpf and g y ∼ cosθ tw /y zpf ] while the coupling betweenŶ (c) andẐ L reduces to −hE d k [the driving amplitude E d and the wave-vector k are defined below Eq. (1)]. It is thus appropriate to identifyŶ (c) as the dark mode (motion orthogonal to the cavity axis), andX (c) as the bright mode (motion parallel to the cavity axis), namely the geometric bright/dark modes. In summary, we have found that the geometric rotation of the reference frame (with the angle of rotation matching the angle between the tweezer polarization and cavity symmetry axis) leads to the bright/dark mode. We thus relabel the mechanical modes in the cavity reference frame as: From the first two lines of Eq. (B14) we can read the transformed mechanical frequencies: Furthermore, we can then rewrite the interaction terms of Eq. (B14) (third and fourth lines) as: where g bd = sinθ tw cosθ tw (ω 2

Dark/bright mode couplings at θ tw = π/4
The analysis in Appendix B 1 and B 2 is valid for any angle θ tw . We now write the couplings for the special case θ tw = π/4 considered in the main text where one has sinθ tw = cosθ tw = 1/ √ 2. From Eqs. (B17) and (B18) we first note that ω b = ω d and define The couplings in Eqs. (B20) and (B21) simplify to respectively.

Nongeometric bright/dark mode
It is instructive to compare the geometric bright/dark modes introduced above with the alternative bright/dark modes introduced in Ref. [29] which are obtained by rotations of the abstract space spanned by the canonical positions and momenta-we will refer to the latter as the nongeometric bright/dark modes. Specifically, in place of Eqs. (B11) and (B12) we now consider a simple rotation of the canonical positions, and of the canonical momenta, where θ ng is the angle of the abstract rotation. We now express Eq. (9) in the transformed coordinates to find:Ĥ (ω x sin 2 θ ng + ω y cos 2 θ ng ) x (c)2 +p (c)2 x + 1 4 (ω x sin 2 θ ng + ω y cos 2 θ ng ) ŷ (c)2 +p (c)2 In other words the nongeometric bright/dark modes cannot be interpreted as the motion along/perpendicular to the cavity axis.
Using trigonometric identities we also find sinθ ng = g x g 2 x + g 2 y , cosθ ng = g y g 2 x + g 2 y . (B28) We can thus recast Eq. (B26) in the form where the frequencies are given by and the couplings reduce to g bd = (ω y − ω x )g x g y 2 g 2 x + g 2 y . (B33)

Comparison with the geometric case
The note that the geometric [Eqs. (B11) and (B12)] and nongeometric transformation [Eqs. (B24) and (B25)] approximately match if the frequency differences are not too large. Let us see how to explicitly recover Eqs. (B32) and (B33) in the limit ω x ∼ ω y from Eqs. (B11) and (B12). We start from the last line of Eq. (B14) and using the fact that it vanishes we can express the angle as: Using trigonometric formulas we also find: cos(θ tw ) = √ ω y g y ω x g 2 x + ω y g 2 y . (B36) We now insert Eqs. (B35) and (B36) in Eqs. (B21) and (B20) [where we use the explicit expressions for the frequencies in Eqs. (B17) and (B18)]. We now finally write ω x = ω y + δω and Taylor expand the couplings in δω: The couplings arising from the nongeometric transformation can be thus seen as a limiting case of the couplings induced by the geometric transformation [compare with in Eqs. (B20) and (B21)]. There are however important differences. First, we note that the expression of g bd in Eq. (B33) is larger by a factor 2 with respect to the coupling obtained in Eq. (B37). Loosely speaking, this difference arises as the coupling g bd in Eq. (B29) couples bothx b ,x d andp b ,p d (i.e., two terms using the nongeometric construction) while in Eq. (B19) it couples onlyx b ,x d (a single term in the geometric construction). Second, we note that g d in Eq. (B32) contains only the lowest-order term of the expansion in Eq. (B38) which contains also a contribution ∼δω.

Two-dimensional cooling formulas for bright/dark modes
In this section we derive the 2D optomechanical cooling rates for the dark/bright mode. We start from the Hamiltonian in Eq. (31) and write Hamilton's equations of motion: where we have introduced g db = g bd to ease the reading of the equations. In the following we will consider also nonconservative terms (damping and input noise) which we have previously omitted for clarity of presentation. We transform Eqs. (B41)-(B44) to second-order differential equations by eliminating the momenta, and express the resulting equations in Fourier space: which can be readily solved. The solutions will be labeled aŝ x 2D b (ω),x 2D d (ω), andẐ L (ω) [each solution in general depends on all three input noisesx b,in (ω),x d,in (ω), andZ L,in (ω)]. The frequency-dependent coupling coefficients are given by where the susceptibilities are given by (B53) To find the self-energy for the bright mode,x b , we need to solve (B46) and (B47) forx d ≡x d (x b ) andẐ L ≡Ẑ L (x b ) and insert the expression in Eq. (B45) for the bright mode. To find the self-energy for the dark mode,x d , we proceed in a completely analogous way-we solve (B45) and (B47) forx b ≡x b (x d ) andẐ L ≡Ẑ L (x d ) and insert the expression in Eq. (B46) for the bright mode. From the imaginary parts of the self-energies we can then readily extract the optomechanical cooling rates: . (B55) We now consider the ideal case θ = π/4 where ω bd ≡ ω b = ω d (see Appendix B 3) and set the detuning to = −ω bd . After some algebra we eventually find opt,b ≡ Im 2ig 2 b η(ω bd ) + 4g 2 bd χ d (ω bd ) ≈ opt,d ≡ Im The corresponding phonon occupancies can be roughly estimated as n j ∼ / opt, j , where is the motional heating rate and j = b, d (as such a procedure gives only a crude estimate for the phonon occupancy we assume the motional heating rate is approximately equal for the bright/dark mode).
To obtain a quantitative estimate for the phonon occupancies one has to consider the PSDs: S x b x b (ω) = |x 2D b (ω)| 2 and S x d x d (ω) = |x 2D d (ω)| 2 . These are related to the corresponding phonon occupancies as the area under the PSD curve [33]: where j = b, d. We finally note that n b can be reliably extracted from the optical (heterodyne) PSD, but extracting n d is nontrivial and requires an indirect inference method (see the last paragraph of Sec. V D).