On the Influence of Earth's Rotation on Light Propagation in Waveguides

We analyse the influence of Earth's rotation (both around its own axis and around the Sun) on the propagation of light in optical media. This is done using both geometrical optics and a perturbative calculation based on Maxwell's equations in rotating coordinates in flat spacetime. Considering light propagation in cylindrical step-index waveguides in particular, the first order correction to electromagnetic modes is computed. The calculation shows that Earth's rotation causes a weak mode coupling, giving rise to sidebands, whose amplitudes are computed as well. The correction to the dispersion relation derived here allows to assess the anisotropy of light propagation due to Earth's rotation. The linearisation of this result is found to agree numerically with a simple formula derived from geometrical optics.


Introduction
The problem of determining the influence of Earth's rotation on optical phenomena was raised in the second half of the 19 th century and was first discussed using the theory of a luminiferous aether, see e.g. [8,9]. Using an argument based on Huygen's principle, Lorentz concluded that the effective refractive index n eff of any medium is related to the ordinary refractive index n by cf. [8, p. 126ff], where v is the velocity of the supposed aether and m is the unit vector in the direction of light propagation. One of the main results of this paper is the derivation of an identical formula, by analysing the Maxwell's equations in a rotating system in Minkowski spacetime. Here v is instead the negative linear velocity due to rotation of Earth, as measured from an inertial system at its center. Moreover, we compute the first order corrections to both electric and magnetic field explicitly. In particular, we derive corrections to the amplitudes of the modes and their sidebands, which are excited due to rotation.
The influence of Earth's rotation (about its own axis) on light was first demonstrated experimentally using a Sagnac interferometer [10], where a beam of light is split into two, which traverse a ring interferometer in opposite sense of rotation. If the apparatus rotates, the different light rays take unequal times to complete a full circle, leading to an interference pattern whose fringes are displaced from the positions they would have if the apparatus were an inertial system. Although Sagnac interpreted this effect as a proof for the existence of aether [11,12], Laue gave an explanation in the framework of Special Relativity [13] and Langevin was able to explain the phenomenon in the language of General Relativity [6,7].
In a 2017 paper [4], Hilweg et al. proposed an experiment to measure the effect of Earth's gravitational potential on single photons. This paper included an estimate for the influence of Earth's rotation, where the photon was modelled as a classical particle traversing a spooled waveguide. In a more recent paper [2], the shift of the wave vector due to Earth's gravitational field was calculated more accurately using Maxwell's equations in a post-Newtonian metric.
In this work, the effect of Earth's rotation is determined from first principles, using Maxwell's equations in a rotating frame in Minkowski spacetime, thereby supplementing the well-known geometrical optics argument by more accurate wave optics calculations.

Outline of this Work
Classical propagation of light can be described either by geometrical optics or wave optics. The latter method describes light by wave solutions of Maxwell's equations whereas geometrical optics is based on the eikonal equation, which can be understood as a limiting case of wave optics where the wavelength becomes infinitely short.
In Section 2, the influence of Earth's rotation on the propagation of light is estimated using geometrical optics. Section 3 is concerned with the derivation of a wave equation from Maxwell's equations in appropriate coordinates. An approximate form of the equation is solved in Section 5. In Section 6, it is shown that more accurate equations essentially lead to the same dispersion relation. Finally, in Section 7 numerical examples of a concrete setup are provided and the effects of Earth's rotation around its own axis and the rotation about the Sun are compared.

Model Assumptions
In order for wave optics to be applicable, we must specify the full geometry of the medium of interest. We consider cylindrical step index waveguides which consist of a cylindrical core of radius a and refractive index n 1 surrounded by a cladding of diameter a a and refractive index n 2 < n 1 . The dielectric will be assumed to be nonmagnetic and its field response will be assumed to be linear. The dependence of the fields D and H on E and B is described covariantly by means of the optical metric, introduced in [3] (not to be confused with the optical metric defined when studying null geodesics in stationary spacetimes).
We will neglect boundary effects which arise at the two ends of the waveguide and at the outer boundary of the cladding. In this sense, the waveguide will be treated as infinitely long and the cladding as infinitely thick.
To describe Earth's rotation we use rigidly rotating coordinates in flat spacetime. The metric tensor will be approximated under the assumption that the rotational velocity is sufficiently slow compared to the speed of light. Such a description applies both to Earth's rotation about its own axis and its orbit around the Sun (neglecting its orbital eccentricity).

Conventions
We use Heaviside-Lorentz units which, in our problem, coincide with Gaussian units, due to the absence of external charges and currents. Furthermore, unless explicitly specified otherwise, we use units where the speed of light in vacuum is equal to one. Compared to SI units, this means that c = µ 0 = ε 0 = 1 .
The Einstein summation convention is used. Greek indices will range from 0 to 4 and Latin indices range from 1 to 3.

Geometrical Optics
In this section we discuss the metric tensor used to describe rotating coordinate systems and introduce the optical metric to describe light propagation in rotating media. A first estimate on the correction to the dispersion relation is derived in the framework of geometrical optics.
For slowly rotating systems, we may neglect terms quadratic in v = Ωρ. We thus arrive at the following metric in "rotating Cartesian coordinates" (x µ ) = (t, x, y, z): where v i = ijk Ω j x k is the local velocity field and Ω is the angular velocity of the rotating system. In the chosen coordinates, the inverse metric tensor is given by In what follows, the components of any tensor will refer to the corotating coordinate system x µ . Furthermore, Ω| x|/c will be assumed to be sufficiently small, so that all terms quadratic in ρΩ can be neglected.
Note that wave equations (in particular those for the electromagnetic field) also depend on derivatives of the metric. By expanding the metric tensor to order v = Ωρ, the wave equation is not expanded in powers of v alone, but also in powers of LΩ, where L is any characteristic length scale of the problem: in our case these are the length and the radius a of the waveguide, as well as the wavelength λ. Thus, by neglecting terms of order O(v 2 ) in the components of the metric tensor, we also neglect terms of order O(a 2 Ω 2 ), O( 2 Ω 2 ) and O(λ 2 Ω 2 ) in the field equations. Since the typical length scales of the dielectric are much smaller than ρ (in our case, this is either Earth's radius or its distance to the Sun), these error terms are expected to be negligible in all practical applications. By slight abuse of notation, we will refer to all of these correction terms by O(Ω 2 ) and thus formally expand all equations in powers of the (dimensionful) angular velocity Ω .

The Optical Metric
The (contravariant) optical metric tensor γ of a linear dielectric with four-velocity vector field u is defined as where g is the spacetime metric and n is the refractive index. In the comoving coordinates the four-velocity of the co-rotating dielectric is given by Using the velocity vector field v = Ω × x, the components of γ (in the comoving coordinates) take the form where terms of order O(v 2 ) have been neglected.

Geometrical Optics Approximation
The eikonal equation may now be used to obtain a first estimate on the influence of Earth's rotation on the propagation of light. One could try to solve (10) in a cylindrical waveguide, using appropriate conditions at the core-cladding interface, which we have not attempted to do. Instead, to simplify the problem, we look for plane wave solutions where the dispersion relation (which relates k 2 and ω 2 ) in the absence of rotation reduces to k = nω, where n is the effective refractive index obtained by solving the full Maxwell's equations in the unperturbed case. Writing where m is the unit vector along the constant vector k, and using the explicit form (9) of the optical metric, the eikonal equation (10) with k = dψ = ωdt − k i dx i yields the quadratic equation where v = m · v. Expanding the positive solution in powers of v one obtains from which one obtains the first order relative correction which approximates unexpectedly well one of the main results of our work below. As stated before, we require the dispersion relation to reproduce β = nω in the unperturbed case, where n is the effective refractive index obtained from the exact solution of the unperturbed problem, which depends on the precise geometry of the waveguide. Comparing the wavelengths in two parallel interferometer arms at different heights, we find that the rotationally induced phaseshift in a Mach-Zehnder interferometer is where A is the enclosed area and ϑ is the angle between the direction of motion and the direction of light propagation. Comparing this with the gravitationally induced phase shift where g is the local gravitational acceleration (see Appendix B for a derivation of this formula), we find that the rotational effect exceeds the gravitational one by a factor of up to Ω/(ng) ≈ 1.6 × 10 3 (depending on the angle ϑ).
In the next sections, we will compute δβ using Maxwell's equations in the metric (5), which requires specialisation to a fixed geometry. We will consider an infinite cylindrical step index waveguide, for which an exact solution in the absence of rotation is known.

Maxwell's Equations in Rotating Coordinates
In this section, we formulate Maxwell's equations in the linearised Born metric (5) and derive a wave equation for the electric and magnetic fields.

Field Equations
In manifestly covariant formulations of electrodynamics, the fields E and B are subsumed by a two-form F , and the fields D and H are combined to form the bivectorF . In the absence of electromagnetic charges or currents, Maxwell's equations can be written in the form where d denotes the exterior derivative and δ is the codifferential. In local coordinates, these equations take the form where the square brackets indicates antisymmetrisation and ∇ denotes the Levi-Civita covariant derivative associated with the spacetime metric g. For linear media of permittivity ε and permeability µ, the optical metric γ allows to write the relation between F andF in the concise form which correctly reduces toF αβ = F αβ in vacuum. We now discuss a "3+1 split" of these equations, which allows us to derive a wave equation for the electric and magnetic fields.

Decomposition of the Field Strength Tensors
Following [2, eq. (3.5)] we decompose the two-form F and the bivectorF as where e and h are spatial one-forms and b and d (not to be confused with the differential operator d) are spatial vectors. More explicitly, this means where, as already pointed out, all field components F µν andF µν refer to the corotating coordinate system. Using (19), we may express the d and h fields in terms of the e and b fields as where terms of order O(v 2 ) have been neglected.

Decomposition of the Field Equations
From [2, eqns. (3.12) and (3.15)] we have the following decomposition of the field equationṡ We need not make any distinction between upper and lower indices since they are manipulated with the spatial metric δ ij . Since b and d have a vanishing spatial divergence, it is natural to work with these two fields and express e and h in terms of them. From (22) one finds to order O(v) where n 2 = µ was used. We note that for n = 1, we may set u := − v/(n 2 − 1) to obtain which coincides with the relations between e, b and d, h in a medium which moves with velocity u relative to an inertial system [5, p. 329, eqns. (76,10) and (76,11)]. This is due to the special form of the metric perturbation, which allows us to write where (u µ ) describes a dielectric at rest and (ũ µ ) describes slow movement of the velocity given above. However, this equivalence does not extend to higher powers of v. Substituting (24) into (23), we obtain the following field equations for d and b.
where [ v, w] denotes the commutator (Lie bracket) of the vector fields v and w.

Electric-Magnetic Symmetry
Note that the field equations are invariant under the following transformation among b and which induces the following transformation of the fields e and h:

The Wave Equation for the Electromagnetic Field
The field equations (27) imply the following wave equations for the electromagnetic field: A derivation of these equations is given in Appendix A.
To get some insight into the relative magnitudes of the two source terms, let us estimate spatial derivatives of the fields by the inverse wavelength 1/λ. Since v is of the order ΩR, we estimate that the first order differential operator is suppressed relative to the second order operator by the factor λ/R, which is clearly negligible. In 6.1, we show that this term does not modify the main result obtained without it.
Furthermore, let us note that the directional derivative along v depends on the position in the dielectric. Decomposing the radius vector x as x = R + r, where R points from centre of rotation to one end of the dielectric, and choosing z to be the distance along the symmetry axis of the waveguide, the coordinate ranges inside the core are Neglecting the dependence of v on x and y is expected to produce relative errors of order O(a/R), while those resulting from neglecting its z dependence are expected to be since the length of a waveguide is typically many orders of magnitude larger than its diameter. As a first approximation we will use v ≈ V , where and in section 6.2 we show that corrections of relative order /R have no influence on the dispersion relation. Instead of considering the full equations we may instead consider the simplified wave equations which yield the same dispersion relation up to terms of order a/R.

Maxwell's Equations in Cylindrical Waveguides
In this section, we give an outline of the calculations in the framework of wave optics and review the solution to the unperturbed problem.

Outline of the Calculation
We seek solutions to Maxwell's equations inside and outside an infinite cylinder of radius a, where all field components depend on time t and z (distance along the symmetry axis) as e i(ωt−βz) . All calculations are performed in corotating cylindrical coordinates (t, r, θ, z). We choose as "fundamental fields" which determine all other field components uniquely. Symbolically, we write F = F (f ; β, ω), where F is an abbreviation for the collection of fields ( E, D, B, H). Due to the linearity of the medium, F depends linearly on f , which satisfies a wave equation of the form where g is a linear differential operator, which vanishes in the unperturbed case. At the core-cladding interface (r = a), the field components D r , E θ , E z , B r , H θ and H z are required to be continuous. As only four of these conditions are linearly independent, it suffices to require the continuity of D r , E z , B r and H z , for which we symbolically write where · is a linear measure of the jump height. The problem at hand thus reduces to a wave equation for f on two disjoint domains (core and cladding), subject to the linear continuity constraint at the interface r = a.

Unperturbed Case
In the unperturbed case, where g vanishes, we consider the excitation of a single mode, i.e.
The two-component radial function h(r; β, ω) has four (complex) degrees of freedom, since it satisfies a homogeneous ordinary differential equation (ODE) of second order on two disjoint domains (core and cladding) and is required to be regular at the origin and to decay at large distances. We choose linear coordinates α on the solution space and write h = h α , such that the four linearly independent continuity conditions are equivalent to a matrix equation of the form To obtain a nontrivial solution (α = 0) one is led to the equation which provides the dispersion relation linking β and ω. Finally, the coefficients α are chosen to be in the kernel of M (0) . An explicit calculation shows that this kernel is one-dimensional.

First Order Corrections
To describe linear corrections due to Earth's rotation, we write where f (1) also depends on θ due to the angular dependence of g. Decomposing f (1) into a Fourier series in the angular variable, every coefficient function f (1) k satisfies an inhomogeneous ODE of second order, subject to the same constraints as above, and is thus an element of a four-dimensional (complex) solution space. Due to the linear independence of the angular modes, continuity conditions hold for all modes separately. Further analysis requires a distinction between the main mode (no further dependence on θ, so k = 0) and the sidebands (k = 0).
For the main mode, f can be chosen to be any particular solution to the ODE, since a homogeneous part can be removed by redefining α. Thus, the radial function of the main mode depends linearly on α, so the continuity condition of this mode is again equivalent to a homogeneous matrix equation which determines the first order correction to the dispersion relation. For the sidebands, the radial functions f (1) k satisfy inhomogeneous ODEs, whose inhomogeneities are linear in the coefficients α. Parameterising homogeneous solutions linearly with parameters χ k , f (1) k becomes a function linear in α and χ k , so the continuity condition is equivalent to an inhomogeneous matrix equation of the form whose solution χ k turns out to be unique.
In the next section, we shall review the unperturbed solution in light of this structure, and in later sections we will carry out the perturbational calculation, following the outline discussed here.

Unperturbed Solution
In this section, we review the solutions of the problem when Ω = 0, using the notation of [2] with φ = 0 (i.e. no gravitational potential) and c = µ 0 = 0 = 1.
Since the dielectric is assumed to be nonmagnetic, we have µ = 1 and thus n 2 = . The relation between D, H and E, B are thus The refractive index is assumed to be piecewise constant, namely where a denotes the radius of the dielectric. We separate the dependence of the fields on t, z and θ via In the absence of a gravitational potential (φ = 0), the constants U and W defined in [2, eqns. (4.17, 4.18)] reduce to with U and W positive. The function is the continuous solution of the equations which is bounded when r → 0 and has finite total energy. The z components of the fields are then given by where the constants A and B are to be determined from continuity conditions. The transverse components of the fields are obtained from the z components by where we have set Together with the relations (45), these equations constitute the map F (f ; β, ω), which was introduced in the previous section. Due to the special form of (45), the continuous fields E z and H z satisfy the same wave equation as the fundamental fields D z and B z . Thus, their continuity can be implemented immediately and the nontrivial continuity conditions reduce to two linearly independent constraints, which can be written in the following form [2, eq. (4.41)]: This is the anticipated matrix formulation of the continuity condition F = 0. The requirement of (54) to admit nontrivial solutions for A and B leads to the dispersion relation which must be solved numerically.

Generalisation
In light of the general formalism of the previous section, it will prove useful to generalise the expressions found above. To allow for discontinuous fundamental fields D z and B z , we define such that the solutions can be written as The original expressions are thus obtained by setting This form of the solution will be used in Chapter 5 to determine first order corrections due to Earth's rotation. In this calculation, it will prove useful to allow for parameters α = (α 1 , . . . , α 4 ), which deviate from the unperturbed values given here.

Main Corrections
In this section, we compute the first order corrections to the fields which arise from the approximate wave equation (35) and the modified relation between e, h and d, b given in (24).
To describe corrections to the unperturbed solution, we write

The Wave Equation for the z Component
In this section, we specialise the vectorial equations (35) to scalar wave equations for the z components of the fields d and b.
From now on, we will use cylindrical coordinates (r, θ, z) and refer all tensorial objects to the orthonormal frame whose coframe is given by Since the wave equation (35) is independent of the spatial coordinate system, we may evaluate it in the coordinate system at hand. As stated in the outline, we are mainly interested in the wave equation for the z components of the fields for which the equation simplifies since We thus obtain the scalar equations where V (b z ) denotes the action of the vector field V on the scalar function b z (i.e. the directional derivative). This is nothing but the scalar wave equation with respect to the optical metric (using the approximation v ≈ V ). Inserting the perturbative expansion (59), neglecting all terms of order O(Ω 2 ), and using the fact that D z and B z satisfy the homogeneous wave equation, we obtain

The Radial Equations
Summarising, we need to analyse the following differential equations for the z components of δ d and δ b: where the fields D z and B z are given by Here, the parameter β and the coefficients α 1 , . . . , α 4 do not necessarily coincide with the unperturbed parameters, which would have been obtained when Ω = 0. Note that in doing so, one modifies the right hand side of the wave equations by terms of order O(Ω 2 ), which we assume to be negligible. Note also that V on the right hand side depends on the angle θ via where V x , V y (and V z ) are constants. Thus, the θ-dependence of the source term is not given by e imθ alone; instead, there are sidebands with an angular dependence of the form e i(m±1)θ . Defining the constants V ± by and decomposing all fields as a Fourier series in θ of the form one has e.g.
where we have set For the perturbed fields δd z , δb z we use the ansatz which leads to the following set of equations for the radial functions: where ζ = n 2 ω 2 − β 2 as before and is either an ordinary or modified Bessel operator, depending on whether ζ > 0 or ζ < 0. We consider only the equations for d, since those for b are identical. Using (57), we obtain equations of the form where the coefficients in front of the source terms have been dropped momentarily -they can be added later due to the linearity of the equations. The powers of a were chosen such that the functions p 0 (α i , α j ; r) and p ± (α i , α j ; r) have the same dimension as f 0 (α i , α j ; r), which, in turn, has the same dimension as the coefficients α i , α j . The homogeneous solutions of interest (we require the functions to be finite at the origin and to decay sufficiently fast at large distances) to the first equation are given by (56). The corresponding homogeneous solutions for the sidebands are then given by The factors ∓U and −W were chosen such that which is due to the following recursion relations Particular solutions p 0 and p ± which also decay at large distances and do not blow up at r = 0 are given by where we have introduced the functions which can be evaluated explicitly using the indefinite integrals (up to additive constants) The solution to the equations (73) are thus The functions p 0 and p ± are plotted in figure 1. As expected, the functions are regular at the origin and decay rapidly in the cladding.

Compatibility of Approximations
Until now we were only concerned with the z components of the fields. In all calculations we have omitted terms quadratic in Ω and we have neglected some terms leading to relative errors of the order of a/R or /R, where a denotes the radius of the dielectric, its length and R denotes Earth's radius.
In the next section we will compute the remaining components using similar approximations. There is however a subtlety: as we will now show, these approximations cannot be discussed term by term, since the resulting expressions must satisfy two equations, which relate three field components each.
We shall later require the continuity of six field components at r = a (the core-cladding interface) due to the following conditions: 1. The non-existence of magnetic surface charges and surface currents is equivalent to the continuity of b r , e θ and e z ; and 2. The absence of electric surface charges and surface currents is equivalent to the continuity of d r , h θ and h z .
Clearly, the continuity conditions must hold for all Fourier modes separately. Note, however, that these conditions constitute six linear equations in the four coefficients which parameterise the solutions in each mode. In order to obtain nontrivial solutions, two continuity conditions must therefore be expressible in terms of the remaining four equations. The existence of such a relation follows in fact from the r-components of (23), One sees that, for solutions as considered in this work, if these two equations hold everywhere, and if the b r , d r , h z , e z are continuous, then h θ and e θ will also be continuous. Note that continuity of b r , d r , h θ , e θ does not imply continuity of h z and e z , since they may be independent of θ.
In order to obtain consistent equations, it is therefore necessary to make approximations which are compatible with these two equations.
The components e z , h z are determined by (24) whereas the components b r , d r , h θ , e θ will be computed from (27). The approximations made in either set of equations thus determines the approximations in the other one.

The Transverse Components of d and b
We use the equations (23) to express the "transverse" r and θ components of the fields in terms of the "longitudinal" z component. Substituting (59), one obtains to first order in Ω where we have used µ = 1, ε = n 2 , D = n 2 E, B = H and fact that n 2 is locally constant. If all field components depend on t and z via exp(i(ωt − βz)), the differential equations for the transverse components reduce to the following set of linear algebraic equations with the solution where ζ = n 2 ω 2 − β 2 . Note that for any vector field w it holds that where ∇ denotes the Levi-Civita covariant derivative with respect to the spatial metric. As a first approximation, we write The approximation ∇ v w ≈ ∇ V w corresponds to x ≈ R, where relative errors are expected to be of order a/R or /R. Moreover, neglecting Ω × w compared to the first term is expected to produce relative errors of the order of λ/R, cf. the discussion in Section 3.5.
Defining the connection one-forms we may decompose the covariant derivative as We find the only non-zero connection forms to be Since the first order perturbation fields have three θ modes, e.g.
we decompose all expressions in a Fourier series in θ using the notation introduced in (68). As before, we use the abbreviation X ± = X ±1 . Note that the cylindrical frame components of E θ , E r , B θ and B r of the unperturbed fields have vanishing m ± 1 components by (52). Using the notation (69) (and a similar notation for Ω), we obtain the following decomposition of the Lie bracket [ w, v], where w is any unperturbed field: where we have used the operator c ± m = ∂ r ∓ m/r, defined in (71). Combining (87) and (94), one obtains the following expressions for the m modes of the electric field d and the magnetic field h and similarly for the m ± 1 modes

The Field Components of e and h
We rewrite equation (24) with µ = 1 and ε = n 2 in the form where v = Ω × x. Substituting the perturbative expansion (59) and using D = n 2 E as well as B = H, one obtains to first order in Ω We compute the θ and z components of these fields, which are required to be continuous. Using the approximation v ≈ V (in accordance with x = R + r ≈ R), one finds for the m modes and for the m ± 1 modes (100)

Consistency of the Approximations
In this section, we verify that the approximations made are consistent with the equations (83)-(84). Due to the electromagnetic symmetry, it suffices to verify the equation Let us start by showing that if all fields depend on z only through e −iβz , and if the field components are related to each other as in the unperturbed case, i.e. the r and θ components are related to the z components by (52) (or equivalently by (87) without the Lie bracket terms) and by h = b, d = n 2 e, then (83)-(84) hold.
Hence, we consider the following subset of the exact Maxwell's equations with Ω ≡ 0: with ζ defined in (53), where all fields depend on z only as e −iβz . We have Equation (101) with the d and h fields replaced by D and H now follows from the last equation in (102), namely H z = B z . In fact, the calculation in (103) shows that if a set of equations with a structure as in (102) holds, then (101) will hold regardless of the precise form of the dependence upon r and θ of the z components of the fields.
Let us return to the problem at hand. As such, (87) can be written as while (98) reads The calculation in (103) applies to those terms in (104)-(105) which have been written out explicitly. Hence, to verify (101) it remains to check that an identity of the form (101) is satisfied by the terms, denoted by ... in (104)-(105), which have not been written in explicit form.
Equivalently, it remains to consider the deviations from the "unperturbed relations" when checking (101). Let us use the notation x y when x (e.g. δd r ) deviates from the "unperturbed expression" by y. For instance, for the m mode we have by equations (95) and (99), suppressing the common factor −V z temporarily Since ζ = n 2 ω 2 − β 2 it follows that and thus which establishes the consistency of the equations for the m mode. For the sidebands, we have from (96) and (100), and after omitting the common factor −V ± , Using the definition c ± m = ∂ r ∓ m/r, one obtains Thus, (111) This completes the proof of the consistency of the approximate equations for the sidebands.

Continuity Conditions
In this section, we discuss a general method to analyse the continuity conditions of the various field components.
According to the discussion in Section 5.3, the requirement that d r , e θ , e z and b r , h θ , h z be continuous at the core-cladding interface constitute only four linearly independent equations. It thus suffices to require d r , e z , b r and h z to be continuous, so we define a vector which summarises the relevant jump heights: where x T denotes the transpose of x. Here, we have introduced the notation f for the jump height of a function f at r = a:

Discontinuity Matrices
Due to linearity of the system, F is a linear homogeneous function of the coefficients which parameterise the fields according to (57) and (82). Thus, there exist (complex) 4 × 4 matrices M 1 , M ± 2 such that We will refer to these matrices as discontinuity matrices. In fact, equations of this form hold for every Fourier mode of the fields separately.
To compute these matrices, we consider the jump height of a field component w as a linear map of the coefficients (e.g. α) to the field of complex numbers. We may thus identify w with a row vector of four complex entries, such that for which we write w .
The subscript α indicates that the form acts on the coefficients α, and not on the coefficients χ ± . We will omit this subscript whenever there is no danger of confusion. Using this notation, the rows of M 1 in the above equation are exactly the row vectors which represent d r , e z , b r and h z , respectively. Using the expressions (56) and (57), we find jump heights of the unperturbed fields to be given by Since the unperturbed fields depend on α 1 , . . . , α 4 only, it is clear that the suppressed subscript is α.
Multiplication of a field with a constant clearly multiplies the corresponding row vector by the same number. If a field component w, whose jump heights are given by is multiplied by a function f , whose left-and right-sided limits exist, then the jump height of f w is given by For example, if f = n 2 or f = ζ, then n 2 w . = n 2 1 w 1 n 2 2 w 2 n 2 1 w 3 n 2 2 w 4 , Similar equations hold for arbitrary powers of n and ζ.

Symmetries of the discontinuity matrices
Due to the electromagnetic symmetry, all discontinuity matrices are of the form Indeed, if the jump height of any component of the d field is given by we may apply the symmetry transformation d → b and b → − d/n 2 to obtain Similarly, if the jump height of e i is given by we may apply the same transformation as above, which entails e → h/n 2 , and find It thus suffices to compute the jump heights of d r and e z , since the first two rows determine the entire matrix.

General Structure of Discontinuity matrices
In general, there are two contributions to these matrices. Some of the terms stem from the z components of the fields, which produce terms in all other field components; while other terms stem from the cross product in (24) which relates d and b to e and h. Let f be any function which enters the expressions for the z components of the fields as We denote by D(f ) the corresponding matrix defined by the condition that the such produced terms in F are This method allows to describe all contributions to the discontinuity matrices which arise due to the z components of the fields alone. To cover the remaining terms which arise from the cross products in (24), we will introduce additional matrices D 1 , D 2 etc. Since, at our level of approximation, the cross product terms are determined by the unperturbed fields D and B, which are parameterised by α alone, these additional matrices will multiply α only.

Continuity of the Main Mode
According to (59), (57) and (82), the main Fourier modes of the z components of d and b are given by Correspondingly, we write the jump height of the fields in the form where D(f 0 ) describes the jump heights due to the functions f 0 alone, D(p 0 ) those due to the presence of the function p 0 , and D 1 the terms which arise form the cross product terms in (24). The factor 2a 2 βωV z was factored from D 1 for later convenience. We start by showing that the matrix D(f 0 ) is given by where all functions whose arguments have been suppressed are evaluated at the core-cladding interface r = a. Thus, all ordinary Bessel functions are evaluated at U and all modified Bessel functions are evaluated at W . Similarly, the functions ∆ ν , T ν with suppressed arguments will be evaluated at r/a = 1.
To compute this matrix, it suffices to consider the equations for the unperturbed fields, since the effects of all perturbations are described by D(p 0 ) and D 1 . From (52) we have Using (118), one obtains and thus D r .
which is exactly the first row of the matrix given in (132). To compute the jump height of E z , we use By virtue of (118), we find which is the second row of D(f 0 ). The third an fourth row are now determined by the electromagnetic symmetry, i.e. by the general form (123).
To compute D(p 0 ), we note that D(f ) is uniquely determined by the limiting behaviour of f and its first derivative at r = a. Comparing the limiting behaviour of p 0 and its first derivatives with the limiting behaviour of f 0 , we find that the matrix D(p 0 ) is obtained from We thus obtain Finally, we compute D 1 explicitly. According to (95), the relevant term for δd z 0 is From (52) we obtain where we have setζ Here, we have introduced the constantsŪ andW according tō Having established that we may use (118) to obtain from which it follows that This determines the first row of D 1 . According to (99) we have n 2 δe z 0 = δd z 0 , so V does not contribute to δe z 0 . Accordingly, the second row of D 1 vanishes. The remaining rows are then then fully determined by the general form (123), so we arrive at the result The continuity conditions of the main mode now take the concise form

Dispersion Relation
We have established that the continuity condition of the main mode takes the form where In the unperturbed case, where Ω = 0, M reduces to M (0) , so δM is the first order perturbation of M due to Earth's rotation. Generically, the kernel of the matrix M is trivial, so we must arrange parameters such that its determinant vanishes, if we wish to obtain a nontrivial solution. Since the only remaining parameters are β and ω, one is determined by the other via the dispersion relation We choose to take the frequency ω as the independent variable, such that the dispersion relation determines β = β(ω). Equation (151) is of major importance since encodes how the wave vector β changes due to Earth's rotation. As in the unperturbed case, this equation can only be solved numerically. If a solution to the unperturbed dispersion relation is known, one may obtain the first order correction to β by expanding the dispersion relation (151) to first order in the velocity V .
Here, we see the importance of letting the parameters α deviate from their unperturbed values α (0) . Had we added homogeneous solutions to the wave equations with parameters α (1) and neglected the term δM α (1) , which is quadratic in the velocity V , instead of (149) we would have obtained since M (0) α (0) = 0. This equation already suggests that the dispersion relation must be modified since it cannot be solved for α (1) due to M being singular when β assumes its unperturbed value. Further analysis would have been necessary to obtain the dispersion relation (151). Linearising the condition det M = 0 in the vicinity of V = 0 using Jacobi's formula Since these matrices are nonsingular (in fact det A 1 = 1 and det A 2 = n 2 1 n 2 2 ) this does not change the roots of det M = 0.
It proves useful to introduce the abbreviations as well as the brackets M (0) then takes the concise form from which one can read off the determinant Since similar structures will appear again in the following calculations, we set such that the defining equation of the unperturbed value of β takes the form Using the fact that as well as the Wronskian the adjugate ofM (0) and the perturbations δM 1 , δM 2 take the form From these matrices, one obtains the following expressions for their traces withM (0) : We thus arrive at the following result for the linear change in β: where This formula is the main result of this section: it determines the linearised deviation from the unperturbed dispersion relation (162). If β solves the unperturbed equation for given ω, (167) determines the first order correction to β for constant ω.
Once the perturbed value of β is known, the coefficients α (which parameterise the main mode) are obtained by finding the kernel of the matrix cf. (148) and (149). From the unperturbed problem it is known that the kernel of M (0) is one-dimensional. Since the dimension of the kernel is an upper semi-continuous function, the same will hold for M , provided that the velocity V is sufficiently small. Note that the determination of α to order O(V ) requires M (0) to be computed using the perturbed value of β. Since δM is already of order O(V ) one may use the unperturbed value of β in this matrix.

Continuity of the Sidebands
In this section, we compute the discontinuity matrices associated to the sidebands. From (82) it is clear that the jump height of the sidebands takes the form The coefficient in front of D ± 2 , which contains the contributions from the usual cross product terms, was factored for convenience.
We may obtain D(f ± ) from D(f 0 ) by noting that is related to f 0 (α 1 , α 2 ; r) by the substitution Applying this to D(f 0 ) and replacing all factors m (which are due to derivatives with respect to the angle θ) by m ± 1, we find Similarly, comparing p 0 with p ± in (79), find that D(p ± ) is related to D(p 0 ) by the substitution so we find where we have set Finally, we compute the matrix D ± 2 . To this end, we note that the relevant terms in the equation for δd r ± are Using (141) and the analogous equation we obtain Simplifying this expression using the identities one obtains Using the explicit expressions for D z and B z , as well as the recursion relations for the Bessel functions (78), we find which implies This determines the first row of D ± 2 .
Next, we compute the second row, which is (up to the factor −iV ± ) given by the jump height of δe z ± . From (100) we find the relevant term to be Substituting this into the above expression for δe z ± and computing the jump height, we find which implies Using the general form (123), we thus arrive at the result Having computed all matrices which enter (170) and having computed the coefficients α of the main mode, the parameters χ ± are determined by Here and in the following sections, we shall assume that the matrices D(f ± ) are invertible.
Recall that the unperturbed dispersion relation is such that det D(f 0 ) = 0. The special case where both D(f 0 ) and either one of the two matrices D(f ± ) are singular for the same value of β would require separate analysis.

Further Corrections
In this section, we briefly discuss additional corrections of order Ω and λΩ, where is the length of the waveguide and λ the light wavelength. Corrections of order λΩ arise from the first order differential operator in equation (34), and corrections of order Ω are obtained by improving the approximation x = R + r ≈ R to x ≈ R + z e z . By adding these two corrections, the only neglected contribution of order Ω is of order O(Ωa), where a is the core radius of the dielectric.
Due to their small amplitudes, we expect that the sidebands will not be experimentally accessible in the near future and thus restrict our attention to the main mode. We show that this main mode is unaffected by the above mentioned effects, i.e. the contributions of order Ω and λΩ vanish.

Corrections of Order λΩ
In the previous analysis, the source term proportional to Ω in (34) was neglected. In this section, we briefly discuss the influence of this term and thus consider the equations The additional terms do not modify the main mode, since the m-th Fourier component vanishes: Thus, this correction does not modify our result for the shift in β but merely modifies the coefficients χ ± 1 , . . . , χ ± 4 , which parameterise the sidebands according to (82).

Corrections of Order Ω
Hitherto, we have approximated to obtain the simplified wave equation (35) and to compute the transverse field components, e.g. in (95) and (99). This approximation is expected to produce relative errors of the order of a/R or /R, where a is the radius of the dielectric, its length and R is Earth's radius. Instead of the above approximation, we thus write The corresponding approximation for the velocity vector field v reads where the relative velocity u = Ω × (z e r ) has the following components when referred to a cylindrical orthonormal frame: But since u depends only on Ω r and Ω θ , but not on Ω z , it depends on the angle θ by sin θ or cos θ, but it does not have a constant θ-independent component. Hence, all terms arising from u only affect the sidebands and leave the main mode unchanged. We thus conclude that the entire analysis of the continuity conditions for the main mode is unaltered by the terms introduced here.

Numerical Examples
We apply our formulae to the same setup as considered in [2]. Table 1 lists both the parameters of the waveguide (refractive indices and core diameter) and the parameters of the light wave (wavelength and angular mode number). The effective refractive index n was found by solving (162) numerically, which yields

Corrections to the Dispersion Relation
The results derived in this work may be applied both to Earth's rotation about its own axis (spin) and its orbital rotation about the Sun, if the eccentricity of Earth's orbit is neglected.
In the first case, we use R = R E ≈ 6371 km and Ω = 2π /d. Since Earth is slightly oblate, more accurate results could be obtained by setting R to be the local distance to Earth's axis. To illustrate the dependence of the effect on the latitude ψ, we compute δβ both at the equator (ψ = π/2) and in Vienna (ψ = 48 • 12 30 ). In both cases, we assume the waveguide to be aligned with the velocity vector due to the rotation, i.e. the z axis points east.
To estimate the effect of Earth's orbital motion, we set Ω = 2π /yr and R = 1 ua. Again, we choose to align the waveguide such that the z axis points along the momentary direction of motion. Note however, that if the alignment relative to Earth is kept fixed, then V z varies over time due to Earth's spin, so δβ will oscillate with a period of one day.
While these numbers quantify the anisotropy in the propagation of light, one may also compare the wavelength of two parallel light rays propagating at different heights, e.g. in a Mach-Zehnder type interferometer whose two arms are displaced by a height difference ∆h. In view of the high quality of the geometrical optics approximation, we discuss this setup using the formula (14). If ϑ denotes the angle between the direction of wave propagation and the local velocity vector field v, the geometrical optics formula implies a difference in wavelengths in the two arms of ∆β = ωΩ∆h cos ϑ .
If denotes the length of the interferometer arms, the phase difference at the recombination point evaluates to ∆Φ rot. = ωA cos ϑ , where A = ∆h is the enclosed area.
We compare this effect with the shift of β caused by Earth's gravitational field. In [2] it was found that a vertical displacement of the dielectric by a heigh δh results in the following first order change in the wavelength β: Since −2g/c 2 ≈ 2.18 × 10 −16 m −1 , this is in good agreement with the following formula, which we derive in Appendix B: where g is the local gravitational acceleration. Even for vertical distances of the order of 100 m, the anisotropy effect due to Earth's spin in Vienna is seven orders of magnitude larger. This suggests that second order effects due to Earth's spin may still be comparable with the first order effect due to Earth's gravitational field. Furthermore, the third order effect due to Earth's orbital motion might still be larger than the first order gravitational effect. At such a level of accuracy, additional effects stemming from the eccentricity of the orbit might become relevant as well.
In the same Mach-Zehnder setup as described above, the phase shift due to the gravitational potential is found to be given by see e.g. [4]. Comparing with the previous formula and restoring the correct powers of c, one finds obtains the ratio ∆Φ rot. ∆Φ grav. = Ωc ng cos ϑ .
In the considered setup, the maximal ratio of the two phase shifts evaluates to Ωc/(ng) ≈ 1.6 × 10 3 , so the rotational effect is three orders of magnitude larger than the gravitational one.
Let us also note that Earth's spin causes a height difference of ∆h = 2R E in the gravitational field of the Sun twice a day. Since the Sun has a local acceleration of roughly g S ≈ 5.93 × 10 −3 m/s 2 , which -by the above formula -corresponds to a relative shift in β of −1.68 × 10 −12 . Even though this effect is much larger than the gravitational effect discussed in [2], it is unlikely to be directly observable since the effect seems to effect light in all arms of an interferometer in the same way.

Main Mode Amplitudes
Having solved for β in terms of ω, the coefficients α must be in the kernel of the continuity matrix for the main mode, as stated in (148). The kernel was computed numerically and the spanning vector α was normalised to α 1 = 1.
Unperturbed Spin (Equator) Spin (Vienna) Orbit  Table 3: Comparison of the main mode coefficients α in the unperturbed and perturbed case. The linearity of the equations was used to normalise the coefficients sucht that α 1 = 1. "Spin" refers to Earth's intrinsic rotation and "Orbit" to the rotation around the Sun. Table 3 shows the numerical results for the coefficients α (normalised to α 1 = 1). As before, the alignment was chosen such that the direction of wave propagation is parallel to the local velocity vector V . We note that the perturbation changes only the modulus of the coefficients, but not their phase.
Comparing, for example, the various results for α 4 , we find a relative change of 3.7% due to Earth's orbit around the Sun and a change of only 0.06% due to its intrinsic rotation at the equator. This ratio is due to the fact that Earth's orbital velocity is roughly 60 times larger than the tangential velocity at the equator.

Conclusion
We have computed the first order corrections to electromagnetic waves in cylindrical waveguides due to Earth's rotation by solving Maxwell's equations in the linearised Born metric. The correction to the dispersion relation was found to be in good agreement with the geometrical optics formula where n is the effective refractive index and v is the projection of the velocity vector v on the direction of wave propagation. The relative change of the wave vector β due to Earth's rotation about the Sun was found to be roughly 7 × 10 −5 , while the effect of Earth's rotation about its own axis causes relative changes of up to 1 × 10 −6 at the equator (the effect depends on the latitude and on the orientation of the waveguide in space). These corrections (relative to the unperturbed value) are many orders of magnitude larger than those due to Earth's gravitational field. Furthermore, considering a Mach-Zehnder type interferometer whose arms are placed at different heights, it was shown that the phase shift due to rotation isdepending on the orientation -up to three orders of magnitude larger than the one due to Earth's gravitational field.
It was found that Earth's rotation couples neighbouring modes of different angular dependence. The relative amplitude of the sidebands caused by Earth's intrinsic rotation was found to be of the order of 10 −5 , whereas Earth's orbital rotation around the Sun induces relative amplitudes of the order of 10 −3 .
We expect that higher order contributions to the dispersion relation cannot be computed using geometrical optics alone, since this method only approximates the result obtained here with an error which we estimate to be much larger than the second order correction.

A Derivation of the Wave Equation
In this section, we derive a wave equation for the electromagnetic field, taking all terms of order O(Ω) into account.
We start by writing (27) as which is equivalent to the original form since Differentiating the first equation in (205) with respect to time and multiplying by n 2 , one obtains Due to the second equation of (205) and n 2 = µ, the second term takes the form Since ∇ · b = 0, the double curl reduces to where ∆ b denotes the vector Laplacian of b. Moreover, due to ∇ · d = 0 and v = Ω × r, we have Together with the intermediate result which is also due to the fact that ∇ · d = 0, we obtain Combining this with we arrive at the wave equation

B Remarks on the Influence of Earth's Gravitational Field
We give a derivation of the first order change in the wave vector due to Earth's gravitational field using geometrical optics. We start from the post-Newtonian metric in the form where φ is Newton's potential. To first order in φ, the components of the inverse metric tensor read g µν = η µν + 2φ δ µν .
A medium at rest in the considered coordinates has four velocity (u µ ) = (u 0 , 0), where u 0 is determined from g(u, u) = −1, which yields Inserting this into (9), one obtains the components of the inverse optical metric The condition of k = ωdt − βdz being null with respect to γ then yields the quadratic equation (1 + 2φ)β 2 − n 2 (1 − 2φ)ω 2 = 0 , whose solution for positive β (ω > 0) is The linearised change in β due to the gravitational field is thus or, using β = nω + O(φ) and δφ = gδh, where g is the local gravitational acceleration and δh is a difference in height δβ β = −2gδh .
This is the same result as equation (1.1) in [2].