Simulations of atomic trajectories near a dielectric surface

We present a semiclassical model of an atom moving in the evanescent field of a microtoroidal resonator. Atoms falling through whispering-gallery modes can achieve strong, coherent coupling with the cavity at distances of approximately 100 nanometers from the surface; in this regime, surface-induced Casmir-Polder level shifts become significant for atomic motion and detection. Atomic transit events detected in recent experiments are analyzed with our simulation, which is extended to consider atom trapping in the evanescent field of a microtoroid.


Introduction
Strong, coherent interactions between atoms and light are an attractive resource for storing, manipulating, and retrieving quantum information in a quantum network with atoms serving as nodes for quantum processing and storage and with photons acting as a long-distance carrier for communication of quantum information [1]. One realization of a quantum node is an optical cavity, where light-matter interactions are enhanced by confining optical fields to small mode volumes. In the canonical implementation, a Fabry-Perot resonator with intracavity trapped atoms enables a panoply of cavity quantum electrodynamics (cQED) phenomena using single photons and single atoms, and thereby, validates many aspects of a cQED quantum node [2,3].
Despite these achievements, high-quality Fabry-Perot mirror cavities typically require significant care to construct and complex experimental instrumentation to stabilize. These practical issues have begun to be addressed by atom chips [4,5], in which atoms are manipulated in integrated on-chip microcavity structures offering a scalable interface between light and matter [6,7,8]. Owing to their high quality factors, low mode volumes, and efficient coupling to tapered optical fibers [9], microtoroidal resonators are a promising example of microcavities well-suited for on-chip cQED with single atoms and single photons [10]. Strong coupling [11,12] and non-classical regulation of optical fields [13,14] have been demonstrated with atoms and the whispering-gallery modes of a silica microtoroidal resonator.
In our experiments with microtoroids, Cs atoms are released from an optical trap and fall near a silica toroid, undergoing coherent interactions with cavity modes as each atom individually transits through the evanescent field of the resonator. In the most recent work of [12], atom transits are triggered in real-time to enable measurement of the Rabi-split spectrum of a strongly-coupled cQED system. Whereas a single atom is sufficient to modify the cavity dynamics, falling atoms are coupled to the cavity for only a few microseconds. Atom dropping experiments necessarily involve a large ensemble of individual atomic trajectories and represent, consequently, a far more complex measurement result.
Interactions between a neutral atom and a dielectric surface modify the radiative environment of the atom resulting in an enhanced decay rate [15] and Casimir-Polder (CP) forces [16,17]. These perturbative radiative surface interactions are usually insignificant in cQED experiments with Fabry-Perot resonators where atoms are far from mirror surfaces, but in microcavity cQED, atoms are localized in evanescent fields with scale lengths λ/2π ∼ 150 nm near a dielectric surface. The experimental conditions for microtoroidal cQED with falling atoms in [12] necessarily involve significant CP forces and level shifts while simultaneously addressing strong coupling to optical cavity modes. Theoretical analysis of this experiment requires addressing both the strong atom-cavity interactions and atom interactions with the dielectric surface of the microtoroid. As reported in [12], spectral and temporal measurements offer signatures of both strong coupling to the cavity mode and the significant influence of surface interactions on atomic motion. The role of these effects is quantified with detailed simulation of the trajectories of falling atoms detected in the real-time at low photon numbers.
In this article, we discuss in detail the approach used to simulate atomic motion near the surface of an axisymmetric dielectric resonator under the influence of strong coherent interactions with cavity modes. The experimental detection method of [12] is implemented stochastically in a semiclassical simulation of atom trajectories. These simulations provide a perspective on the atomic motion of atom transits recorded in our microtoroid experiments, while offering additional insights into the loading of optical evanescent field traps. In section 2, we outline the semiclassical model of a two-level atom coupled to the whispering gallery modes of a microtoroidal resonator. In section 3, we review the optical dipole forces which are a critical factor influencing atomic motion in an optical cavity. Our calculations of modified emission rates and Casimir-Polder surface interactions are detailed in section 4. Section 5 describes the implementation of our model for simulating recent atom-toroid experiments. Finally, section 6 extends our simulation to evanescent field traps around a microtoroid.

Atoms in a microtoroidal cavity
We approach the motion of atoms moving under the influence of surface interactions and coherent cavity dynamics with a semiclassical method to efficiently simulate a large number of atom trajectories. For surface interactions, dispersion forces are calculated perturbatively using the linear response functions of SiO 2 and a multi-level atom. For nearly-resonant non-perturbative coherent interactions between atom and cavity, the atomic internal state and the cavity field are treated quantum mechanically within the two-level and rotating-wave approximations.
Simulations of atomic motion follow the semiclassical method detailed in [18]. Mechanical effects of light are incorporated classically as a force F ( r) on a point particle atom at location r. Trajectories r(t) are calculated with a Langevin equation approach to incorporate momentum diffusion from fluctuations. At each simulation time step t i , the atomic velocity is calculated as: where v i is the velocity at the i time step, m Cs is the atomic mass, and ∆t is the simulation time step t i+1 − t i . The W i are normally distributed with zero mean and standard deviation of 1. Given the force F and diffusion tensor D ij as discussed in section 3, the atom trajectory r(t) and cavity transmission and reflection coefficients, T (t) and R(t) are calculated. A single atom strongly coupled to the cavity mode has a large effect on cavity fields and optical forces, requiring simultaneous solutions of atomic motion and cQED dynamics. Full quantization of atomic motion leads to an unwieldy Hilbert space not conducive to efficient simulation. In contrast, semiclassical methods are well-suited for simulating atomic motion in experiments with falling atoms near resonators. The ratio of the recoil energy to the linewidth of the cesium 6S 1/2 → 6P 3/2 transition is less than 10 −3 . Further, the recoil velocity of ∼ 3.5 mm/s is much less than the typical velocity of falling atoms of order 200 mm/s so that each spontaneous emission event represents a small momentum kick. Cavity fields and internal atomic states respond quickly to environment changes, allowing calculation of optical forces and momentum diffusion in a constant-velocity limit at time t and energy shifts from surface interactions as if atom the atom were stationary. The remainder of this section discusses the quantum mechanical equations of motion for the atom and cavity fields in the low-probe intensity limit, to be followed later by contributions to the force F used in (1).

Modes of a microtoroidal resonator
An idealized microtoroid has axial symmetry, so we work in a standard cylindrical coordinate system r → (ρ, φ, z). The toroid is modeled as a circle of diameter D m with dielectric constant revolved around the z-axis to make a torus of major diameter D M (Fig. 1(a)). The toroid is therefore defined by its minor diameter D m and its principal diameter D p = D M + D m . The fabrication and characterization of high-quality microtoroids are described in detail elsewhere [9].
The axisymmetric cavity modes of interest are whispering-gallery modes (WGM), which lie near the edge of the resonator surface and circulate in either a clockwise or counter-clockwise direction. These modes are characterized by an azimuthal mode number m, whose magnitude gives the periodicity around the toroid and whose sign indicates the direction of propagation. The WGMs for ±m are degenerate in frequency but travel around the toroid in opposite directions. The mode electric fields for the WGM traveling waves are written as E( r) = E max f (ρ, z)e imφ , where f (ρ, z) = E(ρ, z)/E max is the mode function in the ρ − z cross-section normalized by the maximum electric field E max . In general, backscattering couples these two modes so that a more useful eigenbasis for the system consists of the normal, standing wave modes characterized by a phase and the periodicity |m|. This backscattering coupling h is assumed to be real, with the phase absorbed into the definition of the origin of the coordinate φ. In addition, the mode's field decays at a rate κ i through radiation, scattering, and absorption. In our simulations, a cavity mode is fully described by its spatial mode function f ( r), its azimuthal mode number m, its loss rate κ i , and the coupling h to the counter-propagating mode with mode number −m.
We model the microtoroid modes using a commercial finite-element software package (COMSOL) to solve numerically for the vector mode functions f (ρ, z) for modes of a given m [19]. Mode volumes are calculated from, In this notation [10], the coupling of a circularly polarized optical field to an atomic dipole located in the evanescent field of the cavity is calculated as: where d is the dipole operator and ω (0) a = 2πc/λ 0 is the vacuum transition frequency of the two-level atom with free-space wavelength λ 0 . WGMs are predominantly linearly polarized, and so we average over the dipole matrix elements to obtain an effective traveling wave coupling g tw which is approximately ∼ 0.6 of the value for circularly polarized light (see supplementary information of [11]. Travelling wave modes of an axisymmetric resonator are not strictly transverse. For the toroid geometries considered here, with D p , D m λ, the azimuthal component is small and we assume that the optical field is linear outside of the toroid. Since the cavity losses are dominated by absorption and defect scattering rather than the radiative lifetime set by the toroid geometry [10], we let κ i and h be experimental parameters. Fig. 1 shows the lowestorder mode with m = 118 for a toroid with {D p , D m } = {24, 3} µm. The index m is chosen so that the cavity frequency ω c is near the 6S 1/2 → 6P 3/2 transition of Cs with ω (0) a /2π = 351.7 THz.
The local polarization of modes varies throughout the interior and exterior of the toroid. Approximate solutions for constant polarization suggest classifications as quasitransverse modes, labeled transverse electric (TE) and transverse magnetic (TM) modes, although actual solutions are not transverse. A reasonable analytic approximation for the lowest-order mode function with mode number m outside of the toroid is that of a Gaussian wrapped around the toroid's surface that decays exponentially with distance scale set by the free space wavevector 1/λ 0 = 2π/λ 0 , where d(ρ, z) = (ρ − D M /2) 2 + z 2 − D m /2 is the distance to the toroid surface, ψ(ρ, z) = arctan z ρ−D M /2 is the angle around the toroid cross-section (ψ = 0 at z = 0), and ψ 0 is a characteristic mode width (see Fig. 1a). Higher order angular modes are characterized by additional nodes along the coordinate ψ.

Cavity QED in an axisymmetric resonator
We consider a quantum model of a two-level atom at position r(t) coupled to an axisymmetric resonator shown schematically in Fig. 2. The terminology used here follows the supplemental material of [11], [13], and [12], but the general formalism can be found in additional sources (see [20], for example). As described in section 2.1, an axisymmetric resonator supports two degenerate counter-propagating whisperinggallery modes at resonance frequency ω c to which we associate the annihilation (creation) operators a and b (a † and b † ). Each traveling-wave mode has an intrinsic loss rate, κ i ; the modes are coupled via scattering at rate h. External optical access to the cavity is provided by a tapered fiber carrying input fields {a in , b in } at probe frequency ω p . Fiber fields couple to the cavity modes with an external coupling rate κ ex . The output fields of the fiber taper in each direction are the coherent sum of the input field and the leaking cavity field, [11,13]. We specialize to the situation of single-sided excitation, where b in = 0. The input field a in drives the a mode with strength ε p = i √ 2κ ex a in so that the incident photon flux is P in = a † in a in = |ε p | 2 /2κ ex . Experimentally accessible quantities are the transmitted and reflected photon fluxes, P T = a † out a out and P R = b † out b out , respectively. In experiments, data is typically presented as normalized transmission and reflection coefficients defined as T = P T /P in and R = P R /P in . In the absence of an atom, the functions T and R for the bare cavity depend on the detuning ∆ cp = ω c − ω p and the cavity rates h, κ i , and κ ex . At critical coupling, κ ex = κ 2 i + h 2 , the bare cavity T → 0 when ∆ cp = 0 [21].
The cavity modes {a, b} both couple to a single two-level atom with transition frequency ω a at location r. In the context of cQED, the atomic system is described by a single transition with frequency ω a with the associated raising and lowering operators σ + and σ − and an excited state field decay rate γ. The atomic frequency ω a ( r) may be shifted from the free-space value ω (0) a by frequency δ a ( r) due to interactions with the dielectric surface. The coupling of the traveling-wave modes {a, b} to the atomic dipole is described by the single-photon coupling rate g tw ( r) = g max tw f (ρ, z)e ±iθ , where f (ρ, z) is the cavity mode function and θ = mφ. A discussion of f (ρ, z) for the modes of a microtoroid appears in section 2.1. For an atom in motion, ω a ( r), γ( r), and g tw ( r) are spatially-dependent quantities that depend on the atomic position r(t).
To study the atom-cavity dynamics, we write the standard Jaynes-Cummings-style cQED Hamiltonian for coupled field modes [22,11]: Λ 0 Λ + Λ - Figure 2. (a) Schematic of the atom-toroid system. Coherent optical fields in the tapered fiber couple into whispering-gallery cavity modes of an axisymmetric resonator. These fields can couple to an atomic transition with rate g, scatter to the counterpropagating mode (h), escape to the environment (κ i ), or couple in/out of the fiber (κ ex ). The atom is described as a two-level system with transition frequency ω a and spontaneous emission rate γ. (b) Imaginary part of the eigenvalues Λ i of the linearized systems as a function of detuning ∆ = ω c − ω (0) a for a Cs atom at φ = π/4 and g = 60 MHz critically coupled to a cavity with parameters {κ i , h}/2π = {8, 0} MHz (Eqs. (9a)).
Following the rotating-wave approximation, we write the Hamiltonian in a frame rotating at ω p [11,13,20]: where ∆ ap ( r) = ω a ( r) − ω p and ∆ cp = ω c − ω p . Dissipation from coupling to external modes is treated using the master equation for the density operator of the system ρ: Here, κ = κ i + κ ex is the total field decay rate of each cavity mode, and 2γ( r) is the atomic dipole spontaneous emission rate, which is orientation dependent near a dielectric surface (section 4.1). The Hamiltonian (6) can be rewritten in a standing wave basis using normal modes where g A ( r) = g max f (ρ, z) cos θ, g B ( r) = g max f (ρ, z) sin θ, and g max = √ 2g max tw . In the absence of atomic coupling (g tw = 0), these normal modes are eigenstates of (6). With g tw = 0, the eigenstates of the Hamiltonian are dressed states of atom-cavity excitations. With h = 0 and g tw = 0, the atom defines a natural basis in which it couples to only a single standing wave mode. For the modes {A, B} defined above, coupling may occur predominantly, or even exclusively, to one of the two normal modes depending the azimuthal coordinate θ. For such θ, the system can be interpreted as an atom coupled to one normal mode in a traditional Jaynes-Cummings model with dressedstate splitting given by the single-photon Rabi frequency Ω (1) = 2g ≡ 2g max f (ρ, z), along with a second complementary cavity mode uncoupled to the atom. Approximately for g tw h, this interpretation is consistent for any arbitrary atomic coordinate θ. For h = 0 and comparable to κ i with a fixed phase convention (such as Im(h) = 0 used here), this decomposition is not possible for arbitrary atomic coordinate θ; the atom in general couples to both normal modes as a function of φ [11,20].
The master equation (7) can be numerically solved using a truncated number state basis for the cavity modes. Alternatively, the system is linearized by treating the atom operators σ ± as approximate bosonic harmonic oscillator operators with [σ − , σ + ] ≈ 1. For a sufficiently weak probe field, the atomic excited state population is small enough that the oscillator has negligible population above the first excited level and the harmonic approximation is quite good. As part of this linearization, we factor expectation values of normally-ordered operator products into products of operator expectation values [18,23]. Reducing operators to complex numbers suppresses coherence, but numerical calculations confirm that this approximation is accurate when calculating cavity output fields and classical forces for the weak driving power levels considered here. In particular, experiments typically utilize a photon flux P T = a † out a out ∼ 15 cts/µs corresponding to an average cavity photon population of a † a ∼ 0.1. At these photon numbers, cavity expectation values effectively factorize such that a † a ≈ a † a for the semiclassical treatment used here [24]. We use this approximation to write P T = a † out a out ≈ a † out a out , implying that we only need the complex number a out = − a in + √ 2κ ex a and its conjugate to calculate the cavity transmission at these photon numbers. This approximation is not sufficient for calculation of the g (2) (τ ) correlation function where the nonlinearities must be included [12].
The relevant equations of motion for the field amplitudes of the linearized master equation are, Time and spectral dependence of this system of equations are governed by its eigenvalues Λ i . The imaginary part of the eigenvalues as a function of detuning ∆ ≡ ∆ cp − ∆ ap = ω c − ω a are illustrated in Fig. 2b. For large ∆ |g tw |, the three eigenvalues include one atom-like eigenvalue and two cavity-like eigenvalues separated by the mode splitting h. For intermediate ∆, there is an anti-crossing of two dressed-state eigenvalues Λ ± , while the third (cavity-like) Λ 0 is uncoupled to the atom.
For a slowly-moving atom, the mode fields remain approximately in steady state as the parameters evolve with the atom trajectory r(t). Analytic steady-state solutions to (9a) for a ss and b ss are:

Optical forces on an atom in a cavity
Neutral atoms experience forces from the interaction of the atomic dipole moment with the radiation field. These optical dipole forces have a quantum mechanical interpretation as coherent photon scattering [25,26]. For a light field near resonance with the atomic dipole transition, these optical forces can be quite strong, even at the single photon level; cavity-enhanced dipole forces [18,27] have been exploited to trap [28] and localize [29] a single atom with the force generated by a single strongly-coupled photon. In this section, we discuss how the optical forces, their first-order velocity dependence, and their fluctuations are included in our semiclassical simulation.

Dipole forces
In a quantum mechanical treatment of light-matter interactions [26], the eigenstates of the system are dressed states of atom and optical field. The quantum mechanical optical force on the atom at location r can be found from the commutator of the atom momentum p with the interaction Hamiltonian H int consisting of the last two terms from the Hamiltonian (6): The gradient from the position space representation of the momentum operator p only acts on g tw ( r) and not on the field operators [30,31]. The steady-state expectation values of (11) give the dipole force on the atom in the semiclassical approximation: As described in section 2, the steady-state operator expressions are simplified by reducing expectation values of operator products to products of linearized steady-state operator expectation values. Ignoring fiber and spontaneous emission losses, an effective conservative dipole potential U d can be defined by integration of (12).

Velocity-dependent forces on an atom
Non-zero velocity effects on the force (12) are found by including a first-order velocity correction in the steady state expectation values [23,25,31]. Consider a vector of operators O whose expectation values obey a linearized equation system such as (9a). Assuming a small velocity, we expand the operator expectation values O as: where the subscripts denote the order of the velocity v in each term. If an atom is moving through these fields, then the cavity parameters depend in general on atomic position r. As r changes in time, the fields must evolve in response. Consequently, the time derivative of the expectation value evolves not only from explicit time dependence, but from atomic motion as well.
Setting the explicit time derivatives in (14) to zero, the perturbative expansion of the time derivative can be equated to the original linearized equation system. Collecting terms of each order in velocity gives an equation for the first-order term O 1 in terms of the zero-velocity steady-state solution O 0 . This procedure requires the spatial derivative of the zero-order steady-state solutions, where spatial dependence enters through the atomic transition frequency ω a ( r), the spontaneous emission rate γ( r), and the atom-cavity coupling g( r). Only terms linear in velocity are kept in the operator products of the force F ( r) in (12).
In practice, first-order velocity corrections are small in our simulation. For example, Doppler shifts arising from spatial derivatives of the cavity modes are on the order of k · v, where k is the mode wavevector. For typical azimuthal velocities of less than 0.1 m/s, the Doppler shift is less than 1 MHz. The effect becomes more significant as atoms accelerate to high velocities near the surface, but atomic level shifts from surface interactions are more significant in this regime than the Doppler shifts.

Momentum diffusion and the diffusion tensor in a cavity
Quantum fluctuations of optical forces are treated by adding a stochastic momentum diffusion contribution to the atomic velocity in the Langevin equations of motion. We calculate the diffusion tensor components used in (1), D ii , using general expressions for diffusion in an atom-cavity system generalized for the two-mode cavity of a toroid [32]: for i = x, y, z, where γ is the atomic field spontaneous decay rate. The first term represents fluctuations from spontaneous emission, the second term describes a fluctuating atomic dipole coupled to a cavity field, and the third represents a fluctuating cavity field coupled to an atomic dipole. (15) is approximated using steady-state fields calculated from the linearized solutions to the master equation (10a). Although included in the trajectory model, momentum diffusion does not significantly alter averages over ensembles of trajectories at the weak excitation levels and low atomic velocities used in the relevant experiments.

Effects of surfaces on atoms near dielectrics
In the vicinity of a material surface, the mode structure of the full electromagnetic field is modified due to the dielectric properties of nearby objects. These off-resonant radiative interactions modify the dipole decay rate of atomic states and shift electronic energy levels. This surface interaction varies spatially as the relative atom-surface configuration changes. The surface phenomena are dispersive and depend on the multilevel description of the atom's electronic structure; they are calculated using traditional perturbation theory with the full electromagnetic field without focusing on a few select modes enhanced by a cavity in cQED.

Spontaneous emission rate near a surface
When a classical oscillating dipole is placed near a surface, its radiation pattern is modified by the time-lagged reflected field from the dielectric surface. The spontaneous emission rate oscillates with distance d from the surface, which can be interpreted as interference between the radiation field of the dipole and its reflection. The variation of the emission rate depends on whether the dipole vector is parallel or perpendicular to the surface, as intuitively expected from the asymmetry of image dipole orientations of dipoles aligned parallel and perpendicular to the surface normal. For either orientation, the spontaneous emission rate features a marked increase within a wavelength of the surface due to surface evanescent modes that become available for decay for d λ 0 . We calculate the surface-modified dipole decay rates γ ( ) s (d) and γ (⊥) s (d) for a cesium atom near an SiO 2 surface following the methods of Refs. [15,33] (see Fig. 3). This calculation involves an integration of surface reflection coefficients over possible wavevectors of radiated light. The integrand depends on the dielectric function of SiO 2 evaluated at the frequency ω a of the atomic transition. The orientations refer to the alignment of a classical dipole relative to the surface plane.

Calculation of Casimir-Polder potentials
Radiative interactions with a surface are important components of motion for neutral atoms within a few hundred nm of a surface, with the potential for manipulating atomic motion through attractive [16] or repulsive forces [34]. Depending on the theoretical framework, these forces are naturally thought of as radiative self-interactions between two polarizable objects, fluctuations of virtual electromagnetic excitations, or as a manifestation of vacuum energy of the electromagnetic field. These surface interactions, represented by a conservative potential U s , are sensitive to the frequency dispersion of the electromagnetic response properties of the atoms and surfaces. For an atom located a short distance d λ 0 from a dielectric, the fluctuating dipole of the atom interacts with its own surface image dipole in the well-known nonretarded van der Waals interaction. Using only classical electrodynamics with a fluctuating dipole, the surface interaction potential is found to take the Lennard-Jones (LJ) form U LJ s = −C 3 /d 3 , where C 3 is a constant that depends on the atomic polarizability and dielectric permittivity of the surface [35,36,37,38]. At larger separations, virtual photons exchanged between atoms and surfaces cannot travel the distance in time t ∼ 1/ω due to the finite speed of light. Consequently, the interaction potential is reduced, as first calculated in the 1948 paper by Casimir and Polder [39]. The retarded surface potential takes the form U ret s = −C 4 /d 4 for a constant C 4 , where C 4 depends on both c and as this is fundamentally both a relativistic and quantum phenomenon. The full theory of surface forces for real materials with dispersive dielectric functions came with the work of Lifshitz [40,41]. This framework reduces to both the above situations for the proper limits, and, importantly, it accounts for finite temperatures, predicting a U th s ∝ d −3 potential caused by thermal photons dominant at large distances for d c/k B T [42]. In our discussion, we refer to these generalized dispersion forces as Casimir-Polder (CP) forces, whereas U LJ s , U ret s , and U th s refer to the appropriate distance limits.
In microcavity cQED, evanescent field distance scales are set by the scale length of the evanescent field, λ 0 = λ 0 /2π = 136 nm (for the Cs D2 line). The relevant distances (0 < d 300 nm) span both the LJ and retarded regimes, but are much shorter than the thermal regime (d > 5 µm). In the transition region, the limiting power laws do not fully describe U s over the relevant range of d. In our modeling, we utilize a calculation of U s with the Lifshitz approach. The Lennard-Jones, retarded, and thermal limits arise naturally from the Lifshitz formalism [42]. The potential U s enters our simulation in two ways. First, the transition frequency ω a of the two-level atomic system shifts away from the vacuum frequency by δ a = (U ex s ( r)−U g s ( r))/ , where U g s ( r) and U ex s ( r) are the surface potentials for the ground and excited states, respectively. Since the atom transitions between the ground and excited states during its passage through the mode, the average net force used in calculations is found by weighting each contribution by the steady-state atomic state populations, F s = F g s 1 − σ † ss σ ss + F ex s σ † ss σ ss . We calculate U g s and U ex s for a cesium atom near a glass SiO 2 surface using the Lifshitz approach. This calculation depends on the dispersion properties of the response functions of materials, in this case the polarizability of the Cs ground state α(ω) and the complex dielectric function (ω) of the silica surface. Modeling of these functions is discussed in Appendix A. In particular, these response functions must be evaluated on the imaginary frequency axis ω = iξ, as shown in Figure 4.
Following the method of [43], curvature of the toroid surface is implemented by treating the toroid as a cylinder with radius of curvature R = D m /2. The major axis curvature is neglected because for all relevant distances d D M /2. The resulting formula can be interpreted as a sum over discrete Matsubara frequencies ξ n = 2πnk B T / with an integration over transverse wave vectors, which we quote without derivation: [43] Here, α (iξ n ) is the atomic polarizability and r ,⊥ (iξ n , k ⊥ ) are the reflection coefficients of the dielectric material evaluated for imaginary frequency iξ n . The primed summation implies a factor of 1/2 for the n = 0 term. The reflection coefficients for the two orthogonal light polarizations are: where and (iξ n ) is the complex dielectric function evaluated for imaginary frequencies iξ n . Depending on the author, r (r ⊥ ) is sometimes referred to as r TM (r TE ). U g s is calculated by numerical evaluation of (16). U ex s is also calculated using (16), but with an additional contribution accounting for real photon exchange from the excited state with the surface, which is proportional to Re (ωa)−1 (ωa)+1 in the LJ limit [38,44]. The atom-surface potential U g s for the ground state of cesium near a SiO 2 surface is shown in Fig. 5, including calculations for both a planar and a cylindrical surface. Without the cylindrical correction, the potential approaches the LJ, retarded, and thermal limits at appropriate distance scales. For the planar dielectric, our calculation yields C 3 /h = 1178 Hz µm 3 and C 4 /h = 158 Hz µm 4 for the LJ and retarded limits. Note that the transition region between LJ and retarded regimes occurs around d ∼ 100 nm, the relevant distance scale for the experiments we are modeling. For d > D m , the perturbative method accounting for the curvature is no longer accurate [43], but at these distances, the surface forces are insignificant to atomic motion due to their steep power law fall-off. The excited state potential U ex s has a similar form to U g s .

Simulating atoms detected in real-time near microtoroids
In order for the semiclassical model to be applied to our falling atom experiments, we must simulate the atom detection processes. In particular, in [12], falling Cs atoms are detected with real-time photon counting using a field programmable gate array (FPGA), with subsequent probe modulation triggered by atom detection. A microtoroidal cavity with frequency ω c is locked near the 6S 1/2 , F = 4 → 6P 3/2 , F = 5 atomic transition of Cs at ω 0 a at desired detuning ∆ ca = ω c − ω a . Fibercavity coupling is tuned to critical coupling where the bare cavity transmission vanishes, T T min 0.01. For atom detection, a probe field at frequency ω p = ω c and flux P in ∼ 15 cts/µs is launched in the fiber taper and the transmitted output power P T is monitored by a series of single photon detectors. Photoelectric events in a running time window of length ∆t th are counted and compared to a threshold count C th . A single atom disturbs the critical coupling balance so that T /T min > 1, resulting in a burst of photons which correspond to a possible trigger event. Extensive details of the experimental procedure are given in [12] and the associated Supplementary Information.
Whereas only a single atom is required to produce a trigger, spectral and temporal data are accumulated over many thousands of trigger events since each individual atom d (μm)  Figure 5. Atom-surface potentials U g s (red) and U ex s (blue) for a cesium atom at distance d from a SiO 2 surface. The solid lines are for a planar surface whereas the dashed lines are for a curved surface with radius of curvature R = D m /2 = 1.5 µm. The limiting regimes for U g s with a planar surface are shown as dotted lines, each calculated from analytic expressions not using the Lifshitz formalism. The cylindrical surface correction weakens the potential, which is noticeable in the retarded and thermal regimes.
is only coupled to the cavity for a few microseconds. Simulation is a valuable technique to disentangle atomic dynamics from the aggregate data and offer insights into the atomic motion which underlies the experimental measurements.

Simulation procedure
Central to our simulations is the generation of a set of N representative atomic trajectories for the experimental conditions of atoms falling past a microtoroid fulfilling the criteria for real-time detection. Since experimental triggering is stochastic, the trajectory set is generated randomly as well. For each desired collection of experimental parameters P, a set of semiclassical atomic trajectories { r j (t)} P is generated that satisfies the detection trigger criteria. This ensemble is used to extract the cavity output functions T (t, ∆ ap ) and R(t, ∆ ap ). For each individual trajectory, t = 0 is defined to be the time when the trajectory is experimentally triggered by the FPGA. For each set P, N is chosen large enough for a sufficient ensemble average to be obtained for the final output functions, which is typically at least 400 unique triggered trajectories.
Within each simulation, the probe field is fixed to a given ω p . Cavity behavior is determined by the parameters ω c , h, κ i , and κ ex . h and κ i are determined from measurements of the bare cavity with no atoms present. Low-bandwidth fluctuations in κ ex and ω c from mechanical vibration and temperature locking are modeled as normally distributed random variations with standard deviations of 3 MHz and 1.5 MHz, respectively, that are fixed once for the duration of each simulated trajectory. Similar to the experimental procedure, we impose that the bare-cavity output flux is less than 0.4 cts/µs at critical coupling and on resonance. This rate would be identically zero for ∆ cp = 0 and critical coupling in the absence of these fluctuations. If the noise threshold is not met, then the particular trajectory is thrown out as it would have been in experiments.
The atomic cloud is characterized by its temperature, size, and its height above the microtoroid. Its shape is assumed to be Gaussian in each direction with parameters determined by florescence imaging. For each simulation loop, an initial atomic position r in is selected from the cloud and the initial velocity v in is selected from a Maxwell-Boltzmann distribution of temperature T . The trajectory is propagated forward in time under the influence of gravity until it crosses the toroid equatorial plane at z = 0. Only trajectories which pass within 1 µm of the toroid surface at z = 0 are kept as a candidate for detection, as other atoms couple too weakly for triggering. Once an acceptable set of initial conditions is obtained, the trajectory r(t) is calculated over a 50 µs time window around its crossing of z = 0, this time with the gravity, optical dipole forces, and surface interactions included. As the atom moves through the cavity mode, the atom-cavity coupling g, level shifts, decay rates, and forces change with position r, causing deviations of the trajectory from the preliminary free-fall trajectory. If the atom crashes into the surface of the toroid, then the coupling is set to g = 0 onwards and the trajectory effectively ends (except for random 'noise' photon counts arising from the non-zero background transmission).
Using r(t) and the steady-state expressions for the fields (section 2), we find the transmission T (t). The photon count record C i (t j ) on each photodetector i for time step t j is generated from a time-dependent Poisson process with mean count per bin of C i (t j ) = T (t j )P in ∆t, where ∆t = t j+1 − t j = 1 ns and P in is the input flux. Since the typical flux is P in ∼ 10 MHz and the timescale of quantum correlations is ∼ 10 ns, the photon count process is assumed to be Poissonian on the relevant timescale of a few hundred nanoseconds for atom detection. The count record C i (t j ) is compared to the desired threshold of C th in a time window ∆t th [12]. If the trigger condition is met, the initial conditions r in,j , v in,j , the random cavity parameters ω c and κ ex , and the random number seed used to generate W i for diffusion processes are stored for later use. The semiclassical trajectory r j (t) can be fully reconstructed from these parameters. The time coordinate is shifted so that the trigger event occurs at t = 0. This process is repeated to acquire N triggered trajectories.
Cavity output functions such as the experimentally measurable transmission T exp (t, P) for each simulation parameter set P are calculated from the set of trajectories { r j (t)}: Reflection coefficients R exp (t, P) are calculated similarly. Spectra are calculated by averaging output powers over a time window t 1 < t < t 2 for each probe frequency ω p . The times t 1 and t 2 are chosen to be the same as in our experiments, which is typically t 1 = 250 ns and t 2 = 750 ns. The set of triggered trajectories { r j (t)} is valid for a given set of conditions P and detection criteria {C th , ∆t th } until the trigger at t = 0. In experiments, the probe frequency ω p can be changed in power and detuning upon FPGA trigger. Although the same set of trajectories is valid before t = 0 for each detuning, the trajectory set must be recalculated for t > 0 for each probe detuning to mimic experimental conditions for spectral measurements. A numerical solution of the master equation in a number state basis is used for calculation of T (t) in (20); the linearized model is only used to calculate the trajectory r(t) and efficiently generate triggers. Whereas experiments give access only to ensemble averaged output functions, simulations contain the full trajectory paths. Provided that the simulation offers a reasonable approximation of the true ensemble of trajectories, then these results provide a window into the atomic dynamics underlying the cQED measurement of falling atoms which are not readily clear from observations.

Simulation distributions
The experimentally measurable cavity transmission T exp (t) is obtained in (20) as an average over the trajectory set { r j (t)} at each time t. Eq. (20) can formally be written as an integration over the probability distribution of coupling constants at time t, p t (g, θ), for the given experimental parameters P: T exp (t, P) = dg dθ T (g, θ, P)p t (g, θ) The function T (g, θ, P) is shown in Fig. 6 for the parameters P relevant to experiments, specifically with ∆ ca /2π = 0, 60 MHz. For this discussion, we assume all frequencies are fixed and neglect surface shifts. In this perspective, T exp (t) is not directly related to the trajectory set { r j (t)} but rather the probability distribution p t (g, θ) at time t. The time dependence of p t evolves based on the underlying trajectory ensemble.
b.  Figure 6. Plots of T (g tw , θ, P) for (a) and (b) calculated numerically from (7). Atoms with higher g tw generally have higher T and a larger probability for detection. The variation of T with θ is evident, with a different periodicity for the two cavity detunings.
It is instructive to consider the probability distribution p t (g, θ) in more detail since it is the formal output of the simulations. We consider only the distribution p t=0 (g) over the coupling parameter g at the trigger time t = 0 by integrating out the angular dependence. Through a reasonably simple analytic model (detailed in Appendix B), we calculate p t=0 (g) and compare to the results of the semiclassical simulation, which includes dipole and surface forces (Fig. 7). For a cavity on resonance with the atom transition, ∆ ca /2π = 0, the analytic model agrees well with a simulation when dipole and surface forces are not included. In this case, atom trajectories are nearly straight and vertical near the toroid, and the approximations of Appendix B are sufficient. When the full forces are included in the semiclassical model, the additional forces shift the distribution toward lower g. This effect is more significant for ∆ ca /2π = 60 MHz. The corresponding experimental cQED spectra confirm that the semiclassical model with dipole and surface forces is necessary to reproduce spectral features in the real-time experiments (Fig. 7c).  The cavity transmission T varies as a function of the atomic azimuthal coordinate θ = mφ, as evident in Fig. 6. This biases atomic detection towards specific locations around the toroid and leads to a non-uniform angular distribution p t=0 (θ) for atom location at detection. Fig. 8 shows distributions of the atomic angular coordinate at the detection trigger t = 0 for three simulation conditions relevant to the experiments of [12]. Although averaged spectra do not explicitly measure the coordinate θ, these simulation makes clear that trajectories passing through certain regions around the toroid are preferentially detected. The phase of the cavity output field depends on θ, suggesting the possibility for future experiments to measure the distribution of Fig. 8.  Figure 8. Probability distribution p t=0 (θ) of atomic azimuthal angle θ = mφ mod 2π at transit detection time t = 0 presented as histograms of simulation runs. Shown are the cases for cavity detunings (a) ∆ ca = 0 and (b) ∆ ca /2π = ±40 MHz. Normalization is such that the sum across all θ is unity.

Simulated trajectories
We now turn to the simulated trajectories { r j (t)}. In contrast to experiments, in simulations we have the capability of turning certain forces selectively on and off. In particular, we can adjust the surface potential U s and the dipole forces, referred to symbolically as U d (despite them not being strictly derivable from a potential). To investigate the effects these optical phenomena have on atomic trajectories, we run simulations for four cases: the full semiclassical model, the model without surface forces (U s = 0), the model without dipole forces (U d = 0), and the model without any radiative forces (U d = U s = 0).
Considering conditions relevant to [12], we plot simulations for two sets of experimental parameters P 1,2 in Fig. 9. For P 1 , the cavity is detuned to the red, whereas the cavity is blue-detuned in P 2 (∆ ca /2π = −40 MHz for P 1 and +40 MHz for P 2 ). In each set of conditions, the probe field is on resonance with the cavity for high signal-to-noise atom detection (∆ cp = 0) and the average bare-cavity mode population of a is ≈ 0.05 photons. The toroid cavity parameters are those of [12], {g max , h, κ in , κ ex }/2π = {100, 11, 13, 17} MHz. Comparing the full model, we see that trajectories for P 1 primarily crash into the surface, whereas those from P 2 both crash and are repelled from the toroid. This asymmetry is due to the repulsive or attractive dipole force for different probe detunings relative to the atomic transition. The largest effect of turning surface forces off is seen in the blue-detuned trajectories, which have a lower crash rate when U s = 0. With U d = 0, both P 1 and P 2 trajectories look nominally the same; the detuning ∆ ca only affects cQED spectra, with a minor imperceptible effect arising from CP potentials initially shifting the atomic transition either closer to (red) or further from (blue) the cavity field.
In addition to the qualitative differences in detected atom trajectories summarized here, the effects of U d and U s are evident in the experimental quantities T exp (t) and R exp (t). Since here we focus specifically on trajectory calculations, the reader is referred to [12] for detailed comparisons of spectral and temporal simulations to experimental data.  are projected onto the two-dimensional ρ − z plane for all conditions. Magenta trajectories represent un-triggered atoms, blue paths are detected atoms for t < 0 and red paths represent atom trajectories after the trigger for t > 0.

Trapping atoms in the evanescent field of a microtoroid
Our trajectory simulation can be extended to study trapping of atoms in a two-color evanescent far off-resonant trap (eFORT) near a microtoroidal resonator. An evanescent field trap takes advantage of the wavelength dependence of scale lengths for the optical dipole force of two optical fields with frequencies far-detuned from the atomic transition to limit scattering [45,46,47]. The relative powers of the two fields are set so that near the surface, the blue-detuned, repulsive field is stronger than a red-detuned attractive field. As each field falls off with a decay constant of roughly λ = 2π/λ, at some distance, the red, attractive field will dominate and the atom will be attracted to the surface forming a potential minimum. Recently, evanescent fields have been harnessed to trap atoms in a two-color eFORT around a tapered optical fiber [48], where the fiber enables efficient optical access to deliver both high intensity trapping fields and weaker probe fields to the trapped atoms in a single structure. The tapered fiber can be positioned as desired, bringing the trapped atoms near a device for atomic coupling. The tapered nanofiber eFORT is a remarkable achievement toward integrating atom traps with solid-state resonators, but the nanofiber scheme does not allow direct integration with a cavity for achieving strong, coherent coupling between light and trapped atoms. Another disadvantage is that trap depth is limited by the large total power required to achieve trapping with evanescent fields. The high quality factors and monolithic structure of WGM resonators allow evanescent field traps free from these problems while maintaining efficient optical access from tapered fiber coupling. Two-color evanescent field traps in WGM resonators have been analyzed in detail for spheres [49] and microdisks [50]. In this section, we extend our simulations of atoms in the evanescent field of a microtoroid to an eFORT that can capture single falling Cs atoms triggered upon an atom detection event.
Unlike nanofibers, a microtoroid cannot be placed directly in a magneto-optical trap for a source of cold atoms. As shown in [12], we have the experimental capability to detect a single atom falling by a microtoroid and trigger optical fields while that atom remains coupled to the cavity mode. The semiclassical simulations described here are ideal for investigating the capture of falling atoms in a trap triggered upon experimental atom detection.
We add an additional eFORT potential U t to our semiclassical trajectory model in addition to the dipole forces and surface potential U s . For our simulation, U t is formed from a red (blue)-detuned mode near 898 nm (848 nm) with powers ∼ 50 µW to give a trap depth of ∼ 1.5 mK at d ∼ 150 nm from the surface (Fig. 10a). The red (blue) fields interact primarily with the 6S 1/2 → 6P 1/2 (6S 1/2 → 6P 3/2 ) transition. The trap depth is limited by the total power in vacuum that can propagate in the tapered fibers of [12]. Power handling can be improved with specific attention to taper cleanliness, so with experimental care the trap depth can be increased reasonably from the discussion here, although we simulate under the conditions given to illustrate that this trap is already experimentally accessible. c.
d. The difference in vertical scale lengths (ψ 0 in (4)) for modes of different wavelength leads to a trap that is not fully confined if both the red-and blue-detuned trap modes are of the lowest order (as in Fig. 2.1b). As |ψ| increases, the repulsive blue-detuned light weakens faster than the red-detuned field, and atoms can crash into the toroid surface. This problem is alleviated by exciting a higher-order mode for the 898 nm light, as shown in Fig. 10b. The modal pattern confines atoms near z = 0 and prevents trap leakage along ψ. This problem is not present in the microdisk eFORT of [50] because the optical mode extent is determined by structural confinement and not the optical scale length. Use of a higher-order mode was also used to form an atom-gallery in a microsphere [49].
During the detection phase of the simulation, U t = 0. At t = 0 conditioned on an atom detection trigger, U t is turned on. The kinetic energy of an atom with typical fall velocity of v ∼ 0.2 m/s is equivalent to 0.3 mK, so a 1.5 mK trap is sufficiently deep to capture an atom if it is triggered near the trap potential minimum. Defining a trapped trajectory to be one such that the atom has g/2π > 5 MHz at t = 10 µs, approximately 25% of triggered atom trajectories are captured when the trapping potential is turned on. Simulated trapping times exceed 50 µs, limited not by heating from trapping light but by the radiation pressure from the unbalanced traveling whispering-gallery modes of nearly-resonant optical probe field. This probe field can be turned off so that the atoms remain trapped beyond the simulation time.
In contrast to the standing-wave structure of a typical eFORT or Fabry-Perot cavity trap [51], microtoroidal resonators offer the tantalizing possibility of radially confining an atom in a circular orbit around the toroid [49,52]. The U t = 0 outlined here does not confine the atoms azimuthally, forming circular atom-gallery orbits around the microtoroid [52] (Fig. 10c,d). In the same manner as [48], a localized trap can be achieved by exciting a red-detuned standing wave for three-dimensional trap confinement.
This trapping simulation outlines how real-time atom detection can be utilized to trap a falling atom in a microtoroidal eFORT. In practice, microtoroidal traps present some serious practical challenges. Notably, because the trap quality is sensitive to the particular whispering-gallery mode, the excited optical mode must be experimentally controlled. The success of an eFORT for Cs atoms around a tapered nanofiber [48] strongly suggests that similar trap performance might be achieved for an eFORT around a high-Q WGM cavity, localizing atoms in a region of strong coupling to a microresonator.

Conclusion
We have presented simulations of atomic motion near a dielectric surface in the regime of strong coupling to a cavity with weak atomic excitation. As required by experimental distance scales, this simulation includes surface interactions, which manifest through transition level shifts and center-of-mass Casimir-Polder forces. Analysis of the simulated trajectories gives insight into the atomic motion underlying experimental measurements of ensemble-averaged spectral and temporal measurements for single atoms detected in real-time. We have adapted our simulations to investigate the capturing of atoms in an evanescent field far off-resonant optical trap in a microtoroid. Our simulations suggest that falling atoms can be captured into an eFORT around a microtoroid, offering an experimental route towards trapping a single atom in atomchip trap in a regime with both strong cQED interactions and significant Casimir-Polder forces simultaneously. In this system, the sensitivity afforded by coherent cQED can be used not only for atom-chip devices, but also as a tool for precision measurements of optical phenomena near surfaces.
The total atomic polarizability is composed of contributions from valence electron transitions (α v ) and high-energy electron transitions from the core shells to the continuum (α c ), such that α = α v + α c . The valence polarizability α v constitutes 96% of the total static polarizability [54] in Cs, with α c only significant at high frequencies.
For simplicity, all core electron transitions are lumped into a single high-frequency term of the form used in (A.2). This term contains two free parameters, f core and ω core , which are found from the following two conditions. Using the calculation of α v (ω) for the Cs ground state, we enforce that the ground state static polarizability α(ω → 0) matches the known value calculated theoretically [58] α(0) = 5.942 × 10 −23 cm 3 . We also ensure that the ground state LJ constant for a Cs atom near a metallic surface agrees with the known value [54,59] C 3 = − 4πd 3 ∞ 0 α(iξ)dξ = 4.4 · h kHz µm 3 . These conditions are sufficient to fix the two free parameters in α c (ω) for this single oscillator core model, although the high-frequency structure of the core polarizability is lost. For the excited state calculation, we use the same α c (ω).

Appendix B. Analytic model of falling atom detection distributions
Here we develop an analytic model of the distribution p fall (g, θ) of coupling parameters g and azimuthal coordinate θ = mφ. Atoms are assumed to fall at constant vertical velocity with no forces, in contrast to the more complete semiclassical trajectories used in this manuscript to generate p t (g). An abbreviated description of this model appears in the Supplementary Information of [12].
The linearized steady-state cavity transmission T (∆ ap , g( r)) is a known function of ∆ ap and r. We only consider the lowest order mode where the cavity mode function is approximately Gaussian in z and exponential in distance from the surface d. The approximate temporal behavior of the coupling constant g for a single trajectory is, g(ρ, z(t)) = g c (ρ)e −(z(t)/z 0 ) 2 . (B.1) where g c (ρ) is the maximum value of the g at the closest approach of its trajectory (z = 0), z 0 is a characteristic width assumed to be independent of ρ, and z(t) = −vt. g c (ρ) decays exponentially from the maximum g max at the toroid surface, g c (ρ) ∼ g max e −(ρ−Dp)/λ 0 . The transmission T and hence the detection probability depend on θ; in general, if atoms fall uniformly around the toroid, the most numerous trajectories detected will be at the values of θ which maximize T (θ) for the cavity parameters of interest (θ = π/2 for ∆ ca /2π = +40 MHz, for example, as in Fig. 8).
The probability density function for the full ensemble of detected falling atoms p fall (g, θ) can be estimated as the product of the probability of any atom having a particular g and the probability of a trigger event occurring for an atom with coupling g, p fall (g, θ) ∼ p atom (g)p trigger (g, θ). (B.2) An atom transit is triggered when the total detected photon counts exceeds a threshold number, C th , within a detection time window ∆t th . For a probe beam of input flux P in , the mean counts in this window are C = T (g, θ)P in ∆t th . This expression assumes that the atom is moving slowly so that the T (g, θ) at trigger event is the only T (g, θ) that contributes to the detection probability. The detection probability p trigger (g, θ) is estimated from a Poisson distribution of mean count C. From (B.1), p atom (g) can be written as a product of the probability p(g|g c ) of an atom in a trajectory with a given g c to have coupling g and the probability of a trajectory to have that g c , p max (g c ), integrated over all g c , The integral has limits from g to g max since g c cannot be smaller than g. For atoms falling uniformly over the ρ − φ plane, p max (g c ) dg c is proportional to the area of a ring of radius ρ and thickness dρ, p max (g c ) dg c ∼ 2πρ dρ. Using g c (ρ) ∼ e −(ρ−Dp/2)/λ 0 , dgc gc ∼ − dρ λ 0 . Hence, p max (g c ) ∼ 1/g c for (ρ − D p /2) D p /2. To find p(g|g c ) we note that that the probability is proportional to the time an atom in the trajectory is at a particular g. From (B.1) for a constant velocity v, this trajectory is Gaussian and the relative probability must be proportional to dz. Finding the differential as a function of g gives p(g|g c ) ∝ dz ∼ This result diverges as g goes to zero since there are infinite transits with small g c and infinite time for atoms with small g for any transit regardless of g c for t → ±∞. This divergence is not problematic in calculating (B.2) since p trigger (g, θ) cuts off for low g faster than the logarithmic divergence in p atom (g). The spectrum for given experimental parameters as a function of probe detuning ∆ ap = ω p − ω