Polarization transport in optical fibers beyond Rytov's law

We consider the propagation of light in arbitrarily curved step-index optical fibers. Using a multiple-scales approximation scheme, set-up in Fermi normal coordinates, the full vectorial Maxwell equations are solved in a perturbative manner. At leading order, this provides a rigorous derivation of Rytov's law. At next order, we obtain non-trivial dynamics of the electromagnetic field, characterized by two coupling constants, the phase and the polarization curvature moments, which describe the curvature response of the light's phase and its polarization vector, respectively. The latter can be viewed as an inverse spin Hall effect of light, where the direction of propagation is constrained along the optical fiber and the polarization evolves in a frequency-dependent way.


I. INTRODUCTION
Optical fibers are widely used as waveguides for electromagnetic fields, providing a controlled way of transporting light along a given path. They play a fundamental role in many engineering applications, such as telecommunications [1] and sensors [2]. Furthermore, optical fibers play a central role in many areas of physics, such as metrology [3], photonics [4], quantum information [5] and quantum computation [6][7][8], and are planned to be used in future experiments that probe the interplay of quantum optics and Einstein's theory of gravity [9][10][11][12].
The propagation of electromagnetic waves in optical fibers is generally described by Maxwell's equations. For straight optical fibers with circular cross sections, Maxwell's equations can be solved explicitly using certain mode decompositions adapted to the symmetry of the problem [13][14][15]. However, in real-world applications, optical fibers are generally twisted and bent in various ways, leading to optical losses, non-trivial polarization dynamics, and coupling between different modes. In this case, exact solutions to Maxwell's equations are no longer available, which makes it necessary to use approximation schemes.
A geometrical-optics treatment of light rays following nonplanar curves was already given in 1938 by Rytov [16], and later extended by Vladimirskii [17] (see Ref. [18] for a more recent overview of these results). In these papers, the authors derived a transport law for the polarization vector (defined as a unit vector pointing along the direction of the electric field) along the nonplanar curve followed by light rays in inhomogeneous media. This is known as Rytov's law, which states that the polarization vector is Fermi-Walker transported along the curve [19,Ch. 6.1]. In optical fibers, Rytov's law was experimentally observed by Ross [20], followed by a similar experiment by Tomita and Chiao [21] (see also Ref. [22]). A theoretical discussion * thomas.mieling@univie.ac.at; ORCID: 0000-0002-6905-0183 † marius.oancea@univie.ac.at; ORCID: 0000-0002-1242-4041 describing the geometry of this polarization transport law in optical fibers was given by Haldane in Refs. [23,24]. We give a simple illustration of Rytov's law in Fig. 1, where we consider a helical fiber, together with the Fermi-Walker transport of two linear polarization vectors. The standard derivation of Rytov's law relies on the geometrical optics approximation, which breaks down for electromagnetic waves propagating in single-mode optical fibers (SMFs) both because the wavelength is of the same order of magnitude as the diameter of the fiber and because of the rapid formation of caustics. An extension of Rytov's law for light propagating in SMFs was given by Berry [25], where a term proportional to the wavelength was added to the polarization transport law. Similar results were also obtained by Lai et al. [26] using different methods.
However, Berry did not provide an explicit derivation of his result from Maxwell's equations, and the work of Lai et al. does not consider the junction conditions of the electromagnetic field at the core-cladding interface, which renders the relation of their result to optical fibers unclear.
In this paper, we perturbatively solve for the full electromagnetic field in an arbitrarily bent step-index optical fiber under the sole assumption that the radius of curvature of the fiber is much larger than the wavelength of light propagating therein. To do so, we erect cylindrical Fermi coordinates around the baseline of the optical fiber, such that the core-cladding interface is at a constant radial coordinate and the effects of bending are fully encoded in a single component of the metric tensor.
The assumption of weak bending then allows setting up a perturbative scheme based on a multiple-scales method, where the perturbation parameter is given by the ratio of the optical wavelength to the fiber's radius of curvature. We mainly focus on the single-mode regime, but we shall also discuss briefly how our method can be applied in certain multi-mode regimes.
At first order in this perturbative expansion, we recover Rytov's law, which means that the Jones vector is Fermi-Walker transported along the fiber. At next order, we obtain corrections to this transport law, describing FIG. 1. Illustration of Rytov's law in an optical fiber. The two vectors e1 (red) and e2 (green), orthogonal to the tangent vector T = γ ′ (not shown), are Fermi-Walker transported along the curve γ, and they can be viewed as linear polarization vectors aligned with the electric field lines of an electromagnetic wave with horizontal or vertical linear polarization. For example, a linearly polarized electromagnetic wave initially at the lower end of the fiber will experience an counterclockwise rotation of its plane of polarization as it propagates upwards along the helix.
non-trivial dynamics of the electromagnetic phase and polarization due to the bending of the fiber. This correction to the polarization dynamics can be interpreted as an inverse spin Hall effect of light: whereas in the (direct) spin Hall effect of light [27][28][29][30][31][32] the polarization influences the trajectory of light, here the trajectory is prescribed by the optical fiber, which then influences the light polarization. Similar inverse spin Hall effects of light have also been reported in other physical systems [33,34]. This paper is structured as follows. In Section II, we introduce the Fermi coordinate system used in the subsequent calculation and describe the relation of the metric tensor components (in this coordinate system) to the curvature and torsion of the fiber. Section III describes the perturbative ansatz for the electromagnetic potential (using an adapted orthonormal frame along the fiber and the aforementioned multiple-scales approach), and formulates the perturbative equations and their solutions. The expressions that arise there require numerical quadrature. In Section IV we present numerical results for the polarization moments, which describe the coupling of the electromagnetic phase and polarization to the bending of the fiber, and illustrate the polarization transport for the simple case of a helical fiber. Finally, in Section V we compare with previous results obtained by other techniques and provide an outlook on potential future applications of our methods.

II. GEOMETRIC SETUP
Let γ be a curve in three-dimensional Euclidean space representing the baseline of an SMF. In the vicinity of this curve, one can erect a Fermi coordinate system (t, x, y, s), where t is time, s is the arc length of γ, and x and y are transverse coordinates. More details are given in Appendix A. In this coordinate system, the Minkowski metric takes the form where c is the speed of light in vacuum and ν i are the components of the normal vector of the curve in the Fermi coordinate system.
Henceforward, we assume that the fiber is only weakly bent in the following sense. First, we require the radius of curvature 1/κ = 1/ ν 2 x + ν 2 y to be much larger than the SMF's core radius ϱ, i.e. ϱκ ≪ 1. In practice, this requirement is well satisfied, since typical SMFs, bent on the scale of centimeters, have ϱκ ∼ µm/cm = 10 −4 . Second, we require the normal vector ν i to be not only small in norm (of order κ ≪ 1/ϱ), but also to vary slowly when measured in units of the radius ϱ. This means that ν i can be expressed in the form where ε is a small parameter and b i are dimensionless functions whose derivatives are at most of the same order of magnitude as the functions themselves. In practice, since SMF modes have wavelengths comparable to the core radius ϱ, this means that the normal vector is almost constant over an optical wavelength and changes significantly only over a length scale of 1/ε ≫ 1 wavelengths.
To simplify the notation, all subsequent calculations will be carried out in nondimensionalized cylindrical coordinates (t, r, ϑ, σ), defined by t = ϱt/c , x = ϱr cos ϑ , y = ϱr sin ϑ , s = ϱσ . (3) In this coordinate system, the Minkowski metric takes the form where with ς = εσ and b ± being the complex bending functions These complex functions are related to the fiber's curvature κ and torsion τ by where the arbitrary constant describes the freedom to rotate the coordinate system rigidly around the r = 0 axis, cf. Eq. (A19).

III. MAXWELL'S EQUATIONS
The electromagnetic field in optical fibers satisfies the source-free Maxwell equations when using either Gaussian or Lorentz-Heaviside units. For non-magnetic materials, the constitutive relations are with n denoting the refractive index. For step-index optical fibers, n has the form where the refractive indices of the core and the cladding, n 1 and n 2 , are constants satisfying n 1 > n 2 ≥ 1. Here, the operators div and curl are those of flat Euclidean space, which are to be expressed in non-Cartesian coordinates when working with the curvilinear Fermi coordinate system.

A. Gauge-Fixed Field Equations
Here, we work with the electromagnetic potentials ϕ and A, such that Imposing an appropriate generalization of the Lorenz gauge to linear isotropic dielectrics, as in Ref. [12], the field equations reduce to wherever n is constant. Here, ∆ denotes the Laplacian (more precisely: the scalar Laplace-Beltrami operator when acting on ϕ, for the vector Laplace-Beltrami operator when acting on A) of Euclidean space, which we express in curvilinear coordinates below. Equation (13) can also be written in generally covariant form. In regions of constant n, the electromagnetic fourpotential A = −cϕdt + A i dx i satisfies the wave equatioñ (in abstract index notation) where ∇ a is the Levi-Civita derivative of the metric g ab given in Eq. (4) andg ab is the inverse of Gordon's optical metricg ab [35], which here takes the form To decouple the wave equations as much as possible, we use an adapted frame e µ and a coframe e µ : (17) with respect to which the electromagnetic four-potential can be decomposed as To express the junction conditions at the core-cladding interface r = 1 (derived in Ref. [12]) in a form suitable for the following calculations, define the column vector where dots and primes denote partial differentiation with respect tot and r, respectively, and, for any function φ, the expression φ denotes the discontinuity of φ at the interface: The junction conditions can then be expressed as Ψ[A] = 0, which means that all field expressions entering Eq. (19) must be continuous.

B. Perturbative system
In the absence of bending, the fiber modes have constant angular frequency ω and propagation constant β, for which we define the dimensionless counterparts To set up a perturbative scheme for Eq. (13) in the presence of bending, we make the following ansatz. Seeking monochromatic waves of angular frequency ω and propagation constant β, we write where the field's angular dependence is expanded in a Fourier series: with the following perturbative expansion Here, σ and ς = εσ are to be considered as "independent variables" in the sense of the multiple-scales method [36]. The precise form of the last expansion is motivated as follows. In the unperturbed problem, SMFs allow for modes with m = +1 and m = −1 only [13,14]. Hence, we start the series expansion of these Fourier modes at order ε 0 . The bending terms couple neighboring Fourier components, so the series expansions of m = 0 and m = ±2 start at order ε 1 , while Fourier modes with m = ±3 are at most of order ε 2 , and so on. We insert this ansatz into Eq. (13) and consider terms of various powers of ε separately.

Unperturbed system
At order ε 0 , one obtains where the Helmholtz operator H m acts on the fields as follows: with The solution which is regular on the axis r = 0 and decays for r ≫ 1 is given in terms of Bessel functions as The coefficients q (m,k) µ,l may depend on ς, and the functions f (m) are defined in terms of Bessel functions of the first kind, J m , and modified Bessel functions of the second kind, K m , as with For conciseness, we shall abbreviate Eq. (29) as where q (m,0) is a column vector containing all the parameters q that occur in Eq. (29).
The matching conditions at the core-cladding interface r = 1 can then be written in matrix form where M m is a complex 8 × 8 matrix, given explicitly by Eq. (58) in Ref. [12], the precise details of which are not important for our present considerations. For Eq. (33) to admit a non-trivial solution, the determinant of M ±1 must vanish (one can show that the two determinants for m = +1 and m = −1 are proportional, so that they must vanish simultaneously). This determines the standard dispersion relation between β, ω, ϱ and the refractive indices n 1 , n 2 [13,14].
If the dispersion relation is satisfied (SMFs admit only a single solution for β in terms of ω), the matrices M ±1 possess a one-dimensional kernel and co-kernel. Settinĝ q ± to be any vector spanning the kernel of M ±1 (we takê q ± to be ς-independent), the general solution to Eq. (33) can be written as where J ± on ς, we must consider the equations at next order in the perturbative expansion.

First order perturbations
At order ε 1 , one obtains differential equations of the form where the inhomogeneities Σ (m,1) µ , given explicitly in Appendix B, depend on the fields at order ε 0 , their first ς-derivatives, as well as on the bending functions b ± .
These equations can be solved using Green's functions for the Helmholtz operator H m [37,Ch. 1.16]. Setting where Y m and I m denote Bessel functions of the second kind and modified Bessel functions of the first kind, respectively, the general solution can be written in the form The free parameters q (m,1) are needed to satisfy the junction conditions at the core-cladding interface. As in the unperturbed case, these matching conditions can be written in matrix form: Since the matrices M 0 and M ±2 are non-singular, Eqs. (41a) and (41c) can be solved uniquely for the coefficients q (0,1) and q (±2,1) . However, M ±1 is singular due to the dispersion relation. The condition that Eq. (41b) is solvable is equivalent to the right-hand side of Eq. (41b) being in the image of the matrix M ±1 . Denoting by ζ ± any eight-component row vector that spans the co-kernel of M ±1 , the solubility condition can be expressed as Expanding all definitions, one obtains equations of the form where the multiplicative factor is expressible in terms of integrals of Bessel functions. Numerically, we find this factor to be non-zero, leading to As will be explained in more detail below, this equation implies that, to first order in the perturbative expansion, the Jones vector of the SMF is Fermi-Walker transported along the fiber (Rytov's law). Non-trivial dynamics arise only at second order. As Eq. (44) implies that Eq. (41b) reduces to M ±1 · q (±1,1) = 0, it follows that the vectors q (±1,1) are proportional to the vectorsq ± spanning the kernel of M ±1 where the proportionality coefficients may depend on ς. Hence, we obtain where the dynamics of the coefficients J ± is to be determined at the next order.

Second order perturbations
At order ε 2 , one obtains field equations of the form where the inhomogeneities Σ (m,2) µ , given explicitly in Appendix B, depend on the fields at order ε 0 and ε 1 , their ς-derivatives, as well as on the bending functions and their derivatives.
These equations can be solved using the Green's operator defined above: and the interface conditions take the abstract form The structure is similar to Eq. (41). Equations (48a), (48c) and (48d) can be solved uniquely for the coefficients q (0,2) , q (±2,2) , and q (±3,2) . However, as the matrices M ±1 are singular, the solubility of Eq. (48b) provides a non-trivial constraint, which is equivalent to Expanding the definitions, this leads to an equation of the form d dς whereB is a complex 2 × 2 matrix, which we find to be of the general form withξ andη being dimensionless constants that depend on the details of the SMF. (The factors of two on the diagonal were inserted to simplify later calculations using 2b + b − = ϱ 2 κ 2 /ε 2 .) The matrixB is Hermitian since the bending functions b ± are each other's complex conjugates and the parametersξ andη are real (this was verified numerically, see Section IV below).

C. Polarization transport
Combining J (0) ± and J (1) ± into one Jones vector Eqs. (44) and (50) lead to where error terms of order ε 3 have been neglected.
Transforming the complex components J ± to Cartesian components J x and J y according to the evolution equation for the Jones vector J can be written as d ds Extending the Jones vector J to a three-dimensional vector by setting J ∥ = 0, this result can be written in covariant vector notation as where D s is the spatial Fermi-Walker derivative along the fiber's baseline, defined explicitly in Eq. (A4) below.

IV. RESULTS
We have implemented the above calculations in Wolfram Mathematica [38] to numerically compute the phase curvature moment η and the polarization curvature moment ξ for various optical fibers and optical wavelengths. Here, we discuss numerical results for telecommunication SMFs and optical nanofibers, and also present an analytical solution to the transport law (56) for helical optical fibers. Figure 3 shows the dependence of the phase curvature moment η and the polarization curvature moment ξ on the vacuum wavelength λ = 2πc/ω for two typical telecommunication optical fibers with a core radius ϱ = 4.1 µm and refractive indices (i) n 1 = 1.4712 and n 2 = 1.4659 as well as (ii) n 1 = 1.4715 and n 2 = 1.4648.

A. Wavelength-dependence in single-mode fibers
In his 1987 paper [25], Berry predicted the polarization curvature moment to be ξ ≈ 1/2β = λ/(4πn) for weakly guiding fibers, wheren is the effective refractive index. Despite the fact that the fibers considered here are weakly guiding (the relative index differences are ∆ i = 0.36 % and ∆ ii = 0.45 %, respectively), we find a deviation from this prediction. Instead of ξ depending linearly on λ, the curves closely follow an affine dependence on λ. The function ξ(λ) thus has an isolated zero, at which Rytov's law extends to second order in perturbation theory. However, as the following examples show, these curves are not affine throughout.

B. Radius-dependence in nanofibers
In recent years, nanofibers, i.e., optical fibers with submicrometer diameters, have attracted much experimental interest [39][40][41]. In this regime, the two curvature moments η and ξ exhibit a behavior which differs significantly from that in telecommunication SMFs described above. Figure 4 shows the two curvature moments for nanofibers with n 1 = 1.46 (fused silica) and n 2 = 1 (vacuum), operated at an optical vacuum wavelength of λ = 633 nm, as was used experimentally in Ref. [39].
As such fibers cease to be single-mode at ϱ ≳ 230 nm (where modes with m = 0 appear) and even support multiple modes with m = ±1 (labeled by their radial mode index), this plot shows that the polarization curvature moments are continuous functions at the cutoff frequencies of higher modes.

C. Multi-mode regime
In the above calculations, we have mainly focused on SMFs, as the single-mode condition guarantees that the matrices M m with m ̸ = ±1 are non-singular. However, we observe numerically that these matrices are typically also non-singular for fibers which are not single-mode, thus allowing us to extend our calculations into the multimode regime. It should be noted, however, that not all higher-order modes are linearly polarized, in which case Eq. (56) continues to hold on a formal level, but the physical interpretation of J becomes less intuitive (since it merely describes the relative weights of the solutions with m = +1 and m = −1, without a direct relation to a polarization vector).
To investigate the behavior of the curvature moments in the multi-mode regime, we consider optical fibers with n 1 = 1.46 and n 2 = 1 (as in the preceding section) with the normalized frequency V = ϱωc −1 n 2 1 − n 2 2 up to V = 10, where such fibers support five modes with azimuthal mode index m = ±1 (they support further modes for different mode indices, but they are not covered by the calculations presented above). As Fig. 5 shows, the curvature moments grow steeply near the cutoff frequencies, which are marked by vertical dashed lines. This indicates a limit on the range of validity of our perturbative scheme, as first-order perturbations of the form εξ and εη cease to be small near the cutoff.

D. Helical optical fibers
To illustrate the polarization dynamics derived above, we consider the example of an optical fiber whose baseline forms a helix. In this case, the Fermi-Walker transport for the orthonormal frame (A4), as well as the polarization transport equation for the Jones vector (56), can be solved exactly.
We start by defining the curve γ, representing the baseline of the optical fiber, as γ(s) = (R cos(αs), R sin(αs), Hαs) , where s is the arc length, R is the radius of the helix, 2πH is the increase in height of the helix along the z-axis after one turn in the xy-plane, and α = (R 2 + H 2 ) −1/2 . The curve γ has constant curvature and torsion and the Fermi-Walker transported orthonormal frame (T, e 1 , e 2 ) can be constructed explicitly using Eq. (A5). This frame is illustrated in Fig. 1. At some reference point s = 0 on the fiber, the vectors e 1 and e 2 can be aligned with the initial orientation of transverse electric field lines of the linearly polarized fiber modes given in Fig. 2. Then, the Fermi-Walker transport of e 1 and e 2 along γ describes the lowest-order frequency-independent transport law for the polarization of the electromagnetic field along the optical fiber. This is referred to as Rytov's law, and has been experimentally observed in Refs. [20,21].
After traveling along ℓ complete loops of the helical fiber, the plane of polarization of a linearly polarized electromagnetic wave is rotated by the angle This is the angle between the vectors e i (s = 0) and e i (s = 2πℓ/α). Our results also include higher-order frequencydependent corrections to the polarization dynamics given by Fermi-Walker transport, which can be viewed as an inverse spin Hall effect of light. All information about the polarization of the electromagnetic wave traveling along the optical fiber, relative to the Fermi-Walker transported frame, can be described by the Jones vector J = (J + , J − ) T . The dynamics of the Jones vector is given by the transport equation Eq. (53), which, for the helical fiber considered here, takes the form The transport equation can be solved exactly in this case, and the general solution can be written as and c 1 and c 2 are integration constants depending on the initial values J(0). Explicitly, one has with the matrix U (s) being given by 2k sin(ks) cos(ks) + iτ k sin(ks) .
For illustration, consider an electromagnetic wave that is initially right-hand circularly polarized, thus J(0) = (1, 0) T . Then, using the exact solutions given above, one can calculate the probability of observing a beam of opposite circular polarization to be This change of polarization represents an inverse spin Hall effect of light, where the polarization dynamics is controlled by the wavelength and by the bending of the optical fiber. For typical fiber parameters such as R = 10 cm, H = 1 mm and ξ = 10 µm, the above transition probability oscillates with an amplitude 1 − τ 2 /k 2 ≈ 2.5 × 10 −5 and a period πk −1 ≈ 31 m. Furthermore, note that k depends on the polarization curvature moment ξ, which, for telecommunication fibers, is an affine function of the optical wavelength λ (see Fig. 3), but can have a more general behavior for nanofibers (Fig. 4) or for higher-order modes (Fig. 5). As ξ controls the magnitude of the inverse spin Hall effect of light, the results in Fig. 5 suggest that this effect becomes large and potentially observable close to the mode cutoff frequencies or for large values of the normalized frequency V .
When the torsion τ is set to zero, the helix reduces to a circle of curvature κ = R −1 , and the above transition probability becomes Thus, a circularly polarized electromagnetic wave can flip to the opposite state of circular polarization after propagating along a circular fiber. This occurs after the electromagnetic wave has traveled a certain length s along the fiber, depending on the geometry (through κ), the parameters of the fiber, and the wave frequency (ξ depends on the refractive indices in the core and cladding, on the radius of the core, and on wave frequency). More precisely, we will have P +− (s) τ =0 = 1 for s = (2i + 1) π ξκ 2 , with i ∈ Z. A similar transition probability was also obtained in Ref. [26,Eq. 32]. However, our result is different in the sense that ξ is an affine function (at least in the case of the telecommunication fibers presented in Fig. 3) instead of a linear function of λ, as was the case in Ref. [26]. Furthermore, in our case ξ contains additional information about the optical fiber parameters, such as the core radius ϱ, and the refractive indices of the core and cladding.
We can also calculate the transition probability from an initially right-handed circularly polarized wave to a linearly polarized wave: Note that while P +− (s) has a period of πk −1 , P +x (s) involves products of two sinusoidal functions with periods of πk −1 and πτ −1 , respectively. For the fiber parameters mentioned above, we have k ≈ τ , with k − τ ≈ 1.2 × 10 −6 m −1 . In this case, P +x (s) will oscillate on short length scales with frequency k ≈ τ , and we can write In this approximation, the amplitude of the oscillations of P +x (s) is controlled by the wave frequency and fiber parameters (through ξ) and by the geometry (through κ and τ ), while the frequency of the oscillations of P +x (s) is solely controlled by the geometry (through τ ). In particular, a transition from circular to linear polarization is achieved when P +x (s) = 1, although this would require some fine tuning of the geometry, fiber parameters and wave frequency. In general, there will also be a long-scale modulation of P +− (s) due to the small difference between the two frequencies in Eq. (67). However, this additional modulation is only significant on length scales of the order π(k − τ ) −1 ≈ 2.5 × 10 6 m and can thus be neglected.

V. DISCUSSION
We have derived an equation for the transport of the electromagnetic field's polarization and phase in curved optical fibers by perturbatively solving Maxwell's equations, based on a multiple-scales approximation scheme. The response of the electromagnetic field to the bending was found to be characterized by two coupling constants: the polarization curvature moment ξ and the phase curvature moment η. We showed how the polarization curvature moment ξ leads to an inverse spin Hall effect of light, where the state of polarization of the electromagnetic wave is controlled by the wavelength and the bending of the optical fiber.
The equations for the polarization transport derived in this paper are similar in structure to those derived by Berry [25]. In both treatments, the leading-order transport law is given by Fermi-Walker transport, while second-order corrections depend on the fiber's curvature. However, the details of the second-order corrections differ: in Berry's result, second-order corrections depend linearly on frequency and quadratically on curvature, but there is no explicit dependence on fiber parameters (neither on the core radius nor on the refractive indices). Furthermore, Berry's paper does not contain a detailed derivation (the publication of which was announced in the paper, but we are not aware of it). By contrast, our explicit derivation shows that second-order corrections generally can have a more complicated frequency dependence than predicted by Berry, as described by the phase and polarization curvature moments (which depend on the fiber parameters).
Similar results to Berry's were also obtained in the related paper by Lai et al. [26]. While our calculations model the electromagnetic field in both the fiber core and the fiber cladding (with appropriate interface conditions), Lai et al. evaluated the field equations only up to the fiber core radius (i.e. in the core), suggesting that the field was approximated to vanish in the cladding-an approximation which was not necessary in the calculations presented in this paper and is not satisfied in practice for nanofibers.
Moreover, our methods do not rely on the commonly used weak-guidance approximation, nor do they make use of the thin-layer method [26,42]. Instead, the methods presented here implement a perturbative scheme for the full vectorial Maxwell equations, based only on the assumption of weak and slowly varying bending of the fiber (when compared with the wavelength of the fields therein).
We expect the methods described here to be capable also of providing detailed analyses of fiber mode perturbations by other effects, such as variations in the refractive indices due to stresses, variations in the fiber core radius, but also gravitational effects.

ACKNOWLEDGMENTS
We are grateful to Lars Andersson, Piotr Chruściel, and Christopher Hilweg for helpful discussions. T.M. is a recipient of a DOC Fellowship of the Austrian Academy of Sciences at the University of Vienna, Faculty of Physics, and is supported by the Vienna Doctoral School in Physics (VDSP), the research network TURIS, and in part by the Austrian Science Fund (FWF), Project No. P34274, as well as by the European Union (ERC, GRAVITES, 101071779).

Appendix A: Spatial Fermi coordinate system
This section summarizes the construction of Fermi normal coordinates, based on Fermi-Walker transport. For an analogous construction based on the Serret-Frenet frame, see, e.g., Refs. [26,43].
Let (M, g) be an oriented three-dimensional Riemannian manifold with Levi-Civita connection ∇. Let s → γ(s) be a curve parameterized by arc length, and denote by T = γ ′ its tangent vector (with g(T, T ) = 1) and by ν = ∇ T T its normal vector. In regions where the curvature κ = g(ν, ν) is non-zero, one can introduce the orthonormal Serret-Frenet frame where N is the unit normal and B is the binormal. The torsion of the curve γ is then defined as and the evolution of the Serret-Frenet frame along the curve is given by the Serret-Frenet equations [44,Ch. 7.B] The Fermi-Walker derivative D s is defined by its action on any vector v as One readily verifies that D s T = 0, D s N = τ B, D s B = −τ N and D s g = 0. Now, let (T, e 1 , e 2 ) be an orthonormal frame at any "reference point" p = γ(0). Extending e 1 , e 2 along γ by Fermi-Walker transport, i.e. D s e 1 = De 2 = 0, one obtains an orthonormal frame (T (s), e 1 (s), e 2 (s)) along the entire curve γ. The Fermi-Walker transported vectors e 1 (s) and e 2 (s) can be explicitly constructed along any curve γ with nonvanishing curvature κ as where exp : T M → M is the exponential map. The implicit function theorem guarantees that the parameters (s, x, y) define a coordinate system in a neighborhood of the curve γ. As (∂ s , ∂ x , ∂ y ) = (T (s), e 1 (s), e 2 (s)) is an orthonormal frame along γ, one finds that along γ the metric tensor takes the form The Christoffel symbols along γ can be determined using the general formula Using the definition of the Fermi-Walker transport, one finds from which one can read Γ ijs = Γ isj and Γ iss . It thus remains to determine Γ ijk where the last indices are both distinct from s. This is accomplished by considering the curve ξ → ψ(s, c 1 ξ, c 2 ξ), where s, c 1 and c 2 are constant. By Eq. (A7), this is a geodesic, which implies that c 1 ∂ x + c 2 ∂ y is auto-parallel, which means Since the constants c 1 and c 2 are arbitrary, one finds that all covariant derivatives entering this equation vanish, hence Γ ijk vanishes whenever the last two indices are both distinct from s. In total, one thus finds that the only non-zero Christoffel symbols (up to symmetry) along γ are where ν i = (0, ν x , ν y ) are the components of the normal vector ν in the frame (T, e 1 , e 2 ) (they are functions of s alone). Using ∂ i g jk = Γ jki + Γ kij , one finds that the only non-zero derivatives of g ij along γ are Expanding the metric along γ in a Taylor series in x and y, one obtains g ij = δ ij − 2δ s i δ s j (ν x x + ν y y) + O(x 2 ) + O(y 2 ) . (A14) Up to error terms of order r 2 = x 2 + y 2 , one thus has g = dx 2 + dy 2 + (1 − ν i x i ) 2 ds 2 + O(r 2 ) .
While the calculations leading to Eq. (A15) were perturbative, it turns out that the error terms vanish in flat Euclidean space. This can be seen by noting that the metric of Eq. (A15) without the error terms has vanishing Riemann curvature, and that the xy-plane is totally geodesic, so that straight coordinate lines in the xy-plane are exact solutions to the geodesic equation. Error terms thus arise only in the presence of curvature, and are thus absent for the applications in this paper.

Appendix B: Inhomogeneities
Here, we give explicit expressions for the inhomogeneities Σ (m,k) µ arising in Eqs. (35) and (46). To express them in a concise way, we introduce some notation.
First, define the two operators Then, denoting by a µ the column vector of field components defined in Eq. (26), let A, B and C be defined as follows: A − m a µ =     L 2 +m a 0 L 2 +m a ∥ L 2 +(m+1) a + L 2 With these definitions at hand, the inhomogeneities of the first-order equations can then be written as Further, the inhomogeneities of the second-order equations are These source terms do not couple the temporal field component a 0 with the spatial components a ∥ and a ± , which is consistent with the decoupled system given in Eq. (13).