Second-order gravitational self-force in a highly regular gauge: Covariant and coordinate punctures

Gravitational self-force theory is the primary way of modelling extreme-mass-ratio inspirals (EMRIs). One difficulty that appears in second-order self-force calculations is the strong divergence at the worldline of the small object, which causes both numerical and analytical issues. Previous work [Phys. Rev. D 95, 104056 (2017); ibid. 103, 124016 (2021)] demonstrated that this could be alleviated within a class of highly regular gauges and presented the metric perturbations in these gauges in a local coordinate form. We build on this previous work by deriving expressions for the highly regular gauge metric perturbations in both fully covariant form and as a generic coordinate expansion. With the metric perturbations in covariant or generic coordinate form, they can easily be expressed in any convenient coordinate system. These results can then be used as input into a puncture scheme in order to solve the field equations describing an EMRI.


I. INTRODUCTION
Extreme-mass-ratio inspirals (EMRIs) [1] will be a key source of the gravitational waves that will be detected by the Laser Interferometer Space Antenna (LISA), a future space-based gravitational wave detector [2,3].An EMRI features an object of mass m ∼ 1-10 2 M ⊙ slowly spiralling into an object of mass M ∼ 10 5 -10 7 M ⊙ .The smaller object is a compact object, such as a black hole or neutron star, whereas the larger object is a supermassive black hole, existing in the centre of most galaxies [4][5][6].
As the mass ratio, ϵ := m/M ∼ 10 −5 , is very small, the inspiral occurs over a long timescale, with the smaller object expected to complete ϵ −1 ∼ 10 5 intricate orbits before plunging into the central black hole [7,8].Due to the large number of orbits occurring near to the supermassive black hole, the gravitational waves emitted are expected to provide an excellent picture of the geometry of the black hole in the strong-gravity regime.This will allow highly accurate tests of general relativity to be performed [8][9][10][11].

A. Gravitational self-force
The primary method of modelling EMRIs is through a perturbative method known as gravitational self-force theory [12][13][14][15][16].The self-force refers to the process by which changes in an external field caused by an object's dynamics propagate back and affect the motion of the very same object.This method expands the metric describing the geometry of the full spacetime, g µν , around a known, background metric, g µν , with perturbations, h µν , caused by the presence of the small object.In an EMRI, the disparate sizes of the small and large object lead to a natural perturbative parameter, the mass ratio between the two objects, ϵ ≪ 1.One can then write the full spacetime metric as the sum of the background spacetime and these perturbations, where In the case of an EMRI, the background metric describes the geometry of the large black hole if it were isolated in spacetime and is taken to be either the Schwarzschild [17] or Kerr [18] metric.At the leading order in the mass ratio, the small object's worldline, γ, is a geodesic of the background spacetime, g µν .The metric perturbations then alter the motion at higher orders and exert a self-force on the body, moving it away from a background geodesic.This can be written as which reduces to the geodesic equation when ϵ → 0. In Eq. ( 3), z α are coordinates on the accelerated worldline, γ, τ is the proper time in the background metric, g µν , D/ dτ := u µ ∇ µ is the covariant derivative along the worldline and is compatible with g µν , u α := dz α /dτ is the four-velocity and f α n is the nth-order self-force.The self-force (or at least part of it) causes the orbit to evolve at a rate of Ė/E ∼ ϵ, resulting in an inspiral over the radiation reaction time, t rr ∼ E/ Ė ∼ 1/ϵ [13].Here, E is the orbital energy and is one of three constants of motion that completely describe the geodesic of a test particle in the background Kerr spacetime; the other two are the azimuthal angular momentum, L z , and the Carter constant, Q [19].
One challenge is that we are required to go to second order in the mass ratio in order to model the waveforms arXiv:2309.03778v2[gr-qc] 14 Feb 2024 accurately.This is a result of the requirement that for us to extract information from the data gathered by LISA, the phase of the waveform must be accurate to within a fraction of 1 radian.A precise argument for the need for second order was made by Hinderer and Flanagan [20].The orbital parameters, J B = {E, L z , Q}, slowly evolve over the radiation reaction time, t rr ∼ 1/ϵ.This motives the introduction of a "slow time", t = ϵt, so that J B = J B ( t).The orbital frequencies, Ω A = {Ω r , Ω θ , Ω ϕ } in the case of Kerr, are functions of the orbital parameters, J B ( t), and have perturbative expansions, Ω A (J B , ϵ) = Ω (0) where Ω (n≥1) A are the nth order corrections to Ω (0) A due to the conservative part of the self-force.The orbital frequencies evolve with respect to the time, t, as (1) (2) where A is constructed from the nth-order dissipative force.These can then be related to the orbital phases by so that where the adiabatic term, φ (0) A , is constructed from Ω (0) A and F (1) A , and the first post-adiabatic (1PA) term, φ A , is constructed from Ω (1)

A and F (2)
A .One can see this through noting that an integration over t introduces a factor of 1/ϵ through dt = dt d t d t = ϵ −1 d t.Therefore, to calculate the orbital phases with an error much less than order-ϵ 0 requires the entirety of the first-order self-force and the dissipative part of the second-order self-force.
It should be stressed that the conservative piece of the first-order self-force and the dissipative piece of the second-order self-force are on equal footing: even if one has the entirety of the first-order self-force (both dissipative and conservative parts), if one does not have the dissipative piece of the second-order self-force then one cannot correctly track the motion of the small object.
As to the current status of the self-force field, at first order, full inspirals driven by the self-force can be computed for generic orbits in the Schwarzschild spacetime for a spinning small object [21][22][23][24].One can calculate the full first-order self-force for a non-spinning small object on any generic bound orbit in Kerr [25].Adiabatic inspirals in Kerr have been performed for equatorial [26] and generic [27] orbits with Ref. [28] performing an equatorial inspiral using the entirety of the first-order self-force.
B. Local form of the metric perturbations, puncture scheme and infinite mode coupling

Metric perturbations and effective stress-energy tensor
To find the local form of the metric perturbations, one uses the method of matched asymptotic expansions (for a general introduction to matched asymptotic expansions, see, e.g.Refs.[46,47], and for an introduction to their use in self-force, see, e.g.Ref. [13]).When close to the small object, the expansion from Eqs. ( 1)-( 2) breaks down as the gravitational field from the small object dominates over that of the background spacetime.One then introduces a second expansion that focuses in on the small object and then matches this with the external expansion at some appropriate lengthscale.This is then combined with the vacuum Einstein field equations to solve for the metric perturbations, h µν .
The metric perturbation can be split into two fields [48], where h R µν and h S µν are the regular field and singular field, respectively.The regular and singular fields can be expanded in an analogous manner to Eq. ( 2), as The regular field has the form of a Taylor series centred on the worldline of the small object and satisfies the vacuum Einstein field equations, throughout the entire spacetime.When combined with the background metric, it forms a smooth, vacuum effective metric that determines the local geometry that the small object "feels", Through second order, the trajectory of the small object (assuming zero spin) is governed by the equation of motion [49,50] which can be written as a geodesic in the effective spacetime, gµν , as where all quantities with tildes are defined with respect to gµν .This correspondence is known as the generalised equivalence principle [50], which states that (ignoring finite-size effects) a compact object immersed in an external gravitational field will follow a geodesic in some effective metric whose geometry is determined by its own physical mass.
The remaining part of the metric perturbations, the singular field, contains information about the small object's multipole structure [48].Schematically, it has the form where r is the proper spatial distance to γ and M α /S α are the mass/spin dipole terms, respectively.As in previous work, we enforce that the mass dipole and any higherorder corrections to it vanish.This ensures that γ tracks the small object's centre of mass [50][51][52].
In certain classes of gauges, the small object also has the effective stress-energy of a point mass in the effective spacetime [53,54]. 1 Using this effective stress-energy tensor, the field equations can be written in the form where T µν is the Detweiler stress-energy tensor, and all quantities with tildes are defined with respect to the effective metric.The existence of this stress-energy tensor was first postulated by Detweiler [53] and explicitly derived in Ref. [54] (hereafter Paper I).One can also write the left-hand side of Eq. ( 18) in terms of effective quantities as [53,54] δ demonstrating that the system can be described as a linear perturbation of an effective background.
It should be noted that the split into regular and singular fields is not unique [55], but we choose the split to match that of, e.g.Refs.[48][49][50]54], ensuring that the regular and singular fields satisfy the properties listed above.That is, the regular field is smooth on the worldline of the small object, forms the effective metric, gµν , and satisfies the generalised equivalence principle.In addition to the non-uniqueness of the split, it should be emphasised that neither h R µν nor h S µν represent the true physical field; only their sum h µν = h R µν + h S µν does.We stress that the results discussed in this section are all derived from the principle of matched asymptotic expansions.One does not start by assuming that the small object is described by a point-particle stress-energy with some effective equation of motion.Instead, one uses the matching process at each order in ϵ to rigorously derive these properties from first principles.

Puncture scheme
To date, all second-order calculations have involved the use of a puncture scheme [29,30,32]; see, e.g.Refs.[14,16,55,56] for technical details. 2 In this scheme, one introduces a puncture field, h P µν ≈ h S µν , that approximates the singular field to some sufficient order in r away from the worldline, and goes to zero beyond that.From this, one can define a residual field, so that h R µν ≈ h R µν near γ.These fields are then analytically extended down to the worldline, and one solves for the residual field, h R µν , with the puncture field as the source, instead of directly for the physical field, h µν .
We wish to be able to replace h R µν with h R µν in the equation of motion (14).This is possible if h R µν and its first derivatives are identical to h R µν .To ensure this, we impose the conditions lim lim where z µ is a point on the worldline.Explicitly, to calculate the second-order self-force, we need to go to order r in our second-order punctures so that our residual field is once differentiable.Substituting Eq. ( 21) into the field equations and expanding the residual and puncture fields order-by-order, as in Eq. (2), we get These equations can be promoted to the entire domain, including r = 0, provided that the puncture field is known to a sufficiently high order in r; see the discussion after Eq. ( 13) of Paper I. Combining the field equations with the equation of motion (14), one can solve the coupled system of equations and determine how the small object travels in spacetime.

The problem of infinite mode coupling
When implementing the puncture scheme at second order, one encounters the problem of infinite mode coupling [61].To take advantage of the symmetries of the spacetime, one decomposes the metric perturbations into a suitable basis of harmonics. 3For example, in Schwarzschild, one could choose Barack-Lousto-Sago tensor spherical harmonics [62,63], so that the metric perturbations can be decomposed as With the modes written as such, to calculate a single mode of δ 2 G µν [h 1 , h 1 ] requires one to calculate the infinite sum of products of first-order modes [61,64], ] is a certain differential operator [64].From Eq. ( 16), we see that h S1 µν ∼ m/r.This means that, generically, the second-order Einstein tensor diverges as ∼ m 2 /r4 at the worldline of the small object as it has the structural form, 4 .After decomposing into modes and integrating over two of the dimensions, one finds that Eq. ( 28) acts as However, the modes of the first-order field are finite on the worldline [65,66], meaning that Eq. ( 28) is attempting to reconstruct a divergent function through summing up finite modes.Thus to get convergence requires one to calculate an arbitrarily large number of modes of the first-order fields to calculate even one second-order mode.
A way to circumvent this problem was provided by Miller et al. [61].Instead of summing over modes, as in Eq. ( 28), one expands the first-order field into regular and singular pieces.After expanding the first-order field, the second-order Einstein tensor in the source of the secondorder field equations has the form One then replaces the regular and singular fields in Eq. (30) with the residual and puncture fields.The terms are sufficiently well-behaved that one may compute the modes directly from the modes of the first-order residual and puncture fields.As described in Ref. [61], the problem is entirely caused by the slow converge of the modes of δ 2 G iℓm [h P1 , h P1 ] as this is the term that causes the non-mode-decomposed second-order Einstein tensor to diverge as ∼ m 2 /r 4 .Instead of summing up the products of the modes of h P1 µν , Miller et al. [61] directly calculate δ 2 G µν [h P1 , h P1 ] in four dimensions using the four dimensional expression for h P1 µν and then decompose this quantity into modes.Unfortunately, while this makes the calculation of the modes of the source possible, it is incredibly computationally expensive and takes up almost all the code runtime when implemented (such as in Ref. [29]).This is due to having to calculate the modes by numerically integrating the complete four-dimensional expression on a grid of r BL and r values.This will not be efficiently extendible when approaching problems involving more complicated dynamics, such as generic orbits in Kerr.

C. Highly regular gauge
The highly regular gauge was introduced by Pound [50] to ameliorate the strong divergences that occur near the worldline of the small object when in a generic gauge.In this gauge, the most singular piece of the second-order perturbation now has the form ∼ m 2 r 0 instead of the ∼ m 2 /r 2 behaviour previously seen; see Refs.[50,54] for a full discussion.One can divide the second-order singular field into two pieces: a "singular times regular" piece, h SR µν ∼ mh R1 µν /r, and a "singular times singular" piece, h SS µν ∼ m 2 r 0 .By simple order counting of m and h R1 µν , we see that, in the second-order Einstein field equations, h SS µν is sourced by δ 2 G µν [h S1 , h S1 ], as they both feature terms ∼ m 2 , and that h SR µν is sourced by δ 2 G µν [h R1 , h S1 ] as both expressions have terms of the form ∼ mh R1 µν .Although the h SR µν term appears more divergent, as discussed in Paper I, its source, δ 2 G µν [h R1 , h S1 ], is well defined as a distribution.The "singular times singular" term causes FIG. 1. Geometric picture of the gauge conditions for the highly regular gauge.The image features a light cone emanating from the worldline, γ.The null vector, k µ , is tangent to the light cone along radially outgoing curves, and the basis vector, e µ A , is tangent to the light cone along spheres of constant luminosity distance, Sr.Based on Fig. 16 from Ref. [67].
the most issues.Acting on the "singular times singular" piece with the linearised Einstein operator, we see that δG µν [h SS ] ∼ m 2 /r 2 .Therefore, we know that the most singular piece of the second-order Einstein tensor can only act as badly 4 as in a generic gauge.This means that when decomposing into modes, the individual modes of the second-order Einstein tensor can behave, at worst, as While this is still divergent, it is much weaker than in the Lorenz gauge.
The highly regular gauge enforces that the local light cone structure around γ is preserved in the perturbed spacetime.To do so, two gauge conditions are imposed on the singular field.Firstly, the metric perturbations vanish when contracted with k µ , the null vector tangent to the future light cone that emanates from the worldline: Secondly, the perturbations are trace-free with respect to Ω AB , the metric on surfaces of constant luminosity distance: where an upper case Latin letter indicates a quantity defined on those surfaces and e µ A := ∂x µ /∂θ A is the basis vector, where x µ are coordinates in the full spacetime and θ A are coordinates on the surface of constant luminosity distance.These gauge conditions ensure that the local background light cone structure is preserved in the perturbed spacetime and that the background luminosity distance is an affine parameter on the null rays that generate the light cones.An image showing the geometric construction is given in Fig. 1.
When working with a puncture scheme, one can impose different gauge conditions on the residual and puncture fields; see the discussions in Sec.IV A of Ref. [68], Sec. VII A of Ref. [50] and Sec.VI A of Paper I. Therefore, to control the singularity structure, one can impose the highly regular gauge conditions on the puncture.Then, one can impose any convenient gauge conditions on the residual field that simplify the left-hand side of the field equations ( 25)- (26).
Reference [50] only provided the leading-order pieces of the second-order metric perturbations in this gauge.Paper I extended this to include all orders needed to perform a numerical calculation of the self-force.These expressions were provided in Fermi-Walker coordinates, a particular coordinate system that is tethered to an accelerated worldline, γ, and is useful for analysing the properties of fields near to this worldline.However, in order to use the expressions in a puncture scheme, one needs to write them in a coordinate scheme specialised to the problem at hand, such as Boyer-Lindquist coordinates (t BL , r BL , θ BL , ϕ BL ) [69].To avoid a potentially complicated coordinate transformation from Fermi-Walker coordinates to the new coordinate system, one can convert the Fermi-Walker expressions into covariant form.This can then be written in the chosen coordinate system.
To do so, one can use the method given by Pound and Miller [55] (hereafter Paper II).This method was developed in order to transform expressions for the singular field in the Lorenz gauge into covariant form.These expressions, after being written in an appropriate coordinate system and decomposed into a suitable basis of modes, were used as input into the two-timescale expansion [56] that has been used in the only existing calculations of second-order quantities [29,30,32].
The aim in Paper II was the same as the aim here: to convert expressions for the singular field written in Fermi-Walker coordinates into fully covariant expressions.This covariant expression can then be used as input into the previously mentioned puncture scheme.

D. Paper outline
We begin in Secs.II and III by recapping local expansion methods using bitensors; tensorial functions of two spacetime points; the construction of Fermi-Walker coordinates, and the conversion from Fermi-Walker coordinates to covariant form, as introduced by Paper II. Readers familiar with these concepts should feel free to skip directly to Sec.IV, where the covariant punctures for the metric perturbations in the highly regular gauge are derived.These are displayed in an abridged form due to their length, but the full expressions are provided in a Mathematica notebook in the Supplemental Material [70].
Section V then re-expands the covariant expressions from Sec. IV D into a generic coordinate expansion.The method for re-expanding the various covariant quantities is detailed in Sec.V A and, as before, readers familiar with this method can skip directly to Sec.V B where the generic coordinate expansions are presented.As with the covariant expressions, the coordinate punctures are too lengthy to include fully in this paper and are provided in the Supplemental Material [70].
Finally, we sum up the findings of this paper in Sec.VI and discuss potential future avenues for research.

E. Conventions and definitions
We use metric signature (−, +, +, +) and geometric units with c = G = 1.Indices using Greek letters run from 0 to 3 and with lowercase Latin letters run from 1 to 3. Greek/Latin indices are raised and lowered from the background metric, g µν , and the flat-space Euclidean metric, δ ab , respectively.
A primed index on a tensor, A µ ′ , indicates the tensor is evaluated at x ′µ := z µ (τ ), where z µ (τ ) are coordinates on the worldline, γ.An unprimed index on a tensor, A µ , is used for when the tensor is evaluated away from the worldline at x µ .An overset bar on a tensorial index, A μ, is used when a tensor is evaluated at xµ .This is a point on the worldline which is connected to x µ by an orthogonal geodesic.
A hat on a tensor, T a1...ai , refers to the symmetric trace-free (STF) part of the tensor with respect to the flat-space metric, δ ab .The covariant derivative is given by ∇ or a semi-colon and is compatible with the background metric, g µν .The partial derivative is given by ∂ or a comma.
We adopt notation from Ref. [71] for contractions of u µ ′ , σ µ ′ and ∆x µ ′ so that, for example.We use analogous notation for contractions of tensors evaluated at xµ , e.g.

II. LOCAL EXPANSION METHODS
In this section, we recap the methods of performing covariant and coordinate expansions of tensorial quantities near the worldline.We also give an overview of the construction of Fermi-Walker coordinates.

A. Covariant expansions using bitensors
In this section, we outline how one may construct local covariant expansions of tensor fields.Our explanation of the method follows that of Refs.[12,79,80].To do this, we introduce the concept of a bitensor: a tensor which is a function of two spacetime points.One important bitensor that we will make extensive use of is Synge's world function [12,79], where β is the unique geodesic connecting x µ and x µ ′ , s is an affine parameter and ε = ∓1 for time/spacelike geodesics (not to be confused with the mass ratio ϵ).This gives half the geodesic distance squared between the points x µ and x µ ′ .If the two points are connected by a null geodesic, then σ(x, x ′ ) is identically zero.We will use λ as a formal order counting parameter to count powers of spatial distance away from the worldline, γ, so that σ ∼ λ 2 .We denote derivatives of Synge's world function as Note also that we may take derivatives of Synge's world function at the unprimed coordinates as well, giving ).This can be generalised to higher and higher derivatives, e.g.
The indices of σ tell us its tensorial structure at both x µ and x µ ′ , that is, σ µ ′ ν ′ is a rank-2 tensor at x µ ′ but a scalar at x µ .Likewise, σ µ ′ ν is a covector at both x µ and x µ ′ .This property demonstrates that we can always commute primed and unprimed indices as the existence of one does not affect the tensorial rank at the other point.Derivatives of Synge's world function also satisfy the useful identity By taking derivatives of Eq. (37) and then the limit as x µ goes to x µ ′ , one may derive local covariant expansions of σ α ′ ...α... in terms of quantities defined on the worldline.
To see an example, we start by introducing the standard notation for the coincidence limit [79], It immediately follows from Eqs. ( 36)-( 37) that as, if the length of β goes to 0, then the integral in Eq. ( 36) vanishes.Taking primed derivatives of Eq. ( 37), we see which implies that This can be repeated to find higher and higher derivatives of σ(x, x ′ ) [80], Another object we will make use of is the parallel propagator, g µ ′ µ (x, x ′ ) [12,79,80].The parallel propagator parallel transports a tensor from x µ ′ to x µ along β.For instance, the vector A µ (x) can be transported from/to respectively.These expressions hold for covectors as well and tensors with any number of indices with the inclusion of an appropriate number of parallel propagators, e.g.
It also has the properties that when contracted with itself, it returns the Kronecker delta, and is symmetric in indices and arguments, When contracted with Synge's world function, it gives and its derivative contracted with Synge's world function vanishes for all combinations of primed and unprimed indices, e.g.
As we did for Synge's world function with Eq. ( 37), we can calculate different covariant expansions by repeatedly differentiating Eq. ( 52) and taking the coincidence limit.
For example [80], Combining the previous definitions, we can then express an arbitrary tensor A µ ν , evaluated at x, in terms of quantities evaluated at x ′ as where λ is a formal order counting parameter to be set to unity at the end of the calculation.The unknown coefficients, A (N )µ ′ ν ′ α ′ 1...α ′ n , can be found in the same manner as before by repeated differentiation and taking of the coincidence limit.As an example, we seek the covariant expansion of σ µ ′ ν ′ .We first expand, as in Eq. ( 56) but without the need for parallel propagators, as We know from Eq. ( 41), that A (0) Taking primed derivatives and the coincidence limit gives that σ (1) (2) meaning that This can be repeated for any required covariant quantity.
Ref. [81] provides a semi-recursive method for calculating expansions of Synge's world function and the parallel propagator, along with many other covariant quantities.

B. Fermi-Walker coordinates
To analyse the properties of the fields near the worldline of the small object, we introduce Fermi-Walker coordinates, (t, x a ), attached to the accelerated worldline, γ.Our description of Fermi-Walker coordinates summarises that of Refs.[12,82].To begin, we introduce an orthonormal tetrad, (u µ , e µ a ), on γ which is defined at the point z(τ ) so that it satisfies where u µ = dz µ /dτ is the curve's four-velocity, a µ = D 2 z µ / dτ 2 is the acceleration of γ and δ ab = diag(1, 1, 1) is the three-dimensional flat space metric.If γ is a geodesic then a µ vanishes.Equation (61) ensures that the tetrad basis is Fermi-Walker transported along γ, thus keeping it orthogonal to the worldline as it travels along it.This condition reduces to that of parallel transport when the worldline is a geodesic.Equations ( 62)-( 64) then ensure that it is orthonormal at all points on γ.The dual tetrad, (e 0 µ , e a µ ), can be defined as satisfying Equations ( 62)-( 66) then imply that we can write the metric and inverse metric as respectively.
With the orthonormal tetrad constructed, we may now create a local coordinate system so that we may derive the form of the metric near γ.The full technical details are not considered here (see Ref. [12, Chs.9.3-9.5]for more details) but we outline the geometric picture of the coordinate construction.At a point x := z(t) on γ, where t is the proper time, we generate a surface orthogonal to the worldline by emitting spacelike geodesics from z(t) that are orthogonal to γ.We can then label a point on this surface with coordinates x a so that we have coordinates, (t, x a ), that describe points near to the worldline.The tetrad can be written in terms of Synge's world function as As stated previously, Synge's world function gives half the geodesic distance squared between two points (up to a minus sign) meaning that a derivative gives the geodesic distance.This quantity is then contracted with the spatial Fermi-Walker tetrad leg, e a ᾱ, to give the Fermi-Walker spatial distance, x a .The third equation ensures that σ ᾱ is always orthogonal to the worldline.Alternatively, we can write x i = rn i , with r := δ ab x a x b = 2σ(x, x) being the proper distance (along a unique spacelike geodesic orthogonal to γ) from γ to the point being considered and n i being a unit vector giving the direction that the point lies in respective to γ.We note as well that, as with σ α ′ , r ∼ λ and so counts powers of distance from the worldline.A geometric representation of the Fermi-Walker coordinate construction is given in Fig. 2.
Using these coordinates, we can write the metric near γ in the form [50] where all Riemann terms are evaluated on γ at time t.When evaluating Eq. ( 72) on γ, we immediately see that the metric in Fermi-Walker coordinates reduces to the Minkowski metric.However, the Christoffel symbols at FIG. 2. Visualisation of construction of Fermi-Walker coordinates.At the point z(t), we generate an orthogonal surface and label points on that surface with the coordinate x i .The quantity r gives the proper distance to x i and n i picks out the unique orthogonal geodesic that connects x i and γ.Based on Fig. 6 from Ref. [12].
lowest order are not all zero.Instead, Γ t ta | γ = a a and Γ a tt | γ = a a ; both reduce to 0 if γ is a geodesic.As we are looking at a vacuum solution with R µν = 0, we may use the identities from Appendix D3 of Ref. [83] to write and the derivatives as The quantities E and B are the tidal moments felt by an extended body moving on the world line, γ, where two/three indices refer to the quadrupole/octopole moments respectively.They are symmetric and trace-free, with respect to δ ab , over all indices and only depend on the proper time, t.

III. CONVERTING FERMI-WALKER COORDINATES TO COVARIANT FORM
In this section we review the method used in Paper II to derive the covariant Lorenz gauge puncture.While the full technical details containing derivations of the various quantities are contained within that paper, we reproduce the essential results that we will need to produce the highly regular gauge puncture.The final results will be covariant quantities expressed entirely in terms of parallel propagators, the four-velocity, Riemann tensors, and Synge's world function.
The idea behind the method from Paper II is to express the field at a point x in terms of an arbitrary nearby point on the worldline, x ′ = z(τ ′ ).This is done through an intermediary point, x = z(τ ), which lies on γ and is separated from x ′ by the difference in proper time The intermediary point, x, is then connected to x by the unique geodesic that intersects the worldline orthogonally.A visual representation is provided in Fig. 3.As Fermi-Walker coordinates are constructed geometrically, see Sec.II B, there is a very straightforward way to convert them into covariant form.We know from Eqs. ( 69)- (71), that there is a simple correspondence between Fermi-Walker coordinates and covariant quantities.As we saw in the text below Eq. ( 71), we can write the Fermi-Walker radial distance in terms of covariant quantities with where We have added an extra step in Eq. ( 76), where we have rewritten the flat-space metric in terms of the projection operator, which immediately follows from Eq. ( 68).The radial unit vector is then given by Additionally, we must replace the Fermi-Walker basis one-forms, as when written explicitly, the singular field has the standard form These are given in Eqs. ( 82)-( 84) from Paper II by where Finally, the second-order singular field h SR µν features derivatives of the first-order regular field, h R1 µν .Using Eqs. ( 122)-(123) of Paper II, these can be written as where the bar, |, indicates a covariant derivative at x ᾱ and any acceleration terms can be ignored as they would belong to the third-order singular field.These expressions can be derived by taking covariant derivatives of h R1 ᾱ β and calculating the Christoffel symbols constructed from the FW background metric in Eq. (72).
After rewriting all quantities in terms of x, we then reexpand them in powers of ∆τ , the time difference given in Eq. (75).For example, where d dτ ′ = u α ′ ∇ α ′ and the expansion in distance of the difference in proper time is given by ∆τ originally from Eqs. ( 97)-(98) in Paper II.Here, λ is our formal order-counting parameter from Sec. II A, and we have introduced the quantity, and below we will also use the quantity, for notational simplicity. 4This means that the contraction of Synge's world function with itself can be written as Here, r gives a notion of the difference in proper time while ρ denotes a difference in proper distance.
We note that we expand all quantities (such as Eqs.( 89)-( 90)) through four total orders, but we only display the leading two orders here to indicate the forms of the expressions; the full expansions can be found in Paper II.We may do our series expansions as a normal power series as all the Fermi-Walker quantities (including one-forms) are scalars at x.The expansion of Synge's world function is given by Eqs. ( 99)-(101) of Paper II as and expansions of the Fermi-Walker basis one-forms are then given by Eqs. ( 103)-( 106) of Paper II as In the above expressions, we see that acceleration terms have appeared.This is a result of taking the derivatives with respect to τ ′ .As stated, d/dτ ′ = u α ′ ∇ α ′ , so taking multiple τ ′ derivatives results in us taking derivatives of u α ′ along the worldline, providing us with acceleration terms.These can then be differentiated along the worldline, giving us terms like ȧα ′ , where a dot indicates a time derivative in the usual manner.
When accounting for these terms, at first order, we split up h S1 µν into an acceleration-independent and a linear-in-acceleration piece: Recall from Eqs. ( 3) and ( 14) that each acceleration term carries an ϵ.This effectively makes h S1a µν a second-order term and allows us to ignore any non-linear acceleration terms that appear in the expansion of h S1 µν .Additionally, we can ignore any explicit acceleration terms that appear in both h SR µν and h SS µν as these would become third-order terms.

IV. CREATING THE COVARIANT PUNCTURE
With the methods from Paper II recapped, we can now proceed to use them to generate our covariant puncture in the highly regular gauge.We begin in Sec.IV A by reviewing the form of the metric perturbations in the highly regular gauge.Section IV B will provide the components of the highly regular gauge singular field when evaluated at x with each being written in covariant form.We then move to Sec.IV C, which provides the components evaluated at x ′ before combining this with one-form expansions to find the final, fully covariant form in Sec.IV D.

A. Metric perturbations in the highly regular gauge
In this section, we review the main results from Paper I. All results in this section are from there but are reproduced here for convenience.
We write the metric perturbations in the highly regular gauge as where the singular field is given by The second-order singular field is then split as where the "singular times singular" piece, h SS µν , features all terms proportional to m 2 and the "singular times regular" piece, h SR µν , features all terms with the form mh R1 µν .The full expressions for the first-order singular field in the highly regular gauge are given in Eq. ( 56) of Paper I. We reproduce the two leading orders here: Moving to second order, h SR µν is given in full by Eq. ( 130) of Paper I. The two leading orders are Finally, h SS µν is given by Eq. ( 131) of Paper I, which we reproduce here in full as

B. Perturbation components expanded about x
We begin by calculating the form of the components of the first-order singular field, h S1 µν , when expanded around xα .To do so, we substitute the appropriate expressions from Sec. III into Eq.( 101).The components of h S1 μν are then given by We have omitted the highest-order piece of h S1 μν due to its length, but it will be used to calculate the covariant punctures.This can then be continued at second order for the singular fields h SR μν (102) and h SS μν (103).The "singular times regular" piece is given by As in the expression for h S1 μν , we omit the highest-order piece of h SR μν due to length constraints.Finally, the "singular times singular" piece is given by

C. Expansion about x ′
Accounting for the introduction of acceleration terms and splitting up h S1 µν as in Eq. ( 97), we find that the components of h S1 ¡ a µν , when expanded around x ′ α , are given by The acceleration terms that appear as a result of our expansion of the first-order singular field are As h S1a µν is a second-order term, we can neglect any terms of order-λ 2 and higher to match the orders required for h SR µν and h SS µν .
Moving to the second-order field, we calculate the SR components to be where, again, we have omitted the highest order term.
The SS components are calculated to be

D. Final expressions for the covariant punctures
With all of the individual components of the singular field now expressed as functions of x ′ α , we now combine them with the expansions of dt and dx a , given in Eqs. ( 95)-( 96) to find the final form of the covariant punctures.After contracting with the basis vectors, we obtain the covariant form of h S µν dx µ dx ν , as in Eq. ( 80).We then read off the coefficients of dx µ dx ν to obtain h S µν .The first-order singular field is given by We have confirmed that this satisfies the Einstein field equations to the appropriate order, i.e.
At second order, the SS piece of the singular field is given by This again satisfies the appropriate Einstein field equations, The first-order singular field with linear acceleration terms is while the SR piece of the second-order singular field is These need to satisfy We have successfully checked that the covariant punc-tures for h SR µν and h S1,a µν satisfy Eq. ( 117) through the leading two orders, λ −3 and λ −2 .However, we have not been able to verify this at the highest order we have calculated, order λ −1 .This is due to the complexity and length of the expressions when taking multiple different combinations of derivatives.Despite this, we provide all orders of the covariant punctures for the different singular field terms in a Mathematica notebook in the Supplemental Material [70].
Comparing the covariant puncture for h S1 µν from Eq. ( 111) to the Lorenz gauge version of the puncture from Eq. (127) of Paper II, we see that the highly regular gauge puncture has a more complicated form.This continues at higher order with the Lorenz gauge puncture being substantially simpler and shorter at all orders.The more complex form results from the highly regular gauge conditions that seek to preserve the background light cone structure emanating from the worldline in the perturbed spacetime; see Sec.I C for further discussion.This has the knock-on effect that the coordinate expansion in the highly regular gauge will be much more complicated than the Lorenz gauge one as we are introducing more and more terms, and more quantities will need to be expanded.Thus, if we wanted to perform a mode decomposition of the singular field in the highly regular gauge, we would find that the process is likely to be more complicated than in the Lorenz gauge due to an increase in the number of quantities that need to be decomposed into modes.However, we believe that the benefits of the highly regular gauge outweigh any disadvantages that may come from the metric perturbations having a more complicated structure.Merely eliminating the two leading orders of h SS µν in Eq. ( 113) has dramatic consequences as it alleviates the problem of infinite mode coupling [61] that was discussed in the introduction.This should allow one to much more efficiently calculate modes of the second-order source.

V. COORDINATE EXPANSION
In order to implement the covariant expansions in a specific calculation, one must first write them in a chosen coordinate system.This necessitates re-expanding all the covariant quantities in terms of coordinate differences, where ∆x α ′ ∼ λ.A derivative of ∆x α ′ at x µ ′ then gives This leaves us with coefficients evaluated at x µ ′ , as in Eq. ( 56), contracted into certain combinations of ∆x α ′ .
here, we have a slightly different definition for ∆x α ′ and we take the derivatives at x µ ′ instead of x µ .Taking the primed derivative of the appropriate quantities and then substituting these and Eq.(123) into Eq.( 122) gives us the final expression for the coordinate expansion of Synge's world function, where the first four orders are given by σ (1) ) To check these expressions, one can substitute Eq. ( 125) into Eq.( 37) to demonstrate they satisfy the identity for Synge's world function.
We also require the expansions of r and ρ from Eqs. ( 91)-(92) which can be performed by substituting in Eqs. ( 124)-(125).The expression for r is trivial as it just requires us to contract the four-velocity into Eq.( 124), so that, at leading order, where, in analogy with Eq. ( 91), we define the fourvelocity contracted with the coordinate difference as We write the expansion of ρ as a power series, and define We then proceed to substitute our coordinate expansion for σ α ′ from Eq. (124) into the definition for ρ from Eq. (128) and collect terms at each order in λ.The first four orders of the expansion are given by ρ (1) = ρ 0 , (130a) ) To calculate the coordinate expansion of g ν ′ µ , we proceed in a similar way to that of σ α ′ .To begin, we use the ansatz and substitute this into the identity for the derivative of the parallel propagator contracted into a derivative of Synge's world function from Eq. ( 52).We proceed to solve this order-by-order to find As with Eq. ( 123), similar expansions of the parallel propagator have been done previously in Eqs.(3.10)-(3.12) of Ref. [87].We have checked our expressions by substituting them into Eq.( 52) and have verified that they satisfy the identity to the appropriate order in λ.

B. Coordinate expansions of the covariant punctures
With our covariant punctures derived, we can proceed to write them as a generic coordinate expansion using the techniques discussed for the singular scalar field in Sec.V A. This will allow them to be easily written in any desired coordinate system.
To do so, we substitute our coordinate expansion for σ α ′ from Eqs. ( 124 are written in terms of the coordinate difference, ∆x µ ′ and the four-velocity, u µ ′ along with h R1 µ ′ ν ′ , Γ µ ′ ν ′ ρ ′ , and R α ′ β ′ µ ′ ν ′ and their respective derivatives.The final expressions are incredibly long and, as such, we only display them through order λ 0 (except for h SR µν , for which we just display the leading-order term).The higher order terms are available in the Supplemental Material in a Mathematica notebook [70].

VI. CONCLUSION AND APPLICATIONS
The main result of this paper is the conversion of the local coordinate form of the metric perturbations given in Paper I into fully covariant form using the methods of Paper II.These were provided in truncated form in Sec.IV D and in full form in the Mathematica notebook in the Supplemental Material [70].
We have then re-expanded these covariant expressions and written them as a generic coordinate expansion that is valid in any desired coordinate system.As with the covariant expressions, abridged forms were presented in Sec.V B, with the full expressions appearing in the Supplemental Material [70].By providing the metric perturbations in these forms, we have enabled them to be written in any desired coordinate system without necessitating the use of a potentially complicated coordinate transformation from Fermi-Walker coordinates.
One useful immediate extension of this work would be to calculate the modes of the punctures to see how well the highly regular gauge alleviates the problem of infinite mode coupling.For quasicircular orbits in Schwarzschild, for example, one could decompose the punctures into modes using the methods of Ref. [66].From this, one could use the mode coupling formula from Eq. ( 28) to explicitly calculate the behaviour of the second-order Einstein tensor near to the worldline of the small object.
An interesting property of the highly regular gauge to note is that, following from the gauge conditions given in Sec.I C, one can write the singular field metric perturbations in terms of null vectors.For example, if one defines so that k α k α = 0, one can write the first-order singular field from Eq. (111) as One can then write Eq. ( 1) in terms of these null vectors as where V = 2m/(ρ).This has the form of a Kerr-Schild perturbation [88,89] on the background spacetime.However, this correspondence is broken in the singular field at order λ through the introduction of Riemann tidal terms in h S1 µν .Additionally, h R1 µν k µ k ν ̸ = 0 due to the regular field being in a generic gauge.
It would be interesting to further explore the connection between the highly regular gauge and Kerr-Schild gauges, potentially drawing on previous work by Harte [90] and Harte and Vines [91].

FIG. 3 .
FIG. 3. Diagram illustrating the relationship between x, x ′ and x.The two points x ′ and x are points on the worldline, γ, separated by ∆τ while x and x are connected by the geodesic that intersects γ orthogonally.Based on Fig. 1 from Paper II.