{'fraction_non_alphanumeric': 0.04545637179401326, 'fraction_numerical': 0.01375477131185964, 'mean_word_length': 4.055247122545701, 'pattern_counts': {'":': 1, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Compact binary systems in scalar-tensor gravity. II. Tensor gravitational waves to second post-Newtonian order

We derive the tensor gravitational waveform generated by a binary of nonspinning compact objects (black holes or neutron stars) in a general class of scalar-tensor theories of gravity. The waveform is accurate to second post-Newtonian order beyond the leading order quadrupole approximation. We use the direct integration of the relaxed Einstein equations formalism, appropriately adapted to scalar-tensor theories, along with previous results for the equations of motion in these theories. The self-gravity of the compact objects is treated with an approach developed by Eardley. The scalar field causes deviations from the general relativistic waveform that depend only on a small number of parameters. Among the effects of the scalar field are new hereditary terms which depend on the past history of the source. One of these, a dipole-dipole coupling, produces a zero-frequency"gravitational-wave memory"equivalent to the Christodoulou memory of general relativity. In the special case of two black holes, the waveform reduces to the general relativistic waveform. For a mixed (black hole-neutron star) system, the waveform is identical to that of Einstein's theory to first post-Newtonian order, with deviations at higher order depending only on a single parameter. The behavior in these cases matches that found for the equations of motion.

There are many ways by which one may test GR with the GW signal of inspiraling compact binaries. One approach is simply to check the self-consistency of terms in the post-Newtonian sequence [17]. However, such a test can only find deviations from GR, not characterize them. Other tests may involve putting constraints on parameters in specific alternative theories of gravity [18][19][20][21][22][23][24][25][26][27][28]. A third approach parametrizes the waveform in terms of generic, theory-dependent parameters [29][30][31][32], much like the parametrized post-Newtonian formalism did for solar-system tests [16,33]. Most of these analyses have only relied on the dominant, lowest order effects in the waveform model.
A particularly important alternative to GR is the collection of scalar-tensor theories of gravity [34]. They have a long history, dating back over 50 years, and represent one of the simplest possible modifications to Einstein's theory. While solar system and binary pulsar tests put strong constraints on these theories [16], they, like all theories, have not been tested in the strong-field, dynamical regime of inspiraling compact binaries. Furthermore, they remain well motivated. For instance, many so-called f (R) theories, which modify the action of general relativity to allow arbitrary functions of the Ricci scalar, can be expressed in the form of a scalar-tensor theory [35]. These f (R) theories may explain the acceleration of the universe without resorting to dark energy. Scalar-tensor theories are also potential low-energy limits of string theory [34].
This paper is part of a series which seeks to develop the gravitational waveform for inspiraling compact binaries in scalar-tensor theories to high order in the post-Newtonian approximation. Specifically, we are interested in theories described by the action where g µν is the spacetime metric, g is its determinant, R is the Ricci scalar derived from this metric, φ is the scalar field, and ω is the scalar-tensor coupling. Note that ω = ω(φ) is not a constant; that is, we are not restricting our attention to Brans-Dicke theory. We do, however, restrict ourselves only to massless scalar fields (i.e., those without a potential). We have also written S m to represent the matter action. Note that it depends only on the matter fields m and the metric; the scalar field φ does not couple directly to the matter. This means that (1.1) is expressed in the "Jordan" frame, in which standard rods and clocks measure distances and times. All of our work will be done in this frame. An alternative representation is the "Einstein" frame, related to the Jordan frame by a conformal transformation [36]. The first step in the construction of gravitational waveforms is to calculate the equations of motion for the compact objects. This was the subject of the first paper in this series [37]. Mirshekari and Will (hereafter MW) computed the equations of motion to order (v/c) 5 (2.5PN) beyond the leading term. They made use of a method known as direct integration of the relaxed Einstein equations (DIRE), based on the original framework of Epstein and Wagoner [38] and then extended by Will, Wiseman, and Pati [39][40][41][42]. This approach has been shown to give identical results to other methods, including the "post-Minkowskian" method [8], the Hamiltonian approach [43], the "strong-field pointparticle limit" strategy [44], and the "effective field theory" method [45]. It is also easily adapted to scalar-tensor theories.
In the adapted DIRE method, the scalar-tensor field equations are first rewritten in a "relaxed" form: flat-spacetime wave equations for a "gravitational field"h µν and a modified scalar field ϕ. The wave equations are simplified by the choice of a particular coordinate system, represented by a gauge condition onh µν . Together, the wave equations and gauge condition contain all the content of the full field equations.
The wave equations are then solved formally using a retarded Green's function, valid everywhere in spacetime. To convert them to a more useful form, however, these formal solutions are evaluated differently in different regions of spacetime: In the "near zone" close to the source (defined in Sec. II B below), the integrals are expanded using a slow-motion approximation. Far away from the source, in the "radiation zone," a special coordinate transformation is used to evaluate the solutions. The total solution for each field (h µν and ϕ) is then the sum of the two separate solutions; any terms dependent on the arbitrary boundary between zones must cancel out.
In this paper, we calculate the tensor gravitational waveform to order (v/c) 4 (2PN) beyond the leading-order "Newtonian" quadrupole. In doing so, we also make use of the adapted DIRE approach. The difference with MW is fundamentally simple: In that paper, the gravitational fieldh µν was evaluated in the near zone, in order to calculate how it affects the motion of the compact bodies. We, by contrast, evaluateh µν (or more specifically,h ij ) in the radiation zone, where it will be measured by a GW detector. Our procedure very closely follows that of Wiseman and Will [40] (hereafter WW), although the notation is updated to match that of MW. Along the way, we must make use of MW's results for the equations of motion.
As discussed in MW, one key concern in this study is the treatment of the compact bodies' internal gravity. Since scalar-tensor theories do not obey the strong equivalence principle, the motion and gravitational-wave emission of a binary depend on the internal composition of its constituent bodies. To handle this effect, we have adopted the approach of Eardley [46]. We treat the matter stress-energy tensor as a sum of delta functions located at the position of each compact object. However, instead of assigning each body a constant mass, we let the mass be a function of the scalar field, M A = M A (φ). This gives the matter action an indirect dependence on φ, even though we still work in the Jordan frame. In the final waveform, this dependence will appear as the "sensitivity" of the mass to variations in the scalar field, as well as derivatives of this quantity. (The subscript 0 means that the derivative should be evaluated using the asymptotic value of the scalar field, φ 0 .) In the weak-field limit, the sensitivity is proportional to the Newtonian selfgravitational energy per unit mass of the body. For neutron stars, the sensitivity depends on the mass and equation of state of the star, with typical values 0.1-0.3 [47,48]. For black holes, s = 0.5, and all derivatives vanish.
Recently, numerical simulations of compact binaries in scalar-tensor gravity have shown that the sensitivities of neutron stars can change dramatically during the late inspiral [49][50][51]. This "dynamical scalarization" effect is a generalization of the "spontaneous scalarization" of individual stars discovered by Damour and Esposito-Farèse [52,53]. Since we assume that the sensitivities are constant in time, our work does not capture this effect. In any case, dynamical scalarization only becomes relevant during the late portion of the inspiral, when our use of the post-Newtonian approximation also becomes invalid.
We find that, at 0PN and 0.5PN order, the gravitational waves are identical to those in general relativity, except for two changes. First, each explicit factor of total mass m = m 1 + m 2 contained in the expression is modified to αm, where Here ω 0 ≡ ω(φ 0 ). Second, the overall waveform is scaled by (3 + 2ω 0 )/(4 + 2ω 0 ). At higher PN order, the deviations become more complicated. However, they are still described by a relatively small number of parameters, the same as those used in MW to characterize differences between the GR and scalar-tensor equations of motion. Almost every term in the waveform has a counterpart in the GR waveform with the same basic structure. Terms which are entirely new result from the existence of a scalar dipole moment I i s (defined in Sec. IV A 5 below). Unlike the regular dipole moment I i , the scalar dipole moment cannot be made to vanish by choosing center-of-mass coordinates. This is a direct consequence of scalar-tensor theories violating the strong equivalence principle. While the main impact of I i s is in the generation of scalar dipole radiation, its presence also has a strong effect on the tensor waveform.
In general relativity, integration over the radiation zone produces terms which enter the final waveform beginning at 1.5PN order. Some of these terms are "instantaneous"; that is, they depend only on the binary's state at a particular time. (Because the waves travel at the speed of light c, they depend on the binary's state not at the current time t, but at a retarded one τ = t− R/c, where R is the distance to the system.) Other terms are "hereditary": They require an integration of the binary's behavior from an infinite time in the past until the moment the waves are emitted. At 1.5PN order, the hereditary contribution involves a logarithmic factor and is known as the gravitational-wave "tail." It arises from a coupling between monopole and mass quadrupole moments of the source and can be physically described as backscattering of the outgoing radiation on the curved spacetime of the binary. Monopole-mass octupole and monopole-current quadrupole couplings create further tail terms at 2PN order. Higher PN orders include more tails, as well as "tails of tails" arising from three-moment couplings.
The addition of the scalar field to the radiation-zone integrals produces no additional tail terms. However, it does produce another kind of hereditary term, which involves an integral of moments of the system over its entire past history without a logarithmic multiplier. At 1.5PN order, the hereditary term involves a mass dipole-mass dipole coupling, while at 2PN order, the coupling is between the mass dipole and the mass quadrupole. Both of these produce oscillatory terms, like all the other pieces of the waveform. However, the dipole-dipole integral also produces a zero-frequency (DC) term, which grows secularly throughout the inspiral. Such terms are often referred to as "nonlinear gravitational-wave memory," or "Christodoulou memory" [54][55][56][57]. Unlike oscillatory terms, they cause a permanent change in the state of a GW detector. In general relativity, the nonlinear memory does not appear until 2.5PN order, where it is driven by a mass quadrupole-mass quadrupole coupling. Even with the new dipole-dipole effect, the nonlinear memory formally enters the waveform at relatively high post-Newtonian order. However, because it is integrated over the past history of the binary, the memory's actual effect is comparable in magnitude to lower order post-Newtonian terms. It is, therefore, potentially detectable by gravitational-wave detectors [57,58].
One interesting limit of our results is the case where both compact objects are black holes. In scalar-tensor theory, isolated black holes behave identically to those in general relativity. MW verified that a similar statement is true about binary black holes. Specifically, they showed that the 2.5PN equations of motion in scalar-tensor theory are identical to those in general relativity, except for an unmeasurable rescaling of masses. They then conjectured that the same should hold true for the gravitational waves emitted by the binary. Our work shows that this conjecture is correct for tensor gravitational waves, at least to 2PN order. MW also found special results for mixed systems, those containing one black hole and one neutron star. In that case, the equations of motion are identical to GR, with mass rescaling, up to 1PN order. Beyond that order, they deviate from Einstein's theory, but the deviation depends only on a single parameter. Unfortunately, this parameter is the same in Brans-Dicke and more generalized scalar-tensor theories. We have found that all of these properties carry over to the tensor gravitational waves emitted by a mixed system.
To find the gravitational-wave phasing, we require an expression for the energy loss. That, in turn, requires the scalar gravitational waveform. It will be considered in a follow-up paper. While the procedure is the same as for the tensor waveform, the actual calculation is somewhat more complex and lengthy. With both tensor and scalar pieces in hand, we can write down the full response of a detector to the inspiral of a circular binary. Finally, we will be able to use the final waveform in a parameter estimation study. We wish to investigate how well detectors like Advanced LIGO can measure differences between GR and scalar-tensor theories, as well as how the additional terms in the waveform influence the measurement of astrophysical parameters, like those described above.
The outline of the paper is as follows: Section II presents the field equations for the tensor and scalar fields, as derived from the action (1.1). We cast them into a relaxed form and show how they can be solved for field and source points in different regions of spacetime. We also discuss how the tensor and scalar fields affect a gravitationalwave detector. Section III reviews the results from MW in the near zone surrounding the source. Post-Newtonian expansions of the near-zone source are needed to calculate the gravitational waveform.
Section IV describes the calculation of the Epstein-Wagoner moments, the fundamental pieces needed to derive the near-zone contribution to the gravitational waveform. The calculation of the two-index moment I ij EW is described in some detail in order to clarify the major issues involved. Section V shows how the moments can be converted from a generic N -body scenario to the specific two-body case we consider. We also discuss how the equations of motion, taken from MW, are used to expand time derivatives of the moments.
Section VI leaves these pieces aside and presents the other half of the puzzle: the radiation-zone contribution to the gravitational waves. First, we must calculate the tensor and scalar fields in the radiation zone, far from the compact objects. They have both near-zone and radiation-zone sources. Then these fields are used to calculate the GWs. Terms produced here enter the final waveform at 1.5PN and 2PN orders and include the hereditary effects described earlier (tails and memory). Finally, in Sec. VII, we present the full 2PN tensor gravitational waveform for a nonspinning compact binary in massless scalar-tensor theory. We also discuss some features of the results in more detail.
In this paper, we use units in which c = 1. We do not set G = 1; as we shall see, the effective Newtonian gravitational constant depends on the asymptotic value of the scalar field. Greek indices run over four spacetime values (0, 1, 2, 3), while Latin indices run over three spatial values (1,2,3). We use the Einstein summation convention, in which repeated indices are summed over. We use a multi-index notation for products of vector components: x ijk ≡ x i x j x k . A capital letter superscript denotes a product of that dimensionality: x L ≡ x k1 x k2 · · · x k l . Angular brackets around indices denote symmetric, trace-free (STF) products (see Appendix B for details). Finally, we use standard notation for symmetrized and antisymmetrized indices, e.g. x (i y j) ≡ (x i y j +x j y i )/2 and x [i y j] ≡ (x i y j −x j y i )/2. At times we use a bar to separate indices which should be symmetrized from those which should not, e.g. x (i| y j z |k) ≡ (x i y j z k +x k y j z i )/2.

A. Field equations and relaxed form
The field equations for theories described by (1.1) are given by As stated above, we work in the Jordan representation of the theory, so that g µν is the physical metric. The quantity G µν is the Einstein tensor constructed from this metric, while φ is the scalar field, ω(φ) is the coupling function, T µν is the stress energy of matter and nongravitational fields, and T ≡ g αβ T αβ is its trace. We use commas to denote ordinary derivatives. Semicolons denote covariant derivatives (taken using g µν in the usual way), and ✷ g ≡ g αβ ∂ α ∂ β is the d'Alembertian with indices raised by the metric. In the Jordan representation, the derivative ∂T /∂φ is not present for normal matter, but will be present for gravitationally bound bodies, as we will discuss in Sec. III B. We assume that far away from the sources, the metric reduces to the Minkowski metric, η µν , and the scalar field tends to a constant φ 0 . We introduce a rescaled scalar field, and a "gothic" version of that metric with its indices raised, Hereg is the determinant ofg µν . We can then define the gravitational field as We use a tilde to differentiate this field from the gravitational field defined in general relativity, h µν , which has the same definition but with g µν ≡ √ −gg µν in place ofg µν . [See, for instance, WW (2.2).] We impose the Lorenz gauge conditionh µν ,ν = 0 . (2.6) Then the field equation (2.1a) reduces to where ✷ η ≡ η αβ ∂ α ∂ β is the flat-spacetime wave operator and the source is Here T µν is the stress energy of matter and nongravitational fields. In our case, we have no other fields, so this represents the compact piece of the source. The quantity Λ µν represents the gravitational-field contribution to the stress energy: The scalar field equation (2.1b) can also be written as a flat-spacetime wave equation,

B. Solution of the wave equations
The wave equations (2.7) and (2.12) can be solved formally in all spacetime by using a retarded Green's function, The delta function in both these equations restricts the integration to being over the past flat-spacetime null cone C emanating from the field point (t, x). To obtain explicit solutions, we divide the spacetime into two regions. Define the characteristic size of the source as S. We assume the bodies move at velocities v ≪ 1. Then the near zone is defined as the area with |x| = R < R, where R ∼ S/v is the characteristic wavelength of gravitational radiation from the system. (We use capital R to denote the distance from the binary's center of mass to a field point in order to avoid confusion later with r, the orbital radius of the binary.) Everything outside the near zone (R > R) is the radiation zone.
We carry out the integrals (2.14a) and (2.14b) in two pieces: one integral over the near zone and one over the radiation zone. Each is done by a fundamentally different method. In the end, the final solution is the sum of the two pieces. WW and [41] showed explicitly that any terms dependent on the boundary radius R in one-half of the integral would be exactly cancelled by pieces in the other half of the integral, leaving the final answer, as expected, independent of this arbitrary parameter. In our work, we simply assume this property and ignore any terms which depend on R.
The integrals are also evaluated differently depending on what field point x we are interested in. For instance, MW calculatedh µν and ϕ at field points x in the near zone, where the bodies are located. The near-zone fields were then used to calculate the equations of motion for the bodies. In this paper, we are interested in the gravitational waves, so we will ultimately want the fields evaluated in the radiation zone. However, we will still need the near-zone fields as source terms (see Λ ij and Λ ij s above). With two integration regions and two possibilities for field points, there are a total of four distinct ways to evaluate (2.14a) and (2.14b). For a complete description of these methods, refer to WW and [41]. Here we give only a brief summary. For field points in the near zone and integration over the near zone, we can treat the retardation as a small perturbation. Expanding Here N is the three-dimensional hypersurface representing the intersection of the past null cone C and the near-zone world tube. After the expansion, the actual integration takes place over M, the intersection of the hypersurface t = const with the near-zone world tube. Note that in the near zone, the slow motion approximation v ≪ 1 means that each time derivative corresponds to an increase of one-half post-Newtonian order. For field points in the near zone and integration over the radiation zone, we recognize that the source contains only field terms which are themselves retarded. Therefore, it is prudent to change variables and integrate over τ Here, the notation C − N denotes that this is the integration over all pieces of the past null cone which do not intersect the near zone; that is, they are in the radiation zone. For each τ ′ , the inner pieces integrate over the intersection of C with the future null cone emanating from the center of mass of the system at time τ ′ . The τ ′ integration is then a summation over all such future-directed null cones, starting from the infinite past and continuing until the cones no longer overlap. It turns out that the contributions from these integrals only come into play at higher post-Newtonian order than was considered in MW. For field points in the radiation zone and integration over the near zone, we expand the entire integrand in powers of |x ′ |/R and findh µν Here M is again the intersection of the near-zone world tube with a constant-time hypersurface; however, in this case, that time is the retarded time τ = t − R.
For gravitational waves, we can simplify (2.17a) and (2.17b) in two ways: (1) We are only interested in the spatial piece of the tensor,h ij . (2) Gravitational-wave detectors operate in the regime R ≫ R, which we call the far-away zone. (It is, of course, a subset of the radiation zone.) Therefore, we can expand (2.17a) and (2.17b) in powers of 1/R, keeping only the lowest order term. The results arẽ whereN ≡ x/R is the direction from the source to the detector. The Lorenz gauge condition implies a conservation law for the source, Using this, we can rewrite the tensor equation as where the I M+2 EW are known as "Epstein-Wagoner" (EW) moments [38]. They are given by The largest piece of the work in this paper is the calculation of these EW moments to the necessary post-Newtonian order. The use of the conservation law to rearrange the two-and three-index EW moments results in the "surface" moments I ij EW (surf) and I ijk EW (surf) , They are evaluated on ∂M, a sphere of radius R bounding the hypersurface M;n l is a radial unit vector pointing outward from this boundary. Finally, if the field point is in the radiation zone and the integration is over the radiation zone, we can again use a change of variable. The results are identical to (2.16a) and (2.16b), except with a different limit of integration, (2.24b) To find the gravitational waves, we merely ignore all nonspatial components (for the tensor waves) and consider only terms with 1/R dependence. These radiation-zone integrals will produce hereditary terms in the final GW signal, including tail and memory effects.

C. Effect on GW detectors
To find the effect of the (tensor) gravitational field and the scalar field on a GW detector, we first need to convert back to the physical metric. The inverse metric is The determinantg is given byg whereh ≡ η µνh µν is the trace ofh µν . We are not concerned with terms of orderh 2 and higher because we are working in the far-away zone, whereh ∼ 1/R. Putting everything together, we get Here Ψ ≡ ϕ − 1, and we have also made use of the fact that Ψ ∼ 1/R in the far-away zone. The physical metric is therefore whereh µν = η µα η νβh αβ is lowered using the Minkowski metric. A gravitational-wave detector works by measuring the separation ξ between test masses. If the distance between the test masses is small compared to the wavelength of the GWs, and the masses move slowly, the separation obeys the equationξ where dots denote time derivatives and R 0i0j are components of the Riemann curvature tensor. We can use the metric to calculate them, Here "TT" designates the transverse-traceless projection of the gravitational field tensor, which satisfies the conditionŝ We can find the TT part of a tensor by using the projection operator, where P pq = δ pq −N pN q is the transverse projection operator. The fact that only the TT piece ofh ij contributes to the measured GWs will allow us to simplify our calculation by prematurely dropping terms which cannot possibly produce a TT contribution at the end of the day. Note that the scalar field Ψ will contribute a transverse "breathing" mode to the GW signal; it will be treated in a subsequent paper in the series.

III. CALCULATION OF THE NEAR-ZONE SOURCE
In this section, we review results from MW for the source τ µν in the near zone. These results are needed to calculate the Epstein-Wagoner moments and thus the near-zone contribution to the gravitational waveform.
A. Potentials, fields, and field source Following MW, the compact source can be described in terms of densities [59], We can then define a number of Poisson-like potentials. For example, given a generic Poisson integral for a function f (t, x), then the basic "Newtonian" potential is We also have a scalar equivalent, The σ subscript clarifies that these potentials use the densities defined in this section. Later, we will convert to a new density which reflects the specific source (compact binaries) we are studying. We delay defining the rest of the potentials until then. Expressions for all other σ-density potentials can be found in MW (3.12) and (3.13). Note that the generic Poisson integral has the property This will be very useful throughout the calculation. For convenience, we can rewrite the fields ash (3.6e) We use this notation in all spacetime. In the near zone, where the post-Newtonian counting parameter ǫ ∼ v 2 ∼ m/r. Here m is the mass of the system, r is a typical distance scale, and v is its characteristic speed. (Later, r will be the orbital radius of a binary, and v will be the magnitude of its relative velocity.) To obtain expressions for the near-zone fields, we solve (2.15a) and (2.15b) iteratively. At lowest order, we only solve for N and Ψ, the other quantities being intrinsically higher order. The sources only include compact terms: . Using the definitions of the potentials, we find is the effective gravitational constant. The definition is chosen so that, for a perfect fluid with no internal gravitational binding energy, the metric component g 00 matches the result from general relativity, g 00 = −1 + 2GU σ . We do not set G equal to 1, since it depends on the asymptotic value of the scalar field φ, which could potentially vary in time.
The other parameter is For the next order, we need to begin evaluating the field terms of the source, Λ µν and Λ µν s . These are given explicitly in terms of N , K i , B ij , B, and Ψ in MW (3.4) and (3.5). We can plug in (3.7a) for N and (3.7b) for Ψ. This will be enough to get the fields at next order, as shown in MW (4.10). This procedure is iterated until the fields and sources are obtained to the necessary order. In the next section, we will need the expressions for τ µν in order to calculate the Epstein-Wagoner moments.

B. Matter source
The previous section describes a generic matter source characterized by densities σ, σ i , σ ij , and σ s (the "σ densities"). All of the expressions for near-zone fields and sources in MW are written in terms of these generic densities. We now convert to a more realistic source for the situation we are considering: one made up of an arbitrary number of compact objects. (We later specialize to the two-body case, but many steps of the calculation are valid for a general system.) Since a compact object is gravitationally bound, its total mass depends on its internal gravitational energy. This, in turn, depends on the effective local value of the gravitational coupling. In scalar-tensor theory, the coupling is controlled by the value of the scalar field φ in the vicinity of the body.
To deal with this complication, we use the approach of Eardley [46]. In his method, we consider the compact objects to be point masses, with a mass M (φ) that is a function of the scalar field. The stress-energy tensor is then given by (3.10) Here u µ A is the four-velocity of body A, and τ is the proper time measured along its world line. (This is the only instance in which we use the symbol τ for this purpose.) The dependence of mass on φ is what leads to the ∂T /∂φ term in (2.1b), which would not normally be present in the Jordan frame. (Remember that in the Jordan representation, the scalar field does not directly couple to the matter. This indirect coupling is merely a way of treating the complexity of the source in scalar-tensor theory.) We expand M A (φ) about the asymptotic value φ 0 , where m A ≡ M A0 . We define the sensitivity and its derivatives as and so on. Note that s ′ has the opposite sign of the equivalent quantity in [33] and [60]. If we define a new density the stress energy becomes where v µ = (1, v) is the ordinary velocity. The various velocities and the sensitivity s technically should have body labels, but they will each pick one up when multiplied by the delta function in ρ * . We use this convention frequently in the steps to follow. Returning to the definitions of the σ densities, we find and With these quantities, we can also rewrite By substituting the metric into (3.15a)-(3.15d), the σ densities can be written in terms of the ρ * density as a post-Newtonian expansion. These expressions are given in MW (5.12). We will need them to translate τ µν from the expressions in MW to the versions we need. We also define new potentials based on the ρ * density. For instance, More generally, plus natural generalizations of X and Y like X i , X s , etc. All the potentials listed in this paper can be expressed in terms of P (f ) (3.2) or these basic forms. They are listed in Appendix A. MW (5.13)- (5.22) show how to convert between many σ-density and ρ-density potentials (e.g., U σ and U ). These expressions will also be useful in converting τ µν .

IV. EPSTEIN-WAGONER MOMENTS
In this section, we calculate the Epstein-Wagoner moments which are needed to generate the near-zone contribution to the gravitational waveform.
A. Two-index moment I ij

EW
We begin with the two-index moment. Its lowest order piece will produce the lowest order gravitational waves, which we define as 0PN. It is given by For clarity, the moment has been split into three pieces: (1) I ij C , in which the integral is taken over the compact part of τ 00 ; (2) I ij F , in which the integral is taken over the field part of τ 00 ; and (3) I ij S ≡ I ij EW (surf) , the surface moment. To calculate the first two pieces, we need the source τ 00 evaluated to 2PN order, or O(ρǫ 2 ). This is found in MW (4.14a), but it is expressed using the σ density. Converting to the ρ * density, we find Recall that all the potentials are defined in Appendix A. We have introduced the quantity Later, we will also need The compact moment can be written down by inspection, Here body A (for example) has position x A , velocity v A , and acceleration a A . The distance between bodies A and B is r AB , and the unit vectorn AB = x AB /r AB points from body B to body A.
To calculate the field moment I ij F , we must evaluate a series of integrals involving the 1PN and 2PN potentials. Examining (4.2), we see that there are 24 different integrals. However, integrals such as M (∇U ) 2 x ij d 3 x and M (∇U s ) 2 x ij d 3 x are essentially identical: The functional form of the integrand is the same, with the only difference being the addition of sensitivity factors 1 − 2s. If we count such integrals together, there turns out to be 15 unique pieces to compute. (Many of these also share similar fundamental components, but we count them separately due to different factors of velocity or acceleration.)

1PN field integral
The simplest integral is the 1PN term M (∇U ) 2 x ij d 3 x. (When discussing these integrals, we will ignore the constants in front.) We will calculate it explicitly in order to point out some techniques which will be used throughout the calculation. First, we integrate by parts, The surface integral is evaluated on the boundary of the near zone at a constant (retarded) time; therefore, ∂M is a sphere with radius R. We can write x i = Rn i and d 2 S k = R 2nk , wheren i is a unit vector normal to the surface. We then expand U and U ,k in inverse powers of R and look for any R-independent terms in the surface integral. We find that The center of mass (CM) of the system is defined as If we choose to work in coordinates where x CM = 0, then the second term in the above expansion is one post-Newtonian order higher than the others. We use the order parameter ǫ to mark this. With that exception, we only care about the R scaling and the number of unit vectorsn i . The derivative of the potential scales like So the only R-independent terms in the surface integrand scale like ǫn aij . When integrated over the surface, an odd number ofn i vanishes. (See Appendix B.) Therefore, the surface integral makes no contribution. The first volume integral can be simplified using , making it essentially a compact integral. The solution is easily seen to be The second volume integral can be integrated by parts again, The surface integral again vanishes by virtue of integrating an odd number of unit vectors. Meanwhile, the volume integral can be ignored because it will not survive the transverse-traceless projection. So (4.9) alone is the value of the 1PN integral. Throughout this work, we routinely drop terms which will not, in the end, survive the transverse-traceless projection. For convenience, we refer to these terms as "non-TT." Non-TT terms are easy to identify by sight: Any term containing δ ij ,N i , orN j is non-TT. [This can be checked by applying the TT projection (2.32) to such terms.] For the two-index moment, δ ij pieces are the only ones we drop. For higher order moments, terms like δ ik , δ jk , δ il , δ jl , ... are also dropped. This is because the final expression for the waveform (2.21) contracts the EW moments with direction vectorsN k ,N l , ... for all dummy indices (i.e., all those besides i and j). Therefore, a term like δ ik will introduceN i in the final waveform, and that is non-TT. Note that terms like δ kl (or any other involving two dummy indices) must be kept.
As discussed earlier, the integral M (∇U s ) 2 x ij d 3 x is essentially the same as the one we have just calculated. One slight difference does occur in the surface integrals: A does not vanish, even to lowest post-Newtonian order, the 1/R 2 (1/R 3 ) term in the expansion of U s (U ,k s ) will be at the same order as all the other terms. However, this makes no difference: The surface integrals still vanish. The final answer is exactly the same as (4.9), only with sensitivity factors added, (4.11)

2PN two-potential field integrals
The other integrals enter at 2PN order. Eleven of them involve only two potentials; of these, all but the one involving P kl 2 (or P kl 2s ) can be solved using a straightforward procedure. Consider the following example: Here the surface integral vanishes for the same reason as above (i.e., an odd number ofn i ). The first volume integral is evaluated easily using The tricky part is the second volume integral. We can write it as It can be evaluated using techniques developed in WW. We first change integration variables from x to y = x − x A . Looking just at the main piece of the integral, this gives, in our particular example, (4.14) There are two cases to consider: A = B and A = B. In both cases, the infinite series of surface integrals vanishes: The terms either depend on R or average to zero because of an odd number of unit vectors. When A = B, x AB = 0, making the volume integral easy to evaluate. It also vanishes, for the same reasons as the surface integrals. When A = B, we make use of the following expansion: where Y lm are the spherical harmonics, and r <(>) denotes the lesser (greater) of r AB and y. We substitute this expansion into the volume integral and then express all products ofŷ i in terms of symmetric, trace-free (STF) productsŷ L ′ . (See Appendix B.) Here L ′ denotes an l ′ -dimensional STF combination. We can then perform the angular integration using In our example, we find Finally, the radial integral is evaluated using where we have dropped terms dependent on R. For our example, we find where in the last step we have converted back to non-STF notation and discarded the non-TT δ ij terms. This expression can then be used to evaluate the original integral, with the mass and velocity factors added back in. As stated above, ten of the two-potential integrals can be evaluated with this step-by-step procedure: First, integrate by parts so that one potential is undifferentiated. (In some cases, this step is unnecessary because one potential is already undifferentiated.) Next, if any piece can be converted to a compact integral through relations like ∇ 2 Φ 1 = −4πρ * v 2 , evaluate that piece. Then write out the potentials and their derivatives explicitly in terms of masses, positions, velocities, and accelerations. When two derivatives are taken, care must be taken to add appropriate delta functions, where "norm" denotes the derivative computed from the definitions of the potentials. The extra terms are needed to ensure the right answer when the doubly differentiated potentials are integrated in a sphere around the point mass position Next, carefully examine any surface integrals. It turns out that no surface integral contributes to the final answer, for one of three reasons: (1) it has no R-independent terms, (2) the R-independent terms vanish upon integrating over the surface, or (3) the R-independent terms average to something proportional to δ ij . The last type of term does not vanish, but it is non-TT, and we can ignore it.
For the remaining volume integrals, change variables from x to y = x−x A , where A is the label on the differentiated potential. Check that the surface integrals so generated and A = B volume integrals contribute nothing, for one of the three reasons above. Then integrate the A = B volume integral over angle and radius using (4.15), (4.16), and (4.18), keeping only TT terms at the end of the calculation. Many of the A = B integrals which arise in this process appear in several of the ten "main" integrals, and the results can be reused with the appropriate coefficients multiplied in.

P kl 2 field integral
The final two-potential integral is M U ,kl P kl 2 x ij d 3 x (and its counterpart with P kl 2s ), where U ′ is the usual potential written as a function of x ′ . We can integrate over x (unprimed) first, using the techniques of the previous section, except with x B → x ′ . Dropping the primes on the remaining integration variable, the result is Following WW (D5), we can rewrite this as These six terms can be evaluated individually. For the first four, we first integrate by parts and find a vanishing surface integral, a simple-to-compute volume integral (i.e., one involving a Laplacian), and a more difficult volume integral. This last piece is integrated by parts again, leading to another vanishing surface integral and a final volume integral. The final volume integrals in each case seem difficult to evaluate, but fortunately they cancel in pairs when the four terms are combined. The fifth term behaves similarly, except it is possible to evaluate all the integrals eventually. It involves four integrations by parts and three volume integrals which convert to compact integrals by means of a Laplacian. The sixth piece is the most difficult, where For B = C = A, we make a slightly different substitution (since B, not A, labels the differentiated potentials): plus vanishing surface terms. To evaluate this, we use a new expansion in spherical harmonics, The final case, A = B = C, is very complicated. It is worked out in Appendix D of WW. Switching body labels to match their notation, we get where This term is irrelevant for compact binaries, but we keep it for completeness. The final answer for this field integral is There are three field integrals involving the product of three potentials. The simplest is The surface integral is non-TT, and the first volume integral is trivial. Integrating the second volume integral by parts gives another surface integral and another volume integral, both of which are non-TT. The other two field integrals are similar but involve combinations of U and U s . To solve each individually, we would have to compute a three-potential subintegral in the manner of (4.24) above. However, when the two field integrals are appropriately combined in the EW moment, this subintegral cancels and can therefore be ignored. Ignoring non-TT contributions, we find Each of the subintegrals is trivial to evaluate. With the evaluation of the three-body integrals, we have now completed all 15 field integrals necessary to compute the two-index EW moment.

Surface moment
We can rewrite the surface moment (2.23a) as where Λ ij T = Λ ij + Λ ij s . (Recall that the boundary of the near zone is far beyond the compact source: R ≫ S.) We are again only interested in terms which do not depend on R. Since the first term of the integrand multiplies Λ ij T by R 3 , we look for pieces of Λ ij T with R −3 dependence. The second term of the integrand adds an extra factor of R, but the derivative on Λ ij T reduces the overall scaling by R to compensate. To survive the angular integration, a piece of Λ ij T must contain an even number of unit vectors. The two unit vectors in the first term retain then parity of Λ ij T ; the three unit vectors in the second term flip it, but the derivative flips it back.
We begin with the first piece of the source, Λ ij . Because of the two time derivatives, we need to know it to 3PN order [i.e., O(ρǫ 3 )] in order to find the final moment to 2PN order. Unfortunately, MW does not contain an expression for Λ ij to this order. However, because of the way we have defined our quantities, Λ ij will have exactly the same form as the general relativistic version found in Eq. (4.4c) of [41]. The only difference is that the fields N , K i , B ij , and B we plug into (4.4c) are pieces of our new gravitational fieldh µν instead of the GR field h µν .
We must expand Λ ij in the vicinity of r = R and look for terms with R −3 dependence. The first step is to expand the individual fieldsh µν . Forh 00 = N , the lowest order piece is given above by (3.7a), and its expansion is just the expansion of U given in (4.7). Higher order pieces of N are given in MW (4.10a), (4.10e), and (4.15a). The various potentials in these expressions, when converted to ρ * density, can be expanded similarly to U . As with earlier surface integrals, we only care about a term's post-Newtonian order, its dependence on R, and its number of unit vectorŝ n i . With this restriction, we see that all potentials of the same family [i.e., P (or Σ), X, and Y ] share the same basic expansion. The only relevant differences occur for terms containing A m A x i A or its derivatives. As seen in the expansion of U itself, these terms are one post-Newtonian order higher than the others in the same expansion. In the end, we find R 3 + · · · + O(ǫ 7/2 ) .

(4.33)
For simplicity, we have omitted the exact factors ofn. Instead, we just use the superscript (e) to represent a sum of one or more terms with even parityn (e.g., 1 +n ab ), (o) for one or more terms with odd parityn, and (e + o) when terms are present of both parities. Note that the O(ǫ 5/2 ) term is different from the others; as seen in MW (4.10e), it contains no potentials and is, in fact, independent of position x. The other fields can be expanded as The expansion for B ij , which we need only to O(ǫ 2 ), is the same as for B. The O(ǫ 5/2 ) piece of B, like that of N , is independent of position. [See MW (4.10f).] To calculate Λ ij , we need to take spatial and time derivatives of these fields. For spatial derivatives, we merely divide each term (excepting those independent of position) by R and change the parity ofn. Time derivatives add an extra factor of ǫ 1/2 each, while also affecting the coefficients we have chosen to ignore. Only some of the pieces of Λ ij produce terms which scale like R −3 . Ignoring coefficients, these are N ,i N ,j , δ ij (∇N ) 2 , N ,(iK j) , δ ij N ,kK k , N ,(i B ,j) , δ ij ∇N ·∇B, δ ijṄ 2 , NB ij , δ ijṄḂ , and δ ij NṄ 2 . The first six of these scale liken (o) /R 3 , the next three like (δ ij +n (o) )/R 3 , and the last like δ ij /R 3 . (For NB ij , the δ ij factor is hidden inside B ij .) Terms with an odd number of unit vectors will vanish trivially during the angular integration of (4.32). The even parity terms will not vanish. However, plugging into (4.32), we see that the angular integral reduces to the average ofn ij . This produces δ ij , and so these terms are non-TT. We conclude that Λ ij makes no relevant contributions to the surface moment.
We turn next to potential contributions from Λ ij s . We will need to expand it to O(ρǫ 3 ) first. The result is (4.36) The expansion of Ψ is As with N and B, the O(ǫ 5/2 ) term does not contain a potential; in this case, it has a term linearly dependent on the position (so thatn (o) =n a ) in addition to one independent of the position. We again find a number of terms which scale like R −3 . They derive from the Ψ ,i Ψ ,j , δ ij (∇Ψ) 2 , δ ijΨ2 , δ ijΨ2 N , and δ ij ΨΨ 2 pieces of Λ ij s . When inserted into (4.32), all of these but one either vanish because of oddn parity or are ignored because they produce something which is non-TT. The surviving term is (  Plugging into (4.32), we find that the only contribution to the surface moment is (4.39)

Final two-index EW moment
We now add up the results for the 15 field integrals (and their variations), multiplied by the appropriate coefficients, to find the total field moment I ij F . We then add this to the compact moment I ij C and the surface moment I ij S to form the total two-index EW moment. It contains pieces at 0PN, 1PN, 1.5PN, and 2PN order, (4.40) The notation (d/dt) −2 means that two antiderivatives must be taken of the following expression. In practice, we do not need to worry about doing this, since we will eventually take two time derivatives of I ij EW when evaluating (2.21). In the GR limit ζ = 0 (and for completeness, G = 1), this expression reduces to WW (4.17). [Note that (4.17) of WW has a sign error on the term containing ·n AB ]}/2.] There are two important caveats about this expression. First, the surface moment has been written using the lowest order piece of I i s . Note also that the terms − C =B m Cn k BC /r 2 BC from (4.29) do not appear in that form. Similar terms occur in the related integral involving P kl 2s . The two sets of terms can be combined and then simplified using This is the Newtonian equation of motion as defined in MW (6.1). With this substitution, I ij EW is simpler and easier to compare to the WW results (which use the same trick). Because of these two choices, (4.40) is good only for calculations at 2PN order. When going to higher order, the expression should be reverted to its more generic form (not shown) before adding the explicit higher order contributions; otherwise, the final answer will not be accurate.

EW
The rest of the moments are calculated in much the same way as I ij EW . For the three-index moment, it is useful to rewrite it as Here the volume integral can again be split into compact and field components,Ĩ ijk C andĨ ijk F . The surface moment is given by (4.44) To get I ijk EW to 2PN order, we technically need to know the source τ 0i to O(ρǫ 2 ). However, there are no terms at that order, so really we need it only to O(ρǫ 3/2 ). Converting MW (4.12b) to the ρ * density, we find The compact moment is given bỹ There are only two types of field integral, and both can be calculated by the method of Sec. IV A 2. Finally, we must investigate the surface moment. The procedure is similar to that for the two-index surface moment, with three essential differences. First, there is only one time derivative, compared to two in the two-index case. This means that we only need to consider Λ ij T to O(ρǫ 5/2 ). Second, we are interested in terms in Λ ij T which have R −4 dependence. Finally, if a piece of Λ ij T is to survive the angular integration, it must have an odd number of unit vectors. Again, only one term contributes to the surface moment: the one involving the scalar dipole (4.38). The final expression is given by Like the two-index surface moment, this expression uses only the lowest order form of I i s . This is sufficient for our purposes, but care must be taken in any future work to higher post-Newtonian order.
Adding everything up, we get a final expression forĨ ijk EW . It contains pieces at 0.5PN, 1.5PN, and 2PN order, This reduces to WW (4.22) in the GR limit.

EW
To evaluate the four-index moment, we need τ ij to O(ρǫ 2 ), which is not given in MW. We use Eq. (4.4c) of [41] and (4.36) to find We can ignore all terms with δ ij , because they will be non-TT. Substituting for N , K i , B, and Ψ, and changing to the ρ * density, we find an "effective" source, The compact moment is There are eight different field integrals, although some are very similar. They can also be evaluated using the methods of Sec. IV A 2. Adding everything up, we find the final result, which contains 1PN and 2PN terms, This reduces to WW (4.26) in the GR limit.
We can, of course, discard the delta function terms, which are non-TT. There is only one type of field integral to calculate. The final five-index moment is entirely 1.5PN order, This reduces to WW (4.27a) in the GR limit.

EW
With two time derivatives, we only need the source to O(ρǫ), (4.53). There is again only one field integral to calculate, although it is the most difficult by far. The final moment is entirely 2PN order, (4.55) This reduces to WW (4.27b) in the GR limit.

V. TWO-BODY EW MOMENTS
So far we have written down expressions for the EW moments which are fully general and can be applied to a system with any number of compact objects. Now we convert the moments to the two-body case relevant for compact binaries. We define the masses m A of the bodies as m 1 and m 2 . They have positions x 1 and x 2 and velocities v 1 and v 2 .
It will be useful to have the moments expressed in terms of relative variables, r = r 12 , x = x 1 − x 2 ,n =n 12 = x/r, and v = v 1 − v 2 . (Note that we no longer usen to represent unit normals at the surface of M.) In general relativity, the relationship between the individual variables and the relative variables can be fixed by writing down the conserved linear momentum. The loss of momentum to gravitational waves occurs at 2.5PN order, beyond what we need to worry about. However, in scalar-tensor theory, dipole radiation reaction enters at 1.5PN order, so it must be taken into account. From MW (6.9) and (6.10), we find that the individual and relative velocities are related by Here is a scalar-tensor parameter that enters the equations of motion at Newtonian order. We also have defined Finally, we have the usual variables m ≡ m 1 + m 2 , µ ≡ m 1 m 2 /m, η ≡ µ/m, and δm ≡ m 1 − m 2 . The first term in (5.2) is of relative 1PN order and the second term, representing the dipole radiation, is of relative 1.5PN order. We can find an antiderivative and write down the corresponding position relations, It may seem that we would need the O(ǫ 2 ) terms in δ i and δ i x in order to calculate 2PN gravitational waves. As it turns out, the 2PN piece of δ i is not needed anywhere, while the 2PN piece of δ i x is only needed in the first (0PN) term of I ij EW , where it cancels exactly. To simplify our expressions, we introduce other scalar-tensor parameters from MW, We also use the notation where ξ is one ofβ,δ, orχ. (This notation should not be confused with S + and S − .) Note that theχ parameters do not occur in the EW moments; however, they do appear in the equations of motion, which we will need shortly.
Since acceleration terms only appear at 1.5PN (in I ij EW ) or 2PN order (in I ij EW , I ijk EW , and I ijkl EW ), we can replace them with the Newtonian expressions, This will bring the expressions in line with WW (6.6) and allow an easy comparison. Finally, we evaluate (4.42) to get the total three-index moment. The result is a set of simplified, two-body EW moments expressed in relative coordinates, Gαmrn ij x (k δ lm) , (5.10d) Gαmrx ij δ (kl δ mn) .
For the most part, these moments have the same basic forms as those in WW, with the coefficients altered to include the scalar-tensor parameters (5.3), (5.4), and (5.7). As we saw in Sec. IV, both I ij EW and I ijk EW have new terms at 1.5PN and 2PN order, respectively, resulting from the nonvanishing surface moments. The conversion to relative coordinates introduces another new term in I ijk EW at 2PN order. It arises because of the interaction between the lowest order (0.5PN) piece of I ijk EW and the 1.5PN (dipole radiation) terms in (5.2) and (5.6). No such new term is created in the two-index moment: The 1.5PN piece of (5.6) cancels out just like the (unshown) 2PN piece.
The final step is to take time derivatives. The equation forh ij , (2.21), shows that we must take two derivatives of each moment. The five-and six-index moments also contain their own time derivatives. Along the way, we need to substitute the relative equation of motion for each acceleration a i . We take this result from MW, keeping terms up to 2PN order, While the Newtonian, 1PN, and 2PN terms are modified from the GR equations of motion, the 1.5PN terms, caused by dipole radiation reaction, are entirely new. They will cause the introduction of more new terms in the waveform, arising from the two-index EW moment (at 1.5PN order) and the three-index EW moment (at 2PN order). We hold off on presenting the final results until Sec. VII.

VI. RADIATION-ZONE INTEGRALS
So far, our calculation of the gravitational waveform has only considered the contribution from the near zone. We must also calculate the contribution from the radiation zone. To do so, we will need to evaluate (2.24a) forh ij , dropping all terms which fall off faster than 1/R. The first step is to derive an expression for the source τ ij in the radiation zone. Since there are no compact sources at R > R, τ ij is composed purely of field terms. The fields N =h 00 , K i =h 0i , B ij =h ij , and Ψ = ϕ − 1 themselves can be found by summing near-zone and radiation-zone contributions.
A. Radiation-zone fields: Near-zone contributions The near-zone contributions can be found using (2.17a) and (2.17b). Recall that the Epstein-Wagoner construction (2.21) is a special case of (2.17a) forh ij in the far-away zone. Now we want expressions for all of the fields at arbitrary R.
We need to evaluate the various moments defined in (2.18a) and (2.18b). We begin with the moments of τ 00 . We use the expression in (4.2), but we only need it to O(ρǫ). [Actually, all the sources need to be evaluated to O(ρǫ 3/2 ), but only τ 0i has terms at that order.] The integration over the field terms uses the same strategies we used in calculating the EW moments. We find In the second step, we rewrite the moment in terms of m ≡ A m A , the total mass of the system, and E, the (lowest order) conserved energy. [See MW (6.4).] The rest of the moments M 00Q can be renamed I Q , following the definition in MW (3.7a). We have which, to the required order, is equal to zero in the CM frame. (This is the argument used earlier to reduce the order of A m A x i A in surface integrals.) Finally, the next few integrals are only needed to lowest order, The τ 0i moments require the source to O(ρǫ 3/2 ), (4.45). We find This is proportional to the total momentum [see MW (6.5)], and we can set it equal to zero. The other moments are only needed to lowest order, where ǫ ijk is the totally antisymmetric Levi-Civita symbol (ǫ 123 = +1). The current moments J iQ are defined in MW (3.7b) as Note that, unlike the exact equivalence of M 00Q and I Q , the equality between M 0Q , I Q , and J Q−1 is valid only to lowest order in the post-Newtonian expansion.
The τ ij moments require the source to O(ρǫ), (4.53). We need them only to leading order, Here again, equality between M Q , I Q , and (potentially) J Q−1 only holds to the lowest post-Newtonian order. Finally, the scalar moments require τ s to O(ρǫ). Taken from MW (4.9e) and converted to the ρ * density, it is given by (6.8) The lowest order moment is where we have defined a "scalar mass" weighted by sensitivity factors. We also let m s1 be the rest of M s (when calculated to 1PN order). The rest of the M Q s can be renamed I Q s , following MW (3.7d). To the orders we need, they are All of these moments can be plugged back into (2.17a) and (2.17b) to produce the radiation-zone fields, as calculated from near-zone integrals, It is worth pointing out the relative post-Newtonian orders of the terms in these expressions. The lowest order terms, which will serve as the reference, are those involving m and m s . Relative to those, the other terms in N N enter at 1PN (those containing E and I ij ) and 1.5PN (I ijk ) order. The first term in K i N and B ij N is 1PN, while the second in each is 1.5PN beyond the leading order. For Ψ N , the m s1 piece is 1PN. The dipole term has contributions at both 0.5PN and 1.5PN order. The quadrupole piece is relative 1PN order, and the octupole is 1.5PN.

B. Radiation-zone fields: Radiation-zone contributions
We must also calculate the radiation-zone contributions to the radiation-zone fields using (2.24a) and (2.24b). In the radiation zone, there are no compact sources, so τ µν and τ s are made up purely of field terms. For τ µν , we can use MW (3.4) and (3.5), which write Λ µν and Λ µν s generically in terms of N , K i , B ij , B, and Ψ. For these, we substitute (6.12a)-(6.12d).
We start withh 00 . To lowest order, we find This source has the generic form with l = 0 and n = 4. If we restrict ourselves to sources of this form, (2.24a) can be rewritten as where A(s, R) ≡ R+s R P l (ξ) p n−1 dp , (6.16a) B(s, R) ≡ R+s s P l (ξ) p n−1 dp , (6.16b) and ξ ≡ R + 2s In the l = 0 case, this is easy, and we find These terms are 1.5PN beyond the leading order term in (6.12a). Since we kept no higher order terms in that expression, we can also stop here. The total field N =h 00 in the radiation zone is just the sum of (6.12a) and (6.18). The second term in (6.18) is obviously new to scalar-tensor theory. As we shall see, it turns out not to have an impact on the final gravitational waves. Otherwise, N in the radiation zone has exactly the same form as in GR, modified only by the addition of the factors G(1 − ζ), both explicitly and in the definitions of I Q . It turns out that there are no relevant radiation-zone contributions to K i ; to the order we need, (6.12b) is the complete expression in the radiation zone. To calculate B ij C−N , we use We use the same method as for N C−N , except that here we have two different (l, n) combinations, (2,4) and (0, 4), and thus two integrals (6.15) to evaluate. Doing so gives These terms are added to (6.12c) to produce the final expression for B ij . They enter at 1.5PN beyond the leading order. Finally, for the scalar, we need τ s to lowest order. It is given by For a source of the form (6.14), (2.24b) reduces to with A(s, R) and B(s, R) defined just as in (6.16a) and (6.16b). The source has (l, m) = (4, 0), and the result is which adds to Ψ N , (6.12d), at 1.5PN order.

C. Radiation-zone contributions to the GWs
Now that we have the fields in the radiation zone, we need to evaluate the gravitational waveform. Fundamentally, this just means continuing the procedure for calculatingh ij in Sec. VI B to higher order. While doing so, we drop non-TT terms and terms that drop off faster than 1/R. For simplicity, we will treat contributions from Λ ij and Λ ij s separately. We will also do the calculation order by order. We must be careful doing so: The expressions for Λ ij and Λ ij s in MW are sorted by post-Newtonian order in the near zone. In the radiation zone, the ordering can be slightly different. For instance, while K i is at 0.5PN order relative to N in the near zone, it is at 1PN relative order in the radiation zone. Furthermore, time derivatives in the radiation zone do not increase the post-Newtonian order relative to spatial derivatives.
The lowest order term in Λ ij is proportional to N ,i N ,j . (Now that we are dealing with GWs, we can ignore the non-TT δ ij term.) We first plug in the 0PN monopole pieces of N . As seen in Sec. VI B, the resultingh ij scales like 1/R 2 and, as such, dies off in the far-away zone. At the next order, N ,i N ,j generates cross terms between the 0PN piece of N and the 1PN pieces (including both the energy E and the mass quadrupole I ij ). Other terms which contribute at the same order are N ,(iK j) , N ,(i B ,j) , and NB ij , all featuring the 0PN piece of N and the 1PN piece of the other field. The last term is found in (4.4c) of [41]; its two time derivatives place it at higher order in the near-zone counting scheme.
Our expressions for the fields have several unevaluated spatial derivatives; the expression for Λ ij adds more. When we evaluate them, we must remember that the moments are functions of retarded time, τ = t−R, so that, for instance, ∂ c I ab = −İ abN c . Completing all the derivatives and converting products ofN i to STF products, we find (6.24) Here we have ignored δ ij terms, which are non-TT. Careful examination shows that we also should ignore terms proportional toN ij , which will only produce non-TT terms inh ij . This eliminates the term with monopole-monopole coupling (i.e., the one proportional to mE) and leaves only those with monopole-mass quadrupole couplings. In the end, we have to evaluate (6.15) for l = 4, n = 3-6; l = 2, n = 3-6; and l = 0, n = 2-4. Adding everything up, changing fromN L back toN L , and discarding more non-TT terms which arise along the way, we find We have made a slight change of variable from (6.15): s → s/2. This is done to match the first term in WW (5.8); we see that the two expressions are equal except for a factor G 2 (1 − ζ) 2 in the scalar-tensor case. [Recall that the moment I ij contains the factor G(1 − ζ).] Later, we will also bring the first term under the integral, as in WW, but for now we leave it separate. This is to emphasize the difference in the two pieces: The first relies on the instantaneous (but retarded) value of ... I ij , while the second requires a weighted integral of (4) I ij over the entire past history of the source. It is the lowest order "tail" term.
The expression (6.25) turns out to be 1.5PN order beyond the Newtonian quadrupole. Therefore, we only need to go one-half order higher to find the 2PN waveform. In this case, the source will include contributions from the same four terms in Λ ij as before, but they will now be cross terms between the 0PN piece of N and the 1.5PN pieces of N , K i , B, and B ij . Note that all of the 1.5PN pieces coming from radiation-zone integrals (those proportional to m 2 or m 2 s ) will generate only non-TT terms inh ij . Similarly, the term N N ,i N ,j in Λ ij should contribute at this order, but any piece ofh ij it generates will also be non-TT. In the end, the relevant piece of the source will contain only monopole-mass octupole and monopole-current quadrupole couplings. After a great deal of algebra, we find it to be J c(i ǫ j)acN a + · · · . (6.26) We must evaluate (6.15) for l = 5, n = 3-7; l = 3, n = 3-7; and l = 1, n = 2-5. The final result for the GWs is These terms are again identical to the GR results [the remainder of WW (5.8)], except for a factor of G 2 (1 − ζ) 2 . Now we examine the waveform produced by the Λ ij s piece of τ ij . We will mainly be concerned with the term in Λ ij s which is proportional to Ψ ,i Ψ ,j . The lowest order contribution from this term involves the 0PN monopole pieces of Ψ. The resultingh ij was calculated in Sec. VI B and shown to scale like 1/R 2 , making it irrelevant in the far-away zone. Because Ψ contains nonvanishing dipole terms, the next highest order source is only 0.5PN beyond the leading order. It consists of monopole-dipole couplings, However, each piece of this source will generate a term inh ij proportional toN i ,N j , or δ ij , all of which are non-TT. Moving on to the next order, there are three types of terms in τ ij : a monopole-monopole coupling (the cross term of m s and m s1 ), monopole-quadrupole couplings, and dipole-dipole couplings. The first generates only non-TT terms inh ij . The second is analogous to the monopole-quadrupole terms in Λ ij , while the third has no counterpart there.
The total source at this order (without non-TT terms) can be written as (6.29) This requires integrals with l = 4, n = 2-6; l = 2, n = 2-6; and l = 0, n = 2-6. Many of these can be reused from earlier calculations; however, some are brand new. The final expression for the GWs is (Recall that the scalar dipole moment I i s includes a factor of ζ, which will cancel the 1/ζ in front of the second term.) This is an entirely new contribution to the GWs at 1.5PN order, featuring a new type of hereditary integral without a logarithmic factor. We will discuss it further in Sec. VII B.
At the next order, the ΨΨ ,i Ψ ,j term in Λ ij s begins to contribute to the source. However, like the N N ,i N ,j piece of Λ ij , it produces only non-TT terms inh ij and can be ignored. In the Ψ ,i Ψ ,j term, we have monopole-dipole couplings (0PN monopole-1.5PN dipole and 1PN monopole-0.5PN dipole), a triple-monopole coupling (from Ψ C−N ), monopole-octupole couplings, and dipole-quadrupole couplings. The first and second of these produce only non-TT terms in the final waveform. The other two types of terms give This requires integrals with l = 5, n = 2-7; l = 3, n = 2-7; and l = 1, n = 2-7. The final answer is This is a new piece of the waveform at 2PN order.

A. Final tensor waveform
To find the final tensor waveform, we add the contributions from the near and radiation zones. The near-zone contribution is found by inserting the differentiated two-body Epstein-Wagoner moments into (2.21). For the radiationzone pieces (6.25), (6.27), (6.30), and (6.32), we bring the instantaneous terms inside the integrals. Then we write the moments I Q , J Q , I Q s , and their derivatives explicitly in terms of relative two-body variables. Finally, we sum the four pieces to find the complete radiation-zone contribution.
The final tensor waveform can be written as a post-Newtonian expansion, where the superscripts on P denote the PN order of each term. For clarity, we have separated out the 1.5PN and 2PN near-zone terms from the radiation-zone terms at the same order. We find should effectively enter the waveform at 0PN order, with the quadrupole-quadrupole memory at higher order. (The exact PN ordering of the memory terms will depend on the specifics of the scalar-tensor theory and the compact object sensitivities.) Scalar-tensor theory also produces instantaneous dipole-dipole terms at 1.5PN order. Like the hereditary term, they contain DC components; however, since they are not integrated over the binary's history, their effect remains at 1.5PN order. They are equivalent to the nonhereditary, zero-frequency terms Arun et al. discovered at 2.5PN order in general relativity [61]. By contrast, the new 2PN hereditary term does not contain a memory effect. While it has the same basic form as the 1.5PN term (i.e., no logarithm in the integrand), the beating between the mass dipole and mass quadrupole produces no DC component. The same is true for the new instantaneous terms with this coupling. This result is equivalent to the lack of a 0.5PN memory effect (appearing formally at 3PN order) in general relativity: In that case, the beating is between the mass quadrupole and mass octupole, resulting in no zero-frequency terms [6].
It is instructive to examine the waveform in a few special cases. For binary black holes, s 1 = s 2 = 1/2, and all sensitivity derivatives vanish. This means that α = 1 − ζ, and all of the rest of our scalar-tensor parameters (γ,β i , δ i ,χ i , S + , and S − ) vanish. The 2.5PN equations of motion derived in MW then have the exact same form as in general relativity, except for a factor 1 − ζ multiplying each instance of the total mass m. That is, the equations of motion for a binary with masses (m 1 , m 2 ) in general relativity are identical to those for a binary with rescaled masses (m 1 /(1 − ζ), m 2 /(1 − ζ)) in scalar-tensor theory. Since the masses of the bodies are defined by their Keplerian motion, this rescaling is unmeasurable. Therefore, to 2.5PN order, the motion of two black holes in scalar-tensor theory is indistinguishable from the motion in GR. MW predicted that the gravitational waves produced by binary black holes would be similarly indistinguishable from those produced in GR. We see here that the conjecture is correct, at least to 2PN order in the tensor gravitational waves.
As discussed in MW, this result is not surprising. Hawking originally showed that stationary, asymptotically flat black holes in vacuum are identical in both theories [62], leading to conjectures that the same might be true for black hole binaries. Still, this work shows only that the theories are indistinguishable to 2.5PN order in the dynamics and 2PN order in the radiation. It remains for future work to investigate whether indistinguishability holds to all post-Newtonian orders. (See MW for precise details of a conjecture on this point.) There is good evidence that it does: Yunes et al. [27] proved it, but only to lowest order in the mass ratio. Healy et al. [63] used numerical relativity simulations to show that any initial scalar field in the system is quickly radiated away, after which the holes behave identically to the general relativity case. Possible caveats which may break indistinguishability include the introduction of a potential for the scalar field or a time-varying scalar field at infinity [64,65]. It would be interesting to investigate the dynamics and radiation in such scenarios.
For a system containing one neutron star (say, body 1) and one black hole (body 2), s 2 = 1/2, s ′ 2 = 0, and s ′′ 2 = 0. Then α = 1 − ζ, and all other parameters vanish except For the opposite choice of bodies, S − has the opposite sign. Notably, S 2 + = S 2 − = S + S − = Q/(4ζ). Through 1PN order, only α,γ = 0, andβ i = 0 appear in the expressions for the waves, and so the waves are identical to those in GR (after mass rescaling). At 1.5PN order, deviations start to occur. However, the deviations are always parametrized byδ 1 , ζS 2 + , ζS 2 − , or ζS + S − . Therefore, through 2PN order, the tensor waveform for a mixed black hole-neutron star system differs from the general relativity waveform only by the single parameter Q. This is again equivalent to a result found by MW for the equations of motion. Because Q contains no information on the derivatives of the coupling function ω(φ) (i.e., the parameters λ 1 and λ 2 ), we cannot, at 2PN order, formally distinguish the waveform produced in the Brans-Dicke theory [ω(φ) = ω 0 ] from that produced in a general scalar-tensor theory of the type we consider. The only difference will be that for a given neutron star of a certain central density and total number of baryons, different scalar-tensor theories will produce different results for the neutron star mass m 1 and sensitivity s 1 . One can imagine using gravitational waves to measure masses and sensitivities for a wide variety of sources and then producing a mass-sensitivity relation, much like the neutron-star mass-radius relations used to study the nuclear equation of state. This relation could then be used to rule out various models of the coupling function ω(φ).
Before we can completely understand what the measurement of gravitational waves from a compact binary will tell us about scalar-tensor theories of gravity, we must first derive the gravitational-wave phasing. To do so, we will need the rate at which the binary loses energy to gravitational waves, both tensor and scalar. The next paper in our series will derive the scalar waves. The process is identical to that presented in this paper, with two complications. First, we cannot eliminate non-TT terms in the scalar case; indeed, for a scalar, TT is not defined. This will lead to more surviving terms and a need for care when reusing parts of this analysis. Second, and more daunting, the scalar "EW moments" defined by (2.19b) begin with a monopole moment, which has relative order -1PN compared to the tensor two-index moment. The next piece, a dipole moment, has order -0.5PN. Therefore, to obtain the 2PN scalar waveform, we will need to compute the source τ s (τ, x ′ ) to O(ρǫ 3 ), or 3PN order.