Channeling of spontaneous emission from an atom into the fundamental and higher-order modes of a vacuum-clad ultrathin optical fiber

We study spontaneous emission from a rubidium atom into the fundamental and higher-order modes of a vacuum-clad ultrathin optical fiber. We show that the spontaneous emission rate depends on the magnetic sublevel, the type of modes, the orientation of the quantization axis, and the fiber radius. We find that the rate of spontaneous emission into the TE modes is always symmetric with respect to the propagation directions. Directional asymmetry of spontaneous emission into other modes may appear when the quantization axis does not lie in the meridional plane containing the position of the atom. When the fiber radius is in the range from 330 nm to 450 nm, the spontaneous emission into the HE$_{21}$ modes is stronger than into the HE$_{11}$, TE$_{01}$, and TM$_{01}$ modes. At the cutoff for higher-order modes, the rates of spontaneous emission into guided and radiation modes undergo steep variations, which are caused by the changes in the mode structure. We show that the spontaneous emission from the upper level of the cyclic transition into the TM modes is unidirectional when the quantization axis lies at an appropriate azimuthal angle in the fiber transverse plane.


I. INTRODUCTION
Optical fibers can be tapered to a diameter comparable to or smaller than the wavelength of light [1][2][3]. Due to the tapering, the original core almost vanishes and the refractive indices that determine the guiding properties of the tapered fiber are those of the original silica cladding and the surrounding vacuum. Since the diameters of such tapered fibers are in the range of a few hundred nanometers, they are usually called nanofibers. When the radius of the fiber is small enough, it can support only a single mode in the optical region of frequency.
Tapered fibers can also be fabricated with slightly larger diameters or larger refractive indices so that they can support not only the fundamental HE 11 mode but also several higher-order modes. Compared to the HE 11 mode, the higher-order modes have larger cutoff size parameters and more complex intensity, phase, and polar-ization distributions. For ease of reference, the vacuumclad tapered fibers that can support the fundamental mode and several higher-order modes are called ultrathin optical fibers in this paper.
It has been shown that ultrathin optical fibers with higher-order modes can be used to trap, probe, and manipulate atoms, molecules, and particles [28][29][30][31][32][33][34]. The excitation of higher-order modes has been studied [35,36]. The production of ultrathin fibers with higher-order modes [37][38][39] and the experimental studies on the interaction with atoms [40] or particles [41,42] have been reported. The possibility to control and manipulate individual atoms near an ultrathin fiber can also find applications for quantum information.
The interaction between guided light and atoms is of academic and practical interest. Many applications require a deep understanding and an effective control of spontaneous emission of atoms near an ultrathin optical fiber. Radiative decay of an atom in the vicinity of a nanofiber has been studied in the context of a two-level atom [43][44][45] as well as a realistic multilevel atom with a hyperfine structure of energy levels [46,47]. The parameters for the decay of populations [43][44][45][46][47] and cross-level coherences [46][47][48] have been calculated.
It has been shown that spontaneous emission and scattering from an atom with a circular dipole near a nanofiber can be asymmetric with respect to the opposite axial propagation directions [58][59][60][61][62][63]. These directional effects are the signatures of spin-orbit coupling of light [64][65][66][67][68] carrying transverse spin angular momentum [66,69]. They are due to the existence of a nonzero longitudinal component of the nanofiber guided field, which oscillates in phase quadrature with respect to the radial transverse component. The possibility of directional emission from an atom into propagating radiation modes of a nanofiber and the possibility of generation of a lateral force on the atom have been reported [62]. The directiondependent emission and absorption of photons lead to chiral quantum optics [70].
Spontaneous emission from a multilevel atom into the fundamental and higher-order modes of an ultrathin fiber has been studied by Masalov and Minogin [71]. They have found that the decay rates into the higher-order modes can be significantly larger than into the fundamental mode. Their calculations were limited to single transitions and single polarizations. However, all types of transitions and polarizations must be accounted for in a realistic situation. In addition, in Ref. [71] the fiber axis was used as the quantization axis and consequently no direction dependencies of the rates could be observed. Moreover, emission into radiation modes was not considered in Ref. [71].
The aim of the present paper is to investigate directional spontaneous emission from a multilevel atom with an arbitrary quantization axis into an ultrathin fiber. We calculate the rates of spontaneous emission into the fundamental and higher-order guided modes propagating in a given direction. We also calculate the rate of spontaneous emission into radiation modes.
The paper is organized as follows. In Sec. II we describe the interaction of an alkali-metal atom with the electromagnetic field in the presence of an ultrathin optical fiber. Section III is devoted to the basic characteristics of spontaneous emission of the multilevel atom. In Sec. IV we present numerical results. Our conclusions are given in Sec. V.

II. MODEL AND HAMILTONIAN
We consider a multilevel alkali-metal atom trapped in the vicinity of a vacuum-clad ultrathin optical fiber [see Fig. 1(a)]. We use Cartesian coordinates {x, y, z}, where z is the coordinate along the fiber axis, and also cylindrical coordinates {r, ϕ, z}, where r and ϕ are the polar coordinates in the fiber transverse plane xy. The energy levels of the atom are specified in a Cartesian coordinate system {x Q , y Q , z Q }, where z Q is the direction of the quantization axis.
To be concrete, we assume that the atom is 87 Rb. We work with the D 2 line of the rubidium atom, which corresponds to the electric dipole transition from the excited state 5P 3/2 to the ground state 5S 1/2 [see Fig. 1(b)] [72]. We introduce the notations |e = |J ′ F ′ M ′ and |g = |JF M for the magnetic sublevels of the hyperfinestructure (hfs) levels of the excited state and the ground state, respectively. Here, J and J ′ are the total electronic angular momenta, F and F ′ are the total atomic angular momenta, and M and M ′ are the magnetic quantum numbers. We denote the energies of these sublevels as hω e andhω g . The schematic of the hfs levels of the D 2 line of the rubidium-87 atom is illustrated in Fig. 1(b).
We introduce the notation d eg = e|D|g for the dipole matrix element of the transition |e ↔ |g , where D is the electric dipole operator. In the atomic quantization coordinate system {x Q , y Q , z Q }, the spherical components q = 0, ±1 of the dipole matrix element d eg are given by the expression [73] Here, the array in the curly braces is a 6j symbol, the array in the parentheses is a 3j symbol, I is the nuclear spin, and J ′ D J is the reduced electric dipole matrix element in the J basis. Note that d qQ is nonzero only for We assume that the fiber has a cylindrical silica core of radius a and refractive index n 1 and an infinite vacuum cladding of refractive index n 2 = 1. We retain the silica dispersion and at the frequency of the rubidium D 2 line the refractive index n 1 of the fiber is taken as 1.4537. The positive-frequency part E (+) of the electric component of the field can be decomposed into the contributions E (+) g and E (+) r from guided and radiation modes, respectively, as In view of the very low losses of silica in the wavelength range of interest, we neglect material absorption. We follow the continuum field quantization procedures presented in [74]. Regarding the guided modes, we assume that the fiber supports the fundamental HE 11 mode and a few higher-order modes [75] in a finite bandwidth around the central frequency ω 0 = ω e − ω g of the rubidium-87 D 2 line. We label each guided mode in this bandwidth by an index µ = (ω, N, f, p). Here, ω is the mode frequency, the notation N = HE lm , EH lm , TE 0m , or TM 0m stands for the mode type, with l = 1, 2, . . . and m = 1, 2, . . . being the azimuthal and radial mode orders, respectively, and the index f = +1 or −1 denotes respectively the forward or backward propagation direction along the fiber axis z. The HE lm and EH lm modes are hybrid modes. For these modes, the azimuthal order is l = 0, and the index p is equal to +1 or −1, indicating the counterclockwise or clockwise circulation direction of the helical phasefront. The TE 0m and TM 0m modes are transverse electric and magnetic modes. For these modes, the azimuthal mode order is l = 0 and, hence, the mode polarization is single and the polarization index p can take an arbitrary value. For convenience, we assign the value p = 0 to the polarization index p for TE 0m and TM 0m modes. In the interaction picture, the quantum expression for the positive-frequency part E (+) g of the electric component of the field in guided modes is [46] E (+) Here, e (µ) = e (µ) (r, ϕ) is the profile function of the guided mode µ in the classical problem, a µ is the corresponding photon annihilation operator, µ = N f p ∞ 0 dω is the generalized summation over the guided modes, β is the longitudinal propagation constant, and β ′ is the derivative of β with respect to ω. The constant β is determined by the fiber eigenvalue equation [75]. The operators a µ and a † µ satisfy the continuous-mode bosonic commutation rules [a µ , a † where n ref (r) = n 1 for r < a and n 2 for r > a. The explicit expressions for the profile functions e (µ) of guided modes are given in Refs. [75,76] and are summarized in Appendix A. For a hybrid mode N = HE lm and EH lm with the propagation direction f and the phase circulation direction p, the profile function is given in the cylindrical coordinates as (5) where e r , e ϕ , and e z are given by Eqs. (A10) and (A11) for β > 0 and l > 0. For a TE 0m mode with the propagation direction f , the profile function can be written as where the only nonzero cylindrical component e ϕ is given by the second expressions in Eqs. (A17) and (A18). For a TM mode with the propagation direction f , we have where the components e r and e z are given by the first and third expressions in Eqs. (A22) and (A23) for β > 0. An important property of the mode functions of hybrid and TM modes is that the longitudinal component e z is nonvanishing and in quadrature (π/2 out of phase) with the radial component e r .
In the case of radiation modes, the longitudinal propagation constant β for each value of the frequency ω can vary continuously, from −k to k (with k = ω/c). We label each radiation mode by an index ν = (ω, β, l, p), where l = 0, ±1, ±2, . . . is the mode order and p = +, − is the mode polarization. In the interaction picture, the quantum expression for the positive-frequency part E (+) r of the electric component of the field in radiation modes is [46] E (+) Here, e (ν) = e (ν) (r, ϕ) is the profile function of the radiation mode ν in the classical problem, a ν is the corresponding photon annihilation operator, and ν = lp ∞ 0 dω k −k dβ is the generalized summation over the radiation modes. The operators a ν and a † ν satisfy the continuous-mode bosonic commutation rules [a ν , a † ν ′ ] = δ(ω − ω ′ )δ(β − β ′ )δ ll ′ δ pp ′ . In deriving Eq. (8), we have used the normalization condition The explicit expressions for the mode functions e (ν) are given in Refs. [75,76] and are summarized in Appendix B.
Assume that the atom is positioned at a point (r, ϕ, z). The Hamiltonian for the atom-field interaction in the dipole and rotating-wave approximations is given by where the notations α = µ, ν and α = µ + ν stand for the general mode index and the complete mode summation, respectively, and the operators σ ge = |g e| and σ † ge = σ eg = |e g| describe the downward and upward transitions, respectively. The coefficients characterize the coupling of the atomic transition e ↔ g with the guided mode µ and the radiation mode ν. The notation ω eg = ω e − ω g stands for the atomic transition frequency. We note that, for |e = |J ′ F ′ M ′ and |g = |JF M , the scalar product of the atomic dipole vector d eg and the field vector e (α) can be expressed as d eg · e (α) = (−1) q d qQ e zQ , and e (α) Let θ Q be the angle between the quantization axis z Q and the fiber axis z [see Fig. 1(a)]. Assume that the plane (z, z Q ) intersects with the fiber transverse plane xy at a line ζ. Let ϕ Q be the azimuthal angle between ζ and x. We choose the axes x Q and y Q such that x Q is in the plane (z Q , z) and y Q is in the plane (x, y). Then, the transformation for the field vector e (α) from the coordinate system {x, y, z} to the coordinate system {x Q , y Q , z Q } is given by the equations e (α) xQ = (e (α) x cos ϕ Q + e (α) y sin ϕ Q ) cos θ Q − e (α) z sin θ Q , e (α) yQ = −e (α) x sin ϕ Q + e (α) y cos ϕ Q , e (α) zQ = (e (α) x cos ϕ Q + e (α) y sin ϕ Q ) sin θ Q + e (α) z cos θ Q .
The relations between the Cartesian-coordinate vector components e

III. SPONTANEOUS EMISSION OF THE ATOM
In this section, we study spontaneous emission of the multilevel atom. We assume that the field is initially in the vacuum state |0 . In this case, the time evolution of the reduced density operator ρ of the atom is governed by the equations [46] where the coefficients characterize the spontaneous emission process.
In Eqs. (14), the set of coefficients γ (g) ee ′ gg ′ and γ (g) ee ′ describes spontaneous emission into guided modes, and the set of coefficients γ (r) ee ′ gg ′ and γ (r) ee ′ describes spontaneous emission into radiation modes. The expressions for these coefficients are given as [46] γ (g) and γ (r) where G N f peg ≡ G ω0N f peg and G βlpeg ≡ G ω0βlpeg are the coupling coefficients for the resonant guided and radiation modes, respectively. The diagonal decay coefficients γ ee are the rates of spontaneous emission from the magnetic sublevel |e of the atom into guided and radiation modes, respectively. The total decay rate for the population of the sublevel |e is The rate of spontaneous emission from the magnetic sublevel |e of the atom into the guided modes N = HE lm , EH lm , TE 0m , or TM 0m is given by It is clear that We note that the density-matrix equations (13) are in agreement with those used in the treatments for the excitation of a multilevel atom by light of arbitrary polarization [48,[77][78][79][80][81][82]. Equations (13) can, in principle, be used for an arbitrary (degenerate and nondegenerate) multilevel atom. The tensor nature of the Zeeman sublevels and the hfs levels of a realistic alkali-metal atom is expressed by Eq. (1) for the spherical tensor components of the atomic dipole matrix elements d eg . These quantities enter Eqs. (13) through expressions (11) for the coupling coefficients G µeg and G νeg . Unlike the case of the atom-field system in free space [48], the presence of the nanofiber modifies the decay rates γ e and leads to the appearance of the cross-level decay coefficients γ ee ′ (with e = e ′ ) in Eqs. (13) (see [46]).
We introduce the notation which stands for the rate of spontaneous emission into the guided modes N f p via the transition |e → |g . The rate of spontaneous emission from the sublevel |e of the atom into the guided modes N with the propagation direction f is given by The rate of spontaneous emission into all types of guided modes propagating in the direction f is given by For TE modes, the profile function for the electric part of the field does not depend on the propagation direction f [see Eq. (6)]. Therefore, the rates γ For hybrid and TM modes, the longitudinal component e (ωN f p) z of the field is nonvanishing and has opposite signs for opposite propagation directions [see Eqs. (5) and (7)]. Therefore, the rates γ wheref = −f . Hence, we find the relation which yields γ . More generally, we find that γ do not depend on f when the quantization axis z Q lies in the meridional plane containing the position of the atom. In order to show this directional independence, we assume that the atom is on the x axis and the quantization axis z Q lies in the zx plane, that is, ϕ Q = 0. Then, for hybrid and TM modes with the profile functions (5) and (7), Eqs. (12) yield According to Appendix A, for an appropriate choice of the normalization constant, e z and e ϕ are real numbers and e r is an imaginary number. Hence, we can show that the absolute values |e qQ | of the spherical tensor components of the field in the coordinate system {x Q , y Q , z Q } do not depend on f . On the other hand, the dipole matrix element d eg has a single nonzero spherical tensor component d qQ , which is a real number. Consequently, the absolute value of the scalar product d eg · e (µ) is |d eg · e (µ) | = |d qQ ||e and where |ḡ = |F, −M . With the help of the relations (26) and (27), we can also show that Thus, the rates γ In order to get insight into the direction dependencies of the spontaneous emission rates, we consider the rate γ for a given transition |e → |g . When we follow the procedure of Ref. [83], we can decompose this rate as where Here, the notation {A ⊗ B} 2 stands for the irreducible tensor product of rank 2 of arbitrary complex vectors A and B. The quantities γ are called the scalar, vector, and tensor components of the rate γ With the help of the first relation in Eqs. (26), we can show that γ Thus, the direction dependence of the rate γ of the rate can be considered as a result of the interaction between the effective magnetic dipole and the effective magnetic field. Due to spin-orbit coupling of light [64][65][66][67][68], a reverse of the propagation direction leads to a reverse of the spin density of light and, consequently, to a reverse of the vector component γ is a consequence of the fact that the longitudinal component e z of the guided field is not zero.

IV. NUMERICAL RESULTS
In this section, we demonstrate the results of numerical calculations for the decay characteristics of the magnetic sublevels of the excited state 5P 3/2 of a rubidium-87 atom in the presence of an ultrathin optical fiber. The atomic transitions between this state and the ground state 5S 1/2 correspond to the D 2 line and have a wavelength λ 0 = 780 nm. For simplicity, we show only the results of calculations for the spontaneous emission rates γ e of the sublevels |e = |F ′ M ′ = |J ′ F ′ M ′ and their components.
A. Dependencies of the rates on the radial distance In this subsection, we study the dependencies of the rates on the radial distance r. For simplicity, we consider the case where the fiber axis z is used as the quantization axis. In this case, none of the rates depend on the azimuthal angle ϕ. In addition, the decay rates of the sublevels with the magnetic quantum numbers M ′ and −M ′ are the same.
We show in Fig. 2 the radial dependencies of the rates γ (N ) e of spontaneous emission from different magnetic sublevels of the hfs level 5P 3/2 F ′ = 3 of the rubidium atom into different guided modes. The fiber radius is chosen to be a = 400 nm. For the wavelength λ 0 = 780 nm, this fiber can support the HE 11 , TE 01 , TM 01 , and HE 21 modes. According to Fig. 2, the presence of the fiber leads to substantial decay rates into guided modes. Comparison between the different parts of the figure shows that the emission into the HE 21 modes is stronger than into the HE 11 , TE 01 , and TM 01 modes. We observe that different magnetic sublevels have different decay rates, unlike the case of alkali-metal atoms in free space. The rates of spontaneous emission from the outermost magnetic sublevels |F ′ = 3, M ′ = ±3 (red lines) into guided modes are larger than those from the other sublevels. This indicates that the polarization profiles of the guided modes are more favorable to the σ ± transitions than the π transition. The rates of spontaneous emission into guided modes are largest when the atom is positioned on the fiber surface. When the atom is far away from the fiber, γ We note that our results presented in Fig. 2 do not agree quantitatively with the results of Masalov and Minogin [71]. Indeed, the ratio between the rates of emission from the outermost levels into the HE 21 and HE 11 modes at the distance r/a = 1 is equal to about 3 in Fig. 2 but is equal to about 8 in the calculations of Ref. [71]. One of the reasons for the discrepancy is that they considered 85 Rb, while we study 87 Rb. Another reason is that they limited their calculations to atomic transitions and guided modes with a single type of polarization, while we include all atomic transitions and field modes in our treatment. The most important reason for the discrepancy is that Eq. (16) of Ref. [71] is not accurate.
We show in Fig. 3 the radial dependencies of the spontaneous emission rates γ (g) e , γ (r) e , and γ e from different magnetic sublevels of the hfs level 5P 3/2 F ′ = 3 into guided modes, radiation modes, and both types of modes, respectively. We observe from Fig. 3(a) that the rates γ (g) e for the outermost sublevels M ′ = ±3 (red lines) are larger than for the other sublevels. When the radial distance r is not too large, the rates γ constructive and destructive interference due to reflections from the fiber surface [45]. Due to the interference, the total rate γ e can become slightly smaller than γ 0 in some regions outside the fiber [see Fig. 3(c)].
r/a We show in Fig. 4 the radial dependencies of the fractional rates η We show in Fig. 5 the radial dependencies of the fractional rates η e = γ Note that the hfs level 5P 3/2 F ′ = 0 is a singlet state, |F ′ = 0, M ′ = 0 , which is equally coupled to the sublevels |F = 1, M = 0, ±1 of the hfs level F = 1 of the ground state 5S 1/2 . Therefore, the decay rate for the state |F ′ = 0, M ′ = 0 is equal to the average decay rate for an ensemble of two-level emitters with dipoles oriented randomly in space. We show in Fig. 6 the radial dependencies of the spontaneous emission rates γ (g) e , γ (r) e , and γ e from the hfs level 5P 3/2 F ′ = 0 into guided modes, radiation modes, and both types of modes.

B. Dependencies of the rates on the fiber radius
In this subsection, we study the dependencies of the decay rates on the fiber radius a. We again use the fiber axis z as the quantization axis.
In Fig. 7, we show the rates γ (N ) e of spontaneous emission from different magnetic sublevels of the hfs level 5P 3/2 F ′ = 3 into different guided modes as functions of the fiber radius a. We observe from the figure that the rates γ sion from the atom into the fundamental HE 11 modes is strongest when a is around 180 nm. For a given fiber radius a in the range from 330 nm to 450 nm (the sizes that are typically achieved experimentally), the emission into the HE 21 modes is stronger than into the TM 01 , TE 01 , and HE 11 modes. When the atom is positioned on the fiber surface, the rates γ We plot in Fig. 9 the fractional rates η We show in Fig. 10 the fractional rates η e = γ We plot in Fig. 11 the spontaneous emission rates γ (g) e , γ (r) e , and γ e from the hfs level 5P 3/2 F ′ = 0 into guided modes, radiation modes, and both types of modes as functions of the fiber radius. As already noted in the previous subsection, the decay rate for this hfs level is equal to the average decay rate of an ensemble of twolevel emitters with dipoles oriented randomly in space. Figure 11 shows clearly that γ The dipole matrix element d eg is a vector whose spherical tensor components are specified by Eq. (1) in the quantization coordinate system {x Q , y Q , z Q }. It is clear that d eg depends on the orientation of the quantization axis z Q and so do the scalar product d eg ·e (α) and, hence, the spontaneous emission rate for the transition between the sublevels |e and |g . In the previous two subsections, we have studied the case where the quantization axis z Q coincides with the fiber axis z. In this subsection, we examine the dependencies of the rates on the orientation of the quantization axis. For certainty, we assume that the atom is positioned on the axis x.
We plot in Figs. 12 and 13 the dependencies of the fractional rates η e on the radial distance and the fiber radius for the quantization axis z Q = x and y. Com We plot in Figs. 14 and 15 the fractional rates η e as functions of the azimuthal angle ϕ Q and the zenithal angle θ Q of the quantization axis z Q . The figures show that the rates for the magnetic sublevels |F ′ = 3, M ′ = ±2 depend on the orientation of the quantization axis. It is interesting to note that the rate η e for the sublevels |F ′ = 3, M ′ = ±2 (see the green curves) does not depend on ϕ Q and θ Q . This independence is a consequence of the 1/2/0 ratio of the oscillatory strengths of the π/σ ± /σ ∓ transitions from the magnetic sublevels |F ′ = 3, M ′ = ±2 [72]. The symmetry properties of the profile functions with respect to opposite propagation directions and opposite phase circulation directions also play an important role.

D. Directional spontaneous emission rates
It has been shown in Sec. III that, when the quantization axis z Q coincides with the fiber axis z or, more generally, lies in the meridional plane containing the position of the atom, the spontaneous emission rates γ   coincides with the axis y. This axis is perpendicular to the meridional plane containing the position of the atom.
The blue curve in Fig. 20(a), which corresponds to M ′ = 3, N = TM 01 , and θ Q = π/2, indicates that the absolute value of the asymmetry factor ζ (N ) e is equal to 1 at four values ϕ Q = ϕ 0 , π − ϕ 0 , π + ϕ 0 , or 2π − ϕ 0 , where ϕ 0 ≃ 0.108π ≃ 19 • . This means that the spontaneous emission from the outermost sublevel |F ′ = 3, M ′ = 3 of the hfs level 5P 3/2 F ′ = 3 into the TM modes is unidirectional when the quantization axis z Q lies at an appropriate azimuthal angle ϕ Q in the fiber transverse plane xy. This interesting feature arises as a consequence of the properties of the cyclic transition and the TM modes. Indeed, the only allowed electric dipole transition from the sublevel |F ′ = 3, M ′ = 3 of the excited state 5P 3/2 is the σ + transition to the sublevel |F = 2, M = 2 of the ground state 5S 1/2 . The dipole of this transition is coupled to the counterclockwise circular component of the projection of the electric part of the field onto the plane x Q y Q , which is perpendicular to the quantization axis z Q . When the quantization axis lies in the fiber transverse plane xy and is oriented at an azimuthal angle ϕ Q = ϕ 0 , π − ϕ 0 , π + ϕ 0 , or 2π − ϕ 0 , where ϕ 0 = arcsin(|e z |/|e r |), the polarization of the projection of the electric part of a TM mode onto the plane x Q y Q is exactly circular at the position of the atom. The rotation direction of this polarization depends on the propagation direction f . Consequently, spontaneous emission from the sublevel |F ′ = 3, M ′ = 3 into the TM modes is unidirectional.

V. SUMMARY
In this work, we have studied spontaneous emission from a rubidium-87 atom into the fundamental and higher-order modes of a vacuum-clad ultrathin optical fiber. We have shown that the spontaneous emission rate depends on the magnetic sublevel, the type of modes, the orientation of the quantization axis, and the fiber radius. We have found that the rate of spontaneous emission into the TE modes is always symmetric with respect to the propagation directions. Meanwhile, the rates of spontaneous emission into other guided modes do not depend on the propagation direction when the quantization axis lies in the meridional plane containing the position of the atom. Asymmetry of spontaneous emission with respect to the propagation directions may appear when the output modes are not TE modes and the quantization axis does not lie in the meridional plane containing the position of the atom. We have shown that the rate of spontaneous emission into guided modes propagating in a given direction does not change when both the propagation direction and the magnetic quantum number are reversed. This result means that asymmetry of spontaneous emission with respect to the propagation directions leads to asymmetry with respect to the magnetic quantum numbers and vice versa. For the fiber radius in the range from 330 nm to 450 nm, the spontaneous emission into the HE 21 modes is stronger than into the HE 11 , TE 01 , and TM 01 modes. When the quantization axis coincides with the fiber axis and the radial distance is not too large, the rates of spontaneous emission from the outermost magnetic sublevels into guided modes are larger than those from the other sublevels. At the cutoff for higher-order modes, the rates of spontaneous emission into guided and radiation modes undergo steep variations, which are caused by the changes of the mode structure. Due to the mutual compensation of these changes, the variations of the total rate of spontaneous emission into both types of modes are smooth. The total fractional rate of emission into guided modes is most substantial when the fiber radius is around 180 nm, where the fiber supports only the fundamental HE 11 modes, or 340 nm, where the fiber supports not only the HE 11 modes but also the TE 01 , TM 01 , and HE 21 modes. We have shown that the spontaneous emission from the upper level up the cyclic transition into the TM modes is unidirectional when the quantization axis lies at an appropriate azimuthal angle in the fiber transverse plane. Our results lay the foundations for future research on manipulating and controlling the coupling of atoms, molecules, and dielectric particles to higher-order modes of ultrathin optical fibers.

ACKNOWLEDGMENTS
We acknowledge support for this work from the Okinawa Institute of Science and Technology Graduate University. S.N.C. and T.B. are grateful to JSPS for partial support from a Grant-in-Aid for Scientific Research (Grant No. 26400422).

Appendix A: Guided modes of a step-index fiber
Consider the model of a step-index fiber that is a dielectric cylinder of radius a and refractive index n 1 and is surrounded by an infinite background medium of refractive index n 2 , where n 2 < n 1 . We use the Cartesian coordinates {x, y, z}, where z is the coordinate along the fiber axis. We also use the cylindrical coordinates {r, ϕ, z}, where r and ϕ are the polar coordinates in the fiber transverse plane xy.
For a guided light field of frequency ω (free-space wavelength λ = 2πc/ω and free-space wave number k = ω/c), the propagation constant β is determined by the fiber eigenvalue equation [75] J ′ l (ha) haJ l (ha) Here, we have introduced the parameters h = (n 2 1 k 2 − β 2 ) 1/2 and q = (β 2 − n 2 2 k 2 ) 1/2 , which characterize the scales of the spatial variations of the field inside and outside the fiber, respectively. The integer index l = 0, 1, 2, . . . is the azimuthal mode order, which determines the helical phasefront and the associated phase gradient in the fiber transverse plane. The notations J l and K l stand for the Bessel functions of the first kind and the modified Bessel functions of the second kind, respectively.
The notations J ′ l (x) and K ′ l (x) stand for the derivatives of J l (x) and K l (x) with respect to the argument x.
For l ≥ 1, the eigenvalue equation (A1) leads to hybrid HE and EH modes [75]. The eigenvalue equation is given, for HE modes, as and, for EH modes, as Here, we have introduced the notation We label HE and EH modes as HE lm and EH lm , respectively, where l = 1, 2, . . . and m = 1, 2, . . . are the azimuthal and radial mode orders, respectively. Here, the radial mode order m implies that the HE lm or EH lm mode is the mth solution to the corresponding eigenvalue equation (A2) or (A3), respectively.
For l = 0, the eigenvalue equation (A1) leads to TE and TM modes [75]. The eigenvalue equation is given, for TE modes, as and, for TM modes, as We label TE and TM modes as TE 0m and TM 0m , respectively, where m = 1, 2, . . . is the radial mode order. The subscript 0 implies that the azimuthal mode order of TE and TM modes is l = 0. The radial mode order m implies that the TE 0m or TM 0m mode is the mth solution to the corresponding eigenvalue equation (A5) or (A6), respectively. According to [75], the fiber size parameter V is defined as V = ka n 2 1 − n 2 2 . The cutoff values V c for HE 1m modes are determined as solutions to the equation The electric component of the field can be presented in the form where E is the envelope. For a guided mode with a propagation constant β and an azimuthal mode order l, we can write where e is the mode profile function. In Eq. (A8), the parameters β and l can take not only positive but also negative values. We decompose the vectorial function e into the radial, azimuthal and axial components denoted by the subscripts r, ϕ and z, respectively. We summarize the expressions for the mode functions of hybrid modes, TE modes, and TM modes in the below [75].

Hybrid modes
We consider hybrid modes N = HE lm or EH lm . It is convenient to introduce the parameter s = l 1 h 2 a 2 + 1 q 2 a 2 J ′ l (ha) haJ l (ha) + K ′ l (qa) qaK l (qa) Here, the parameter A is a constant that can be determined from the propagating power of the field. Without loss of generality, we take A to be a real number.
In the cylindrical coordinates, the mode profile function of the electric component of a quasicircularly polarized hybrid mode N with a propagation direction f = ± and a phase circulation direction p = ± is given by e (ωN f p) =re r + pφe ϕ + fẑe z , where the mode function components e r , e ϕ , and e z are given by Eqs. (A10) and (A11) for β > 0 and l > 0. These components depend explicitly on the azimuthal mode order l and implicitly on the radial mode order m.
An important property of the mode functions of hybrid modes is that the longitudinal component e z is nonvanishing and in quadrature (π/2 out of phase) with the radial component e r . In addition, the azimuthal component e ϕ is also nonvanishing and in quadrature with the radial component e r . We note that the full mode function of the quasicircularly polarized hybrid mode is E (ωN f p) = e (ωN f p) e if βz+iplϕ , where β > 0 and l > 0. We have the following symmetry relations: e

TE modes
We consider transverse electric modes N = TE 0m . For r < a, we have Without loss of generality, we take A to be a real number. The mode profile function of the electric component of a TE 0m mode with a propagation direction f = ± can be written as which yield e (ωTE0mf ) = −e (ωTE0mf ) * .
Without loss of generality, we take A to be a real number. The mode profile function of the electric component of a TM mode with a propagation direction f = ± can be written as (B14)