Motion of localized sources in general relativity: gravitational self-force from quasilocal conservation laws

An idealized"test"object in general relativity moves along a geodesic. However, if the object has a finite mass, this will create additional curvature in the spacetime, causing it to deviate from geodesic motion. If the mass is nonetheless sufficiently small, such an effect is usually treated perturbatively and is known as the gravitational self-force due to the object. This issue is still an open problem in gravitational physics today, motivated not only by basic foundational interest, but also by the need for its direct application in gravitational-wave astronomy. In particular, the observation of extreme-mass-ratio inspirals by the future space-based detector LISA will rely crucially on an accurate modeling of the self-force driving the orbital evolution and gravitational wave emission of such systems. In this paper, we present a novel derivation, based on conservation laws, of the basic equations of motion for this problem. They are formulated with the use of a quasilocal (rather than matter) stress-energy-momentum tensor---in particular, the Brown-York tensor---so as to capture gravitational effects in the momentum flux of the object, including the self-force. Our formulation and resulting equations of motion are independent of the choice of the perturbative gauge. We show that, in addition to the usual gravitational self-force term, they also lead to an additional"self-pressure"force not found in previous analyses, and also that our results correctly recover known formulas under appropriate conditions. Our approach thus offers a fresh geometrical picture from which to understand the self-force fundamentally, and potentially useful new avenues for computing it practically.


I. INTRODUCTION
A. Gravitational waves and extreme-mass-ratio inspirals The recent advent of gravitational wave astronomypropelled by the ground-based direct detections achieved by the LIGO/Virgo Collaboration (see Ref. [1] for the detections during the O1 and O2 observing runs), the success of the LISA Pathfinder mission as a proof of principle for future space-based interferometric detectors [2,3], and the subsequent approval of the LISA mission for launch in the 2030s [4,5]-has brought a multitude of both practical and foundational problems to the foreground of gravitational physics today. While a plethora of possibilities for gravitational wave sources are actively being investigated theoretically and anticipated to become accessible observationally, both on the Earth as well as in space, the most ubiquitous class of such sources has manifestly been-and foreseeably will remain-the coalescence of compact object binaries [6,7]. These are two-body systems consisting of a pair of compact objects, say of masses M 1 and M 2 , orbiting and eventually spiraling into each other. Each of these is, usually, either a black hole (BH) or a neutron star. There are also more general possibilities being investigated, including that of having a brown dwarf as one of the objects [8].
The LIGO/Virgo detections during the first scientific runs [1], O1 and O2, have all involved binaries of stellarmass compact objects (SCOs) located in our local neighborhood. These have comparable masses, of the order of a few tens of solar masses each (M 1 ∼ M 2 ∼ 10 0−2 M ⊙ ). In addition second-and third-generation terrestrial detectors can also eventually see intermediate-mass-ratio inspirals, binaries consisting of an intermediate-mass BH, of 10 2−4 M ⊙ , and a SCO. While there is as yet no direct evidence for the existence of the former sorts of objects, there are good reasons to anticipate their detection (through gravitational waves) most likely at the centers of globular clusters, and their study provides an essential link to the strongly perturbative regime of compact object binary dynamics.
It is even further in this direction that future space-based detectors such as LISA are anticipated to take us. In particular, LISA is expected to see extreme-mass-ratio inspirals (EMRIs) [9], compact binaries where M 1 ≫ M 2 . An elementary sketch is depicted in Fig. 1. The more massive object could be a (super) massive black hole (MBH) of mass M 1 ¼ M ∼ 10 4−7 M ⊙ located at a galactic center, with the significantly less massive object-effectively orbiting and eventually spiraling into the MBH-being a SCO: either a stellar-mass black hole or a neutron star, with M 2 ¼ m ∼ 10 0−2 M ⊙ .
Average estimates indicate that LISA will be able to see on the order of hundreds of EMRI events per year [10], with an expectation of observing, for each, thousands of orbital cycles over a period on the order of one year before the final plunge [11]. The trajectories defining these cycles and the gravitational wave signals produced by them will generally look much more complex than the relatively generic signals from mergers of stellar-mass black holes of comparable masses as observed, for example, by LIGO/Virgo. EMRIs will therefore offer an ideal experimental milieu for strong gravity: the complicated motion of the SCO around the MBH will effectively "map out" the geometry-that is, the gravitational field-around the MBH, thus presenting us with an unprecedented opportunity for studying gravity in the very strong regime [10,12]. In particular, among the possibilities offered by EMRIs are the measurement of the mass and spin of the MBH to very high accuracy, testing the validity of the Kerr metric as the correct description of BHs within general relativity (GR), and testing GR itself as the correct theory of gravity.
Yet, the richness of the observational opportunities presented by EMRIs comes with an inexorable cost: that is, a significant and as yet ongoing technical challenge in their theoretical modeling. This is all the more pressing as the EMRI signals expected from LISA are anticipated to be much weaker than the instrumental noise of the detector. Effectively, what this means is that extremely accurate models are necessary in order to produce the waveform templates that can be used to extract the relevant signals from the detector data stream. At the theoretical level, the problem of EMRI modeling cannot be tackled directly with numerical relativity (used for the LIGO/Virgo detections), simply due to the great discrepancy in (mass/length) scales; however, for the same reason, the approach that readily suggests itself is perturbation theory. See Fig. 2 for a graphic depicting the main methods used for compact object binary modeling in the different regimes. In particular, modeling the strong gravity, extreme-mass-ratio regime turns out to be equivalent to a general and quite old  The main approaches used in practice for the modeling of compact object binaries as a function of the mass ratio (increasing from 1) and the spacetime curvature involved. For low curvature (high separation between the bodies), post-Newtonian and post-Minkowskian methods are used. For high curvature (low separation) and low mass ratio, numerical relativity is used. For high curvature and extreme mass ratios, as the scale of a numerical grid would have to span orders of magnitude thus rendering it impracticable, perturbation theory must be used-in particular, self-force methods. problem which can be posed in any (not just gravitational) classical field theory: the so-called self-force problem.

B. Self-force problem
Suppose we are dealing with a theory for a field ψðxÞ in some spacetime. If the theory admits a Lagrangian formulation, we can usually assume that the field equations have the general form where L is a (partial, possibly nonlinear and typically second-order) differential operator, and we refer to S as the source of the field ψ. Broadly speaking, the problem of the self-force is to find solutions ψðxÞ satisfying (1) when S is "localized" in spacetime. Intuitively, it is the question of how the existence of a dynamical (field-generating) "small object" (a mass, a charge, etc.) backreacts upon the total field ψ, and hence in turn upon its own future evolution subject to that field. Thus, an essential part of any detailed self-force analysis is a precise specification of what exactly it means for S to be localized. In standard approaches, one typically devises a perturbative procedure whereby S ends up being approximated as a distribution, usually a Dirac delta, compactly supported on a worldline-that is, the background (zeroth perturbative order) worldline of the small object. However, this already introduces a nontrivial mathematical issue: if L is nonlinear (in the standard partial differential equation sense), then the problem (1) with a distributional source S is mathematically ill defined, at least within the classical theory of distributions [13] where products of distributions do not make sense [14]. 1 One might therefore worry that nonlinear physical theories, such as GR, would a priori not admit solutions sourced by distributions, and we refer the interested reader to Ref. [20] for a classic detailed discussion of this topic. The saving point is that, while the full field equation (in this case, the Einstein equation) may indeed be generally nonlinear, if we devise a perturbative procedure (where the meaning of the perturbation is prescribed in such a way as to account for the presence of the small object itself), then the first-order field equation is, by construction, linear in the (first-order) perturbation δψ of ψ. Thus, assuming the background field is a known exact solution of the theory, it always makes sense to seek solutions δψ to for a distributional source S, where δL indicates the firstorder part of the operator L in the full field equation (1). As this only makes sense for the (linear) first-order problem, such an approach becomes again ill defined if we begin to ask about the (nonlinear) second-or any higher-order problem. Additional technical constructions are needed to deal with these, the most common of which for the gravitational self-force has been the so-called "puncture" (or "effective source") method [11,[21][22][23]; similar ideas have proven to be very useful also in numerical relativity [24,25]. For work on the second-order equation of motion for the gravitational self-force problem, see e.g., Refs. [26][27][28]. For now, we assume that we are interested here in the first-order self-force problem (2) only. Now, concretely, in GR, our physical field ψ is simply the spacetime metric g ab (where latin letters from the beginning of the alphabet indicate spacetime indices), and following standard convention, we denote a first-order perturbation thereof by δg ab ¼ h ab . The problem (2) is then just the first-order Einstein equation, where G ab is the Einstein tensor, κ ¼ 8π (in geometrized units c ¼ G ¼ 1) is the Einstein constant, and T PP ab is the energy-momentum tensor of a "point particle" (PP) compactly supported on a given worldline C ∘ . We will return later to discussing this more precisely, but in typical approaches, C ∘ turns out to be a geodesic-that is, the "background motion" of the small object, which is in this case a small mass. 2 Thus, simply solving (3) for h ab assuming a fixed C ∘ for all time, though mathematically well defined, is by itself physically meaningless; it would simply give us the metric perturbations caused by a small object eternally moving on the same geodesic. Instead, what we would ultimately like is a way to take into account how h ab modifies the motion of the small object itself. Thus, in addition to the field equation (3), any self-force analysis must be supplemented by an equation of motion (EoM) telling us, essentially, how to move from a given 1 Nonlinear theories of distributions are being actively investigated by mathematicians [15][16][17]. Some work has been done to apply these to the electromagnetic self-force problem [18] and to study their general applicability in GR [19]; however, at this point, to our knowledge, their potential usefulness for the gravitational self-force problem has not been contemplated to any significant extent. 2 The problem of deriving geodesic motion for appropriately defined nondynamical "test particles" from the Einstein equation (in lieu of postulating it as an independent axiom of GR additional to the Einstein equation) is a long-standing and interesting issue in its own right. Einstein was involved in some of the earliest work on this [29], and over the decades, various proofs have been put forward outside of the context of the gravitational self-force problem. See Refs. [30,31] for some of the most famous such proofs. See also Ref. [32] for a recent general review of the most widely used approaches as well as an interesting novel proposal. We consider later in this paper in detail one approach to the gravitational self-force which also proves geodesic motion as the background motion of point particles in GR.

MOTION OF LOCALIZED SOURCES IN GENERAL …
PHYS. REV. D 101, 064060 (2020) 064060-3 background geodesic C ∘ at one step in the (ultimately numerical) time evolution problem to a new background geodesic C ∘ 0 at the next time step-with respect to which the field equation (3) is solved anew, and so on. This is sometimes referred to as a "self-consistent" approach. See Fig. 3 for a visual depiction. The first proposal for an EoM for the gravitational selfforce (GSF) problem was put forward in the late 1990s, since known as the MiSaTaQuWa equation after its authors [34,35]. On any C ∘ , it reads The lhs is a second (proper) time derivative of a deviation vector Z a on C ∘ pointing in the direction of the "true motion" (away from C ∘ ), to be defined more precisely later.
On the rhs, E ∘ ab is the electric part of the Weyl tensor on C ∘ , such that the first term is a usual "geodesic deviation" term. The second term on the rhs is the one usually understood as being responsible for self-force effects: F a ½·; · is a fourvector functional of a symmetric rank-2 contravariant tensor and a vector, to which we refer in general (for any arguments) as the GSF functional. In any spacetime with a given metric g ∘ ab and compatible derivative operator ∇ ∘ a , it is explicitly given by the following simple action of a first-order differential operator: While this is easy enough to calculate once one knows the arguments, the main technical challenge in using the MiSaTaQuWa equation (4) lies precisely in the determination thereof: in particular, h tail ab is not the full metric perturbation h ab which solves the field equation (3), but instead represents what is called the "tail" integral of the Green functions of h ab [36]. This quantity is well defined, but difficult to calculate in practice and usually requires the fixing of a perturbative gauge-typically the Lorenz gauge, Physically, h tail ab can be thought of as the part of the full perturbation h ab which is scattered back by the spacetime curvature. (In this way, h ab can be regarded as the sum of h tail ab and the remainder, what is sometimes called the "instantaneous" or "direct" part h direct ab , responsible for waves radiated to infinity [37].) An alternative, equivalent GSF EoM was proposed by Detweiler and Whiting in the early 2000s [38]. It relies upon a regularization procedure for the metric perturbations, i.e., a choice of a decomposition for h ab [the full solution of the field equation (3)] into the sum of two parts: one which diverges-in fact, one which contains all divergent contributions-on C ∘ , denoted h S ab (the so-called singular field, related to the direct part of the metric perturbation), and a remainder which is finite, h R ab (the so-called regular field, related to the tail part), so that one writes h ab ¼ h S ab þ h R ab . An analogy with the self-force problem in electromagnetism gives some physical intuition behind how to interpret the meaning of this decomposition [11], with h S ab ∼ m=r having the heuristic form of a "Coulombian self-field." However, no procedure is known for obtaining the precise expression of h S ab in an arbitrary perturbative gauge, and moreover, once a gauge is fixed (again, usually the Lorenz gauge), this splitting is not unique [11]. Nevertheless, if and when such an h S ab is obtained (from which we thus also get h R ab ¼ h ab − h S ab ), the Detweiler-Whiting EoM for the GSF reads The EoM's (4) and (6) are equivalent in the Lorenz gauge and form the basis of the two most popular methods used today for the numerical computation of the GSF. Yet a great deal of additional technical machinery is required for handling gauge transformations. This is essential because, FIG. 3. A depiction of the perturbative problem for the GSF. In particular, this represents one of the most popular conceptions of a so-called self-consistent approach [33]: at a given step (on a given Cauchy surface) in the time evolution problem, one computes the "correction to the motion" away from geodesic (C ∘ ) in the form of a deviation vector Z a , determined by the GSF. Then, at the next time step, one begins on a new ("corrected") geodesic (C in the EMRI problem, the background spacetime metricthat of the MBH-is usually assumed to be Schwarzschild-Droste or Kerr. Perturbation theory for such spacetimes has been developed and is most easily carried out in, respectively, the so-called Regge-Wheeler and radiation gauges; in other words, in practice, it is often difficult (though not infeasible-see, e.g., Ref. [39]) to compute h ab directly in the Lorenz gauge for use in (4) or (6).
A proposal for an EoM for the GSF problem that is valid in a wider class of perturbative gauges was presented by Gralla in 2011 [40]. It was therein formulated in what are called "parity-regular" gauges, i.e., gauges satisfying a certain parity condition. This condition ultimately has its origins in the Hamiltonian analysis of Regge and Teitleboim in the 1970s [41], wherein the authors introduce it in order to facilitate the vanishing of certain surface integrals and thus to render certain general-relativistic Hamiltonian notions, such as multipoles and "center of mass," well defined mathematically. In parity-regular gauges (satisfying the Regge-Teitleboim parity condition), the Gralla EoM-mathematically equivalent, in the Lorenz gauge, to the MiSaTaQuWa and the Detweiler-Whiting EoM's-is The GSF (last) term on the rhs is obtained in this approach by essentially relating the deviation vector Z a (the evolution of which is expressed by the lhs) with a gauge transformation vector and then performing an "angle average" over a "small" r-radius two-sphere S 2 r , with ϵ S 2 the volume form of the unit two-sphere, of the so-called bare GSF, F a ½h bc ; U ∘ d . The latter is just the GSF functional [Eq. (5)] evaluated directly using the full metric perturbatiuon h ab (i.e., the tail plus direct parts, or equivalently the "regular" plus "singular" parts), around (rather than at the location of) the distributional source. The point therefore is that this formula never requires the evaluation of h ab on C ∘ itself, where it is divergent by construction; instead, away from C ∘ it is always finite, 3 and (7) says that it suffices to compute the GSF functional (5) with h ab directly in the argument (requiring no further transformations), and integrate it over a small sphere.
The manifest advantage of (7) relative to (4) or (6) is that no computations of tail integrals or regularizations of the metric perturbations are needed at all. Yet, to our knowledge, there has thus far been no attempted numerical computation of the GSF using this formula. One of the issues with this remains that of the perturbative gauge: depending upon the detailed setup of the problem, one may still not easily be able to compute h ab directly in a parityregular gauge (although manifestly, working in the parityregular "no-string" radiation gauge [42] may be useful for a GSF calculation in Kerr), i.e., a gauge in which (7) is strictly valid, and so further gauge transformations may be needed. Aside from the practical issues with a possible numerical implementation of this, there is also a conceptual issue: this formula originates from an essentially mathematical argument-by a convenient "averaging" over the angles-so as to make it well defined in a Hamiltonian setting via a relation to a canonical definition of the center of mass. Yet its general form as a closed two-surface integral suggestively hints at the possibility of interpreting it not merely as a convenient mathematical relation, but as a real physical flux of (some notion of) "gravitational momentum." We contend and will demonstrate in this paper that indeed an even more general version of (7), not restricted by any perturbative gauge choice (so long as one does not construct it in such a way that produces divergences in h away from C ∘ ), results from the consideration of momentum conservation laws in GR.

C. Self-force problem via conservation laws
The idea of using conservation laws for tackling the selfforce problem was appreciated and promptly exploited quite early on for the electromagnetic self-force. In the 1930s [43], Dirac was the first to put forward such an analysis in flat spacetime, and later on in 1960 [44], DeWitt and Brehme extended it to nondynamically curved spacetimes. 4 In such approaches, it can be shown 5 that the EoM for the electromagnetic self-force follows from local conservation expressions of the form where the lhs expresses the flux of matter four-momentum P a (associated with T ab ) between the "caps" of (i.e., closed spatial two-surfaces delimiting) a portion (or "time interval") of a thin worldtube boundary ℬ (topologically R × S 2 ), with natural volume form ϵ ℬ and (outwarddirected) unit normal n a (see Fig. 4). In particular, one takes a time derivative of (8) to obtain an EoM expressing 3 This is true provided one does not transform to a perturbative gauge wherein any of the multipole moments of h diverge away from C ∘ as a consequence of the gauge definition. Examples of gauges leading to such divergences are the "half-string" and "fullstring" radiation gauges of Ref. [42], which exhibit stringlike singularities in h along a radial direction. Nevertheless, in this work, it was also shown that one can define a no-string radiation gauge which is in the parity-regular class, and where the singularities in h remain only on C ∘ , thus rendering the integral (5) well defined. 4 By this, we mean spacetimes with nonflat but fixed metrics, which do not evolve dynamically (gravitationally) in response to the matter stress-energy-momentum present therein. 5 See Ref. [45] for a basic and more contemporary presentation.
the time rate of change of momentum in the form of a closed spatial two-surface integral (by differentiating the worldtube boundary integral). For the electromagnetic selfforce problem, the introduction of an appropriate matter stress-energy-momentum tensor T ab into Eq. (8) and a bit of subsequent argumentation reduces the integral expression to the famous Lorentz-Dirac equation; on a spatial three-slice in a Lorentz frame and in the absence of external forces, for example, this simply reduces to _ P i ¼ 2 3 q 2 _ a i for a charge q.
Along these lines, see also Ref. [46] for a treatment of the scalar and electromagnetic self-forces based on the notion of a limit of an extended charge. Formulations of the scalar and electromagnetic self-forces using generalized Killing fields have more recently been put forward in Refs. [47,48]. Then in Ref. [49], the idea of generalized Killing fields was also used to study the gravitational self-force problem by employing definitions of momenta as bulk (volume) integrals (over the worldtube interior), involving the internal matter stress-energy-momentum tensor of the moving object (mass). This in fact proved the validity of the Detweiler-Whiting approach, i.e., viewing the "point mass" as moving on a geodesic in an effective (regularized) metric, to all multiple orders and also for the angular momentum. As we will shortly elaborate, we posit that gravitational momentum in GR cannot fundamentally be regarded as a bulk (local) density, and a more detailed analysis to reveal how the (locally formulated) results of Ref. [49] relate to our work here would be very interesting.
The success of conservation law approaches for formulating the electromagnetic self-force in itself inspires hope that the same may be done in the case of the GSF problem. In particular, Gralla's EoM (7) strongly hints at the possibility of understanding the rhs not just as a mathematical (angle averaging) device, but as a true, physical flux of gravitational momentum arising from a consideration of conservation expressions.
Nevertheless, to our knowledge, there has thus far been no proposed general treatment of the GSF following such an approach. This may, in large part, be conceivably attributed to the notorious conceptual difficulties surrounding the very question of the basic formulation of conservation laws in GR. Local conservation laws, along the lines of Eq. (8) that can readily be used for electromagnetism, no longer make sense fundamentally once gravity is treated as dynamical. The reason has a simple explanation in the equivalence principle [50]: one can always find a local frame of reference with a vanishing local "gravitational field" (metric connection coefficients), and hence a vanishing local "gravitational energy-momentum," irrespective of how one might feel inclined to define the latter. 6 A wide variety of approaches have been taken over the decades toward formulating sensible notions of gravitational energy-momentum, with still no general consensus among relativists today on which to qualify as "the best" [52,53]. Often the preference for employing certain definitions over others may simply come down to context or convenience, but in any case, there exist agreements between the most typical definitions in various limits. A very common feature among them is the idea of replacing a local notion of gravitational energy-momentum, i.e., energy-momentum as a volume density, with what is known as a quasilocal energy-momentum, i.e., energymomentum as a boundary density. The typical Hamiltonian definitions of the (total) gravitational energy-momentum for an asymptotically flat spacetime, for example, are of such a form. Among the most commonly used generalizations of these definitions to arbitrary (finite) spacetime FIG. 4. A worldtube boundary ℬ (topologically R × S 2 ) in ℳ, with (outward-directed) unit normal n a . The change in matter four-momentum between two constant time slices of this worldtube is given by the flux of the normal projection (in one index) of the matter stress-energy-momentum tensor T ab through the portion of ℬ bounded thereby. 6 It is worth remarking here that, in a perturbative setting, an approach that is sometimes taken is to work with an "effective" local gravitational stress-energy-momentum tensor, defined as the rhs of a suitably rearranged (first-order) Einstein equation. This is a common tactic often used for studying, for example, the energy-momentum of gravitational waves, with some applicable caveats (see, e.g., Chapter 35 of Ref. [50]). In fact, one of the first formulations of the gravitational self-force-in particular, the derivation of the MiSaTaQuWa EoM (4) presented in Sec. III of Ref. [34]-made use of the local conservation (vanishing of the spacetime divergence) of a suitably defined local tensor of such a sort (in analogy with the approach of DeWitt and Brehme to the electromagnetic self-force [44]). We elaborate in the remainder of this subsection and at greater length in Sec. II on why such a notion of gravitational conservation principles, while demonstrably useful for operational computations in some situations, cannot in general be expected to capture the fundamentally quasilocal (boundary density) nature of gravitational energymomentum. See Ref. [51] for a detailed discussion and comparison between these two (local and quasilocal) views of gravitational energy-momentum. regions was proposed in the early 1990s by Brown and York [54], and follow from what is now eponymously known as the Brown-York stress-energy-momentum tensor. It is a quasilocal tensor, meaning it is only defined on the boundary of an arbitrary spacetime region. For example, using this, the total (matter plus gravitational) energy inside a spatial volume is given up to a constant factor by the closed two-surface (boundary) integral of the trace of the boundary extrinsic curvature-precisely in agreement with the Hamiltonian definition of energy for the entire spacetime in the appropriate limit (where the closed two-surface approaches a two-sphere at asymptotically flat spatial infinity) but, in principle, applicable to any region in any spacetime.
The formulation of general energy-momentum conservation laws in GR from the Brown-York tensor has been achieved in recent years with the use of a construction called quasilocal frames [51], a concept first proposed in Ref. [55]. Essentially, the idea is that, in order to describe a (finite) gravitational system of interest and its associated boundary fluxes of gravitational energy-momentum, it does not suffice to merely specify, as in the local matter conservation laws of the form of Eq. (8), a worldtube boundary ℬ (the interior of which is to be regarded as containing that system, and across which its energymomentum flux is measured) as an embedded submanifold of ℳ. What is in fact required is the specification of a congruence making up this worldtube boundary, i.e., a twoparameter family of timelike worldlines with some chosen four-velocity field representing the motion of a topological two-sphere's worth of quasilocal observers. We will motivate this construction in greater amplitude shortly, but the reason for needing it is basically to be able to meaningfully define "time-time" and "time-space" directions on ℬ for our conservation laws. A congruence of this sort is what is meant by a quasilocal frame.
The enormous advantage in using these quasilocal conservation laws over other approaches lies in the fact that they hold in any arbitrary spacetime. Thus, the existence of Killing vector fields-a typical requirement in other conservation law formulations-is in no way needed here. This idea has been used successfully in a number of applications so far [51,[56][57][58][59][60][61]. These include the resolution of a variation of Bell's spaceship paradox 7 in which a box accelerates rigidly in a transverse, uniform electric field [56], recovering under appropriate conditions the typical (but more limited) local matter conservation expressions of the form of Eq. (8) from the quasilocal ones [51] and application to post-Newtonian theory [58] and to relativistic geodesy [60,61].
A similar idea to quasilocal frames, called gravitational screens, was proposed more recently in Refs. [64,65].
There, the authors also make use of quasilocal ideas to develop conservation laws very similar in style and form to those obtained via quasilocal frames. A detailed comparison between these two approaches has thus far not been carried out, but it would be very interesting to do so in future work. In particular, the notion of gravitational screens has been motivated more from thermodynamic considerations, and similarly casting quasilocal frames in this language could prove quite fruitful. For example, just as these approaches have given us operational definitions of concepts like the "energy-momentum in an arbitrary spacetime region" (and not just for special cases such as an entire spacetime), they may help to do the same for concepts like "entropy in an arbitrary spacetime region" (and not just for known special cases such as a black hole).

D. Executive summary of the paper
We now summarize the structure and main result of this paper.
Section II is entirely devoted to an overview of quasilocal frames and quasilocal conservation laws, in complete self-contained technical detail for our purposes here.
In Sec. III, we prove the main general result of this paper from the quasilocal momentum conservation law: in particular, we show that for any localized gravitational system, the change between some initial and final time in its total linear momentum 8 p, along a spatial direction determined by a vector ϕ (the precise meaning of which is to be elaborated later), is given by the following flux through the worldtube boundary interval Δℬ bounded by these times (see Fig. 4): Here, F a ½·; · is an extended GSF functional. In particular, it is the usual GSF functional F a [see Eq. (5)] plus a novel piece to which we refer as the gravitational self-pressure force ℘ a , arising from a quasilocal pressure effect (also to be elaborated later): The first argument of this functional F a , just as in Gralla's formula [our Eq. (7)], is the metric perturbation h ab on ℬ. This avoids any potential singularities in inside the worldtube, i.e., within the spacetime volume of which ℬ is the (exterior) boundary, and therefore the need for performing regularizations or any further transformations. The second argument of F a in our result [see Eq. (10)], i.e., u ∘ a , is the four-velocity not of any background geodesic contained inside the worldtube, as in the typical GSF EoM's discussed earlier-indeed, strictly speaking, the main result [Eq. (10)] holds without necessarily having to even introduce any such geodesic, or more generally, without having to say anything specific about the content of the worldtube interior-but instead that of the background quasilocal observers on ℬ itself. Manifestly, our formula [Eq. (10)] bears significant resemblance to that of Gralla [Eq. (7)], and it is the scope of Sec. IV to compare the two in the appropriate setting. For this, we introduce the general setup of the Gralla-Wald approach to the GSF, and apply our conservation law formula for two choices of quasilocal frames: one which is inertial with the SCO in the perturbed spacetime (and hence not inertial with the geodesic-following point particle in the background) and one which is inertial with the geodesicfollowing point particle in the background (and hence not inertial with the SCO in the perturbed spacetime). We derive the EoM's in both of these cases and discuss their correspondence to the known GSF EoM's.
Finally, Sec. V concludes the paper.

E. Notation and conventions
We work in the ð−; þ; þ; þÞ signature of spacetime. Script uppercase letters (A; ℬ; C; …) are reserved for denoting mathematical spaces (manifolds, curves, etc.). The n-dimensional Euclidean space is denoted as usual by R n , the n-sphere of radius r is denoted by S n r , and the unit n-sphere is denoted by S n ¼ S n 1 . For any two spaces A and ℬ that are topologically equivalent (i.e., homeomorphic), we indicate this by writing A ≃ ℬ.
The set of ðk; lÞ-tensors on any manifold U is denoted by J k l ðUÞ. In particular, TU ¼ J 1 0 ðUÞ is the tangent bundle, and T Ã U ¼ J 0 1 ðUÞ is the dual thereto. Any ðk; lÞ-tensor in any ð3 þ 1Þ-dimensional (Lorentzian) spacetime ℳ is equivalently denoted either using the (boldface) index-free notation A ∈ J k l ðℳÞ following the practice of, e.g., Refs. [50,66], or the abstract index notation A a 1 ÁÁÁa k b 1 ÁÁÁb l ∈ J k l ðℳÞ following that of, e.g., Ref. [67]; that is, depending upon convenience, we equivalently write with latin letters from the beginning of the alphabet (a; b; c; …) being used for spacetime indices (0,1,2,3). The components of A in a particular choice of coordinates fx α g 3 α¼0 are denoted by A α 1 ÁÁÁα k β 1 ÁÁÁβ l , using greek (rather than latin) letters from the beginning of the alphabet (α; β; γ; …). Spatial indices on an appropriately defined (three-dimensional Riemannian spacelike) constant time slice of ℳ are denoted using latin letters from the middle third of the alphabet in roman font: in lowercase (i; j; k; …) if they are abstract and in uppercase (I; J; K; …) if a particular choice of coordinates fx I g 3 I¼1 has been made. For any n-dimensional manifold ðU; g U ; ∇ U Þ with metric g U and compatible derivative operator ∇ U , we denote its natural volume form by Let I ≃ S 2 be any (Riemannian) closed two-surface that is topologically a two-sphere. Latin letters from the middle third of the alphabet in fraktur font (i; j; k; …) are reserved for indices of tensors in J k l ðIÞ. In particular, for S 2 itself, S ij is the metric, D i is the associated derivative operator, and ϵ S 2 ij is the volume form; in standard spherical coordinates fθ; ϕg, the latter is simply given by Contractions are indicated in the usual way in the abstract index notation: e.g., U a V a is the contraction of U and V. Equivalently, when applicable, we may simply use the "dot product" in the index-free notation, e.g., U a V a ¼ U · V, A ab B ab ¼ A∶ℬ, etc. We must keep in mind that such contractions are to be performed using the metric of the space on which the relevant tensors are defined. Additionally, often we find it convenient to denote the component (projection) of a tensor in a certain direction by simply replacing its pertinent abstract index therewith: e.g., we equivalently write For any (0,2)-tensor A ab that is not a metric, we indicate its trace (in nonboldface) as A ¼ A a a ¼ trðAÞ. (If it is a metric, this notation is reserved, as usual, for its determinant rather than its trace.) Finally, let U and V be any two diffeomorphic manifolds and let f∶U → V be a map between them. This naturally defines a map between tensors on the two manifolds, which we denote by f Ã ∶J k l ðUÞ → J k l ðVÞ and its inverse ðf −1 Þ Ã ¼ f Ã ∶J k l ðVÞ → J k l ðUÞ if it exists. We generically refer to any map of this sort as a tensor transport [68]. It is simply the generalization to arbitrary tensors of the pushforward f Ã ∶TU → TV and pullback f Ã ∶T Ã V → T Ã U, the action of which is defined in the standard way-see, e.g., Appendix C of Ref. [67]. (Note that our convention for sub-/superscripting the star is the generally more common one used in geometry [68,69]; it is sometimes opposite to and sometimes congruous with that used in the physics literature, e.g., Refs. [67] and [70], respectively).
holds. In what follows, we introduce the concept of quasilocal frames [51,[55][56][57][58][59][60][61] and describe the basic steps for their construction, as well as the energy and momentum conservation laws associated therewith. In Sec. II A, we offer an heuristic idea of quasilocal frames before proceeding in Sec. II B to present the full mathematical construction. Then, in Sec. II C, we motivate and discuss the quasilocal stress-energy-momentum tensor used in this work, that is, the Brown-York tensor. Finally, in Sec. II D, we review the formulation of quasilocal conservation laws using these ingredients.

A. Quasilocal frames: Heuristic idea
Before we enter into the technical details, we would like to offer a heuristic picture and motivation for defining the concept of quasilocal frames.
We would like to show how the GSF arises from generalrelativistic conservation laws. For this, we require first the embedding into our spacetime ℳ of a worldtube boundary ℬ ≃ R × S 2 . The worldtube interior contains the system of which the dynamics we are interested in describing. In principle, such a ℬ can be completely specified by choosing an appropriate radial function rðxÞ on ℳ and setting it equal to a non-negative constant [such that the rðxÞ ¼ const > 0 Lorentzian slices of ℳ have topology R × S 2 ]. This would be analogous to defining a (Riemannian, with topology R 3 ) Cauchy surface by the constancy of a time function tðxÞ on ℳ.
However, this does not quite suffice. As we have briefly argued in the Introduction (and will shortly elaborate upon in greater technicality), the conservation laws appropriate to GR ought to be quasilocal in form, that is, involving stress-energy-momentum as boundary (not volume) densities. One may readily assume that the latter are defined by a quasilocal stress-energy-momentum tensor living on ℬ, which we denote-for the moment, generally-by τ ab . (Later, we give an explicit definition, namely, that of the Brown-York tensor, for τ.) To construct conservation laws, then, one would need to project this τ into directions on ℬ, giving quantities such as energy or momenta, and then to consider their flux through a portion of ℬ (an interval of time along the worldtube boundary). But in this case, we have to make clear what is meant by the energy (time-time) and momenta (time-space) components of τ within ℬ, the changes in which we are interested in studying. For this reason, additional constructions are required.
In particular, what we need is a congruence of observers with respect to which projections of τ yield stress-energymomentum quantities. Since τ is only defined on ℬ, this therefore needs to be a two-parameter family of (timelike) worldlines, the union of which is ℬ itself. This is analogous to how the integral curves of a time flow vector field (as in canonical GR) altogether constitute ("fill up") the entire spacetime ℳ, except that there we are dealing with a three-(rather than two-)parameter family of timelike worldlines.
We refer to any set of observers, the worldlines of which form a two-parameter family constituting ℬ ≃ R × S 2 , as quasilocal observers. A specification of such a two-parameter family, equivalent to specifying the unit four-velocity u a ∈ Tℬ of these observers (the integral curves of which "trace out" ℬ), is what is meant by a quasilocal frame.
With this, we can now meaningfully talk about projections of τ into directions on ℬ as stress-energy-momentum quantities. For example, τ uu may appear immediately suggestible as a definition for the (boundary) energy density. Indeed, later, we take precisely this definition, and we will furthermore see how momenta (the basis of the GSF problem) can be defined as well.

B. Quasilocal frames: Mathematical construction
Concordant with our discussion in the previous subsection, a quasilocal frame (see Fig. 5 for a graphical illustration of the construction) is defined as a twoparameter family of timelike worldlines constituting the worldtube boundary (topologically R × S 2 ) of the history FIG. 5. A portion of a quasilocal frame ðℬ; uÞ in a spacetime ℳ, bounded by constant t two-surfaces I i and I f . In particular, ℬ ≃ R × S 2 is the union of all integral curves (two-parameter family of timelike worldlines), depicted in the figure as dotted red lines, of the vector field u ∈ Tℬ which represents the unit fourvelocity of quasilocal observers making up the congruence. The unit normal to ℬ (in ℳ) is n, and the normal to each constant t slice I of ℬ isũ (not necessarily coincidental with u). Finally, H (with induced metric σ) is the two-dimensional subspace of Tℬ consisting of the spatial vectors orthogonal to u. Note that unlike I, H need not be integrable (indicated in the figure by the failure of H to make a closed two-surface). of a finite (closed) spatial three-volume in ℳ. Let u a denote the timelike unit vector field tangent to these worldlines. Such a congruence constitutes a submanifold of ℳ that we call ℬ ≃ R × S 2 . Let n a be the outward-pointing unit vector field normal to ℬ; note that n is uniquely fixed once ℬ is specified. There is thus a Lorentzian metric γ [of signature ð−; þ; þÞ] induced on ℬ, the components of which are given by We denote the induced derivative operator compatible therewith by D. To indicate that a topologically R × S 2 submanifold ðℬ; γ; DÞ of ℳ is a quasilocal frame (that is to say, defined as a particular congruence with four-velocity u as detailed above, and not just as an embedded submanifold) in ℳ, we write ðℬ; γ; D; uÞ or simply ðℬ; uÞ.
Let H be the two-dimensional subspace of Tℬ consisting of the "spatial" vectors orthogonal to u. Let σ denote the two-dimensional (spatial) Riemannian metric [of signature ðþ; þÞ] that projects tensor indices into H, and is induced on ℬ by the choice of u (and thus also n), given by The induced derivative operator compatible with σ is denoted by D. Let fx i g 2 i¼1 (written using fraktur indices from the middle third of the latin alphabet) be spatial coordinates on ℬ that label the worldlines of the observers, and let t be a time coordinate on ℬ such that surfaces of constant t, to which there exists a unit normal vector that we denote byũ a ∈ Tℬ, foliate ℬ by closed spatial twosurfaces I (with topology S 2 ). Letting N denote the lapse function of g, we have u ¼ N −1 ∂=∂t.
Note that in general, H need not coincide with the constant time slices I. Equivalently, u need not coincide withũ. In general, there will be a shift between them, such thatũ where v a represents the spatial two-velocity of fiducial observers that are at rest with respect to I as measured by our congruence of quasilocal observers (the four-velocity of which is u), is the Lorentz factor. The specification of a quasilocal frame is thus equivalent to making a particular choice of a two-parameter family of timelike worldlines comprising ℬ. There are, a priori, three degrees of freedom (DoF's) available to us for doing this. Heuristically, these can be regarded as corresponding to the three DoF's in choosing the direction of u-from which n and all induced quantities are then computable. (Note that u has four components, but one of the four is fixed by the normalization requirement u · u ¼ −1, leaving three independent direction DoF's.) Equivalently, we are in principle free to pick any three geometrical conditions (along the congruence) to fix a quasilocal frame. In practice, usually it is physically more natural, as well as mathematically easier, to work with geometric quantities other than u itself to achieve this.
Yet, it is worth remarking that simply writing down three desired equations (or conditions) to be satisfied by geometrical quantities on ℬ does not itself guarantee that, in general, a submanifold ðℬ; γ; DÞ obeying those three particular equations will always exist-and, if it does, that it will be the unique such submanifold-in an arbitrary ðℳ; g; ∇Þ. Nevertheless, one choice of quasilocal frame that is known to always exist (a claim we will qualify more carefully in a moment) is that where the two-metric σ on H is "rigid" (or "time independent")-these are called rigid quasilocal frames.
Most of the past work on quasilocal frames has in fact been done in the rigid case [51,[55][56][57][58][59]. We know, however, that other quasilocal frame choices are also possible, such as geoids-dubbed geoid quasilocal frames [60,61]; these are the general-relativistic generalization of "constant gravitational potential" surfaces in Newtonian gravity. Regardless, the quasilocal frame choice that we will mainly consider in this paper is the rigid one (and we will be clear when this choice is explicitly enacted).
Intuitively, the reason for this preference is that imposing in this way the condition of "spatial rigidity" on ðℬ; uÞ-a two-dimensional (boundary) rigidity requirement, which, unlike three-dimensional rigidity, is permissible in GReliminates from the description of the system any effects arising simply from the motion of the quasilocal observers relative to each other. Thus, the physics of what is going on inside the system (i.e., the worldtube interior) is essentially all that affects its dynamics.
Technically, there is a further reason: a proof of the existence of solutions-i.e., the existence of a submanifold ℬ ≃ R × S 2 in ℳ that is also a quasilocal frame ðℬ; uÞfor any spacetime ðℳ; g; ∇Þ has up to now only been fully carried out for rigid quasilocal frames. 9 While, as we have commented, other quasilocal frame choices may be generally possible in principle (and may be shown to be possible to construct, case-by-case, in specific spacetimes-as we have done, e.g., with geoid quasilocal frames [60,61]), they are as yet not rigorously guaranteed to exist in arbitrary spacetimes.
The quasilocal rigidity conditions can be stated in a number of ways. Most generally, defining to be the strain rate tensor of the congruence, they amount to the requirement of vanishing expansion θ ¼ trðθÞ and shear θ habi ¼ θ ðabÞ − 1 2 θσ ab , i.e., 9 The idea of the proof is to explicitly construct the solutions order by order in an expansion in the areal radius around an arbitrary worldline in an arbitrary spacetime [57].
In the adapted coordinates, these three conditions are expressible as the vanishing of the time derivative of the two-metric on H, i.e., 0 ¼ ∂ t σ. Both of these two equivalent mathematical conditions, θ ðabÞ ¼ 0 ¼ ∂ t σ, capture physically the meaning of the quasilocal observers moving rigidly with respect to each other (i.e., the "radar-ranging" distances between them does not change in time).

C. Quasilocal stress-energy-momentum tensor
Before we consider the formulation of conservation laws with the use of quasilocal frames (from which our analysis of the GSF will eventually emerge), we wish to address in a bit more detail an even more fundamental question: what are conservation laws in GR actually supposed to be about? At the most basic level, they should express changes (over time) in some appropriately defined notion of energymomentum. As we are interested in gravitational systems (and specifically, those driven by the effect of the GSF), this energy-momentum must include that of the gravitational field, in addition to that of any matter fields if present.
Hence, we may assert from the outset that it does not make much sense in GR to seek conservation laws based solely on the matter stress-energy-momentum tensor T, such as Eq. (8). It is evident that these would, by construction, account for matter only-leaving out gravitational effects in general (which could exist in the complete absence of matter, e.g., gravitational waves), and thus the GSF in particular. What is more, such conservation laws are logically inconsistent from a general-relativistic point of view: a nonvanishing T implies a nontrivial gravitational field (through the Einstein equation) and thus a necessity of taking into account that field along with the matter one(s) for a proper accounting of energy-momentum transfer. A further technical problem is also that the formulation of conservation laws of this sort is typically predicated upon the existence of Killing vector fields or other types of symmetry generators in ℳ, which one does not have in general-and which do not exist in spacetimes pertinent for the GSF problem in particular.
We are therefore led to ask how we can meaningfully define a total-gravity plus matter-stress-energy-momentum tensor in GR. It turns out that the precise answer to this question, while certainly not intractable, is unfortunately also not unique-or at least, it lacks a clear consensus among relativists, even today. See, e.g., Refs. [52,53] for reviews of the variety of proposals that have been put forward toward addressing this question. Nonetheless, for reasons already touched upon and to be elaborated presently, what is clear and generally accepted is that such a tensor cannot be local in nature (as T is) and for this reason is referred to as quasilocal.
Let τ ab denote this quasilocal, total (matter plus gravity) stress-energy-momentum tensor that we eventually seek to use for our conservation laws. It has long been understood [50] that whatever the notion of gravitational energymomentum (defined by τ) might mean, it is not something localizable; in other words, there is no way of meaningfully defining an "energy-momentum volume density" for gravity. This is, ultimately, due to the equivalence principle: locally, one can always find a reference frame in which all local gravitational fields (the connection coefficients), and thus any notion of energy-momentum volume density associated therewith, disappear. The remedy is to make τ quasilocal, meaning that, rather than volume densities, it should define surface densities (of energy, momentum, etc.)-a type of construction which is mathematically realizable and physically sensible in general.
The specific choice we make for how to define this total (matter plus gravity), quasilocal energy-momentum tensor τ is the so-called Brown-York tensor, first put forward by the authors in Ref. [54]; see also Ref. [71] for a detailed review. This proposal was based originally upon a Hamilton-Jacobi analysis; here, we will offer a simpler argument for its definition, sketched out initially in Ref. [51].
Consider the standard gravitational action S G for a spacetime volume V ⊂ ℳ such that ∂V ¼ ℬ ≃ R × S 2 is a worldtube boundary as in the previous subsection (possibly constituting a quasilocal frame, but not necessarily). This action is given by the sum of two terms, a bulk and a boundary term, respectively: In particular, the first is the Einstein-Hilbert bulk term, and the second is the Gibbons-Hawking-York boundary term [72,73], Here, Additionally, the matter action S M for any set of matter fields Ψ described by a Lagrangian L M is The definition of the total (quasilocal) stress-energymomentum tensor τ for gravity plus matter can be obtained effectively in the same way as that of the (local) stressenergy-momentum tensor T for matter alone-from the total action in Eq. (21) rather than just, respectively, the matter action in Eq. (24). In particular, T is defined by computing the variation δ (with respect to the spacetime metric) of the matter action: In other words, one defines the matter stress-energymomentum tensor as the functional derivative, The definition of the Brown-York tensor follows completely analogously, except that now gravity is also included. That is, for the total action of gravity (minimally) coupled to matter, we have that the metric variation is In the equality of Eq. (28), Π is the canonical momentum of ðℬ; γ; DÞ, given by Π ¼ Θ − Θγ. It follows from direct computation using Eqs. (21), (24), and (25); for a review of this derivation carefully accounting for the boundary term see, e.g., Chapter 12 of Ref. [74]. In the equality of Eq. (29), the Einstein equation G ¼ κT has been invoked (in other words, we impose the Einstein equation to be satisfied in the bulk), thus leading to the vanishing of the bulk term; meanwhile in the boundary term, a gravity plus matter stress-energy-momentum tensor τ (the Brown-York tensor) has been defined in direct analogy with the definition of the matter energy-momentum tensor T in Eq. (25). Hence, just as Eq. (25) implies Eq. (26), Eq. (29) implies Henceforth, τ refers strictly to this (Brown-York) quasilocal stress-energy-momentum tensor of Eq. (30), and not to any other definition. It is useful to decompose τ in a similar way as is ordinarily done with T, so we define as the quasilocal energy, momentum, and stress, respectively, with units of energy per unit area, momentum per unit area, and force per unit length. Equivalently,

D. Conservation laws
The construction of general conservation laws from τ was first achieved in Refs. [51,56], and proceeds along the following lines. Let ψ ∈ Tℬ be an arbitrary vector field in ℬ. We begin by considering a projection of Π in the direction of ψ (in one index), i.e., Π ab ψ b , and computing its divergence in ℬ. By using the Leibnitz rule, we simply have Next, we integrate this equation over a portion Δℬ of ℬ bounded by initial and final constant t surfaces I i and I f , as depicted in Fig. 5. On the resulting lhs, we apply Stokes's theorem, and on the first term on the rhs, we use the Gauss-Codazzi identity: D a Π ab ¼ n a γ b c G ac . Thus, using the notation for tensor projections in certain directions introduced in Sec. I E for ease of readability (e.g., G ab n a ψ b ¼ G nψ and similarly for other contractions), we obtain Z where ϵ I denotes the volume form on the constant time closed two-surfaces I, and we have used the notation We also remind the reader that u represents the unit normal to each constant time closed two-surface, which in general need not coincide with the quasilocal observers' four-velocity u but is related to it by a Lorentz transformation, Eq. (18); see also Fig. 5.
We stress that, so far, Eq. (36) is a purely geometrical identity, completely general for any Lorentzian manifold ℳ; in other words, thus far, we have said nothing about physics. Now, to give this identity physical meaning, we invoke the definition of the Brown-York tensor in Eq. (30) (giving the boundary extrinsic geometry its meaning as stressenergy-momentum) as well as the Einstein equation [Eq. (15)], giving the spacetime curvature its meaning as the gravitational field. With these, Eq. (36) turns into On the lhs, we have inserted the relationũ ¼γðu þ vÞ, with v a representing the spatial two-velocity of fiducial observers that are at rest with respect to I (the hypersurfaceorthogonal four-velocity of which isũ) as measured by our congruence of quasilocal observers (the four-velocity of which is u), is the Lorentz factor. Observe that Eq. (37) expresses the change of some component of the quasilocal stress-energy-momentum tensor integrated over two different t ¼ const closed two-surfaces I as a flux through the worldtube boundary Δℬ between them. The identification of the different components of τ as the various components of the total energy-momentum of the system thus leads to the understanding of Eq. (37) as a general conservation law for the system contained inside of Δℬ. Thus, depending on our particular choice of ψ ∈ Tℬ, Eq. (37) will represent a conservation law for the total energy, momentum, or angular momentum of this system [51].
Let us now assume that ðℬ; uÞ is a rigid quasilocal frame. If we choose ψ ¼ u, then Eq. (37) becomes the energy conservation law, Z where α a ¼ σ ab a b is the H projection of the acceleration of the quasilocal observers, defined by a a ¼ ∇ u u a . Now suppose, on the other hand, that we instead choose ψ ¼ −ϕ where ϕ ∈ H is orthogonal to u (with the minus sign introduced for convenience), and represents a stationary conformal Killing vector field with respect to σ. This means that ϕ is chosen such that it satisfies the conformal Killing equation, L ϕ σ ¼ ðD · ϕÞσ, with L the Lie derivative and D the derivative on H (compatible with σ). A set of six such conformal Killing vectors always exists: three for translations and three for rotations, respectively generating the action of boosts and rotations of the Lorentz group on the two-sphere [51]. The idea, then, is that the contraction of these vectors with the quasilocal momentum integrated over a constant-time topological two-sphere boundary expresses, respectively, the total linear and angular momentum (in the three ordinary spatial directions each) at that time instant. Thus, Eq. (37) becomes the (respectively, linear and angular) momentum conservation law, Z where ν ¼ 1 2 ϵ ab H D a u b is the twist of the congruence (with ϵ H ab ¼ ϵ ℳ abcd u c n d the induced volume form on H) and P ¼ 1 2 σ∶S is the quasilocal pressure (force per unit length) between the worldlines of ℬ. We remark that the latter can be shown to satisfy the very useful general identity (which we will expediently invoke in our later calculations): An analysis of the gravitational self-force problem should consider the conservation law in Eq. (39) for linear momentum. Thus, we will use the fact, described in greater detail in Appendix A, that the conformal Killing vector ϕ ∈ H for linear momentum admits a multipole decomposition of the following form: with the dots indicating higher harmonics. Here, r is the area radius of the quasilocal frame (such that ℬ is a constant r hypersurface in ℳ), r I denotes the standard direction cosines of a radial unit vector in R 3 , and B i I ¼ ∂ i r I are the boost generators on the two-sphere. See Appendix A for a detailed discussion regarding conformal Killing vectors and the two-sphere. In spherical coordinates fθ; ϕg, we have r I ¼ ðsin θ cos ϕ; sin θ sin ϕ; cos θÞ. Thus, Eq. (41) gives us a decomposition of ϕ in terms of multipole moments, with the l ¼ 1 coefficients Φ I simply representing vectors in R 3 in the direction of which we are considering the conservation law.

III. GENERAL DERIVATION OF THE GRAVITATIONAL SELF-FORCE FROM QUASILOCAL CONSERVATION LAWS
In this section, we will show how the GSF is a general consequence of the momentum conservation law in Eq. (39) for any system which is sufficiently localized. By that, we mean something very simple: taking the r → 0 limit of a quasilocal frame around the moving object which is treated as small, i.e., as a formal perturbation about some background. No further assumptions are for the moment needed. In particular, we do not even need to enter into the precise details of how to specify the perturbation family for this problem; that will be left to the following section, where we will carefully define and work with the family of perturbed spacetimes typically employed for applications of the GSF.
We first review the basic formulation of perturbation theory in GR in Sec. III A. While this material is well known, we find it useful to include it here both for establishing notation as well as carefully defining the concepts that we need to work with at an adequate level of rigour. Then, in Sec. III B, we show that the first-order perturbation of the momentum conservation law in Eq. (39) always contains the GSF, and that it dominates the dynamics for localized systems.

A. Perturbation theory in GR
Our exposition of perturbation theory in this subsection follows closely the treatment of Ref. [75]. (See also Chapter 7 of Ref. [67] for a simpler treatment of this topic but following the same philosophy.) Perturbation theory in GR is best made sense of from the point of view of "stacked" manifolds off some known background. To be more precise, let λ ≥ 0 represent our perturbation parameter. It is a purely formal parameter, in the sense that it should be set equal to 1 at the end of any computation and serves only to indicate the order of the perturbation. The idea, then, is to define a one-parameter family of spacetimes fðℳ ðλÞ ; g ðλÞ ; ∇ ðλÞ Þg λ≥0 , where ∇ ðλÞ is the connection compatible with the metric g ðλÞ in ℳ ðλÞ , ∀ λ ≥ 0, such that ðℳ ð0Þ ; g ð0Þ ; ∇ ð0Þ Þ ¼ ðℳ Þ is a known, exact spacetime-the background. See Fig. 6 for a visual depiction. For notational convenience, any object with a subscripted "(0)" (from a one-parameter perturbative family) is equivalently written with an overset "∘" instead.
For the GSF problem, ℳ ∘ is usually the Schwarzschild-Droste or Kerr spacetime. Then, one should establish a way of smoothly relating the elements of this one-parameter family (between each other) such that calculations on any ℳ ðλÞ for λ > 0-which may be, in principle, intractable analytically-can be mapped to calculations on ℳ ∘ in the form of infinite (Taylor) series in λ-which, provided ℳ ∘ is chosen to be a known, exact spacetime, become tractable, order by order, in λ.
Thus, one begins by defining a (five-dimensional, Lorentzian) product manifold, the natural differentiable structure of which is given simply by the direct product of those on ℳ ðλÞ and the non-negative real numbers (labeling the perturbation parameter), R ≥ ¼ fλ ∈ Rjλ ≥ 0g. For any one-parameter family of ðk; lÞtensors fA ðλÞ g λ≥0 such that A ðλÞ ∈ J k l ðℳ ðλÞ Þ, ∀ λ ≥ 0, we define A ∈ J k l ðNÞ by the relation A α 1 ÁÁÁα k β 1 ÁÁÁβ l ðp; λÞ ¼ A α 1 ÁÁÁα k ðλÞ β 1 ÁÁÁβ l ðpÞ, ∀ p ∈ ℳ ðλÞ and ∀ λ ≥ 0. Henceforth, any such tensor living on the product manifold will be denoted in serif font-instead of roman font, which remains reserved for tensors living on ð3 þ 1Þ-dimensional spacetimes. Furthermore, any spacetime tensor (except for volume forms) or operator written without a sub-or superscripted ðλÞ lives on ℳ ∘ . Conversely, any tensor (except for volume forms) or operator living on ℳ ðλÞ , ∀ λ > 0, is indicated via a sub-or (equivalently, if notationally more convenient) superscripted ðλÞ; e.g., A ðλÞ ¼ A ðλÞ ∈ J k l ðℳ ðλÞ Þ is always tensor in ℳ ðλÞ . The volume form of any (sub)manifold U is always simply denoted by the standard notation ϵ U (and is always understood to live on U), as in Eq. (13).
Let Φ X ðλÞ ∶N → N be a one-parameter group of diffeomorphisms generated by a vector field X ∈ TN. (That is to say, the integral curves of X define a flow on N which connects any two leaves of the product manifold.) For notational convenience, we denote its restriction to maps from the background to a particular perturbed spacetime (identified by a particular value of λ > 0) as p ↦ φ X ðλÞ ðpÞ: ð45Þ The choice of X-equivalently, the choice of φ X ðλÞ -is not unique; there exists freedom in choosing it, and for this reason, X-equivalently, φ X ðλÞ -is referred to as the perturbative gauge. We may work with any different gauge do not need to render the issue of gauge specification explicit, we may drop the superscript and, instead of φ X ðλÞ , simply write φ ðλÞ .
Consider now the transport under φ X ðλÞ of any tensor A ðλÞ ∈ J k l ðℳ ðλÞ Þ from a perturbed spacetime to the background manifold. We always denote the transport of any such a tensor by simply dropping the ðλÞ sub-or superscript and optionally including a superscript to indicate the gauge-that is, ∀ A ðλÞ ∈ J k l ðℳ ðλÞ Þ, and similarly the transport of ∇ ðλÞ to ℳ ∘ is ∇. We know, moreover, that we can express any such A as a Taylor series around its background value, A ð0Þ ¼ A ∘ , in Lie derivatives along X (see Ref. [75]): where L denotes the Lie derivative; in the last equality, we have defined δ n A ¼ ð1=n!Þð∂ n λ AÞj λ¼0 , and so the (gaugedependent) first-order perturbation is Note that the symbol δ n , ∀ n, can be thought of as an operator δ n ¼ ð1=n!Þ∂ n λ j λ¼0 that acts upon and extracts the Oðλ n Þ part of any tensor in ℳ ∘ . So now, in particular, we have that the background value of g ¼ ðφ X ðλÞ Þ Ã g ðλÞ is g ∘ , and we denote its first-order perturbation for convenience and according to convention as h ¼ δg. Thus, we have where we have omitted explicitly specifying the gauge (X) dependence for now. Let us define one further piece of notation that we shall need to use: let Γ is in fact a tensor. Note that C ∘ ¼ 0, i.e., C ¼ λδC þ Oðλ 2 Þ. In particular, it is given by B. Gravitational self-force from the general momentum conservation law Let fðℬ ðλÞ ; u ðλÞ Þg λ≥0 be an arbitrary one-parameter family of quasilocal frames (defined as in Sec. II) each of which is embedded, respectively, in the corresponding element of the family of perturbed spacetimes fðℳ ðλÞ ; g ðλÞ ; ∇ ðλÞ Þg λ≥0 described in the previous subsection (Fig. 77). Consider the general geometrical identity (36)  For λ ¼ 0, this gives us our conservation laws in the background for λ > 0, and in a perturbed spacetime for λ > 0. It is the latter that we are interested in, but since we do not know how to do calculations in ℳ ðλÞ ∀ λ > 0, we have to work with Eq. (51) transported to ℳ ∘ . This is easily achieved by using the fact that for any diffeomorphism f∶U → V between two (oriented) smooth n-dimensional manifolds U and V and any (compactly supported) n-form ω in V, we have that as the inverse image of the worldtube boundary (quasilocal frame) in the background manifold, and using the fact that the tensor transport commutes with contractions, the above can simply be written in the notation we have established as Z So far, we have been completely general. Now, let us restrict our attention to the momentum conservation law (ψ ¼ −ϕ ∈ H) given by Eq. (53), and let us assume that we do not have any matter on Δℬ (hence, by the Einstein equation, G nϕ j Δℬ ¼ κT nϕ j Δℬ ¼ 0), or even simply that any matter if present there is subdominant to the linear perturbation, i.e., Tj Δℬ ¼ Oðλ 2 Þ. The lhs then expresses the change in momentum of the system (inside the worldtube interval in the perturbed spacetime) between some initial and final time slices; for notational ease, we will simply denote this by Δp ðϕÞ . [Note that we prefer to use typewriter font for the total quasilocal momentum, so as to avoid any confusion with matter four-momentum defined in the typical way from T ab and traditionally labeled by P a , as, e.g., in Eq. (8).] Then, inserting also the definition of the Brown-York tensor [Eq. (30)] on the rhs and replacing D with ∇ since it does not affect the contractions, Eq. (53) becomes We claim, and will now demonstrate, that the OðλÞ part of this always contains the GSF. Let us consider Eq. (54) term by term. First, we have the transport-in this case, the pullback-under φ ðλÞ of the volume form of ℬ ðλÞ . Now, we know that the pullback under a diffeomorphism of the volume form of a manifold is, in general, not simply the volume form of the inverse image of that manifold under the diffeomorphism. However, it is always true (see, e.g., Chapter 7 of Ref. [76]) that they are proportional, with the proportionality given by a smooth function called the Jacobian determinant and usually denoted by J. That is, in our case, we have φ Ã ðλÞ ϵ ℬ ðλÞ ¼ Jϵ ℬ , with J ∈ C ∞ ðℬÞ. In particular, this function is given by JðpÞ ¼ detðT p φ ðλÞ Þ, ∀ p ∈ ℬ, where T p φ ðλÞ ¼ ðφ ðλÞ Þ Ã ∶T p ℬ → T φ ðλÞ ℬ ðλÞ is the pushforward, and the determinant is computed with respect to the volume forms ϵ ℬ ðpÞ on T p ℬ and ϵ ℬ ðλÞ ðφ ðλÞ ðpÞÞ on T φ ðλÞ ðpÞ ℬ ðλÞ . Now, it is clear that we have J ¼ 1 þ OðλÞ, as φ ð0Þ is simply the identity map. Therefore, we have As for the other terms in the integrand of Eq. (54), we simply have Hence, we can see that there will be three contributions to the OðλÞ rhs of Eq. (54). Respectively, from Eqs. (55)-(57), these are the OðλÞ parts of the volume form pullback, which may not be easy to compute in practice; the Brown-York tensor τ, which may be computed from its definition [Eq. (30)]; and the derivative of the conformal Killing vector ϕ, which may be readily carried out and, as we will presently show, always contains the GSF. Thus, we denote this contribution to the OðλÞ part of Δp ðϕÞ as Δp ðϕÞ self , Now, we proceed with the computation of Eq. (58). In particular, let us consider the series expansion of Eq. (58) in the areal radius r of ℬ. This can be defined for any time slice by r ¼ ð 1 4π R I ϵ I Þ 1=2 , such that a constant r slice of ℳ ∘ defines ℬ (and n ¼ M∇ ∘ r for some positive function M on ℬ). It has been shown [57] that the Brown-York tensor has, in general, the following expansion in r, where are called the vacuum energy and vacuum pressure, respectively. Some remarks regarding these are warranted before we move on. In particular, these are terms which have sometimes been argued to play the role of "subtraction terms" (to be removed from the quasilocal energymomentum tensor); see, e.g., Ref. [71]. From this point of view, the definition of the Brown-York tensor [Eq. (30)] may be regarded as carrying a certain amount of freedom, inasmuch as any freedom may be assumed to exist to define a "reference" action S 0 to be subtracted from the total (gravitational plus matter) action S GþM in the variational principle discussed in Sec. II C. Such a subtraction of a reference action, while common practice in gravitational physics, has the sole function of shifting the numerical value of the action such that, ultimately, the numerical value of the Hamiltonian constructed from the modified action S GþM − S 0 may be interpreted as the Arnowitt-Deser-Misner (ADM) energy. However, this essentially amounts to a presumption that we are free to pick the zero of the energy-in other words, that the vacuum energy may be freely subtracted away without affecting the physics. Though we refrain from entering into much further detail here, it has been shown [51] that these vacuum terms, Eqs. (60) and (61), are in fact crucial for our conservation laws to yield physically reasonable answers and to make mathematical sense-evidencing that the vacuum energy/ pressure should be taken seriously as having physically real significance. We will now lend further credibility to this by showing that they are precisely the energy (and pressure) associated with the momentum flux that are typically interpreted as the GSF. Actually, we argue in this paper that the term implicating the vacuum energy yields the standard form of the GSF, and the vacuum pressure term is novel in our analysis.
Collecting all of our results so far-inserting Eqs. (59)- (62) into Eq. (58)-we thus get Let us now look at the contractions in the integrand. For the first (energy) term, inserting the connection coefficient (50), we have by direct computation where the functional F is precisely the GSF four-vector functional defined in the introduction [Eq. (5)], and to write the final equality, we have used the orthogonality property ϕ u ∘ ¼ 0. Thus, we see that this is indeed the term that yields the GSF. For the second (pressure) term in Eq. (63), we similarly obtain by direct computation where, in expressing the rhs, it is convenient to define a general functional of two (0,2)-tensors similar to the GSF functional: We call this novel term the gravitational self-pressure force. Now, we can collect all of the above and insert them into (63). Before writing down the result, it is convenient to define a total functional F as the sum of F and ℘, We refer to this as the extended GSF functional. Note that for F we write only the functional dependence on h and u ∘ since the two-metric σ ∘ is determined uniquely by u ∘ . With this, and setting the perturbation parameter to unity, Eq. (63) becomes This is to be compared with Gralla's formula [40] discussed in the Introduction, Eq. (7). While the equivalence thereto is immediately suggestive based on the general form of our result, we have to do a bit more work to show that indeed Eq. (69), both on the lhs and the rhs, recoversthough in general will, evidently at least from our novel gravitational self-pressure force, also have extra terms added to-Eq. (7). We leave this task to the following section, the purpose of which is to consider in detail the application of our conservation law formulation to a concrete example of a perturbative family of spacetimes defined for a self-force analysis, namely, the Gralla-Wald family.
Concordantly, we emphasize that the result above [Eq. (69)] holds for any family of perturbed manifolds fℳ ðλÞ g λ≥0 and is completely independent of the internal description of our system, i.e., the worldtube interior. In other words, what we have just demonstrated-provided only that one accepts a quasilocal notion of energymomentum-is that the (generalized) GSF is a completely generic perturbative effect in GR for localized systems; it arises as a linear-order contribution of any spacetime perturbation to the momentum flux of a system in the limit where its areal radius is small.
This view of the self-force may cast fresh conceptual light on the old and seemingly arcane problem of deciphering its physical origin and meaning. In particular, recall the common view that the GSF is caused by the backreaction of the "mass" of a small object upon its own motion. Yet, what we have seen here is that it is actually the vacuum mass, or vacuum energy that is responsible for the GSF. We may still regard the effect as a "backreaction," in the sense that it is the boundary metric perturbations of the system-the h on ℬ-which determine its momentum flux, but the point is that this flux is inexorably present and given by Eq. (69) regardless of where exactly this h is coming from. Presumably, the dominant part of h would arise from the system itself-if we further assume that the system itself is indeed what is being treated perturbatively by the family fℳ ðλÞ g λ≥0 , as is the case with typical selfforce analyses-but in principle, h can comprise absolutely any perturbations; i.e., its physical origin does not even have to be from inside the system.
In this way, we may regard the GSF as a completely geometrical, purely general-relativistic backreaction of the mass (and pressure) of the spacetime vacuum-not of the object inside-upon the motion of a localized system (i.e., its momentum flux). This point of view frees us from having to invoke such potentially ambiguous notions as "mass ratios" (in a two-body system for example), let alone "Coulombian m=r fields," to make basic sense of self-force effects. They simply-and always-happen from the interaction of the vacuum with any boundary perturbation, and are dominant if that boundary is not too far out.

IV. APPLICATION TO THE GRALLA-WALD APPROACH TO THE GRAVITATIONAL SELF-FORCE
In this section, we will consider in detail the application of our ideas to a particular approach to the self-force, that is

MOTION OF LOCALIZED SOURCES IN GENERAL …
PHYS. REV. D 101, 064060 (2020)  to say, a particular specification of fðℳ ðλÞ ; g ðλÞ Þg via a few additional assumptions aimed at encoding the notion of a small object being "scaled down" to zero "size" and mass as λ → 0. In other words, we now identify the perturbation (which has up to this point been treated completely abstractly) defined by fðℳ ðλÞ ; g ðλÞ Þg as actually being that caused by the presence of the small object; that could mean regular matter (in particular, a compact object such as a neutron star) or a black hole. The assumptions (on fg ðλÞ g) that we choose to work with here are those of the approach of Gralla and Wald [33]. Certainly, the application of our perturbed quasilocal conservation laws could just as well be carried out in the context of any other self-force analysis-such as, e.g., the self-consistent approximation of Pound [77] (the mathematical correspondence of which to the Gralla-Wald approach has, in any case, been shown in Ref. [78]).
Our motivation for starting with the Gralla-Wald approach in particular is twofold. On the one hand, it furnishes a mathematically rigorous and physically clear picture (which we show in Fig. 8)-arguably more so than any other available GSF treatment-of what it means to scale down a small object to zero size and mass (or, equivalently, of perturbing any spacetime by the presence of an object with small size and mass-we will be more precise momentarily). On the other hand, it is within this approach that the formula for the GSF has been obtained (in Ref. [40]) as a closed twosurface (small two-sphere) integral around the object (in lieu of evaluating the GSF at a spacetime point identified as the location of the object), in the form of the Gralla angle averaging formula [Eq. (7)]-with which our extended GSF formula (69) is to be compared.
In Sec. IVA, we provide an overview of the assumptions and consequences of the Gralla-Wald approach to the GSF.
Afterward, in Sec. IV B, we describe the general embedding of rigid quasilocal frames in the Gralla-Wald family of spacetimes, and then in Sec. IV C, we describe their detailed construction in the background spacetime in this family. Having established this, we then proceed to derive equations of motion in two ways. In particular, we carry out the analysis with two separate choices of rigid quasilocal frames ("frames of reference"): first, inertially with the point-particle approximation of the moving object in the background in Sec. IV D and, second, inertially with the object itself in the perturbed spacetime in Sec. IV E.

A. Gralla-Wald approach to the GSF
The basic idea of Gralla and Wald [33] for defining a family fðℳ ðλÞ ; g ðλÞ Þg λ≥0 such that λ > 0 represents the inclusion of perturbations generated by a small object is the following one. One begins by imposing certain smoothness conditions on fg ðλÞ g λ≥0 corresponding to the existence of certain limits of each g ðλÞ . In particular, two limits are sought corresponding intuitively to two limiting views of the system: first, a view from "far away" from which the "motion" of the (extended but localized) object reduces to a worldline and, second, a view from "close by" the object from which the rest of the Universe (and, in particular, the MBH it might be orbiting as in an EMRI) looks "pushed away" to infinity. A third requirement must be added to  Fig. 1 of Ref. [33].) The lined green region that "fills in" ℳ ðλÞ for r ≤ Cλ is the small object which scales down to zero size and mass in the background ℳ ∘ . The solid black lines represent taking the ordinary limit (the far away view where the motion appears reduced to a worldline), and the dashed black lines represent the scaled limit (the close by view where the rest of the Universe appears pushed away to infinity). The worldline C ∘ , which can be proven to be a geodesic, is parametrized by z ∘a ðτ ∘ Þ and has four-velocity U ∘ . The deviation vector Z on C ∘ is used for formulating the first-order correction to the motion. this, namely, that both of these limiting pictures nonetheless coexist in the same spacetime; i.e., the two limits are smoothly related (or, in other words, there is no pathological behavior when taking these limits along different directions). While in principle this may sound rather technical, one can actually motivate each of these conditions with very sensible physical arguments as we shall momentarily elaborate further upon. From them, Gralla and Wald have shown [33] that it is possible to derive a number of consequences, including geodesic motion in the background at zeroth order and the MiSaTaQuWa equation [34,35] for the GSF at first order in λ.
Let us now be more precise. Let fðℳ ðλÞ ; g ðλÞ Þg λ≥0 be a perturbative one-parameter family of spacetimes as in the previous section. We assume that fg ðλÞ g λ≥0 satisfies the following conditions, depicted visually in Fig. 8: (i) Existence of an "ordinary limit".-There exist coordinates fx α g in ℳ ðλÞ such that g ðλÞ βγ ðx α Þ is jointly smooth in ðλ; x α Þ for r > Cλ where C > 0 is a constant and r ¼ ðx i x i Þ 1=2 . For all λ ≥ 0 and r > Cλ, g ðλÞ is a vacuum solution of the Einstein equation.
Furthermore, g ∘ βγ ðx α Þ is smooth in x α including at r ¼ 0, and the curve C (iii) Uniformity condition.-Define A ¼ r, B ¼ λ=r, and n i ¼ x i =r. Then, each g ðλÞ βγ ðx α Þ is jointly smooth in ðA; B; n i ; tÞ. Mathematically, the first two conditions respectively ensure the existence of an appropriate Taylor expansion (in r and λ) of the metric in a "far zone" (on length scales comparable with the mass of the MBH in an EMRI, r ∼ M) and a "near zone" (on length scales comparable with the mass of the object, r ∼ m). Meanwhile, the third is simply a consistency requirement ensuring the existence of a "buffer zone" (m ≪ r ≪ M) where both expansions are valid. (This idea is in many ways similar to the method of "matched asymptotic expansions" [34].) From a physical point of view, what is happening in the first ("ordinary") limit is that the body is shrinking down to a worldline C ∘ with its mass (understood as defining the perturbation) going to zero at least as fast as its radius. (As we increase the perturbative parameter λ from zero, the radius is not allowed to grow faster than linearly with λ; viewed conversely, this condition ensures that the object does not collapse to a black hole if it was not one already before reaching the point-particle limit.) In the second (scaled) limit, the object is shrinking down to zero size in an asymptotically self-similar manner (its mass is proportional to its size, and its "shape" is not changing). Finally, the uniformity condition ensures that there are no "bumps of curvature" in the one-parameter family. (Essentially, this guarantees that there are no inconsistencies in evaluating the limits along different directions.) From these assumptions alone, Gralla and Wald [33] are able to derive the following consequences: (a) Background motion.-The worldline C ∘ is a geodesic in ℳ ∘ ; writing its parametrization in terms of proper time τ ∘ as C ∘ ¼ fz ∘ a ðτ ∘ Þg τ ∘ ∈R and denoting its four- i.e., the field equation is where Here, m is a constant along C ∘ and is interpreted as representing the mass of the object-or, more precisely, the mass of the point particle which approximates the object in the background. (This is a subtle point that should be kept in mind, and which will be better elucidated in our analysis further on.) (d) First-order equation of motion.-At OðλÞ, the correction to the motion in the Lorenz gaugecorresponding to the choice of a certain gauge vector L ∈ TN defined by the condition where h ¼ trðhÞ-is given by the MiSaTaQuWa equation [34,35], is the electric part of the Weyl tensor and h tail is a tail integral of the retarded Green's functions of h. The above is an equation for a four-vector Z called the deviation vector; the lhs is the acceleration associated therewith, and the rhs is a geodesic deviation term plus the GSF. This deviation vector is defined on C correction needed to move off C ∘ and onto the worldline representing the center of mass of the perturbed spacetime, defined as in the Hamiltonian analysis of Regge and Teitelboim [41]. Let us make a few comments on these results, specifically concerning (a) and (c). On the one hand, it is quite remarkable that geodesic motion can be recovered as a consequence 10 of this analysis-i.e., without having to posit it as an assumption-just from smoothness properties (existence of appropriate limits) of our family of metrics fg ðλÞ g, and on the other, this analysis offers sensible meaning to the usual "delta function cartoon" (ubiquitous in essentially all self-force analyses) of the matter stressenergy-momentum tensor describing the object in the background spacetime. The point is that the description of the object is completely arbitrary inside the region that is not covered by the smoothness conditions of the family fg ðλÞ g, i.e., for r ≤ Cλ when λ > 0. (Indeed, this region can be filled in even with exotic matter, e.g., failing to satisfy the dominant energy condition, or a naked singularity, as long as a well-posed initial value formulation exists.) Regardless of what this description is, the smoothness conditions essentially ensure that its "reduction" to ℳ ∘ (or, more precisely, the transport of any effect thereof with respect to the family fg ðλÞ g) simply becomes that of a point particle sourcing the field equation at OðλÞ. In this way, the background "point-particle cartoon" is justified as the simplest possible idealization of a small object.
What we are going to do, essentially, is to accept consequences (a)-(c) [in fact, we will not even explicitly need (b)], the proofs of which do not rely upon any further limiting conditions such as a restriction of the perturbative gauge, and to obtain, using our perturbed momentum conservation law, a more general version of the EoM, i.e., consequence (d). For the latter, Gralla and Wald [33] instead rely on the typical but laborious Hadamard expansion techniques of DeWitt and Brehme [44], wherein the "mass dipole moment" of the object is set to zero. It is possible [41] to have such a notion in a well-defined Hamiltonian sense by virtue of (b). While mathematically rigorous and conducive to obtaining the correct known form of the MiSaTaQuWa equation, their derivation and final result suffer not only from the limitation of having to fix the perturbative gauge but also from the (as we shall see, potentially avoidable) technical complexity of arriving at the final answer-including the evaluation of h tail (or otherwise taking recourse to a regularization procedure).
The link between this approach and our conservation law derivation of the EoM which we are about to carry out is established by the work of Gralla [40], who discovered that Eq. (74) can be equivalently written as Here, the GSF term F½h tail ; U ∘ in the MiSaTaQuWa equation [Eq. (74)] is substituted by an integral expression-an average over the angles-of F. In particular (as, strictly speaking, one cannot define integrals of vectors as such), this is evaluated by using the exponential map based on C ∘ to associate a flat metric, in terms of which the integration is performed over a two-sphere of radius r, S 2 r , with ϵ S 2 denoting the volume form of S 2 .
Observe that, here, the functional dependence of F is on h itself (and not on h tail or any sort of regularized h) and for this reason is referred to as the "bare" GSF. Moreover, this formula is actually valid in a wider class of gauges than just the Lorenz gauge, in particular, it holds in what are referred to as parity-regular gauges [40]. We refrain from entering here into the technical details of exactly how such gauges are defined, except to say that the eponymous "parity condition" that they need to satisfy has its ultimate origin in the Hamiltonian analysis of Regge and Teitleboim [41] and is imposed so as to make certain Hamiltonian definitions-and in particular for Gralla's analysis [40], the Hamiltonian center of mass-well defined. These, however, are not limitations of our quasilocal formalism, where we know how to define energy-momentum notions more generally than any Hamiltonian approach. Thus, in our result, there will be no restriction on the perturbative gauge. This may constitute a great advantage, as the parity-regular gauge class-though an improvement from being limited to the Lorenz gauge in formulating the EoM-still excludes entire classes of perturbative gauges convenient for formulating black hole perturbation theory (e.g., the Regge-Wheeler gauge in Schwarzschild-Droste) and hence for carrying out practical EMRI calculations.
We proceed to apply our quasilocal analysis to the Gralla-Wald family of spacetimes, beginning with a general setup in this family of rigid quasilocal frames.

B. General setup of rigid quasilocal frames in the Gralla-Wald family
Let ðℬ ðλÞ ; u ðλÞ Þ be a quasilocal frame in ðℳ ðλÞ ; g ðλÞ ; ∇ ðλÞ Þ, for any λ > 0, constructed just as described in Sec. II: with unit four-velocity u ðλÞ , unit normal n ðλÞ , induced metric γ ðλÞ , and so on. Using the fact that the tensor transport is linear and commutes with tensor products, we can compute the transport (in the five-dimensional stacked manifold N ¼ ℳ ðλÞ × R ≥ used in our perturbative setup, as in Sec. III A) of any geometrical quantity of interest to the background. For example, 10 See again footnote 2 and the references mentioned therein for more on this topic. ¼ g ab − n a n b Similarly, Now, let us assume that ðℬ ðλÞ ; u ðλÞ Þ is a rigid quasilocal frame, meaning that the congruence defining it has a vanishing symmetrized strain rate tensor in ℳ ðλÞ , θ ðλÞ ðabÞ ¼ 0: ð85Þ ing the transport of the quasilocal observers' four-velocity, n ¼ n ∘ þ λδn þ Oðλ 2 Þ the unit normal, and so on. In other words, ðℬ; uÞ is the background mapping of the perturbed congruence ðℬ ðλÞ ; u ðλÞ Þ, and so will itself constitute a congruence (in the background), i.e., a quasilocal frame defined by a two-parameter family of worldlines with unit four-velocity u in ℳ ∘ . However, although ðℬ ðλÞ ; u ðλÞ Þ is a rigid quasilocal frame in ℳ ðλÞ , ðℬ; uÞ is not in general a rigid quasilocal frame in ℳ ∘ (with respect to the background metric g ∘ ). One can see this easily as follows. Let ϑ ∈ J 0 2 ðℳ ∘ Þ be the strain rate tensor of ðℬ; uÞ, so that it is given by The rhs is a series in λ, owing to the fact that u (and therefore σ, the two-metric on the space H orthogonal to u in ℬ) is transported from a perturbed congruence in ℳ ðλÞ . Upon expansion, we obtain is just the strain rate tensor of the background congruence-i.e., the congruence defined by u ∘ -and is the first-order term in λ. Note that we are abusing our established notation slightly in writing Eq. (87) Since 0 ¼ θ ðλÞ ðabÞ identically in ℳ ðλÞ [as we demand that ðℬ ðλÞ ; u ðλÞ Þ is a rigid quasilocal frame], Eq. (93) must vanish order by order in λ. That implies, in particular, that the zeroth-order congruence (defined by u ∘ ) is a rigid quasilocal frame, and that the symmetrized strain rate tensor of the background-mapped perturbed congruence (defined by u) is given by This tells us that the deviation from rigidity of ðℬ; uÞ in ℳ ∘ occurs only at OðλÞ (and, in particular, is caused by the same perturbed connection coefficient term that is responsible for the GSF). In other words, we can treat ðℬ; uÞ as a rigid quasilocal frame at zeroth order. This zeroth-order congruence actually makes up a different worldtube via the exponential map), and σ ∘ ¼ r 2 S; i.e., it is the metric of S 2 r . This is the most trivial possible rigid quasilocal frame: at any instant of time, a two-sphere worth of quasilocal observers moving with the same four-velocity as the point at its center (parametrizing the given worldline).
At first order, the equation 0 ¼ δθ ðabÞ can be regarded as the constraint on the linear perturbations (δu) in the motion of the quasilocal observers in terms of the metric perturbations guaranteeing that the perturbed congruence is rigid in the perturbed spacetime. (So, presumably, going to nth order in λ would yield equations for every term up to the nth-order piece of the motion of the quasilocal observers, δ n u.) Now, recall the momentum conservation law for rigid quasilocal frames, Eq. (39). This holds for ðℬ ðλÞ ; u ðλÞ Þ in ℳ ðλÞ . Just as we did in the previous section with the general conservation law, we can use φ ðλÞ to turn this into an Let us now further assume that we can ignore the Jacobian determinant (discussed in the previous section) as well as the shift v of the quasilocal observers (relative to constant time surfaces). Then, dividing the above equation by Δt, where t represents the adapted time coordinate on ℬ, and taking the Δt → 0 limit, we get the time rate of change of the momentum, where _ p ðϕÞ ¼ dp ðϕÞ =dt, and we must keep in mind that the derivative is with respect to the adapted time on (the inverse image on the background of) our congruence.

C. Detailed construction of background rigid quasilocal frames
Let G be any timelike worldline in ℳ ∘ . Any background metric g ∘ on ℳ ∘ in a neighborhood of G admits an expression in Fermi normal coordinates [36,50], which we label by fX α g ¼ fT ¼ X 0 ; X I g 3 I¼1 , as a power series in the areal radius. Denoting by A K ðTÞ and W K ðTÞ the proper acceleration and proper rate of rotation of the spatial axes (triad) along G (as functions of the proper time T along G), respectively, this is given by where R 2 ¼ δ IJ X I X J is the square of the radius in these coordinates (not the square of the Ricci scalar) and P KL ¼ δ KL − X K X L =R 2 projects vectors perpendicular to the radial direction X I =R. Here, we have to remember that the Riemann tensor R ∘ IJKL (along with A and W) is understood to be evaluated on G.
For all cases that we will be interested in, we will ignore the possibility of rotation, so we set W I ¼ 0 from now on.
Let us now assume that our background rigid quasilocal frame ðℬ ∘ ; u ∘ Þ is constructed around G; that is to say, into this coordinate system there is embedded a two-parameter family of worldlines representing a topological two-sphere worth of observers, i.e., a fibrated timelike worldtube ℬ ∘ surrounding G. This may be conveniently described, as detailed in Sec. II B, by defining a new set of coordinates fx α g ¼ ft; r; x i g 2 i¼1 given simply by the adapted coordinates ft; x i g 2 i¼1 on ℬ ∘ supplemented with a radial coordinate r. Then, denoting fx i g ¼ fθ; ϕg, we introduce, as done in previous calculations with rigid quasilocal frames in Fermi normal coordinates [57], the following coordinate transformation, Tðt; r; θ; ϕÞ ¼ t þ Oðr 2 ; RÞ; ð102Þ X I ðt; r; θ; ϕÞ ¼ rr I ðθ; ϕÞ þ Oðr 2 ; RÞ; ð103Þ where r I ðθ; ϕÞ ¼ ðsin θ cos ϕ; sin θ sin ϕ; cos θÞ ð104Þ are the standard direction cosines of a radial unit vector in spherical coordinates in R 3 and R here represents the order of the perturbations of the quasilocal frame away from the round two-sphere due to the background curvature effects. In particular, for rigid quasilocal frames, we know that this is in fact simply the order of the Riemann tensor on G, i.e., R ∘ IJKL ¼ OðRÞ. Thus, one may ultimately desire to take OðRÞ effects into account for a full calculation, but for the moment-since, in principle, this R is unrelated to λ and we can assume it to be subdominant thereto-we simply omit them. Thus, we can simply take I ∘ ¼ S 2 r , and we can assume that there is no shift, so thatγ ¼ 1.
Applying the coordinate transformation in Eqs. (102) and (103) to the background metric given by Eqs. (99)-(101) with W ¼ 0, and then using all of the definitions that we have established so far, it is possible to obtain by direct computation all of the quantities appearing in the integrand of the conservation law [Eq. (98)] as series in r. We display the results only up to leading order in r, including the possibility of setting A ¼ 0: Here, j G are, respectively, the electric and magnetic parts of the Weyl tensor evaluated on G. Also, B I i ¼ ∂ i r I and R I i ¼ ϵ S 2 i j B I j are, respectively, the boost and rotation generators of S 2 . See Appendix A for more technical details on this. We remind the reader that E vac and P vac are, respectively, the vacuum energy and pressure, Eqs. (60) and (61), respectively.
The way to proceed is now clear: we expand Eq. (98) as a series in λ, using the zeroth-order parts of the various terms written above. We need only to specify the worldline G in ℳ ∘ about which we are carrying out the Fermi normal coordinate expansion (in r). We will consider two cases: G ¼ C ∘ (the geodesic, such that ℬ is inertial with the point particle in ℳ ∘ ) and G ¼ C (an accelerated worldline such that ℬ ðλÞ is inertial with the object in ℳ ðλÞ , i.e., it is defined by a constant r > Cλ in ℳ ðλÞ ). These will give us equivalent descriptions of the dynamics of the system, from two different "points of view," or (quasilocal) frames of reference.
Before entering into the calculations, we can simplify things further by remarking that the zeroth-order expansions in Eqs. (105)-(110) will always make the twist (ν) term in the conservation law [Eq. (98)] appear at OðrÞ or higher, both in ð _ p ðϕÞ Þ ð0Þ and δ _ p ðϕÞ , regardless of our choice of G. Hence, we can safely ignore it, as we are interested (at least for this work) only in the part of the conservation law which is zeroth order in r. Thus, we simply work with Into this, we furthermore have to insert the multipole expansion of the conformal Killing vector ϕ given by Eq. (41). We correspondingly write such that for any l ∈ N, we have Explicitly, the first two terms are D. Equation of motion inertial with the background point particle Then, A ¼ 0. We will take this to be the case for the rest of this subsection-corresponding, as discussed, to a rigid quasilocal frame the inverse image in the background of which is inertial with the point-particle approximation of the moving object in the background spacetime. This situation is displayed visually in Fig. 9.
Let us first compute the zeroth-order (in λ) part of _ p ðϕÞ . Inserting (105)-(110) into the zeroth-order part of (115) and (116), and making use of the various properties in Appendix A, we find by direct computation: ð _ p ðϕ l¼2 Þ Þ ð0Þ ¼ Oðr 2 Þ: We provide the steps of the calculation in Appendix B. Let us now compute the OðλÞ, l ¼ 1 part of _ p ðϕÞ , i.e., the OðλÞ part of Eq. (115), which as usual we denote by δ _ p ðϕ l¼1 Þ . One can see that this will involve contributions from five OðλÞ terms, respectively containing δN, δE, δα, δP, and δD. For convenience, we will use the notation ð _ p ðϕ l Þ ðQÞ Þ ðnÞ to indicate the term of δ n ð _ p ðϕ l Þ Þ that is linear in Q, for any l, n. Thus, we write All of the computational steps are again in Appendix B. We find If δN does not vary significantly over S 2 r , the Oðr 0 Þ part of the above would be negligible owing to the fact that R S 2 r ϵ S 2 r I ¼ 0. Next, let us consider the δE and δP terms. For this, we find it useful to depict the instantaneous quasilocal frame ðS 2 r ; r 2 S; DÞ embedded in a constant-time three-slice of ℳ ∘ in Fig. 10.
The δE term can be easily determined by realizing that in our current choice of quasilocal frame, the only background matter is the point particle which is always at the center of our present coordinate system; i.e., it is always on C ∘ (on which we are here centering our Fermi normal coordinates). Interpreting the constant m as in the Gralla-Wald approach [33] to be the mass of this point particle, this simply means that so that when this is integrated (as a surface energy density) over S 2 r , we simply recover the mass: We remark that, by definition, it is possible to express the quasilocal energy as E ¼ u a u b τ ab ¼ − 1 κ k with k ¼ σ∶Θ the trace of the two-dimensional boundary extrinsic curvature. Notice that the integral of this over a closed twosurface in the r → ∞ limit is in fact the same as the usual ADM definition of the mass/energy; thus, δE ¼ − 1 κ δk, and so it makes sense to interpret m as the ADM mass of the object. So, now, using Eq. (121), we can find that the δE contribution to δ _ p ðϕ l¼1 Þ is also at most quadratic in r: 10. An instantaneous rigid quasilocal frame ðS 2 r ; r 2 S; DÞ (where S and D, respectively, are the metric and derivative compatible with the unit two-sphere) inertial with the background point particle. This means that the latter is located at the center of our Fermi normal coordinate system.
To compute the δP term, we now employ the useful identity in Eq. (40), which tells us that Using this, into which we insert the δE from Eq. (121), we find that the δP contribution to δ _ p ðϕ l¼1 Þ is at most quadratic in r as well, Note that the above results may in fact be higher order in r than quadratic. We have only explicitly checked that they vanish up to linear order inclusive.
Finally, we are left with the δα and δD contributions to δ _ p ðϕ l¼1 Þ . By direct computation, it is possible to show that their sum is in fact precisely what we have referred to as the extended GSF in our general analysis of the preceding section, i.e., it is the l ¼ 1 part of Eq. (69), In particular, they respectively contribute the usual GSF (from δα) and the gravitational self-pressure force (from δD). Thus, we have found that the total OðλÞ, l ¼ 1 part of the momentum time rate of change is given at leading (zeroth) order in r by nothing more than the generalized GSF. In other words, where we have defined Without loss of generality, let us now pick Φ I ¼ ð0; 0; 1Þ to be the unit vector in the Cartesian X 3 ¼ Z direction, and denote its corresponding conformal Killing vector as ϕ l¼1 ¼ ϕ Z l¼1 . (Alternately, pick the Z axis to be oriented along Φ I .) We know S ij B Z j ¼ ð−1= sin θ; 0Þ; moreover, by the coordinate transformation Inserting these into Eq. (126), we get The first line is precisely in the form of the GSF term from the Gralla formula, Eq. (7) [40], except here in the integrand we have (the Z-component of) our extended GSF F [Eq. (11)]: the usual GSF F (the only self-force term in Gralla's formula) plus our self-pressure term, ℘. The second line contains additional terms involving the extended GSF in the other two (Cartesian) spatial directions. Notice, however, that R r dθ ∧ dϕ cos θ sin ϕ, so if F X and F Y do not vary significantly over S 2 r , their contribution will be subdominant to that of F Z .
Thus, we have shown that our EoM (126) always contains Gralla's angle average of the bare (usual) GSF. However, the form of (126) (expressing the perturbative change in the quasilocal momentum) still cannot be directly compared, as such, with Gralla's EoM (75) (expressing the change in a deviation vector representing the perturbative correction to the motion). In the following subsection, we clarify the correspondence by repeating the calculation using a quasilocal frame inertial with the moving extended object in the perturbed spacetime (rather than with the geodesic in the background, as here). Furthermore, we conjecture that a careful imposition of the parity condition on the perturbative gauge-of which we have made no explicit use so far-would make the contribution from our self-pressure term vanish, but a detailed proof is required and remains to be carried out. is the inverse image of the rigid quasilocal frame ðℬ ðλÞ ; u ðλÞ Þ defined by r ¼ Cλ þ ε ¼ const, ∀ ε > 0, in ℳ ðλÞ . The meaning of the r coordinate in the latter is as given in the Gralla-Wald assumptions (Sec. IVA). This situation is displayed in Fig. 11.
We now proceed to calculate, in the same way as we did for the point-particle-inertial case, the various terms in the expansion of the momentum conservation law, Eqs. (115) and (116). At zeroth order, we obtain The steps of all these computations are again shown in Appendix B.
Let us now compute the OðλÞ, l ¼ 1 part of _ p ðϕÞ .
First, we find that δ _ p ðϕ l¼1 Þ ðδNÞ is the same as in the pointparticle-inertial case, so if δN does not vary significantly over S 2 r , the Oðr 0 Þ part thereof is negligible. Next, let us look at the δE and δP parts. Again, it is useful to consider in this case the visual depiction of the instantaneous quasilocal frame, shown in Fig. 12.  ð136Þ Since Φ I is arbitrary, we thus get the EoM r ; r 2 S; DÞ (where S and D, respectively, are the metric and derivative compatible with the unit two-sphere) inertial with the moving object in the perturbed spacetime. This means that the pointparticle approximation of this object in the background spacetime is not located at the center of our Fermi normal coordinate system. Instead, it is displaced in some direction ρ I , which must be OðλÞ. than in terms of the proper acceleration A of C, we use the generalized deviation equation (as the name suggests, the deviation equation between arbitrary worldlines, not necessarily geodesics), Eq. (37) of Ref. [79]. In our case, this Combining this with Eq. (137), we finally recover the OðλÞ EoM Note that the factor of 3 multiplying the self-force term is in fact present in the EoM in Gralla's Appendix B, that is Eq. (B3) of Ref. [40]. The latter, in this case, expresses the time evolution not of the deviation vector itself but of the change in this deviation vector due to a gauge transformation, possibly including extra terms in case that transformation is out of the parity-regular class. We conjecture that a detailed analysis of the precise correspondence between our deviation vector definition and that of Gralla-Wald (which, while encoding the same intuitive notion of a perturbative correction to the motion, may not be completely identical in general), together with a relation of their gauge transformation properties, would make it possible to relate these EoM exactly.

V. CONCLUSIONS
In this paper, we have used quasilocal conservation laws to develop a novel formulation of self-force effects in general relativity, one that is independent of the choice of the perturbative gauge and applicable to any perturbative scheme designed to describe the correction to the motion of a localized object. In particular, we have shown that the correction to the motion of any finite spatial region, due to any perturbation of any spacetime metric, is dominated when that region is small (i.e., at zeroth order in a series expansion in its areal radius) by an extended gravitational self-force; this is the standard gravitational self-force term known up to now plus a new term, not found in previous analyses and attributable to a gravitational pressure effect with no analog in Newtonian gravity, which we have dubbed the gravitational self-pressure force. Mathematically, we have found that the total change in momentum Δp ðϕÞ ¼ p ðϕÞ final − p ðϕÞ initial between an initial and final time of any (gravitational plus matter) system subject to any metric perturbation h is given, in a direction determined by a conformal Killing vector ϕ (see Sec. II D), by the following flux through the portion of the quasilocal frame (worldtube boundary) ðℬ; u ∘ Þ delimited thereby, where we have restored units, r is the areal radius, and F is the extended self-force functional. In particular, F ¼ F þ ℘ where F is the usual bare self-force [determined by the functional in Eq. (5)], and ℘ is our novel self-pressure force [determined by the functional in Eq. (67)]. The most relevant practical application of the self-force is in the context of modeling EMRIs. Ideally, one would like to compute the correction to the motion at the location of the moving object (SCO). Yet, once a concrete perturbative procedure is established, the latter usually ends up being described by a distribution (Dirac delta function), rendering such a computation ill defined unless additional tactics (typically in the form of regularizations or Green's functions methods) are introduced. However, if one takes a step back from the exact point denoting the location of the "particle" (the distributional support), and instead considers a flux around it, any singularities introduced in such a model are avoided by construction.
We have, moreover, shown that our formulation, when applied in the context of one particular and very common approach to the self-force-namely, that of Gralla and Wald [33]-yields equations of motion of the same form as those known up to now; in particular, they always contain, in the appropriate limit, the angle average self-force term of Gralla [40]. We conjecture that a more rigorous study of these equations of motion and their gauge transformation properties would prove their exact correspondence under appropriate conditions. We would like here to offer a concluding discussion on our results in this paper in Sec. VA, as well as an outlook toward future work in Sec. V B.

A. Discussion of results
From a physical point of view, our approach offers a fresh and conceptually clear perspective on the basic mechanism responsible for the emergence of self-force effects in general relativity. In particular, we have demonstrated that the self-force may be regarded as nothing more than the manifestation of a physical flux of gravitational momentum passing through the boundary enclosing the small moving object. This gravitational momentum, and gravitational stress-energy-momentum in general, cannot be defined locally in general relativity. As we have argued at length in this paper, such notions must instead be defined quasilocally, i.e., as a boundary rather than as a volume densities. This is why the self-force appears mathematically as a boundary integral around the moving object [Eq. (69)], dominant in the limit where the areal radius is small.
The interpretation of the physical meaning of the selfforce as a consequence of conservation principles leads to many interesting implications. As we have seen, the mass of the moving object-e.g., the mass m of the SCO in the EMRI problem-seems to have nothing to do fundamentally with the general existence of a self-force effect. Indeed, according to our analysis, the self-force is in fact generically present as a correction to the motion-and dominant when the moving region is small-whenever one formulate the problem at the "particle location," i.e., the support of the distributional matter stress-energymomentum tensor modeling the moving object (SCO) in the background spacetime. Of course, due to the distributional source, h actually diverges on its support, and so regularization or Green's function methods are typically employed in order to make progress. However, in principle, no such obstacles are encountered (nor are the aforementioned technical solutions needed) if the self-force is evaluated on a boundary around-very close to, but at a finite distance away from-the particle, where no formal singularity is ever encountered; h remains everywhere finite over the integration, and therefore so does the (extended) self-force functional [Eq. (11)] with it directly as its argument.

B. Outlook to future work
A numerical implementation of a concrete self-force computation using the approach developed in this paper would be arguably the most salient next step to take. To our knowledge, no numerical work has been put forth even using Gralla's angle average integral formula [40] (which would further require gauge transformations away from parity-regular gauges).
We stress here that our proposed equation of motion involving the gravitational self-force is entirely formulated and in principle valid in any choice of perturbative gauge. To our knowledge, this is the first such proposal bearing this feature. This may provide a great advantage for numerical work, as black hole metric perturbations h are often most easily computed (by solving the linearized Einstein equation, usually with a delta-function source motivated as in or similarly to the Gralla-Wald approach [33] described in Sec. IVA) in gauges that are not in the parity-regular class restricting Gralla's formula [40]. In other words, we claim that one may solve the linearized Einstein equation [Eq. (71)] for h X in any desired choice of gauge X, insert this h X into our extended GSF functional [Eq. (11)] to obtain F X ½h X ; u ∘ X (for some choice of background quasilocal frame with four-velocity u ∘ ), and then integrate this over a "small radius" topological two-sphere surrounding the particle (so that u ∘ can be approximated by the background geodesic four-velocity of the particle, U ∘ ), to obtain the full extended gravitational self-force (or correction to the motion) directly in that gauge X. It is easy to speculate that this may simplify some numerical issues tremendously vis-à-vis current approaches, where much technical machinery is needed to handle (and to do so in a sufficiently efficient way for future waveform applications) the necessary gauge transformations involving distributional source terms.
Nevertheless, further work is needed to bring the relatively abstract analysis developed in this paper into a form more readily suited for practical numerics. The most apparent technical issue to be tackled involves the fact that h is usually computed (in some kind of harmonics) in angular coordinates centered on the MBH, while the functional F ½h; u ∘ is evaluated in angular coordinates (on a small topological two-sphere) centered on the moving particle, i.e., the SCO. A detailed understanding of the transformation between the two sets of angular coordinates is thus essential to formulate this problem numerically. This issue is discussed a bit further in Gralla's paper [40], but a detailed implementation of such a computation remains to be attempted. The abstraction and generality of our approach may, on the other hand, also provide useful ways to address some other technical issues surrounding the self-force problem. For example, all the calculations in this paper may be carried on to second order (in the formal expansion parameter λ). This is conceptually straightforward given our basic perturbative setup, but of course requires a detailed analysis in its own right. Nonetheless, one may readily see that any higher-order correction to the motion manifestly remains here in the form of a boundary fluxonly now involving nonlinear terms in the integrand. Thus, any sort of singular behavior is avoided at the level of the equations of motion in our approach, up to any order.
As another example, if ever desired (e.g., for astrophysical reasons), linear or any higher-order in r (the areal radius of the SCO boundary) effects on the correction to the motion can also be computed using our approach. Moreover, any matter fluxes (described by the usual matter stress-energy-momentum tensor, T) can also be accommodated thanks to our general (gravity plus matter) conservation laws [Eq. (39)].
Furthermore, while we have applied our ideas in this paper to a specific self-force approach-that of Gralla and Wald [33]-our general formulation (Sec. III) can just as well be used in any other approach to the gravitational selfforce, i.e., any other specification of a perturbative procedure (of a family of perturbed spacetimes fðℳ ðλÞ ; g ðλÞ g) for this problem. In other words, our approach permits any alternative specification of what is meant by a (sufficiently) "localized source" in general relativity, as our conservation expressions always involve fluxes on their boundaries and are not conditioned in any way by the exact details of their interior modeling. Thus, our equation of motion [Eq. (139)] could be used not only for a self-consistent computation (using, e.g., an approach such as that of Refs. [80,81] for solving the field equations in this context) within the Gralla-Wald approach, but also, for example, in the context of the (mathematically equivalent) self-consistent formulation of Pound [77].
Beyond the gravitational self-force, another avenue to explore from here-of interest at the very least for conceptual consistency-is how our approach handles the electromagnetic self-force problem. Although undoubtedly some conceptual parallels may be drawn between the gravitational and electromagnetic self-force problems (see, e.g., Ref. [11]), foundationally they are usually treated as separate problems. Indeed, shortly after the paper of Gralla and Wald [33] detailing the self-force approach used in this work, Gralla et al. [82] put forth a similar analysis, with an analogous approach and level of rigor, of the electromagnetic self-force. It would be of great interest to apply our quasilocal conservation laws in this setting, as they can be used to account not just for gravitational but also (and in a consistent way) matter fluxes as well. It may thus prove insightful to study how the transfer of energymomentum is actually accounted for (between the gravitational and the matter sector), as in our approach we are not restricted to fixing a nondynamical metric in the spacetime. In other words, the conservation laws account completely for fluxes due to a dynamical geometry as well as matter.