Competitive effects between gravitational radiation and mass variation for two-body systems in circular orbits

This work provides, at lower order, general analytical solutions for the orbital separation, merging time, and orbital frequency of binary systems emitting gravitational waves while being submitted to mass variations. Specific features, depending on the exponent of the mass derivative, are investigated in details. Two phenomenologically interesting cases are explicitly considered : i) binaries formed by two light primordial black holes submitted to Hawking evaporation and ii) bodies driven by a Bondi accretion of phantom dark energy. It is shown that three different regimes arise, including an intricate non-monotonic behaviour of the system. We study subtle imprints that could be associated with those phenomena. A careful analysis of the conditions of validity of the different hypotheses performed is finally carried out.


I. INTRODUCTION
Newtonian orbits are stable for objects with constant masses.This is obviously not the case anymore if the considered masses become time-dependent, whatever the reason for this variation.This also becomes wrong in General Relativity (GR) when the system radiates energy through gravitational waves.The very existence of this energy has been intensively debated, even quite recently [1], as one might wonder if objects following geodesics -thus being "force-free" -can dissipate energy.It is however well settled since the works performed by Penrose, Bondi, and Sachs in the 1960s, although the conceptual arguments are somehow subtle and the technicalities quite involved [2].This work aims at investigating -at lowest order -the detailed evolution of the orbits of systems submitted both to the emission of gravitational waves (GWs) and to mass variations.
One of the most striking sources of GWs is the coalescence of black holes (BHs), which has been exhaustively studied theoretically and has proved to be experimentally fruitful since most of the GW signals detected so far by LIGO and Virgo collaborations come from this phenomenon (see, e.g., [3,4]).Although the most simple and straightforward cases are those of binaries with constant masses, inspiral binaries with time-varying masses and their consequences on the emitted gravitational radiation have also been studied in the literature [5][6][7][8][9][10][11][12][13][14][15].It is wellknown that the emission of GWs tend to make the two bodies of the binary come closer to each other, until they ultimately coalesce [16].However, the time-varying mass can have the converse effect, depending on whether the bodies gain or lose mass.If they gain mass, for instance through matter accretion, they inspiral, thus enhancing the effect of the emission of GWs.On the other hand, if they lose mass, they outspiral: in that case, mass loss produces an antagonistic effect to the one of gravitational radiation.
In this note, we focus on presenting a comprehensive view of the dynamics induced by the competitive effects of gravitational radiation and mass loss, using Newtonian analysis and treating the binary system as Keplerian, as well as the characteristic features of the imprint it produces on the emitted GWs.In section II, we review how one can theoretically deal with the coupling between mass variation and the backreaction produced by the emission of gravitational waves.However, unlike previous works -see Refs.[7][8][9][10][11][12] for accreting bodies or Refs.[13][14][15] for evaporating BHs in the context of braneworld models and extra-dimensions scenarioswe present a generic analytic solution of the differential equation satisfied by the orbital separation when the mass variation is taken into account.This leads in turn to analytic solutions for the frequency and strain of the emitted GWs.
For cosmology and phenomenology, we focus specifically on two situations that are driven by such a competition.In section III, we investigate the case of inspiral binaries of light primordial black holes (PBHs) submitted to Hawking evaporation.Due to their low mass, they are good candidates as sources of high-frequency GWs [17] in addition to standing as plausible dark matter candidates and valuable probes of the early universe.
The second case of interest is the one of binaries of BHs accreting phantom energy (i.e.violating energy dominance conditions) which makes their mass decrease [18].Although the accretion of a fluid onto a black hole is a long-standing problem in astrophysics [19], the study of this so-called phantom energy has proved particularly relevant in the context of dark energy [20].Assuming that the latter can be described by a perfect fluid of density ρ and pressure p, an equation of state parameter w = p/ρ smaller than -1 would make dark energy behave as phantom energy: accretion would decrease the mass.Coining down the value of w is a hot spot in modern cosmology [21] and doing so by examining small-scale systems like BHs binaries was proposed in [7,22] and questioned in [8].In section IV, we refine the analysis and draw a clear conclusion.
Beyond the specific examples chosen for their physical relevance, this study aims at offering an exhaustive classification of possible behaviors of binary systems emitting gravitational waves and submitted to mass variation effects (either accretion or radiation) described by a (positive or negative) power law.These generic results are presented in section V Eventually in section VI, we show how accounting for a time-varying mass changes the formal expression of the radiated gravitational power and in section VII we discuss how the results of section V are changed when considering non-identical masses presenting a high hierarchy.

II. COUPLING MASS VARIATION AND GWS EMISSION
The two-body problem with variable mass, for pointlike entities, was extensively studied in Ref. [23,24].When the variation of mass is isotropic, the equation of motion is given by where r = r 2 -r 1 , thus R ≡ |r| is the orbital separation between the two bodies of respective mass m 1 and m 2 , the total mass being m tot (t) ≡ m 1 (t)+m 2 (t).For orbital separations much larger than the Schwarzschild radii of the black holes, R ≫ R s = 2Gm/c 2 , we can describe the binary system using Keplerian dynamics.If µ ≡ m 1 m 2 /(m 1 +m 2 ) stands for the reduced mass, the orbital angular momentum of the system, for a circular orbit, whose generic expression is given by is conserved, even in the case of variable masses [23].If we consider that the two bodies have the same mass m(t), then the conservation of the orbital angular momentum Jorb = 0 provides showing that in the case of mass loss ( ṁ < 0) the two bodies outspiral, i.e. they drive away from each other, which is consistent with Refs.[14,15,23,25].From now on, we shall always assume, unless otherwise stated, that the two bodies have same mass m(t).This is the most interesting situation and it is phenomenologically sound -at least as a first approximation -if the mass spectrum of PBHs is nearly monochromatic (see [26] for a review of formation mechanisms).Only in section VII, will we drop this assumption in order to examine how certain results presented in the following are modified when the binary system is made of bodies of different masses, but presenting a strong hierarchy (typically one mass dominates over the other).The mass loss, in addition to separating the two bodies also induces a change in the orbital energy E orb = −Gm 2 /(2R), which is the total energy of the system, sum of kinetic and potential energy of the orbit.For variable masses, the variation of orbital energy − dE orb dt must be equal to the power associated to mass loss P ml , thus leading to On the other hand, the emission of GWs costs energy which is taken from the orbital energy of the system, carried away at a rate [16] P gw (t) = 64 5 In fact, m having an explicit time-dependence, it affects the form of the power radiated by GWs thus modifying Eq. ( 5), as already noted in [8], through several corrective terms.We discuss this point in further detail in section V, and show how the choice of the form (5) for P gw is connected to the condition of circularity of the orbit.
To consider the concomitant effect of gravitational radiation and mass variation, one can use the conservation of energy and write [8] − dE orbit dt = P ml + P gw , with P ml and P gw respectively given by Eqs. ( 4) and (5).
Taking the derivative of the orbital energy leads to Ṙ = − 128 5 which is the general differential equation relevant for our problem.This agrees with the result of Ref. [8] -up to the corrected prefactor 3 in the last term.This deserves a brief specific discussion.Equation (3) holds only if the masses are identical.There are two other limit cases allowing for simple formulas (as can immediately be seen from the expression of the angular momentum).They correspond to a strong mass hierarchy with asymmetrical losses.We specifically investigate those cases in the last section of this work.We will show there that if the varying mass is the small one, the prefactor 2, inappropriately used in [8] indeed appears.Hence the probable origin of the mistake.
If the functions m(t) and ṁ(t) are explicitly known, this differential equation, although non-linear, can be integrated since it is a Bernoulli differential equation of the form y ′ + P (t)y = Q(t)y n with n = −3 whose general solution is known.In the following two sections, we shall focus on mass losses described by rates of the form ṁ ∝ −1/m 2 (see section III) and ṁ ∝ −m 2 (see section IV).More generic results shall be found in section V.
Throughout all this work, we call "merging" or "coalescence" the situation corresponding to the precise vanishing of the orbital separation (R = 0), and not to the "contact" of horizons (R = 4Gm/c 2 ).For most observables this makes nearly no difference at all, but it is important when considering the end of the process.

III. INSPIRAL BINARIES WITH EVAPORATING BHS
Because of their small masses, primordial black holes can be particularly sensitive to the Hawking evaporation process (see [27,28] for pioneering works and [26] for a recent review).Quantum field theory in curved spacetimes predicts that black holes evaporate with a temperature [29,30] for a BH of mass m, ℏ being the Planck constant, k B the Boltzmann constant, c the speed of light, and G the gravitational constant.Even though for typical solarmass BHs, Hawking evaporation is too weak to play any significant role in the dynamics of binaries, for two-body systems made of PBHs, it could produce non-negligible effects on the overall dynamics, and leave a possible imprint on the emitted gravitational waves.The Hawking process leads to a mass loss at the rate [27,28] where α H accounts for the degrees of freedom of each particle contributing to the evaporation.For simplicity, we shall assume that α H is a constant, which is a good approximation.In this case, the differential equation is separable and integrates into where m 0 represents the initial mass of the BH whereas t ev ≡ m 3 0 /(3α H ) corresponds to the time of evaporation, i.e. the typical time it takes for the BH to evaporate completely.The other time scale of the problem under scrutiny is the time of coalescence which accounts for the limited life duration of the binary system.However, contrary to t ev which only depends on constants fixed during the initial setup (namely the initial mass m 0 ), the time of coalescence has an explicit dependence on the mass m (and not only on the initial mass m 0 ).Otherwise stated, if t ev is not to be altered by the dynamics of the system, the time of coalescence will necessarily be affected by the mass loss.Consequently, we will call t cc the time of coalescence of the binary system if the two BHs were of constant mass m 0 , with an initial orbital separation R 0 , i.e. [16] t cc = 5 512 while t coal will denote the real time of coalescence, that is to say the one taking into account the mass loss, which is likely to be different from t cc and yet to be determined.
A. Evolution of the orbital separation Specifying (7) to the case of Hawking evaporation provides Ṙ = − 128 5 Evaluating the above at an initial time If one assumes that the system is prepared -which boils down to giving ourselves m 0 and an initial separation R 0 -such that t ev < 4t cc , then the evaporation process initially dominates.Since the evaporation increases R and since the GW part is in m 3 /R 3 -see Eq. ( 12) -when m diminishes and R increases, both these variables contribute to the dwindling of this term, that will thus never be able to grow again to counterbalance the outspiraling effect.In conclusion, when t ev < 4t cc , the system initially outspirals and continues to do so until the PBHs eventually evaporate completely, leading to the disappearance of the binary.In the other case, no definite conclusion can be stated by the sole inspection of Eq. (12).When 4t cc < t ev , we are in a regime where, initially, the radiation of GWs dominates, thus R decreases.But as R and m diminish, there might be a counterbalancing effect due to the term of relative mass loss.Using the explicit expression (10) and solving the Bernoulli differential equation provides When t cc → ∞, the second bracket boils down to unity and we recover the solution where no gravitational waves are emitted with the sole effect of mass loss leading to an outspiralling dynamics.Equation ( 14) is the exact analytical solution to the problem under consideration.
The equation R(t) = 0 has for solution which corresponds to the real time of coalescence.Mathematically, this solution is well-defined only if 1 − 6t cc /t ev > 0 which leads to the condition If Eq. ( 16) is satisfied, the equation R(t) = 0 always admits a solution and thus the two PBHs do coalesce, in a similar fashion to what is observed in the case of constant masses.However, as shown in Fig. 1, because of the mass loss which renders the system less tightly bound, t coal is always larger than t cc .
0  15) divided by the time of evaporation tev (purple curve) in function of the ratio tcc/tev.At tcc/tev = 1/6, this function is finite equal to 1, and for higher values, it ceases to be mathematically well-defined.The grey dashed line is the identity function.
This immediately raises the question of what happens in the small window of parameters for which 4t cc < t ev < 6t cc .Equation ( 13) indicates that we are then in a regime where the emission of gravitational waves initially dominates, thus leading, in a first step, to the decrease of the orbital separation.However, since the condition ( 16) is not fulfilled, the equation R(t) = 0 does not admit any solution, which indicates that the system will never merge.This is due to the evaporation term that gets progressively more important and eventually dominates the GW emission.Once the evaporation process leads, there is no coming back and it ultimately makes the two PBHs outspiral, until they completely disappear.This is an interesting and highly non-trivial case where the orbital separation evolution is non-monotonic.
To sum up, there are three distinct regimes (see Fig.

2) :
• If t ev < 4t cc , the system outspirals due to the domination of the evaporation process and R increases with a final divergence corresponding to the full evaporation of the two PBHs; • If 4t cc < t ev < 6t cc , the two PBHs come closer together in a first step due to the emission of GWs and then, as the evaporation takes the lead back, they ultimately outspiral as in the first case; and • If t ev > 6t cc the emission of GWs dominates entirely, leading to the merger of the system, in a similar fashion to what is observed in the case of constant masses.The only real influence of the mass loss is to increase the time of coalescence.
Evolution of the orbital separation R (normalized to the initial orbital separation R0) between two evaporating black holes with initial mass m0 = 10 12 kg, as a function of time (in seconds).Three distinct regimes arise : for tev < 4tcc (green curve) the binary outspirals whereas for tev > 6tcc (blue curve) it inspirals following a trend similar to the standard case of constant masses.For 4tcc < tev < tcc (orange curve), one has an intermediate regime mixing the two behaviours, with an unusual non-monotonic evolution.
The above inequalities can be conveniently expressed as conditions on the initial orbital separation R 0 for a given initial mass m 0 by using Eq.(11).If one defines and the ouspiralling regime is reached when This might be relevant when evaluating the merging rate of PBHs [31][32][33].

B. Analysis of the frequency
Since the system is assumed to be Keplerian at every step in its evolution, using Kepler's third law ω 2 = 2Gm/R 3 along with Eq. ( 14) provides an analytical expression for the orbital frequency ω: with ω 2 0 ≡ 2Gm 0 /R 3 0 the initial orbital frequency.As for the orbital separation, it admits three distinct regimes (see Fig. 3).The frequency of gravitational waves is simply twice the orbital frequency.
As expected, when R → ∞ in the outspiralling case, ω → 0. This basically means that the two bodies are not sufficiently tightly bound to be considered anymore as a binary system, hence the very notion of orbital frequency becomes irrelevant.When the system actually merges (blue curve in Fig. 3), we recover the familiar chirping trend.To determine the chirping frequency, we place ourselves at the innermost stable circular orbit (ISCO) which corresponds to an orbital separation R ISCO = 12Gm ISCO /c 2 where m ISCO is the time-varying mass evaluated at t = t ISCO .Thus one can solve the equation R(t) = R ISCO using Eq. ( 14), find the corresponding time t ISCO and then evaluate ω(t ISCO ) with Eq. (19).Since R(t) = R ISCO is an implicit equation that doesn't admit any simple analytical solution, this should be done numerically.
Nonetheless, two important features can be noticed.First, as for the standard case, t ISCO is very close to the effective time of coalescence t coal .Solving R(t) = R ISCO for various values of m 0 (typically from 10 9 kg to 10 20 kg) shows that, if γ is the order of magnitude of the initial mass, i.e. if m 0 ≃ 10 γ kg, then t coal − t ISCO ≃ 10 γ−34 s.For PBHS, one can thus legitimately approximate t ISCO ∼ t coal .Conveniently, we readily have an explicit expression for t coal , see Eq. (15).However, by construction, ω(t coal ) = ∞.Although tiny in itself, the difference between t ISCO and t coal is, in principle, relevant to determine the final frequency.
Analytically, the expression ( 19) is not of great use and one should simply go back to Kepler's third law.In fact, whenever the binary system actually merges, the relative variation of mass during the process is small ; more specifically, we shall demonstrate that it is perturbatively small at ISCO.Let us introduce where the last equality comes from Eq. (10).Using that t ISCO ∼ t coal and Eq. ( 15), we are led to Looking at Fig. 1, unless we put ourselves right at the saturated condition t ev = 6t cc -this case is treated at the end of the discussion -we have t ISCO /t ev ≲ 0.4 which leads to ε ≲ 0.2 ≪ 1, i.e. the variation of mass is always small whenever the binary system merges.Consequently, one can rely on the usual expression for the frequency at ISCO, but evaluated for bodies that have lost a fraction ε of their initial masses such that m(t denotes the frequency at ISCO for the evaporating binary system (resp.for a binary system formed of two constant masses m 0 ), then hence for an evaporating binary system, we expect (slightly) higher frequencies.For numerical estimates, plugging Eq. ( 20), supplemented by Eq. ( 21), into Eq.( 22), one is led to: where M ⊙ is the stellar mass.As for the strain amplitude at the ISCO, it is basically given by the usual formula [34]: where M c is the chirp mass -simply given by M c = m/2 1/5 in the case of bodies of identical masses -and D is the distance to the observer.Applying the same reasoning than for the frequency associated with the ISCO, one easily obtains that h H max /h cc max = 1 − ε < 1, i.e. using again Eqs. ( 20) and ( 21): Let us now examine the specific case for which t ev = 6t cc .According to Fig. ( 1) -or equivalently Eq. ( 15) -it means that the system merges concomitantly to the full evaporation of the two BHs, i.e. t coal = t ev .Furthermore, plugging the condition t ev = 6t cc into Eq.( 14) provides which proves that in that case, the system does actually merge.On the other hand, using Eq. ( 19), it is easy to see that the frequency is simply given by We recover that at coalescence, the frequency formally diverges.The analysis of the frequency reached at the ISCO can however be here conducted in a fully analytical way.Solving R(t) = R ISCO indeed leads to One can plug the above expression into Eq.( 27) and use the explicit form of ω 0 , as well as the condition t ev = 6t cc , which enables to express the initial orbital separation R 0 as a function of the initial mass m 0 and some constants (or vice versa), to obtain that the frequency associated at the ISCO is given by: Nicely, this formula does depend neither on m 0 , nor on R 0 .As the ISCO is reached very late in the process (one should keep in mind that, in this case, the mass decreases during the inspiral and vanishes at merging), the frequency is huge, not far from the Planck frequency.This is obviously a purely academic situations but it enlightens the behaviour of the system in the most extreme (merging) case.Interestingly, one should notice that the binary system always reaches its ISCO, which was not a priori obvious, as R ISCO = 12Gm/c 2 decreases as R decreases.The two curves (R(t) and R ISCO (t)) however necessarily intersect each other.

IV. BONDI ACCRETION OF PHANTOM DARK ENERGY
We now turn to the study of binaries formed of two black holes accreting phantom dark energy.As in Ref. [7,8], to describe the accretion of matter from interstellar medium by a compact object, we shall use the standard Bondi approximation corresponding to a rate of mass change ṁ ∝ −m 2 , which integrates into the typical time of evolution τ reading, in the context of dark energy [8,18], where M ⊙ is the Solar mass, w = p/ρ is the equation of state parameter, ρ c ∼ 10 −26 kg/m 3 is the critical density of the universe, and ρ d (∞) is the dark energy density (which is of the same order of magnitude as shown by observations [35]).For numerical estimates, we assume in the following 1 + w ∼ −0.1.Several remarks are in order at this point.First, let us emphasize that we focus here on phantom energy (w < −1), that is on "anti-accretion" (accretion inducing a mass decrease) as the case w > −1 leads to a standard accretion during which the mass increases.Although interesting in itself, this situation does not bring any significant new features to the evolution of the orbital separation: both the effects of gravitational waves and of accretion play in the same direction.The case of a pure cosmological constant, that is w = −1, leads to no mass variation at all.
Second, it should be emphasized that in the case of a standard accretion, as it will be discussed later in this article, Eq. ( 30) would lead to a mass divergence in a finite amount of time.Although not directly related to this work, this opens interesting phenomenological features [36].The situation is, in a sense, formally closealthough reversed -to the one of the Hawking evaporation.The deep reason for this is quite simple.In the case of Hawking evaporation, as in the case of standard accretion, the mass variation gets amplified by the evolution it generates.When it evaporates, a BH becomes smaller and smaller, hence hotter and hotter.The phenomenon diverges in finite time.Exactly as what happens for standard accretion: the more a BH absorbs usual matter (with w > −1), the larger the cross section gets and the faster its mass grows, leading to a divergence in finite time.Both cases are fundamentally unstable and the relevant question is basically to wonder if the coalescence (if any) happens before or after the singularity associated with the mass variation.On the other hand, the case of phantom energy accretion is stable, hence the regular behaviour of Eq. ( 30), whatever the value of t.This is because, when w < −1, as the anti-accretion takes place, the BH gets smaller and smaller and, therefore, sees its cross section and, consequently its mass loss rate, decrease with time.It is a negative feedback whereas the diverging cases correspond to positive feedbacks.

A. Evolution of the orbital separation
Focusing on the case of phantom dark energy, the same method as used in section III can be applied.The differ-ential equation now reads which clearly shows that, depending on an initial setup favouring either inspiral or outspiral, there is always one of the two terms which presents a competitive behaviour between the time evolution of the mass and its orbital separation counterpart (in addition to the obvious competition between the two terms themselves, simply due to the opposite sign).Explicitly: if R decreases, the evolution of the amplitude of the first term is not obvious as m also decreases, making the fate of m 3 /R 3 a priori non-trivial.If, the other way round, R increases, the evolution of the amplitude of the second term is now not obvious as m still decreases, making the evolution of mR possibly intricate.At an initial time t 0 = 0, one has Ṙ and the integration of Eq. ( 32) provides Solving the equation R(t) = 0 gives the following effective time of coalescence which is well-defined mathematically only for τ > 14t cc .As a consequence, the system outspirals for τ < 12t cc , inspirals when τ > 14t cc and for 12t cc < τ < 14t cc we recover this so-called intermediate regime beginning with an initial inspiral but with an ultimate outspiral due to the time-varying mass.Contrary to the case of Hawking evaporation though, the outspiralling dynamics does not lead to any divergence of the orbital separation at finite time, as it can be seen in Fig. 4. As the evolution law for the mass -Eq.( 30) -does not present any pole, the orbital radius gently tends to infinity for t → ∞.

B. Evolution of the frequency
It is straightforward, using Kepler's third law, to obtain an analytical expression for the orbital frequency.Its evolution is shown in Fig. 5.The ISCO analysis performed in the case of the Hawking evaporation is still applicable.One can indeed demonstrate that t ISCO ∼ t coal and plotting the curve t coal /τ using Eq. ( 4) shows that the previously introduced ε parameter is such that ε ≲ 0.23.Then an approximate formula for the frequency associated to the ISCO when the system is in the inspiralling regime is .We recover the outspiralling dynamics for τ < 12tcc (green curve), the inspiralling one for τ > 14tcc (blue curve), and the intermediate non-monotonic regime for 12tcc < τ < 14tcc (orange curve).
with the time of coalescence t coal being given by Eq. ( 35).Again, the correction brought by the time-varying mass is meager.Carrying out the same analysis, but with a non-phantom dark energy, i.e. with w > −1, shows, as expected, that the binary system always merges.In addition, because the mass growth enhances the inspiralling effect, it does so with a shorter time of coalescence.However, as we shall show later on, for all ranges of mass, the difference between these two times of coalescence is too small in practice to constitute a reliable experimental criteria to differentiate the case w > −1 from the case w < −1.

C. Strain -observational distance relation
We now consider the possibility, in the light of our results and in the same vain as in Ref. [8], to use observations of a binary system to measure local dark energy properties.To get an order of magnitude one can safely assume that, for the accretion of dark energy not to be negligible, the two terms entering the differential equation (32) must be of the same order of magnitude.This implies that the two bodies should be then separated by (at least) a distance R ∼ 128 15 Using Kepler's third law and the above relation to express the frequency solely in function of the mass, and injecting the result in the strain formula, Eq. ( 24), shows that, given a minimum detectable strain amplitude h min , and demanding h > h min , constrains the distance to the observer D to be less than where the normalizing value 10 −23 corresponds to the order of magnitude of the minimum measurable strain amplitude currently reached by gravitational wave detectors (see e.g.[37]).This expression depends only on the mass as the orbital separation is chosen so that the anti-accretion term plays a significant role.It should be stressed that it is very different -not only in prefactors but also in the mass functional dependence -from the usual formula encountered in gravitational wave physics, D < 1.6(m/M ⊙ )/(h × 10 23 ) Gpc, as the latter assumes that the system is seen at the merging time.This is far from being the case for the situation we consider here.
Beyond the distance D max , even if the radial change due to the accretion of dark energy were big when compared to the one due to the emission of GWs, the strain produced would be too small to be detected.For solarmass BHs, with h min ∼ 10 −23 , the corresponding order of magnitude is approximately the size of the solar system.However, for supermassive black holes with m = 10 8 M ⊙ , Eq. ( 38) gives D max ∼ 29 Gpc which is more than the cosmological horizon radius corresponding to the limit of the observable universe.It therefore seems that wherever a binary system of 10 8 M ⊙ BHs, seen early enough in the inspiralling process, exists, the accretion of phantom dark energy could be measured on Earth.This is however not that simple.
Whenever phantom dark energy dominates over gravitational radiation, the binary system is in an outspiralling regime, the orbital frequency thus decreases (see the orange and green curves in Fig. 5).This is an appealing feature as it never happens in standard cases corresponding to a dynamics driven only by the emission of GWs.It could be a "smoking gun" for an effect beyond the standard model (BSM).For the strain to be high enough to be detected without requiring the system to violate obvious bounds (as it would, e.g., be the case if one assumes the existence solar mass BHs within the solar system), the best candidates are supermassive BHs.Equation (37) shows that they should be separated by a distance of order R ∼ 10 18 m, corresponding to a time of merging above the Hubble time.However, although the strain would actually be measurable and the sign of the frequency drift would be favourable (that is negative and, therefore, distinct from any usual phenomenon), the amplitude of the variation of frequency -which stands as the relevant parameter -would be ridiculously small.Fixing orders of magnitude by assuming constant masses, one would have | ḟ /f | ∼ 10 7 f 8/3 which leads to | ḟ /f | ∼ 10 −30 for m = 10 8 M ⊙ .By inspection of Fig. 5, it is obvious that taking into account the variation of the mass and the precise form of the frequency would only worsen the effect.In addition, the absolute value of the frequency would also be way below anything measurable in the next decades.
It should be noticed that, in the context of binaries of evaporating BHs, the question of the maximum distance of observation was already considered in Ref. [38]: BHs undergoing the Hawking effect (with masses typically below M * ≡ 10 12 kg) would have fully evaporated in a time smaller than the age of the Universe.There is no chance to observe them in our necessarily confined scope of observation.In this case, the phenomenological relevance of our work is not about the detection of a single event but, rather, about the statistical features of the expected merging rate.

V. GENERIC RESULTS
Beyond the examples of (potential) phenomenological interest given in the previous sections, we now focus on offering an exhaustive classification of the possible behaviours of binaries emitting GWs submitted to mass variation effects described by either a positive or a negative power law.The analytical solutions we give -not previously known in the relevant literature to the best of our knowledge -can be used in many different situations.
We consider mass loss rates of the form: with β > 0 in all cases.For now, we ignore the cases for which the above differential equations do not lead to a mass evolution following a power law -that is corresponding to ṁ ∝ m which implies an exponential trend.Those specific types of evolution are dealt with at the end of this section.
Let us comment about the physical signification of some particular cases.Hawking radiation corresponds to (i) with k = 3.A constant mass loss rate -as it would be the case for a star during most of its life -corresponds to (i) with k = 1.A standard Bondi accretion by a black hole -that is to say a spherically symmetric, adiabatic, steady-flow gas accretion -corresponds to (iii) with k = 1.A Bondi anti-accretion -that is with a phatom fluid -corresponds to (ii) with k = 1.A naive mass loss rate proportional to the area of the compact (non-black hole) object corresponds to (ii) with k = 1/3.A hypothetical situation where the accretion rate of a BH would be solely driven by its surface gravity corresponds to (iv) with k = 2.
Each of these differential equations for β constant are separable and easily integrate into : It should be noticed that (ii) and (iv) exhibit behaviours qualitatively different from (i) and (iii).In the latter two cases, there is a critical point in the evolution of the mass: for (i) it corresponds to t = t ev , where the mass suddenly plunges to zero whereas for (iii) it corresponds to the divergence at finite time t = τ * .As previously explained, those cases are physically related to unstable processes driven by positive feedbacks.Conversely, the cases (ii) and (iv) correspond, respectively, to a smooth decrease and a smooth increase of the mass.They are related to stable processes with negative feedbacks.Then, τ and τ can be understood as typical time scales on which the considered body has lost or gained a significant amount of mass.
The differential equations satisfied by the orbital separation R are respectively given by : (i) Ṙ = − 128 5 which are all Bernoulli equations.The solutions to these four differential equations are given by : The times of coalescence are given by (i) The different situations are deeply nonequivalent.
• In the case (i), the system outspirals whenever t ev < 12 k t cc .It first inspirals and, then, outspirals, for 12  k t cc < t ev < 15+k k t cc .It only inspirals when t ev > 15+k k t cc .• The case (ii) has to be dived into two sub-cases: -If k < 3, the system behaves as previously, but with different bounds.It outspirals for τ < k t cc .This new behaviour was not encountered in the phenomenological cases presented before.It is illustrated in Fig. 6 with k = 4 (which would correspond to a mass loss rate proportional to T 4 , with the temperature T proportional to m 5/4 as expected in some stellar models).
-If k = 3, there is no intermediate regime.The system outspirals for 4t cc < τ and inspirals for 4t cc > τ .
• In the cases (iii) and (iv), there is only one regime corresponding to a pure inspiralling behavior leading to the coalescence of the binary system.This was expected since a gain of mass enhances the effect due to gravitational waves.In the case (iii), the mass diverges in a finite amount of time but the system coalesces anyway, either before (when k < 15) or precisely at this time (when k ≥ 15).This is not trivial and this results from the fact that the very rapidly growing mass also speeds up the inspiral.
The richness of the solutions encountered lies in the asymmetrical role played by the different terms entering the differential equation: the mass variation can either tend to decrease or increase the orbital separation (depending on the sign of ṁ) and can either be smooth of catastrophic (depending on the exponent of m) but the GW term always plays in the same direction, although with a strength coupled to the mass loss.
Let us now turn to the case where ṁ = ±βm, β > 0 which leads to a mass evolution exponentially decreasing or increasing, depending on the sign.The above differential equation integrates into m ± (t) = m 0 e ±t/t * , t * = β −1 . (44) The differential equation for the orbital separation simply reads Ṙ = − 128 5 tcc the binary system outspirals (green curve) whereas for τ > 3tcc it inspirals (blue curve).The regime 11  4 tcc < τ < 3tcc presents a non-monotonic evolution with an initial outspiral followed by an inspiral when the gravitational radiation takes the lead and makes the binary coalesce (orange curve).and integrates into these two solutions (again, depending on the initially chosen sign): In the case R + , The dynamics is always outspiralling.The exponential accretion overwhelms any effect coming from gravitational radiation.It can be clearly seen analytically by taking the limit t → ∞ in Eq. ( 46) for which the term in brackets vanishes and at large t, the overall dynamics is dictated by e 3t/t * .Furthermore, it is straightforward to notice that the equation R(t) = 0 does not admit any mathematically sound solution.Conversely, for R − , the binary system always merges with a time of coalescence given by As it could have been expected, in both cases, the decreasing or increasing exponential evolution of the mass completely overwhelms the effects coming from gravitational radiation.Although mathematically rigorous, this case of exponential evolution of the mass might suffer from the fact that the formal expression of the radiated power by gravitational waves should be modified (beyond the simple accounting for the m(t) term) when one allows a time-dependent mass, notably including terms related to the rate of variation of the mass.This point is discussed in the following section.

VI. RADIATED GW POWER AND CIRCULARITY CONDITION
In section II, we have taken into account the radiated power of emitted gravitational waves by substituting in the standard formula a varying mass.One should however be careful about the fact that the standard formula for P gw is derived with the assumption of a constant mass.For circular orbits and time-varying masses, with orbital frequency ω, the second mass moments M ij = µx i (t)x j (t) read: since for identical masses, the reduced mass is given by µ = m/2.The total radiated power due to the emission of gravitational waves involves the third temporal derivatives of the second mass moments (48) since it is given, in the quadrupole approximation, by [16] where the average ⟨•⟩ is in the time domain over several characteristic periods of the gravitational waves.Using Eq. ( 48), we get If only the last term in the average in Eq. ( 50) is taken into account, one recovers the formula (5).All the other terms involve derivatives (up to order three) of the mass and are thus the direct consequence of allowing timevarying masses.If m (n) ≡ d n m dt n denotes the n-th derivative of the mass, with n a positive integer, then these corrective terms can be neglected if which imposes that the mass of the bodies is slowly varying in comparison to the typical orbital evolution of the binary system.Then, over the temporal window on which the average is performed, that is to say on a few characteristic periods of the gravitational waves, the mass remains nearly constant equal to m(t).
The second assumption used to compute the radiated GW power thanks to Eqs. (48) and (49) relies on the fact that the orbit remains circular.In the standard case of bodies with constant mass, the circularity is ensured by demanding | ω| ≪ ω 2 [16].We investigate if this condition is changed by the variability of the mass.By taking the derivative of Kepler's third law, one has Making use of the triangle inequality Since | Ṙ| ≡ v ∥ is the radial velocity of the system, while ωR ≡ v ⊥ represents its tangential velocity, Eq. (53) shows that the orbits are (quasi)-circular, i.e. v ∥ ≪ v ⊥ , if | ω| ≪ ω 2 and | ṁ| ≪ mω.This last condition exactly corresponds to the case n = 1 of Eq. (51).Under those two conditions, the approximation of a circular orbit with slowly varying radius is still applicable, even in the case of a variable mass system.In conclusion, the slowly varying mass condition, encapsulated by the conditions (51), which enables to simplify the radiated GW power (50) to Eq. ( 5), is linked to the condition of (quasi-)circularity of the orbital trajectories during the inspiral phase.
To conclude this section, one should emphasize that the generic results of section V, although perfectly correct from a mathematical point of view, might suffer some lack of physical relevance in certain cases.Indeed, for mass evolutions presenting a stiff variation, the derivatives ṁ, m, ... m cannot be legitimately neglected in the radiated power P gw and Eq. ( 50) should be used instead, thus providing supplementary terms to the differential equation on R. Investigating the precise effect of these additional terms would constitute a work in itself, as it might as well, as explained just above, have some consequences on the form of the orbital trajectories.Let us however note that, for instance, in the context of BHs submitted to Hawking evaporation, the slowly varying mass assumption is accurately satisfied -except for the higly fine-tuned case for which t coal = t ev which might have an influence on the result of Eq. ( 29).

VII. CASE OF A HIGH HIERARCHY OF MASS
Henceforth, we have examined the case of a binary system constituted by two identical masses m(t), following the same rate of variation.We now turn to the case where one of the bodies forming the binary has constant mass M whereas the other body has a varying mass m(t), assuming that the masses obey the hierarchy ∀t, M ≫ m(t).In this case, the orbital angular momentum reads and its conservation leads to As Kepler's third law involves the sum of the masses of the two bodies, it now simply reads ω 2 ≈ GM/R 3 .Furthermore, the orbital energy is so that the power radiated by the mass loss is One must also be careful about the expression of P gw in this case.We established in the previous analysis that, in full generality, one must take into account the variable mass through the computations of the derivatives of the second mass moments.It is still the case for a high hierarchy because the second mass moments M ij depend on the reduced mass µ ≈ m for M ≫ m.Since Kepler's third law shows that ω 2 ∝ M , the circularity condition of the orbit is satisfied only by demanding | ω| ≪ ω 2 , as in the standard case of constant masses.Consequently, if one still wants to get rid of the higher derivatives terms of m in the expression of P gw (i.e.ṁ, m and ... m), it boils down to an additional assumption put by hand, which has -contrary to the case of identical evaporating masses -no link whatsoever with the circularity of the orbit.Under this hypothesis, one can use The generic differential equation satisfied by the orbital separation (i.e. the counterpart of Eq. ( 7)) is consequently Ṙ = − 64 5 which is again a Bernoulli differential equation which can be analytically integrated given an explicit expression of the mass.If one takes the same mass evolutions as the one given in Eq. ( 40), then Eqs. ( 42) and ( 43) are still valid up to the following systematic substitution (in the exponents as well as in some prefactors, the various signs being unchanged) In addition, the typical time t cc now reads t cc ≡ 5 256 One can also investigate the converse situation for which it is the biggest mass that varies in time i.e. ∀t, M (t) ≫ m and m is constant.Let us note that for Hawking radiation for instance, this situation is physically irrelevant since bigger BHs are less sensitive to the effect.For completeness, however, we shall treat all the cases.The conservation of the orbital angular momentum, using Eq.(54), leads this time to In the expression of the orbital energy, and thus in the power radiated by the variation of mass, m and M play symmetric roles, so Eq. ( 56) is unchanged and we have The conservation of the energy leads to the following (Bernoulli) differential equation for the orbital separation Ṙ = − 64 5 If one takes the same mass evolutions as the ones given by Eq. (40), then Eqs.(42) and (43) are still valid up to the following systematic substitution (in the exponents as well as in some prefactors, the various signs being unchanged) and the typical time t cc is still given by Eq. (61).

VIII. CONCLUSION
Somehow surprisingly, even at lowest order, the dynamics of binary systems with varying mass is highly non-trivial when the emission of gravitational waves is taken into account.
In this work, we have investigated in details the evolution of binaries composed of evaporating primordial black holes with equal masses.We have shown that, depending on the initial conditions, there exist three different regimes.Both the orbital separation and the frequency of emitted gravitational waves are studied under the assumption of circularity.The critical case corresponding to an evaporation time exactly equal to the coalescence time was also considered.
Then, we focused on the case of a Bondi accretion of phantom dark energy.We took this opportunity to correct a mistake made in the literature and to settle a controversy about possible observations of this effect.This required to consider the strain, the frequency and the frequency drift of emitted gravitational waves.
Building on those results, a general study of all possible (power law) cases, including both mass losses and mass gains, was proposed.We gave full analytical solutions and, interestingly, suggested a taxonomy based only on variables easy to determine a priori -that is the typical variation time-scale and the time the coalescence would take with constant masses (not the actual coalescence time).We have shown that quite a few different situations appear, depending on initial conditions, on the sign of the mass variation, and on the regular or "explosive" nature of this variation.The existence of new non-monotonic regimes is underlined.
Finally, we have carefully stated the domain of validity of the different assumptions performed throughout the work and considered also the case of a high mass hierarchy between the inspiralling (or outspiralling) black holes.
In the future, it would be welcome to extend this work to the case of eccentric orbits and to investigate deeper the possible phenomenological consequences.

FIG. 3 .
FIG.3.Evolution of the orbital frequency ω (normalized to the initial orbital frequency ω0) between two evaporating black holes with initial mass m0 = 10 12 kg, as a function of time (in seconds).The color of the three different curves correspond to the ones used in Fig.2.

FIG. 4 .
FIG.4.Evolution of the orbital separation R (normalized to the initial orbital separation R0) between two black holes submitted to a Bondi-type accretion of phantom dark energy with initial mass m0 = 10 12 kg, as a function of time (in seconds).We recover the outspiralling dynamics for τ < 12tcc (green curve), the inspiralling one for τ > 14tcc (blue curve), and the intermediate non-monotonic regime for 12tcc < τ < 14tcc (orange curve).

FIG. 5 .
FIG.5.Evolution of the orbital frequency ω (normalized by the initial orbital frequency ω0) between two black holes submitted to Bondi-type accretion of phantom dark energy with initial mass m0 = 10 12 kg, as a function of time (in seconds).The color of the three different curves correspond to the ones used in Fig.4.

45) 0 1 × 10 1 × 10 FIG. 6 .
FIG.6.Illustration of the behaviour for the case (ii) with k > 3. The plot corresponds to k = 4, with initial mass m0 = 10 12 kg.Panel (a) shows the orbital separation R (normalized to the initial orbital separation R0) while panel (b) displays the corresponding orbital frequency.For τ <11  4 tcc the binary system outspirals (green curve) whereas for τ > 3tcc it inspirals (blue curve).The regime11  4 tcc < τ < 3tcc presents a non-monotonic evolution with an initial outspiral followed by an inspiral when the gravitational radiation takes the lead and makes the binary coalesce (orange curve).