Stability and dynamics of optically levitated dielectric disks in a Gaussian standing wave beyond the harmonic approximation

Forces and torques exerted on dielectric disks trapped in a Gaussian standing wave are analyzed theoretically for disks of radius $2~\mu\text{m}$ with index of refraction $n=1.45$ and $n=2.0$ as well as disks of radius 200 nm with $n=1.45$. Calculations of the forces and torques were conducted both analytically and numerically using a discrete-dipole approximation method. Besides harmonic terms, third order ro-translational coupling terms in the potential energy can be significant and a necessary consideration when describing the dynamics of disks outside of the Rayleigh limit. The coupling terms are a result of the finite extension of the disk coupling to both the Gaussian and standing wave geometry of the beam. The resulting dynamics of the degrees of freedom most affected by the coupling terms exhibit several sidebands as evidenced in the power spectral densities. Simulations show that for Gaussian beam waists of $1-4~\mu\text{m}$ the disk remains stably trapped.

Forces and torques exerted on dielectric disks trapped in a Gaussian standing wave are analyzed theoretically for disks of radius 2 µm with index of refraction n = 1.45 and n = 2.0 as well as disks of radius 200 nm with n = 1.45. Calculations of the forces and torques were conducted both analytically and numerically using a discrete-dipole approximation method. Besides harmonic terms, third order ro-translational coupling terms in the potential energy can be significant and a necessary consideration when describing the dynamics of disks outside of the Rayleigh limit. The coupling terms are a result of the finite extension of the disk coupling to both the Gaussian and standing wave geometry of the beam. The resulting dynamics of the degrees of freedom most affected by the coupling terms exhibit several sidebands as evidenced in the power spectral densities. Simulations show that for Gaussian beam waists of 1 − 4 µm the disk remains stably trapped.

I. INTRODUCTION
The choice of particle used in levitated optomechanics is an important factor that depends on the goal of application. The most widely used particle in the field is a silica sphere with radius small compared to the wavelength. The dynamics of spheres trapped in cavities and focused laser beams are well understood and used for cooling to the motional ground state as well as force sensing [1][2][3][4]. This is owing to the simple harmonic translational and free rotational dynamics making it an ideal system to handle for both experimentalists and theorists. Particles with decreased particle symmetry allow rotational degrees of freedom to enter into the potential energy. A nanorod has large differences in moments of inertia and polarizability which allows rotations to be described as decoupled librations about the laser polarization axis. The motion of nanodumbells or generally anisotropic materials requires rigid-body dynamics since these particles have moments of inertia of similar magnitude [5,6]. Increasing the size of the particle relative to the wavelength of the laser further complicates the motion for any particle shape [7]. Still, terms necessary to describe nanorods and nanodumbbells have been investigated and the motion is also well understood [8][9][10][11].
Dielectric disks also have a relatively simple shape, but have not seen as much attention as other particle geometries. Several studies point to thin nanodisk scattering being more realistically described in a Rayleigh-Gans rather than a Rayleigh approximation for index of refraction n ∼ 1 [12][13][14]. This generally leads to an orientational dependent shape function in the form of a Bessel function. From studies investigating the applications of disks for various purposes, it is unclear whether there is consensus on the necessity of including the shape function or other non-harmonic terms in the dynamics [15][16][17]. There are few experimental studies involving disks, however two such studies suggest terms of higher order may be necessary for describing the motion [18,19].
In this paper it is shown that higher order terms of at least third order in the potential energy are necessary for describing the dynamics of disks outside the Rayleigh regime in a Gaussian standing wave. A disk experiences restoring forces in all three translational degrees of freedom and torques in two rotational degrees of freedom. Similar to rods and nanodumbells, the rotation about the disk's symmetry axis is unaffected by light coupling and is a constant of the motion. Focus is given to the effects due to the third order terms which provide unique rotranslational couplings that have not yet been discussed in levitated optomechanics. The coupling terms are a result of the finite extension of the disk coupling to both the Gaussian and standing wave geometry of the beam. Inclusion of the coupling terms results in dynamics with several different modes of oscillation for each degree of freedom which are evident in the power spectral density. Simulations show no evidence of instability.
An analytical as well as numerical approach using a discrete-dipole approximation method is used to identify the forces and torques on disks of radius 2 µm with index of refraction n = 1.45 and n = 2.0 as well as disks of radius 200 nm with n = 1.45. The Gaussian standing wave is constructed with a wavelength λ = 850 nm and various waists w 0 = 2, 2.5, 3, 4 µm.
The coupling terms presented in this paper may hinder or benefit applications for levitated disks. Disks have been proposed as potential accelerometers for gravitational wave detection [16]. The third order coupling terms may complicate determining which degree of freedom experienced a force or torque. On the other hand, it may be used as a means for indirectly detecting the motion of several degrees of freedom with a single detection scheme and therefore an efficient force/torque detector. Another common application is cooling the motion of the disk in attempt to study macroscopic quantum mechanics [20][21][22]. As energy from one degree of freedom can be transferred to another through the couplings, it may have potential for sympathetically cooling several degrees of freedom by performing a cooling method on only one of the degrees of freedom. Preliminary results show that this is indeed possible for both radii studied using parametric feedback cooling or cold damping.
The coupling terms are found to scale as the square ratio of the radius to the beam waist, a 2 /w 2 0 , and may therefore have less of an impact on the dynamics for particles of smaller radii compared to the wavelength. It is also found that the influence of the coupling may be reduced by sufficiently separating each degree of freedom's harmonic frequency.
This paper is organized as follows. Section II illustrates the analytical calculation of the potential energy of thin dielectric disks in a Gaussian standing wave. The potential energy is approximated to reveal a term third order in displacements and rotations. In Sec. III, the procedure for calculating forces and torques on a disk using the discrete-dipole approximation is outlined. The corresponding coefficients/frequencies are presented for various Gaussian beam waists. Lastly, Sec. IV examines the resulting dynamics due to the harmonic and coupling terms described in the previous sections.
FIG. 1. Coordinate system in the lab frame (x, y, z) and the particle frame (x , y , z ). The disk's symmetry axis is aligned with the particle frame z axis. The disk's center of mass as measured from the lab frame r0 is shown in red. The location to a point on the thin disk is given by the polar coordinates (ρ , φ ) in the particle frame which are shown in purple.

II. APPROXIMATE ANALYTICAL POTENTIAL ENERGY
This section outlines the analytical calculation of the potential energy of a thin dielectric disk in a Gaussian standing wave in the Rayleigh-Gans approximation. In this approach the disk thickness is taken to be very thin so that the Rayleigh approximation holds along that direction [14,15]. The approximated results verify the existence and helps elucidate the origin of the terms responsible for the dynamics seen in the following sections.
The disk is described with radius a, thickness T λ, index of refraction n, and susceptibilities χ = n 2 −1 and χ ⊥ = χ /n 2 corresponding to the susceptibility parallel and perpendicular to the disk symmetry axis (z axis), respectively [15]. The principal moments of inertia are I z = ma 2 /2 and I x = I y = m(3a 2 + T 2 )/12. The disk's center of mass is located at r 0 = x 0 , y 0 , z 0 and rotations are described in terms of the Euler angles (α, β, γ) in the z − y − z convention [5,23,24].
The Gaussian standing wave is formed by two counterpropagating Gaussian waves with non-zero longitudinal components so that they satisfy Maxwell's equations [25]. Each traveling wave has the symmetric waist w 0 , wavenumber k = (±x)2π/λ, and is polarized in theẑ direction. Around the focus, x = 0, each wave takes the form where x R = kw 2 0 /2 is the Rayleigh range and + (−) stands for the the right (left) traveling wave. The incident fields used for the numerical calculations in Secs. III and IV are found by propagating Eq. (1) throughout all space using the angular spectrum representation [25]. For the analytical calculations performed in this section and in Appendix B the approximated Gaussian standing wave is used, which is valid in the space |x| x R . The mechanical potential energy associated with the interaction between the light and the dielectric is where the integral is over the volume of the disk, χ 0 is the diagonal susceptibility matrix in the nanoparticle frame, and ↔ R is the rotation matrix. The rotation matrix in terms of the Euler angles explicitly can be found in Appendix A. The potential energy in the Rayleigh-Gans approximation with the incident field, Eq. (2), becomes where χ 1 = ∆χ cos 2 β + χ ⊥ , χ 2 = ∆χ sin β cos β sin α, ∆χ = χ − χ ⊥ , and a higher order term proportional to z 2 (r ) sin 2 kx(r ) was dropped. To evaluate Eq. (4) the coordinates of the disk must be projected onto each lab frame coordinate (x, y, z) and it is favorable to move to polar coordinates. First, in the limit T λ the functions in Eq. (4) are independent of the thickness leaving the functions in the integral dependent only on the disk's radial and angular coordinates (r ) = (ρ , φ ) (see Fig.  1). In terms of the center of mass and disk coordinates, x(r ), y(r ), and z(r ) in Eq. (4) are and the R ij are matrix components in the rotation matrix ↔ R (see Appendix A). Insertion of Eqs. (5) and (6) into Eq. (4) leads to analytic solutions in terms of Bessel functions. In the limit of small radius w 0 a, r 0 a w 2 0 where the zeroth order approximation to the exponentials (∼ 1) can be used, a Bessel function of the first kind is obtained as was found in Ref. [15]. However, this approximation misses the coupling of the disk to the Gaussian standing wave and a fourth order expansion in the coordinates is required to resolve it. Practical parameters in levitated optomechanics are in the range (λ, w 0 ) ∼ 1 µm and (r 0 , a) ∼ 1 − 0.1 µm. For the derivation, the limits a 2 w 2 0 , r 2 0 w 2 0 are used to expand each function in Eq. (4) to fourth order in the coordinates and terms O(a 6 /w 6 0 ) as well as O(x n 0,i π m j ), where n + m ≥ 4, π j = (α, β), are dropped which retains terms up to third order in the coordinates. Due to the symmetry of the disk, the potential energy is independent of the angle γ. Further, the disk's symmetry axis is primarily aligned along the lab framex direction and, as will be justified in the next section, rotates at angles that justify the small angle approximation α → 0 + θ z , β → π/2+θ y with θ z , θ y , small. Here θ z represents small angle rotations about the lab frame z axis while θ y is a small rotation about the lab frame y axis. The resulting potential energy is of the form Explicit expressions for the ω i may be found in Appendix B. The terms in the first row in the above potential describe simple harmonic motion for the three translational and two rotational degrees of freedom. The last term is a coupling between the translational and rotational degrees of freedom that is of third order in the coordinates. The coupling terms arise due to the finite radius of the disk and the Gaussian and standing wave geometry of the beam. An asymmetric electric field gradient across the disk produces a stronger force on the section of the disk with greater laser intensity. That section of the disk is pulled into the region of the trap with greater laser intensity more strongly than the section of the disk with less field intensity. As the radius increases and the trap becomes more confining, the greater the electric field gradient across the disk and the more influential the coupling terms are. With reference to Eq. (4), it is a result of the ro-translational coupling in the Gaussian together with the x 0 dependence in cos 2 kx(r ) describing the standing wave. The asymmetry in the (ω 1 , ω 2 ) coefficients is due to thex component of the incident electric field propor- To garner an idea of the dynamics that arise due to the coupling, consider the x 0 y 0 θ z term in Eq. (7). A disk displaced by r 0 = x 0 , y 0 , 0 in Fig. 1 experiences a torque about the −z axis due to a greater electric field intensity on the side of the disk nearest the focus. These terms are therefore a gradient force/torque as a consequence of the electric field gradient along the finite extension of the disk.

III. NUMERICAL EVALUATION OF THE FORCES AND TORQUES
A. System and Procedure The optical scattering problem for finite sized dielectric objects is generally difficult to solve analytically. As was done in the previous section, approximations are often required to glean insight into the dynamics. Another rigorous approach is to numerically solve for the scattered electromagnetic waves and use the resulting Maxwell stress tensor to obtain the forces and torques. This section details the results from performing the latter method by numerically implementing the discrete-dipole approximation (DDA) to calculate the scattered fields of the disk [26,27].
In the DDA, the disk is composed of N discrete spherical dipoles each with polarizibility α and the internal fields of the dielectric are solved for self-consistently to retrieve the scattered fields outside the particle. In the implementation of the DDA used for this paper, each dipole that composed the spherical dipole had a polarizibility α = 4π 0 R 3 n 2 − 1 n 2 + 2 . The method developed has been shown to be accurate to within 1% by comparing the scattered fields from a discretized sphere to the exact Mie scattering solutions [28]. The scattered fields that are generated from the DDA are then added to the incident field and inserted into the Maxwell stress tensor [29] T in order to obtain the forces and torques where The surface over which the integration is performed was taken to be a sphere centered at the disk center with radius 1.5× that of the disk. The surface integration was performed using Gaussian quadrature with increasing number of points until convergence was demonstrated.
The above procedure was performed for dielectric disks located near the intensity maximum of a Gaussian standing wave. To construct the standing wave, a righttraveling wave, E R (x, y, z) is found by propagating Eq.
(1) throughout all space using the angular spectrum representation with no paraxial approximation. A lefttraveling wave, E L (x, y, z) = E R (−x, −y, z), is added to the right-traveling wave to form the standing wave. The wavelength of each wave is λ = 850 nm and is fixed throughout this paper. A range of Gaussian beam waists were explored w 0 = 2, 2.5, 3, 4 µm to define the optical trap.
Most of the calculations performed were for disks of radius a = 2 µm, thickness T = λ/4n to achieve maximum light coupling, and index of refraction n = 1.45 or n = 2.0. The indices of refraction correspond to materials composed of silica and silicon nitride, respectively. Unless otherwise stated, the data and discussions that follow will refer to this set of parameters.
The following example outlines the steps for how a calculation is performed: the disk's symmetry axis is aligned with the axial direction (x axis), the disk is displaced a distance y 0 from the focus of the standing wave, the scattered waves are calculated using the DDA, the forces and torques are computed using Eqs. (9) and (10). The process is identical for rotations: the disk is initially situated at r 0 = 0, 0, 0 and (α = 0, β = π/2), a rotation is made α = 0 + θ z , the scattered waves are calculated using the DDA, the forces and torques are calculated. The baseline for the calculations is when the disk is placed symmetrically at the focus of the standing wave, r 0 = 0, 0, 0 , (α = 0, β = π/2) which should be a potential minimum. Indeed, a force or torque due to a displacement generally gives a value at least ten orders of magnitude greater than the baseline.

B. Forces and Torques
As is expected in levitated optomechanics, small displacements in one direction reveals a spring force in that same direction F i = −k i x i,0 and torque τ i = −κ i π i , π i = (α, β). The spring constants for each degree of freedom, (k i , κ i ), are determined by direct division, k i = −F i /x i,0 . At the harmonic level, no coupling of the different degrees of freedom through the potential energy were found.
Being that there are 6 degrees of freedom (including γ), there are 15 different second order couplings possible in the forces and torques. Of these possibilities, only terms similar to that in Eq. (7) were found to be above the baseline. These terms were found to be significant for disks of large and small radii. For example, a displacement of the center of mass by r 0 = x 0 , 0, z 0 produces a torque about the y axis, suggesting a term in the potential energy U ∝ D 1 y 0 z 0 θ y , with D 1 a proportionality constant. A similar coupling of the same order was found U ∝ D 2 x 0 y 0 θ z , with D 2 = D 1 necessarily. The coefficients D 1 and D 2 are also determined by division, i.e. D 1 = F z /(y 0 θ x ). Interestingly, the coefficients computed in this way generally gives different values for the force in theŷ andẑ directions for the first coupling term, and for the second coupling term, with A ≈ C 1 and B ≈ C 2 .
The coefficients A and B can differ from C 1 and C 2 by 2% using a waist of w 0 = 2 µm and 20% using a waist of w 0 = 4 µm. Although the discrepancy is suspected to be due to higher order terms, we are only interested in the dynamics due to this term and the average values D 1 = (A + 2C 1 )/3 and D 2 = (B + 2C 2 )/3 will be used from here on so that potential energy can be written in the form of Eq. (7). The consequences of using the average values is insignificant and will be discussed in Sec. IV. The spring and coupling constants (k i , κ i , D i ) have the same units and are most useful when written in terms of frequencies for translational harmonic motion, for rotational harmonic motion, and for the coupling terms.
w0 (µm) ωx (kHz) ωy (kHz) ωz (kHz) ω θy (kHz) ω θz (kHz) ω1 (kHz) ω2 (kHz)  2  394  38  38  537  390  46  39  3  264  17  17  361  263  21  17   TABLE I. Frequencies for a silica disks of radius a = 200 nm and thickness T = λ/(40n) for two beam waists w0 = 2, 3 µm. The disk has dimensions that are reduced by a factor of ten from the a = 2 µm, T = λ/(4n) disks. A fixed total power of 100 mW is used for the calculations. The number of points used to compose the disk was N = 37488 and the thickness of the disk was 4 points. Values for the frequencies as a function of beam waist are shown in Fig. 2 for silica and Fig. 3 for silicon nitride using a fixed total laser power of 100 mW. The general trend identified from the figures is that each frequency decreases as the waist increases. This feature is not unexpected, however, for particles in the Rayleigh regime λ a, ω i ∝ 1/w 2 0 while for a = 2 µm disks the dependence is nearly linear.
For both materials, the frequency in the axial direction is in the 150 − 200 kHz range while the radial degrees of freedom oscillate in the 1 − 10 kHz range. The axial frequency is most strongly affected by the standing wave which is independent of the waist. However, the radial frequencies are dominantly due to the Gaussian geometry. To leading order (see Appendix B), for fixed power the axial frequencies depend inversely on the wavelength and waist ω x ∝ 1/λw 0 while the radial frequencies depend on the waist as ω y,z ∝ 1/w 2 0 , hence the disparity between the axial and radial frequencies. Note that part of the waist dependence on each frequency is due to the dependence of the laser intensity on the waist. Each frequency therefore shares a 1/w 0 dependence from the power.
For a 2 µm radius disk at T = 300 K, these frequencies correspond to translational oscillation amplitudes of x 0 ∼ 1 nm and (0, y 0 , z 0 ) ∼ 20 nm. The rotational frequencies are closer to the axial frequency and in the range 190−125 kHz. The rotational frequencies differ by 20% between the two materials at the same waist. Using the average frequency, this corresponds to angular displacements of ∼ 1 mrad. Displacements of this size justify some of the approximations made in Sec. II since r 0 w 0 and sin α ≈ θ z .
Also shown in Figs. 2 and 3 are the coupling coefficients (ω 1 , ω 2 ). The coefficients being in the 50−200 kHz range are comparable to both the rotational and axial frequencies. Due to the large coupling frequencies combined with the relatively large oscillation amplitude in the radial degrees of freedom, the resulting forces/torques due to the coupling terms have an impact on the dynamics as shown in Sec. IV.
Force and torque calculations were also performed for silica disks of radius a = 200 nm and thickness T = λ/(40n) for the two beam waists w 0 = 2, 3 µm. The dimensions are 10× smaller than the a = 2 µm, T = λ/(4n) disk. The resulting frequencies are shown in Table I. From the table, each frequency scales as ω i ∼ 1/w 0 except for the radial frequencies (ω y , ω z ) ∼ 1/w 2 0 . This dependence on the waist is consistent with the analytical frequencies given in Appendix B. Also from the table, each harmonic frequency is larger, and the coupling frequencies reduced, compared to its a = 2 µm and T = λ/(4n) counterpart in Fig. 2. The dependence of each frequency on the radius is also consistent with that found analytically in Appendix B. The harmonic frequencies increase as the radius decreases since the disk has greater field intensity per volume. The coupling frequencies scale as ∼ a/w 2 0 due to the electric field gradient across the disk. This dependence provides a factor of ten between the a = 200 nm and a = 2 µm coupling frequencies.
C. Accuracy of the DDA

FIG. 4.
Frequencies obtained for a a = 2 µm silica disk using DDA for varying number of points that the disk was composed of relative to the frequency obtained using 299744 points, ωi,0. The legend describes the various frequencies for the x, y, z, θy, θz degrees of freedom as well as the ω1, ω2 coupling frequencies. The data points along the x-axis are 4680, 15804, 37488, and 299744 points. Comparing the left and rightmost data points in the figure shows that using 64 times more points changes the frequencies by less than 2%.
The frequencies shown in Sec. III B were obtained through several numerical operations such as integrations and the implementation of the DDA. One of the major questions regarding convergence of these values is how many points (i.e. number of discrete dipoles), N , should be used to discretize the disk. Figure 4 shows the relative change of the various frequencies discussed in the previous subsections as a function of the number of points used to compose the disk. Here, ω i,0 is the frequency calculated using the largest number of points shown in the plot, N = 299744. The frequency calculated using N points is ω i . The change in the frequency ω i compared to ω i,0 points is then ∆ω i = ω i − ω i,0 . The plot is shown for all of the various frequencies discussed above using a a = 2 µm silica disk with a w 0 = 2 µm waist. Increasing the number of points by a factor of 64 from N = 4680 to N = 299744 changes the frequency by less than 2%. On the other hand, the time complexity of the DDA method used to calculate the scattered light from the disk scales as N ln N .

IV. DYNAMICS
The previous two sections have illustrated that disks levitated in Gaussian standing waves experience simple harmonic motion as well as non-harmonic forces and torques involving second order couplings. This section discusses the resulting dynamics due to these forces and torques as well as the natural torques that arise in rigid body dynamics.
Thus far the focus has been on identifying terms in the potential energy. For translational motion the kinetic energy is trivial and leads to the equations of motion for small angle oscillations. As was shown in Ref. [5], for a symmetric top-like rigid body the rotational kinetic energy naturally involves coupling between the α,α, β, andβ degrees of freedom. Whether these terms are significant or not depends on the geometry. For a = 200 nm disks each non-linear coupling term is significant and must be considered. For a = 2 µm disks, the term responsible for precession about the x axis is the largest, but is still 10 −4 times smaller than the harmonic term and is therefore negligible. The equations of motion for a = 2 µm disks are then written asθ for small angle oscillations.
FIG. 5. Example trajectories of the x and two rotational degrees of freedom as well as the power spectral density of the axial motion for a a = 2 µm silica disk in a w0 = 3 µm waist trap. The influence of the second order coupling term produces several amplitude modulations at different frequencies, but the disk remains stable. The frequencies of modulation in the x degree of freedom can be seen in the power spectral density. Note that the rotational amplitudes remain in the ∼ mrad range, justifying the small angle approximation. Figure 5 shows sample trajectories of the θ y , θ z , and x motions of a a = 2 µm silica disk in a w 0 = 3 µm waist Gaussian standing wave by simulating Eqs. (20) to (24) at T = 300 K. The influence of the second order coupling terms are seen to be significant for the three degrees of freedom with each trajectory containing modulations at various frequencies. Without the couplings the oscillations would be at the same amplitude for all times. In a gaseous environment these modulations might be mistaken for noise in an experiment.
The bottom-rightmost plot in Fig. 5 shows the power spectral density (PSD) of the x motion. The harmonic frequency ω x /2π = 163 kHz is the largest and rightmost peak in the PSD. The other frequencies in the figure are the harmonic frequency plus the sums and differences of the various y, z, θ y , and θ z frequencies. Whereas sidebands due to coupling typically appear symmetrically on each side of the harmonic frequency, the frequency structure seen in Fig. 5 is such that all significant modes have smaller frequency than the harmonic frequency. This is not a general feature of the coupling term and depends on the degree of freedom that is being observed and the various levels of degeneracy.
The influence of the coupling term on each degree of freedom has two factors: the size of the coefficients ω 1 and ω 2 , and the level of degeneracy of the coupled degrees of freedom. First, the ω 1 and ω 2 coupling coefficients are relatively large ∼ 100 kHz. Second, strong coupling is achieved when the frequencies are nearly degenerate. Because ω x , ω θy , and ω θz are close in frequency the coupling term produces a larger effect on these degrees of freedom. Since the radial degrees of freedom oscillate 10× slower, the influence of the coupling term is significantly reduced, but not absent.
The question of stability is one of the most important for applications using levitated nanodisks. Despite the seemingly complicated motion, simulations have shown no evidence that this motion is unstable. The disk remains stable in the trap after several thousand oscillations for all of the beam waists explored w 0 = 2, 2.5, 3, 4 µm. The a = 200 nm disk was found to be stable at all frequencies, even with inclusion of the nonlinear coupling terms in the rotational kinetic energy [5]. Recall from Sec. III B the differing coefficients in Eqs. (11) to (16) as produced from the DDA calculations. Simulating the equations of motion with different coefficients attached to each degree of freedom's coupling term causes no issue for stability.
A common application in levitated optomechanics is cooling the motion of the levitated particle in attempt to reach the ground state, or to reach lower pressures [4]. The couplings in this paper offer a possibility of cooling one or more degrees of freedom sympathetically by actively cooling only one degree of freedom. The full dynamics of cooling using the couplings is beyond the scope of this paper, but we note some preliminary findings. Through simulations of the equations of motion Eqs. (20) to (24), results show that sympathetic cooling is indeed possible. For both radii, parametric feedback or cold damping [30,31] is an effective method for cooling multiple degrees of freedom. By inserting artificial numbers for the frequencies in the simulation, two relations were found for optimal cooling. Frequencies tailored within a few kHz of the relations ω x = ω θy ± ω z and/or ω x = ω θz ± ω y , can achieve significant sympathetic cooling to at least the mK regime. From Fig. 2, a = 2 µm silica disks are naturally in this regime. From Appendix B, each frequency depends on several parameters and has the possibility to be tuned to achieve optimal cooling experimentally.

V. CONCLUSION
The forces and torques exerted on dielectric disks trapped in a Gaussian standing wave were analyzed for disks of radius 2 µm with index of refraction n = 1.45 and n = 2.0 as well as disks of radius 200 nm with n = 1.45. Calculations of the forces and torques were conducted both analytically and under numerical simulation using a discrete-dipole approximation method.
Similar to nanodumbbells, a nanodisk experiences restoring forces in all three translational degrees of freedom, restoring torques in two rotational degrees of freedom, and has constant spin about the symmetry axis. Due to the finite geometry of the disk, third order, rotranslational coupling terms in the potential energy are found to be a necessary consideration when describing the dynamics of disks. The coupling terms are the result of an electric field gradient across the disk and depend on the ratio of the radius to the beam waist and on the temperature.
The ro-translational coupling produces several modes of oscillation in the coupled degrees of freedom which are evident in the power spectral density. While the restoring forces are dominant, the coupling terms can become sizable through strong coupling, which manifests when the coupled degrees of freedom are nearly degenerate. Despite the couplings, simulations show no evidence that the motion is unstable, which is of utmost importance for applications such as gravitational wave detection, force sensing, and ground state cooling.

ACKNOWLEDGMENTS
Supported by the Laboratory Directed Research and Development program at Sandia National Laboratories, a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia LLC, a wholly owned subsidiary of Honeywell International Inc. for the U.S. Department of Energy's National Nuclear Security Administration under contract de-na0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do no necessarily represent the views of the U.S. Department of Energy or the United States Government.
We would like to acknowledge Alejandro Grine, Darwin Serkland, Justin Schultz, Michael Wood, Peter Schwindt, and Tongcang Li for motivation of pursuit on the topic and useful discussions.