Gravitational Radiation from Close Binaries with Time-Varying Masses
A. Miguel Holgado Paul M. Ricker



















































A. Miguel Holgado
holgado2@illinois.edu

0000-0003-4143-8132]A. Miguel Holgado
0000-0002-5294-0630]Paul M. Ricker
Department of Astronomy and National Center for Supercomputing Applications, University of Illinois at Urbana-Champaign, Urbana IL, 61801, USA
















In the quadrupole approximation of General Relativity in the weak-field limit, a time-varying quadrupole moment generates gravitational radiation. 
Binary orbits are one of the main mechanisms for producing gravitational waves and are the main sources and backgrounds for gravitational-wave detectors across the multi-band spectrum. 
In this Paper, we introduce additional contributions to the gravitational radiation from close binaries that arise from time-varying masses along with those produced by orbital motion.
We derive phase-dependent formulae for these effects in the quadrupolar limit for binary point masses, which reduce to the formulae that Peters and Mathews (1963) derived when the mass of each component is taken to be constant. 
We show that gravitational radiation from mass variation can be orders of magnitude greater than that of orbital motion. 









§ INTRODUCTION
 
Einstein's theory of General Relativity (GR) predicted the existence of gravitational waves (GWs) <cit.>. 
The existence of GWs was later indirectly confirmed from observations of the orbital decay of a pulsar binary <cit.>. 
Recently, the Laser Interferometer Gravitational Wave Observatory (LIGO) has ushered in the era of GW astronomy with detections of GWs from compact-object mergers <cit.>.
The Laser Interferometer Space Antenna (LISA) will then open up the millihertz frequency band <cit.>. 
The combination of ground-based and space-based detectors holds promise for advances in multi-band GW astronomy, standard-siren cosmology, and tests of GR <cit.>. 
At nanohertz frequencies, pulsar timing arrays are currently operating and will eventually detect the stochastic GW background from the cosmic population of inspiraling supermassive black-hole binaries <cit.>.

Even though a detection has yet to be made, the upper limits are providing astrophysically relevant constraints on the cosmic
supermassive black-hole binary population <cit.>.

The decihertz band may be opened with future space-based detectors. 
The main astrophysical GW sources and GW backgrounds that these detectors across the multi-band GW spectrum aim to detect are from the orbits of compact binaries. 
In this Paper, we consider an additional effect: time-varying masses in a binary due to wind mass loss or accretion from mass transfer, for example, generate time-varying quadrupole moments that generate gravitational radiation in addition to those produced by orbital motion. 

The analytic expressions for gravitational radiation from point-mass binaries in the quadrupole approximation of GR in the weak-field limit <cit.> have been used widely in the astrophysics community for modeling the evolution of close stellar binaries, generating compact-binary distributions from binary population synthesis codes <cit.>, and estimating the strain from sub-parsec supermassive black-hole binary candidates.
In this Paper, we consider the more general treatment of the two-body problem in which the mass of the binary is time-dependent. 
This is motivated by the variety of binary systems that interact, i.e., low-mass X-ray binaries, cataclysmic variables, and supernovae in close binaries.
Binary interactions include Roche-lobe mass transfer <cit.>, mass loss, and common-envelope evolution <cit.>. 


A number of previous works in the astrophysics and fundamental-physics communities have considered close binaries with time-varying masses in different contexts <cit.>. 
<cit.> considered the effect of the time-varying masses on the quadrupole moment in the context of neutron-star spindown, but did not consider the phase-dependent effects that we present here. 
We derive new expressions for gravitational radiation from binaries with variable masses and we quantify the differences between treatments with and without the additional contributions for an explosive mass-loss case.
The expressions that arise from this new analysis reduce to the Peters and Mathews (PM herein) formulae <cit.> when the mass of each component is taken to be constant.




§ THE PERTURBED TWO-BODY PROBLEM

If two constant point masses are in an orbit bound by their mutual gravitational interaction, their motion can be described as a Keplerian orbit in the absence of any internal and external perturbations. 
The six orbital elements (a,e,i,ϖ,ω,ν) fully define the Keplerian orbit, where a is the semi-major axis, e is the eccentricity, i is the inclination, ϖ is the longitude of the ascending node, ω is the argument of periapsis, and ν is the true anomaly.
We consider a gravitationally bound point-mass binary with time-dependent component masses m_1(t) and m_2(t) with a reduced mass μ = m_1 m_2/M(t) ≡μ(t), where M(t) =m_1+m_2 is the total mass.

Allowing for time-dependent masses introduces perturbative forces on the two bodies. 

We consider the case where the primary's mass remains constant and the secondary undergoes isotropic mass loss: ṁ_̇1̇ = 0 and ṁ_2 < 0. 
The resulting total mass-loss rate of the binary is Ṁ = ṁ_2. 
For the isotropic case, only the self-gravity between the two components governs the equation of motion for the relative separation vector r


    r̈ = - GM/r^3 r ,


where G is the gravitational constant, r is the magnitude of r, and where we take the same coordinate system as <cit.> with the origin located at the barycenter.
There is an implicit perturbing force that induces the following acceleration


    f = -1/2Ṁ/Mṙ ,


which is implicit in the isotropic variation of the mass-losing body. 
This has been derived in several previous works <cit.>. 
If the mass variation is anisotropic, i.e., either from non-conservative mass exchange and asymmetric ejecta, then the equation of motion will have an explicit perturbing force <cit.>. 
The orbital evolution can be described as an osculating orbit, which is the Keplerian orbit the binary would follow if the perturbation ceased immediately.
The full set of equations for the time evolution of the osculating orbital elements from the isotropically varying mass are <cit.>



    da/ dt   = - (1+2ecosν + e^2/1-e^2) Ṁ/M a  , 
     de/ dt   = - (e+cosν) Ṁ/M , 
     dω/ dt   = -sinν/eṀ/M , 
     dν/ dt   = [G M a (1-e^2)]^1/2(1+ecosν)^2/a^2 (1-e^2)^2 + sinν/eṀ/M ,



where the specific angular momentum is conserved, i.e., the quantity r^2 θ̇ = [G M a (1-e^2)]^1/2 is a constant of the motion. 
Since there is no net force in the direction perpendicular to the orbital plane, i̇ = ϖ̇ = 0, and the orbital motion thus remains coplanar. 
The mass variation causes the argument of periapsis to precess at a rate ω̇ (see Eq. <ref>).




The upcoming LISA mission will be observing mass-transferring binaries, including low-mass X-ray binaries, cataclysmic variables, ultra-compact X-ray binaries, etc. over a ∼5-yr baseline.
For a binary with an orbital period of 1   hr, the 5-yr observation baseline would correspond to ∼3×10^5   rad in phase, which is obtained by taking the product of the angular frequency of the binary and the observing duration. 
Investigating the importance of periapsis precession for mass-transferring binaries relevant to the LISA mission is left for future work. 
For now, we take ω̇≈ 0.




§ QUADRUPOLE GRAVITATIONAL RADIATION

For eccentric Keplerian orbits, the orbital separation r between the binary masses obeys


    r = a(1-e^2)/1+ecosν .


The components of the mass quadrupole tensor I_jk for an eccentric binary are



    I_xx   = μ r^2 cos^2 θ = μ a^2 (1-e^2)^2 cos^2θ/(1+ecosν)^2 , 
    
    I_yy   = μ r^2 sin^2θ = μ a^2 (1-e^2)^2 sin^2θ/(1+ecosν)^2 , 
    
    I_xy   = μ r^2 sinθcosθ = μ a^2 (1-e^2)^2 sinθcosθ/(1+ecosν)^2 ,

 

where θ = ν + ω, which simplifies to θ = ν if the precession of the argument of periapsis is taken to be negligible.
We can choose the initial θ to begin where ω = 0.  
For isotropic wind mass loss, the periapsis distance r_ p = a(1-e) is constant, so ṙ_ p = 0.
When including gravitational radiation reaction, the periapsis distance will indeed decrease, though one can take ṙ_ p≈ 0 on timescales less than the mass-loss timescale M/|Ṁ| when the GW coalescence timescale is much longer than the mass-loss timescale. 



The third time derivatives of the mass-quadrupole tensor components are



 

    ⃛I_xx   = 2 G^3/2μ M^3/2/a^5/2(1-e^2)^5/2(1+ecosθ)^2 (2sin2θ + 3 e sinθcos^2θ) + 6 G μ̇ M /a (1-e^2)(1-2cos^2θ - e cos^3 θ) 
       =-6 μ̈[G M a (1-e^2)]^1/2sinθcosθ/1+ecosθ + ⃛μ a^2 (1-e^2)^2 cos^2θ/(1+ecosθ)^2 , 
    ⃛I_yy   = - 2 G^3/2μ M^3/2/a^5/2(1-e^2)^5/2(1+ecosθ)^2 (2sin  2θ + esinθ(1+3cos^2θ)) + 6 G μ̇ M/a (1 - e^2)(e cos^3 θ + 2 cos^2 θ + ecosθ + e^2 - 1) 
       =+  6 μ̈[GMa(1-e^2)]^1/2esinθ +sinθcosθ/1+ecosθ + ⃛μ a^2 (1-e^2)^2 sin^2θ/(1+ecosθ)^2 , 
    ⃛I_xy   = -2 G^3/2μ M^3/2/a^5/2(1-e^2)^5/2 (1+ecosθ)^2(2cos2θ - ecosθ(1 - 3 cos^2θ)) - 6 G μ̇ M/a (1-e^2)sinθ(e + 2 cosθ + e cos^2θ) 
       = - 3 μ̈[G M a (1-e^2)]^1/21-2 cos^2θ - ecosθ/1+ecosθ + ⃛μ a^2 (1-e^2)^2 sinθcosθ/(1+ecosθ)^2 .






 

    ⃛I_xx   = - 2 G μ M/r^2ṙcos^2θ + 4 G μ M/rθ̇sin 2θ - 2 G μṀ/rcos^2θ - 6 G μ̇ M/rcos^2θ + 6 μ̇ r^2 θ̇^2 (e^2 sin^2νcos^2 θ/(1+ecosν)^2 - e sinνsin 2θ/1+ecosν + sin^2θ) 
       =+ 3 μ̈ r^2 θ̇(2 e sinνcos^2θ/1 + ecosν - sin 2θ) + ⃛μ r^2 cos^2θ , 
    ⃛I_yy   = - 2Gμ M/r^2ṙsin^2θ - 4 G μ M/rθ̇sin 2θ - 2 G μṀ/rsin^2θ - 6 G μ̇ M/rsin^2 θ + 2 μ̇ r^2 θ̇^2 (e^2 sin^2 νsin^2θ/(1+ecosν)^2 + esinνsin 2θ/1+ecosν + cos^2θ) 
       =+ 3 μ̈ r^2 θ̇(2 esinνsin^2θ/1+ecosν - sin 2θ) + ⃛μ r^2 sin^2 θ , 
    ⃛I_xy   = - Gμ M/r^2ṙsin 2θ - 4 G μ M/rθ̇cos 2θ - G μṀ/rsin 2θ - 3 G μ̇ M/rsin 2θ + 3 μ r^2 θ̇(e^2 sin^2 νsin 2θ/(1+ecosν)^2 + 2 e sinνcos 2θ/1+ecosν - sin 2θ) 
       =+ 3 μ̈ r^2 θ̇(e sinνsin 2θ/1+ecosν + cos 2θ) + 1/2⃛μ r^2 sin 2θ .




If the periapsis precession is taken to be negligible, ω̇ = 0, then the first two terms of Eqs. (<ref>) simplify to the canonical PM expressions. 
The additional terms in Eqs. (<ref>) contain time derivatives of the masses of the binary, which may also be thought of as time derivatives of the binary's chirp mass M = (m_1 m_2)^3/5/ M^1/5. 




Gravitational radiation carries away the orbital energy and angular momentum of the binary at the following rates



    .  dE/ dt|_ GW   = 1/5G/c^5⟨⃛ I_jk⃛ I_jk⟩ , 
    .  dJ_i/ dt|_ GW   = 2/5G/c^5ϵ_ijk⟨Ï_jm⃛ I_km⟩ ,



where c is the speed of light and I_jk is the reduced (trace-free) quadrupole moment tensor I_jk = I_jk - 1/3δ_jkδ_ℓ m I_ℓ m.


The secular evolution equations for the semi-major axis and eccentricity are



    da/ dt   = … , 
     de/ dt   = … .




If we take the orbit to be circular, then the GW luminosity simplifies to


    L_ GW = 32/5G^4/c^5μ^2 M^3/a^5 + 72/5G^3/c^5μ̇^2 M^2/a^2 - 48/5G^3/c^5μ̈μ M^2 /a^2
     
    - 12/5G^2/c^5⃛μμ̇ M a + 18/5G^2/c^5μ̈^2 M a + 1/5G/c^5⃛μ^2 a^4  ,


where the positive terms in Eq. (<ref>) correspond to GW energy loss from the orbit and the negative terms correspond to GW energy gain in the orbit.
From the strain tensor in the transverse-traceless gauge,


    h_jk^ TT = 2 G/c^4 D_ LÏ_jk ,


where D_ L is the luminosity distance to the source, the strain amplitude for a circular binary is


    h_0 = 4 G^5/3/c^4μ M^2/3Ω^2/3/D_ L + 2 G^5/3/c^4μ̇ M^2/3/D_ LΩ^1/3 + G^5/3/c^4μ̈ M^2/3/D_ LΩ^4/3 ,


where Ω is the angular frequency of the binary, Ω = 2 π f_ orb = 2π / P_ orb.
In the transverse-traceless gauge, the plus and cross polarizations of the strain are h_+ = h_xx^ TT = - h_yy^ TT and h_× = h_xy^ TT = h_yx^ TT.

From Eq. (<ref>), we can define the following parameters



    ζ_1 ≡μ̇/μ/2Ω , 
    ζ_2 ≡μ̈/μ/(2Ω)^2 ,



which quantify the relative importance of mass-variation effects on the GW amplitude relative to the PM term of a given source.
For ζ_1 = 1, the mass-transfer timescale needs to be comparable to the orbital period.
The true trajectory may differ significantly from that of an osculating orbit since the perturbing force in this regime is comparable to the self-gravity. 

We emphasize that even though the above analysis is solely applicable to a binary whose mass is varying isotropically, it is illustrative of how one includes time-varying masses in the PM analysis and eventually extensions to the post-Newtonian analysis. 
In the limit that the time-derivatives of the masses are zero, then the canonical PM expressions remain. 



§ EXPLOSIVE MASS LOSS IN A CLOSE BINARY

Supernova explosions in close binaries result in extreme mass loss. 
If enough mass is ejected during the explosion, then the binary can be disrupted, i.e., the final energy of the binary is ≥ 0 and the orbit becomes unbound. 
We again consider the case of isotropic mass loss such that the kick imparted to the remnant of the exploding body is solely due to the acceleration from the ejected mass (i.e., a “Blaauw kick” <cit.>), such that the equation of motion remains the same as Eq. (<ref>). 
In this regime, we choose to work in Cartesian coordinates such that r = (x,y,z), ṙ = (ẋ,ẏ,ż), and r̈ = (ẍ, ÿ, z̈).
We also ignore any interaction between the companion and the ejecta. 
When the mass-loss timescale is shorter than the orbital period, the perturbing force can become greater than the gravitational attractive force and the motion can no longer be accurately described as an osculating Keplerian orbit.
One then needs to integrate the equation of motion directly. 
The mass quadrupole tensor of the binary in these coordinates is  


    I_jk =
    [ μ x^2  μ xy  μ xz;  μ xy μ y^2  μ yz;  μ xz  μ yz μ z^2 ] ,


and the third derivatives of each component are


 

    ⃛I_xx   = ⃛μ x^2 + 6 μ̈ x ẋ + 6 μ̇ẋ^2 + 6 μ̇ x ẍ + 6 μẋẍ + 2 μ x ⃛x , 
    ⃛I_yy   = ⃛μ y^2 + 6 μ̈ y ẏ + 6 μ̇ẏ^2 + 6 μ̇ y ÿ + 6 μẏÿ + 2 μ y ⃛y , 
    ⃛I_zz   = ⃛μ z^2 + 6 μ̈ z ż + 6 μ̇ż^2 + 6 μ̇ z z̈ + 6 μżz̈ + 2 μ z ⃛z , 
    ⃛I_xy   = ⃛I_yx = ⃛μ xy + 3 μ̈ẋ y + 3 μ̇ẍ y + μ⃛x y + 3 μ̈ x ẏ
       yx = = + 6 μ̇ẋẏ + 3 μẍẏ + 3 μ̇ x ÿ + 3 μẋÿ + μ x ⃛y , 
    ⃛I_xz   = ⃛I_zx = ⃛μ xz + 3 μ̈ẋ z + 3 μ̇ẍ z + μ⃛x z + 3 μ̈ x ż
       zx = = + 6 μ̇ẋż + 3 μẍż + 3 μ̇ x z̈ + 3 μẋz̈ + μ x ⃛z , 
    ⃛I_yz   = ⃛I_zy = ⃛μ yz + 3 μ̈ẏ z + 3 μ̇ÿ z + μ⃛y z + 3 μ̈ y ż
       zy = = + 6 μ̇ẏż + 3 μÿż + 3 μ̇ y z̈ + 3 μẏz̈ + μ y ⃛z .




The GW luminosity and strain waveform can then be calculated with these expressions using Eqs. (<ref>) and (<ref>), respectively. 
Eqs. (<ref>) are general for any point-mass binary with a given equation of motion and time-dependent masses.
These equations can also be translated into any choice of coordinate system. 
Since there are no forces acting in the z direction, then the z components will be zero and the orbital motion will remain coplanar in the x-y plane. 
In binary population synthesis codes, supernova mass loss in close binaries is often treated as being effectively instantaneous and conservative of the total orbital energy and orbital angular momentum. 
When treating the gravitational radiation, the mass-loss timescale should be bounded by the light-crossing timescale of the binary.
The total energy and total angular momentum can be treated as conserved as long as the amount of energy and angular momentum radiated away in GWs is included in the balance equations.
Supernovae in close binaries are the most extreme cases of mass loss and are thus expected to be the loudest sources of GWs from mass variation.



 §.§ Numerical Comparison

We consider a close binary consisting of a neutron star and a helium star in an initially circular orbit with an initial semi-major axis of a_0 = 10^-2  au. 
We fix the neutron-star mass to be m_ NS,1 = 1.4 M_⊙ and the helium-star mass to be m_ He = 2.4 M_⊙. 
The helium star may be large enough to be tidally deformed, but we ignore this effect and are solely interested in the GWs sourced by the supernova mass loss since this process operates on dynamical timescales, while the tidal effects on the orbit operate on secular timescales.
The helium star will undergo a supernova explosion, which we take to be isotropic for simplicity and will form a neutron star with the same mass as the primary m_ NS,2 = m_ NS,1 = 1.4 M_⊙. 
The supernova explosion itself will be a source of GWs in the LIGO band, but here, we solely focus on the GWs sourced by the binary mass loss, whose characteristic frequency corresponds to the mass-loss timescale.
We choose the following smooth function as a model for the mass-losing body such that the time derivatives of the mass are continuous:


    m_2(t) = m_2, f + 1/2(m_2, i - m_2,  f) (1 + 2/πtan^-1[(t-t_0)/τ]) ,


where m_2, i is the initial mass, m_2, f is the final mass, τ is a given timescale for the mass loss and t_0 is the the time at which the mass loss peaks. 
Given that the only constraint for τ is the light-crossing time τ_ c of the binary, we consider two cases: τ = [τ_ c, 100 τ_ c]. 
The light-crossing time relative to the orbital period for this particular system is Ωτ_ c≈ 0.002, so the mass-loss timescale for both cases are less than the orbital period. 
Other smooth functions may be used to model the sharp decrease in mass, but such choices will not affect our overall results. 
We plot the orbital trajectory for the two cases in the top panel Figure <ref>. 
For each case, we integrate the equation of motion (see Eq. <ref>) with the mass of the secondary obeying Eq. (<ref>).
The peak of the mass loss is taken to occur at t_0 = P_0, where P_0 is the initial orbital period of the binary (obtained from Kepler's Third Law). 
The GW luminosity evaluated from Eqs. (<ref>) and (<ref>) and the GW luminosity predicted from the PM formalism [i.e., excluding terms in Eqs. (<ref>) with O (≥μ̇)] are plotted in the bottom panel of Figure <ref>. 



A spike in the GW luminosity occurs during the explosive mass loss and the amplitude of the spike increases with a decreasing mass-loss timescale.
The dashed lines in the bottom panel of Figure <ref> correspond to the GW luminosity computed with the PM formalism, i.e., ignoring the time-derivatives of the reduced mass in Eqs. (<ref>). 
The maximum difference between our predicted GW luminosity and that of PM is about 4.18 and 1.37 for the τ_ c and 100 τ_ c, respectively, in units of log_10(L_ GW/L_ GW,0).
Even for the longer mass-loss timescale case, our predicted GW luminosity is about an order of magnitude larger than the PM prediction.  
The blue dashed line, corresponding to the PM case with τ = τ_ c exhibits a peak which comes from the acceleration on the binary from the mass loss, i.e., the Blauuw kick.
Since the mass loss is isotropic, the final orbital configuration is in the same orbital plane as the initial one. 
A more systematic and comprehensive survey of the parameter space is left for future study, especially since additional physics such as natal kicks from asymmetric supernovae will affect the GW luminosity and final state of the binary <cit.>.



§ DISCUSSION AND CONCLUSIONS
 
In this Paper, we have shown that time-varying masses in close binaries contribute additional gravitational radiation along with the radiation produced by orbital motion. 
We have also shown that gravitational radiation from time-varying masses can be orders of magnitude greater than that of the orbital motion.
This analysis is performed in the quadrupolar approximation of GR in the weak-field limit and the derived formulae reduce to the PM formulae when the mass of each binary component is taken to be constant. 
Our analysis may be useful for implementations in binary population synthesis and generating approximate waveforms.
It may also provide new testable predictions for the evolution of compact binaries. 
Applying the analysis of this Paper to higher-order post-Newtonian orders for higher-compactness binary systems and to the gravitational self-force for high to extreme mass-ratio binaries is the next extension of this work.
These extensions may reveal richer physical behavior than the quadrupole approximation that is presented here. 
Future work will also include a more comprehensive survey of the parameter space including natal kicks from asymmetric supernovae. 
Additional future studies will investigate how these new contributions to the gravitational radiation will affect the formation of compact binaries and the rates and distributions of compact-object mergers. 
The analysis we present can also be extended to hierarchical triple systems.  
A detection of GWs from mass transfer would not only validate this effect, but may also provide new and complementary tests of GR at different regimes of curvature and field strengths. 




AMH is supported by the DOE NNSA Stewardship Science Graduate Fellowship under grant number DE-NA0003864.
AMH thanks Charalampos Markakis, Eliu Huerta, Gabrielle Allen, Xin Liu, Charles Gammie, Nico Yunes, Fred Rasio, and Vicky Kalogera for helpful comments and fruitful discussions.
PMR acknowledges support from NSF under AST 14-13367.




















 Matplotlib <cit.>, Numpy <cit.> .

yahapj


