Higher-order modes of vacuum-clad ultrathin optical fibers

We present a systematic treatment of higher-order modes of vacuum-clad ultrathin optical fibers. We show that, for a given fiber, the higher-order modes have larger penetration lengths, larger effective mode radii, and larger fractional powers outside the fiber than the fundamental mode. We calculate, both analytically and numerically, the Poynting vector, propagating power, energy, angular momentum, and helicity (or chirality) of the guided light. The axial and azimuthal components of the Poynting vector can be negative with respect to the direction of propagation and the direction of phase circulation, respectively, depending on the position, the mode type, and the fiber parameters. The orbital and spin parts of the Poynting vector may also have opposite signs in some regions of space. We show that the angular momentum per photon decreases with increasing fiber radius and increases with increasing azimuthal mode order. The orbital part of angular momentum of guided light depends not only on the phase gradient but also on the field polarization, and is positive with respect to the direction of the phase circulation axis. Meanwhile, depending on the mode type, the spin and surface parts of angular momentum and the helicity of the field can be negative with respect to the direction of the phase circulation axis.


I. INTRODUCTION
Near-field optics using optical fibers is currently a highly active and productive area of research that has implications for optical communication, sensing, computing, and even quantum information. Its main tool are so-called nanofibers, which are optical fibers that are tapered to a diameter comparable to or smaller than the wavelength of light [1][2][3]. The essence of the tapering technique is to heat and pull a single-mode optical fiber to a very small thickness, while maintaining the taper condition adiabatically [1][2][3][4]. 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. Thus, these fibers can be treated as very thin vacuum-clad silica-core fibers.
In a vacuum-clad nanofiber, the guided field penetrates an appreciable distance into the surrounding medium and appears as an evanescent wave carrying a significant fraction of the power and having a complex polarization pattern [5][6][7]. These fibers offer high transmission and strong confinement of guided light in the transverse plane of the fiber. This confinement allows one to efficiently couple guided light to emitters placed on or near the fiber surface. Such fibers are therefore versatile tools for coupling light and matter and have a wide range of potential practical applications [8,9]. For example, they have been used for trapping atoms [10][11][12], for probing atoms [13][14][15][16][17][18][19], molecules [20], quantum dots [21], and color centers in nanodiamonds [22,23], and for mechanical manipulation of small particles [24][25][26]. Due to the lack of cutoff as well as the possession of a small mode area and a simple mode structure, the fundamental HE 11 mode has been exploited in most studies to date.
However, tapered fibers can also be fabricated with slightly larger diameters and/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 polarization distributions. In addition, the higher-order modes can have larger angular momentum compared to the HE 11 mode. For ease of reference, the micro-and nanofibers that can support the fundamental mode and several higher-order modes are called ultrathin fibers in this paper.
Despite increased interest in higher-order modes of ultrathin fibers, systematic treatments for the basic prop-erties of light fields in such modes do not exist. Although the full and exact fiber theory [42] is applicable also to ultrathin fibers, deep understanding can only be reached only by combining a systematic and comprehensive analysis with detailed numerical calculations for fibers with parameters in the range of experimental interest. The purpose of this work is to present such a systematic treatment. We show that, for a given fiber, the higher-order modes have larger penetration lengths, larger effective mode radii, and larger fractional powers outside the fiber than the fundamental mode. We calculate analytically and numerically the Poynting vector, propagating power, energy, angular momentum, and helicity (or chirality) of guided light.
The paper is organized as follows. In Sec. II we review the theory of guided modes of optical fibers and present the results of numerical calculations for the propagation constants and penetration lengths of the modes of fibers with the parameters in the range of experimental interest. Section III is devoted to the study of the electric intensity distribution and the effective mode radius. In Sec. IV we calculate the Poynting vector, propagating power, and energy per unit length, and examine the orbital and spin parts of the Poynting vector. Section V is devoted to the study of angular momentum of guided light and its orbital, spin, and surface parts. In Sec. VI we calculate the helicity and the associated chirality of guided light. Our conclusions are given in Sec. VII.

II. GUIDED MODES OF OPTICAL FIBERS
In this section, we first briefly review the theory of guided modes of optical fibers and then calculate the propagation constants and evanescent-wave penetration lengths of the fundamental mode and higher-order modes of ultrathin fibers with parameters in the range of experimental interest.
For this we consider the model of a step-index fiber that is a dielectric cylinder of radius a and refractive index n 1 , surrounded by an infinite background medium of refractive index n 2 , where n 2 < n 1 . 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.
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 [42] J ′ l (ha) haJ l (ha) + K ′ l (qa) qaK l (qa) n 2 1 J ′ l (ha) haJ l (ha) + n 2 2 K ′ l (qa) qaK l (qa) = l 2 1 h 2 a 2 + 1 q 2 a 2 2 β 2 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 (1) leads to hybrid HE and EH modes [42], for which the eigenvalue equations are given by Eqs. (A1) and (A2) in Appendix A. We label these modes as HE lm and EH lm , where l = 1, 2, . . . is the azimuthal and m = 1, 2, . . . the radial mode orders. The radial mode order m implies that the HE lm or EH lm mode is the m-th solution to the corresponding eigenvalue equation.
For l = 0, the eigenvalue equation (1) leads to TE and TM modes [42], for which the eigenvalue equations are given by Eqs. (A4) and (A5) in Appendix A. We label these modes as TE 0m and TM 0m , where again m = 1, 2, . . . is the radial mode order and the subscript 0 implies that the azimuthal mode order of each mode is l = 0. We are interested in vacuum-clad ultrathin fibers, which can support not only the fundamental HE 11 mode but also several higher-order modes in the optical region. For this we plot in Fig. 1 the propagation constant β for the HE 11 mode and several higher-order modes as a function of the fiber radius a for a wavelength of light that is chosen to be λ = 780 nm. The fiber is assumed to be made of silica, with a refractive index n 1 = 1.4537, and the surrounding medium is air or vacuum, with a refractive index n 2 = 1. One can see that the first two higherorder modes, TE 01 and TM 01 , appear when a ≃ 283 nm, and the next higher-order mode, HE 21 , appears when a ≃ 325 nm. It is clear that the number of modes supported by the fiber increases with increasing fiber radius a. The numerical results presented in Fig. 1  Outside the fiber, the guided modes are evanescent waves in the radial direction r. The penetration depth is characterized by the parameter Λ = 1/q, which we show for the HE 11 mode and several higher-order modes as a function of the fiber radius a in Fig. 2. One can see that, near to the cutoffs, the penetration length Λ is large, that is, the field is not tightly confined inside the fiber. Furthermore, when the fiber radius a increases, the penetration length decreases to the limiting value Λ min = 1/(k n 2 1 − n 2 2 ). In general, for a given fiber, the higher-order modes have larger penetration lengths than the HE 11 mode.
We will next discuss the mode functions [42]. For this we will write the electric and magnetic components of the field in the form E = (Ee −iωt + c.c.)/2 and H = (He −iωt + c.c.)/2, where E and H are spatial envelope functions, which obey the Helmholtz equation. They are the mode functions we are interested in and, for a guided mode with a propagation constant β and an azimuthal mode order l, we can write them as E = ee iβz+ilϕ and H = he iβz+ilϕ . Here, e and h are the reduced mode profile functions of the electric and magnetic components of the field, respectively, and β and l can take not only positive but also negative values. In the following we will consider hybrid modes, TE modes, and TM modes separately.
A. Hybrid modes

Quasicircularly polarized hybrid modes
In this section, we consider quasicircularly polarized hybrid HE and EH modes. For convenience, we use the notations β > 0 and l > 0 for the propagation constant and the azimuthal mode order, respectively. We introduce the index f = +1 or −1 (or simply f = + or −) for the positive (+ẑ) or negative (−ẑ) propagation direction, which leads to the corresponding propagation phase factor of e iβz or e −iβz . We also introduce the index p = +1 or −1 (or simply p = + or −) for the counterclockwise or clockwise phase circulation, corresponding to the azimuthal phase factor of e ilϕ or e −ilϕ . The index p = + or − also indicates that the central phase circulation axis is +ẑ or −ẑ. We can label quasicircularly polarized hybrid modes by the mode index µ = (f lp), which can also be extended to include the mode type, HE or EH, and the radial mode order, m, when necessary.
We choose a notation in which we decompose an arbitrary vector V =rV r +φV ϕ +ẑV z into the radial, azimuthal, and axial components denoted by the subscripts r, ϕ, and z. The notationsr =x cos ϕ +ŷ sin ϕ, ϕ = −x sin ϕ+ŷ cos ϕ, andẑ stand for the unit basis vectors of the cylindrical coordinate system {r, ϕ, z}, withx andŷ being the unit basis vectors of the Cartesian coordinate system for the fiber transverse plane xy. The position vector in the fiber transverse plane is given by r = rr = xx + yŷ.
In the cylindrical coordinates, the reduced mode profile functions e (f lp) (r) and h (f lp) (r) of the electric and magnetic components of a quasicircularly polarized hybrid mode with the propagation direction f , the azimuthal mode order l, and the phase circulation direction p are then given by where the electric mode function components e r , e ϕ , and e z and the magnetic mode function components h r , h ϕ , and h z are given by Eqs. (A9)-(A12) for β > 0 and l > 0 in Appendix A. These mode function components depend explicitly on the azimuthal mode order l and are implicitly dependent on the radial mode order m. An important property of the mode functions is that the longitudinal components e z and h z are nonvanishing and in quadrature (π/2 out of phase) with the radial components e r and h r , respectively. In addition, the azimuthal components e ϕ and h ϕ are also nonvanishing and in quadrature with the radial components e r and h r , respectively. The electric and magnetic polarizations of hybrid modes are not of the TE and TM types. Note that the full mode functions for quasicircularly polarized hybrid modes are given by (3)

Quasilinearly polarized hybrid modes
Quasilinearly polarized hybrid modes are linear superpositions of counterclockwise and clockwise quasicircularly polarized hybrid modes. The full mode functions of the electric and magnetic components of the guided field in a quasilinearly polarized hybrid mode (f, l, ϕ pol ) are given by [42] , Here, the phase angle ϕ pol determines the orientation of the symmetry axes of the mode profile in the fiber transverse plane. In particular, the specific phase angle values ϕ pol = 0 and π/2 define two orthogonal polarization profiles, one being symmetric with respect to the x axis and the other being the result of the rotation of the first one by an angle of π/2l in the fiber transverse plane xy. We can write where e (f lϕ pol ) and h (f lϕ pol ) are the reduced mode profile functions of quasilinearly polarized hybrid modes and are given as , Inserting Eqs. (2) into Eqs. (6) yields In particular, we find, for ϕ pol = 0, e (f l,0) = √ 2(re r cos lϕ + iφe ϕ sin lϕ + fẑe z cos lϕ), and, for ϕ pol = π/2, e (f l,π/2) = √ 2(re r sin lϕ − iφe ϕ cos lϕ + fẑe z sin lϕ), h (f l,π/2) = √ 2(−ifrh r cos lϕ + fφh ϕ sin lϕ − iẑh z cos lϕ).

B. TE modes
We again label the propagation directions of TE modes by the index f = + or −. The reduced mode profile functions of the electric and magnetic components of TE modes with the propagation directions f can be written as where the mode function components e ϕ , h r , and h z are given by Eqs. (A13)-(A16) for β > 0 in Appendix A. They depend implicitly on the radial mode order m. It is clear from Eqs. (10) that, for TE modes, we have e The electric polarization of a TE mode is therefore linear and aligned along the azimuthal direction. Meanwhile, since h r is π/2 out of phase with respect to h z , the magnetic polarization of the mode is elliptical in the meridional rz plane, which contains the radial r axis and the fiber z axis. The full mode functions of TE modes are given by

C. TM modes
We also label the propagation directions of TM modes by the index f = + or −. The reduced mode profile functions of the electric and magnetic components of TM modes with the propagation directions f can be written as where the mode function components e r , e z , and h ϕ are given by Eqs. (A17)-(A20) for β > 0 in Appendix A. They depend implicitly on the radial mode order m. It is clear from Eqs. (11) that, for TM modes, we have e The magnetic polarization of a TM mode is therefore linear and aligned along the azimuthal direction. Meanwhile, since e r is π/2 out of phase with respect to e z , the electric polarization of the mode is elliptical in the meridional rz plane. The full mode functions of TM modes are given by

III. SPATIAL INTENSITY DISTRIBUTIONS
In this section, we study the electric intensity distributions |e| 2 = |e r | 2 + |e ϕ | 2 + |e z | 2 of the fields in the fundamental HE 11 and several higher-order modes, namely the TE 01 , TM 01 , HE 21 , and EH 11 modes. In the cases of the hybrid HE 11 , HE 21 , and EH 11 modes, we examine both quasicircular and quasilinear polarizations. The inner (r/a < 1) and outer (r/a > 1) parts are distinguished by the blue and cyan colors, respectively. In all four cases The distributions are normalized to the same power. The fiber radius is chosen to be a = 400 nm. All other parameters are as for Fig. 1.
The cross-sectional profiles of the electric intensity distributions |e| 2 of the fields in the quasicircularly polarized HE 11 mode, the TE 01 mode, the TM 01 mode, and the quasicircularly polarized HE 21 mode are shown in Fig. 3. One can note that all of them are azimuthally symmetric. To show the spatial dependencies of these distributions more clearly, we display in Fig. 4 cuts in the radial direction. We note that the group of the TE 01 , TM 01 , and HE 21 modes corresponds to the first higher-order LP 11 mode of weakly guiding fibers [42]. Meanwhile, the fundamental HE 11 mode corresponds to the lowest LP 01 mode of weakly guiding fibers [42].
For hybrid modes, we can use quasilinear polarization instead of quasicircular polarization. In order to illustrate fields in hybrid modes with quasilinear polariza-Radial distance r (nm) Electric intensity (arb. units)  tion, we display in Fig. 5 the cross-sectional profiles of the electric intensity distributions |e| 2 of the fields in the quasilinearly polarized HE 11 mode and the quasilinearly polarized HE 21 mode. To show the spatial dependencies more clearly, radial cuts for both modes in two different directions are shown in Fig. 6.
In addition to the TE, TM, and HE modes, there is an-Radial distance r (nm) Electric intensity (arb. units) other type of guided modes, namely the EH modes. The cross-sectional profiles of the electric intensity distributions |e| 2 of the fields in the quasicircularly and quasilinearly polarized EH 11 modes are shown in Fig. 7. Similarly, Fig. 8 depicts cuts in the radial direction. Inside the fiber, the field intensity is not a fast reducing function of the radial distance r. A discontinuity of the field intensity is observed at the position of the fiber surface. This discontinuity is due to the boundary condition for the normal (radial) component e r of the electric field. Since the difference between the refractive indices of the silica core and the vacuum cladding is large, the discontinuity of the field at the position of the fiber surface is dramatic.
Outside the fiber, the field intensity monotonically and quickly reduces with increasing radial distance r. This behavior is a consequence of the evanescent-wave nature of guided fields, which do not propagate along the radial direction. Comparison between the figures shows that, for the parameters used, the fraction of the field intensity distribution outside the fiber for higher-order modes is larger than that for the HE 11 mode.
We observe from Figs. 3 and 7(a) that, for quasicircularly polarized hybrid HE and EH modes, TE modes, and TM modes, the spatial distribution of the field intensity is cylindrically symmetric. In these cases, the outer parts of the electric intensity distributions of the different modes look very similar to each other as they exhibit the evanescent wave behavior. Meanwhile, the inner parts of the electric intensity distributions of different modes look very different from each other. Indeed, the inner parts of the electric intensity profiles have the shape of a cone in Figs. 3(a) and 3(c), the shape of a doughnut in Figs. 3(b) and 3(d), and the shape of a combination of a cone and a doughnut in Fig. 7(a).
We observe from Figs. 5 and 7(b) that, for quasilinearly polarized hybrid modes, the spatial distribution of the field intensity is not cylindrically symmetric. In the inner and outer vicinities of the fiber surface, the field intensity strongly varies with varying azimuthal angle.
Finally, from Figs. 4(b), 4(d), and 6(b) one can see that, in the cases of the TE 01 and HE 21 modes, the electric field intensity is exactly equal to zero at the center of the fiber. Figure 8(b) shows that, for the quasilinearly polarized EH 11 mode, the electric field intensity is exactly equal to zero at two centrally symmetric off-center positions along the y axis inside the fiber.
The spatial profiles of the fields presented in Figs The effective mode area can be defined as A eff = ( |e| 2 dr) 2 / |e| 4 dr, where we use the notation dr = 2π 0 dϕ ∞ 0 r dr. This allows us to define an effective mode radius as r eff = A eff /π. The parameters A eff and r eff characterize the confinement of the field mode in the fiber transverse plane. We show in Fig. 9 the effective mode radius r eff as a function of the fiber radius a for the fundamental mode and several higher-order modes. It is clear that the effective radii of the higher-order modes are larger than that of the fundamental mode. In addition, different modes have different minimum effective radii. These minimum values are achieved at different points corresponding to different values of a. For the light wavelength λ = 780 nm used in our numerical calculations, the smallest value of r eff is about 353 nm and is achieved for the fundamental HE 11 mode of a fiber with the radius a = 275 nm.

IV. POYNTING VECTOR, POWER, AND ENERGY PER UNIT LENGTH
Next, we calculate the Poynting vector, propagating power, and energy per unit length. We show that the axial and azimuthal components of the Poynting vector can be negative with respect to the direction of propagation and the direction of phase circulation, respectively. In order to get deeper insight into the connection between linear and angular momenta of light, we also study the decomposition of the Poynting vector into the orbital and spin parts. We show that the orbital and spin parts of the Poynting vector can have opposite signs in some regions of space.

A. Poynting vector
An important characteristic of light propagation is the cycle-averaged Poynting vector We introduce the notations S z , S ϕ , and S r for the axial, azimuthal, and radial components of the vector S in the cylindrical coordinates. For guided modes of fibers, we have S r = 0 and The explicit expressions for S z and S ϕ are given by Eqs. (B1)-(B8) in Appendix B. We note that the existence of a nonzero azimuthal component S ϕ of the Poynting vector for guided fields leads to a force transverse to the direction of propagation. This is similar to the situation for light beams with a transverse phase gradient, for which transverse optical forces have been experimentally observed [43]. It is worth nothing that, due to the interference between different terms associated with different Bessel functions, the sign of the azimuthal component S ϕ of the Poynting vector of a quasicircularly polarized hybrid mode and the sign of the axial component S z of the Poynting vector of a quasilinearly polarized hybrid mode can vary in space. The details are given in Appendix B.
For the quasicircularly polarized HE 11 mode, the TE 01 mode, the TM 01 mode, and the quasicircularly polarized HE 21  The axial component S z and the azimuthal component S ϕ of the Poynting vector for the quasicircularly polarized HE 12 and EH 11 modes are displayed in Fig. 11. The dashed blue curves of the figure shows that S ϕ is negative in a localized region of space. For the HE 12 mode [ Fig. 11(a)], this region is inside the fiber. However, for the EH 11 mode [ Fig. 11(b)], part of this region is the outside and part is in the inside of the fiber.
Figures 10(a), 10(d), 11(a), and 11(b) and additional numerical calculations, which are not shown here, confirm that, outside the fiber, the azimuthal component S ϕ of the Poynting vector is positive for quasicircularly polarized HE modes but negative for quasicircularly polarized EH modes.
It is not surprising that a component of the Poynting vector can have different signs in different regions of space [44,45]. Similar results have been obtained for the axial component of the Poynting vector of a guided mode [44] and for the axial and azimuthal components of the Poynting vector of a Bessel beam [45]. In fact, we have confirmed that the axial component S z of the Poynting vector of the quasilinearly polarized HE 11 mode can become negative when the refractive index n 1 of the fiber is large enough (n 1 /n 2 > 2.71 for the HE 11 mode) [44]. We show in Fig. 12 a similar result for the quasilinearly polarized higher-order HE 21 mode. One can see that in this case the Poynting vector is negative in four regions around the fiber surface at azimuthal angles around the Radial distance r (nm) Poynting vector (arb. units) values ϕ = 0, π/2, π, and 3π/2.

B. Propagating power
The optical power carried by the fiber is given by It can be split as P = P in + P out , where P in and P out are the propagating powers inside and outside the fiber and explicit expressions for both are given by Eqs. (C1)-(C6) in Appendix C. The fractional power outside the fiber η P is defined as η P = P out /P . We display η P in Fig. 13 as a function of the fiber radius a for the HE 11 mode and several higherorder modes. We observe from the figure that η P reduces with increasing a and that the fractional powers outside the fiber for higher-order modes are larger than that for the fundamental mode. It is interesting to note that, near the cutoffs for the EH 11 and HE 31 modes, the factor η P is significantly smaller than unity, unlike the cases of the HE 11 and HE 12 modes. We show in Appendix C that, for the EH lm modes with l = 1, 2, . . . and the HE lm modes with l = 3, 4, . . . , the limiting values of the factor η P in the cutoff regions are smaller than unity, in agreement with the aforementioned numerical results. We also show in Appendix C that, for the TE 0m and TM 0m modes and the HE lm modes with l = 1 or 2, the limiting values of the factor η P in the cutoff regions are equal to unity. Despite this prediction, we observe from Fig. 13   TM 01 , and HE 21 modes are slightly smaller than unity. These numerical deviations are due to the steep slopes of the curves that make it difficult to approach the cutoffs.
Note that the numerical results presented in Fig. 13 are in agreement with the results presented in Ref. [46].

C. Energy per unit length
The cycle-averaged energy per unit length is given by where n(r) = n 1 for r < a and n 2 for r > a. The first and second terms on the right-hand side of expression (15) correspond to the electric and magnetic parts, respectively, of the energy of the field. For guided modes, these parts are equal to each other. We can split U as U = U in + U out , where U in and U out are the energies per unit length inside and outside the fiber, and their explicit expressions are given by Eqs. (D1)-(D6) in Appendix D. The fractional energy outside the fiber η U = U out /U is shown as a function of the fiber radius a for the HE 11 mode and several higher-order modes in Fig. 14. One can see that the behavior of η U is very similar, but not identical, to that of η P .

D. Orbital and spin parts of the Poynting vector
It is known that the Poynting vector of the field can be decomposed into two parts, the orbital part and the spin part [47][48][49][50][51][52][53]. In the dual-symmetric formalism, the decomposition takes the form [47][48][49][50][51][52][53] where S orb = S e-orb + S h-orb is the orbital part, with its electric and magnetic components and S spin = S e-spin + S h-spin is the spin part, with its electric and magnetic components In Eq. (17), the dot product applies to the field vectors, that is, A · (∇)B ≡ i=x,y,z A i ∇B i for arbitrary field vectors A and B.
In general, we have the equality S e = S h , where S e = S e-orb + S e-spin and S h = S h-orb + S h-spin are the electric and magnetic components of the Poynting vector. However, we may observe the inequalities S e-orb = S h-orb and S e-spin = S h-spin . The explicit expressions for the electric and magnetic components of the orbital and spin parts of the Poynting vector of guided light are given in Appendix E. It is worth noting that the orbital part S orb of the Poynting vector is proportional to the canonical momentum of light, which determines the radiation pressure force upon a small dipole Rayleigh particle [49][50][51][52][53].
We show in Appendix E that the orbital parts S orb z and S orb ϕ of the axial and azimuthal components, respectively, of the Poynting vector are positive with respect to the direction of propagation and the direction of phase circulation, respectively. Meanwhile, the signs of the spin parts, S spin z and S spin ϕ , of the Poynting vector can vary in the fiber transverse plane and can hence be negative with respect to the direction of propagation and the direction of phase circulation, respectively, in some regions of space. Thus, the orbital and spin parts of the Poynting vector can have opposite signs in certain regions of space. We show numerical results confirming this in Figs. 15-17. We show in Appendix E that the orbital part S orb z of the axial component S z of the Poynting vector is determined by the local density of energy. Meanwhile, the orbital part S orb ϕ of the azimuthal component S ϕ of the Poynting vector of a quasicircularly polarized hybrid mode depends on not only the local phase gradient but also the local polarization, unlike the case of uniformly polarized paraxial beams [53].
The radial dependencies of the orbital part S orb z and the spin part S spin z of the axial component S z of the Poynting vector for the quasicircularly polarized HE 11 mode, the TE 01 mode, the TM 01 mode, and the quasicircularly polarized HE 21 mode are shown in Fig. 15. The radial dependencies of the orbital part S orb ϕ and the spin part S spin ϕ of the azimuthal component S ϕ of the Poynting vector for the quasicircularly polarized HE 11 and HE 21 modes are shown in Fig. 16. Additionally, the radial dependencies of the orbital and spin parts of the axial and azimuthal components of the Poynting vector for the quasicircularly polarized EH 11 mode are displayed in Fig. 17.
These figures show that the orbital parts S orb  Figures 16 and 17(b) show that, outside the fiber, the spin part S spin ϕ of the azimuthal component S ϕ of the Poynting vector is positive for HE modes and is negative for EH modes. These features are also observed for S ϕ (see Figs. 10 and 11).

V. ANGULAR MOMENTUM OF GUIDED LIGHT
In this section, we calculate the angular momentum of guided light and also study its orbital, spin, and surface parts. We show that the orbital part of angular momentum depends not only on the phase gradient, but also on the field polarization, and is always positive with respect to the direction of the phase circulation axis. Meanwhile, the spin and surface parts of angular momentum and the helicity (chirality) of light can be negative with respect to the direction of the phase circulation axis. We find that the signs of the spin and surface parts of the transverse angular momentum density of the fundamental and higher-order modes depend on the direction of propagation.

A. Angular momentum of guided light
For the electromagnetic field in free space, the linear momentum density is given by p local = S/c 2 [54]. For the field in a dielectric medium, several formulations for the linear momentum density can be found in the literature [55]. The Abraham formulation [56] takes p local = [E × H]/c 2 , which is sometimes interpreted as the field-only contribution to the momentum of light. On the other hand, the Minkowski formulation [57] takes p local = [D × B]. While the appropriate form remains contentious because the debate has not been settled by experiments, the Abraham formulation is generally accepted [54,58]. Therefore, in our basic calculations, we adopt the Abraham formulation for the field linear momentum density inside and outside the fiber.
With the above definition of the linear momentum density, the angular momentum density of the electromagnetic field is given by j local ≡ (R × p local ) = (R × S)/c 2 . Here, R = xx + yŷ + zẑ is the position vector in the three-dimensional space. Integrating j local over the crosssectional plane of the fiber then yields the angular momentum per unit length Note that TE and TM modes and quasilinearly polarized HE and EH modes have no angular momentum. We consider the cycle-averaged angular momentum per unit length J of quasicircularly polarized HE and EH modes. The only nonzero component of J is aligned along the fiber axis and is given by Thus, the axial angular momentum per unit length J z is determined by the azimuthal component S ϕ of the Poynting vector. We can write J z = J in z + J out z , where J in z and J out z are the parts of the angular momentum of light inside and outside the fiber. The explicit analytical expressions for J in z and J out z are given by Eqs. (F1) and (F2) in Appendix F. According to these expressions, the axial angular momentum per unit length J z depends on the direction of phase circulation, specified by the index p, but does not depend on the direction of propagation, specified by the index f .
The angular momentum per photon j z =hωJ z /U for quasicircularly polarized hybrid modes with the positive (counterclockwise) phase circulation direction p = + is shown as a function of the fiber radius a in Fig. 18. One can see that j z decreases with increasing a and increases with increasing l. Comparison between HE and EH modes shows that, for a given set of l and m, the angular momentum per photon j z for an EH lm mode is smaller than that for the corresponding HE lm mode. This feature is related to the fact that, outside the fiber, the azimuthal component S ϕ of the Poynting vector is positive for HE modes and is negative for EH modes (see Fig. 11). Figure 18 also shows that the EH 11 mode has a (nm) the lowest angular momentum per photon. It is clear that the angular momentum per photon in a higher-order hybrid mode is large when the azimuthal mode order l is large.

B. Orbital and spin parts of angular momentum
The angular momentum per unit length J of a light beam can be decomposed into orbital, spin, and surface parts as J = J orb + J spin + J surf [59][60][61][62][63]. However, the identification of terms as orbital, spin, and surface components is not unique [64,65].
In the dual-symmetric formalism, the orbital and spin parts of angular momentum per unit length are given as [50][51][52][53]66] and In Eq. (21), the dot product applies to the field vectors, that is, A · (R × ∇)B ≡ i=x,y,z A i (R × ∇)B i for arbitrary field vectors A and B. Detailed discussions of various aspects of optical orbital angular momentum can be found in Ref. [67].
Meanwhile, the surface part of angular momentum per unit length is, the dual-symmetric formalism, given as [59][60][61][62][63] where we have used the notation The orbital part of angular momentum per unit length is related to the orbital part S orb of the Poynting vector via the formula J orb = (1/c 2 ) (R × S orb ) dr. The spin and surface parts of angular momentum per unit length are related to the spin part S spin of the Poynting vector via the formula J spin + J surf = (1/c 2 ) (R × S spin ) dr.
The surface part of angular momentum is usually omitted in the literature [63]. The reason is that, when the field vanishes sufficiently quickly in the limit of large distances, the surface part is, due to the Gaussian theorem, identical to zero. For example, for a bullet-like light wave packet with a finite transverse and longitudinal extent, the surface part of angular momentum can be neglected. However, for a pencil-like light beam whose span along the direction of propagation is virtually infinite, the surface part is not vanishing [63].
For TE and TM modes, we have J orb = J spin = J surf = 0. For quasicircularly polarized HE and EH modes, the only nonzero components of the vectors J orb , J spin , and J surf are the axial components and Here, we have introduced the notations a ± 0 = lim ε→0 (a ± ε). According to Eqs. An important point to note here is that the above results, derived for angular momentum of light in guided modes, are different from the results for angular momentum of light in scalar Laguerre-Gaussian beams [68,69]. The main reason is that a guided light beam is a vector beam [70], whose polarization is not uniform in the crosssectional plane. Another important reason is that the guided mode has two parts: one inside the fiber, where the medium is a dielectric, and the other one outside the fiber, where the medium is the vacuum. In addition, the discontinuity of the refractive index at the fiber surface leads to the appearance of the surface part of the angular momentum of light.
It follows from Eq. (24) that the orbital part J orb z of angular momentum of light in an arbitrary hybrid mode is positive or negative when the phase circulation direction index p is positive or negative, respectively. Note that p = + or− means that the phase circulation direction in the xy plane is counterclockwise or clockwise, that is, the phase circulation axis is +ẑ or −ẑ, respectively. Thus, the orbital part J orb z of angular momentum of light in a hybrid mode is positive with respect to the direction of the phase circulation axis. Meanwhile, the expression on the right-hand side of Eq. (24) contains not only the terms l|e| 2 and l|h| 2 , which result from the local phase gradient, but also the terms Im(e * r e ϕ ) and Im(h * r h ϕ ), which result from the local polarization. Thus, the orbital part J orb z of angular momentum depends not only on the phase gradient but also on the field polarization.
Equation (25) shows that the spin part J spin z of angular momentum is determined by the polarization of the field, whereas, according to Eq. (26), the surface part J surf z of angular momentum is associated with the discontinuity of the spin density at the fiber surface. It is clear that the discontinuity of the spin density is induced by the discontinuity of the refractive index of the medium at the fiber surface. Unlike the orbital part J orb z z of angular momentum per photon is always positive with respect to the direction of the phase circulation axis, which is in agreement with Eq. (24). Note that j orb tially smaller thanh in the cases of the HE 11 and HE 12 modes but is comparable to or larger thanh in the cases of higher-order HE and EH modes. From Fig. 20 one can see that the spin part j spin z of angular momentum per photon is positive with respect to the direction of the phase circulation axis for the HE modes and negative for the EH modes. Furthermore, for the HE lm modes with the azimuthal mode order l = 1, the spin part j spin z is dominant to the orbital part j orb z .
However, for the HE lm modes with l ≥ 2 and the EH lm modes, the orbital part j orb . The above result is in agreement with the results of Ref. [69], where it has been shown for Laguerre-Gaussian beams that the local spin density can be positive in some regions and negative in others.
Although the transverse component of angular momentum of guided light is zero, the local density of this component is not zero. Indeed, the local density of the azimuthal component of angular momentum is given by ρ Jϕ = −rS z /c 2 . It can be decomposed as The azimuthal orbital angular momentum density ρ J orb ϕ originates from the orbital part S orb z of the axial component of the Poynting vector. The azimuthal spin and surface angular momentum densities ρ J spin ϕ and ρ J surf ϕ result from the spin part S spin z of the axial component of the Poynting vector. Note that the first and second terms in the expressions on the right-hand side of Eqs. (27) correspond to the electric and magnetic parts, respectively. Equations (27) can be used not only for quasicircularly polarized hybrid modes but also for TE and TM modes.
According to Eq. (27), the signs of ρ J orb ϕ , ρ J spin ϕ , and ρ J surf ϕ depend on the direction of propagation f . The dependence of the local transverse spin density ρ J spin ϕ on the direction of propagation is a signature of spin-orbit coupling of light [51][52][53][72][73][74][75]. Note that both ρ J spin ϕ and ρ J surf ϕ appear as a result of the facts that the longitudinal field components e z and h z are nonvanishing and in quadrature with the radial field components e r and h r , respectively. It has been shown that, due to spin-orbit coupling of light, 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 [76][77][78][79][80][81].
The helicity of a light beam is closely related to its chirality. Indeed, according to [87], the cycle-averaged optical chirality density of a monochromatic light field can be characterized by the quantity [82,[87][88][89][90][91] ρ chir = n 2 2c Im[E · H * ], (29) so that ρ chir = n 2 ωρ hlcy . Thus, in the frequency domain, the chirality density is proportional to the helicity density and the proportionality factor is n 2 ω. Note that the optical chirality density (29) can be measured from the asymmetry in the rates of excitation between a small chiral molecule and its mirror image [87,90]. However, according to [92], there is no single measure of chirality. The optical chirality measure (29) is appropriate to chiral effects arising from interference between electric and magnetic dipole transitions [87,90], whereas for chiral effects in spontaneous emission and scattering from atoms with rotating electric dipoles, an appropriate measure of optical chirality is the ellipticity vector of the field polarization [76][77][78][79][80][81]93].
For quasicircularly polarized hybrid modes, we find the following expression for the helicity density: The optical helicity per unit length is It is clear from Eqs. (30) and (31) that when we reverse the propagation direction f or the phase circulation direction p, the sign of helicity per unit length is reversed. The explicit expression for the helicity per unit length J hlcy in terms of the fiber parameters is given by Eq. (G3) in Appendix G. Note that the optical helicity of TE and TM guided modes is zero. The helicity per photon j hlcy =hωJ hlcy /U is shown as a function of the fiber radius a in Fig. 22. One can see that, for the propagation direction f = + and the phase circulation direction p = +, the helicity per photon j hlcy is positive for the HE modes and negative for the EH modes. We note that the magnitude of the helicity per photon j hlcy in a guided mode does not exceed the valueh, which is the value of the helicity per photon of circularly polarized light in free space.

VII. SUMMARY
In this work, we have presented a systematic treatment of higher-order modes of vacuum-clad ultrathin optical fibers. We have shown that, for a given fiber, the higherorder modes have larger penetration lengths, larger effective mode radii, and larger fractional powers outside the fiber than the fundamental mode. We have calculated analytically and numerically the Poynting vector, propagating power, energy, angular momentum, and helicity of the field. In doing so we have shown that the axial component S z and the azimuthal component S ϕ of the Poynting vector can be negative with respect to the direction of propagation and the direction of phase circulation, respectively, depending on the position, the mode type, and the fiber parameters. The occurrence of such a negative axial or azimuthal component of the Poynting vector indicates the possibility of the occurrence of a negative force upon an atom or a small particle. We have also found that the orbital and spin parts of the Poynting vector may have opposite signs in some regions of space. We have shown that, for the EH lm modes with l = 1, 2, . . . and the HE lm modes with l = 3, 4, . . . , the limiting values of the fractional power outside the fiber in the cutoff regions are significantly smaller than unity. Meanwhile, for the TE 0m and TM 0m modes and the HE lm modes with l = 1 or 2, the limiting values of the fractional power outside the fiber in the cutoff regions are equal to unity. Our calculations have shown that the angular momentum per photon decreases with increasing fiber radius and increases with increasing azimuthal mode order, and the angular momentum per photon of an EH lm mode is smaller than that of the corresponding HE lm mode. We have found that the orbital part of angular momentum of guided light depends on not only the phase gradient but also the field polarization, and is positive with respect to the direction of the phase circulation axis. Meanwhile, the spin and surface parts of angular momentum and the helicity (chirality) of light in an EH mode are negative with respect to the direction of the phase circulation axis. We have shown that the signs of the spin and sur-face parts of the transverse angular momentum density of the fundamental and higher-order modes depend on the direction of propagation. The directional dependence of the local transverse spin and surface angular momentum densities is a signature of spin-orbit coupling of light and appears as a result of the facts that the longitudinal field components are nonvanishing and in quadrature with the radial field components. Our results lay the foundations to future research on manipulating and controlling the motion of atoms, molecules, and dielectric particles using higher-order modes of ultrathin 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
For l ≥ 1, the eigenvalue equation (1) leads to hybrid HE and EH modes [42]. For the HE modes, the respective eigenvalue equation is given as and, for the EH modes, as Here, we have introduced the notation For l = 0, the eigenvalue equation (1) leads to TE and TM modes [42], with the eigenvalue equation for the TE modes given as and for the TM modes as According to [42], 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 J 1 (V c ) = 0. For HE lm modes with l = 2, 3, . . . , the cutoff values are obtained as nonzero solutions to the equation (n 2 1 /n 2 2 + 1)(l − 1)J l−1 (V c ) = V c J l (V c ). The cutoff values V c for EH lm modes, where l = 1, 2, . . . , are determined as nonzero solutions to the equation J l (V c ) = 0. For TE 0m and TM 0m modes, the cutoff values V c are obtained as solutions to the equation J 0 (V c ) = 0.
The electric and magnetic components of the field can be presented in the form where E and H are the envelope functions. For a guided mode with a propagation constant β and an azimuthal mode order l, we can write where e and h are the reduced mode profile functions and β and l can take not only positive but also negative values. The mode profile functions e and h can be decomposed into 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 [42].

Hybrid modes
It is convenient to introduce the parameters s, Then, we find, for r < a, and, for r > a, e r = iA β 2q and Here, the parameter A is a constant that can be determined from the propagating power of the field.

TE modes
For r < a, we have e r = 0, and For r > a, we have e r = 0, and h r = i β q J 0 (ha) K 0 (qa) AK 1 (qr),

TM modes
For r < a, we have and h r = 0, For r > a, we have and h r = 0,

Appendix B: Poynting vector
In this appendix, we calculate the Poynting vector S for the different mode families. First we note that, for guided modes of fibers, the radial component of the Poynting vector is always zero, that is, S r = 0.

Hybrid modes
For quasicircularly polarized hybrid modes, we find that the axial and azimuthal components of the Poynting vector are given, for r < a, by and, for r > a, by Note that the expressions for S ϕ in Eqs. (B1) and (B2) contain cross terms of the types J l±1 (hr)J l (hr) and K l±1 (qr)K l (qr). These terms appear as a result of the interference between different terms associated with different Bessel functions. Due to the interference, the azimuthal component S ϕ of the Poynting vector of a quasicircularly polarized hybrid mode may have different signs in different regions of space.
For quasilinearly polarized hybrid modes, we find that the axial component of the Poynting vector is given, for r < a, by and, for r > a, by In both regions, we have S ϕ = 0.
It is worth noting that expressions (B3) and (B4) for S z contain cross terms of the types J l−1 (hr)J l+1 (hr) and K l−1 (qr)K l+1 (qr). A similar argument as the one above confirms that the axial component S z of the Poynting vector of a quasilinearly polarized hybrid mode can have different signs in different regions of space.

TE modes
For TE modes, we find that the axial component of the Poynting vector is given, for r < a, by and, for r > a, by The azimuthal component is zero, that is, S ϕ = 0. It is clear from Eqs. (B5) and (B6) that, for TE modes, the axial component S z of the Poynting vector is positive with respect to the direction of propagation, in agreement with the results of Ref. [44].

TM modes
For TM modes, we find that the axial component of the Poynting vector is given, for r < a, by S z = f |A| 2 ωǫ 0 n 2 1 β 2h 2 J 2 1 (hr) (B7) and, for r > a, by S z = f |A| 2 ωǫ 0 n 2 2 β 2q 2 J 2 0 (ha) K 2 0 (qa) The azimuthal component is zero, that is, S ϕ = 0. It is clear from Eqs. (B7) and (B8) that, for TM modes, the axial component S z of the Poynting vector is positive with respect to the direction of propagation, in agreement with the results of Ref. [44].

Appendix C: Power
The propagating power of a guided mode is P = P in +P out , where P in and P out are the propagating powers inside and outside the fiber.
For EH lm modes, the cutoff value V c is determined as a nonzero solution to the equation J l (V c ) = 0. In the limit V → V c for EH lm modes, the parameters s and s 2 tend to a limiting value s c , where s c = −1. Consequently, the term in the last line of Eq. (C2) is dominant. On the other hand, in the limit qa → 0, we have K l (qa) ≃ (1/2)(l − 1)!(2/qa) l for l ≥ 1 [94]. With the use of the boundary conditions for the field at the fiber surface, we can show that J l (ha) = O(q 2 a 2 ). When we use the aforementioned asymptotic expressions, we can show that P out tends to a finite value. Meanwhile, P in tends to a nonzero finite value. Consequently, in the limit V → V c for EH lm modes, the fractional power outside the fiber η P tends to a limiting value that is smaller than unity.
For HE lm modes, the cutoff value V c is not a solution to the equation J l (V c ) = 0 except for the case of l = 1. In the limit V → V c for HE lm modes, we have s = −1 + O(q 2 a 2 ) and s 2 = −1 + O(q 2 a 2 ). When we use these asymptotic expressions and the approximate expression K l (qa) ≃ (1/2)(l − 1)!(2/qa) l for l ≥ 1 [94], we can show that, for l ≥ 3, the power outside the fiber P out tends to a finite value. Meanwhile, P in tends to a nonzero finite value. Consequently, in the limit V → V c for HE lm modes with l ≥ 3, the fractional power outside the fiber η P tends to a limiting value that is smaller than unity.
The analysis in the above paragraph is not valid for the HE lm modes with l = 1, 2. Indeed, for l = 1, 2, the expression on the right-hand side of Eq. (C2) contains the modified Bessel function K 0 (qa). The asymptotic expression for this function with a small argument qa is K 0 (qa) ≃ − ln(qa/2) − γ, where γ ≃ 0.5772 is the Euler-Mascheroni constant [94]. In addition, for l = 1, the cutoff value V c is a solution to the equation J 1 (V c ) = 0, and the corresponding magnitude of J 1 (ha) in the limit V → V c is found to be on the order of 1/| ln qa|. Then, we can show that, in the limit V → V c for the HE lm modes with l = 1, 2, we have P out → ∞. Meanwhile, P in tends to a finite value. Consequently, in the limit V → V c for the HE lm modes with l = 1, 2, the fractional power outside the fiber η P tends to unity.

TE modes
For TE modes, the explicit expressions for the powers inside and outside the fiber are [42] P in = f |A| 2 πa 2 ωµ 0 β 2h 2 [J 2 1 (ha) − J 0 (ha)J 2 (ha)] (C3) and P out = f |A| 2 πa 2 ωµ 0 β 2q 2 J 2 0 (ha) K 2 0 (qa) [K 0 (qa)K 2 (qa) − K 2 1 (qa)]. (C4) We calculate the fractional power outside the fiber η P for a TE mode in the limit where the fiber size parameter V tends to the cutoff value V c . This cutoff value V c is determined as a solution to the equation J 0 (V c ) = 0. In the limit V → V c , we have qa → 0 and ha → V c . When we use the eigenvalue equation (A4) for TE modes and the asymptotic expressions for the modified Bessel functions K 0 (qa) and K 1 (qa) with a small argument qa, we find from Eq. (C4) that P out → ∞. Meanwhile, P in tends to a finite value. Consequently, in the limit V → V c for TE modes, the fractional power outside the fiber η P tends to unity.
(C6) We calculate the fractional power outside the fiber η P for a TM mode in the limit where the fiber size parameter V tends to the cutoff value V c . The cutoff value V c is determined as a solution to the equation J 0 (V c ) = 0. In the limit V → V c , we have qa → 0 and ha → V c . When we use the eigenvalue equation (A5) for TM modes and the asymptotic expressions for the modified Bessel functions K 0 (qa) and K 1 (qa) with a small argument qa, we find from Eq. (C6) that P out → ∞. Meanwhile, P in tends to a finite value. Consequently, in the limit V → V c for TM modes, the fractional power outside the fiber η P tends to unity.

TE modes
For TE modes, the explicit expressions for the energies per unit length inside and outside the fiber are found to be U in = |A| 2 πa 2 µ 0 4 J 2 0 (ha) + 2n 2 1 k 2 h 2 J 2 1 (ha) and U out = |A| 2 πa 2 µ 0 4