Local Continuity of Angular Momentum and Noether Charge for Matter in General Relativity

Conservation laws have many applications in numerical relativity. However, it is not straightforward to define local conservation laws for general dynamic spacetimes due the lack of coordinate translation symmetries. In flat space, the rate of change of energy-momentum within a finite spacelike volume is equivalent to the flux integrated over the surface of this volume; for general spacetimes it is necessary to include a volume integral of a source term arising from spacetime curvature. In this work a study of continuity of matter in general relativity is extended to include angular momentum of matter and Noether currents associated with gauge symmetries. Expressions for the Noether charge and flux of complex scalar fields and complex Proca fields are found using this formalism. Expressions for the angular momentum density, flux and source are also derived which are then applied to a numerical relativity collision of boson stars in 3D with non-zero impact parameter as an illustration of the methods.

Throughout this work the metric has sign {−, +, +, +} and physical quantities will be expressed as a dimensionless ratio of the Planck length L pl , time T pl and mass M pl unless stated otherwise; for example Newtons equation of gravity would be written as where F pl = M pl L pl T −2 pl is the Planck force. Consequently c, G and take the numerical value of 1. Additionally, tensor fields will be denoted using bold font for index free notation and normal font for the components. The dot product between two vector fields will be written interchangeably as A · B ↔ A µ B µ for readability. Additionally, ∇ µ denotes the covariant derivative and ∂ µ is the partial derivative, both with respect to coordinate x µ . Finally, unless stated otherwise, Greek indices such as {α, β, ..., µ, ν, ...} label four dimensional tensor components whereas late Latin indices such as {i, j, k, ...} label three dimensional tensor components and early Latin indices such as {a, b, ...} label two dimensional ones.

I. INTRODUCTION
Conservation laws play an important role in many areas of physics. For a general Lagrangian density L, dependant on fields φ i and derivatives ∂ k φ i for i ∈ {1, 2, ..., m} and k ∈ {1, 2, ..., n}, if a field transformation φ i → φ i + δφ i leaves the Lagrangian constant the Euler-Lagrange equations imply there is a conserved current J , * rc634@cam.ac.uk with zero divergence, given by In curved space a conserved current J satisfies ∇ µ J µ = 0. A charge Q within 3-volume V and a flux F though ∂V , the boundary of V , can be associated with J as described later in Eqs. (20) and (21). If J is conserved then Q is a conserved charge satisfying This says the rate of change of a charge in a volume V is equal to the flux across the boundary ∂V of V . In the case that J has a non-zero divergence, ∇ µ J µ = 0, Eq. (3) generalises to the continuity equation where S is defined in Eq. (22); S is the source of J in V which can be understood as the destruction or creation of charge Q. Eq. (4) is a simplified version of Eq. (16), later referred to as the QFS system. Evaluation of the continuity equations above, and their corresponding charges Q, have many uses in the study of fundamental fields in Numerical Relativity. One such use is the measurement of the Noether current J of a complex scalar/vector field which arises from a U (1) gauge symmetry of the matter fields ψ j of the form ψ j → ψ j e ia ∼ ψ j + iaψ j for some small constant a and j ∈ {1, 2, ..., n}. The total charge Q, also called Noether charge in this case, is useful to track during numerical simulations as it gives insight into the numerical quality of a simulation. A violation of Noether charge conservation can arise from insufficient resolution in some region of the simulation or due to boundary conditions in a finite volume simulation. In the case of Sommerfeld (outgoing wave) boundary conditions [1] we might expect charge to be transported out of a finite computational domain and the total charge Q in the simulation should decrease.
Monitoring only Q within some volume V , it is impossible to tell whether Noether charge violation is due to a flux F through the surface ∂V or undesirable numerical inaccuracies such as dissipation. It is more useful to check whether the continuity Eq. (3) (or equivalently Eq. (4) if there were a non-zero source term) is obeyed for a finite domain V ; if this fails there is likely a problem as the continuity equations should be exactly observed for general spacetimes.
Another use of the continuity equations is to measure the amount of energy-momentum belonging to matter fields within a volume V . This has many possible applications such as calculating the total energy or momentum of compact objects such as boson stars and neutron stars. The energy-momentum of matter obeys a conservation law in General Relativity as given by Penrose [2] where the considered spacetime is assumed to admit a Killing vector. In the case a Killing vector exists then a conserved current J associated with the energy-momentum tensor T can be identified. The current is J µ = T µ ν ξ ν for some Killing vector ξ and satisfies ∇ µ J µ = 0. If ξ is a Killing vector then Eq. (3) is the correct continuity equation and the charge Q is conserved. In General Relativity the existence of Killing vectors is rare, reserved for spacetimes with special symmetries. Generic dynamic spacetimes with no symmetries, such as inspirals and grazing collisions of compact objects, have no easily identifiable Killing vector fields. If there is no Killing vector the divergence of J becomes ∇ µ J µ = T µν ∇ µ ξ ν and the source term S is non-zero. Now Eq. (4) is the correct continuity equation and the charge Q is no longer conserved. In section V we will show how the choice of ξ affects the type of current J , and therefore charge Q, obtained. While measures of energy or momentum are interesting in their own right, the measure of Eq. (4) within some volume V can be a good measure of numerical quality of a simulation in a similar fashion to the measure of Noether charge mentioned already.
When dealing with black hole spacetimes resolution requirements typically become very strict towards the singularity and lead to a local violation of Eqs. (3) and (4). This might not doom a simulation as for most physical applications in GR singularities are contained by an event horizon and are therefore causally disconnected from the rest of the simulation; a resolution problem in the vicinity of a singularity therefore may not propagate to the exterior. It could be helpful instead to consider a volumē V equal to V but removing a set of finite volumesṼ i which surround any singularities. Testing Eqs. (3) and (4) in volumeV would then give a measure of the simulation resolution untainted by the resolution issues at a singularity.
Currently in Numerical Relativity it is common to measure energy-momentum in a localised region with with Eq. (20) for the charge Q where the charge density Q is given in section V. Examples of this can be seen in [3], [4], [5] and [6]. While this is a good measure it neglects any radiation and the transfer of energy-momentum between matter and spacetime curvature; if the spacetime does not contain the corresponding Killing vector the charge cannot be treated as a conserved quantity. Instead the combination of variables Q − S, where S is defined in Eq. (22) and S in section V, should be treated as a conserved quantity. Other popular methods to obtain the energy-momentum of a system include integrating asymptotic quantities such as the ADM mass and momentum, however these can not be used locally as they are defined in the limit of large radii only.
Recent work by Clough [7] evaluates Q, F and S for energy and linear momentum with the assumption that the approximate Killing vector ξ is a coordinate basis vector satisfying ∂ i ξ j = 0. Successful numerical tests of Eq. (4) are given for fixed and dynamic background simulations.
This paper builds on the work of [7] and generalises the system to measure angular momentum conservation and the conservation of Noether charges of complex scalar fields and spin-1 complex Proca fields. The assumption that the approximate Killing vector ξ is a basis vector satisfying ∂ i ξ j = 0 is dropped and leads to a more general source term S. The QFS system for angular momentum is also tested using fully non-linear numerical Relativity simulations of a spacetime consisting of two boson stars colliding in a grazing fashion. This paper is organised as follows. In section II the QFS system for a general non-conserved current is derived and section III explicitly expands the results for use with a spherical extraction surface. Even though no other extraction surfaces are considered, the results of III are easily adaptable to other shapes. Section IV is a standalone derivation of the well known Noether charge density from the QFS perspective and goes on to find the flux variable; results for complex scalar fields and complex Proca fields are given. The application of the QFS system to energy momentum currents, angular momentum and energy are given in section V. A fully non-linear test of the QFS system for angular momentum, using GRChombo [8,9] to perform Numerical Relativity simulations, is presented in section VI along with a convergence analysis.

II. DERIVATION OF THE QFS SYSTEM
For a spacetime (M, g) we start by defining a vector field J and subjecting it to the following continuity equation, where S is a source term and describes the nonconservation of J . In the case S = 0 the current is conserved. We are interested in the charge density Q and source density S associated with J in a spatial 3volume V ∈ Σ. Here Σ is the usual 3-dimensional spacelike manifold Σ consisting of the set of all points with constant time coordinate t, equipped with metric γ. We are also interested in the flux density F through ∂V , the boundary of V with metric σ. Σ is spanned by spatial coordinates x i related to the full spacetime coordinates x µ by x µ = {t, x i }. The normal to Σ is the unit co-vector n defined as, where (α, β i ) are the usual lapse and shift from the ADM 3+1 spacetime decomposition [10]. The reader is directed to [11] for a comprehensive introduction to the 3 + 1 decomposition. In Eq. (7), t µ = (1, 0, 0, 0) is the future directed vector and is distinct from n µ . Time vector t is useful as its integral curves form lines of constant spatial coordinates. With this knowledge we can define the 4volume M , the spatial 3-volume V evolved along integral curves of t between times t 0 ≤ t ≤ t 0 +δt in the limit δt → 0. Finally we define the 3-dimensional volume H, with metric h. H is the evolution of ∂V along integral curves of t between times t 0 ≤ t ≤ t 0 + δt and is the 3-volume the flux crosses; clearly our definition of H will affect our definition of flux density. There is no reason to choose the timelike vector t, rather than n, to evolve V and ∂V in time and both will result in a different definition of flux density. However, it is shown in appendix B that these two choices result in the same total integrated flux. A diagram summarising the relevant geometry can be found in Fig. 1.
With the relevant geometry discussed we can derive the QFS system, Eq. (4), by integrating Eq. (5) over M ; Let us start by using Gauss' theorem for curved space [12] on the left hand side, with (3) g being the volume element of a generic 3surface andŝ being the corresponding unit normal. Note thatŝ is outward directed when spacelike and inward directed when timelike. The integral of ∇·J over M is now transformed to a surface integral of J over the compound 3-volume ∂M . This surface is split into an integral over H and two integrals over V at times t 0 and t 0 + δt. The integrals over V give, where we made use of the fact that the coordinate volume V is constant for all times due to it evolving in time with t µ . Here √ γ is the volume element on the spacelike manifold Σ. Now let us evaluate the integral over H, with metric h of signature {−, +, +}, where N is the unit normal to H as shown in Fig. 1. Given that n is not tangent to H (but the time vector t is) we have N · n = 0 and N · t = 0. This means we must normalise N with metric g and not γ as N is not tangent to Σ and g(N , N ) = γ(N , N ). In other words, N is a 4-vector in the case we time evolve V with time vector t but N is a 3-vector if we time evolve with n. Finally the right hand side source integral from (8) becomes, Combining Eqs. (11), (13) and (15) transforms Eq. (8) into, where the density, flux density and source density (Q, F, S) of angular momentum are defined as, The integrated version of these quantities can be written as For later sections it is useful to split the normal vector N into its spacelike and timelike parts, where ⊥ N µ =⊥ ν µ N ν is the projected part of N onto Σ with n · ⊥ N = 0 and ⊥ is the projection operator onto Σ, Using Eqs. (23) and (24) the flux term becomes, The term on the left arises from flux through the surface ∂V and the term on the right is a consequence of the coordinate volume V moving with respect to a normal observer with worldline traced by n. Writing the flux term as above makes it obvious how the definition of flux depends on the 3-volume H which determines √ −h, √ σ and N . Equivalently it can be seen that the density (17) and source (19) terms do not depend on the extraction surface.

III. APPLICATION TO SPHERICAL EXTRACTION
The numerical application of the QFS system in section VI chooses a spherical volume V to extract the angular momentum flux. Using standard Cartesian and spherical polar coordinates, x i cart = {x, y, z} and x i polar = {r, θ, φ} respectively, we can define H as the coordinate volume r = r 0 , t 0 ≤ t ≤ t + δt. Thus, the normal N to H is proportional to ∇(r − r 0 ) = ∇( x 2 + y 2 + z 2 − r 0 ). Explicitly calculating the components N µ , and normalising to unity, with spherical polar spacelike coordinates gives, Note that if we had chosen H, the future evolution of ∂V , to be evolved along the unit vector n rather than time vector t then we would have obtained a different definition of N perpendicular to n rather than t. The consequences of the alternate choice of H are explored in appendix B. The density and source terms Q and S do not depend on the integration domain V , but the flux term F does. The calculation of the Flux term requires the evaluation of the volume element √ −h of H. Due to the choice that H is the surface of constant radial coordinate, finding the metric of this surface is straightforward. Here we define spherical polar coordinates x µ = {t, r, θ, φ} on M and X m = {t, θ, φ} spanning H. Projecting the 4-metric g onto H we can write (4) where (4) h belongs to M. The line element of a curve residing in H can be equivalently evaluated in M or H; the pullback of (4) h from M| r=r0 to H gives the 3-metric h belonging to H, A similar argument can be made for ∂V , the set of all points satisfying r = r 0 and t = t 0 , with metric σ. The metric components and volume element are Using Cramer's rule for the inverse of a matrix with Eqs. (31) and (32) we get and reading from Eq. (29) gives Similarly to Eq. (30), the pushforward of h on H to (4) h on M| r=r0 gives which shows that h tt = (4) h tt . Combining Eqs. (34) and (36) it can be shown that Using this with Eq. (28) we can expand Eq. (18) for the flux term F, but for practical purposes it is helpful to decompose this in terms of 3+1 variables as in Eq. (26), It is straightforward to re-derive the results of this section for other extraction volume shapes, such as cylinders or cubes/rectangles, with a redefinition of spacelike volume V giving rise to different N , H and ∂V . It is wise to pick suitable coordinates adapted to the symmetry of the problem.

IV. NOETHER CURRENTS
In this section we apply the previous results of the QFS system Eq. (16) to the continuity of Noether charge for both a complex scalar and the Proca field. The charge Q represents the number density of particles. Since the total particle number minus antiparticle number is always conserved the conservation law is exact and the source term S vanishes.
Globally the total Noether charge should be conserved, however in numerical simulation this might not always be the case. Two common ways for non-conservation to occur are for the matter to interact with the simulation boundary conditions (often unproblematic) or some region of the simulation being insufficiently resolved (often problematic). Without knowledge of the Noether flux it is difficult to know what a change in Noether charge should be attributed to. If a large extraction volume containing the relevant physics shows a violation of Eq. (16) then the change in total Noether charge being due to boundary conditions can be ruled out and resolution is likely the culprit.
When considering black hole spacetimes it is common for Noether charge to be dissipated as matter approaches the singularity; this is due to resolution requirements typically becoming very high in this region. The violation of Eq. (16) inside a black hole horizon might not cause any resolution problems for the black hole exterior however due to causal disconnection. If the extraction volume is modified to exclude finite regions containing any black hole singularities then Eq. (16) could be used to monitor the conservation of Noether charge away from troublesome singularities. This would be a good way of checking the resolution of a matter field in situations such as boson/Proca stars colliding with black holes or scalar/vector accretion onto a black hole.

Complex Scalar Fields
In this section we consider the conserved Noether current associated with a complex scalar field ϕ with Lagrangian where V is some real potential function. There is a U(1) symmetry where a complex rotation of the scalar field ϕ → ϕe ia , for constant a, leaves the action unchanged. The associated Noether current J can be found in [13], and satisfies ∇ · J = 0. The conservation is exact here which tells us the source term vanishes. In this case the Noether charge density (17) and flux density (26) are, where Π = −n · ∇ϕ is the conjugate momentum of the scalar field. Using Eq. (43) for a spherical extraction surface, and spherical polar spacelike coordinates {r, θ, φ}, this explicitly becomes

Complex Vector Fields
The Complex vector field A, also called a Proca field, has Lagrangian Again V is some real potential function. The action is invariant under a similar U (1) complex rotation of the vector field A µ → A µ e ia for constant a. Following [14] this leads to the following Noether current J , which again satisfies ∇ · J = 0 and the source term vanishes. Defining a 3 + 1 decomposition compatible with [15] gives, where φ, E, a and B all belong to Σ. Additionally D is the covariant 3-derivative of Σ. Note that F has no time-time component as n µ n ν F µν = 0 from the antisymmetry of F . Using these, the Noether charge (17) becomes, = i n µāν n µ E ν − n µ a ν n µĒν , Using Eq. (26), the Noether flux is, where the Christoffel symbols from D µ in B µν cancel out. Putting this into the expression for the Proca Noether flux we get, and equation (43) gives the flux term for a spherical extraction surface, where spherical polar spacelike coordinates {r, θ, φ} used.

V. ENERGY-MOMENTUM CURRENTS
To find the current associated with energy-momentum we consider a vector field J defined with respect to a second vector field ξ and the stress tensor T by Calculating the divergence of this vector leads to the following continuity equation, where (72) shows the divergence vanishes if ξ is a Killing vector of the spacetime; a vanishing divergence corresponds to a conserved current with a zero source term. For more general spacetimes where ξ is not Killing, the right hand side of (72) leads to a non-zero source term accounting for the transfer of energy-momentum between matter and spacetime curvature [7]. The choice of ξ dictates the type of energy-momentum current retrieved; for instance ξ = ∂ t will correspond to an energy current J µ = T µ ν (∂ t ) ν = T µ t and the spacelike choice ξ = ∂ i gives a momentum current J µ = T µ ν (∂ i ) ν = T µ i corresponding to the coordinate x i . For an account of energy and linear momentum continuity see [7].

Angular Momentum
The numerical test of the QFS system (16) in section VI measures the conservation of angular momentum. To do this we choose ξ = ∂ φ , used in Eq. (70), which is the coordinate basis vector of some azimuthal coordinate φ. The angular momentum current is Any spacetime with azimuthal symmetry (e.g. the Kerr spacetime) will have a vanishing source term as ∂ φ is a Killing vector. This includes numerical simulations of matter in a fixed background. The example simulation in section VI is the fully nonlinear grazing collision of two boson stars and ξ is not a Killing vector for finite distances from the collision centre. In this case the source term is non-zero. Using the standard 3+1 decomposition of the stress tensor [11], [16] and explicitly expanding the density term from Eq. (17) gives, where x and y are Cartesian coordinates related to spherical polar coordinates in the usual way. Combining Eqs. (26), (73) and (77) we can get the angular momentum flux through a spherical extraction surface, Using Eq. (43), for a spherical extraction surface, the flux term becomes, in spherical polar coordinates. The explicit expansion of the source term S is left for the appendix A, but the result is given here, As noted in appendix A, when choosing a coordinate system to evaluate S, if ξ is a coordinate basis vector then the ∂ i ξ j terms vanish. It would be simple to re-derive these results for linear momentum, by using ξ = ∂ i for momentum in the x i direction for example, where x i is some Cartesian spatial coordinate. Results for linear momentum can be found in [7].

Energy
A local conservation system can also be applied to energy with the choice of an approximate Killing vector ξ, ξ µ = (∂ t ) µ = t µ = (1, 0, 0, 0), and energy current Using the standard 3+1 decomposition of the stressenergy tensor from [11] or [16] the energy density Q is, from Eq. (17). Similarly, combining Eqs. (18) and (23), the energy flux is, and Eqs. (28) and (39) can be used for a spherical extraction surface, The source term is omitted here as the expression derived in appendix A assumes that ξ is spacelike. For energy continuity a timelike approximate Killing vector ξ is used and leads to a different expression for the source term that can be found in [7] along with the above density Q and flux term F.

VI. NUMERICAL APPLICATION
To numerically test the QFS system, given in Eq. (16), for angular momentum an example spacetime consisting of colliding boson stars is simulated in 3D using GRChombo [8,9]. GRChombo is a modern, open source, Numerical Relativity code with fully Adaptive Mesh Refinement (AMR) using the Berger-Rigoutsos blockstructured adaptive mesh algorithm [17]. The CCZ4 constraint damping formulation [17,18] is used with the moving puncture gauge [19,20]. Time integration is done with 4th order Runge-Kutta method of lines.

Numerical Setup of Simulations
The boson stars are composed of a complex scalar field ϕ, minimally coupled to gravity with the Lagrangian given in Eq. (44). Boson stars are stable self-gravitating spherically symmetric solutions of the Einstein-Klein-Gordon system in curved space; for a detailed review see [13]. In this work, the Klein-Gordon potential is chosen to be V = m 2 ϕφ, where m is the mass of a bosonic particle, leading to so called mini boson stars. The Kaup limit for the maximum mass of a mini boson star can be found numerically as approximately where the physical constants are included for completeness, but have numerical value 1 in Planck units. Notably, the maximum mass of a mini boson star scales inversely with the boson particle mass m.
The Lagrangian in Eq. (44) with potential V = m 2 ϕφ is unchanged up to an overall constant under a rescaling of the boson mass like m → bm, for some dimensionless constant b, while simultaneously rescaling x µ → b −1 x µ for coordinates with dimension length/time. Consequently a mini boson star solution, categorised by the central scalar field amplitude ϕ c , represents a one parameter family of solutions with ADM mass and radius inversely proportional to m. To keep the choice of m arbitrary the coordinates used in the simulation are mx µ , which are exactly Planck units in the case m = 1 (i.e. the Planck mass).
To measure the charge associated with angular momentum, the following angular momentum measures are considered, where Q, F and S are defined in Eqs. (74), (83) and (84) respectively.Q is the angular momentum modified by δQ S ; this is equivalent to absorbing the source term into Q.Q is the initial angular momentum modified by the time integrated total flux. Equation (16) impliesQ =Q exactly, and we define the relative numerical error measure e 1 by which converges to zero in the continuum limit. We can alternatively define a different relative error where the source term is not absorbed into Q. Again, Eq. (16) implies that Q =Q, or e 2 = 0, in the continuum limit.
The initial data of the numerical simulations consists of two boson stars, each with mass M = 0.395(0) m −1 , boosted towards each other in a grazing configuration. The data for two single boosted stars are superposed as in Ref. [21] to minimise errors in the Hamiltonian and momentum constraints and spurious oscillations in the scalar field amplitudes of the stars. The physical domain  1, 0, 0) along the x axis. The stars travel towards each other and undergo a grazing collision to form a short lived dense object at time t ∼ 375 m −1 . Afterwards, much of the scalar field (and angular momentum) leaves the extraction radii as it is ejected to spatial infinity. Figures. 2, 3 and 4 show the angular momentum within radii r = {20, 40, 60} m −1 . The Newtonian angular momentum for this configuration is which is in close agreement with Q andQ in Figs. 2, 3 and 4 while the matter is contained by the extraction radii. Given that we are dealing with a fully non-linear spacetime in general relativity there is no reason why the naive Newtonian angular momentum should agree so well with the numerically integrated values Q orQ; this could be due to the mass of the stars being M = 0.395(0) m −1 , well below the Kaup limit M Kaup ∼ 0.633 m −1 and the mild boost velocities v = 0.1. In the case that the star masses/densities and velocities tend to zero we expect general relativity to approach the Newtonian limit; conversely for large masses/densities and boost velocities the Newtonian estimate likely becomes less accurate. Finally we note in Figs. 2, 3 and 4 that the sourcecorrected density variableQ is less prone to oscillations than Q and is closer to being constant at early times when no angular momentum flux is radiated.Q has another advantage over Q; at extraction radii sufficiently far from any matterQ will remain constant due to the flux F vanishing. Q will only remain constant if the source term integral δQ S also remains constant which does not happen in general dynamic spacetimes, even for large extraction radii.

Convergence Analysis
Three numerical simulations are used to test the convergence of the angular momentum measures as the continuum limit is approached. They have N ∈ {320, 384, 448} gridpoints on the coarsest level, named level 0 with grid spacing ∆x 0 = L/N . Each finer level, named level n, has grid spacing ∆x n = 2 −n ∆x 0 . Any gridpoints that fall inside radius r = 200 m −1 are forced to be resolved by at least AMR level 1. Similarly, any points within radius r < 60 m −1 are resolved by at least AMR level 3; this modification quadruples the default resolution for r < 60 m −1 compared to level 1. These two radii have a 20% extra buffer zone to ensure that AMR boundaries are outside and away from the desired radii. On top of this the AMR is triggered to regrid when a tagging criterion is exceeded; a description of the algorithm can be found in section 2.2.2 of [22]. The tagging criteria used in this paper involve gradients of the scalar field and spatial metric determinant; this loosely means as a region of spacetime becomes more curved, or matter becomes denser, the region is resolved with higher resolution. Figs. 5 and 6 show the relative errors e 1 and e 2 for the convergence sequence; it can be seen that e 1 , the relative error ofQ, is less prone to oscillations than e 2 , the relative error of Q. The choice of enforcing AMR regridding to level 3 within r < 60 m −1 is very problem specific and the grid structure has been chosen carefully for the particular physical scenario to give higher resolution around the late time scalar field configuration at the origin; this enables accurate simulation of the extended object after merger. Simulations prior to this modification showed approximately five times higher relative error e 1 and much worse Noether charge conservation. The highest resolution simulation, with N = 448, shows that the relative error e 1 is 3% after 8000 m −1 time units.
We now obtain the order of convergence ω of e 1 . It is convenient to express e 1 as three functions {f 1 , f 2 , f 3 } corresponding to the three different resolution simulations with N = {320, 384, 448} and f ∞ to denote the continuum limit solution. A traditional convergence analysis, as in [23], assumes that the numerical error of a function (i.e. difference from f ∞ ) is dominated by a term proportional to ∆x ω i for an order of convergence ω; thus we can write for some constant coefficient E for all resolutions i. Equation (104) with i = {1, 2, 3} can be used to eliminate both E and f ∞ giving the well known result for ideal convergence. Figure 7 shows for three orders of convergence ω = {2, 3, 4}; the two expressions should be equal for an ideal order of convergence ω. It can be seen by eye that ω = 3 is the best estimate. The black curve shows the difference between f3 and f2; the relative error e1 of the two highest resolution simulations in Section VI 2. The three coloured curves show the difference between the two lowest resolution simulations f2 and f1, but modified by (∆x ω 3 − ∆x ω 2 )/(∆x ω 2 − ∆x ω 1 ) in accordance with Eq. (105), for three idealised orders of convergence ω = {2, 3, 4}. The black curve is in best agreement with the green curve giving an estimate of ω = 3 for the order of convergence.
To quantify the order of convergence, rather than guessing, we define the deviation factor D as which averages the violation of Eq. (105) between times t 0 ≤ t ≤ t 1 . Figure 8 plots D versus ω with a red curve and the order of convergence can be estimated by minimising D(ω) with respect to ω. As can be seen in Fig. 8, the traditional order of convergence is approximately 3.2.
Given that e 1 vanishes in the continuum limit we can set f ∞ = 0 to find the order of convergence to zero. Using Eq. (104) with the two highest resolutions i = {2, 3}, and setting f ∞ = 0, E can be eliminated to give Similarly to before, we can define a deviation factorD, which time averages the violation of Eq. (107). The black curve in Fig. 8 plots Eq. (108) versus ω and the order of convergence to zero can be estimated by minimisingD(ω) with respect to ω. As can be seen in Fig. 8, the order of convergence convergence to zero is approximately 1.9. The red curve gives the deviation from ideal traditional convergence D, defined by Eq. (106), as a function of ω. The black curve shows the deviation from ideal convergence to zeroD, using the definition given in Eq. (108), as a function of ω. For both curves the estimated order of convergence is found by minimisation with respect to ω; this gives ω = 3.2 for traditional convergence and ω = 1.9 for convergence to zero.

VII. CONCLUSION
A derivation of the QFS system (16) for continuity equations, valid locally for general spacetimes, is derived and applied to spherical integration surfaces. Although spherical extraction surfaces are used, the methods of section III can be applied to general extraction surfaces with minor adjustments. The QFS system is used to calculate the well known Noether charge densities for complex scalar and complex vector (Proca) fields along with novel expressions for the flux variable F in section IV. Next the QFS system for energy momentum currents associated with matter are found and the main result of this paper is the explicit derivation of the angular momentum QFS variables Q, F and S. The three variables can be used to measure the angular momentum of matter within a region, the flux of angular momentum of matter through the boundary of that region and the transfer of angular momentum between matter and curvature; they can also be used with Eq. (16) to determine the numerical quality of a simulation as the QFS system is exactly satisfied in the continuum limit. In section VI the combination of variables Q and S is shown to be a superior measure of angular momentum than integrals of only the charge density Q in two ways; firstly its measurement is less prone to oscillations and secondly it is conserved in the large radius limit.
The QFS system for angular momentum is numerically tested on a dynamic non-linear spacetime consisting of two colliding boson stars; the collision has a small impact parameter giving rise to a non-zero total angular momentum. The stars promptly collide and form a highly perturbed, localised scalar field configuration partially retaining angular momentum. The total angular momentum of the spacetime is measured using the QFS variables (Eqs. (77), (83) and (84)) and is shown to agree well with the Newtonian approximation. This is a good check on the normalisation of the QFS variables as they should return the Newtonian calculation in the low energy limit; even though we simulate a fully non-linear spacetime the density and boost velocity of the stars are mild. The final numerical result is the convergence test of the QFS system which measures the relative error described in VI. The relative error converges to zero with order ω ≈ 1.9 in the continuum limit and the highest resolution simulation gives a fractional error of approximately 3% in the total angular momentum after 8000 time units.
The QFS system is straightforward to implement and it is hoped these results will be useful to the Numerical Relativity community for better measurement of local energy-momentum of matter and Noether charge aswell as powerful check on simulation resolution. This gives us our final form for the angular momentum source density, If we pick a coordinate basis vector as our approximate Killing vector, for example with components ξ µ = (∂ φ ) µ = (0, 0, 0, 1) µ in polar coordinates, then the ∂ µ ξ ν terms will vanish. However if we wish to work in Cartesian coordinates, which is very common for numerical codes, then the vector componentsξ µ become, wherex µ are Cartesian coordinates and x µ are spherical polar coordinates.

Appendix B: Generality of Result
Here we demonstrate that the choice of 4-volume M integrated in Eq. (8) does not change the resulting QFS system (4). We start by defining the extraction 3-volume V 1 ∈ Σ at time t = t 0 . The boundary of V 1 is the 2volume ∂V 1 with metric σ. As in section II, Σ is the 3-manifold defined by the set of all points with constant time coordinate t, equipped with metric γ and unit normal n like Eq. (6). We now choose to define a 4-volumẽ M , different to M , as the evolution of V 1 along integral curves of n between times t 0 ≤ t ≤ t 0 + δt in the limit δt → 0. The boundary ofM , ∂M , is composed of three coordinate 3-volumes, V 1 ∈ Σ t , V 2 ∈ Σ t+δt andH. Here V 2 = V 1 + δV and is the future of V 1 , at time t = t 0 + δt, found by following integral curves of n.H is the 3-volume defined by the time evolution of 2-volume ∂V 1 with n. A diagram showing the differences between the choices of time evolution vectors n and t is given in Fig. 9.
We start again by using Gauss' theorem like in Eq. (9) which results in three surface integrals over V 1 ∈ Σ t , V 2 ∈ Σ t+δt , andH, whereÑ is the unit normal toH. Starting with the integrals over V 1 and V 2 = V 1 + δV we get, FIG. 9. Comparison of two possible geometries for derivation of QFS system in section II on manifold M. Σt is the spatial hypersurface at time t and Σ t+δt is the the spatial hypersurface at a later time t + δt. V1 is the coordinate volume, with surface ∂V1 (not labelled), that we wish to use as an extraction volume on Σt. The red cylinder, defined by ∂V1 evolved along integral curves of t = ∂t, is the same as in Fig. 1. Evolving ∂V1 forward in time with n, as demonstrated with the blue cylinder, gives a different coordinate volume ∂V2 (not labelled) on Σ t+δt . Similarly to the red cylinder, the blue cylinder has sides labelled byH and an interiorM . The difference in the volumes V1 and V2 on Σ t+δt is denoted by δV.
with the new integral over δV appearing because V 1 = V 2 as V 1 is evolved along integral curves of n rather than time basis vector t = ∂ t ; this is demonstrated in Fig. 9.
In the limit that δt → 0 it can be seen that, where s i are the components of the unit normal to the coordinate surface ∂V 1 in R 3 rather than Σ t . The overall negative sign in Eq. (B4) comes from the defined direction of the shift vector β as seen Fig. 9. Addressing the integral overH gives, which is in the same form as Eq. (16) with definitions, where we used n · J = Q for the flux term. The density term Q and source term S are agnostic to our choice of extraction surface and its time evolution so have turned out the same as Eqs. (17) and (19). At first glance the flux term F seems different to Eq. (26) but evaluating this term in a coordinate basis will show otherwise.
Choosing a spherical extraction surface as in Section III and using spherical polar coordinates, x µ = {t, r, θ, φ}, V 1 becomes the coordinate 3-volume r ≤ r 0 , t = t 0 . The unit normalÑ satisfiesÑ · n = 0 sõ N µ = (0,Ñ i ) where, Cramer's rule for matrix inverse, it can be shown that, where it should be noted thath < 0 and σ > 0. Deriving Eq. (B13) uses the fact thatH intersects Σ on ∂V 1 and therefore must have the same line element for variations in angular coordinates; hence g θθ =h θθ , g θφ =h θφ and g φφ =h φφ . The final component we need is to calculatẽ h tt which can be done by projecting the 4-metric g ontõ H as, = −α −2 , (B16) where (4)h is a 4-tensor belonging to M andÑ t = 0. Using the pushforward ofh onH to (4)h on M| p∈H , similarly to Sec. III, it can be shown that (4) h tt = h tt . Equations (B13) and (B16) combine to give, again notingh < 0. Now we can re-write the flux (B8) term as, F = αγ rν J νÑr + 1 √ γ rr β r Q, and this is identical to Eq. (43) found earlier.