The Response of Laser Interferometric Gravitational Wave Detectors Beyond the Eikonal Equation

The response of Michelson interferometers to weak plane gravitational waves is computed at one order of accuracy beyond the eikonal equation. The modulation of the electromagnetic field amplitude and polarisation are taken into account by solving the transport equations of geometrical optics with boundary conditions adapted to laser interferometry. Considering both DC and balanced homodyne readout schemes, explicit formulae for the interferometer output signals are derived. These signals comprise perturbations of the optical path length, frequency and amplitude, and are shown to be insensitive to polarisation perturbations.


Introduction
Laser interferometric detectors of gravitational waves (gw's) have been studied extensively using a variety of methods. The perhaps simplest models are based on the Jacobi equation [1; 2], assuming the mirrors to follow nearby geodesics whose separation is read out interferometrically.
The assumption of small arm lengths can be removed by computing round-trip times of light rays, which are commonly modelled either as null curves with prescribed spatial trajectories [3][4][5][6][7] or as null geodesics [8][9][10][11]. A comparison of these three methods can be found in Ref. [12], where the underlying assumptions are assessed and criteria are formulated under which circumstances these methods yield the same result.
While the ray round-trip time determines the optical phase via the eikonal equation, it does not provide information about gw effects on the light amplitude [13] and polarisation [14; 15], the description of which requires methods beyond the eikonal equation, i.e. the Einstein-Maxwell equations or, in the high-frequency limit, the transport equations of geometrical optics. The influence of gravity on electromagnetic (em) polarisation has become known as the gravitational Faraday effect and was studied extensively in various scenarios, i.a. for light coming from distant stars and being perturbed by the presence of massive objects [16][17][18] or by gravitational waves, see e.g. Refs. [19][20][21] and references therein. A recent analysis quantifying the gw effect on light polarisation using Stokes parameters can be found in Ref. [22]. Naturally, such calculations were mainly concerned with one-way light propagation and used emission conditions different from those suitable to describe laser interferometers.
To the best of our knowledge, the first description of laser interferometry beyond the eikonal equation was given in Ref. [23] to compute the gw response of standing em waves in Fabry-Pérot cavities. The analysis, however, was limited to special alignments, where the electromagnetic and gravitational waves propagate either in parallel or orthogonal directions. The general case of arbitrary alignment was discussed much later in Ref. [24]. This calculation was criticised in Ref. [25] for not implementing appropriate boundary conditions, but the alternative analysis presented there was again limited to a special alignment, where the gw propagates orthogonally to both interferometer arms. General alignments were reconsidered in Ref. [26], but the analysis of interferometers provided there was restricted again to the phase shift. Although giving explicit expressions for the full em field, the authors did not describe polarisation reflection at mirrors but instead considered other setups to determine the polarisation perturbation. The analysis was then extended to describe the em field in a cavity [27], but the applications were restricted to specific alignments and polarisations of the gravitational wave.
Subsequent analyses of gw effects beyond the eikonal equation were provided for single light rays [28], optical cavities [29] (where the problem is typically reduced to a single dimension for simplicity), and Ref. [30] provided examples of polarisation effects in Kerr-Schild gw metrics.
However, we are not aware of comprehensive analyses of laser interferometers taking all these effects into account and allowing for arbitrary alignments relative to the gw.
The aim of this paper is to provide a systematic description of Michelson interferometers in the presence of weak gw's, based on the transport equations of geometrical optics. (A similar analysis based on the full Einstein-Maxwell equations is in preparation.) To this end, we model the emission of laser radiation, light propagation in interferometer arms, reflection at mirrors and the behaviour of light at beam splitters -all in the presence of a weak plane gravitational wave. In particular, light emission is modelled by prescribing the field on a suitably chosen timelike hypersurface, where we impose the em field to have a definite frequency and constant intensity, but allow for general perturbations of the em polarisation. The reflection is assumed to take place at perfect mirrors, and the partial reflection at the beam splitter is described using a simple model which is motivated by the description of light reflection at perfectly reflecting surfaces. The resulting expression for the light intensity at the detector, which is derived from the energy-momentum tensor, thus contains gw effects on the optical phase, amplitude and polarisation. The calculations presented here allow for arbitrary incidence angles, polarisations and waveforms of the gravitational wave and do not involve unidimensional approximations.
The structure of this paper is as follows. In Section 2 we describe the setup of simple Michelson interferometers used for gravitational wave detection, and in Section 3 we review the fundamental equations of geometrical optics. The emission of monochromatic plane waves by a laser source and the propagation of such radiation in the gw background is described in Section 4. Reflection of such radiation by perfect mirrors is modelled in Section 5, leading to a description of light rays as they travel back and forth in an interferometer arm. In Section 6 we give a model for the behaviour of the fields at (non-polarising) beam splitters, allowing to compute the radiation sent into the two arms from the laser source, and to compute the radiation reaching the detector via the output port. The results are combined in Section 7, where the output intensity of the interferometer is computed and discussed. Finally, approximate formulae for the interferometer output are derived and plotted in Section 8.
After finalising this paper, we became aware of the recent preprint [31], where the same setup is analysed.

Setup: Michelson Interferometers and Plane Gravitational Waves
Consider a linearised plane gravitational wave of amplitude ε 1 in transverse, traceless, and synchronous gauge g µν = η µν + εh µν (κ ρ x ρ ) , (2.1) where η is the background Minkowski metric (with the sign convention that spacelike vectors have positive norm), and κ is the gw wave vector where ω g is the gw frequency and n i is a unit vector in the unperturbed geometry. Note that we take κ to be past pointing so that for a wave propagating along the x-axis, the argument κ ρ x ρ coincides (up to the factor ω g ) with the standard null coordinate u = t − x. The gauge conditions then require i.e. the metric perturbation is purely spatial, transverse and traceless (where indices of h are raised and lowered with the background Minkowski metric).
This paper uses the following notation. For any two vectors v and w we write v.w ≡ η µν v µ w ν , and similarly we write the argument of the metric perturbation as κ.x ≡ κ µ x µ . Further, for quadratic forms such as h and g we use the notation h(v, w) = h µν v µ w ν .
Let us now consider Michelson interferometers in the metric (2.1). Here, we follow the standard practice of describing the geometry of the interferometer in terms of the (unperturbed) background metric. Physically, this can be justified by assuming the apparatus to be aligned before the gw arrives, and assuming all material points of the interferometer to follow geodesic motion as the gw passes by. For extended bodies, such as the beam splitter or the mirrors, deviations from this model could be described using elasticity theory, but this is beyond the scope of this paper. Note, however, that such arising corrections to the boundary conditions will be of order ε and can thus be described using the unperturbed equations of geometrical optics. Ref. [32,Sect. 5], but here we restrict ourselves to this simplified setup.) In typical applications for gravitational wave detection, the distances between the laser, beam splitter, and detector are much shorter than the two arm lengths of the interferometer. Hence, we will neglect gw effects in the segments (laser → beam splitter) and (beam splitter → detector), and focus on the propagation of light in the two interferometer arms.

Review of Geometrical Optics
This section briefly reviews the transport equations of geometrical optics in a notation useful for the following analysis. The equations presented here coincide with those of Ref. [30], see also Ref. [33] for the relation to other formalisms based on 3+1 decompositions such as those of Refs. [34; 35].
Geometrical optics at leading order 1 describes the electromagnetic field strength in the form (3.1) or sums of such expressions, where the eikonal ψ satisfies the eikonal equation and the amplitude f (a two-form) satisfies A motivation for these equations and the extension to higher order expansions is sketched in Appendix A. Writing the amplitude as where A satisfies the scalar amplitude transport equation the equations for the em polarisation tensor f µν reduce to so that f µν is parallel transported along the rays (integral curves of ∇ψ). This system of equations is consistent because ∇ψ satisfies the geodesic equation. The first two equations here entail that f can be written as where the em polarisation one-form E is orthogonal to dψ (and one is always free to add any multiple of dψ). The factor 1/ω was introduced here so that f has the same units as E. As made precise in Appendix B, the field E encodes the polarisation of the electromagnetic wave.
Solutions to the system (3.6) are readily constructed by choosing E to be parallel transported along the rays, i.e. ∇ ∇ψ E = 0, which is consistent with the condition E µ ∇ µ ψ = 0. Since the norm of E is preserved under parallel transport, the amplitude of the field is fully characterised by A , so that one may assume E to be normalised according to g(E, E) = 1. (E is necessarily spacelike, since it is orthogonal to the lightlike vector ∇ψ and linearly independent of it.) We are thus led to the system of equations Because ∇ψ is geodetic, the last two conditions propagate, i.e. if they hold at any point of a light ray, they hold along the entire ray.
The geometrical optics description of the field (at leading order) is thus obtained by first solving the eikonal equation (3.2), and then determining the amplitude A from (3.5), and the polarisation E from (3.8). In all cases, one must prescribe boundary values on a non-characteristic hypersurface Σ. In the next section, we prescribe such boundary data as to model the emission of light by a monochromatic source and solve for the eikonal, amplitude and polarisation.

Emission of Monochromatic Plane Waves
To model the emission of light from a radiating surface, we prescribe boundary data on a timelike hypersurface Σ, for which we take as the simplest model a coordinate-hyperplane Here, the m i 's are constants which may be taken to be normalised with respect to the unperturbed spatial metric, i.e. δ ij m i m j = 1, so that m.m = 1. The unit conormal (normal one-form) to this surface is given by For the eikonal, we take ψ = −ωt on Σ (corresponding to coherent emission with frequency ω).
This entails that on Σ we have dψ = −ωdt + αν, where α is determined by the eikonal equation up to a sign: α = ±ω. We shall assume ν(ψ)| Σ > 0, corresponding to emission along ν (instead of emission in the opposite direction).
For the amplitude A and the polarisation E, however, one cannot directly carry over the unperturbed boundary data to the perturbed problem: as we will see, the boundary data on Σ must satisfy a set of equations (specifically, the algebraic conditions of Eq. (3.8) and certain constraint equations discussed below) which differ from those of flat space, so that perturbations of the boundary data are unavoidable.
The algebraic conditions state that the prescribed polarisation, E, on Σ must be orthogonal to dψ = ω(ν − dt), and normalised (both in the perturbed metric), as dictated by Eq. (3.8).
Without loss of generality, E can be chosen to be purely spatial on Σ (this can be accomplished by adding to E a suitable multiple of dψ, cf. Appendix B), so that the algebraic conditions reduce to Seeking solutions of the form E = E (0) + O(ε), one obtains the perturbed polarisation on Σ where The first correction term in Eq. (4.4) ensures that E is orthogonal to ν in the perturbed geometry, the second correction term fixes the normalisation, and the last term allows for small variations of the emitted polarisation.
Additionally, as shown in Appendix C, the emitted polarisation must satisfy certain differential constraints, analogous to the divergence constraints in the usual spacelike Cauchy problem.
Denoting byẼ the pull-back of the overall field A E to Σ at any instant of time, they read where A, B are spatial indices associated to the surface Σ (at any instant of time),∇ is the induced Levi-Civita connection, and˜ denotes the induced volume form.
In flat space, these equations imply that the componentsẼ A are harmonic functions, so that the requirement that the field be bounded implies thatẼ (0) is spatially constant on Σ. Since one is free to prescribe the time dependence ofẼ (0) , the solutions are parameterised by two real functions of time alone. Plane polarised waves of constant energy density then correspond to A (0) = 1 and E (0) being a constant spatial vector of unit norm, tangent to Σ.
In the perturbed case, the solutions of Eq. (4.6) are still parameterised by two real functions of time t alone, but contrary to the unperturbed case it is not possible in general to arrange forẼ to have constant norm everywhere on Σ (see Appendix C). However, it turns out that the precise details of the perturbations are not relevant for the applications considered here, as one can nonetheless imposeẼ AẼ A = 1 to hold at the spatial coordinate origin with significant deviations arising only at length scales comparable to the gw wavelength. Hence, as the field is only read out on much smaller regions, such perturbations will not contribute to the final output signal.
Consistent boundary data can thus be obtained by construction a solutionẼ A of the form A to Eq. (4.6) (using the methods described in Appendix C) and separating amplitude and polarisation information by writing the perturbation in the form whereẼ is orthogonal toẼ (0) in the background metric. δA then corresponds to the amplitude perturbation, and extending the formẼ trivially to a space-time object by defining the temporal and normal components to be zero, one obtains a polarisation of the form (4.4), which satisfies all constraints.
We stress that δA ,Ẽ, and thus also E are defined on Σ only (the propagation of the fields away from this surface will be computed later). To emphasise this point and to simplify the notation in the following calculations, we write and similarly forẼ and E , where τ, ξ, ζ and m form an orthonormal basis in the unperturbed geometry (for example, one could take τ = dt and (m, ξ, ζ) an orthonormal basis of R 3 ).
With this notation, the boundary conditions for the amplitude and polarisation can be written as where we find it useful to decompose the polarisation perturbation into normal and tangential parts according to As mentioned above, the precise details of the functions δA and E are not important for the analysis here: the only assumption on the amplitude is so that the field at the emission point has unit energy density, and the assumptions on E are to be found in Eq. (4.5).
To summarise, we consider the equations of geometrical optics with the following data prescribed on the surface Σ = {m.x = 0}:

Eikonal Equation
In the unperturbed case, the eikonal is For the perturbed eikonal, the ansatz ψ = k.x + εψ (1) leads to which is solved by integrating the right-hand side along the unperturbed rays. To this end, we write x µ = x µ 0 + sk µ , where x µ 0 lies on the emission surface Σ. Contracting with m i , one finds s = m.x/ω, so that the emission point (for any given x away from Σ) is seen to be Note that by construction one has k ν ∂ ν x µ 0 = 0. Defining the integrated waveform one obtains the eikonal This clearly satisfies the boundary conditions: since κ.x 0 coincides with κ.x on Σ, the first-order correction vanishes there.
Taking the derivative, the local wave vector is found to be where ν is as in Eq. (4.2) (but evaluated away from the emission surface Σ), and where we have Amplitude Transport Equation Consider now the transport equation (3.5) for the scalar amplitude A . In the absence of a gravitational wave (unperturbed case), a plane wave has A = 1, so for the perturbed case we write A = 1 + εA (1) and obtain at first order Since we are working in harmonic coordinates, the scalar wave operator is simply = g µν ∂ µ ∂ ν .
Because the first term in Eq. (4.19) has vanishing second derivatives, only the second term contributes, where the unperturbed expression for the wave operator suffices.
where h ρσ (u) = ∂ u h ρσ (u). Consequently, the transport equation for A (1) takes the form Since k µ ∂ µ (κ.x 0 ) = 0, this leads to an amplitude which grows linearly with the distance from the emission surface. Indeed, imposing the boundary condition (4.23), the solution is The second term here, δA (x 0 ), corresponds to the boundary value perturbations of Eq. (4.23), evaluated at the emission point x 0 (x), as defined in Eq. (4.17). Due to the linear growth of the first term, we expect the leading order geometrical optics approximation to be only applicable for short distances from the emission surface Σ.

Polarisation Transport Equation
In the unperturbed case, E µ is a constant vector on Σ with i.e. the emitted polarisation is purely spatial and orthogonal to m. The transport equation then µ is constant everywhere (also away from Σ). In the perturbed case, we write E = E (0) + εE (1) and impose on E the boundary values given in Eq. (4.14). At first order, the equations for the polarisation are thus where Γ (1)ν µρ denotes the Christoffel symbols at first order in ε. These equations are readily integrated to yield As for the scalar amplitude, δE µ (x 0 ) corresponds to boundary values (4.29) evaluated at the emission point x 0 (x), as defined in Eq. (4.17).
Explicitly, using the standard formula for the Christoffel symbols, Γ (k, u 2 , u 1 ) is found to be given by Summary All in all, the leading order geometrical optics approximation for the electromagnetic field satisfying the emission conditions (4.12)-(4.14) was found to be with

34)
where m and n are unit vectors specifying the propagation direction of the electromagnetic and gravitational waves, respectively, and where the subscript of h indicates the value of the parameter u where the metric perturbation is evaluated.

Reflection at Mirrors
Consider the reflection of the wave (4.33) at a mirror described by the timelike hypersurfacẽ i.e. a surface which -in our coordinate system -is parallel to the emission surface Σ and separated from it by a coordinate distance .
As the simplest model, we assume the mirror to be perfectly reflecting. (More realistic models would use generalisations of the Fresnel equations to the general relativistic setting, as described in Ref. [38].) As explained in Appendix D, the corresponding boundary conditions relating the where R β α is the reflection along the surface normal to ν as defined in Eq. (4.2): where the index of ν is raised with the full metric g. These equations say that the phase and amplitude of the reflected wave coincide with those of the incident wave, and that the reflected polarisation is obtained from the incident polarisation by applying the reflection operator R and changing the sign. (One could, alternatively, include the negative sign inǍ , or add ±π to the reflected phase, but this does not alter the overall fieldF .) To describe the reflected eikonal, we setǩ to be the unperturbed reflected wave vectoř and analogously to the definition of x 0 in Eq. (4.17), let With this notation, the eikonal of the reflected wave is readily found to bě where the function H is defined in Eq. (4.18). For the scalar amplitudeǍ = 1 + εǍ (1) , one obtains the transport equatioň A direct calculation shows that so that the perturbation of the reflected amplitude is found to bě To compute the reflected polarisation, we writě f = (dψ/ω) ∧Ě , (5.10) in complete analogy to Eq. (3.7). Since, according to Eq. (D.2), the reflected eikonal satisfies the last condition of Eq. (5.2) can be solved by imposinǧ The solution to the transport equatioň 13) with this boundary value is then found to be given by where, in the last two lines, we have used the decomposition of the emitted polarisation into normal and tangential parts, given in Eq. (4.10).
Finally, we evaluate the field back at the emission surface Σ, for which we introduce some notation. Given any point x on Σ, denote by x R the point obtained by following the unperturbed incoming light ray back in time (along −ǩ) until one reaches the mirror surface. Following the unperturbed light ray further back in time (along −k), one reaches the point x E on Σ, where the ray was initially emitted a time 2 ago. Explicitly, the coordinates of these points are given by With this notation, the returning field on Σ can be written aš where the subscripts indicate the values of the parameter u in the metric perturbation h, and the last expression is understood in matrix notation (i.e. 1 is the identity matrix, Γ is the matrix defined in Eq. (4.31), and R (0) Γ denotes the matrix product of R (0) and Γ as defined in Eqs. (4.31) and (5.3)).

Non-Polarising Beam Splitters
Considering now a non-polarising 50:50 beam splitter, we model the two "output ports" by timelike hypersurfaces Σ i and Σ ii (as sketched in Fig. 1), whose unperturbed normals m i and m ii are orthogonal in the background geometry. The associated unit one-forms (in the perturbed geometry) are . Radiation sent into the beam splitter partially passes through, and is partially deflected. This deflection can be described using the same methods as discussed in the previous section, namely by reflecting the incoming field strength tensor along the normal cf. the description of reflection in Appendix D. To this end, define the reflection operator which acts by interchanging ν i and ν ii , and leaving all vectors orthogonal to them unchanged. In With this, we can now formulate a simple model of the beam splitter. Consider incoming radiation at the surface Σ in (coming from a laser) of the form where E in is as in Eq. (4.4) and the normalisation is such that the emitted radiation has unit energy density.
We then assume the fields at the output ports (facing the mirrors) to be The outgoing fields have only half the amplitude of the incoming field, in accordance with energy conservation. This means that half the radiation passes through the beam splitter unchanged, while the other half is reflected along the normalν defined above.
To obtain an expression for the amplitude everywhere on Σ, one must then construct a solution to the constraint equations with this prescribed behaviour at the spatial origin. The method of doing so is described in Appendix C, but the precise details are irrelevant for the applications where.
These fields propagate towards the mirrors, are reflected there, and return to the beam splitter.
In the course of this propagation, the fields undergo changes as described by the equations (5.17)-(5.20). Irrespective of the precise form of this transformation, the returning radiation (restricted to Σ i and Σ ii ) can be written in the form Finally, these fields pass through the beam splitter to reach the output port. To model this, we assume the radiation impinging at Σ i to be deflected as above, while the radiation impinging at Σ ii is inverted because of the π phase shift at a beam splitter. This translates to the following formula for the output field at the spatial origin:

Interferometer Output Signal
Having developed a description of light propagation, reflection at mirrors and at beam splitters, all in the presence of a plane gravitational wave, we can now compute the output of a simple Michelson interferometer as sketched in Fig. 1.
Ray I One ray is first transmitted by the beam splitter, then reflects at the mirror surface {m i .x = i }, and propagates back to the beam splitter where it is finally deflected towards the detector.
The radiation sent into the interferometer arm is precisely of the form (4.12)-(4.14), rescaled only by a factor 1/ √ 2 in accordance with Eq. (6.6). Up to this global factor, the field returning to the beam splitter is thus obtained from Eq. (5.16) by the substitution where x E i and x R i are defined as in Eq. (5.15), except for this same substitution: According to Eq. (6.9), the corresponding radiation reaching the output port of the beam splitter is obtained by applying R BS to both the local wave vector and the polarisation, and changing the overall sign of the latter. A direct calculation shows that so that the contribution of this first ray to the output field (omitting the háček) is where

5)
where However, due to the deflection at the beam splitter before the emission in the detector arm, we must also substitute as follows from Eq. (7.3). This leads to the following expression for the contribution of the second ray to the interferometer output: where 14) Here, we have decided to put the π phase shift at the beam splitter into the eikonal so that all amplitudes are positive.
These results can now be used to compute the final output field.
Interferometer Output So far, we have worked with a complex electromagnetic field for simplicity. The physical field, however, is obtained by taking the real part. Assuming the emitted polarisation E to be real (corresponding to linear polarisation), the real fields are obtained simply by replacing e iψ by cos(ψ) in Eqs. (7.4) and (7.13). The overall real field at the output port is There are multiple readout schemes to measure the gw perturbations of this resulting field, notably the homodyne, heterodyne, and dc readout schemes [39][40][41]. While the homodyne readout scheme was used in the first generation of gw detectors, the Advanced ligo, Advanced Virgo, kagra and geo 600 detectors now use the dc readout scheme (dcr) [42][43][44][45]. However, as the planned ligo a+ upgrade will transition to the balanced homodyne readout scheme (bhr) [46], we will describe both the dcr and bhr schemes in the following analysis.
In all cases, the signal is read out using photo diodes, which are sensitive to the energy density of the em field. As shown in Appendix E, the time-averaged intensity of a field of the form (7.18) at first order in ε is (7.19) where, notably, no inner products of the wave vectors and/or polarisation vectors enter. This is because linear perturbations of these normalised vectors produce only quadratic corrections to the relevant inner products and are thus negligible. This means that although the gravitational wave causes perturbations of the eikonal, amplitude, local wave vector, and the polarisation, only the first three perturbations produce a signal in the detectors considered.
Note that we have chosen the input field to be normalised to unit energy density. In the general case, the expression given here is the output power normalised to the input power. Equation (7.19), as well as its generalisation to more than two interfering rays, can be used to compute the signals in the dcr and bhr schemes. Here and in the following, all expressions are evaluated at the spatial coordinate origin, where the detector is positioned.

DC Readout Scheme
In the dcr scheme, one directly measures the intensity of the interferometer output, yielding the signal In the unperturbed case, this reduces to which vanishes if both arm lengths are chosen to be equal. The gravitational wave produces three correction terms, corresponding to perturbations of the eikonal, amplitude, and wave vector: Using Eqs. (7.5)-(7.7), Eqs. (7.14)-(7.16), as well as Eq. (4.11), these first-order perturbations evaluate to all of which vanish for equal arm lengths. Hence, the dcr scheme requires a slight asymmetry ∆ = i − ii to obtain a linear response to the gw perturbation.
The three terms appearing here can be interpreted as follows. S dcr,ψ arises due to the perturbation of the optical path lengths along the two interferometer arms. The multiplicative factor sin(2ω( i − ii )) entails that the phase response vanishes when the unperturbed output power is extremal, and is maximal in between. The second term, S dcr,A , arises from amplitude perturbations, which can be regarded as scintillations of the light signal, which grow with the length of the interferometer arms (for very long distances, corrections from higher order terms of the geometrical optics expansion might become relevant). The strength of this effect is directly proportional to the output intensity in the unperturbed configuration, as can be seen from the multiplicative factor sin 2 (ω( i − ii )). Finally, the last term, S dcr,K , arises because the gravitational wave modulates the frequency of the light rays, which translates to a modulation of the output power, similar to the amplitude modulation.
Hence, the first contribution, S dcr,ψ , describes the standard phase response (as it is determined from the eikonal alone), while the remaining terms describe the corrections beyond the eikonal equation.

Balanced Homodyne Readout Scheme
In the bhr scheme, one mixes the output field (7.18) with a further ray -typically derived from the same laser source -which is not sent through the interferometer. For this ray the gw perturbations are negligible so that one may set where ϕ denotes the readout angle (an experimentally controlled parameter), and the amplitude normalisation is the same as in Eq. (6.5). Sending this ray and the output field (7.18) orthogonally through a further 50:50 beam splitter, the resulting intensities at the two output ports evaluate to where the difference in sign of the interference terms involving the third ray arises because of the π phase shift at the beam splitter. The signal in the bhr scheme is then obtained by measuring the difference of the two intensities, i.e.
The unperturbed signal in this readout scheme is which vanishes for all values of ϕ if the arm lengths are chosen to be equal.
As for the dcr scheme, it is convenient to separate the first-order perturbations of the phase, amplitude, and frequency by writing bhr + εS bhr,ψ + εS bhr,A + εS bhr,K + O(ε 2 ) , (7.31) which, for equal arm lengths, evaluate to While the phase response in the bhr scheme (7.32) is the same as in the dcr scheme (7.23), the amplitude and frequency responses differ. Specifically, Eqs. (7.24) and (7.25) contain the sum of the perturbations in the two interferometer arms, while Eqs. (7.33) and (7.34) contain their difference.

Detector Pattern Functions in the Low Frequency Limit
In the low frequency approximation, which applies if the gw wavelength is much longer than the interferometer arms, the above formulae can be simplified significantly by assuming the gw waveform to vary slowly along the light ray trajectories. Formally, this corresponds to a Taylor expansion of h near x = 0 to leading order (i.e. including the first non-trivial term).
To set up the notation, we first consider the well-known phase perturbations in Eqs. (7.23) and (7.32). To analyse the phase response of the interferometer in both readout schemes, set Assuming (almost) equal arm lengths, the low-frequency limit of S ψ evaluates to where the metric perturbation, h, is evaluated at the coordinate origin. Decomposing h in terms of the two polarisations "+" and "×" (described in more detail below) as one can write where The functions F + ψ and F × ψ are known as the detector pattern functions). To describe their behaviour for various incidence angles of the gravitational wave, one must parameterise the gw wave vector as well as the gw polarisation tensors A + and A × . Specifically, one must relate the For comparison, we note that the transformation matrix given in Ref. [2] uses the "ZXZ" convention for Euler angles and is thus obtained from the one here by the substitution Φ → Φ − π/2, Ψ → Ψ + π/2. To parameterise the gw polarisation tensors, denote by e 1 , e 2 , e 3 the canonical orthonormal basis of R 3 . The X and Y basis vectors of the gw frame then correspond to the following vectors in the detection frame: These quantities depend implicitly on all Euler angles Φ, Θ, Ψ, but it suffices to consider Ψ = 0 since the general case is obtained from the particular using the formula With this parametrisation of the gw polarisation, and choosing m i = e 1 as well as m ii = e 2 , the pattern functions (8.6) evaluate to (8.11) in agreement with standard results [47,Eq. (58)]. The absolute values of these functions are plotted in Fig. 3, where it is seen that the gw perturbation of the phase is maximal if the gw propagates orthogonally to both interferometer arms.
Let us now consider the terms which are beyond the eikonal equation.

DC Readout Scheme
The first corrections beyond the eikonal equation in the dcr scheme are given by Eqs. (7.24) and (7.25). Expanding these expressions to first order in the low-frequency limit, one obtains where, as above, the metric perturbation is evaluated at the spatial coordinate origin. Since the expressions coincide up to an overall factor, the corrections beyond the eikonal equation can be summarised as Separating the two gw polarisations as in Eq. (8.4), this can be written as where the antenna pattern functions describing the response beyond the eikonal equation in the dcr scheme are Setting m i = e 1 , m ii = e 2 and inserting the parametrisation of A + and A × from Eq. (8.9), they evaluate (for Ψ = 0) to The angular dependence of F + dcr (in absolute value) is plotted in Fig. 4a. Clearly, the effect is maximal when the gw wave vector lies in the xy plane (spanned by the interferometer arms), and a gw propagating orthogonally to both arms produces no amplitude perturbations (contrary to the phase perturbation, which is maximal for orthogonal propagation).

Balanced Homodyne Readout Scheme
The analysis of the signals in the bhr scheme proceeds analogously to those in the dcr scheme. Expanding Eqs. (7.33) and (7.34) to leading order in the low-frequency limit, one obtains Also here, the first corrections beyond the eikonal equation can be summarised succinctly as (8.20) and decomposing the gw waveform into its two polarisation components (as above), this can be written as where the antenna pattern functions for the response beyond the eikonal equation in the bhr scheme (for Ψ = 0) are Up to a factor −1/2, these functions coincide with the eikonal pattern functions given in Eq. (8.11).
Their absolute values are plotted in Figs. 4b and 4c, for direct comparison with the pattern function of the dcr scheme.

Angular Efficiency Factors
To compare the sensitivity of the two readout schemes, it is customary to define the angular "sky-average" of a function f as (8.23) and to define the angular efficiency factor as For the eikonal response, and the dcr and bhr responses to corrections beyond the eikonal equation, these factors evaluate to This shows that the dcr scheme is more susceptible to amplitude and frequency perturbations than the bhr scheme by a factor of F 2 dcr /F 2 bhr = 4/3.

Waveforms
To illustrate the signals produced by the first corrections beyond the eikonal equation, we plot the detector responses for the gw waveform of peak amplitude ε ≈ 4.75 × 10 −21 , shown in Fig. 5. This waveform was generated from the SEOBNRv4 model [49] using the PyCBC software package [50]. The inspiral phase of the waveform shown here starts with a frequency of Defining the dimensionless functions (8.26) the signals at leading order in the low-frequency limit can be written as

27)
δS dcr ≈ ω g f dcr sin 2 (ω∆ ) , δS bhr ≈ ω g f bhr cos(ϕ − 2ω ) .   both interferometer arms. In this case, S ψ and δS bhr are maximal, while δS dcr vanishes. For comparison, Figure 6b shows the signals for the case where the gw propagates parallelly to one of the interferometer arms (in which case the phase response is sensitive only to the + polarisation mode). In this configuration, both the dcr and the bhr schemes are susceptible to beyond-eikonal corrections and, in fact, the dcr scheme produces a stronger response.
However, since the corrections shown here are suppressed relatively to the phase response by the frequency ratio ω g /ω ≈ 10 −13 , the experimental detection of such effects in the foreseeable future seems highly unlikely.

Discussion
We have determined the response of laser interferometric gravitational wave detectors at one order of accuracy beyond the eikonal equation by solving also the transport equations for the scalar amplitude and the polarisation vector. To this end, we have modelled the emission of monochromatic plane waves, described their propagation and reflection, as well as their behaviour at beam splitters. The final result for the interferometric output intensity encompasses gwinduced perturbations of the optical path length (via the eikonal ψ), modulations of the optical frequency (as given by the time-component of the local wave vectors K), and scintillations (described by the amplitudes A ), but was found to be independent of the perturbation of the polarisation vector E arising from a non-trivial holonomy integral along the light rays.
The analysis presented here goes beyond previous results which have either allowed for arbitrary gw incidence angles but were restricted to the level of the eikonal equation [26], or have described amplitude and polarisation effects but were restricted to special alignments [25; 27; 31]. Moreover, the analysis here does not rely on unidimensional approximations such as in Refs. [24; 29], and the result is well-behaved also for small angles between the gw and em wave vectors, where calculations based on the vector potential often find difficulties [22].
Finally, we comment on the model assumptions made here. The geometry of the interferometer was assumed to be well-described in terms of the background geometry (when using transversetraceless coordinates), which relies on the assumption that all material points of the interferometer follow geodesic motion. This will not hold exactly true in experiment, and deviations from this model would require detailed descriptions of the elastic properties of the respective materials.
However, since the mirrors and beam splitters are much smaller than the interferometer arm lengths, we expect such corrections to be suppressed by the ratio of the object size compared to the arm length and thus negligible for most applications. Further, we have made specific assumptions of the light sent into the interferometer. While the assumptions on the eikonal and amplitude are simple to interpret (emission occurs with definite and fixed frequency and the emitted radiation has constant energy density), the emitted polarisation is more difficult to model, leading to the unspecified term E in Eq. (4.4). However, the final result is independent of E due to it being orthogonal to the unperturbed polarisation E (0) . Hence, we expect most other models to differ mainly in the choice of the emission surfaces on which the boundary conditions are prescribed, rather than in the data prescribed thereon. However, provided the surfaces differ from the ones used here only at first order in ε, such arising corrections can be described entirely using the unperturbed transport equations since corrections to the latter would produce terms quadratic in ε.
We expect that the calculations presented here can be generalised without much complication to Fabry-Pérot cavities, while optical cavities whose mirrors are not plane-parallel (but rather concentric or confocal, for example) could be described using similar techniques, but would require different boundary conditions adapted to the problem.

A Derivation of the Polarisation Transport Equations
The equations of geometrical optics can be motivated by writing the electromagnetic field strength tensor as and inserting for f µν a formal expansion in inverse powers of (iω), see e.g. Refs. [36,Sect. 4 where a dot indicates contraction of adjacent indices. Together, these equations imply Assuming ψ to satisfy the eikonal equation one is left with the intermediate result where HdR f = −d(∇ · f) − ∇ · (df) is the Hodge-de Rham wave operator, which is related to the connection-induced wave operator by the Weitzenböck identity where R α βµν denotes the Riemann curvature tensor. A direct calculation shows which can be thought of as constraint equations limiting the number of degrees of freedom to two. The geometrical optic scheme described here is in accord with the one given in Ref. [52].
We emphasise that the series (3.5) is first and foremost a formal series, which is not necessarily convergent, but rather asymptotic in many typical applications, see e.g. Ref. [37,Sect. 2.7] for a discussion of the asymptotic nature of such series for initial data prescribed on spacelike hypersurfaces. In the concrete setup here, our analysis of the full Einstein-Maxwell equations (in preparation) will assess the quality of the leading order geometrical optics approximation and estimate the errors made in using this approximation.

B Interpretation of the Polarisation Vector
In this section, we show how the one-form E in the field strength tensor is related to the polarisation of the electromagnetic wave. Here, E is actually only defined up to addition of a multiple of dψ.
Consider an observer O with four-velocity u (i.e. a unit timelike vector). Since the eikonal ψ describes the rapidly oscillating phase of the electromagnetic wave, the measured frequency is Regarding the electric field E and the magnetic field B as one-and two-forms, respectively, the decomposition of the field strength can be written as where E and B are such that their contractions with u vanish (i.e. they are spatial for O). In particular, they are determined by We now apply this decomposition to the field given in Eq. (B.1). Using the freedom to add any multiple of dψ to E, we write which is chosen such that u µĒ µ = 0. The above formula for E then yields so thatĒ determines the direction of the electric field and A ω O /ω its amplitude (note thatĒ has the same norm as E). To obtain an expression for the magnetic field B, we first decompose the eikonal gradient as where n µ is a unit spacelike vector orthogonal to u. Using the second equation in (B.4), one which is the standard formula stating that E, B and n form a right-handed orthogonal system, and that the electric and magnetic fields have the same norm.

C Constraint Equations for Boundary Data
As is well-known from the 3 + 1 decomposition of Maxwell's equations, one is not free to prescribe the electromagnetic field on spacelike hypersurfaces arbitrarily. Instead, the initially prescribed electric and magnetic fields E and B (in Heaviside-Lorentz units) must satisfy D · E = ρ and D · B = 0, where ρ is the charge density on the initial hypersurface, and D is the spatial covariant derivative.
Similarly, we show here that boundary values prescribed on a timelike hypersurface Σ must also satisfy certain constraint equations, which we compute explicitly for the case of normal emission, as considered in the main body of the text.
To this end, consider a hypersurface-orthogonal unit timelike vector field u. With respect to this field, one may decompose the field strength, F , as well as its Hodge dual, * F , in terms electric and magnetic fields as where E and B are both spatial in the sense that their contraction with u vanishes, cf. Eq. (B.3).
Maxwell's equations in vacuum then take the form where L u denotes the Lie derivative along the vector field u.
Consider now a timelike hypersurface Σ with unit normal n, such that u and n are orthogonal.
This entails that n is spatial with respect to u, and also that u is tangent to Σ. Contracting Eqs. (C.2) and (C.3) with the surface normal n and imposing normal emission in the sense where a = ∇ u u is the acceleration vector. These equations are constraints on boundary data, as they only involve tangential derivatives, but no derivatives normal to Σ.
It turns out that the first terms involving L u of E and B give no contribution. Indeed, writing the Lie derivative in terms of the Lie bracket L u E = [u, E], and noting that u and E being tangent to Σ implies that their Lie bracket is also tangent to Σ, one finds L u E µ n µ = L u B µ n µ = 0.
Hence, the constraint equations reduce to which are the same for both E and B. Assuming further (as is the case for plane waves) that n, E, B form a right-handed orthogonal frame with |E| = |B|, such that Eq. (B.9) applies, one may eliminate one of the fields from the constraints. In this case, one has so that the first equation of (C.7) takes the form whereδ is the divergence operator of the intersection of Σ with a time-slice, and similarly denoting the exterior derivative intrinsic to such a surface byd, one obtains the constraint equations in the form where a = a − n(n · a) is the tangential part of the acceleration, and where we have identified the vector a with its metric-equivalent one-form. In particular, for geodesic u, this simplifies tõ Here, ∆ denotes the spatial scalar Laplacian on Σ. Note that χ is not required to remain bounded: instead, since the physically relevant field is E = dχ, we require χ to have bounded first derivatives.
To construct perturbative solutions, it is convenient to write the Laplace equation using the "densitised contravariant metric" in the form whereg AB = √ detgg AB , withg being the metric induced by g on Σ (at any instant of time).
Inserting the perturbative expansion one obtains where ∆ (0) = δ AB ∂ A ∂ B is the unperturbed Laplacian on Σ. Sinceg is a function of κ.x alone, one finds a particular solution to be given by (1)BC (u) du , (C. 17) and requiring the overall function χ to have bounded first derivatives, one is only free to add homogeneous solutions of the form χ A x A , where the χ A 's are constant. Taking the derivative of and choosing the constants χ A such that E (1) assumes a prescribed value at the spatial origin, one obtains

D Boundary Conditions for Perfect Reflection
In this section, we give a generally covariant description of perfect reflection at the level of geometrical optics. For this, we consider a field of the form where F 1 describes a wave impinging on a reflecting surface Σ, and F 2 is the reflected wave. We first consider how dψ 2 relates to dψ 1 (on Σ), and then we determine E 2 from E 1 .
For the eikonal functions, we assume the standard matching condition that ψ 1 and ψ 2 coincide on Σ, cf. e.g. Ref. [36,Sect. 2.5]. This entails that dψ 2 = dψ 1 + αν (on Σ), where ν is the unit conormal to Σ, and α is some function on Σ which can be determined as follows. Since we are concerned with waves propagating in different directions, we must have α = 0. Further, taking the norm on both sides and using the fact that both functions ψ 1,2 satisfy the eikonal equation, one obtains α = −2ν σ ∇ σ ψ 1 . Using the definition (5.3), this can be written as i.e. the wave vector is simply reflected along the normal to Σ.
We now turn to the relationship between E 1 and E 2 . For this, consider an arbitrary point on Σ (which we shall keep fixed from now on) and let u denote any (future-pointing) unit timelike vector at this point, which is orthogonal to ν (the final result will be independent of this choice).
By the argument of Appendix B, we may assume without loss of generality that the covectors E 1,2 are orthogonal to u. Using local inertial coordinates adapted to u (at the considered point), we may write where ω O is the frequency measured by observers with four-velocity u, cf.
Now, since we are working in a local inertial frame, we may use the standard equations ν × E = 0 , ν · B = 0 , (D. 6) which are equivalent to the following equations, in exterior algebra notation: The first equation here entails that E 2 = −E 1 + βν for some constant β, and the condition of perfect reflection requires that E 1 and E 2 have the same norm, yielding The second matching condition ν ∧ B = 0 is then automatically satisfied. Indeed, since ν ∧ m 2 = ν ∧ m 1 , one has ν ∧ B = ν ∧ m 1 ∧ E , (D.9) which vanishes because ν ∧ E = 0.
So far, we have shown that in local inertial coordinates (using a decomposition of F adapted to the four-velocity u) the conditions for perfect reflection are which, together, imply But since this last equation (which is independent of the decomposition of F into E and B) is tensorial, it holds in all coordinate systems, and since the considered point on Σ was arbitrary, it holds everywhere on this surface.

E Output Power of an Interferometer
In Eq. (7.18), the field at the output port of the interferometer was found to be of the form where k i and k ii are null vectors, and E i and E ii are spacelike unit vectors (with respect to the full space-time metric g), which satisfy g(k i , E i ) = 0 and g(k ii , E ii ) = 0. In this section, we compute the time-averaged energy density for such a field at first order in the gw amplitude ε.
The intensity measured by a photodetector is given by the temporal component of the energy-momentum tensor Considering the first term only for the moment, we have which can be simplified considerably. Since E i and E ii have unit norm, the first line simplifies to (A i k i0 ) 2 cos(ψ i ) 2 + (A ii k ii0 ) 2 cos(ψ ii ) 2 . For the remaining expression, we make use of the fact that in the unperturbed problem the k's coincide, and so do the E's.
Moreover, the inner products g(E i , k ii ) and g(E ii , k i ) in the third line are of order ε since the unperturbed wave vectors are orthogonal to the unperturbed polarisation vectors. But since the unperturbed polarisation vectors are purely spatial, the multiplying factors E i0 and E ii0 are of order ε as well, so that both summands in the last line are of second order and thus negligible.
Hence, to first order in ε one has Now consider the second term in the energy-momentum tensor: By a similar argument as before, one finds this expression to be of order ε 2 . Hence, at first order in ε, the energy density reduces to Finally, we compute the time average of this expression since the read-out cannot resolve the fast oscillations of the electromagnetic field itself but only the slow modulation arising from the gravitational wave. (Recall that optical frequencies are in the order of 500 THz while observed gw frequencies are roughly in the range of 10 to 500 Hz). Formally, such an average is given by In particular, the amplitudes A , the polarisation vectors E and the wave vectors k oscillate slowly and may be taken outside the average, leaving us only with the eikonal terms, which oscillate rapidly in time as −ωt. Consequently, we have cos(ψ i ) 2 = cos(ψ ii ) 2 = 1 2 , cos(ψ i ) cos(ψ ii ) = 1 2 cos(ψ i − ψ ii ) , (E.10) so that the averaged energy density at the output is To conclude this section, we verify explicitly that the conditions (E.5) are satisfied for the fields given in the main body of the text. Writing Eqs. (7.7) and (7.16) as BS κ + m ii (κ.ǩ i ))] , (E.12) K ii = −dt + ν ii + 1 2 β[ωκ + m ii (κ.ǩ ii ))] , (E. 13) where the precise form of the factors α and β is irrelevant for the argument here, one finds g(K i , K ii ) = 1 2 α ω(−dt − m i ).κ + κ.ǩ i + 1 2 β ω(−dt − m ii ).κ + κ.ǩ ii + O(ε 2 ) , (E.14) which vanishes sinceǩ i = ω(−dt − m i ) andǩ ii = ω(−dt − m ii ).
Similarly, using the expressions for the polarisation vectors from Eqs. (7.8) and (7.17) together with the explicit form of δE given in Eq. (4.10), as well the fact that R ii and R BS are orthogonal matrices, their inner product evaluates to BS E (0) ) κ.x + 1 2 εh(R