Extraction of gravitational waves in numerical relativity

A numerical-relativity calculation yields in general a solution of the Einstein equations including also a radiative part, which is in practice computed in a region of finite extent. Since gravitational radiation is properly defined only at null infinity and in an appropriate coordinate system, the accurate estimation of the emitted gravitational waves represents an old and non-trivial problem in numerical relativity. A number of methods have been developed over the years to “extract” the radiative part of the solution from a numerical simulation and these include: quadrupole formulas, gauge-invariant metric perturbations, Weyl scalars, and characteristic extraction. We review and discuss each method, in terms of both its theoretical background as well as its implementation. Finally, we provide a brief comparison of the various methods in terms of their inherent advantages and disadvantages. Electronic supplementary material The online version of this article (doi:10.1007/s41114-016-0001-9) contains supplementary material, which is available to authorized users.


Introduction
With the commissioning of the second generation of laser interferometric gravitationalwave detectors, and the recent detection of gravitational waves (Abbott 2016), there is considerable interest in gravitational-wave astronomy. This is a huge field, covering the diverse topics of: detector hardware construction and design; data analysis; astrophysical source modeling; approximate methods for gravitational-wave calculation; and, when the weak field approach is not valid, numerical relativity.
Numerical relativity is concerned with the construction of a numerical solution to the Einstein equations, so obtaining an approximate description of a spacetime, and is reviewed, for example, in the textbooks by Alcubierre (2008), , Baumgarte and Shapiro (2010), Gourgoulhon (2012) and Rezzolla and Zanotti (2013). The physics in the simulation may be only gravity, as is the case of a binary black hole scenario, but it may also include matter fields and/or electromagnetic fields. Thus numerical relativity may be included in the modeling of a wide range of astrophysical processes. Often (but not always), an important desired outcome of the modeling process will be a prediction of the emitted gravitational waves. However, obtaining an accurate estimate of gravitational waves from the variables evolved in the simulation is normally a rather complicated process. The key difficulty is that gravitational waves are unambiguously defined only at future null infinity (J + ), whereas in practice the domain of numerical simulations is a region of finite extent using a "3+1" foliation of the spacetime. This is true for most of the numerical codes, but there are also notable exceptions. Indeed, there have been attempts towards the construction of codes that include both null infinity and the central dynamic region in the domain, but they have not been successful in the general case. These attempts include the hyperboloidal method (Frauendiener 2004), Cauchy characteristic matching (Winicour 2005), and a characteristic code (Bishop et al. 1997b). The only successful application to an astrophysical problem has been to axisymmetric core collapse using a characteristic code (Siebel et al. 2003).
In the linearized approximation, where gravitational fields are weak and velocities are small, it is straightforward to derive a relationship between the matter dynamics and the emission of gravitational waves, the well-known quadrupole formula. This can be traced back to work by Einstein (1916Einstein ( , 1918 shortly after the publication of general relativity. The method is widely used to estimate gravitational-wave production in many astrophysical processes. However, the strongest gravitational-wave signals come from highly compact systems with large velocities, that is from processes where the linearized assumptions do not apply. And of course, it is an event producing a powerful signal that is most likely to be found in gravitational-wave detector data. Thus it is important to be able to calculate gravitational-wave emission accurately for processes such as black hole or neutron star inspiral and merger, stellar core collapse, etc. Such problems cannot be solved analytically and instead are modeled by numerical relativity, as described in the previous paragraph, to compute the gravitational field near the source. The procedure of using this data to measure the gravitational radiation far from the source is called "extraction" of gravitational waves from the numerical solution.
In addition to the quadrupole formula and full numerical relativity, there are a number of other approaches to calculating gravitational-wave emission from astrophysical sources. These techniques are not discussed here and are reviewed elsewhere. They include post-Newtonian methods (Blanchet 2014), effective one-body methods (Damour and Nagar 2016), and self-force methods (Poisson et al. 2011). Another approach, now no-longer pursued, is the so-called "Lazarus approach", that combined analytical and numerical techniques (Baker et al. 2000b(Baker et al. , 2002a. In this article we will review a number of different extraction methods: (a) Quadrupole formula and its variations (Sect. 2.3); (b) methods using the Newman-Penrose scalar ψ 4 evaluated on a worldtube (Γ ) (Sect. 3.3); (c) Cauchy Perturbative methods, using data on Γ to construct an approximation to a perturbative solution on a known curved background (Sects. 4, 5; Evans 1988, 1990); and (d) Characteristic extraction, using data on Γ as inner boundary data for a characteristic code to find the waveform at J + (Sects. 6, 7). The description of the methods is fairly complete, with derivations given from first principles and in some detail. In cases (c) and (d), the theory involved is quite lengthy, so we also provide implementation summaries for the reader who is more interested in applying, rather than fully understanding, a particular method, see Sects. 5.6 and 7.8.
In addition, this review provides background material on gravitational waves (Sect. 2), on the "3+1" formalism for evolving the Einstein equations (Sect. 3), and on the characteristic formalism with particular reference to its use in estimating gravitational radiation (Sect. 6). The review concludes with a comparison of the various methods for extracting gravitational waves (Sect. 8). This review uses many different symbols, and their use and meaning is summarized in "Appendix 1". Spin-weighted, and other, spherical harmonics are discussed in "Appendix 2", and various computer algebra scripts and numerical codes are given in "Appendix 3".
Throughout, we will use a spacelike signature (−, +, +, +) and a system of geometrised units in which G = c = 1, although when needed we will also indi-cate the speed of light, c, explicitly. We will indicate with a boldface any tensor, e.g., V and with the standard arrow any three-dimensional vector or operator, e.g., v and ∇. Four-dimensional covariant and partial derivatives will be indicated in general with ∇ μ and ∂ μ , but other symbols may be introduced for less common definitions, or when we want to aid the comparison with classical Newtonian expressions. Within the standard convention of a summation of repeated indices, Greek letters will be taken to run from 0 to 3, while Latin indices run from 1 to 3.
We note that some of the material in this review has already appeared in books or other review articles. In particular, we have abundantly used parts of the text from the book "Relativistic Hydrodynamics", by Rezzolla and Zanotti (2013), from the review article "Gauge-invariant non-spherical metric perturbations of Schwarzschild black-hole spacetimes", by Nagar and Rezzolla (2006), as well as adaptations of the text from the article "Cauchy-characteristic matching", by Bishop et al. (1999a).

A quick review of gravitational waves 2.1 Linearized Einstein equations
When considering the Einstein equations as a set of second-order partial differential equations it is not easy to predict that there exist solutions behaving as waves. Indeed, the concept of gravitational waves as solutions of the Einstein equations written as linear and homogeneous wave equations is valid only under some rather idealised assumptions, such as a vacuum and asymptotically flat spacetime, a linearised regime for the gravitational fields and suitable gauges. If these assumptions are removed, the definition of gravitational waves becomes much more difficult, although still possible. It should be noted, however, that in this respect gravitational waves are not peculiar. Any wave-like phenomenon, in fact, can be described in terms of homogeneous wave equations only under simplified assumptions, such as those requiring a uniform "background" for the fields propagating as waves. These considerations suggest that the search for wave-like solutions to the Einstein equations should be made in a spacetime with very modest curvature and with a line element which is that of flat spacetime but for small deviations of nonzero curvature, i.e., where the linearised regime is guaranteed by the fact that |h μν | 1. Before writing the linearised version of the Einstein equations (1) it is necessary to derive the linearised expression for the Christoffel symbols. In a Cartesian coordinate basis (such as the one we will assume hereafter), we recall that the general expression for the affine connection is given by where the partial derivatives are readily calculated as As a result, the linearised Christoffel symbols become Note that the operation of lowering and raising the indices in expression (5) is not made through the metric tensors g μν and g μν but, rather, through the spacetime metric tensors η μν and η μν . This is just the consequence of the linearised approximation and, despite this, the spacetime is really curved! Once the linearised Christoffel symbols have been computed, it is possible to derive the linearised expression for the Ricci tensor which takes the form where is the trace of the metric perturbations. The resulting Ricci scalar is then given by R := g μν R μν η μν R μν .
Making use of (6) and (8), it is possible to rewrite the Einstein equations (1) in a linearised form as Although linearised, the Einstein equations (9) do not yet seem to suggest a wavelike behaviour. A good step in the direction of unveiling this behaviour can be made if we introduce a more compact notation, which makes use of "trace-free" tensors defined ash where the "bar-operator" in (10) can be applied to any symmetric tensor so that, for instance,R μν = G μν , and also iteratively, i.e.,h μν = h μν . 1 Using this notation, the linearised Einstein equations (9) take the more compact form − ∂ α ∂ αhμν − η μν ∂ α ∂ βh αβ + ∂ α ∂ μhνα = 16π T μν , where the first term on the left-hand side of (11) can be easily recognised as the Dalambertian (or wave) operator, i.e., ∂ α ∂ αh μν = h μν . At this stage, we can exploit the gauge freedom inherent in general relativity (see also below for an extended discussion) to recast Eq. (11) in a more convenient form. More specifically, we exploit this gauge freedom by choosing the metric perturbations h μν so as to eliminate the terms in (11) that spoil the wave-like structure. Most notably, the coordinates can be selected so that the metric perturbations satisfy Making use of the gauge (12), which is also known as the Lorenz (or Hilbert) gauge, the linearised field equations take the form that, in vacuum reduce to the desired result Equations (14) show that, in the Lorenz gauge and in vacuum, the metric perturbations propagate as waves distorting flat spacetime.
The simplest solution to the linearised Einstein equations (14) is that of a plane wave of the typeh μν = A μν exp(iκ α x α ), where of course we are interested only in the real part of (15), with A being the amplitude tensor. Substitution of the ansatz (15) into Eq. (14) implies that κ α κ α = 0 so that κ is a null four-vector. In such a solution, the plane wave (15) travels in the spatial direction k = (κ x , κ y , κ z )/κ 0 with frequency ω := κ 0 = (κ j κ j ) 1/2 . The next step is to substitute the ansatz (15) into the Lorenz gauge condition Eq. (12), yielding A μν κ μ = 0 so that A and κ are orthogonal. Consequently, the amplitude tensor A, which in principle has 16 −6 = 10 independent components, satisfies four conditions. Thus the imposition of the Lorenz gauge reduces the independent components of A to six. We now investigate how to reduce the number of independent components to match the number of dynamical degrees of freedom of general relativity, i.e., two. While a Lorenz gauge has been imposed [cf. Eq. (12)], this does not completely fix the coordinate system of a linearised theory. A residual ambiguity, in fact, is preserved through arbitrary gauge changes, i.e., through infinitesimal coordinate transformations that are consistent with the gauge that has been selected. The freedom to make such a transformation follows from a foundation of general relativity, the principle of general covariance. To better appreciate this matter, consider an infinitesimal coordinate transformation in terms of a small but otherwise arbitrary displacement four-vector ξ Applying this transformation to the linearised metric (2) generates a "new" metric tensor that, to the lowest order, is so that the "new" and "old" perturbations are related by the following expression or, alternatively, byh new μ ν =h old μν − ∂ ν ξ μ − ∂ μ ξ ν + η μν ∂ α ξ α .
Requiring now that the new coordinates satisfy the condition (12) of the Lorenz gauge ∂ αhnew μα = 0, forces the displacement vector to be solution of the homogeneous wave equation ∂ β ∂ β ξ α = 0.
As a result, the plane-wave vector with components generates, through the four arbitrary constants C α , a gauge transformation that changes arbitrarily four components of A in addition to those coming from the condition A · κ = 0. Effectively, therefore, A μν has only 10 − 4 − 4 = 2 linearly independent components, corresponding to the number of degrees of freedom in general relativity (Misner et al. 1973). Note that these considerations are not unique to general relativity and similar arguments can also be made in classical electrodynamics, where the Maxwell equations are invariant under transformations of the vector potentials of the type A μ → A μ = A μ + ∂ μ Ψ , where Ψ is an arbitrary scalar function, so that the corresponding electromagnetic tensor is F new μ ν = ∂ ν A μ − ∂ μ A ν = F old μ ν . Similarly, in a linearised theory of general relativity, the gauge transformation (18) will preserve the components of the Riemann tensor, i.e., R new αβμν = R old αβμν + O(R 2 ). To summarise, it is convenient to constrain the components of the amplitude tensor through the following conditions: (a) orthogonality condition: four components of the amplitude tensor can be specified since the Lorenz gauge implies that A and κ are orthogonal, i.e., A μν κ ν = 0. (b) choice of observer: three components of the amplitude tensor can be eliminated after selecting the infinitesimal displacement vector ξ μ = iC μ exp(iκ α x α ) so that A μν u μ = 0 for some chosen four-velocity vector u. This means that the coordinates are chosen so that for an observer with four-velocity u μ the gravitational wave has an effect only in spatial directions. 2 (c) traceless condition: one final component of the amplitude tensor can be eliminated after selecting the infinitesimal displacement vector ξ μ = iC μ exp(iκ α x α ) so that A μ μ = 0.
Conditions (a), (b) and (c) define the so-called transverse-traceless (TT) gauge, which represents a most convenient gauge for the analysis of gravitational waves. To appreciate the significance of these conditions, consider them implemented in a reference frame which is globally at rest, i.e., with four-velocity u α = (1, 0, 0, 0), where the amplitude tensor must satisfy: i.e., the spatial components of h μν are divergence-free.
i.e., only the spatial components of h μν are nonzero, hence the transverse character of the TT gauge.
i.e., the spatial components of h μν are trace free hence the trace-free character of the TT gauge. Because of this, and only in this gauge,h μν = h μν .

Making sense of the TT gauge
As introduced so far, the TT gauge might appear rather abstract and not particularly interesting. Quite the opposite, the TT gauge introduces a number of important advantages and simplifications in the study of gravitational waves. The most important of these is that, in this gauge, the only nonzero components of the Riemann tensor are However, since the use of the TT gauge indicates that a travelling gravitational wave with periodic time behaviour h TT jk ∝ exp(iωt) can be associated to a local oscillation of the spacetime, i.e., To better appreciate the effects of the propagation of a gravitational wave, it is useful to consider the separation between two neighbouring particles A and B on a geodesic motion and how this separation changes in the presence of an incident gravitational wave (see Fig. 1). For this purpose, let us introduce a coordinate system xα in the Schematic diagram for the changes in the separation vector n between two particles A and B moving along geodesic trajectories produced by the interaction with a gravitational wave propagating along the direction κ neighbourhood of particle A so that along the worldline of the particle A the line element will have the form The arrival of a gravitational wave will perturb the geodesic motion of the two particles and produce a nonzero contribution to the geodesic-deviation equation. We recall that the changes in the separation four-vector n between two geodesic trajectories with tangent four-vector u are expressed through the geodesic-deviation equation (see Fig. 1) where the operator D Dτ is the covariant time derivative along the worldline (in this case a geodesic) of a particle.
Indicating now with nˆj B := xˆj B − xˆj A = xˆj B the components of the separation threevector in the positions of the two particles, the geodesic-deviation equation (29) can be written as A first simplification to these equations comes from the fact that around the particle A, the affine connections vanish (i.e., Γˆĵ αβ = 0) and the covariant derivative in (31) can be replaced by an ordinary total derivative. Furthermore, because in the TT gauge the coordinate system xα moves together with the particle A, the proper and the coordinate time coincide at first order in the metric perturbation [i.e., τ = t + O((h TT μν ) 2 )]. As a result, equation (31) effectively becomes and has solution Equation (33) has a straightforward interpretation and indicates that, in the reference frame comoving with A, the particle B is seen oscillating with an amplitude proportional to h TT jk .
Note that because these are transverse waves, they will produce a local deformation of the spacetime only in the plane orthogonal to their direction of propagation. As a result, if the two particles lay along the direction of propagation (i.e., if n κ), then h Let us now consider a concrete example and in particular a planar gravitational wave propagating in the positive z-direction. In this case where A + , A × represent the two independent modes of polarization, and the symbol refers to the real part. As in classical electromagnetism, in fact, it is possible to decompose a gravitational wave in two linearly polarized plane waves or in two circularly polarized ones. In the first case, and for a gravitational wave propagating in the z-direction, the polarization tensors + ("plus") and × ("cross") are defined as The deformations that are associated with these two modes of linear polarization are shown in Fig. 2 where the positions of a ring of freely-falling particles are schematically represented at different fractions of an oscillation period. Note that the two linear polarization modes are simply rotated of π/4. In a similar way, it is possible to define two tensors describing the two states of circular polarization and indicate with e R the circular polarization that rotates clockwise (see Fig. 3) and with e L the circular polarization that rotates counter-clockwise (see Fig. 3) The deformations that are associated to these two modes of circular polarization are shown in Fig. 3.

The quadrupole formula
The quadrupole formula and its domain of applicability were mentioned in Sect. 1, and some examples of its use in a numerical simulation are presented in Sect. 8. In practice, the quadrupole formula represents a low-velocity, weak-field approximation to measure the gravitational-wave emission within a purely Newtonian description of gravity. 3 In practice, the formula is employed in those numerical simulations that either treat gravity in an approximate manner (e.g., via a post-Newtonian approximation or a conformally flat metric) or that, although in full general relativity, have computational domains that are too small for an accurate calculation of the radiative emission.
In what follows we briefly discuss the amounts of energy carried by gravitational waves and provide simple expressions to estimate the gravitational radiation luminosity of potential sources. Although the estimates made here come from analogies with electromagnetism, they provide a reasonable approximation to more accurate expressions from which they differ for factors of a few. Note also that while obtaining such a level of accuracy requires only a small effort, reaching the accuracy required of a template to be used in the realistic detection of gravitational waves is far more difficult and often imposes the use of numerical relativity calculations on modern supercomputers.
In classical electrodynamics, the energy emitted per unit time by an oscillating electric dipole d = qx, with q the electrical charges and x their separation, is easily estimated to be where the number of "dots" counts the order of the total time derivative. Equally simple is to calculate the corresponding luminosity in gravitational waves produced by an oscillating mass dipole. In the case of a system of N point-like particles of mass m A (A = 1, 2, . . . , N ), in fact, the total mass dipole and its first time derivative are respectively. However, the requirement that the system conserves its total linear momentumd :=ṗ = 0, forces to conclude that L mass dip. = 0, i.e., that there is no mass-dipole radiation in general relativity (This is equivalent to the impossibility of having electromagnetic radiation from an electric monopole oscillating in time.). Next, consider the electromagnetic energy emission produced by an oscillating electric quadrupole. In classical electrodynamics, this energy loss is given by where is the electric quadrupole for a distribution of N charges (q 1 , q 2 , . . . , q N ). In close analogy with expression (44), the energy loss per unit time due to an oscillating mass quadrupole is calculated to be L mass quad. := 1 5 ( ... where I − jk is the trace-less mass quadrupole (or "reduced" mass quadrupole), defined as and the brackets indicate a time average [Clearly, the second expression in (47) refers to a continuous distribution of particles with rest-mass density ρ.].
A crude estimate of the third time derivative of the mass quadrupole of the system is given by where v is the mean internal velocity. Stated differently, ...
where L int is the power of the system flowing from one part of the system to the other. As a result, the gravitational-wave luminosity in the quadrupole approximation can be calculated to be (we here restore the explicit use of the gravitational constant and of the speed of light) (50) The second equality has been derived using the virial theorem for which the kinetic energy is of the same order of the potential one, i.e., M v 2 ∼ G M 2 /R, and assuming that the oscillation timescale is inversely proportional to the mean stellar density, i.e., τ ∼ (1/G ρ ) 1/2 ∼ (R 3 /G M) 1/2 . Similarly, the third equality expresses the luminosity in terms of dimensionless quantities such as the size of the source relative to the Schwarzschild radius R S = 2G M/c 2 , and the source speed in units of the speed of light. Note that the quantity c 5 G has indeed the units of a luminosity, i.e., erg s −1 = cm 2 g s −3 in cgs units.
Although extremely simplified, expressions (48) and (50) contain the two most important pieces of information about the generation of gravitational waves. The first one is that the conversion of any type of energy into gravitational waves is, in general, not efficient. To see this it is necessary to bear in mind that expression (46) is in geometrized units and that the conversion to conventional units, say cgs units, requires dividing (46) by the very large factor c 5 /G 3.63 × 10 59 erg s −1 . The second one is contained in the last expression in Eq. (50) and that highlights how the gravitationalwave luminosity can also be extremely large. There are in fact astrophysical situations, such as those right before the merger of a binary system of compact objects, in which v 2 ∼ 0.1 c and R ∼ 10 R S , so that L mass−quad ∼ 10 51 erg s −1 ∼ 10 18 L , that is, 10 18 times the luminosity of the Sun; this is surely an impressive release of energy.

Extensions of the quadrupole formula
Although valid only in the low-velocity, weak-field limit, the quadrupole-formula approximation has been used extensively in the past and still finds use in several simulations, ranging from stellar-core collapse (see, e.g., Zwerger and Müller 1997 for some initial application) to binary neutron-star mergers (see, e.g., Oechslin et al. 2002 for some initial application). In many of the simulations carried out to study stellar collapse, one makes the additional assumption that the system remains axisymmetric and the presence of an azimuthal Killing vector has two important consequences. Firstly, the gravitational waves produced in this case will carry away energy but not angular momentum, which is a conserved quantity in this spacetime. Secondly, the gravitational waves produced will have a single polarization state, so that the transverse traceless gravitational field is completely determined in terms of its only nonzero transverse and traceless (TT) independent component. Following Zwerger and Müller (1997) and considering for simplicity an axisymmetric system, it is useful to express the gravitational strain h TT (t) observed at a distance R from the source in terms of the quadrupole wave amplitude A 20 (Zanotti et al. 2003) h where F + = F + (R, θ, φ) is the detector's beam pattern function and depends on the orientation of the source with respect to the observer. As customary in these calculations, we will assume it to be optimal, i.e., F + = 1. The = 2, m = 0 wave amplitude A 20 in Eq. (51) is simply the second time derivative of the reduced mass quadrupole moment in axisymmetry and can effectively be calculated without taking time derivatives numerically, which are instead replaced by spatial derivatives of evolved quantities after exploiting the continuity and the Euler equations (Finn and Evans 1990;Blanchet et al. 1990;Rezzolla et al. 1999b). The result in a spherical coordinate system is where z := cos θ , k = 16π 3/2 / √ 15, is the Newtonian gravitational potential, and I [ax] is the appropriate component of the Newtonian reduced mass-quadrupole moment in axisymmetry Of course it is possible to consider more generic conditions and derive expressions for the strain coming from more realistic sources, such as a an astrophysical system with equatorial symmetry. In this case, focussing on the lowest = 2 moments, the relevant multipolar components for the strain are Baiotti et al. (2009) where are the more general Newtonian mass and mass-current quadrupoles. Expressions (52)-(58) are strictly Newtonian. Yet, these expression are often implemented in numerical codes that are either fully general relativistic or exploit some level of general-relativistic approximation. More seriously, these expressions completely ignore considerations that emerge in a relativistic context, such as the significance of the coordinate chosen for their calculation. As a way to resolve these inconsistencies, improvements to these expressions have been made to increase the accuracy of the computed gravitational-wave emission. For instance, for calculations on known spacetime metrics, the gravitational potential in expression (52) is often approximated with expressions derived from the metric, e.g., as = (1 − g rr )/2 (Zanotti et al. 2003), which is correct to the first Post-Newtonian (PN) order. Improvements to the mass quadrupole (53) inspired by a similar spirit have been computed in Blanchet et al. (1990), and further refined and tested in Shibata and Sekiguchi (2004), Nagar et al. (2005), Cerdá-Durán et al. (2005), Pazos et al. (2007), Baiotti et al. (2007), Dimmelmeier et al. (2007), Corvino et al. (2010). A systematic comparison among the different expressions of the quadrupole formulas developed over the years was carried out in Baiotti et al. (2009), where a generalization of the mass-quadrupole formula (57) was introduced. In essence, following previous work in Nagar et al. (2005), Baiotti et al. (2009) introduced a "generalized" mass-quadrupole moment of the form where the generalized rest-mass density was assumed to take a number of possible expressions, namely, Clearly, the first option corresponds to the "standard" quadrupole formula, but, as remarked in Baiotti et al. (2009) none of the alternative quadrupole formulas obtained using these generalized quadrupole moments should be considered better than the others, at least mathematically. None of them is gauge invariant and indeed they yield different results depending on the underlining choice made for the coordinates. Yet, the comparison is meaningful in that these expressions were and still are in use in many numerical codes, and it is therefore useful to determine which expression is effectively closer to the fully general-relativistic one.
Making use of a fully general-relativistic measurement of the gravitational-wave emission from a neutron star oscillating nonradially as a result of an initial pressure perturbation, Baiotti et al. (2009) concluded that the various quadrupole formulas are comparable and give a very good approximation to the phasing of the gravitationalwave signals. At the same time, they also suffer from systematic over-estimate [expression (61)] or under-estimates of the gravitational-wave amplitude [expressions (60), and (62)- (63)]. In all cases, however, the relative difference in amplitude was of 50 % at most, which is probably acceptable given that these formulas are usually employed in complex astrophysical calculations in which the systematic errors coming from the microphysical modelling are often much larger.

The 3+1 decomposition of spacetime
At the heart of Einstein's theory of general relativity is the equivalence among all coordinates, so that the distinction of spatial and time coordinates is more an organisational matter than a requirement of the theory. Despite this "covariant view", however, our experience, and the laws of physics on sufficiently large scales, do suggest that a distinction of the time coordinate from the spatial ones is the most natural one in describing physical processes. Furthermore, while not strictly necessary, such a distinction of time and space is the simplest way to exploit a large literature on the numerical solution of hyperbolic partial differential equations as those of relativistic hydrodynamics. In a generic spacetime, analytic solutions to the Einstein equations are not known, and a numerical approach is often the only way to obtain an estimate of the solution.
Following this principle, a decomposition of spacetime into "space" and "time" was already proposed in the 1960s within a Hamiltonian formulation of general relativity and later as an aid to the numerical solution of the Einstein equations in vacuum. The basic idea is rather simple and consists in "foliating" spacetime in terms of a set of non-intersecting spacelike hypersurfaces Σ := Σ(t), each of which is parameterised by a constant value of the coordinate t. In this way, the three spatial coordinates are split from the one temporal coordinate and the resulting construction is called the 3+1 decomposition of spacetime (Misner et al. 1973).
Given one such constant-time hypersurface, Σ t , belonging to the foliation Σ, we can introduce a timelike four-vector n normal to the hypersurface at each event in the spacetime and such that its dual one-form Ω := ∇t is parallel to the gradient of the coordinate t, i.e., with n μ = {A, 0, 0, 0} and A a constant to be determined. If we now require that the four-vector n defines an observer and thus that it measures the corresponding fourvelocity, then from the normalisation condition on timelike four-vectors, n μ n μ = −1, we find that n μ n μ = g μν n μ n ν = g tt A 2 = − 1 where we have defined α 2 := −1/g tt . From the last equality in expression (65) it follows that A = ±α and we will select A = −α, such that the associated vector field n μ is future directed. The quantity α is commonly referred to as the lapse function, it measures the rate of change of the coordinate time along the vector n μ (see Fig. 4), and will be a building block of the metric in a 3+1 decomposition [cf. Eq. (72)]. The specification of the normal vector n allows us to define the metric associated to each hypersurface, i.e., γ μν := g μν + n μ n ν , γ μν := g μν + n μ n ν , where γ 0μ = 0, γ i j = g i j , but in general γ i j = g i j . Also note that γ ik γ k j = δ i j , that is, γ i j and γ i j are the inverse of each other, and so can be used for raising and lowering the indices of purely spatial vectors and tensors (that is, defined on the hypersurface Σ t ).
The tensors n and γ provide us with two useful tools to decompose any fourdimensional tensor into a purely spatial part (hence contained in the hypersurface Σ t ) and a purely timelike part (hence orthogonal to Σ t and aligned with n). Not surprisingly, the spatial part is readily obtained after contracting with the spatial projection operator (or spatial projection tensor) Fig. 4 Schematic representation of the 3+1 decomposition of spacetime with hypersurfaces of constant time coordinate Σ t and Σ t+dt foliating the spacetime. The four-vector t represents the direction of evolution of the time coordinate t and can be split into a timelike component αn, where n is a timelike unit normal to the hypersurface, and into a spacelike component, represented by the spacelike four-vector β. The function α is the "lapse" and measures the proper time between adjacent hypersurfaces, while the components of the "shift" vector β i measure the change of coordinates from one hypersurface to the subsequent one while the timelike part is obtained after contracting with the time projection operator (or time projection tensor) N μ ν := −n μ n ν , and where the two projectors are obviously orthogonal, i.e., We can now introduce a new vector, t, along which to carry out the time evolutions and that is dual to the surface one-form Ω. Such a vector is just the time-coordinate basis vector and is defined as the linear superposition of a purely temporal part (parallel to n) and of a purely spatial one (orthogonal to n), namely The purely spatial vector β [i.e., β μ = (0, β i )] is usually referred to as the shift vector and will be another building block of the metric in a 3+1 decomposition [cf. Eq. (72)]. The decomposition of the vector t into a timelike component nα and a spatial component β is shown in Fig. 4.
Because t is a coordinate basis vector, the integral curves of t μ are naturally parameterised by the time coordinate. As a result, all infinitesimal vectors t μ originating at a given point x i 0 on one hypersurface Σ t would end up on the hypersurface Σ t+dt at a point whose coordinates are also x i 0 . This condition is not guaranteed for translations along Ω μ unless β μ = 0 since t μ t μ = g tt = −α 2 + β μ β μ , and as illustrated in Fig. 4.
In summary, the components of n are given by and we are now ready to deduce that the lapse function and the shift vector can be employed to express the generic line element in a 3+1 decomposition as Expression (72) clearly emphasises that when β i = 0 = dx i , the lapse measures the proper time, dτ 2 = −ds 2 , between two adjacent hypersurfaces, i.e., while the shift vector measures the change of coordinates of a point that is moved along n from the hypersurface Σ t to the hypersurface Σ t+dt , i.e., Similarly, the covariant and contravariant components of the metric (72) can be written explicitly as from which it is easy to obtain an important identity which will be used extensively hereafter, i.e., where g := det(g μν ) and γ := det(γ i j ). When defining the unit timelike normal n in Eq. (65), we have mentioned that it can be associated to the four-velocity of a special class of observers, which are referred to as normal or Eulerian observers. Although this denomination is somewhat confusing, since such observers are not at rest with respect to infinity but have a coordinate velocity dx i /dt = n i = −β i /α, we will adopt this traditional nomenclature also in the following and thus take an "Eulerian observer" as one with four-velocity given by (71).
When considering a fluid with four-velocity u, the spatial four-velocity v measured by an Eulerian observer will be given by the ratio between the projection of u in the space orthogonal to n, i.e., γ i μ u μ = u i , and the Lorentz factor W of u as measured by n (de Felice and Clarke 1990) − n μ u μ = αu t = W.
As a result, the spatial four-velocity of a fluid as measured by an Eulerian observer will be given by Using now the normalisation condition u μ u μ = −1, we obtain so that the components of v can be written as where in the last equality we have exploited the fact that γ i j u j = u i − β i W/α.

The ADM formalism: 3+1 decomposition of the Einstein equations
The 3+1 decomposition introduced in Sect. 3.1 can be used not only to decompose tensors, but also equations and, in particular, the Einstein equations, which are then cast into an initial-value form suitable to be solved numerically. A 3+1 decomposition of the Einstein equations was presented by Arnowitt et al. (2008), but it is really the reformulation suggested by York (1979) that represents what is now widely known as the ADM formulation (see, e.g., Alcubierre 2008; Gourgoulhon 2012 for a detailed and historical discussion). As we will see in detail later on, in this formulation the Einstein equations are written in terms of purely spatial tensors that can be integrated forward in time once some constraints are satisfied initially.
Here, we only outline the ADM formalism, and refer to the literature for the derivation and justification. Further, it is important to note that the ADM formulation is, nowadays, not used in practice because it is only weakly hyperbolic. However, the variables used in the ADM method, in particular the three-metric and the extrinsic curvature, are what will be needed later for gravitational-wave extraction, and are easily obtained from the output of other evolution methods (see discussion in Sects. 5, 7).
Instead of the ADM formalism, modern simulations mainly formulate the Einstein equations using: the BSSNOK method (Nakamura et al. 1987;Shibata and Nakamura 1995;Baumgarte and Shapiro 1999); the CCZ4 formulation (Alic et al. 2012), which was developed from the Z4 method (Bona et al. 2003(Bona et al. , 2004) (see also Bernuzzi and Hilditch 2010 for the so-called Z4c formulation and Alic et al. 2013 for some comparisons); or the generalized harmonic method (Pretorius 2005) (see also Baumgarte and Shapiro 2010;Rezzolla and Zanotti 2013 for more details).
We start by noting that once a 3+1 decomposition is introduced as discussed in Sect. 3.1, it is then possible to define the three-dimensional covariant derivative D i . Formally, this is done by projecting the standard covariant derivative onto the space orthogonal to n μ , and the result is a covariant derivative defined with respect to the connection coefficients where we will use the upper left index (3) to mark a purely spatial quantity that needs to be distinguished from its spacetime counterpart. 4 Similarly, the three-dimensional Riemann tensor (3) R i jk associated with γ has an explicit expression given by In a similar manner, the three-dimensional contractions of the three-dimensional Riemann tensor, i.e., the three-dimensional Ricci tensor and the three-dimensional Ricci scalar, are defined respectively as their four-dimensional counterparts, i.e., The information present in R μ νκσ and missing in (3) R i jk can be found in another symmetric tensor, the extrinsic curvature K i j , which is purely spatial. Loosely speaking, the extrinsic curvature provides a measure of how the three-dimensional hypersurface Σ t is curved with respect to the four-dimensional spacetime. For our purposes, it is convenient to define the extrinsic curvature as (but note that other definitions, which can be shown to be equivalent, are common) where L n is the Lie derivative relative to the normal vector field n. Expression (84) provides a simple interpretation of the extrinsic curvature K i j as the rate of change of the three-metric γ i j as measured by an Eulerian observer. Using properties of the Lie derivative, it follows that Note that Eq. (85) is a geometrical result and is independent of the Einstein equations. The next step is to note further purely geometric relations, how the spacetime curvature is related to the intrinsic and extrinsic curvatures of the hypersurface Σ t . These formulas are known as the Gauss-Codazzi equations and the Codazzi-Mainardi equations. They are We now have enough identities to rewrite the Einstein equations in a 3+1 decomposition. After contraction, we can use the Einstein equations to replace the spacetime Ricci tensor with terms involving the stress-energy tensor, and then after further manipulation the final result is: (3) The following definitions have been made for the "matter" quantities that is, for contractions of the energy-momentum tensor that would obviously be zero in vacuum spacetimes. Overall, the six equations (89), together with the six equations (85) represent the time-evolving part of the ADM equations and prescribe how the three-metric and the extrinsic curvature change from one hypersurface to the following one. In contrast, Eqs. (90) and (91) are constraints that need to be satisfied on each hypersurface. This distinction into evolution equations and constraint equations is not unique to the ADM formulation and is indeed present also in classical electromagnetism. Just as in electrodynamics the divergence of the magnetic field remains zero if the field is divergence-free at the initial time, so the constraint equations (90) and (91), by virtue of the Bianchi identities (Alcubierre 2008;Baumgarte and Shapiro 2010;Gourgoulhon 2012;Rezzolla and Zanotti 2013), will remain satisfied during the evolution if they are satisfied initially (Frittelli 1997). Of course, this concept is strictly true in the continuum limit, while numerically the situation is rather different. However, that issue is not pursued here.
Two remarks should be made before concluding this section. The first one is about the gauge quantities, namely, the lapse function α and the shift vector β i . Since they represent the four degrees of freedom of general relativity, they are not specified by the equations discussed above and indeed they can be prescribed arbitrarily, although in practice great care must be taken in deciding which prescription is the most useful. The second comment is about the mathematical properties of the time-evolution ADM equations (89) and (85). The analysis of these properties can be found, for instance, in Reula (1998) or in Frittelli and Gómez (2000), and reveals that such a system is only weakly hyperbolic with zero eigenvalues and, as such, not necessarily wellposed. The weak-hyperbolicity of the ADM equations explains why, while an historical cornerstone in the 3+1 formulation of the Einstein equations, they are rarely used in practice and have met only limited successes in multidimensional calculations (Cook et al. 1998;Abrahams et al. 1998). At the same time, the weak hyperbolicity of the ADM equations and the difficulty in obtaining stable evolutions, has motivated, and still motivates, the search for alternative formulations.

Gravitational waves from ψ 4 on a finite worldtube(s)
The Newman-Penrose scalars are scalar quantities defined as contractions between the Weyl, or conformal, tensor (in four dimensions), and an orthonormal null tetrad α , n α [N P] , m α ,m α (Newman and Penrose 1963). The null tetrad is constructed from an orthonormal tetrad, and we use the notation n α [N P] , rather than the usual n α , because n α [N P] is obtained in terms of the hypersurface normal n α . Supposing that the spatial coordinates (x, y, z) are approximately Cartesian, then spherical polar coordinates are defined using Eqs. (399) and (400). However, in general these coordinates are not exactly spherical polar, and in particular the radial coordinate is not a surface area coordinate (for which the 2-surface r = t = constant must have area 4πr 2 ). We reserve the notation r for a surface area radial coordinate, so the radial coordinate just constructed will be denoted by s. Then the outward-pointing radial unit normal e s is An orthonormal basis (e s , e θ , e φ ) of Σ t is obtained by Gram-Schmidt orthogonalization, and is extended to be an orthonormal tetrad of the spacetime by incuding the hypersurface normal n. Then the orthonormal null tetrad is (The reader should be aware that some authors use different conventions, e.g., without a factor of √ 2, leading to different forms for various equations). For example, in Minkowski spacetime there is no distinction between s and r , and in spherical polar The null tetrad satisfies the orthonormality conditions (97) The Newman-Penrose, or Weyl, scalars (Newman and Penrose 1963) are defined as For our purposes, the most important of these quantities is ψ 4 since in the asymptotic limit it completely describes the outgoing gravitational radiation field: far from a source, a gravitational wave is locally plane and ψ 4 is directly related to the metric perturbation in the TT gauge In an asymptotically flat spacetime using appropriate coordinates (these issues are discussed more formally in Sect. 6), the peeling theorem (Penrose 1965a;Geroch 1977;Hinder et al. 2011) shows that ψ 4 falls off as r −1 , and more generally that ψ n falls off as r n−5 . Thus, gravitational waves are normally described not by ψ 4 but by r ψ 4 which should be evaluated in the limit as r → ∞ (which in practice may mean evaluated at as large a value of r as is feasible). Often lim r →∞ r ψ 4 is denoted by ψ 0 4 (t, θ, φ), but that notation will not be used in this section. These issues are discussed further in Sect. 6, but for now we will regard gravitational waves, and specifically r ψ 4 , as properly defined only in a spacetime whose metric can be written in a form that tends to the Minkowski metric, and for which the appropriate definition of the null tetrad is one that tends to the form Eq. (96), as r → ∞.
Equation (102) for ψ 4 involves spacetime, rather than hypersurface, quantities, and this is not convenient in a "3+1" simulation. However, the expression for ψ 4 can be manipulated into a form involving hypersurface quantities only (Gunnarsen et al. 1995) (there is also a derivation in the textbook (Alcubierre 2008), but note the sign difference in the definition of ψ 4 used there): The proof is not given here, but in summary is based on using an arbitrary timelike vector, in this case the hypersurface normal n, to decompose the Weyl tensor into its "electric" and "magnetic" parts.
The above procedures lead to an estimate ψ 4 , but results are rarely reported in this form. Instead, ψ 4 is decomposed into spin-weighted spherical harmonics (see section "The spin-weighted spherical harmonics s Y m " in "Appendix 2"), and the r ψ m 4 are evaluated and reported. Although, normally, the dominant part of a gravitational-wave signal is in the lowest modes with = 2, the other modes are important to gravitational-wave data analysis, recoil calculations, etc.

Extracting gravitational waves using ψ 4 on a finite worldtube
"3+1" numerical simulations are restricted to a finite domain, so it is not normally possible to calculate exactly a quantity given by an asymptotic formula (but see Sects. 6, 7). A simple estimate of r ψ 4 can be obtained by constructing coordinates (s, θ, φ) and an angular null tetrad vector m as discussed at the beginning of Sect. 3.3. Then r ψ 4 can be evaluated using Eq. (104) on a worldtube s = constant, and the estimate is r ψ 4 = sψ 4 or alternatively r ψ 4 = ψ 4 √ A/4π , where A is the area of the worldtube at time t. This approach was first used in Smarr (1977), and subsequently in, for example, Pollney et al. (2007), Pfeiffer (2007) and Scheel et al. (2009). This method does not give a unique answer, and there are many variations in the details of its implementation. However, the various estimates obtained for r ψ 4 should differ by no more than O(r −1 ).
The quantity ψ 4 has no free indices and so tensorially is a scalar, but its value does depend on the choice of tetrad. However, it may be shown that ψ 4 is first-order tetradinvariant if the tetrad is a small perturbation about a natural tetrad of the Kerr spacetime. This result was shown by Teukolsky (1972Teukolsky ( , 1973; see also Chandrasekhar (1978), and Campanelli et al. (2000). Briefly, the reasoning is as follows. The Kinnersley null tetrad is an exact null tetrad field in the Kerr geometry (Kinnersley 1969). It has the required asymptotic limit, and the vectors α , n α [N P] are generators of outgoing and ingoing radial null geodesics respectively. In the Kerr geometry C [Kerr] αβμν = 0, but using the Kinnersley tetrad all ψ n are zero except ψ 2 . Thus, to first-order, ψ 4 is evaluated using the perturbed Weyl tensor and the background tetrad; provided terms of the form C [Kerr] αβμν n α where three of the tetrad vectors take background values and only one is perturbed, are ignorable. Allowing for those ψ n that are zero, and using the symmetry properties of the Weyl tensor, all such terms vanish. This implies that the ambiguity in the choice of tetrad is of limited importance because it is a secondorder effect; see also ; . These ideas have been used to develop analytic methods for estimating ψ 4 Baker et al. 2000a;Baker and Campanelli 2000;Baker et al. 2001). Further, the Kinnersley tetrad is the staring point for a numerical extraction procedure.
In practice the spacetime being evolved is not Kerr, but in many cases at least far from the source it should be Kerr plus a small perturbation, and in the far future it should tend to Kerr. Thus an idea for an appropriate tetrad for use on a finite worldtube is to construct an approximation to the Kinnersley form, now known as the quasi-Kinnersley null tetrad Nerozzi et al. 2005). The quasi-Kinnersley tetrad has the property that as the spacetime tends to Kerr, then the quasi-Kinnersley tetrad tends to the Kinnersley tetrad. The method was used in a number of applications in the mid-2000s (Nerozzi et al. 2006;Campanelli et al. 2006;Fiske et al. 2005;Nerozzi 2007).
Despite the mathematical attraction of the quasi-Kinnersley approach, nowadays the extrapolation method which assumes the simpler Schwarzschild background (see next section) is preferred since, at a practical level, and as discussed in Sect. 8, extrapolation can give highly accurate results. Modern simulations typically extract on a worldtube at between 100 and 1000 M where the correction due to the background being Kerr rather than Schwarzschild is negligible. More precisely, an invariant measure of curvature is the square root of the Kretschmann scalar, which for the Kerr geometry (Henry 2000) takes the asymptotic form where a := J/M 2 . The curvature is already small in the Schwarzschild (a = 0) case, and the effect of ignoring a is a small relative error of order a 2 /r 2 .

Extracting gravitational waves using ψ 4 in practice: the extrapolation method
The method most commonly used at present is an adaptation of a simple estimate on a finite worldtube, and has become known as the extrapolation method. A schematic illustration of the method is given in Fig. 6. A preliminary version of extrapolation was used in 2005 (Baker et al. 2006). However, the method, as used at present, was developed in 2009 by two different groups Boyle and Mroué 2009), and a recent description is given by Taylor et al. (2013); see also Pollney et al. (2011). The essential idea is that ψ 4 is estimated on worldtubes at a number of different radii, and then the data is fitted to a polynomial of form ψ 4 = N n=1 A n /r n so that lim r →∞ r ψ 4 is approximated by A 1 . However, there are some subtleties that complicate the procedure a little. The expected polynomial form of ψ 4 is applicable only on an outgoing null cone; and further r should be a surface area coordinate (although often requiring this property is not important). We assume that the data available is ψ m 4 (t, s) obtained by decomposing ψ 4 into spherical harmonic components on a spherical surface of fixed coordinate radius s at a given coordinate time t. Then the first step in extrapolation is to obtain ψ m 4 (t * , r ), where t * is a retarted time coordinate specified in Eq. (108) below, and where Schematic illustration of ψ 4 extrapolation. The Cauchy evolution is shown with green slices, with an outer boundary in blue subject to a boundary condition that excludes incoming gravitational waves. The light blue lines are approximations to outgoing null slices. ψ 4 is evaluated where the Cauchy slices meet the innermost worldtube r Γ ; and also at fixed values of r > r Γ on each Cauchy slice, and then interpolated onto the black dots shown on the null slices. Values of ψ 4 at the black dots on a given null slice are then with A the area of the coordinate 2-surface t = constant, s = constant. Because the spacetime is dynamic, it would be a complicated process to construct t * exactly. Instead, it is assumed that extraction is performed in a region of spacetime in which the geometry is approximately Schwarzschild with (t, s) approximately standard Schwarzschild coordinates. Then In Eq. (108), M is an estimate of the initial mass of the system, usually the ADM mass, and g tt , g ss are averaged over the 2-sphere t = constant, s j = constant. It is straightforward to check that if (t, s) are exactly Schwarzschild coordinates, then t * is null.
In this way we obtain, for fixed t * , ψ m 4 at a number of different extraction radii; that is, we have data of the form ψ m 4 (t * , r k ), k = 1, . . . , K . 5 In practice, the real and imaginary parts of ψ m 4 may vary rapidly, and it has been found to be smoother to fit the data to the amplitude and phase. 6 For each spherical harmonic component two data-fitting problems are solved using the least-squares method, and where the A n , φ n are all real. Note that φ(t * , r ) must be continuous in r , and so for certain values of r k it may be necessary to add ±2π to arg(ψ m 4 (t * , r k )). Then the estimate is The remaining issue is the specification of N , and of the extraction spheres r k (or more precisely, since extraction is performed on spheres of specified coordinate radius, of the s k ). The key factors are the innermost and outermost extraction spheres, i.e., the values of r 1 and r K , and of course the requirement that K > N + 1. Essentially, the extrapolation process uses data over the interval 1/r ∈ [1/r K , 1/r 1 ] to construct an estimate at 1/r = 0. Polynomial extrapolation can be unreliable, or even divergent, as N is increased; it can also be unreliable when the distance from the closest data point is larger than the size of the interval over which the data is fitted. As a result of this latter condition, it is normal to require r K > 2r 1 . On the other hand, increasing r K increases the computational cost of a simulation, and decreasing r 1 could mean that a higher order polynomial is needed for accurate modelling of the data at that point. A compromise is needed between these conflicting factors. Values commonly used are that N is between 3 and 5, r 1 is normally of order 100 M, and r K is 300 M with values as large as 1000 M reported. Typically, K is about 8, with the 1/r k evenly distributed over the interval [1/r K , 1/r 1 ].
If the desired output of a computation is a waveform (to be used, say, in the analysis of LIGO detector data), then ψ 4 needs to be translated into its wave strain components where the constants of integration A m , B m need to be fixed by the imposition of some physical condition, for example that the strain should tend to zero towards the end of the computation. While this procedure is simple and straightforward, in practice it has been observed that the double time integration may lead to a reduction in accuracy, and in particular may introduce nonlinear drifts into the waveform. (The presence of a linear drift is easily corrected by means of an adjustment to the integration constants A m , B m ). It was shown in Reisswig and Pollney (2011) that the cause of the problem is that ψ m 4 includes random noise, and this can lead to noticeable drifts after a double integration. The usual procedure to control the effect is via a transform to the Fourier domain. The process to construct the wave strain from ψ m 4 , without any correction for drift, isψ where F is the Fourier transform operator, ω denotes frequency in the Fourier domain, and˜denotes a Fourier transformed function. The division by ω 2 in the second line of Eq. (114) is clearly potentially problematic for small ω, and an obvious strategy is to apply a filter to modify this equation. A number of such filters have been proposed, based on reducing those frequency components that are lower than ω 0 -the lowest frequency expected, on physical grounds, in the waveform. The simplest choice is a step function (Campanelli et al. 2009;Aylott 2009), but it has the drawback that it leads to Gibbs phenomena. To suppress this effect, a smooth transition is needed near ω 0 and various filters have been investigated (Santamaría et al. 2010;McKechan et al. 2010;Reisswig and Pollney 2011). A particularly simple choice of filter, yet effective in many cases , is

Energy, momentum and angular momentum in the waves
Starting from the mass loss result of Bondi et al. (1962), the theory of energy and momentum radiated as gravitational waves was further developed in the 1960s (Penrose 1963(Penrose , 1965aTamburino and Winicour 1966;Winicour 1968;Isaacson 1968) and subsequently (Geroch 1977;Thorne 1980a;Geroch and Winicour 1981). Formulas for the radiated angular momentum were presented in Campanelli and Lousto (1999), Lousto and Zlochower (2007) based on earlier work by Winicour (1980); formulas were also obtained in Ruiz et al. (2007), Ruiz et al. (2008) using the Isaacson effective stress-energy tensor of gravitational waves (Isaacson 1968). The result is formulas that express the energy, momentum and angular momentum content of the gravitational radiation in terms of ψ 4 . Strictly, all the quantities should be evaluated in the limit as r → ∞ and using an appropriate null tetrad. The energy equation is The linear momentum equations are wherer i is a unit radial vector. If the angular coordinate system being used is spherical polars, then from Eq. (399)r i = (sin θ cos φ, sin θ sin φ, cos θ), whereas if the coordinates are stereographicr i would by given by Eqs. (405) and (406). The angular momentum equations are (119) where theĴ i are operators given, in spherical polar coordinates, bŷ In practice, the above formulas are rarely used directly, and instead ψ 4 is first decomposed into spin-weighted spherical harmonics using Eq. (105). Then the energy equation is The momentum flux leaving the system is where The angular momentum equations become where the symbol refers to the imaginary part and

Gravitational waves in the Cauchy-perturbative approach
Black-hole perturbation theory has been fundamental not only for understanding the stability and oscillations properties of black hole spacetimes (Regge and Wheeler 1957), but also as an essential tool for clarifying the dynamics that accompanies the process of black hole formation as a result of gravitational collapse (Price 1972a, b).
As one example among the many possible, the use of perturbation theory has led to the discovery that Schwarzschild black holes are characterised by decaying modes of oscillation that depend on the black hole mass only, i.e., the black hole quasinormal modes (Vishveshwara 1970b, a;Press 1971;Chandrasekhar and Detweiler 1975). Similarly, black-hole perturbation theory and the identification of a powerlaw decay in the late-time dynamics of generic black-hole perturbations has led to important theorems, such as the "no hair" theorem, underlining the basic black-hole property of removing all perturbations so that "all that can be radiated away is radiated away" (Price 1972a, b;Misner et al. 1973).
The foundations of non-spherical metric perturbations of Schwarzschild black holes date back to the work of Regge and Wheeler (1957), who first addressed the linear stability of the Schwarzschild solution. A number of investigations, both gauge-invariant and not, then followed in the 1970s, when many different approaches were proposed and some of the most important results about the physics of perturbed spherical and rotating black holes established (Price 1972a, b;Vishveshwara 1970b, a;Chandrasekhar and Detweiler 1975;Zerilli 1970a, b;Moncrief 1974;Cunningham et al. 1978Cunningham et al. , 1979Teukolsky 1972Teukolsky , 1973. Building on these studies, which defined most of the mathematical apparatus behind generic perturbations of black holes, a number of applications have been performed to study, for instance, the evolutions of perturbations on a collapsing background spacetime Sengupta 1979b, a, 1980;Karlovini 2002;Seidel et al. 1987Seidel et al. , 1988Seidel 1990Seidel , 1991. Furthermore, the gauge-invariant and coordinate independent formalism for perturbations of spherically symmetric spectimes developed in the 1970s by Sengupta (1979b, a, 1980), has been recently extended to higher-dimensional spacetimes with a maximally symmetric subspace in Kodama et al. (2000), Kodama and Ishibashi (2003), Ishibashi and Kodama (2003), Kodama and Ishibashi (2004), for the study of perturbations in brane-world models.
Also nowadays, when numerical relativity calculations allow to evolve the Einstein equations in the absence of symmetries and in fully nonlinear regimes, black hole perturbative techniques represent important tools. 7 Schwarzschild perturbation theory, for instance, has been useful in studying the late-time behaviour of the coalescence of compact binaries in a numerical simulation after the apparent horizon has formed (Price and Pullin 1994;Abrahams and Cook 1994;Abrahams et al. 1995). In addition, methods have been developed that match a fully numerical and threedimensional Cauchy solution of Einstein's equations on spacelike hypersurfaces with a perturbative solution in a region where the components of three-metric (or of the extrinsic curvature) can be treated as linear perturbations of a Schwarzschild black hole [this is usually referred to as the "Cauchy-Perturbative Matching"] Rupright et al. 1998;Camarda and Seidel 1999;Allen et al. 1998;Rezzolla et al. 1999a;Lousto et al. 2010;Nakano et al. 2015). This method, in turn, allows to "extract" the gravitational waves generated by the simulation, evolve them out to the wave-zone where they assume their asymptotic form, and ultimately provide outer boundary conditions for the numerical evolution.
This section intends to review the mathematical aspects of the metric perturbations of a Schwarzschild black hole, especially in its gauge-invariant formulations. Special care is paid to "filter" those technical details that may obscure the important results and provide the reader with a set of expressions that can be readily used for the calculation of the odd and even-parity perturbations of a Schwarzschild spacetime in the presence of generic matter-sources. Also, an effort is made to "steer" the reader through the numerous conventions and notations that have accompanied the development of the formalism over the years. Finally, as mentioned in the Introduction, a lot of the material presented here has already appeared in the Topical Review by Nagar and Rezzolla (2006).

Gauge-invariant metric perturbations
It is useful to recall that even if the coordinate system of the background spacetime has been fixed, the coordinate freedom of general relativity introduces a problem when linear perturbations are added. In particular, it is not possible to distinguish an infinitesimal "physical" perturbation from one produced as a result of an infinitesimal coordinate transformation (or gauge-transformation). This difficulty, however, can be removed either by explicitly fixing a gauge (see, e.g., Regge and Wheeler 1957;Price 1972a, b;Vishveshwara 1970b, a;Zerilli 1970a, b), or by introducing linearly gaugeinvariant perturbations (as initially suggested by Moncrief 1974 and subsequently adopted in several applications Cunningham et al. 1978Cunningham et al. , 1979Seidel et al. 1987Seidel et al. , 1988Seidel 1990Seidel , 1991.
More specifically, given a tensor field X and its infinitesimal perturbation δ X, an infinitesimal coordinate transformation where L ξ is the Lie derivative along ξ . We will then consider δ X to be gauge-invariant if and only if L ξ X = 0, i.e., if δ X = δ X. In particular, since gravitational waves are metric perturbations, we will consider the case that X is the background metric g 0 , and then metric perturbations are gauge invariant if and only if L ξ g 0 = 0.
Stated differently, the possibility of building gauge-invariant metric perturbations relies on the existence of symmetries of the background metric. In the case of a general spherically symmetric background spacetime (i.e., one allowing for a time dependence) and which has been decomposed in multipoles (see Sect. 4.2), the construction of gauge-invariant quantities is possible for multipoles of order ≥ 2 only (Gerlach and Sengupta 1979b, a;Martín-García and Gundlach 1999;Gundlach and Martín-García 2000). In practice, the advantage in the use of gauge-invariant quantities is that they are naturally related to scalar observables and, for what is relevant here, to the energy and momentum of gravitational waves. At the same time, this choice guarantees that possible gauge-dependent contributions are excluded by construction.
Of course, this procedure is possible if and only if the background metric has the proper symmetries under infinitesimal coordinates transformation; in turn, a gaugeinvariant formulation of the Einstein equations for the perturbations of a general spacetime is not possible. Nevertheless, since any asymptotically flat spacetime can in general be matched to a Schwarzschild one at sufficiently large distances, a gaugeinvariant formulation can be an effective tool to extract physical information about the gravitational waves generated in a numerically evolved, asymptotically flat spacetime Rupright et al. 1998;Camarda and Seidel 1999;Allen et al. 1998;Rezzolla et al. 1999a) (see also Sect. 5.6 for additional implementational details). The following section is dedicated to a review of the mathematical techniques to obtain gauge-invariant perturbations of a the Schwarzschild metric.

Multipolar expansion of metric perturbations
Given a spherically symmetric Schwarzschild solution with metric g 0 and line element where e 2a = e −2b = (1 − 2M/r ), we generically consider small non-spherical perturbations h μν such that the new perturbed metric is where |h μν |/|g 0 μν | 1. Although we have chosen to employ Schwarzschild coordinates to facilitate the comparison with much of the previous literature, this is not the only possible choice, nor the best one. Indeed, it is possible to formulate the perturbations equations independently of the choice of coordinates as discussed in Martel and Poisson (2005), or in horizon-penetrating coordinates when the perturbations are in vacuum (Sarbach and Tiglio 2001). Mostly to remain with the spirit of a review and because most of the results have historically been derived in these coordinates, we will hereafter continue to use Schwarzschild coordinates although the reader should bear in mind that this is not the optimal choice. Because the background manifold M is spatially spherically symmetric, it can be written as the product M = M 2 × S 2 , where M 2 is a Lorentzian 2-dimensional manifold of coordinates (t, r ) and S 2 is the 2-sphere of unit radius and coordinates (θ, φ). As a result, the perturbations can be split "ab initio" in a part confined to M 2 and in a part confined on the 2-sphere S 2 of metric γ . Exploiting this, we can expand the metric perturbations h in multipoles referred to as "odd" or "even-parity" according to their transformation properties under parity. In particular, are odd (or axial) multipoles those that transform as (−1) +1 , under a parity transformation (θ, φ) → (π −θ, π+φ), while are even (or polar) those multipoles that transform as (−1) . As a result, the metric perturbations can be written as m=− , and the upper indices (o) and (e) distinguish odd and even-parity objects, respectively. Adapting now a notation inspired by that of Gerlach and Sengupta Sengupta 1979b, a, 1980) and recently revived by Gundlach and Martín-García (Martín-García and Gundlach 1999;Martín-García 2000, 2001), we use lower-case indices a, b . . . = 0, 1 to label the coordinates of M 2 and upper-case indices C, D . . . = 2, 3 to label the coordinates of S 2 . Using this notation, the scalar spherical harmonics Y m are then simply defined as where ∇ C indicates the covariant derivative with respect to the metric γ := diag(1, sin 2 θ) of S 2 , and where It is now convenient to express the odd and even-parity metric functions in (133) in terms of tensor spherical harmonics. To do this we introduce the axial vector S m C defined as S m where C D is the volume form on S 2 as defined by the condition C D C E = γ E D and such that ∇ C AB = 0. In this way, each odd-parity metric function in (133) can be written as where h, h Proceeding in a similar manner, each even-parity metric function can be decomposed in tensor spherical harmonics as where 1 K , and G (with the indices , m omitted for clarity) are the coefficients of the even-parity expansion, are also functions of (t, r ) only.
Note that we have used the Regge-Wheeler set of tensor harmonics to decompose the even-parity part of the metric in multipoles (Regge and Wheeler 1957). Despite this being a popular choice in the literature, it is not the most convenient one since the tensor harmonics in this set are not linearly independent. An orthonormal set is instead given by the Zerilli-Mathews tensor harmonics (Zerilli 1970;Mathews 1962) and the transformation from one basis to the other is given by defining the tensor Z m C D confined on the 2-sphere S 2 (see also section "Vector and tensor spherical harmonics" in "Appendix 2") and then replacing in Eq. (139) the second covariant derivative of the spherical harmon- This transformation has to be taken into account, for instance, when developing gauge-invariant procedures for extracting the gravitational-wave content from numerically generated spacetimes which are "almost" Schwarzschild spacetimes (Abrahams et al. 1992;Anninos et al. 1993Anninos et al. , 1995aAbrahams and Price 1996b, a).
Besides vacuum tensor perturbations, the background Schwarzschild spacetime can be modified if non-vacuum tensor perturbations are present and have a nonzero mass-energy, but much smaller than that of the black hole. In this case, the generic stress-energy tensor t μν describing the matter-sources can be similarly decomposed in odd and even-parity parts that are naturally gauge-invariant since the background is the vacuum Schwarzschild spacetime and are given explicitely by for the odd-parity part and by for the even-parity one. Note that we have now used the Zerilli-Matthews set of harmonics for the expansion, that the ten coefficients L m a , L m , T m ab , T m a , T m 2 , T m 3 are gauge-invariant, and that explicit expressions for them will be presented in the following sections.
Let us now consider the Einstein field equations that, in the static vacuum background, take the simple form R where R 0 is the Ricci tensor built from the background metric g 0 . At first order in the perturbations, the field equations reduce to where R is now the Ricci tensor built from the metric perturbations h. Note that while a generic perturbation will be a mixture of odd and even-parity contributions, we will exploit the linearity of the approach to handle them separately and simplify the treatment. In the following two sections we will discuss the form the Einstein equations (145) assume in response to purely odd and even-parity perturbations over a Schwarzschild background. In particular, we will show how the three odd-parity coefficients of the expansion in harmonics of the metric, i.e., h (o) a , h, and the seven even-parity ones, i.e., H 0 , can be combined to give two gauge-invariant master equations, named respectively after Regge and Wheeler (1957) and Zerilli (1970), each of which is a wave-like equation in a scattering potential. 8 Although our attention is here focussed on the radiative degrees of freedom of the perturbations (i.e., those with ≥ 2) because of their obvious application to the modelling of sources of gravitational waves, a comment should be made also on lowerorder multipoles. In particular, it is worth remarking that the monopole component of the metric for a vacuum perturbation (i.e., with = 0) is only of even-parity type and represents a variation in the mass-parameter of the Schwarzschild solution. On the other hand, the dipole component of the even-parity metric for a vacuum perturbation (i.e., with = 1) is of pure-gauge type and it can be removed by means of a suitable gauge transformation (Zerilli 1970b). This is not the case for a dipolar odd-parity metric perturbation, which can instead be associated to the introduction of angular momentum onto the background metric.

Gauge-invariant odd-parity perturbations
Before discussing the derivation of the odd-parity equation, we should make a choice for the odd-parity master function. Unfortunately, this choice has not been unique over the years and different authors have made different choices, making comparisons between different approaches less straightforward. Here, we will make a choice which highlights the relation with the gravitational-wave amplitude measured by a distant observer and, in particular, we construct the gauge-invariant combination of multipoles (Gerlach and Sengupta 1979b;Martín-García and Gundlach 1999;Harada et al. 2003) where, we recall, ∇ a represents the covariant derivative with respect to the connection of the submanifold M 2 . If ab is the antisymmetric volume form on M 2 , then the function is gauge-invariant and will be our choice for the Regge-Wheeler master function (Gerlach and Sengupta 1979b;Martín-García and Gundlach 1999;Gundlach and Martín-García 2000;Harada et al. 2003). A slight variation of the master function (147) has been introduced by Cunningham et al. (1978) in terms of the function ψ := Λ (o) and this has been used so extensively in the literature (Seidel et al. 1987(Seidel et al. , 1988Seidel 1991) that it is now commonly referred to as the Cunningham-Price-Moncrief (CPM) convention. We partly follow this suggestion and introduce a new master function for the odd-parity perturbations defined as With the choice (148), the Einstein field equations (145) with odd-parity perturbations lead to the inhomogeneous "Regge-Wheeler" equation where is the "tortoise coordinate" (Misner et al. 1973) and V (o) is the odd-parity potential, defined as The right-hand side of Eq. (149) represents the generic odd-parity "matter-source" and is given by with the components of the odd-parity matter-source vector defined as and where dΩ = sin θ dφdθ is the surface element on the 2-sphere S 2 .
Another choice for the gauge-invariant odd-parity master variable is possible and indeed was originally proposed by Moncrief (1974). This function, which hereafter we will refer to as the odd-parity Moncrief function, is defined as where the first expression is coordinate independent (Martel and Poisson 2005), while the second one is specialized to Schwarzschild coordinates with h 2 = −2h (Moncrief 1974). In the Regge-Wheeler gauge, i.e., for h 2 = h = 0, the definition (154) coincides with the variable used by Regge and Wheeler (1957). Historically, the choice of (154) as master variable has been the most common in the literature to describe odd-parity perturbations of a Schwarzschild spacetime and we will refer to it as "Regge-Wheeler" (RW) convention. It should be noted that while (154) is a solution of the Regge-Wheeler equation, the corresponding source term differs from expression (152). A general expression of the source in the RW convention can be found in Martel and Poisson (2005), Martel (2004) together with its specification for a point-particle (see also Andrade and Price 1999;Tominaga et al. 1999;Ferrari and Kokkotas 2000). The two master functions Q (o) and Ψ (o) are intimately related through the variational formalism employed by Moncrief (1974), and through the explicit expression (Martel and Poisson 2005) Note that Eq. (155) highlights an important difference between the two master functions which is not just a dimensional one (i.e., Ψ (o) has the dimensions of a length, while Q (o) is dimensionless) and this will have consequences on the asymptotic expressions for the gravitational waveforms when these are expressed in one or in the other convention. A detailed discussion of this will be made in Sects. 5.1, 5.3 and 5.4.

Gauge-invariant even-parity perturbations
Also in the case of even-parity perturbations, it is possible to express the evolution of the even-parity perturbations in terms of a wave-like equation in a scattering potential [cf. the Regge-Wheeler equation (149)]. In particular, following Moncrief (1974), we define the gauge-invariant functions and where their linear combination is also a gauge-invariant function. Strictly related to expression (158) is the gaugeinvariant function most frequently used in the literature (Seidel 1990;Gundlach and Martín-García 2001;Martel 2004;Lousto and Price 1997;Ruoff 2001;Ruoff et al. 2001;Martel and Poisson 2002;Nagar et al. 2004;Poisson 2004) which is also the solution of the inhomogeneous even-parity master equation or "Zerilli" equation and, again, is a wave-like equation in the scattering Zerilli potential (Zerilli 1970a) (161) The even-parity matter-source has a rather extended expression given by Martel (2004), Nagar et al. (2004Nagar et al. ( , 2005 Note that the expressions of the even-parity vector and tensor spherical-harmonics for the matter-source needed in (162) can be obtained from the orthogonality properties of the harmonics and are where the angular functions W m (θ, φ) and X m (θ, φ) are defined as Regge and Wheeler (1957) and where the overbar stands for complex conjugation. We should note that even-parity functions can be found in the literature under different notations. A particularly common choice is that proposed in Moncrief (1974) for the even parity gauge-invariant master function Q (e) which is related to Zerilli function (159) simply as Q (e) = ΛΨ (e) , while other authors use instead a master function defined as Z := 2Ψ (e) (Martel 2004;Lousto and Price 1997). Another evenparity function can be introduced in terms of two new gauge-invariant metric functions k and χ are defined as Martín-García (2000, 2001) and such that In this case, the Zerilli function (159) can be equivalently defined as Gundlach and Martín-García (2001) Finally, the homogeneous odd and even-parity master equations (149) and (160) can be transformed into each other by means of differential operations (Chandrasekhar 1983), and that they are connected to the master equation that Bardeen and Press have derived via the Newman-Penrose formalism (Bardeen and Press 1973).

Numerical implementations of the Cauchy-perturbative approach
In the previous section we have reviewed the derivation of the equations describing the evolution of perturbations of nonrotating black holes induced, for instance, by a nonzero stress-energy tensor. These perturbations have been assumed to be generic in nature, needing to satisfy only the condition of having a mass-energy much smaller than that of the black hole. The solution of these equations with suitable initial conditions completely specifies the reaction of the black hole to the perturbations and this is essentially represented by the emission of gravitational waves.
As mentioned in Sect. 4.1, the importance of the gauge-invariant variables used so far is that they are directly related to the amplitude and energy of the gravitational-wave signal measured at large distances. The purpose of this Chapter is to review the steps necessary to obtain the relations between the master functions for the odd and evenparity perturbations and the "plus" and "cross" polarisation amplitudes h + , h × of a gravitational wave in the TT gauge. In practice, and following the guidelines tracked in Cunningham et al. (1978Cunningham et al. ( , 1979, we will derive an expression for the perturbation metric h equivalent to that obtained in the standard TT gauge on a Minkowski spacetime and relate it to the odd and even-parity master functions Ψ (o) and Ψ (e) .
To obtain this result a number of conditions need to be met. First, we need to evaluate each multipole of the decomposed metric perturbations in the tetrad eν of stationary observers in the background Schwarzschild spacetime, i.e., hμν = e and where the indicesμ refer to the locally "flat" coordinates. Second, all of the quantities need to be evaluated far away from the source (i.e., in the "wave-zone") and in the so-called radiation gauge. In practice, this amounts to requiring that that components hθθ , hφφ and hθφ are functions of the type f (t − r )/r (i.e., they are outgoing spherical waves), while all the other components have a more rapid decay of O(1/r 2 ). Finally, we need to impose the condition that the metric is traceless modulo higher order terms, i.e., hθθ + hφφ = 0 + O(1/r 2 ).
In the following sections we will discuss the asymptotic expressions from odd-and even-parity perturbations, and how to implement the Cauchy-perturbative approach to extract gravitational-wave information within a standard numerical-relativity code.

Asymptotic expressions from odd-parity perturbations
We first consider odd-parity perturbations and recall that from the radiation-gauge conditions and since for large r the metric asymptotes that of a flat spacetime, i.e., e b ∼ e −b ∼ 1, we have Note that since h has the dimensions of a length squared, both h + and h × are dimensionless. Next, we need to relate the perturbation h to the odd-parity master function Ψ (o) . To do so, we follow the procedure outlined in Cunningham et al. (1978), and note that (cf. Eq. (III-20) 9 ) Equation (178)  Exploiting now the outgoing-wave behaviour of h at large distances we can write so that asymptotically Eq. (178) simply becomes and its integration yields 9 We recall that in the notation of Cunningham et al. (1978)  As a result, the "+" and "×" polarisation amplitudes of the gravitational wave can be calculated from Eqs. (177) Expressions (182) and (183) can be written in a compact form using the s = −2 spin-weighted spherical harmonics (see also section "Vector and tensor spherical harmonics" in "Appendix 2") so that expressions (182) and (183) can be combined into a single complex expression given by where, for clarity, we have explicitly restored the multipole indices , m.

The master function Q (o)
As discussed in Sect. 4.3, the odd-parity metric perturbations are sometimes expressed in terms of the odd-parity Moncrief function Q (o) [cf. Eq. (154)]; indeed it is not unusual to find in the literature the gravitational-wave amplitudes expressed in terms of this quantity. However, great care must be paid to the asymptotic relation between the master function Q (o) and the gravitational-wave amplitudes and, indeed, this is sometimes a source of confusion (Kawamura and Oohara 2004;Kawamura et al. 2003). To clarify this point, we recall that the derivation of the asymptotic relation between Q (o) and h proceeds in a way similar to the one discussed above. In the radiation gauge and at large distances from the black hole, we can use relation (179) in the definition (154) with h 2 = −2h, so that which is also a dimensionless quantity. Since h ∼ O(r ), the function Q (o) does not depend on r at leading order and Eq. (186) can be integrated to give where the integration constant can be defined as and it can be set to zero in the case of asymptotically flat metric perturbations (h = 0) at earlier times. Combining now expressions (167), (168) and (187), the gravitationalwave amplitudes in the two polarisations and with the new master function read (189) Note that in expressions (185) and (189) (149), but they yield two different asymptotic expressions for the gravitational-wave amplitudes. This difference, which is consistent with Eq. (155) when evaluated in a an asymptotic region of the spacetime where L m 1 = 0, is subtle but important and, as mentioned above, it has led to some inconsistencies in the literature both for the determination of the asymptotic gravitational-wave amplitudes and for the energy losses. This will be further discussed in Sects. 5.3 and 5.4.

Asymptotic expressions from even-parity perturbations
A calculation conceptually analogous to the one carried out in Sect. 5.1 leads to the relation between the gravitational-wave amplitude and the even-parity master function. In particular, after projecting Eq. (139) along the stationary tetrad, the asymptotic wave amplitudes in the two polarisation states are so that we essentially need to relate the metric perturbation G with the even-parity function Ψ (e) . Firstly, it is easy to realize that the even-parity metric projected onto the tetrad, h μν , is such that so that the terms proportional to these multipoles are of higher order for large r and can be neglected. Furthermore, from the transverse traceless condition we obtain an asymptotic relation between the gauge-invariant functions K and G where we have used the definition (135) to derive the right-hand side of expression (194). As a result, the asymptotic relation between the two components of the evenparity part of the perturbation metric simply reads Using now the definitions (156)- (157), we have that asymptotically and their linear combination (157) becomes Finally, the asymptotic gauge-invariant even-parity master function reads so that, modulo higher-order terms, the even-parity gravitational-wave amplitudes measured by a distant observer can be written in the compact form

Asymptotic general expressions
It is often convenient to combine the expressions for the asymptotic gravitational-wave amplitudes related to odd and even-parity perturbations into the single expression or, equivalently where we have defined h + := h × . Note that X 0 = 0 for any value of , so that in the case of axisymmetry the gravitational-wave signal is proportional to W 0 only.
It is also useful to underline that while expression (201) resembles the corresponding expression (10) of Kawamura and Oohara (2004), it is indeed different. Firstly, because in Kawamura and Oohara (2004) the Moncrief function is adopted for the odd-parity part of the perturbations and hence, modulo a normalisation factor, the function Ψ (o) appearing there corresponds to our function Q (o) [cf. expression (154)]. Secondly, because with this choice for the odd-parity perturbations a time derivative is needed in the asymptotic expression for the gravitational-wave amplitudes [cf. the discussion in the derivation of Eq. (189)]. As a result, expression (10) of Kawamura and Oohara (2004) (which is also missing the distinction between the real and imaginary parts) should really be replaced by expression (202). A similar use of the Moncrief function for the odd-parity part is present also in , Shibata and Sekiguchi (2003), Shibata and Sekiguchi (2005), where it is employed to calculate the gravitational-wave content of numerically simulated spacetimes.

Energy and angular momentum losses
Using the expressions derived in the previous sections we can now estimate the energy and angular momentum losses due to gravitational waves propagating outwards to spatial infinity. More specifically, this can be done by using expression (201) and the definition of Isaacson's stress-energy pseudo-tensor τ μν for the gravitational-wave field h propagating in the curved background field g 0 and in a Lorentz gauge (Isaacson 1968;Landau and Lifshitz 1975) where the brackets · · · refer to a spatial average over a length much larger than the typical gravitational wavelength (Isaacson 1968;Misner et al. 1973). The averaged expression (203) is gauge-invariant (Isaacson 1968) and holds in the "limit of high frequency" (or short-wave approximation), i.e., whenever the wavelength of the gravitational-wave field is small when compared to the local radius of curvature of the background geometry. In practice, gravitational radiation from isolated systems is of high frequency whenever it is far enough away from its source. Expression (203) accounts for the amount of energy and momentum carried by the gravitational wave over a certain region of spacetime, but since we are interested in the energy flux as measured by an inertial observer, we need to project the pseudo-tensor on the observer's locally orthonormal tetrad, where it becomes withhμν := hμν − 1 2 hημν and h being now the trace of hμν. As a result, the energy per unit time and angle carried by the gravitational waves and measured by a stationary observer at large distance is given by where the total derivative is made with respect to the asymptotic observer's time. Integrating (205) over the solid angle, the total power emitted in gravitational waves is then given by where expression (207) was first presented in Cunningham et al. (1978Cunningham et al. ( , 1979. Note that as discussed at the end of Sect. 5.3, these expressions need to be suitably modified when the energy losses are expressed in terms of the odd-parity Moncrief function Q (o) , in which case the energy-loss rate needs to be modified as Similarly, the angular momentum flux carried away in the form of gravitational waves can also be calculated in terms of the energy-momentum tensor (204). In particular, using spherical coordinates and assuming that the rotation is parametrised by the angle φ, we have Since the metric components in the radiation-gauge behave like outgoing spherical waves and since hθθ = h + and hθφ = h × , the angular momentum carried away in the form of gravitational waves (209) is then expressed as Proceeding in a way similar to the one followed in the calculation of the emitted power, the total angular momentum lost per unit time to gravitational wave reads (Martel and Poisson 2005) or, using the Moncrief master function (154) for the odd-parity perturbations (Poisson 2004) To conclude, we report the expression for the energy spectrum d E/dω, which is readily calculated from Eq. (206) after performing the Fourier transform of the odd and even-parity master functions, i.e., where we have indicated with f (ω, r ) the Fourier transform of the timeseries f (t, r ). Similarly, when using the odd-parity Moncrief function one obtains

A commonly used convention
A rather popular choice for the gauge-invariant master functions has found successful application in the extraction of the gravitational-wave content of numerically simulated spacetimes (Abrahams and Price 1996b, a;Abrahams et al. 1998;Rupright et al. 1998;Rezzolla et al. 1999a). For instance, the convention discussed below has been implemented in the Cactus computational toolkit (Camarda and Seidel 1999; Allen et al. 1998), a diffused and freely available infrastructure for the numerical solution of the Einstein equations (Allen et al. 1999;Cactus 2016). Numerous tests and applications of this implementation have been performed over the years and we refer the reader to Camarda and Seidel (1999), Allen et al. (1998), Font et al. (2002, Baiotti et al. (2005) for examples both in vacuum and non-vacuum spacetimes. The reference work for this convention is the one by Abrahams and Price (1996b, a), although a similar approach for the even-parity part of the perturbations was also adopted in previous works (Abrahams et al. 1992;Anninos et al. 1995b). We first note that the coefficients c 0 , c 1 and c 2 introduced in Abrahams and Price (1996b, a) are related simply to the multipolar coefficients of the odd-parity part introduced in Sect. 4.2. More specifically, considering that c 2 = −2h = h 2 , c 0 = h (o) 0 , and 1 , it is then easy to realise that the odd and even-parity master functions Q × m and Q + m defined in Abrahams and Price (1996b, a) are related to the master functions discussed so far through the simple algebraic expressions so that the asymptotic expression for the gravitational-wave amplitudes in the two polarisations are given by Similarly, expressions (217) and (218) can be combined into a single one which closely resembles expression (202) and that in its compactness highlights the advantage of the normalisation (215)-(216). We should remark that the notation in Eq. (219) could be misleading as it seems to suggest that h × is always of odd-parity and h + is always of even-parity. Indeed this is not true in general and in the absence of axisymmetry, i.e., when m = 0, both h × and h + are a superposition of odd and even parity modes. It is only for axisymmetric systems, for which only m = 0 modes are present, that Q × m and Q + m are real numbers, that h + is only even-parity and h × is only odd-parity.
Also very compact is the expression for the emitted power that, with this convention, simply reads Finally, the flux of linear momentum emitted in gravitational waves in the i-direction can be computed from the Isaacson's energy-momentum tensor and can be written in terms of the two polarization amplitudes as Favata et al. (2004) where n i = x i /r is the unit radial vector that points from the source to the observer. The calculation of this flux in terms of Q + m and Q × m can be computed after inserting Eq. (219) in Eq. (221), decomposing n i in spherical harmonics and performing the angular integral. This procedure goes along the lines discussed in Thorne (1980b), where all the relevant formulae are essentially available [cf. Eq. (4.20) there, but see also Sopuerta et al. (2006), Pollney et al. (2007)], so that we only need to adapt them to our notation. More specifically, in Thorne (1980b) the even-parity (or electric) multipoles are indicated with I m and the odd-parity (or magnetic) ones with S m , and are related to our notation by where ( ) where the asterisk indicates complex conjugation, we rewrite Eq. (4.20) of Thorne (1980b) in a more compact form. Following Damour and Gopakumar (2006) where the lowest multipolar contribution was explicitly computed in this way, it is convenient to combine the components of the linear momentum flux in the equatorial plane in a complex number as F x + iF y . The multipolar expansion of the flux vector can be written as Pollney et al. (2007) where δ m = 1 if m = 0 and δ m = 1/2 if m = 0. A more extended representation in terms of the various multipoles reads where a ± m := ( ± m)( ∓ m + 1), Note that here both F m x and F m y are real numbers and are obtained as the real and imaginary part of the right-hand side of Eq. (226).

Implementation summary
All of the material presented in the previous sections about the gauge-invariant description of the perturbations of a Schwarzschild black hole has laid the ground for the actual implementation of the Cauchy-perturbative extraction method in numerical-relativity calculations. We recall that the goal of the Cauchy-perturbative method is that of replacing, at least in parts of the three-dimensional numerical domain, the solution of the full nonlinear Einstein's equations with the solution of a set of simpler linear equations that can be integrated to high accuracy with minimal computational cost. In turn, this provides an unexpensive evolution of the radiative degrees of freedom, the extraction of the gravitational-wave information, and, if needed, the imposition of boundary conditions via the reconstruction of the relevant quantities at the edge of the three-dimensional computational domain.
In order to do this, it is necessary to determine the region of spacetime where a perturbative approach can be applied. In general, the three-dimensional numerical grid (indicated as N in Fig. 7) will comprise an isolated region of spacetime where the gravitational fields are strong and highly dynamical. In this region, indicated as S in Fig. 7, the full nonlinear Einstein equations must be solved. Outside of S, however, in what we will refer to as the perturbative region P, a perturbative approach is not only possible but highly advantageous. Anywhere in the portion of P covered by N we can place a two-dimensional surface, indicated as Γ in Fig. 7, which will serve as the surface joining numerically the highly dynamical strong-field region S and the perturbative one P. In practice, it is easier to choose this surface to be a 2-sphere of radius r Γ , where r Γ can either be the local coordinate radius, the corresponding Schwarzschild radial coordinate, or some more sophisticated radial coordinate deduced from the local values of the metric (cf. discussion in Sect. 4.2). 10 It is important to emphasize that the 2-sphere Γ need not be in a region of spacetime where the gravitational Fig. 7 Schematic illustration of a perturbative Cauchy extraction on a single spacelike hypersurface. Schematic picture of the Cauchy-perturbative matching procedure for a spacelike slice of spacetime (one spatial dimension has been suppressed). N is the three-dimensional numerical grid on the spacelike hypersurface Σ t in which the full Einstein equations are solved, and B its two-dimensional outer boundary. The interior (dark shaded) region S shows the strong-field highly dynamical region of spacetime fully covered by N. P is the region of spacetime where a perturbative solution can be performed and extends from the 2-sphere Γ (of radius r Γ ) to the 2-sphere A (of radius r A ) located in the asymptotically flat region of spacetime. P is covered entirely by a one-dimensional grid L and partially by the three-dimensional grid N fields are weak or the curvature is small. In contrast to approaches which matched Einstein's equations onto a Minkowski background Evans 1988, 1990), the matching is here made on a Schwarzschild background, so that the only requirement is that the spacetime outside of S approaches a Schwarzschild one. Of course, even in the case of a binary black-hole merger, it will be possible to find a region of spacetime, sufficiently distant from the black holes, where this requirement is met to the desired precision (Price and Pullin 1994;Abrahams and Cook 1994;Abrahams and Price 1996b, a;Abrahams et al. 1995).
In a practical implementation of the Cauchy-perturbative approach Rezzolla et al. 1999a), a numerical code provides the solution to the full nonlinear Einstein equations everywhere in the three-dimensional grid N except at its outer boundary surface B. At the extraction 2-sphere Γ , a different code (i.e., the perturbative module) "extracts" the gravitational wave information and transforms it into a set of multipole amplitudes which are chosen to depend only on the radial and time coordinates of the background Schwarzschild metric Rezzolla et al. 1999a).
In this way, two of the three spatial dimensions of the problem are suppressed and the propagation of gravitational waves on a curved background is reduced to a one-dimensional problem. During each timestep, information about the gravitational field read-off at Γ is propagated by the perturbative module out to the 2-sphere A in the asymptotic flat region of spacetime. This is done by solving a set of coupled one-dimensional linear differential equations (one for each of the multipoles extracted at Γ ) on the one-dimensional grid L covering the perturbative region P and ranging between r Γ and r A r Γ . From a computational point of view, this represents an enormous advantage: with a few straightforward transformations, the computationally Fig. 8 Schematic illustration of a perturbative Cauchy extraction. The Cauchy evolution is shown with green slices, comprising hypersurfaces Σ t on each of which is constructed a three-dimensional grid N The outer-boundary surface B of the three-dimensional grid is shown in dark blue, and is subject to a boundary condition that excludes incoming gravitational waves. Data from the Cauchy evolution on the worldtube Γ supplies boundary data to the perturbative equations, whose solution leads to the gravitational waves on the asymptotic boundary A. Note the difference between the asymptotic values of the gravitational waves extracted at A (filled blue circles) with the boundary values that can instead be injected on B expensive three-dimensional evolution of the gravitational waves via the nonlinear Einstein equations is replaced with a set of one-dimensional linear equations that can be integrated to high accuracy with minimal computational cost. Although linear, these equations account for all of the effects of wave propagation in a curved spacetime and, in particular, automatically incorporate the effects of backscatter off the curvature.
Note that as a result of this construction, (and as shown in Fig. 7), the perturbative region P is entirely covered by a one-dimensional grid L and only partially by a three-dimensional grid in the complement to S in N. The overlap between these two grids is essential. In fact, the knowledge of the solution on P allows the perturbative approach to provide boundary conditions at the outer boundary surface B and, if useful, Dirichlet data on every gridpoint of N outside the strong region S. This is also illustrated in Fig. 8, which represents a one-dimensional cut of Fig. 7, and highlights the difference between the asymptotic values of the gravitational waves extracted at the boundary A of the one-dimensional grid (filled blue circles) with and the boundary values that can be instead specified (i.e., "injected") on the outer boundary surface B of the three-dimensional grid.
The freedom to specify boundary data on a 2-surface of arbitrary shape as well as on a whole three-dimensional region of N represents an important advantage of the perturbative approach over similar approaches to the problem of gravitational-wave extraction and imposition of boundary conditions.
In what follows we briefly review the main steps necessary for the numerical implementation of the Cauchy-perturbative approach in a numerical-relativity code solving the Einstein equations in a 3+1 split of spacetime. This approach, which follows closely the discussion made in Rupright et al. (1998), Rezzolla et al. (1999a), basically consists of three steps: (1) extraction of the independent multipole amplitudes on Γ ; (2) evolution of the radial wave equations (247)-(249) on L out to the distant wave zone; (3) reconstruction of K i j and ∂ t K i j at specified gridpoints at the outer boundary of N. We next discuss in detail each of these steps.

Perturbative expansion
The first step is to linearize the Einstein equations around a static Schwarzschild background by separating the gravitational quantities of interest into background (denoted by a tilde) and perturbed parts: the three-metric γ i j = γ i j + h i j , the extrinsic curvature K i j = K i j + κ i j , the lapse N = N + α, and the shift vector β i = β i + v i . Note that the large majority of modern numerical-relativity codes implement the BSSNOK (Nakamura et al. 1987;Shibata and Nakamura 1995;Baumgarte and Shapiro 1999) or the CCZ4 (Alic et al. 2012) formulation of the Einstein equations. As mentioned in Sect. 3.2, in these formulations, the extrinsic curvature tensor is not evolved directly, but rather a traceless tensor extrinsic curvature tensor related to a conformal decomposition of the three-metric (Alcubierre 2008;Baumgarte and Shapiro 2010;Gourgoulhon 2012;Rezzolla and Zanotti 2013). Of course, also in these formulations it is possible to reconstruct the physically related extrinsic curvature tensor K i j and we will therefore continue to make use of K i j hereafter.
In Schwarzschild coordinates (t, r, θ, φ), the background quantities are given by g i j dx i dx j = N −2 dr 2 + r 2 (dθ 2 + sin 2 θ dφ 2 ), while the perturbed quantities have arbitrary angular dependence. The background quantities satisfy the dynamical equations ∂ t γ i j = 0, ∂ t N = 0, and thus remain constant for all time. The perturbed quantities, on the other hand, obey the following evolution equations where κ := κ i i and, as mentioned above, the tilde denotes a spatial quantity defined in terms of the background metric, γ i j . Note that the wave equation for κ i j involves only the background lapse and curvature.
Next, it is possible to simplify the evolution equation (237) by separating out the angular dependence, thus reducing it to a set of one-dimensional equations. This is accomplished by expanding the extrinsic curvature in Regge-Wheeler tensor spherical harmonics (Regge and Wheeler 1957) and substituting this expansion into (237). Using the notation of Moncrief (1974) we express the expansion as where (ê 1 ) i j , . . . , (f 4 ) i j are the Regge-Wheeler harmonics, which are functions of (θ, φ) and have suppressed angular indices ( , m) for each mode. Explicit expressions for these tensors are given in section "Regge-Wheeler harmonics" in "Appendix 2". The odd-parity multipoles (a × and b × ) and the even-parity multipoles (a + , b + , c + , and d + ) also have suppressed indices for each angular mode and there is an implicit sum over all modes in (238). The six multipole amplitudes correspond to the six components of κ i j . However, using the linearized momentum constraints we reduce the number of independent components of κ i j to three. An important relation is also obtained through the wave equation for κ, whose multipole expansion is simply given by κ = h(t, r )Y m . Using this expansion, in conjunction with the momentum constraints (239), we derive a set of radial constraint equations which relate the dependent amplitudes for each ( , m) mode.

Extraction
Taking the extraction 2-sphere Γ as the surface joining the evolution of the highly dynamical, strong field region (dark shaded area of Fig. 7) and the perturbative regions (light shaded areas), at each timestep, K i j and ∂ t K i j are computed on N as a solution to Einstein's equations. Assuming that N uses topologically Cartesian coordinates, 11 the Cartesian components of these tensors are then transformed into their equivalents in a spherical coordinate basis and their traces are computed using the inverse background metric, i.e., H = γ i j K i j , ∂ t H = γ i j ∂ t K i j . From the spherical components of K i j and ∂ t K i j , the independent multipole amplitudes for each ( , m) mode are then derived by an integration over the 2-sphere: Their time derivatives are computed similarly. Rather than performing the integrations (244)-(246) using spherical polar coordinates, it is useful to cover Γ with two stereographic coordinate "patches". These are uniformly spaced two-dimensional grids onto which the values of K i j and ∂ t K i j are interpolated using either a three-linear or a three-cubic polynomial interpolation scheme. As a result of this transformation, the integrals over the 2-sphere in (244)-(246) are computed avoiding polar singularities (see discussion in section "Stereographic coordinates" in "Appendix 2").

Perturbative evolution
Substituting (238) into (237) and using the constraint equations (240), we obtain a set of linearized radial wave equations for each independent amplitude. For each ( , m) mode we have one odd-parity equation and two coupled even-parity equations, These equations are related to the standard Regge-Wheeler and Zerilli equations (Regge and Wheeler 1957;Zerilli 1970a).
Once the multipole amplitudes, (a × ) m , (a + ) m , (h) m and their time derivatives are computed on Γ in the timeslice t = t 0 , they are imposed as inner boundary conditions on the one-dimensional grid. Using a suitably accurate integration scheme, the radial wave equations (247)-(249) can be evolved for each ( , m) mode forward to the next timeslice at t = t 1 . The outer boundary of the one-dimensional grid is always placed at a distance large enough that background field and near-zone effects are unimportant, and a radial Sommerfeld condition for the wave equations (247) can also be integrated using the data for K i j computed in this region. Note also that because h i j and α evolve along the coordinate time axis, these equations need only be integrated in the region in which their values are desired, not over the whole region L.
Of course, the initial data on L must be consistent with the initial data on N, and this can be determined by applying the aforementioned extraction procedure to the initial data set at each gridpoint of L in the region of overlap with N. In the latter case, initial data outside the overlap region can be set by considering the asymptotic fall-off of each variable.

Reconstruction
An important side product of the evolution step discussed above is that outer boundary values for N can now be computed, although, to the best of our knowledge, this procedure has not been implemented yet as a way to obtain outer boundary conditions. In particular, for codes using the BSSNOK (Nakamura et al. 1987;Shibata and Nakamura 1995;Baumgarte and Shapiro 1999) or the CCZ4 (Alic et al. 2012) formulation of the Einstein equations, it is sufficient to provide boundary data only for K i j , since the interior code can calculate γ i j at the outer boundary by integrating in time the boundary values for K i j .
In order to compute K i j at an outer boundary point of N (or any other point in the overlap between N and P), it is necessary to reconstruct K i j from the multipole amplitudes and tensor spherical harmonics. The Schwarzschild coordinate values (r, θ, φ) of the relevant gridpoint are first determined. Next, (a × ) m , (a + ) m , and (h m ) for each ( , m) mode are interpolated to the radial coordinate value of that point. The dependent multipole amplitudes (b × ) m , (b + ) m , (c + ) m , and (d + ) m are then computed using the constraint equations (240). Finally, the Regge-Wheeler tensor spherical harmonics (ê 1 ) i j -(f 4 ) i j are computed for the angular coordinates (θ, φ) for each ( , m) mode and the sum in Eq. (238) is performed. This leads to the reconstructed component of κ i j (and therefore K i j ). A completely analogous algorithm can be used to reconstruct ∂ t K i j in formulations in which this information is needed.
It is important to emphasize that this procedure allows one to compute K i j at any point of N which is covered by the perturbative region. As a result, the numerical module can reconstruct the values of K i j and ∂ t K i j on a 2-surface of arbitrary shape, or any collection of points outside of Γ .

Gravitational waves in the characteristic approach
The formalism for expressing Einstein's equations as an evolution system based on characteristic, or null-cone, coordinates is based on work originally due to Bondi (1960) and Bondi et al. (1962) for axisymmetry, and extended to the general case by Sachs (1962). The formalism is covered in the review by Winicour (2005), to which the reader is referred for an in-depth discussion of its development and the associated literature.
Most work on characteristic evolution uses, or is an adpatation of, a finite difference code that was originally developed at the University of Pittsburgh and has become known as the PITT null code. The early work that eventually led to the PITT code was for the case of axisymmetry (Isaacson et al. 1983;Bishop et al. 1990;Gómez et al. 1994), and a general vacuum code was developed in the mid-1990s (Bishop et al. 1996b(Bishop et al. , 1997bLehner 1998Lehner , 1999Lehner , 2001. Subsequently, the code was extended to the non-vacuum case (Bishop et al. 1999b, and code adaptations in terms of variables, coordinates and order of accuracy have been investigated (Gómez 2001;Gómez and Frittelli 2003;Reisswig et al. 2007Reisswig et al. , 2013a. Spectral, rather than finite difference, implementations have also been developed, for both the axially symmetric case (de Oliveira and Rodrigues 2009) and in general (Handmer and Szilágyi 2015). One potential difficulty, although in practice it has not been important in characteristic extraction, is the development of caustics during the evolution, and algorithms to handle the problem have been proposed (Stewart and Friedrich 1982;Corkill and Stewart 1983). There are also approaches that use outgoing null cones but for which the coordinates are not Bondi-Sachs (Bartnik 1997;Bartnik and Norton 2000).
Shortly after the publication of the Bondi and Bondi-Sachs metrics and formalism, the idea of conformal compactification was introduced. This led to the well-known asymptotic description of spacetime, and the definitions of asymptotic flatness, past, future and spacelike infinity (I + , I − , I 0 ), and of past and future null infinity (J − , J + ) (Penrose 1963); see also Penrose (1964Penrose ( , 1965b and Tamburino and Winicour (1966); and the reviews by Adamo et al. (2012) and Frauendiener (2004). The key result is that gravitational radiation can be defined unambiguously in an asymptotically flat spacetime only at null infinity. The waves may be expressed in terms of the Bondi news N (see Eq. (271) below), the Newman-Penrose quantity ψ 4 , or the wave strain (h + , h × ).
After a characteristic code has been run using a compactified radial coordinate as in Eq. (259), the metric is known at J + , and so it would seem to be straightforward to calculate the emitted gravitational radiation. Unfortunately, this is not in general the case because of gauge, or coordinate freedom, issues. The formulas do take a very simple form when expressed in terms of coordinates that satisfy the Bondi gauge condition in which the asymptotic flatness property is obviously satisfied, and for which conditions set at J + are propagated inwards along radial null geodesics. However, in a numerical simulation that is not the case: coordinate conditions are fixed on an extraction worldtube (in the case of characteristic extraction), or perhaps on a worldline (Siebel et al. 2003) or ingoing null hypersurface, and then propagated outwards to J + . The result is that the geometry at and near J + may appear very different to one that is foliated by spherical 2-surfaces of constant curvature. Of course, the Bondi gauge and the general gauge are related by a coordinate transformation, and formulas for N and ψ 4 are obtained by constructing the transformation.
An explicit formula in the general gauge for the news was obtained in Bishop et al. (1997b) ("Appendix 2"); and a calculation of ψ 4 was reported in Babiuc et al. (2009), but the formula produced was so lengthy that it was not published. These formulas have been used in the production of most waveforms calculated by characteristic codes. An alternative approach, in which the coordinate transformation is explicit, rather than partially implicit, was suggested (Bishop and Deshingkar 2003) but has not been further used or developed. Recently, a formula for the wave strain (h + , h × ), which is the quantity used in the construction of templates for gravitational-wave data analysis, was derived (Bishop and Reisswig 2014). An important special case is that of the linearized approximation, in which deviations from the Bondi gauge are small. The resulting formulas for N , ψ 4 and (h + , h × ), are much simpler and so much easier to interpret than in the general case. Further these formulas are widely used because the linearized approximation often applies to the results of a waveform computation in a realistic scenario.
We set the context for this section by summarizing the Einstein equations in characteristic coordinates, and outlining the characteristic evolution procedure. The focus of this section is formulas for gravitational waves, and we next present the formulas in the simplest case, when the coordinates satisfy the Bondi gauge conditions. Much of the remainder of the section will be devoted to formulas for gravitational waves in the general gauge, and will include a discussion of conformal compactification. This section makes extensive use of spin-weighted spherical harmonics and the eth formalism, which topics are discussed in "Appendix 2".

The Einstein equations in Bondi-Sachs coordinates
We start with coordinates based upon a family of outgoing null hypersurfaces. Let u label these hypersurfaces, φ A (A = 2, 3) be angular coordinates labelling the null rays, and r be a surface area coordinate. In the resulting x α = (u, r, φ A ) coordinates, the metric takes the Bondi-Sachs form where h AB h BC = δ A C and det(h AB ) = det(q AB ), with q AB a metric representing a unit 2-sphere; W c is a normalized variable used in the code, related to the usual Bondi-Sachs variable V by V = r + W c r 2 . It should be noted here that different references use various notations for what is here denoted as W c , and in particular (Bishop et al. 1997b) uses W with W := r 2 W c . As discussed in Sect. 1, we represent q AB by means of a complex dyad q A , then h AB can be represented by its dyad component J := h AB q A q B /2. We also introduce the fields K := 1 + JJ and U := U A q A . The spin-weight s of a quantity is defined and discussed in section "Spin-weighted fields" in "Appendix 2"; for the quantities used in the Bondi-Sachs metric We would like to emphasize two matters: (a) The metric Eq. (250) applies quite generally, and does not rely on the spacetime having any particular properties. (b) There are many different metrics of the form Eq. (250) that describe a given spacetime, and changing from one to another is known as a gauge transformation (about which more will be said later).
The form of the Einstein equations for the general Bondi-Sachs metric has been known for some time, but it was only in 1997 (Bishop et al. 1997b) that they were used for a numerical evolution. [see also Gómez and Frittelli 2003 for an alternative semifirst-order form that avoids second angular derivatives (ð 2 ,ð 2 ,ðð)]. The equations are rather lengthy, and only the hypersurface and evolution equations are given in that paper, in an "Appendix". 12 See also section "Computer algebra" in "Appendix 3". Here, in order to make the discussion of the Einstein equations precise but without being overwhelmed by detail, we give the equations in vacuum in the linearized case, that is when any second-order term in the quantities J, β, U, W c can be ignored. The Einstein equations are categorized into three classes, hypersurface, evolution, and constraint. The hypersurface equations are The evolution equation is The constraint equations are Reisswig et al. (2007) R 00 : 12 There is a misprint in Eq. (A3) of the journal version of the reference, which has been corrected in the version on the arXiv, and also in Reisswig et al. (2013a).

Fig. 9
Schematic illustration of the boundary data required for the characteristic code. The data required is J at u = 0 and β, J, U, ∂ r U, W c on the worldtube Γ An evolution problem is normally formulated in the region of spacetime between a timelike or null worldtube Γ and future null infinity (J + ), with (free) initial data J given on u = 0, and with boundary data for β, U, ∂ r U, W c , J satisfying the constraints given on Γ (Fig. 9). (In characteristic extraction, the data satisfies the Einstein equations inside Γ , and so the issue of ensuring that the boundary data must satisfy the characteristic constraint equations does not arise). The hypersurface equations are solved to find β, U, W c , and then the evolution equation gives ∂ u J and thence J on the "next" null cone. See Kreiss and Winicour (2011) and Babiuc et al. (2014) for a discussion of the well-posedness of the problem.
We extend the computational grid to J + by compactifying the radial coordinate r by means of a transformation r → x = f (r ) where lim r →∞ f (r ) is finite. In characteristic coordinates, the Einstein equations remain regular at J + under such a transformation. In practice, in numerical work the compactification is usually However, for the purpose of extracting gravitational waves, it is more convenient to express quantities as power series about J + , and so we compactify using (Common practice has been to use the notation for 1/r , but since we will have expressions involving the compactified radial coordinate and spherical harmonics such a notation would be confusing). Starting from the Bondi-Sachs metric Eq. (250), we make the coordinate transformation (260) to obtain In contravariant form, (262) Later, we will need to use the asymptotic Einstein equations, that is the Einstein equations keeping only the leading order terms when the limit r → ∞, or equivalently ρ → 0, is taken. We write the metric variables as J = J (0) + J (1) ρ, and similarly for β, U and W c . Each Einstein equation is expressed as a series in ρ and only leading order terms are considered. There is considerable redundancy, and instead of 10 independent relations we find (see section "Computer algebra" in "Appendix 3") The above are for the fully nonlinear case, with the linearized approximation obtained by setting K (0) = e −2β (0) = 1, and ignoring products of J and U terms.

The Bondi gauge
In the Bondi gauge, the form of the Bondi-Sachs metric is manifestly asymptotically flat since it tends to Minkowskian form as r → ∞. In order to see what conditions are thus imposed, the first step is to write the Minkowskian metric in compactified Bondi-Sachs coordinates. Starting from the Minkowski metric in spherical coordinates (t, r, φ A ), we make the coordinate transformation (t, r ) → (u, ρ) where to obtain We use the notation˜to denote quantities in the Bondi gauge. The metric of Eqs. (261) and (262) still applies, with the additional properties asρ → 0, The undifferentiated conditions can be regarded as defining the Bondi gauge, being motivated by the geometric condition that the metric Eq. (261) should take the form Eq. (268) in the limit as ρ → 0. The conditions on ∂ρβ, ∂ρŨ and ∂ρK follow from the asymptotic Eqs. (263), (264), and (423) respectively; and the condition on ∂ρW c is obtained from the asymptotic Einstein equation h AB R AB = 0 to second order in ρ and applying the Bondi gauge conditions already obtained. The null tetrad in the Bondi gauge will be denoted by˜ α ,ñ α [N P] ,m α , with components to leading order inρ [obtained by applying the coordinate transformation (267) to Eq. (96)] (270) The gravitational news was defined by Bondi et al. (1962) and is evaluated in the limitρ → 0, and is related to the strain in the TT gauge by where the rescaled strain H is which result is a straightforward consequence of the relationJ = h + + ih × discussed in section "Spin-weighted representation of deviations from spherical symmetry" in "Appendix 2". When using the Newman-Penrose formalism to describe gravitational waves, it is convenient to introduce since it will be important, when considering conformal compactification (Sect. 6.5), to have a quantity that is defined atρ = 0. In the Bondi gauge, as shown in section "Computer algebra" in "Appendix 3", ψ 0 4 simplifies to evaluated in the limitρ → 0. Thus ψ 0 4 , N and H are related by

General gauge
We construct quantities in the general gauge by means of a coordinate transformation to the Bondi gauge, although this transformation is largely implicit because it does not appear in many of the final formulas. The transformation is written as a series expansion in ρ with coefficients arbitrary functions of the other coordinates. Thus it is a general transformation, and the requirements that g αβ must be of Bondi-Sachs form, and thatg αβ must be in the Bondi gauge, impose conditions on the transformation coefficients. The transformation is ( Fig. 10) where the transformation coefficients u 0 , A u , ω, A ρ , φ A 0 , A A are all functions of u and φ A only. Conditions on the coefficients are found by applying the tensor transformation for specific cases of α, β, using the form of the metric in Eqs. (261) and (262) and also applying the conditions in Eq. (269) tog αβ andg μν (Bishop et al. 1997b;Bishop and Deshingkar 2003;Bishop and Reisswig 2014). The procedure is shown in some detail for one case, with the other cases being handled in a similar way. The actual calculations are performed by computer algebra as discussed in section "Computer algebra" in "Appendix 3". . Now consider a coordinate transformationθ → θ with dθ/dθ > 1 near whereθ is 0 and π , and with dθ/dθ < 1 near whereθ is π/2. In these coordinates, the surface of constant surface area r = 1 will be as shown in the right panel, and will not be a spherical surface of constant curvature From Eqs. (262) and (269),g 01 =ρ 2 + O(ρ 4 ). Then using the contravariant transformation in Eq. (278) with α = 0, β = 1, we havẽ Evaluating the right hand side to O(ρ 2 ), the resulting equation simplifies to give The remaining conditions follow in a similar way In the next equations, X 0 = q A φ A 0 , A = q A A A ; the introduction of these quantities is a convenience to reduce the number of terms in the formulas, since φ A 0 , A A do not transform as 2-vectors. As a result, the quantity ζ = q + i p also appears, and the formulas are specific to stereographic coordinates. We find so that to O(ρ 2 ) 0 = 2 Aω + 2 A u X 0 Uζ e −2β + K ðu 0 (2 +ðX 0 + 2X 0ζ ) so that at ρ = 0 and so that at ρ = 0 J = (1 + q 2 + p 2 ) 2 2(1 +q 2 +p 2 ) 2 ω 2 ðX 0 (2 + ðX 0 + 2X ζ ). (289) Explicit expressions for the null tetrad vectors n α [N P] and m α (but not α ) will be needed. n α [N P] is found by applying the coordinate transformation Eq. (278) toñ [N P] α (Eq. (270)) and then raising the index, giving to leading order in ρ The calculation of an expression for m α is indirect. Let F A be a dyad of the angular part of the general gauge metric, so it must satisfy Eq. (384), then Bishop et al. (1997b); Babiuc et al. (2009), with F A undetermined up to an arbitrary phase factor e −iδ(u,x A ) . We then define m α The suffix [G] is used to distinguish the above form from that defined in Eq. (270) since m α [G] = m α . However, it will be shown later (see Sects. 6.5.1, 6.5.2 and section "Computer algebra" in "Appendix 3") that the value of gravitational-wave descriptors is unaffected by the use of m α [G] rather than m α in its evaluation; thus it is permissible, for our purposes, to approximate m α by m α [G] . We now transform m α [G] in Eq. (292) to the Bondi gauge, The componentm 1 [G] is of the same order as˜ 1 , and som α [G] andm α are not equivalent. It can be checked (see section "Computer algebra" in "Appendix 3") (Bishop and Reisswig 2014) that the angular part ofm α [G] is equivalent tom α , sincem α [G]m α = 0 andm where |ν| = 1.
Since we actually require ν = 1, Eq. (294) can be used to set the phase factor δ explicitly. The result is An alternative approach (Bishop et al. 1997b;Babiuc et al. 2009), to the phase factor δ uses the condition that m α [G] is parallel propagated along J + in the direction n α [N P] , yielding the evolution equation where ∇ A is the covariant derivative with respect to the angular part of the metric h AB .

The gravitational-wave strain
An expression for the contravariant metricg αβ in the Bondi gauge is obtained from Eqs. (262) and (269), and each metric variable is expressed as a Taylor series about ρ = 0 (e.g.,J = 0 +ρ∂ρJ + O(ρ 2 )). Applying the coordinate transformation (277) we find g αβ , then useρ = ωρ + A ρ ρ 2 to express each component as a series in ρ; note that the coefficients are constructed from terms in the Bondi gauge, e.g., ∂ρJ . Then both sides of with g AB given in Eq. (262), are expressed as series in ρ, and the coefficients of ρ 1 are equated. This leads to an equation in which ∂ ρ J depends linearly on ∂ρJ (= H = limr →∞r (h + + ih × ), the rescaled strain defined in Eq. (273) (Bishop and Reisswig 2014), which may be inverted to give where the coefficients are These results are obtained using computer algebra as discussed in section "Computer algebra" in "Appendix 3". The above formula for the wave strain involves intermediate variables, and the procedure for calculating them is to solve equations for the variable indicated in the following order: Eq. (281) for x A 0 and thus X 0 , Eq. (287) for ω, Eq. (280) for u 0 , Eq. (283) for A u , and Eq. (285) for A.

Conformal compactification
Here we give only a brief introduction to this topic, as these matters are discussed more fully in many standard texts and reviews, e.g., Wald (1984) and Frauendiener (2004). We have made a coordinate compactification, resulting in the metric and null tetrad being singular at ρ = 0, which is therefore not included in the manifold. Thus, quantities are not evaluated at ρ = 0, but in the limit as ρ → 0. Introducing a conformal transformation has the advantage that this technical issue is avoided and J + at ρ = 0 is included in the manifold; but also that the resulting formulas for N and ψ 0 4 are simpler. (Of course, it should be possible to use the asymptotic Einstein equations to simplify expressions derived in physical space, but due to the complexity of the formulas this approach has not been adopted).
We use the notationˆfor quantities in conformal space. In the general gauge, the conformal metricĝ αβ is related to the metric Eq. (261) by g αβ = ρ −2ĝ αβ so that In a similar way, the Bondi gauge conformal metricĝ αβ is related to the Bondi gauge metricg αβ byg αβ =ρ −2ĝ αβ . Thusĝ αβ andĝ αβ are regular at ρ = 0 (or equivalently atρ = 0) and so in the conformal picture ρ = 0 is included in the manifold. The conformal metricsĝ αβ andĝ αβ are related bŷ which at ρ = 0 is the usual tensor transformation law with an additional factor ω −2 . A quantity that obeys this property is said to be conformally invariant with weight n where n is the power of ω in the additional factor; thus the metric tensor is conformally invariant with weight −2. In practice, it is not necessary to establish a relation of the form Eq. (303) to prove conformal invariance. The key step is to be able to show that a tensor quantity T a··· b··· satisfiesT a··· b··· = ρ −n T a··· b··· , then conformal invariance with weight n easily follows. In Eq. (303) the error term O(ρ) is shown explicitly, although it turns out to be irrelevant since the relation is evaluated at ρ = 0. This is generally the case, so from now on the error terms will not be taken into account; the one exception will be in the News calculation Sect. 6.5.1 which in places involves off-J + derivatives (since ∂ ρ O(ρ) = O(1)).
It is important to note that not all tensor quantities are conformally invariant, and in particular this applies to the metric connection and thus covariant derivativeŝ and to the Ricci scalar in n-dimensionŝ The Weyl tensor, however, is conformally invariant, so that the forms C α βγ δ and C αβγ δ are conformally invariant with weights 0 and −2 respectively. The construction of the conformal null tetrad vectors is not unique. It is necessary that orthonormality conditions analogous to Eq. (97) be satisfied, and it is desirable that the component of leading order in ρ should be finite but nonzero at ρ = 0. These conditions are achieved by defininĝ Thusn α [N P] andn a [N P] are related by the usual tensor transformation law, and With these definitions, n α [N P] , m α [G] are conformally invariant with weights 0 and 1 respectively.
Considering the conformally compactified metric of the spherical 2-surface described by the angular coordinates (φ A or φ A ) at J + , we have since ω =ρ/ρ. The curvature of J + is evaluated in two different ways, and the results are equated. The metric on the LHS is that of a unit sphere, and therefore has Ricci scalar R(x A ) = 2; and the metric on the RHS is evaluated using Eq. (305)with n = 2. Thus, leading to where the relationship between R(φ A ) and J, K is derived in Gómez et al. (1997), and where h AB ∇ A ∇ B log(ω) is given in terms of the ð operator in Eq. (B1) of Bishop et al. (1997b). Eq. (312) is a nonlinear elliptic equation, and in practice is not actually solved. However, it will be used later, when considering the linearized approximation, since in that case it has a simple analytic solution.

The news N
A difficulty with evaluating the gravitational radiation by means of Eq. (271) is that it is valid only in a specific coordinate system, so a more useful approach is to use the definition (Penrose 1963;Bishop et al. 1997b;Babiuc et al. 2009) At first sight Eqs. (271) and (313) do not appear to be equivalent, but the relationship follows by expanding out the covariant derivatives in Eq. (313) then expressing the metric coefficientsJ etc. as power series inρ as introduced just before Eq. (263). Using the Bondi gauge conditions Eq. (269), it quickly follows that −q A q BΓ 1 (0)AB /2 = ∂ũJ /2, which is zero, and the result follows since −q A q BΓ 1 (1)AB /2 = ∂ũ∂ρJ /2. Computer algebra has been used to check that replacingm α in Eq. (313) bym α [G] (withm α [G] defined in Eq. (293)) has no effect. 13 Because covariant derivatives are not conformally invariant, transforming Eq. (313) into the general gauge is a little tricky. We need to transform to physical space, where tensor quantities with no free indices are invariant across coordinate systems, and then to conformal space in the general gauge. From Eq. (304), and usingρ = ρω and Eq. (315). Using the conditions that (a) scalar quantities are invariant in physical space so thatm α it follows that Since m α [G] and∇ α ρ are orthogonal, we may write∇ α∇β (ρω) = ω∇ α∇β ρ +ρ∇ α∇β ω, so that This expression is simplified by (a) expanding out the covariant derivatives, (b) expressing the metric as a power series in ρ and using m α Finally, Eq. (292) is used to replace m α [G] in terms of F A , and the whole expression is divided byρ 3 , yielding (Bishop et al. 1997b;Babiuc et al. 2009) The limit is evaluated by expressing each metric coefficient as a power series in ρ, e.g., J = J (0) + ρ J (1) , and then writing . Direct evaluation combined with use of the asymptotic Einstein Eq. (266) shows that F A F BΓ 1 AB (0) = 0 (see section "Computer algebra" in "Appendix 3"), so that the limit evaluates to F A F BΓ 1 AB (1) . Further evaluation of Eq. (321) into computational ð form is handled by computer algebra, as discussed in section "Computer algebra" in "Appendix 3".
The attentive reader may have noticed that the derivation above usedρ = ρω rather thanρ = ρ(ω + ρ A ρ ), so that ∂ ρ ω should not be taken as 0 but as A ρ .

However, the corrections that would be introduced remain O(ρ) since (a)m 1
[G] = 0, (b) in Eq. (319) the termĝ 11 A ρ contained inĝ 1γ ∂ γ ω is O(ρ)A ρ , and (c) in Eq. (321) the term 6.5.2 The Newman-Penrose quantity ψ 0 4 The Newman-Penrose quantity ψ 4 , and its re-scaled version ψ 0 4 , were introduced in Sect. 3.3, and defined there for the case of physical space. Because the Weyl tensor is conformally invariant, it is straightforward to transform the earlier definition into one in the conformal gauge. Thus, in the conformal Bondi gauge, and again, as in the case of the news N , the limiting process means that the metric variables need to be expressed as power series inρ. Calculating the Weyl tensor is discussed in section "Computer algebra" in "Appendix 3", and the result is ψ 0 4 = ∂ 2 u ∂ρJ as given in Eq. (275). The "Appendix" also checks that replacingm α in Eq. (322) bym α [G] (withm α [G] defined in Eq. (293)) does not affect the result for ψ 0 4 . In this case, transforming Eq. (322) to the conformal general gauge is straightforward, since the tensor quantities are conformally invariant and the net weight is 0. The result is wherem α [G] = e −iδ (0, 0, F A ), and is further evaluated, by directly calculating the Weyl tensor, in section "Computer algebra" in "Appendix 3" (Babiuc et al. 2009) (but note that this reference uses a different approach to the evaluation of ψ 0 4 ).

Linearized case
In the linearized case the Bondi-Sachs metric variables β, J, U, W c and the coordinate transformation variables u 0 , A u , (ω − 1), A ρ , x A 0 , A A are regarded as small. Algebraically, the approximation is implemented by introducing a parameter = max(|β|, |J |, |U |, |W c |) in a neighbourhood of J + . Then, the metric variables are re-written as β → β etc., and quantities such as N , ψ 0 4 are expressed as Taylor series in with terms O( 2 ) ignored, leading to considerable simplifications. It is common practice to assume that the error in the approximation is about 2 . While computational results do not contradict this assumption, a word of caution is needed: no work on establishing a formal error bound for this problem has been reported.
Equations (280)-(282) and Eq. (285), simplify to It will also be useful to note the linearized form of Eq. (289), In the linearized case, Eq. (312) takes the form (Bishop 2005) 2 + 4(ω − 1) which is solved by decomposing ω and J into spherical harmonic components leading to [recall that Λ = ( + 1)] so that Evaluating Eq. (321) for the news, and using Eq. (325), is discussed in section "Computer algebra" in "Appendix 3". The result is Bishop (2005) Now from the linearized asymptotic Einstein equations, ∂ u J + ðU = 0 and ∂ ρ U − 2ðβ = 0, so we get The result is more convenient on decomposition into spherical harmonics, N = In the linearized case, the evaluation of ψ 0 4 is straightforward, because the Weyl tensor is a first-order term so the null tetrad vectors need be correct only to zeroth order. As discussed in section "Computer algebra" in "Appendix 3", we find (Babiuc et al. 2009) ψ 0

Numerical implementations of the characteristic approach
The idea of combining the "3+1" and characteristic approaches to extract the gravitational-wave signal from a numerical simulation was introduced in the 1990s (Bishop 1992(Bishop , 1993. However, this early work focused on using the combination for the whole spacetime and was called Cauchy-characteristic matching (CCM). Subsequently, implementation difficulties with CCM led to the development of something less ambitious known as characteristic extraction (CE). Under certain conditions (which in practice are achievable), CE is just as accurate as CCM. The advantage of CCM, if it can be implemented, is that it has the potential to make a significant contribution to overall code efficiency (Bishop et al. 1996b). An outline of what is meant by CCM and CE, and the differences between the two approaches, is illustrated in Fig. 11 and described in the caption.
As first steps towards CCM in relativity, it was implemented for the model problem of a nonlinear scalar wave equation (Bishop et al. 1996a(Bishop et al. , 1997a without any symmetries, and for the Einstein equations with a scalar field under the condition of spherical symmetry (Gómez et al. 1996;Lehner 2000). There has been a series of papers on CCM under axial symmetry (Clarke and d'Inverno 1994;Clarke et al. 1995;d'Inverno and Vickers 1997;d'Inverno et al. 2000;d'Inverno and Vickers 1996;Dubal et al. 1995Dubal et al. , 1998. A detailed algorithm for CCM in relativity in the general case was Fig. 11 Schematic illustration of Cauchy-characteristic matching (CCM) and characteristic extraction(CE). In both cases there is a Cauchy evolution green slices, and a characteristic evolution light blue slices between r Γ and J + . The difference is that in CCM the outer boundary of the Cauchy evolution is at the worldtube r Γ with boundary data supplied by the characteristic evolution; and in CE the outer boundary of the Cauchy evolution is as shown in dark blue and is subject to a boundary condition that excludes incoming gravitational waves presented in Bishop et al. (1999a). The stable implementation of matching is quite a challenge, and this goal has not yet been achieved Szilágyi 2000); although a stable implementation without symmetry has been reported with the Einstein equations linearized and using harmonic "3+1" coordinates (Szilágyi and Winicour 2003;Szilágyi et al. 2002). The issue of progress towards CCM is much more fully discussed in the review by Winicour (2005).
As a consequence of the difficulties with a stable implementation of CCM, in the 2000s attention shifted to the issue of developing CE, for which stability is not expected to be an issue. Further, although CCM has the advantage of high computational efficiency (Bishop et al. 1996b), it was realized that CE can be as accurate as CCM, provided the outer boundary of the Cauchy evolution is sufficiently far from the worldtube Γ that the two are not causally related, as indicated in Fig. 12. The implementation of CE for a test problem was described in 2005 (Babiuc et al. 2005). Subsequently, codes have been developed that yield useful results for the astrophysical problem of the inspiral and merger of two black holes Reisswig 2010;Reisswig et al. 2010;Babiuc et al. 2011a, b). Work that uses, rather than develops, characteristic extraction includes (Ott 2011;Reisswig et al. , 2013b. There is an alternative approach (Helfer 2010) that yields the emitted energy, momentum, and angular momentum, although it has not been implemented numerically.
The key feature of characteristic extraction is just a coordinate transformation, from Cauchy to characteristic coordinates in a neighbourhood of the worldtube Γ . However, it is a more complicated procedure than it might appear, because the Bondi-Sachs radial coordinate r is a surface area coordinate, so it cannot be expressed explicitly in terms of the "3+1" coordinates. This complicates the matter in two ways Fig. 12 Portions of the Kruskal diagram that are determined numerically. The green horizontal lines indicate the region of spacetime that is determined by the Cauchy evolution and which has finite spatial extent with artificial outer boundary at r B . The light blue diagonal lines indicate the region that is determined by characteristic evolution, and which starts off from a worldtube Γ located at r Γ using boundary data from the Cauchy evolution. The future Cauchy horizon of the Cauchy initial data is indicated by the dotted diagonal line L parallel to J + . As long as the worldtube Γ is located within the future Cauchy horizon, the numerically evolved subset of the spacetime is consistently determined 1. The coordinate transformation has to be made in two steps, firstly to a null coordinate system in which the radial coordinate is an affine parameter on the outgoing null radial geodesics, and secondly to Bondi-Sachs coordinates. 2. In general Γ is not a worldtube of constant r , so setting data at the innermost radial grid point of the Bondi-Sachs system requires special care.
Implementations of characteristic extraction mainly follow Bishop et al. (1999a), but do differ in certain aspects. Further, the field clearly needs to develop since most implementations are second-order accurate. This is a significant limitation since codes with higher order accuracy have existed for some time on the Cauchy side; and recently the first fourth order characterictic code has been reported (Reisswig et al. 2013a), as well as a spectral characteristic code (Handmer and Szilágyi 2015). An important recent development is the implementation of spectral characteristic extraction (Handmer et al. , 2016, so that the whole computation is spectrally convergent. Below we outline the characteristic extraction procedure including some of the variations currently in use.

Worldtube boundary data
Characteristic extraction is conceptually a post-processing procedure, as opposed to CCM in which the "3+1" and characteristic codes must run in step with each other. The transformation to Bondi-Sachs coordinates will require of the "3+1" data the full four-metric and all its first derivatives on the extraction worldtube Γ . Thus the first issue is to consider precisely what data the "3+1" code should dump to file, for subsequent processing by the characteristic extraction code. The simplistic solution of just dumping everything is not practical, because the data set is too large, and this is even the case should the data dump be restricted to those grid-points that are close to Γ . Some form of data compaction is required. On a given timeslice the extraction surface is spherical, so the natural compaction procedure is decomposition in terms of spherical harmonics. It turns out that this procedure also has some beneficial sideeffects 1. It filters out high frequency noise. 2. It greatly simplifies the process of interpolation onto a regular angular grid.
Assuming that the "3+1" spacelike coordinates are approximately Cartesian x i [C] = (x, y, z), the extraction worldtube Γ is defined by for some fixed radius R. Of course the above is a simple coordinate-specific, rather than a geometric, definition; but in practice the definition has worked well and Γ has not exhibited any pathologies such as becoming non-convex. As indicated above, the extraction code will need the full four-metric and its first-derivatives, and it is a matter of choice as to whether conversion from the "3+1" variables (such as lapse, shift and three-metric) to four-metric is performed in the "3+1" code or in the extraction routine; for simplicity, this discussion will be on the basis that the conversion is performed in the "3+1" code. The conversion formulas are given in Eq. (75). The time derivatives of the four-metric could be found by finite differencing, but the results are likely to be less noisy if they can be expressed in terms of other variables in the "3+1" codee.g., the "1+log" slicing condition and the hyperbolicΓ -driver condition , if being used, would mean that time derivatives of the lapse and shift are known directly, and the time derivative of the three-metric may be obtainable from the extrinsic curvature. For the spatial derivatives of the four-metric, it is sufficient to calculate and write to file only the radial derivative, since ∂ x , ∂ y , ∂ z can later be reconstructed from the radial derivative and angular derivatives of the spherical harmonics which are known analytically. The radial derivative in terms of the Cartesian derivatives is Having calculated the above variables and derivatives at "3+1" grid points in a neighbourhood of Γ , they must each be interpolated onto points on the coordinate sphere Γ using (at least) fourth-order interpolation. Then each quantity A, whether scalar, vector or tensor, is decomposed as and the A m are written to file. The decomposition is performed for ≤ max , and in practice max ≈ 8. A variation of the above procedure was introduced by Babiuc et al. (2011a). Instead of calculating the radial derivative of a quantity A, the idea is to decompose A into a product of spherical harmonics and Chebyshev polynomials in r . More precisely, we consider a "thick" worldtube R 1 < R < R 2 and the idea is to seek coefficients A k, ,m such that we may write where U k is a Chebyshev polynomial of the second kind, and so that, within the thick worldtube, the argument of U k has the required range of −1 to 1. The coefficients A k, ,m are then determined by a least squares fit to the data at each "3+1" grid point within the thick worldtube. The decomposition is carried out for k ≤ k max , and in practice (Babiuc et al. 2011a) takes k max = 6. Thus, this procedure involves writing three times as much data to file compared to that of calculating the radial derivative, but is probably more accurate and has the flexibility of being able to reconstruct data, including radial derivatives, at any point within the thick worldtube. Babiuc et al. (2011a) also introduced the option of calculating time derivatives via a Fourier transform process, so being able to filter out high frequency noise.

Reconstruction from spectral modes
The variables are reconstructed via in the case of decomposition only into angular modes, or via Eq. (342) in the case of decomposition into both angular and radial modes. The radial derivatives at the extraction worldtube R = R Γ are obtained either directly, or by analytic differentiation in the case that the reconstruction is in terms of Chebyshev polynomials. We then need to obtain the Cartesian derivatives in terms of radial and angular derivatives, and by the chain rule where φ A = φ 2 , φ 3 are the angular coordinates. The angular derivatives of the Y m may be re-expresed in terms of spin-weighted spherical harmonics, and ∂ i φ A expressed explicitly in terms of the Cartesian coordinates. The details depend on the specific angular coordinates being used, and in the common case of stereographic coordinates the formulas are Reisswig et al. (2010) where where the upper sign is valid for the north patch and the lower sign is valid for the south patch.

Transformation to null affine coordinates
In this section we construct the coordinate transformation from the Cartesian like "3+1" coordinates to a null coordinate system in which the radial coordinate is an Fig. 13 Schematic illustration of the (first stage) construction of characteristic coordinates and metric affine parameter rather than the Bondi-Sachs surface area coordinate. As already mentioned, we need this intermediate step because the surface area coordinate cannot be expressed as a function of only the "3+1" coordinates, but would also need terms involving the three-metric γ i j . The procedure is illustrated schematically in Fig. 13.
On the extraction worldtube Γ we can simply define the coordinate transformation, but off Γ it will need to be calculated. The null affine coordinates are x α [N ] = (u, λ, q, p), and the relation to the "3+1" coordinates x α [C] = (t, x, y, z) on Γ is defined to be u = t, λ = 0, with q, p given by Eq. (406) for r = R Γ . Although the above is given in terms of stereographic angular coordinates (q, p), rather than general angular coordinates φ A , the formulas that follow will not be specific to stereographic coordinates.
The unit normal n μ to the hypersurface Σ t is determined from the lapse and shift as stated in Eq. (71). Let s α = (s i , 0) be the outward pointing unit normal to the section S t of the worldtube at time t n . By construction, s i lies in the slice Σ t , and is given by Eq. (94). The generators α of the outgoing null cone through S t are given on the worldtube by which is normalized so that α t α = −1, where t α = (1, 0, 0, 0) is the Cauchy evolution vector.
We may now build the coordinate transformation between the "3+1" Cartesian coordinates x α and the (null) affine coordinates y α . As already discussed, we need this in a neighborhood of the worldtube, not just on the worldtube. Along each outgoing null geodesic emerging from S t , angular and time coordinates are defined by setting their values to be constant along the rays, and equal to their values on the worldtube. Geometrically, the idea is that we define (u, q, p) to be constant on each null geodesic generator. Algebraically, given (u, λ, q, p), the "3+1" coordinates are given by where x (0) [C] α is given by Eq. (350), and where α is given by Eq. (351). This expression determines x α [C] (u, λ, q, p) to O(λ 2 ). Consequently, the calculation of any quantity off-Γ is restricted to be be second-order accurate. If higher order is required, we would need to take into account how the geodesic generators, i.e., the α , vary off-Γ , which would mean using information provided by the geodesic equation.
Then the metric g [N ]αβ in null affine coordinates x α [N ] = (u, λ, q, p) is expressed in terms of the "3+1" metric as The Jacobian of the coordinate transformation is now expressed as a series expansion in the parameter λ. We do not need the ∂ λ x μ [C] because the coordinate λ is an affine parameter of the null geodesics, which fixes the g [N ]λμ : with the numerical value of the last condition a consequence of the normalization condition t α α = −1. The relevant part of the coordinate transformation is then (u, q, p).  (u, q, p); since (0)μ is known on the worldtube analytically in the angular directions from the spherical harmonic decomposition, and analytically or on a regular (q, p, u) grid in the time direction, ∂ ξ (0)μ can easily be found by analytic differentiation or finite differencing.

Null affine metric
The λ-derivative of the Cauchy 4-metric at the worldtube can be expressed as so that the null affine metric takes the form In the above Eq., the (0) coefficients are and the λ derivative coefficients are and where the λ-derivatives of the Cauchy metric are evaluated as The contravariant null affine metric, g αβ [N ] , is also expressed as an expansion in λ, [N ] The coefficients are obtained from the conditions as well as the requirement that certain components are fixed (which follows from Eq. (354)) g λu Thus the contravariant null affine metric and its λ derivative are

Metric in Bondi-Sachs coordinates
We are now able to construct the surface area coordinate r (u, λ, φ A ) In order to make the coordinate transformation x α , we need expressions for ∂ λ r , ∂ A r and ∂ u r . From Eq. (365), where in which the ∂ A ∂ C x i [C] are evaluated analytically in terms of φ A . An expression for ∂ u r will be required later but only on the worldtube Γ , so that Eq. (358) may be used when simplifying ∂ u applied to Eq. (365); further on Γ , ∂ u = ∂ t , and by construction ∂ A x i [C] is independent of time. Thus, The metric g αβ [B] in Bondi-Sachs coordinates is obtained from the coordinate transformation Note that the spherical part of the metric is unchanged, i.e., g AB [B] = g AB [N ] , and only the components g 11 [B] , g 1A [B] and g 01 [B] on Γ need to be determined. From Eq. (363), As discussed in Sect. 6.1, the characteristic Einstein equations are not formulated directly in terms of the metric components, but in terms of quantities derived from the metric, specifically J, β, U and W c . Explicitly, the relations between these quantities and the contravariant metric components are 2r 2 , β = − 1 2 log(g 01 [B] ), U = g 1A [B] g 01 [B] q A , 7.6 Starting up the null code at the worldtube As already mentioned, a difficulty faced is that Eq. (372) gives metric quantities on the worldtube Γ , which is not in general a hypersurface at a constant value of the r -coordinate. The original method for tackling the problem makes use of a Taylor series in λ (Bishop et al. 1999a), and has been implemented in Szilágyi et al. (2000), Szilágyi (2000), Babiuc et al. (2005), Reisswig et al. (2009) and Reisswig et al. (2010).
Recently, a method that uses a special integration algorithm between the worldtube and the first characteristic grid-point, has been proposed and tested (Babiuc et al. 2011a, b). Both approaches are outlined below.

Taylor series method
The Taylor series method is based on writing, for some quantity A, where A represents J, β, U and W c . A needs to be written as a function of null affine coordinates, so that ∂ λ A can be evaluated. Also, using Eqs. (365) and (367) evaluated on the worldtube, we need to find the value of λ at which r (λ) = r i , where r i is an r -gridpoint near the worldtube; this needs to be done for each grid-point in both the angular and time domains. The derivation of the Taylor expansions is straightforward, with second λ-derivatives eliminated using the Einstein equations (Bishop et al. 1999a). The results are

Special evolution routine between the worldtube and the first radial grid-point
In this approach, on a null cone say u = u n , we need only the values of the Bondi-Sachs metric variables at the angular grid-points on the worldtube. We also suppose that the value of J is known at all grid-points of the Bondi-Sachs coordinate system on the given null cone, either as initial data or from evolution from the previous null cone. A mask is set to identify those radial grid-points for which x i − x Γ < Δx, and these points will be called "B points". The special algorithm is concerned with setting data at the points i = B + 1, called "B+1 points". The first hypersurface equation in the hierarchy is the one for β, and is the simplest one to hadle. The algorithm is where Δ r = r B+1 − r Γ . The local truncation error associated with this algorithm is O(Δ 3 r ). The remaining hypersurface equations involve angular derivatives, which cannot be evaluated on the worldtube because it is not, in general, a hypersurface of constant r . Consequently, the right hand sides of these equations are evaluated at the B+1 points rather than at the points mid-way between r B+1 and r Γ . Schematically, the hypersurface equations are of the form (r n A) ,r = f , and the algorithm is The result is that the local truncation error for these equations is reduced to O(Δ 2 r ). Even so, one start-up step with error O(Δ 2 r ) is consistent with the global error of O(Δ 2 x ). Since the value of r varies on the worldtube, it may happen that the angular neighbour of a B+1 point is a B point. Thus, the code must also set data for the metric variables at the B points, even though much of this data will not be needed.

Initial data
The above discussion has shown how data should be set at, or on a neighbourhood of, the inner worldtube Γ , but in order to run a characteristic code data for J is also required on an initial null cone u = constant. Earlier work has adopted the simplistic but unphysical approach of just setting J = 0, assuming that the error so introduced would quickly be eliminated from the system. Babiuc et al. (2011a) and Bishop et al. (2011) investigated the matter. It was found that the error due to simplistic initial data is usually small, but it can take a surprisingly long time, up to 800 M, until saturation by other effects occurs. In terms of observations by a gravitational-wave detector, the effect of the error in search templates is not relevant. However, if a signal is detected, the effect would be relevant for accurate parameter estimation at large SNR (signal to noise ratio), but no quantitative estimates have been given.
Two methods for setting physically realistic initial data for a characteristic evolution have been proposed and tested. In Bishop et al. (2011) the initial data is set by means of fitting the boundary data to a general form of a linearized solution to the vacuum Einstein equations. On the other hand, Babiuc et al. (2011a) sets the initial data by means of the simple condition as in this case there should be no incoming radiation since the Newman-Penrose quantity ψ 0 = 0.

Implementation summary
The issues summarized here are: (1) setting up a characteristic code that starts from the output of a "3+1" code; (2) estimating the gravitational waves from metric data in a compactified domain output by a characteristic code; (3) estimating quantities derived from the gravitational waves, i.e., the energy, momentum and angular momentum.

Setting worldtube boundary data for the characteristic code
The coding of characteristic extraction is a complex process, and is not simply a matter of implementing a few of the formulas derived earlier in this section. Below, we outline the key steps that are required. The reader is also referred to section "Numerical codes" in "Appendix 3" for information about computer code that implements characteristic extraction.
1. Within the "3+1" code, write a routine that uses Eq. (342) to perform a spectral decomposition of the three-metric, lapse and shift, and outputs the data to file. 2. In a front-end to the characteristic code, write a routine that reads the data from the file created in the previous step, and reconstructs the four-metric and its first derivatives at the angular grid-points of the extraction worldtube. 3. Construct the generators α of the outgoing null cone using Eq. (351), and then the Jacobian ∂ x μ [C] /∂ x α [N ] as a series expansion in the affine paramenter λ, for each angular grid-point on the worldtube. 4. As described in Sect. 7.4, construct the null affine metric g [N ]αβ and its first λderivative at the angular grid-points of the extraction worldtube; then construct the contravariant forms g αβ [N ] and ∂ λ g αβ [N ] . 5. From Eq. (365), determine the surface area coordinate r and its first derivatives at the angular grid-points of the extraction worldtube. 6. Construct the Jacobian ∂ x μ [B] /∂ x α [N ] , and thus the Bondi-Sachs metric g αβ [N ] and then the metric coefficients β, J, U, W c at the angular grid-points of the extraction worldtube. 7. Implement either of the special start-up procedures described in Sect. 7.6. 8. The construction of a characteristic code is not described in this review, but see section "Numerical codes" in "Appendix 3" for information about the availability of such codes.

Estimation of gravitational waves
In Sect. 6, many formulas used ρ (=1/r ) as the radial coordinate, but that is unlikely to apply in practice. In the case that the radial code coordinate is x given by Eq. (259), the relation between ∂ x and ∂ ρ at J + is If all the metric coefficients near J + are small (which, in practice, is often but not always the case), then the linearized formulas apply, and: ψ 0 4 is evaluated using Eq. (334). -The news N is evaluated, decomposed into spherical harmonics, using Eq. (332).
In the general (nonlinear) case, it is first necessary to evaluate the coordinate transformation functions φ A 0 (u, x A ), ω(u, x A ) and u 0 (u, x A ). The reason for doing so is that it is then possible to determine, at each (u, φ A ) grid point, the corresponding values of the Bondi gauge coordinates (ũ,φ A ). Thus the gravitational-wave quantities can be expressed as functions of the (physically meaningful) Bondi gauge coordinates, rather than as functions of the code coordinates. This issue did not arise in the linearized case because it is a second-order effect and thus ignorable. The procedure for evaluating these functions is: Solve the evolution problem Eq. (281) with initial data φ A 0 (0, x A ) = 0. This initial condition assumes that the initial data for J has been set with J = 0 at J + .
ω(u, x A ). Either solve the evolution problem Eq. (282) with initial data ω(0, x A ) = 1, or evaluate the explicit formula Eq. (287). u 0 (u, x A ). Solve the evolution problem Eq. (280). In this case, there is a gauge freedom to set the initial data u 0 (0, x A ) arbitrarily. In the cases of N and ψ 0 4 , the phase factor δ(u, x A ) also needs to be evaluated. This can be done either explicitly, Eq. (295), or by solving the evolution problem Eq. (296) with initial data δ(0, x A ) = 0. Then: ψ 0 4 is evaluated using Eq. (323). -The news N is evaluated using Eq. (321).

Energy, momentum and angular momentum in the waves
The formulas for the energy, momentum and angular momentum have already been given in terms of ψ 4 in Sect. 3.3.3, and these formulas are directly applicable here on substituting ψ 4 by ψ 0 4 /r . The resulting formulas involve one or two time integrals of ψ 0 4 , and it is useful to note that here all such integration can be avoided by using

A comparison among different methods
This review has described the following methods for extracting the gravitational-wave signal from a numerical simulation -The quadrupole formula, including various modifications, leading to the wave strain (h + , h × ); -ψ 4 (fixed radius) and ψ 4 (extrapolation), leading to the Newman-Penrose quantity ψ 4 ; -Gauge-invariant metric perturbations, leading to the wave strain (h + , h × ); -Characteristic extraction, leading to the wave strain (h + , h × ), the gravitational news N , or the Newman-Penrose quantity ψ 4 .
There are a number of factors that need to be taken into account in deciding the appropriate method for a particular simulation. In outline, these factors are: -Physical problem motivating the simulation. The most appropriate method for extracting gravitational waves is affected by how the result is to be used. It may be that only moderate accuracy is required, as would be the case for waveform template construction for use in searches in detector data; on the other hand, high accuracy would be needed for parameter estimation of an event in detector data at large SNR. Further, the purpose of the simulation may be not to determine a waveform, but to find the emitted momentum of the radiation and thus the recoil velocity of the remnant. -Domain and accuracy of the simulation. The domain of the simulation may restrict the extraction methods that can be used. All methods, except that using the quadrupole formula, require the existence of a worldtube, well removed from the domain boundary, on which the metric is Minkowskian (or Schwarzschild) plus a small correction. As discussed in Sect. 3.3.2, extrapolation methods need these worldtubes over an extended region. Further, the accuracy of the simulation in a neighbourhood of the extraction process clearly limits the accuracy that can be expected from any gravitational-wave extraction method. -Ease of implementation of the various extraction methods. All the methods described in this review are well understood and have been applied in different contexts and by different groups. Nevertheless, the implementation of a new gravitational-wave extraction tool will always require some effort, depending on the method, for coding, testing and verification. -Accuracy of the various extraction methods. Theoretical estimates of the expected accuracy of each method are known, but precise data on actual performance is more limited because suitable exact solutions are not available. In a simulation of a realistic astrophysical scenario, at least part of the evolution is highly nonlinear, and the emitted gravitational waves are oscillatory and of varying amplitude and frequency. On the other hand, exact solutions are known in the linearized case with constant amplitude and frequency, or in the general case under unphysical conditions (planar or cylindrical symmetry, or non-vacuum). One exception is the Robinson-Trautman solution (Robinson and Trautman 1962), but in that case the gravitational waves are not oscillatory and instead decay exponentially.
Thus, in an astrophysical application, the accuracy of a computed waveform is estimated by repeating the simulation using a different method; then the difference between the two waveforms is an estimate of the error, provided that it is in line with the theoretical error estimates. In some work, the purpose of comparing results of different methods is not method testing, but rather to provide validation of the gravitational-wave signal prediction. The only method that is, in principle, free of any systematic error is characteristic extraction, but the method was not available for general purpose use until the early 2010s. It should also be noted that there remains some uncertainty about factors that could influence the reliability of a computed waveform (Boyle 2016).

Comparisons of the accuracy of extraction methods
A number of computational tests have been reported, in which the accuracy of various extraction methods is compared. Such tests are, of course, specific to a particular physical scenario, and to the choice of "3+1" evolution code, initial data, gauge conditions, etc. Some of the tests reported are now outlined, together with the results that were obtained. While it is natural to want to generalize these results, a word of caution is needed since the testing that has been undertaken is quite limited. Thus any generalization should be regarded as providing only a guide to which there may well be exceptions.
- Nagar et al. (2005) investigates various modifications of the standard quadrupole formula in comparison to results obtained using gauge-invariant metric perturbations for the case of oscillating accretion tori. Good results are obtained when back-scattering is negligible, otherwise noticeable differences in amplitude occur. - Balakrishna et al. (2006) computes ψ 4 (fixed radius) and gauge-invariant metric perturbations for gravitational waves from boson star perturbations, but detailed comparisons between the two methods were not made. - Pollney et al. (2007) compares gauge-invariant metric perturbations to ψ 4 (fixed radius) extraction for the recoil resulting from a binary black hole merger. It was found that results for the recoil velocities are consistent between the two extraction methods. -  and Baiotti et al. (2009) compare gauge-invariant metric perturbations, modified quadrupole formula and ψ 4 (fixed radius) extraction for a perturbed neutron star. While the results are generally consistent, each method experienced some drawback. The gauge-invariant method has a spurious initial junk component that gets larger as the worldtube radius is increased. In ψ 4 extraction, fixing the constants of integration that arise in obtaining the wave strain can be a delicate issue, although such problems did not arise in this case. The generalized quadrupole formula led to good predictions of the phase, but to noticeable error in the signal amplitude. - Reisswig et al. (2009Reisswig et al. ( , 2010 compare ψ 0 4 from characteristic extraction and from ψ 4 -extrapolation for Binary Black Hole (BBH) inspiral and merger in spinning and non-spinning equal mass cases. The "3+1" evolution was performed using a finite difference BSSNOK code . A comparison was also made in Babiuc et al. (2011a) for the equal mass, non-spinning case. Recently, a more detailed investigation of the same problem and covering a somewhat wider range of BBH parameter space, was undertaken (Taylor et al. 2013) using SpEC for "3+1" evolution ). These results lead to two main conclusions. (1) The improved accuracy of characteristic extraction is not necessary in the context of constructing waveform templates to be used for event searches in detector data. (2) Characteristic extraction does provide improved accuracy over methods that extract at only one radius. The ψ 4 -extrapolation method performs better, and there are results for ψ 0 4 that are equivalent to characteristic extraction in the sense that the difference between the two methods is less than an estimate of other errors. However, that does not apply to all modes, particularly the slowly varying m = 0 "memory" modes.
-A study of gravitational-wave extraction methods in the case of stellar core collapse ) compared characteristic extraction, ψ 4 extraction (fixed radius), gauge-invariant metric perturbations, and the quadrupole formula. In these scenarios, the quadrupole formula performed surprisingly well, and gave results for the phase equivalent to those obtained by characteristic extraction, with a small under-estimate of the amplitude. However, quadrupole formula methods fail if a black hole forms and the region inside the horizon is excised from the spacetime. The gauge-invariant metric perturbation method gave the poorest results, with spurious high frequency components introduced to the signal. In characteristic extraction and ψ 4 extraction the waveform was obtained via a double time integration, and the signal was cleaned up using Fourier methods to remove spurious low frequency components. -It is only very recently (Bishop and Reisswig 2014) that a method was developed in characteristic extraction to obtain the wave strain directly instead of via integration of N or ψ 0 4 . That work also compared the accuracy of the waveform obtained to that found from integration of ψ 0 4 using ψ 4 -extrapolation, in two cases-a binary black hole merger, and a stellar core collapse simulation. When comparing the wave strain from characteristic extraction to that found by time integration of ψ 0 4 , good agreement was found for the dominant (2,2) mode, but there were differences for ≥ 4.

Appendix 2: Spin-weighed spherical harmonics and the ð operator
A convenient way to represent vector and tensor quantities over the sphere, including their angular derivatives, is to use spin-weighted quantities and the ð operator. The formalism was introduced by Newman and collaborators in the 1960s (Newman and Penrose 1966;Goldberg et al. 1967), and has been described in textbooks such as Rindler (1984, 1986) and Stewart (1990). Even so, the theory is not well known and there are variations in notation and conventions (which topic is discussed further in Sect. 1), so the theory will be presented here in some detail based on the conventions of Gómez et al. (1997). We will describe the theory, both in general terms and with specific reference to the coordinates commonly used, i.e., spherical polar and stereographic. Further, the spin-weighted spherical harmonic functions will be introduced, as well as the vector and tensor spherical harmonics (Newman and Silva-Ortigoza 2006). All of these can be used as basis functions on the sphere.

The complex dyad
Let φ A , q AB and q AB (2 ≤ A, B ≤ 3) be coordinates and the associated metric of a unit sphere. For example, in standard spherical polars, φ A = (θ, φ) and The first step is to introduce a complex dyad q A . Geometrically, the dyad is a 2vector that can be written as q A = (q A ) + i (q A ) where (q A ), (q A ) are real and orthonormal. In other words, the real and imaginary parts of q A are unit vectors at right-angles to each other. From this definition, it is straightforward to verify the following properties Clearly, the dyad q A is not unique being arbitary up to a rotation and/or reflection, so that p A = e iγ q A and p A =q A are also dyads. Even so, it is convenient to write the dyad in a way that is natural to the coordinates being used, which for a diagonal metric means that the real part of the dyad should be aligned to the φ 2 direction, and the imaginary part to the φ 3 direction. Thus, the dyad usually used in spherical polar coordinates is q A = 1, i sin θ , and q A = (1, i sin θ).
The possibility of different parities, i.e., of a reflection, introduces an additional complication, which can be avoided by convention: All the coordinate systems used in this review are right-handed, and we will not consider dyads that are related by complex conjugation. This is equivalent to embedding the sphere into Euclidean space with coordinates (r, φ 2 , φ 3 ) that are always right-handed, i.e., the vector product (q) × (q) points in the positive r -direction.
The dyad q A resembles the angular components of the tetrad vector m α given in Eq. (96) when evaluated on a unit sphere, but there is a difference of a factor of √ 2. This is an example of different conventions used by different authors, which matter is discussed further in Sect. 1.

Spin-weighted fields
Having defined the dyad q A , we are now in a position to construct spin-weighted fields from vector and tensor fields. In the simplest case, suppose that μ(φ 2 , φ 3 ) is a scalar field. Then μ is also a spin-weighted field with spin-weight s = 0. Given a vector field v A (φ 2 , φ 3 ), we define to be a field with spin-weight s = 1. Since V is a complex quantity, it contains two independent fields and thus uniquely represents the two components of v A . We can also define the quantityV =q A v A with spin-weight s = −1, but of course it is not independent of V . From a second-rank tensor t AB we can construct two independent fields with spin-weights s = 2 and s = 0 respectively. Together T and τ have 4 independent fields to represent the 4 components of t AB . In general, given a tensor field w A···D , the quantity W = w A···BC···D q A · · · q Bq C · · ·q D , with m factors q and n factorsq is defined to be a quantity with spin-weight s = m −n. Spin-weighted quantities may also be defined in terms of their transformation properties. However, here, we use the definition above and later derive the transformation rule Eq. (398).

Differentiation and the ð operator
We would now like to define derivative operators ð andð that act on spin-weighted fields, and that are consistent with covariant differentiation. This means that if W is defined as in Eq. (388), then we would like ðW = ∇ E w A···BC···D q A · · · q Bq C · · ·q D q E ,ðW = ∇ E w A···BC···D · · · q Bq C · · ·q Dq E .
(389) We see immediately that, if W has spin-weight s, then ðW has spin-weight s + 1 and ðW has spin-weight s − 1. We achieve the desired effect by the definition where (Normally, the quantity denoted here by ϒ is given the notation Γ ; but we use the notation ϒ because Γ is used to represent the extraction worldtube). We demonstrate the above for a spin-weight 1 field V = q A v A . Starting from the definition Eq. (389) and the ansatz Eq. (390) with ϒ to be determined, we obtain and after renaming some dummy indices. This must be true for all v C , and the bracketed term is just a covariant deriavtive, so we have We lower the free superscript C and then contract withq C to obtain the desired result.
It should be noted that, in general, ð andð do not commute. The commutator is so that the operators commute only in the case of a quantity with spin-weight s = 0.

Coordinate transformation of spin-weighted quantities: Rotation factors exp(iγ )
Spin-weighted quantities are defined in a way that they have no free tensorial indices so it would appear that they are scalars, but this is misleading because different dyads are used in the different coordinate systems. Suppose that we have two coordinate systems S (q) and S ( p) with natural dyads q and p respectively. Each dyad has components in each of the coordinate systems, and so we define q A (q) , q A ( p) to mean the components of q in S (q) , S ( p) respectively; and similarly for p A (q) and p A ( p) . Assuming that S (q) and S ( p) have the same parity, then their dyads are related by a rotation and, as discussed just after Eq. (384), p Suppose that v is a vector and that V (q) , V ( p) are the corresponding spin-weighted quantities with respect to the dyads q, p respectively. Thus Generalizing to the case where V is defined with m factors q A and m − s factorsq A so that the spin-weight is s, we find

Specific coordinate systems
It is very convenient to define coordinate systems on the unit sphere in terms of a coordinate transformation from a Cartesian system. The Cartesian coordinates are denoted by x a [C] = (x, y, z), and the spherical coordinates by x a = (r, φ 2 , φ 3 ). We use computer algebra to construct the Jacobian of the transformation bewteen spherical and Cartesian coordinates, and then to find the components of the dyad with respect to the Cartesian system. We use the notation Q i to denote the components of a dyad with respect to Cartesian coordinates. The computer algebra also evaluates the quantity ϒ used in the definition of the ð operator Eq. (391), see section "Computer algebra" in "Appendix 3".

Spherical polar coordinates
We use coordinates (r, θ, φ) for standard spherical polar coordinates. Only one coordinate patch is required, but the coordinate system is singular at the poles θ = 0 and θ = π . The relation to Cartesian coordinates is x = r sin θ cos φ, y = r sin θ sin φ, z = r cos θ and the inverse transformation is r = x 2 + y 2 + z 2 , θ = arccos z x 2 + y 2 + z 2 , φ = arctan y x .
(400) Transforming the Cartesian metric to (r, θ, φ) coordinates, we find, as is well known, ds 2 = dr 2 + r 2 (dθ 2 + sin 2 θ dφ 2 ), which is normally represented by the dyad The components of the dyad with respect to the Cartesian coordinates are and the quantity ϒ is ϒ = − cot θ.

Stereographic coordinates
In stereographic coordinates, the sphere is described by means of two patches, called North and South, with local coordinates x a [N ] = (r, q, p) and x a [S] = (r, q, p) defined on each patch. Where it is necessary to distinguish between (q, p) on the North and  (x, y, z) coordinates. Construct a straight line from the north pole N of the sphere (at x = y = 0, z = 1) to a point P on the sphere, and let the line meet the plane z = 0 at X . Let the Cartesian coordinates of X be (x, y, z) = (q [S] , p [S] , 0), then the southern patch stereographic coordinates of P are (q [S] , p [S] ). The northern patch stereographic coordinates are constructed in a similar way, with the straight line in this case starting from the south pole South patches, we will use the suffix [N ] or [S] , but otherwise the suffix will be omitted. The relation to Cartesian coordinates is (Fig. 14) North: x = 2qr 1 + q 2 + p 2 , y = 2 pr 1 + q 2 + p 2 , z = r (1 − q 2 − p 2 ) 1 + q 2 + p 2 , South: x = 2qr 1 + q 2 + p 2 , y = −2 pr 1 + q 2 + p 2 , z =− r (1 − q 2 − p 2 ) 1 + q 2 + p 2 . (405) The inverse transformation is North: r = x 2 + y 2 + z 2 , q = x x 2 + y 2 + z 2 + z , p = y x 2 + y 2 + z 2 + z , South: r = x 2 + y 2 + z 2 , q = x x 2 + y 2 + z 2 − z , p = − y It is then straightforward to construct the Jacobian and to transform the Cartesian metric into the metric in (r, q, p) coordinates. We find, on both patches, ds 2 = dr 2 + 4r 2 (1 + q 2 + p 2 ) 2 dq 2 + dp 2 , which is normally represented by the dyad q A = 1 + q 2 + p 2 2 (1, i), or q A = 2 1 + q 2 + p 2 (1, i).
where in the above formula m > 0, and where The symbol means truncation to an integer; for example, 4/2 = 5/2 = 2. The (spin-weight zero) spherical harmonics in angular coordinates are then obtained by simply substituting the appropriate coordinate transformation x a [C] = x a [C] (x a [S] ) into the above formulas for Y m , and it will be found that r cancels out of the result, leaving a formula in terms only of the angular coordinates. The constant factors in the definitions are chosen so that the Y m satisfy the orthonormality condition where integration is over the unit sphere, and where dΩ = det(q AB )dφ 2 dφ 3e.g., in spherical polars dΩ = sin θ dθ dφ. The spin-weighted spherical harmonics are found by repeated application of the operators ð orð: The s Y m are defined only in the cases that |s| ≤ and |m| ≤ , and in the spinweight zero case it is usual to omit the s , i.e., 0 Y m = Y m . The s Y m satisfy the same orthonormality condition as the Y m in Eq. (418) and this is the origin of the square root factors in Eq. (419). However orthonormality does not fix the phase, and this raises issues that are pursued in Sect. 1. The s Y m constitute a large number of different cases, and are obtained from computer algebra scripts, see section "Computer algebra" in Appendix 3. The Y m satisfy the property so that ðð is the Laplacian operator on the sphere. Using Eqs. (420) and (395), it follows that ð 2ð2 Y m =ð 2 ð 2 Y m = (Λ 2 − 2Λ)Y m .
Then transforming from (θ, φ) to (x, y) coordinates, and in a neighbourhood of (x, y) = (0, 0), gives which is the expected form in the TT gauge.

Z m , the "real" Y m
Since metric quantities are real, a decomposition in terms of the Y m may introduce mode mixing between ±m modes. This can be avoided by making use of the formalism described in Zlochower et al. (2003) and Bishop (2005), and using basis functions which, in the spin-weight 0 case, are purely real; following Bishop (2005), these are denoted as s Z m .
The s Z m obey orthonormal properties similar to those of the Y m in Eq. (418), and have a relationship with ð s Z m similar to that for the s Y m in Eq. (419).

Vector and tensor spherical harmonics
Particularly within the context of gauge invariant perturbation theory, it is common practice to use vector and tensor spherical harmonics rather than spin-weighted spherical harmonics (Nagar and Rezzolla 2006). The vector spherical harmonics are defined, in the even parity case and in the odd parity case S m where C D is the Levi-Civita completely antisymmetric tensor on the 2-sphere; for example, in spherical polars θθ = φφ = 0, θφ = − φθ = sin θ . The tensor spherical harmonics are defined, in the even parity case and in the odd parity case The vector and tensor spherical harmonics are related to the ð operator and thereby to the spin-weighted spherical harmonics

Regge-Wheeler harmonics
We report below the explicit expressions of the Regge-Wheeler harmonics, (ê 1 ) i j , . . . , (f 4 ) i j , which have been introduced in Sect. 5.6.2 when discussing the numerical implementation of the Cauchy-perturbative method. In particular, the tensor spherical harmonics (ê 1 ) i j and (ê 2 ) i j in have the rather lengthy but otherwise straightforward expressions and Similarly, the tensor spherical harmonics (f 1 ) i j − (f 4 ) i j which enter in the decomposition of even-parity perturbations have the form and

Issues of convention in the definitions of spin-weighted quantities
The reader needs to be aware that different authors use different conventions in the definitions of quantities discussed in this section, and that it is therefore inadvisable to use expressions from different sources without first carefully checking the conventions used. Fortunately, the definitions of the Y m do seem to be standard throughout the mathematical-physics community. However, there are differences in the definition of the complex dyad, here denoted by q A and normalized so that q Aq A = 2; much other work uses m A with normalization m Am A = 1 so that q A = √ 2m A . Clearly the dyad definition-q A or m A -affects the definition of spin-weighted quantities. The ð operator is defined so that for a spin-weight 0 scalar V , ððV = ∇ 2 V where ∇ 2 is the Laplacian operator on the unit sphere, and the various definitions of ð all have the same magnitude. However, there can be a variation in sign. For example, the definition used in Alcubierre (2008) is −1 times that used here. This implies that the definition of s Y m is −1 times that used here for any odd s, positive or negative. However, for even s definitions of s Y m are consistent. the metric, dyad and ϒ (used in the evaluation of the operator ð) in the spherical coordinates. Each file passes these quantities to the procedure "C2P" in the file procs.map, which checks that all the relations given in Eqs. (384) and (391) are satisfied; this procedure also calculates and outputs the dyad transformed into Cartesian coordinates, i.e., Q a . The procedure "sYlm" in procs.map requires as input the various quantities defined in the coordinate-specific files, together with values for s, and m; it then calculates s Y m in the appropriate coordinate patch using the equations given in Appendix 2. The file sYlm.map reads each of the driver files in turn, and for each coordinate patch calculates all the s Y m for ≤ max , |s| ≤ min( max , s max ), with default values max = 3, s max = 2. The output is written to polars-sYlm.out, stereoNorth-sYlm.out and stereoSouth-sYlm.out.

Numerical codes
The Einstein toolkit (http://einsteintoolkit.org) contains a wide range of numerical relativity codes using the Cactus code framework (http://cactuscode.org). In particular, it contains the following thorns relevant to gravitational-wave extraction: "Extract" which implements the Cauchy-perturbative method, "WeylScal4" which implements ψ 4 extraction. The process of characteristic extraction is started in "NullSHRExtract" which constructs the worldtube boundary data, and then the various thorns listed under "PITTNullCode" are used for the characteristic evolution and the determination of gravitational-wave descriptors at J + .