Einstein vs Hawking: black hole binaries and cosmological expansion

This note aims at investigating two different situations where the classical general relativistic dynamics compete with the evolution driven by Hawking evaporation. We focus, in particular, on binary systems of black holes emitting gravitational waves and gravitons, and on the cosmological evolution when black holes are immersed in their own radiation bath. Several non-trivial features are underlined in both cases.


I. INTRODUCTION
An uncharged and non-rotating black hole (BH) in vacuum is fully described by the Schwarzschild solution to the Einstein's field equations.According to Birkhoff's theorem, this solution is actually unique [1].The main expected standard physics addition to this general.relativisticpicture is Hawking evaporation [2], which leads to an explosive mass loss when the BH is light enough.Neglecting higher-order effects, the entire geometric picture is therefore encoded in the treatment provided by general relativity supplemented by the Hawking process.This is obviously well-known [3].However, we shall show in this work that quite unexpected behaviors do emerge in some circumstances, in particular when both effects are simultaneously at play.
In this article, we intend to focus on two specific situations where the evolution driven by classical effects competes with the one induced by the Hawking mechanism.The aim of this study is not to study new physics effects, nor to provide a framework for explicit phenomenological investigations.We simply try to clarify, at the theoretical level, the subtle interplay between classical (i.e.described by Einstein gravity) and semi-classical (i.e.accounted for by the Hawking mechanism) effects.As we will demonstrate in the following, many complicated situations can be encountered.Our objective is to put the emphasis on cases where the dominant contribution could not have been easily guessed and to provide a detailed quantitative analysis of the subtleties entering the game.Although assuming that the considered system is isolated is obviously a strong hypothesis, we insist that within this restriction our analysis is complete.The relative amplitude of the considered effects depends upon the initial values of the masses and on the orbital separation.We provide an exhaustive investigation considering all cases.
The first case we consider consists of a binary system of two BHs.Energy is inevitably radiated away by gravitational waves (although unquestionably correct let us underline that this apparently simple statement is actually not simple [4]).This makes the system inspiral.However, simultaneously, the Hawking radiation induces a mass decrease and this makes the system outspiral.Among other particles, spin-2 bosons -referred to as gravitons in the following -are emitted.We investigate how spacetime excitations produced by gravitational waves compare with those triggered by gravitons.In particular, we show that the total energy radiated in the form of gravitons is always greater than the one carried away by gravitational waves.We also compare the associated frequencies and derive general results.We take this opportunity to provide some clarifications about the highly non-trivial trajectories.
The second situation we consider corresponds to the cosmological expansion if the Universe is filled with a gas of black holes.If space is maximally symmetric, the full general relativistic geometrical picture is encoded in Friedmann's equation for a matter fluid.However, if the BHs evaporate, things get more involved.On the one hand, the cosmic expansion favors matter over radiation (as the dilution is faster for radiation).But, on the other hand, matter gets converted -at an increasing rate -into radiation by the Hawking process.Who wins and why?We study in details the situation and provide new analytical results.Throughout all this work, we put the emphasis on clarity.We also keep the full SI units, making numerical estimates, if required, easier than in Planck units.
In principle, the dynamics of an isolated binary system of Schwarzschild black holes is fully determined by both the emission of gravitational waves and the mass loss induced by the evaporation.Depending on the mass and on the distance between the BHs, one of these effects can, of course, be very strongly dominant over the other.Considering a priori both phenomena is however the rigorous approach to get an exhaustive picture of what happens in vacuum.In addition, we shall precisely show in the following that both phenomena can simultaneously play an important role.The relevant mass range will be precisely calculated.From a phenomenological perspective, only primordial black holes -for which accretion is anyway negligible -are relevant in this case.We, however, keep our study at the general theoretical level.As gravitational waves tend to induce an inspiralling trajectory whereas the mass decrease makes the system outspiral, this obviously triggers a competition between both effects.This has been investigated recently in [5], assuming circular orbits and performing the calculations at the lowest meaningful order.The same hypotheses will be made in this work.
Let us make a few comments about the circularity hypothesis.First, it is well known that back-reaction induces a circularization of the orbit.The eccentricity decreases quite fast and, for all realistic situations, becomes negligible long before the coalescence phase [6].Second, the third time derivatives of the second mass moments M ij ≡ µx i (t)x j (t) can be computed for elliptic orbits [7], including the additional contributions arising from the time-varying masses along with those produced by the orbital motion.These additional terms involve an angular dependency upon the true anomaly ψ and the eccentricity e of the elliptic orbit but these contributions are all bounded to numerically evolve between 0 and 1.What is interesting is that, at the dimensional level, when the two masses of the binary system are equal (let us denote them by m), these terms scale as ṁΩ 2 0 a 2 , mΩ 0 a 2 , and ... ma 2 , with a the semi-major axis of the ellipse and Ω 2 0 ≡ 2Gm/a 3 the fundamental frequency of the emitted gravitational waves for elliptic orbits.For a generic value of the eccentricity e, all harmonics contribute to the frequency spectrum, i.e. there are waves radiated at all frequencies Ω l = lΩ 0 for all integer values of l.
On the other hand, the dimensional prefactor appearing in the classic formula of Peter and Mathews [8] when the masses are taken constant (which is obviously also present in the formulae that generalize this situation to time-varying masses) is mΩ 3 0 a 2 .Consequently, as long as with n an integer, the corrective terms can be neglected (see appendix).In fact, this condition is the exact analog of the one derived in [5], the so-called slowly varying mass condition, which imposes that the mass of the bodies is slowly varying when compared with the typical orbital evolution Ω 0 of the elliptic 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 and equal to m(t).All in all, this means that if this hypothesis holds during the Hawing process (which is the case except during the dramatic end), if one were wishing to take into account the non-circularity at the beginning of the process, the power radiated should simply be corrected by the classic factor [8] f (e) = 1 This, however, does not change the order of magnitude of our conclusions.
The two-body problem with variable mass was initially considered in [9,10].Using the conservation of the angular momentum, it is straightforward to show that the equation of motion reads where R is the orbital separation and m is the (variable) mass of each black hole.We assume that the BHs have the same mass as this is the most formally interesting case and as several models naturally lead to nearly monochromatic mass spectra (see [11] for a review of formation mechanisms of primordial black holes).This work is mostly conceptual but if it were to have any phenomenological relevance, the focus should clearly be on primordial black holes (PBHs) as the mass loss rate would otherwise be too small to change substantially the trajectory.Hence, the reference to a mass spectrum, as PBHs are expected to be formed in the early universe by phenomena possibly involving a continuum of scales.Should the BHs making the binary system present a high mass hierarchy, the factor 3 entering Eq. ( 3) would have to be replaced by a factor 2 or 1 (depending on which BH is actually loosing mass).
On the other hand, the equation of motion for a binary system (with fixed masses) emitting gravitational waves reads [6] ṘR 3 = − 128 5 where c is the speed of light, and G is the gravitational constant.Both equations easily combine [12], leading to: which is the general differential equation relevant for our problem.
In [5], we have shown that this quite simple equation leads to a surprisingly rich landscape of behaviors.
We now focus on a specific question: which process releases the highest amount of power in the form of spacetime vibrations?Otherwise stated: how does the gravitational wave power, emitted by the orbiting black holes, compare with the graviton power, due to the Hawking evaporation?The answer is not a priori obvious and deserves a specific study.The Hawking process converts directly mass into radiation which represents a kind of ultimate efficiency.However, when the system merges the emission of gravitational waves can also extract a huge amount of energy from the orbit.Who wins?Why and when?Incidentally, how does the frequency of gravitons compare with the one of gravitational waves?
A fundamental remark is in order at this stage.Equation ( 5) is deeply asymmetric as the evaporation process is unaffected by the emission of gravitational waves whereas the orbit is drastically influenced by the mass loss induced by the Hawking effect.This makes the situation quite tricky.
Basically, Hawking evaporation is characterized by a (nearly) thermal spectrum with temperature [2] where ℏ is the Planck constant, and k B is the Boltzmann constant.As the BH evaporates it gets lighter and lighter, hence hotter and hotter: the process is explosive and characterized by a negative heat capacity.Importantly, being gravitational, the evaporation is nearly democratic in the sense that the BH couples equally to all the fundamental degrees of freedom.In principle, the spin does play a role both in the accurate expression of the spectrum and in the so-called greybody factors which account for slight modulations with respect to the standard blackbody behaviour [13].Those effects however remain weak and will be neglected in this study.
In addition, we assume that all particles are emitted at the same energy E = ζk B T .This approximates the thermal distribution by its maximum.For integer spins, ζ ≈ 1.6 [14].As gravitons represent 2 degrees of freedom among the 128 of the Standard Model (SM), they will constitute ξ = 1/64 of the total emitted flux.This assumes that the temperature of the black hole is higher than all the SM masses.If the BH is large enough so that this is not true anymore, this assumption will slightly underestimate the relative importance of gravitons.This however does not change any conclusion.
The dynamics of the evaporation can easily be obtained by integrating the instantaneous Hawking spectrum.This leads to: where α H accounts for the number of particles and their spins (see, i.e., [15] for numerical values).Consistently with the previously mentioned hypotheses, we take α H to be constant, which is a good approximation.This integrates into where m 0 is the initial mass and t ev ≡ m 3 0 /(3α H ) corresponds to the time it takes to the BH to fully evaporate.
The dynamics of the orbit is obviously more involved as it depends both on the emission of gravitational waves and on the mass loss.We have shown in [5] that for any initial mass m 0 , it exists three different regimes, depending on the initial orbital separation R 0 .Let us define and When R 0 < R 1 , the system is inspiralling: the emission of gravitational waves dominates the dynamics.When R 0 > R 2 , the other way around, the system is outspiralling as the mass decrease plays the dominant role in the orbital evolution.For R 1 < R 0 < R 2 , the dynamics is non-monotonic: the system first inspirals and, then, outspirals.
With α H = 7.8 × 10 17 kg 3 s −1 , one can numerically express R 1 and R 2 under the form and where M * ≡ 10 12 kg is the typical mass of black holes whose time of evaporation is equal to Hubble time.
In the following, we consider the evolution from the initial time to the evaporation time if the system outspirals.We instead consider the coalescence -or innermost stable circular orbit (ISCO), which is nearly the sametime as the final point if the system inspirals.At the very end of the section, we extend the analysis.

B. Frequencies
Intuition is easier when considering wavelengths.For an evaporating BH, the Schwarzschild radius R S = 2Gm/c 2 is the only length scale of the problem.The wavelength of the emitted gravitons should therefore be of the same order of magnitude: λ gr ∼ R S .
On the other hand, for gravitational waves, another scale enters the game: the orbital separation.Kepler's third law reads, for the orbital frequency, ω 2 = 2Gm/R 3 , which leads to an associated gravitational wavelength λ GW ∼ R 3 /R S at all times.
As the orbital separation is always greater than the Schwarzschild radius (otherwise the system would have merged) it is obvious that λ gr < λ GW .
The order of magnitude of the ratio is λ GW /λ gr ∼ (R/R S ) 3/2 and its evolution is therefore, as expected, only governed by the relative evolution of the two scales of the problem.
Let us now be more quantitative.It is possible to calculate the exact analytical expression for both frequencies as a function of time.Relying on the R(t) expression derived in [5], one can easily obtain for gravitational waves: where ω 0 ≡ 2Gm 0 /R 3 0 is the initial orbital frequency of the system and t cc ≡ 5 512 corresponds to the time of coalescence if the masses were constant (which is therefore always smaller than the actual value).
On the other hand, the frequency for gravitons is simply given by: When the system outspirals (which is always the case when R 0 > R 2 and which happens soon after the initial time for R 1 < R 0 < R 2 ), it is clear that the frequency given by Eq. ( 13) is always smaller than the one of Eq. (15).As the BHs drift apart, the ratio ω gr /ω GW indeed increases (starting with an initial value greater than 1) as the graviton frequency only increases with time whereas the gravitational wave frequency always decreases with time.When the system approaches t ev , that is when the mass of the BHs approaches zero, the graviton frequency tends to infinity -or, at least, to the Planck value -whereas the gravitational wave frequency vanishes.
However, when the system is inspiralling (which is always the case when R 0 < R 1 and which holds briefly after the initial time for R 1 < R 0 < R 2 ), both frequencies increase with time.This is the non-trivial case.This time, the ratio ω gr /ω GW does decrease (obviously still starting from a value greater than one).
The naive dimensional analysis presented at the beginning of this section suggests that both frequencies converge when the system coalesces.As the orbital frequency cannot even be defined beyond the ISCO, we evaluate the final quantities at this point.In practice, going from R = 12Gm/c 2 , that is the ISCO, to R = 4Gm/c 2 does not change the orders of magnitude.At ISCO, we obtain: where m f ≡ m(t ISCO ).The ratio is It is, as expected, of order unity and does not depend on the mass.This means that both frequencies converge at the ISCO (or merging) with a wavelength of the order of the only scale of the problem which is then the Schwarzschild radius.As shown in Fig. 1, the gravitational wave frequency, which starts from a much lower value, generically increases much faster than the graviton frequency.
To summarize, the frequency of gravitons is always higher than the one of gravitational waves.The ratio ω gr /ω GW decreases if the orbital separation decreases and both frequencies become roughly equal at the ISCO.The other way around, the ratio diverges at finite time if the system outspirals.This is not surprising, although the details might not have been guessed so easily.

Gravitational waves
A more subtle question is the one of the comparison between the power radiated by the system as gravitational waves and as gravitons.Let us first focus on the gravitational wave power alone and ask a preliminary question: how does it vary with time?In the approximation scheme previously mentioned, the energy radiated by unit time reads [6] P GW = 64 5 where both R and m are time-dependent (the evolution of the first being, in addition, affected by the second).It can be seen that the gravitational wave frequency ωGW varies much faster than the graviton frequency ωgr.The curve is extended after the red dot corresponding to the intersection but is not physical anymore.
When the system is outspiralling, it is obvious that the gravitational wave power decreases.However, when the system inspirals, the situation gets tricky: the decrease of R clearly tends to increase the power but the decrease of m plays in the opposite direction.Gathering all terms involved in the power, one is led to: with P 0 = 64G 4 m 5 0 /(5c 5 R 5 0 ). Figure 2 shows the typical evolution of P GW and exhibits three main regimes, similarly to the trends that were presented in Ref. [5] for the orbital separation and frequency.The condition Ṗ > 0 is equivalent to This translates into R < R c with This can be numerically expressed as ) Interestingly, this lies between R 1 and R 2 .It means that when R 0 > R 2 the system is only outspiralling and P GW is always monotonically decreasing.When R 0 < R 1 the system inspirals and the power monotonically increases.However when R c < R 0 < R 2 , which is a small but nonvanishing interval, the gravitational wave power does decrease at the beginning of the trajectory while the orbital separation also decreases.There is no one-to-one correspondence between the sign of Ṗ and the sign of Ṙ (see Fig. 3 as an illustration).Basically, this -perhaps surprising -behavior can be understood as follows.For most of the parameter space, the evolution of the gravitational wave power is determined by the orbital evolution.This is simply due to the fact that if the graviton power was much higher than the gravitational wave power, the system would anyway be outspiralling.It is not possible for the binary system to inspiral while having P gr ≫ P GW .However, the power is not what directly dictates the dynamics and this is why, in a narrow window, the power evolution might not follow the evolution of the orbit.
At the ISCO and merging the gravitational wave power reads This raises two important comments.First, those values do not depend on the mass.The power released at the final stage of a merger of black holes is independent of their masses: lighter BHs come closer to one another and this "compensates" for the smaller source term.Of course, the derivative of the power depends on the mass (as 1/m) and so does the full radiated energy.Secondly, the power value at merging is huge, reaching a nearly Planckian value.

Comparison between gravitational waves and gravitons
While the time evolution of the gravitational wave power is given by Eq. ( 19), the graviton power reads: .
In principle, all the needed material is therefore at disposal.The situation is however quite involved and FIG. 3. Normalized derivative of the orbital separation and of the power radiated as gravitational waves.It can be seen that they can be both negative at the same time, which corresponds to a situation where the power released as gravitational wave decreases whereas the system inspirals.
the best way to get a clear idea of the main trends is to focus on the initial and final behaviors of the relevant powers.
Let us first mention that the ultimate fate of the system just depends upon the relative values of R 0 and R 1 .The relative position with respect to R 2 only changes the initial behavior.
If R 0 < R 1 , the final power in gravitational waves is given by Eq. ( 23) while the final power in gravitons is As the gravitational wave power is nearly Planckian whereas the graviton power is "cut" at a mass larger than the Planck mass, the inequality P ISCO GW > P ISCO gr always holds.It approaches equality only if the mass is small at the ISCO.
If R 0 > R 1 , the final power in gravitational waves vanishes while the final power in gravitons diverges (or becomes Planckian if one introduces a regulator in the final stages of the evaporation process).It is then ensured that P gr (t ev ) ≫ P GW (t ev ).
This completely determines the powers at the end of the process.
On the other hand, the initial power radiated as gravitational waves and as gravitons are given by This defines a critical radius such that if R 0 < R G , then P 0 GW > P 0 gr while if R 0 > R G , then P 0 GW < P 0 gr .This can be numerically expressed as (28) At this stage, the initial and final conditions are known.This, however, immediately raises a question: how does R G compare with R 1 ?Interestingly, they do not coincide in general, which means that the initial domination of powers does not determine the final one -although they are obviously not uncorrelated.As expected, the mass dependence in both the formulae of R G and R 1 exhibit very close but not exactly identical exponents: the scaling appears with a power 1.5 for R 1 , and with power 1.4 for R G .This is not so surprising but not a priori obvious.The numerical evaluation shows that both radii coincide for M ∼ 1000 M Pl where M Pl is the Planck mass.The case R G > R 1 is therefore mostly formal and corresponds to a physical situation which, although possible, remains unlikely for phenomenological purposes.
To fully understand the situation, it is interesting to introduce the following mass: This allows a refinement of the previous taxonomy.
• If m 0 > m 1 : if R 0 < R G , gravitational wave power dominates at the initial and final stages; if R G < R 0 < R 1 , graviton power initially dominates but gravitational wave power finally dominates; if R 0 > R 1 graviton power dominates at the initial and final stages.
if R 0 < R 1 , gravitational wave power dominates at the initial and final stages; gravitational wave power initially dominates but graviton power finally dominates; if R 0 > R G , graviton power dominates at the initial and final stages.
Figure 4 illustrates some typical behaviors.Although it is not of great phenomenological relevance, a subtlety should be underlined at this stage: one might assume that if either gravitational wave or graviton power dominates both at the initial and final stages, it should always dominate.This is however not the case.Figure 5 exhibits a double crossing situation where graviton power dominates at the beginning and at the end but not during As a conclusion, except for exceptional cases, the preceding list captures the diversity of behaviors for both powers as a function of the initial conditions.Obviously, the answer to the competition between gravitons and gravitational waves as the main source of space-time vibration did not have a simple intuitive answer and required a somehow extensive investigation.

Comparison between gravitational waves and gravitons
Beyond the comparison of instantaneous powers, it makes sense to evaluate the full amount of energy radiated as gravitons and as gravitational waves during the entire process from the remote past to the ISCO.
The energy emitted as gravitons is simply given by ∆E gr = ξ∆mc 2 which leads to . ( The typical behavior is shown in the bottom right plot of Fig. 6.The energy emitted as gravitational waves is obviously less simple to evaluate as soon as the mass varies.Once again, for t in the interval [t 0 , t ISCO ] we define with P GW given by Eq. ( 19).The three regimes associated with the behavior of P GW are present for ∆E GW in Fig. 6.When the system outspirals, the total energy rapidly reaches a nearly constant value as the radiated power become negligible (both the increase of the orbital separation and the mass loss tend to decrease P GW ).
When the system inspirals, one recovers the usual trend but with slight differences.Figure 7 shows the evolution of the total energy radiated as gravitational waves for the actual system with a varying mass together with its approximation by either a constant mass equal to the initial mass or by a constant mass equal to the final one.As expected, the actual behavior lies in between.
This raises an interesting comment.We have seen previously that the gravitational wave power necessarily increases when the system inspirals (except, in some very specific situations, for a brief amount of time).This was not obvious as the mass loss competes with the decrease of the orbital separation.We now ask a different question: what about the total energy emitted as gravitational waves?Let us assume that the mass of the BHs varies from m 0 initially to m f at the ISCO.Is ∆E GW better approximated by a constant mass equal to m 0 of to m f ?Otherwise stated: which part of the dynamics does release more energy in the surrounding space?Even if the final power dominates it could be that the time integration favors the initial stages where the system spends more time (in addition of being more massive).We do not address this question frontally and the main conclusion will be given at the end of the point 3 of this subsection.t (s) FIG. 6.On the top and bottom left plots, the temporal evolution of the total energy radiated as gravitational waves (divided by the initial value of the power P0) is depicted in each of the three regimes: blue for outspiralling, green for inspiralling, and orange for the intermediate case.These curves correspond to the results from the numerical integration of Eq. (31) while the black dashed ones correspond to analytical formulas given by Eq. (34) and Eq. ( 35).The bottom right plot displays the radiated energy as gravitons given in Eq. (30).
To conclude the present subsection, let us make a brief digression that might help the reader wishing to use explicit formulas.It is indeed possible to compute analytically the integral given by Eq. (31) once the expression of the radiated power P GW given by Eq. ( 19) is plugged into it.To this purpose, we recognize incomplete Beta functions by performing simple changes of variables.For this to work, one quickly finds out that certain mathematical manipulations are allowed if and only if one specifies in which regime -whether outspiralling or inspiralling -we are, which provides a constraint on the sign of the term 1 − tev 6tcc .Consequently, we expect two different expressions for the integrated power depending on the regime under consideration, which reflects the fact that ∆E GW inherits from the behaviour of P GW .Furthermore, for most purposes, it is generally easier to manipulate hypergeometric functions rather than incomplete Beta functions.The link between them, more specifically between B(z; a, b) and the Gaussian (or ordinary) hypergeometric function 2 F 1 (a, b, c; z) is ensured by the following relation [16] Therefore, in the outspiralling case, the change of variable u = (1 − t ′ /t ev ) 6 leads to ∆E GW = 18 123 ; On the top plots and the bottom left plot of Fig. 6, it is shown that these analytical formulas fit perfectly the curves resulting from the numerical integration of Eq. (31).Using the main properties of hypergeometric functions, one can recover the main characteristic behaviors exhibited in Fig. 6, for instance the convergence of ∆E GW to a finite value for t → t ev in the outspiralling case using that lim z→0 2 F 1 (a, b, c; z) = 1.

Fully time-integrated process
In the previous analysis, we focused on a process assumed to take place between an arbitrary initial time and a final moment defined as the merging (or ISCO) time.This led to an impressive diversity of conclusions.It is however also possible to consider the fully time-integrated process without interrupting the reasoning at merging.This leads to a (maybe surprising) generic result: the total energy released as spacetime vibrations by gravitational ways is always of the same order (or smaller) than the energy released as gravitons.
The total energy radiated as gravitational waves between a remote initial time and the ISCO is, if the mass does not vary substantially, On the other hand, the total energy emitted as gravitons over an unlimited amount of time is trivially where the superscript "full" means here that we also consider the evaporation even after merging.If the black hole is hot enough to emit all standard model degrees of freedom, ξ = 1/64.This elegantly leads to meaning that as soon as the dynamics of the system is such that it inspirals -that is as soon as the evaporation is not too fast at the initial stage -the gravitational wave energy and the graviton energy are of the same order.If the system outspirals, gravitons obviously dominate strongly.
Let us now question the validity of this result by considering several subtleties.
First, the energy emitted during the merging and ringdown phase can, obviously, not be captured by Eq. (36) which is simply obtained by evaluating the Newtonian orbital energy at the ISCO.The measurements [17] however suggest that the "missing mass" between the initial black holes and the final one is of the same order as ∆E GW .As a crude approximation, one therefore gets where the superscript "full" means that the gravitational waves emitted during the non-perturbative phase are also accounted for.This clearly increases the amount of GW energy but does not change the orders of magnitude.Second, one could instead consider very large black holes with a temperature (during most of the process) smaller than the mass of all massive particles.In that case, ξ = 1/2.This increases the number of gravitons but does not radically alter the conclusion.Anyway, the energy released as gravitons is of the same order as the energy released as gravitational waves and one has: Third, one can question the fact that the gravitational wave energy is here approximated assuming a constant mass.As we shall see later, there indeed exist very special cases where the system inspirals with a strongly varying mass.Then, gravitational waves are suppressed and ∆E GW < ∆E gr .This however represents a tiny part of the parameter space requiring strong fine-tuning.For the very vast majority of initial conditions, ∆E GW ∼ ∆E gr .
Although mathematically very simple, the fact that both energies are comparable was not a priori so obvious.A binary system of black holes basically "excites" spacetime as much through gravitational waves than through gravitons.
The only important case where a strong hierarchy appears is when the mass is small enough (or the orbital separation is large enough).If it is such that the system outspirals, the amount of gravitational waves becomes obviously negligible.The very specific situation where the system still merges but with a mass at the ISCO which is much smaller than the initial mass will be analytically studied in the next section.
In conclusion, the energy released by gravitational waves is always smaller than the one emitted as gravitons.The bound is saturated for an inspiralling system with slowly varying mass and the hierarchy becomes large if the mass loss becomes sizable during the process.It should be noticed that the amount of emitted gravitons is always (nearly) the same, it is only the power radiated as gravitational waves which highly depends upon the orbit.

Extreme case
The analytical solutions are textbook when the mass is nearly constant.There is however another case, on the extreme other side, for which a simple analytical solution can be obtained.It can indeed be shown [5] that a system whose full evaporation time is exactly equal to the merging time does actually exist (which is not obvious as it might have been guessed to outspiral) and is easy to solve.In this situation, the power radiated by the emission of gravitational waves is given by (41) Using Eq. ( 31), a simple analytic expression for the energy dissipated can be found: 45 256 . (42) This allows to evaluate the ratio of emitted energies at the ISCO: Interestingly, this shows that as soon as the initial mass is not negligible -if it was the case, it would anyway mean that the process would be nearly instantaneousthe gravitational wave energy is highly suppressed when compared to the one associated with gravitons.Although quite expected on the one hand, this is, on the other hand, a non-trivial result as the gravitational wave power emitted at the ISCO is still Planckian.
In the extreme case considered here, the BHs have nearly fully evaporated when they reach the ISCO.This implies that ∆E gr ≈ ξm 0 c 2 .In more details, Eq. ( 43) therefore shows that This answers the question that we previously asked.When the mass loss is substantial during the merging process, the total amount of energy emitted as gravitational waves is neither approximated by assuming a constant mass equal to the initial one nor a constant mass equal to the final one.

III. SPACE EXPANSION
We now focus on a second situation where a competition takes place between the dynamics driven by the classical BH behavior and the evolution due to Hawking evaporation.Let us consider a gas of BHs in the early universe.Those BHs evaporate and produce radiation.On the one hand, the expansion of space tends to favor matter -that is BHs -over radiation as the energy density scales as a −3 for the former and as a −4 for the latter, a being the scale factor.On the other hand, the evaporation process tends to favor radiation over matter as it converts the BH mass into relativistic quanta.It is not a priori obvious to determine the evolution of the system.Is the energy density of the universe dominated by matter or radiation?If a transition takes place, when does this happen and why?
In most previous studies, authors focused mostly on the ratio between the BH energy density and the total energy density (thus comprising both matter and radiation) [18,19], which somehow blurs the competitive effects between radiation due to the evaporation and the contribution of BHs themselves.The first articles to address the question analytically were [20,21].In [20], the focus was on a power law mass spectrum for BHs while we will consider here the monochromatic case, as it makes the situation much clearer.We provide analytical solutions, completing and refining the works of [21,22] as well as correcting certain mistakes in [21].

A. Numerical results
The gas of BHs is modeled as a dust-like perfect fluid (i.e.p m = 0) whose constituents are Schwarzschild black holes, with an energy density ρ m .If we denote by f PBH (m, t) the distribution of BHs of mass m at time t -the so-called mass spectrum -the BH energy density is given by The energy density for radiation is denoted by ρ r .The initial total density is ρ 0 and the parameter β fixes the initial ratio between matter and radiation, namely ρ m (t 0 ) = βρ 0 and ρ r (t 0 ) = (1 − β)ρ 0 .For an evaporating distribution of BHs, the continuity equations are modified by an interaction term (see for instance [22]), leading to the so-called Friedmann-Boltzmann equations.In the case of a purely monochromatic spectrum, we thus have the following set of coupled ordinary differential equations for the energy density of BHs (thus matter) and radiation ρ r : supplemented by the Friedmann equation for two components universe: This set, knowing the evolution law for m, as given by Eq. ( 8), can easily be numerically integrated.The result is displayed by the blue curve in Fig. 8 for β = 0.9.The analytical solution discussed in the following subsections is represented in dashed orange.Several interesting features should be noticed.
First, the ratio ρ m /ρ r is initially increasing.This means that, at the beginning of the process, the differential dilution (favouring matter over radiation) is more efficient than Hawking evaporation.Although BHs lose mass, their contribution to the energy budget of the Universe increases with time.
Second, there exists an inflection point -corresponding to the first vertical line and defining the time t invwhere the derivative of ρ m /ρ r vanishes.It means that at some point, the evaporation becomes efficient enough to overcome the dilution.Very interestingly, this does not happen when the Hawking evaporation becomes deeply explosive, that is when the mass plunges to zero.This can be understood by looking carefully at the mass evolution given by Eq. ( 8).It is correct to think that this behavior, driven by a positive feedback and a negative heat capacity, is somehow catastrophic and leads to a singularity (in the sense that dm / dt diverges) at t = t ev .However, the actual evaporation time differs only by a factor of 1/3 from the one naively extrapolated from the initial evaporation.This basically means that contrary to what is usually assumed, the evaporation is not that slow at the initial time (when compared to the actual evaporation time).
Third, although both the previously mentioned regimes are obviously described by power laws, the absolute value of the slope of the second one is larger.When the evaporation process begins to dominate, it does so more efficiently.Let us underline that by evaporation "domination", we mean here that the derivative of ρ m /ρ r is negative.However, during nearly all the evolution (except, of course, if β is chosen tiny), matter dominates and the scale factor behaves as a ∝ t 2/3 .
Fourth, there exists a final regime -corresponding to the second vertical line -where the considered ratio decreases very fast before vanishing.This corresponds to the "explosion" of the BHs.For most -but not all -of the parameter space, this happens when ρ m /ρ r becomes close to its initial value.
Fifth, one should also be careful not to be fooled by the logarithmic scale used in Fig. 8: in cosmic time, the inflection point happens very soon after the beginning of the process.However, ρ m /ρ r has already considerably increased in this small amount of time.to the time of evaporation tev.The blue curve represents the results from numerical integration while the orange dotted one is our analytical solution.The dashed horizontal line corresponds to y = 1 while the two vertical lines respectively correspond to t = tinv and t = tev.

B. Analytical solutions
In the case of a purely monochromatic spectrum, f PBH remains monochromatic at all times so that Eq. ( 45) straightforwardly becomes with . As displayed in Fig. 9 and as expected, the analytical solution for ρ m matches perfectly the one resulting from the numerical integration.More interestingly, an approaching analytical formula can also be derived for ρ r .As shown in Eq. ( 46), it obeys a differential equation of the form y ′ +p(x)y = b(x) whose solution is given by y(x) = e −P (x) e P (s) b(s) ds + Ce −P (x) where P is the primitive of the function p and where C is a constant fixed by specifying initial conditions.In our case, p(t) = 4H which is straightforwardly integrated into P (t) = a(t 0 ) 4 /a(t) 4 .Furthermore, the right-hand side of the ODE is now completely determined using Eq. ( 8) and Eq.(48).Enforcing that at the initial time t 0 , ρ r takes the value (1 − β)ρ 0 , we obtain: The only difficulty to overcome is to compute the integral involved in the above equation, which necessarily means specifying an expression for the scale factor a(t).At this step, we are coerced into resorting to numerical resolutions of the set of coupled differential equations which show that, as could be expected, a(t) nearly behaves as a power law, but with an exponent that can, in principle, shift between 1/2 (radiation domination) and 2/3 (matter domination).Actually, even when radiation dominates initially -that is to say β is close to 0 -a(t) quickly recovers the matter-dominated trend as it can be seen in Fig. 10.As a first approximation, we shall thus consider that a(t) ≈ At 2/3 with A a constant.This is natural as we have previously shown that ρ m /ρ r strongly increases at the beginning of the process and still dominates even when the evaporation begins to become significant.
FIG. 10.Scale factor a(t) as a function of time, resulting from numerical integration and represented in blue.For comparison, in dashed orange (resp.green), is represented a(t) ∝ t 2/3 (resp.a(t) ∝ t 1/2 ).Even for β = 0.1, which is the value chosen for this plot, it can be seen that the slope 2/3 is quickly recovered.As expected, the fit a(t) ∝ t 2/3 is increasingly better as β approaches unity.
In this case, the integral is analytical and can be split into two parts in order to recognize incomplete Beta functions (see Eq. ( 32)).Using the fact that they can be expressed in terms of hypergeometric functions (see Eq.  (50) As it can be seen in Fig. 11, the above analytical form matches very well the numerical results.However, one should notice that as β gets closer to 0, that is to say as the initial amount of matter becomes negligible, the fit of the first branch of the numerical integration becomes less accurate.This is the direct consequence of having assumed a(t) ∝ t 2/3 during the whole evolution of the system (obviously, when β ∼ 0, the universe is initially radiation-dominated and, at the very beginning of the evolution, a(t) ∝ t 1/2 ).A possible refinement would consist in finding the value t * for which the power-law exponent for the scale factor shifts from 1/2 to 2/3, and splitting the integral in Eq. ( 49) in order to plug the correct expressions of a(t) for t < t * .However, such manipulations only make the computations heavier, bringing neither new physics nor substantial improvements in accuracy.is always negligible when compared to its time-dependent counterpart.This, of course, is not fully correct as the first term in Eq. (50) does not vanish initially at t = t 0 .However, when focusing on most of the temporal evolution of ρ r , this approximation is completely legitimate.Furthermore, as the behaviour of the hypergeometric function is such that [23] lim z→0 2 F 1 (a, b, c; z) = 1, if the system is not close to the typical time of evaporation t ev , it does not contribute.Only for t ∼ t ev , does the hypergeometric function play a role, essentially modeling the fact that the BHs have entered the stiff part of the evaporation process, thus leading to a kink in the radiation energy density contribution.
Consequently, for t < t ev , Eq. ( 50) can be simplified to ρ r (t) ≈ Aρ 0 PBH 5a(t 0 )t ev a(t 0 ) a(t) (51) As highlighted by Fig. 12, the above formula is sufficient to fit the trend of a 3 ρ r except for t very close to t ev .It also gives a mathematical explanation of the reason why the curve reaches a minimum before reaching the second cutoff associated to the full evaporation as a consequence of the competition between the functions Aρ 0 PBH 5a(t0)tev t 5/3 and (1 − β)ρ 0 t −2/3 , the second being pure radiation while the first is modeling the smooth evaporation.Both these quantities present monotonic behaviours and the minimum reached by a 3 (t)ρ r (t) corresponds to their equality.Equating them straightforwardly gives the so-called time of inversion .
Since t ev ∝ m 3 0 , studying t inv as a function of t ev shows that, as expected, the smaller the initial mass of black holes, the closer to the initial time of formation t 0 will t inv be.In the extreme limit where β = 1, this competitive effect is not present anymore and the radiation energy density starts by increasing from 0 to a non-vanishing value thanks to the radiation created by the evaporation process.
It can also be noticed that for any generic initial condition different from the very specific β = 1 case, this inversion always occurs.This detail was missed in the analysis presented in [21] and could be guessed from the plots of [22] but was neither explained nor proved.Although it does not affect the big picture of a matter dominated regime followed by a radiation dominated one at the end of the evaporation process, it highlights that even in its early stage, the evaporation process plays a nonnegligible role that is not restrained to what occurs near t ev .Furthermore, one could have assumed that for light enough black holes -that are evaporating fast enough, as dm / dt ∝ 1/m 2 -radiation immediately dominates.This is not the case.There is no critical initial mass below which this happens.
Let us finally analyze the end behaviour of ρ r .As expected, it can be seen that it converges to a finite value at t = t ev .Our analytical solutions also exhibit these behaviours, through the fact that 2 F 1 ( 5 3 , 2 3 , 8 3 , 1) = Γ( 13 )Γ( 83 ).From here, one can plot ρ r for a range of initial conditions and see that there are some instances, depending upon the initial value of β, for which the final value of ρ r at t ev will exceed the initial value of ρ r , and others such that ρ r converges to a value that is less than that of the initial conditions.The critical value β c of β for which it occurs, is obtained by solving the equation a 3 (t 0 )ρ r (t 0 ) = a 3 (t ev )ρ r (t ev ), using Eq. ( 11) and the convergence of the hypergeometric function.It can be well approximated by C. Relics In this work, we have resisted the temptation to bring any exotic physics into the game and we stood with purely standard processes.Let us make here an exception.Based on several distinct arguments (see, e.g., references in [24]), it is possible -if not probable -that stable relics are formed at the end of the evaporation process, with masses close to the Planck mass m rel ∼ m Pl .This makes the evolution even more interesting.
Starting with BHs with masses larger than the Planck mass, the dynamics begin, as usual, by an increase of the matter contribution over the radiation one (in a sense, space expansion overcomes the evaporation).At some point (t = t inv ), the trend reverses and the relative contribution of matter begins to decrease (evaporation overcomes space expansion) -although it remains dominant for a while.At t ≈ t ev , the evolution becomes dramatic and ρ m /ρ r plunges.If it was not already the case, radiation soon dominates.However, in this case, the matter density does not reach exactly zero as the evaporation stops at m = m rel .Due to the differential dilution, mat- ter (that is relics) contribution slowly begins to increase again.At some point, it eventually drives the cosmological dynamics.Once relics dominate, they obviously keep dominating.
A simple gas of BHs can generate a quite involved cosmological expansion history, as shown in Fig. 13.

D. Cosmological constant
Although the previous discussion would mostly apply to the primordial Universe, it makes sense to add a cosmological constant (assumed to be positive, as measured).This, indeed, remains compatible with the core hypothesis of this work, which consists in studying theoretically the full dynamics in vacuum.The Friedmann equation now reads where ρ Λ ≡ Λc 2 8πG .
Should the ratios of the BH and radiation densities be expressed as a function of the scale factor, the addition of a cosmological constant would change nothing at all, whatever its value.As ρ m and ρ r are decreasing functions of the scale factor whereas ρ Λ is constant, the latter will inevitably dominate at some point.When using a as the evolution variable, this anyway makes strictly no difference.
The situation is obviously different when using the usual cosmic time.At some stage, the cosmological constant dominates.When this happens, the scale factor begins to grow exponentially.As matter is less diluted than radiation, this favors black holes over photons.This is precisely what can be seen in Fig. 14 where a vacuum energy density, comparable to the total matter and radiation density at t = 10 −5 s, was added.As expected, when the cosmological constant dominates -that is for t ≳ 10 −5 s -the ratio ρ m /ρ r decreases more slowly, and eventually increases, until the black holes "explode".In this final stage, taking place around t = 3 × 10 −4 s for the initial conditions chosen in Fig. 14, the ratio necessarily plunges to zero, whatever the background dynamics.
Interestingly, the addition of a cosmological constant can therefore enrich the previously studied evolution.If its value is such that is begins to dominate after the (first) maximum of ρ m /ρ r but before the full disappearance of the black holes, it will lead to the second local maximum before the the unavoidable final decrease.
The influence of the Λ-term on the individual densities of BHs and radiation, ρ m and ρ r , is much less interesting and simply leads to an exponential dilution at the time of dominance of the vacuum energy.

IV. CONCLUSION
In this work, we have considered several situations where Hawking evaporation induces a dynamic competition with classical evolution.Our aim was not to provide "new physics" results but to highlight some non-trivial conclusions that are not intuitively straightforward.
When considering binary systems of Schwarzschild black holes, we have, in particular, studied exhaustively the evolution of the gravitational wave and graviton frequencies, powers, and cumulative energies for all possible scenarios (inspiraling, outspiraling, and non-monotonic).Most conclusions were not a priori obvious.We have also established that the full energy transferred to space-time as gravitons is often of the same order and in any case higher than the one transferred as gravitational waves.
When dealing with the expansion of a space filled with a gas of Schwarzschild black holes, we have proven that the matter and radiation energy densities evolve in an intricate way.The main features are universal and do not depend on the initial conditions.Some useful analytical approximations were also provided.
Several developments could be considered in future works.First, it would be interesting to also consider hyperbolic trajectories [25].The emission of gravitational waves is very different in this case [26].Second, it could make sense to also focus on the graviton emission enhancement taking place if black holes are spinning [27].Finally, possible modifications of the Hawking evaporation happening in the quantum gravity regime (see references in [24]) should be taken into account.

V. ACKNOWLEDGEMENTS
The authors would like to thank Cyril Renevey who provided us with the first version of the code initially used for numerical simulations.

VI. APPENDIX
Let us review here in details the influence of the noncircularity of the orbits.To this aim, we parameterize the elliptic trajectory by its eccentricity e, such that 0 ≤ e < 1, its semi-major axis a, and by the polar coordinates (r, ψ) in the plane of the orbit, with the origin located at the position of the center-of-mass.where the functions g 1 , g 2 , etc. are only depending on the true anomaly ψ and on the eccentricity e.They are bounded to be smaller that 3. Explicitly, they read, in the approximation that that the precession of the argument of periapsis is negligible: (61) For isotropic 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/| ṁ| that is when the gravitational waves coalescence timescale is much longer than the mass-loss timescale.This is the only hypothesis that was made in order to derive these results.It should also be noticed that when enforcing constant masses in to the above formulae, one recovers the canonical expressions of [8].If one doesn't make the assumption of identical time-varying masses, the corrective terms involve derivatives of the reduced mass µ and of the total mass m tot , 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 tot .
From the inspection of the above formulae, it is clear that the corrective terms are negligible in comparison to the PM canonical expression when with n an integer.This constraint is simply the slowlyvarying mass condition, which is satisfied during the whole evaporation process except very close to the time of evaporation for which the mass quickly decreases.

FIG. 1 .
FIG.1.Typical evolution of the frequencies close to the merging.It can be seen that the gravitational wave frequency ωGW varies much faster than the graviton frequency ωgr.The curve is extended after the red dot corresponding to the intersection but is not physical anymore.

FIG. 2 .
FIG.2.Typical time evolution of the power radiated as gravitational waves for three different regimes (in green, the purely inspiralling case, in blue the purely outspiralling case, and in brown an intermediate case).

FIG. 4 .FIG. 5 .
FIG.4.Some sketches of the typical time evolution (from 0 to tev) of the gravitational wave radiated power PGW and its graviton counterpart Pgr, showing the diversity of situations and illustrating the taxonomy presented in the text.The top left plot depicts an outspiralling system, with PGW dominating at the initial and final stages.The top right and bottom left plots depict situations where PGW and Pgr intersect, respectively in the outspiralling regime and in the non-monotonic regime.Eventually, the bottom right plot represents an outspiralling system for which Pgr dominates at all times.

FIG. 8 .
FIG.8.Ratio of matter over radiation energy densities y = ρm/ρr, as a function of time, from an initial formation time t0 to the time of evaporation tev.The blue curve represents the results from numerical integration while the orange dotted one is our analytical solution.The dashed horizontal line corresponds to y = 1 while the two vertical lines respectively correspond to t = tinv and t = tev.

3 ρ m /ρ 0 FIG. 9 .
FIG. 9. Evolution of the comoving matter energy density (normalized to the reference energy scale ρ0) as a function of time from t0 to tev, with the same color codes as in Fig.8.

3 ρ r /ρ 0 FIG. 11 .
FIG. 11.Evolution of the comoving radiation energy density (normalized to the reference energy scale ρ0) as a function of time from t0 to tev, with the same color codes as in Fig.8.The horizontal dashed line corresponds to y = a 3 (t0)ρr(t0) = (1 − β)ρ0 while the two vertical lines respectively correspond to t = tinv and t = tev.

FIG. 13 .
FIG.13.Evolution of ρm/ρr as a function of time when BHs do not fully evaporate but, instead, form stable relics when approaching the Planck mass.The new regime starts here for t ≈ 2 × 10 −4 s.
These two curves exactly intersect at the so-called time of inversion tinv given by Eq. (52).