Cosmological evolution of primordial black holes
Jared R. Rice Bing Zhang






jrice@physics.unlv.edu


zhang@physics.unlv.edu

Department of Physics and Astronomy, University of Nevada, Las Vegas, NV 89154, USA


The cosmological evolution of primordial black holes (PBHs) is considered. A comprehensive view of the accretion and evaporation histories of PBHs across the entire cosmic history is presented, with focus on the critical mass holes. The critical mass of a PBH for current era evaporation is M_cr∼ 5.1×10^14 g. Across cosmic time such a black hole will not accrete radiation or matter in sufficient quantity to hasten the inevitable evaporation, if the black hole remains within an average volume of the universe. The accretion rate onto PBHs is most sensitive to the mass of the hole, the sound speed in the cosmological fluid, and the energy density of the accreted components. It is not easy for a PBH to accrete the average cosmological fluid to reach 30M_⊙ by z∼0.1, the approximate mass and redshift of the merging BHs that were the sources of the gravitational wave events GW150914 and GW151226. A PBH located in an overdense region can undergo enhanced accretion leading to the possibility of growing by many orders of magnitude across cosmic history. Thus, two merging PBHs are a plausible source for the observed gravitational wave events. However, it is difficult for isolated PBHs to grow to supermassive black holes (SMBHs) at high redshift with masses large enough to fit observational constraints.









§  INTRODUCTION

Primordial black holes (PBHs) are among the most intriguing ghosts in the universe. A singular PBH of sufficient mass can navigate the history of the universe without detectable clues to its existence; a true cosmic ghost. Low mass PBHs evaporate before the current epoch and the radiation signature of an isolated high mass PBH is too weak to detect. The last moments of a PBH evaporation reveal the hole through a burst of high-energy radiation that is distinguishable from that of short gamma-ray bursts (GRBs) <cit.>.

The upper limits on the number density of PBHs across a wide range of masses is discussed extensively in <cit.>. To date there are no confirmed PBH burst signals, but these compelling ghosts are ripe cosmological messengers that will enhance our understanding of the universe if observed. The PBHs evolving through cosmic history could be used as a proxy for understanding the conditions in the early universe. PBHs of significant mass may gain a dark matter (DM) halo, e.g. <cit.>. Since the PBH evaporation rate depends only on the mass of the hole and the assumed particle physics model <cit.>, PBHs in similar astrophysical environments should produce similar radiation signatures; the ultimate “standard candles.”

This study explores the evolution of PBHs through accretion and evaporation across the entire cosmic history. Special attention is paid to the changes in the density, temperature, and sound speed in the cosmological fluid because of their influence on the accretion rate of that fluid onto the PBHs. In <ref> the concordance cosmological model of ΛCDM is discussed. In <ref> the PBH accretion and evaporation models are discussed and formulae are given for the accretion rates in the various cosmological eras. In <ref> the results of the study are discussed. Finally in <ref> the conclusions are presented and a discussion of astrophysical implications is made.






§  COSMOLOGICAL MODEL

The concordance cosmology assumed throughout this study is the six parameter ΛCDM model, implementing the most recent Planck Collaboration results <cit.>. The model consists of the homogeneous and isotropic Friedmann-Robertson-Walker (FRW) geometry dynamically evolving according to the Einstein field equations. The Einstein equations, also called the Friedmann equations <cit.> in this case, describe the evolution of the curvature and energy content of the universe as

    (   ȧ/a)^2 +kc^2/a^2 = 8π G/3ρ,
       ä/a = -4π G/3(ρ+3P/c^2),

where the scale factor is a≡ a_0(1+z)^-1, with the scale factor today a_0≡1 and z the cosmological redshift, k=0,±1 indicating zero, positive, or negative spatial curvature respectively, G is the universal gravitation constant, and c is the speed of light. The term ρ is the sum of the proper inertial mass densities of the cosmological fluid and the contribution from spatial curvature, and P is the pressure contribution from matter, radiation, and vacuum energy (or cosmological constant Λ).

The equation of state of each cosmological fluid can be expressed P_i = w_iρ_i c^2 (no sum over i) with equation of state parameter w_i. The equation of state parameters for matter, radiation, cosmological constant, and spatial curvature are 0, 1/3, -1, and -1/3 respectively. Note that the equation of state of a baryonic gas P∝ρ^γ, where γ is the adiabatic index, is relevant when calculating accretion rates onto a compact object. The approximation w_i=b≈0 for a baryonic gas holds on cosmological scales.

The Hubble parameter H is a measure of the temporal (extrinsic) curvature of the FRW geometry and is defined

    (ȧ/a)^2    ≡ H^2 = H_0^2ℰ^2,

where subscript-0 implies evaluation of a quantity today. The Hubble constant is H_0 = 100h km s^-1Mpc^-1 where the dimensionless Hubble parameter is h=0.6774 from Planck <cit.>. Before defining ℰ^2, it is convenient to introduce the dimensionless density parameters today

    Ω_i,0   ≡ρ_i,0/ρ_cr,0 = 8π G/3H_0^2ρ_i,0,

where i indicates baryonic matter, dark matter, radiation, and Λ. Dividing Eq. (<ref>) by H^2 and evaluating the quantities today gives an expression for the `effective' dimensionless density parameter for spatial curvature

    Ω_k,0   = 1 - Ω_0,

where Ω_0 ≡Ω_r,0+Ω_m,0+Ω_Λ,0. The combined Planck and baryon acoustic oscillation data <cit.> are consistent with Ω_k,0 = 0.000±0.005, i.e. the universe has zero spatial curvature to within 0.5% accuracy. The term ℰ^2 in Eq. (<ref>) is a function of the dimensionless density parameters with their redshift dependencies

    ℰ^2    ≡∑_iΩ_i,0(1+z)^3(1+w_i)
       =Ω_r,0(1+z)^4 + Ω_m,0(1+z)^3
       ==+(1-Ω_0)(1+z)^2 + Ω_Λ,0.

In the FRW geometry, proper time is related to redshift through the differential ż = -H(z)(1+z). Therefore the time Δ t≡ t_2-t_1 elapsed between any two redshifts z_1 and z_2 is given by the integral

    Δ t    ≡∫_t_1^t_2 dt = H_0^-1∫_z_2^z_1dz/(1+z)ℰ(z),

which has no tractable analytic solution ordinarily, but may be calculated analytically in simple cases or numerically in general. The age of the universe calculated numerically from Eq. (<ref>) is t_0=13.8 Gyr, which was reported in the 2015 Planck results <cit.>.

The spatially-averaged inertial mass densities of the various components of the cosmological fluid decrease as power-laws with decreasing redshift according to their equation of state, i.e.

    ρ_i    ∝ (1+z)^3(1+w_i)

for matter, radiation, curvature, and the cosmological constant. The average matter density in the universe evolves as ρ_m ∝ (1+z)^3. The effective mass density of radiation evolves as ρ_r ∝ (1+z)^4. In the early universe the redshift dependence of this term is more complicated due to the presence of radiation in the form of neutrinos and other relativistic Standard Model (SM) particles in addition to photons. Thus ρ_r  evolves <cit.> according to the expression

    ρ_r(T)    = π^2/30g_⋆(T)k_B^4T^4/c^5ħ^3,

where g_⋆(T) is the effective number of relativistic degrees of freedom, k_B is Boltzmann's constant, and T is the temperature. Thus the effective radiation mass density evolves with redshift as ρ_r∝ g_⋆(1+z)^4.


A list of the important particle mass and energy thresholds is given in Table <ref>. Shown in Fig. <ref> is the corresponding plot of g_⋆ as a function of temperature.

The factor g_⋆ increases by up to a factor 106.75/3.38∼31.6 at high redshift when all SM particles are relativistic. The factor g_⋆(T) for all relativistic particle species in thermal equilibrium can be calculated as the sum <cit.>

    g_⋆(T)    = ∑_i=bg_i(T_i/T)^4 + 7/8∑_i=fg_i(T_i/T)^4,

where the first term is the sum over all bosons and the second term is the sum over all fermions.

Each component of the cosmological fluid has an associated temperature whose value depends on redshift. At high redshift after inflation the universe is dominated by radiation and the components of the cosmological fluid are in equilibrium. The temperature of radiation evolves simply as T_r = T_0(1+z) where T_r is the radiation temperature at any redshift lower than the neutrino decoupling redshift z_dec,ν and the temperature today is T_0 = 2.72548 K <cit.>. The baryonic matter is coupled to the radiation through Compton scattering prior to the thermalization redshift <cit.>

    1 + z_th   ∼ 800(Ω_b,0h^2)^2/5∼ 174,

and in this redshift regime the baryonic matter temperature evolves as T_b = T_0(1+z). After z_th, the baryonic matter temperature evolves adiabatically until reionization as

    T_b = T_0(1+z)^2/1+z_th,

so that at thermalization, T_b = T_r(z=z_th).

It is assumed that the dark matter is mostly weakly interacting massive particles (WIMPs) denoted χ. For simplicity other dark matter models are not considered in this study; see <cit.> for more details on all the dark matter models. For a discussion on the possibility of the PBHs themselves being the dark matter see <cit.>. If the dark matter is composed of WIMPs the DM temperature decouples from the radiation temperature at high redshift when thermal freeze out of the dark matter particles occurs <cit.>. The model forces the DM thermal freeze out to occur when k_BT_fr∼ m_χ c^2/20. If m_χ c^2 = 100 GeV the freeze out redshift is 1+z_fr = T_fr/T_0∼2×10^13.

The root-mean-square velocity v_rms is approximated by comparing the relativistic kinetic energy of the dark matter particles to their thermal energy:

    (γ_χ-1)m_χ c^2    = 3/2k_BT_χ,

where the Lorentz factor associated with β_χ≡ v_rms/c of the dark matter particles is given by

    γ_χ   ≡ (1-β_χ^2)^-1/2.

It is useful to define the dimensionless quantity

    Θ_χ   ≡3k_BT_χ/2m_χ c^2

such that Eq. (<ref>) becomes γ_χ = 1+Θ_χ and therefore

    β_χ^2    = Θ_χ^2+2Θ_χ/(1+Θ_χ)^2,

which is friendly to numerical evaluations for all possible physical values of Θ_χ. The value of β_χ becomes close to unity quite rapidly and for the assumed freeze-out temperature, β_χ≃ 0.367, i.e. the dark matter particles are mildly relativistic at freeze-out.






§  PRIMORDIAL BLACK HOLE ACCRETION AND EVAPORATION

The history of the universe may be divided into redshift regimes to simplify analysis. The relevant physical processes in the very early universe are distinct from those acting in the current era and thus it is important to summarize the physics in each regime.

Though important to the dynamics of the universe in general, the history of the very early universe (prior to inflation) is not considered in detail in this study. The number density of a pre-inflation cosmological relic, e.g. any pre-inflation PBHs, is negligible after the inflationary epoch. The number density of a pre-inflation relic depends on the amount of inflation. The amount of inflation is calculated by finding the number of e-foldings during the inflationary epoch

    N    ≡∫_a_i^a_fdln a = ∫_t_i^t_f H dt,

where t_i denotes the start of inflation, t_f denotes the end of inflation, and H is the Hubble parameter given in Eq. (<ref>). A successful inflation model requires the number of e-foldings to be at least N_min∼50 in order to solve the horizon problem <cit.>. The number density of a relic which formed prior to inflation will thus decrease by a factor e^3N_min∼ e^150∼ 10^65. Thus any PBHs which formed prior to inflation are unlikely to be located in the observable universe today.

After inflation, PBHs may form through a variety of mechanisms including collapse of primordial inhomogeneities, phase transitions, and cosmic string or domain wall collisions <cit.>. If the energy density fluctuations have a strength δρ/ρ∼ 1 in a particular spacetime volume, the region will likely collapse to a black hole. In this study it is assumed that the collapse to a black hole occurs on a time scale much shorter than the Hubble time so that the expansion is irrelevant to PBH formation. A black hole forming at a time t after the Big Bang will have a mass less than or equal to the Hubble mass at that time, i.e.

    M_H   = c^3t/G∼ (4.0×10^14 g)t_-24,

using the useful notation f = 10^nf_n. At t=1.0 s after the Big Bang, the Hubble mass is M_H ∼ 2.0×10^5M_⊙. A derivation of Eq. (<ref>) can be found in <ref>. Models for the mass function of PBHs are discussed in <cit.> with emphasis on the behavior of the mass function near the critical mass regime.

Particles with spin s between energy E and E+dE are emitted near the horizon of a Schwarzschild black hole of mass M at a rate <cit.>

    dṄ   = Γ_sdE/2πħ[exp(8π GME/ħ c^3)-(-1)^2s]^-1,

where Γ_s is the absorption probability for a mode with spin s <cit.>. This is the so-called Hawking radiation which resembles emission from a blackbody with radius R_s = 2GM/c^2 of temperature

    k_BT_BH   = ħ c^3/8π GM = (10.6 MeV)M_15^-1.

The BH can only radiate when the temperature of the hole is greater than that of the radiation bath of the early universe. The temperature of radiation in the early universe evolves as T_r∝ 1+z and is less than 10.6 MeV when z<4.5×10^10, i.e. when t>0.01 s. This will have a negligible effect on the evolution of PBHs near 10^15 g. The absorption probability asymptotes to the geometric optics limit

    Γ_s    = 27G^2M^2E^2/ħ^2c^6

when the particle energy is  E ≫ k_BT_BH. The functional form of Γ_s is much more complicated for lower energy E∼ k_BT_BH interactions as discussed in <cit.>.

The mass loss rate due to the Hawking emission from a Schwarzschild black hole of mass M requires a sum over all particle types and an integration over the particle energies leading to the simple equation 

    dM_15/dt   = (-5.34×10^-5 gs^-1)f(M)M_15^-2,

where f(M) is a function <cit.> allowing for the emission of particles other than photons and f(M)=1 for M≫10^17 g. The function f(M) increases when the mass of the PBH crosses a particle mass threshold (see Table <ref>), after which the PBH may emit that particle. A good approximation is f(M) ∼ f(M_i) because for the majority of its lifetime the mass of a PBH remains near its formation mass M_i <cit.>.

For supermassive black holes or stellar mass black holes the evaporation rate in Eq. (<ref>) is negligibly small. The evaporation rate becomes important on cosmological time scales for black holes with mass M∼10^15 g. This is seen by integrating Eq. (<ref>) to get the evaporation timescale

    t_evap   = (6.24× 10^18 s)f(M_i)^-1M_i,15^3.

Assuming t_evap = t_0 = 13.8 Gyr, the critical formation mass for evaporation today is

    M_cr   = 5.1× 10^14 g,

where the parameter f for the critical mass is f(M_cr)∼1.9 as assumed in <cit.>.

In every cosmic era the Hawking evaporation of a near-critical mass PBH will compete with accretion of the cosmological fluid onto the hole. For the PBHs of M∼ M_cr the accretion turns out to be irrelevant if the hole accretes the cosmological fluid at spatially-averaged densities. For PBHs much smaller than M_cr accretion is completely unimportant. For PBHs larger than M_cr the accretion becomes ever more important and the evaporation rate becomes ever smaller. Thus it is important to quantify the various accretion rates at the relevant cosmic epochs. The accretion rates are dependent on the physical parameters of the cosmological fluid, which change dramatically with redshift. The full equation to be solved is the first-order nonlinear ordinary differential equation in M

    dM/dt   = Ṁ_evap(f,M;z) + Ṁ_acc(ρ,c_s,M;z),

where Ṁ_evap is given by Eq. (<ref>) and the mass accretion term Ṁ_acc will be calculated in the following sections. The equation is integrated from the formation time t_i to any desired final time (or the evaporation time for small holes) using Eq. (<ref>), with the concordance cosmological model accounted for at all times. The Hawking evaporation term has an explicit dependence on f(M) and the mass M of the hole and an implicit dependence on z (or t). The accretion term has explicit dependence on ρ, c_s, and M and an implicit dependence on z due to the evolution of those quantities across cosmic time.

The mass accretion term in Eq. (<ref>) is split into its component parts

    Ṁ_acc(ρ,c_s,M;z)    = Ṁ_r + Ṁ_b + Ṁ_χ,

where r indicates radiation (γ, ν, and other SM particles), b indicates baryonic matter, and χ the dark matter particles. When the universe is cool enough (i.e. T_r < 0.511 MeV), the radiation term consists only of photons and neutrinos. At higher redshift, the other SM particles become relativistic and can be accreted. When the baryonic matter is coupled to the radiation the two accretion rates become coupled and are written Ṁ_b+r. In the sections following, the mass accretion term is calculated explicitly for the different cosmic eras.






 §.§  Late universe accretion

In the late universe at z ≲ z_th the relevant cosmic scales are set by the formation and evolution of structure, i.e. the distribution of dark and baryonic matter in the cosmic web. The details of cosmic structure formation are ripe with rich and complicated physics and are not included in this study; see <cit.>.

To set a bound on late universe accretion, all accretion terms here are set by the spatially-averaged fluid quantities. The PBHs in our universe will likely form and evolve within overdense regions, so the use of spatially-averaged quantities gives a good idea what to expect with relatively isolated holes. The accretion of radiation in this redshift regime is unimportant for ∼ M_cr PBHs because it is a horizon-limited growth given in the approximate form

    Ṁ_r    = 4π R_S^2cρ_r
       = 16π G^2/c^3ρ_r,0(1+z)^4M^2
       = (6.5×10^-48 gs^-1)(1+z)^4M_15^2,

where ρ_r is the equivalent mass density in radiation. The accretion rate in Eq. (<ref>) is comparable to the magnitude of the Hawking evaporation rate when the mass of the PBH is

    M    = (6.3×10^25 g)(1+z)^-1.

Thus the accretion of background radiation in this redshift regime is unimportant to critical mass holes. A PBH of the mass given in Eq. (<ref>) will not evaporate until long after the current era.

The accretion of baryonic matter is more complicated as it is governed by gas dynamics in the vicinity of the PBH. If the PBH is in an `average' region of the universe, i.e. of average baryonic matter density and temperature, the accretion of baryons will be a competition between Bondi and Eddington-limited accretion <cit.>. The accretion rate found in this manner will inform a lower bound for any relevant PBH accretion activity. Below z∼ 30 the details of cosmic structure will change this simplified picture, but it is useful to set a first approximation. A complete picture of the baryonic accretion has not been properly solved and is the subject of intense study from both theoretical and observational perspectives. With the complicated gas dynamics removed from the analysis in this simplified calculation, the accretion rate can be expressed as

    Ṁ_b    = min(Ṁ_b,B,Ṁ_b,E)
       = min(4πλ_sG^2/c_s,b^3ρ_b,0(1+z)^3M^2,4π Gm_p/σ_TcM).

where λ_s = 1/4 for a γ = 5/3 baryonic gas <cit.>, m_p is the proton mass, and σ_T is the Thomson scattering cross section for electrons. It is clear that the Eddington limit is redshift-independent and is equal to

    Ṁ_E    = (7.03×10^-2 gs^-1)M_15.

The Bondi rate is also redshift-independent in this redshift regime. The temperature of the baryonic gas in this regime is

    T_b    = T_r(1+z)/1+z_th = T_0(1+z)^2/1+z_th,

and therefore the sound speed in the baryonic gas (assuming it is entirely hydrogen) is

    c_s,b   = (5k_BT_b/3m_p)^1/2
       = (1+z)(5k_BT_0/3m_p(1+z_th))^1/2
       = (1.5×10^3 cms^-1)(1+z).

Thus the redshift dependence of the Bondi rate goes away and Eq. (<ref>) becomes

    Ṁ_b    = min[π G^2(5k_BT_0/3m_p(1+z_th))^-3/2ρ_b,0M^2,4π Gm_p/σ_TcM]
       = 
    (1.9×10^-24 gs^-1)M_15^2,     M<M_cr,1
    
    (7.0×10^-2 gs^-1)M_15,     M>M_cr,1,

where M_cr,1 = 3.8×10^37 g is the mass of a PBH that gives an equivalence in the Bondi and Eddington rates in this redshift regime. Comparing Eq. (<ref>) to Eq. (<ref>) it is clear that an isolated near-critical mass PBH cannot accrete sufficiently to beat the Hawking evaporation rate. The Bondi rate is comparable to the magnitude of the Hawking evaporation rate when the PBH has the characteristic mass

    M_ch,1   = 7.3×10^19 g.

Any relevant growth of a near-critical mass PBH in this redshift regime will have to come from enhanced accretion if the hole is located within a significant density perturbation such as an individual galaxy or galaxy cluster.






 §.§ Post-recombination accretion

In the post-recombination universe (z_th < z < z_rec) the matter temperature is coupled to the radiation temperature via Compton scattering, i.e. T_b = T_0(1+z). The recombination redshift is listed in <cit.> as z_rec = 1089.90. Starting 381,000 yr after the Big Bang and until thermal decoupling, the sound speed in the baryonic gas can be expressed as

    c_s,b   = (5k_BT_b/3m_p)^1/2
       = (1+z)^1/2(5k_BT_0/3m_p)^1/2
       = (1.9×10^4 cms^-1)(1+z)^1/2.

In this redshift regime, the accretion of radiation is still horizon limited and given by Eq. (<ref>). The accretion of baryonic matter is the Bondi accretion rate at lower mass and is Eddington-limited growth if the mass is large enough. Using the same arguments as before

    Ṁ_b    = min(Ṁ_b,B,Ṁ_b,E)
       = min(π G^2/c_s,b^3ρ_b,0(1+z)^3M^2,4π Gm_p/σ_TcM)
       = 
    (8.1×10^-28 gs^-1)(1+z)^3/2M_15^2,     M<M_cr,2
    
    (7.0×10^-2 gs^-1)M_15,     M>M_cr,2,

where M_cr,2 = (8.7×10^40 g)(1+z)^-3/2 is the PBH mass that gives an equivalent Bondi and Eddington rate. The Bondi accretion rate in Eq. (<ref>) is comparable to the magnitude of the Hawking evaporation rate when the PBH has a characteristic mass

    M_ch,2   = (5.1×10^20 g)(1+z)^-3/8.

Thus in the post-recombination era until thermal decoupling, the relevant process for near-critical mass PBHs is Hawking evaporation.






 §.§  Pre-recombination accretion

In the pre-recombination era (z_rec < z < z_mr) after matter-radiation equality the baryonic matter and radiation are fully coupled and cannot accrete independently. Thus the assumptions present in the Bondi accretion formula fail <cit.> and the accretion of the coupled fluid is horizon-limited. The temperature of the baryonic gas is coupled to the radiation temperature and the sound speed in the fluid can be written (see <ref>)

    c_s^2    = c^2/34ρ_r/4ρ_r+3ρ_b.

At the matter-radiation equality the sound speed in Eq. (<ref>) is a few percent below the asymptotic value c/√(3).

The accretion rate of this coupled fluid onto a PBH is the horizon-limited rate

    Ṁ_b+r   = 4π R_S^2c_s(ρ_r+ρ_b)
       = 16π G^2/√(3)c^3[4ρ_r,0(1+z)^4/4ρ_r,0(1+z)^4+3ρ_b,0(1+z)^3]^1/2
       ==×[ρ_r,0(1+z)^4+ρ_b,0(1+z)^3]M^2.

This rate has a complicated dependence on redshift so it is useful to expand the right hand side of Eq. (<ref>) near the boundaries of this redshift regime. Defining the intermediary terms ρ_r'≡4ρ_r,0 and ρ_b'≡3ρ_b,0 near z_rec the rate takes the form

    Ṁ_b+r   = (6.8×10^-36 gs^-1)[0.26ρ_r'(1+z/1+z_rec)^4
       ==+0.17ρ_b'(1+z/1+z_rec)^3]M_15^2.

The redshift dependence of the sound speed in Eq. (<ref>) is included in the expansion above and in the expansion that follows. The rate in Eq. (<ref>) becomes comparable in magnitude to the Hawking evaporation rate when the PBH has a characteristic mass

    M_ch,3a   ≃ (6.2×10^22 g)[0.26ρ_r'(1+z/1+z_rec)^4
       ==+0.17ρ_b'(1+z/1+z_rec)^3]^-1/4.

So again the Hawking evaporation is most important for critical mass, i.e. Eq. (<ref>), PBHs. Closer to z_mr the rate in Eq. (<ref>) takes the form

    Ṁ_b+r   = (5.5×10^-34 gs^-1)[0.25ρ_r'(1+z/1+z_mr)^4
       ==+0.21ρ_b'(1+z/1+z_mr)^3]M_15^2,

so the accretion of the baryonic radiation fluid occurs slowly for near critical PBHs. This rate becomes comparable in magnitude to the Hawking evaporation rate when the hole is of characteristic mass

    M_ch,3b   ≃ (2.1×10^22 g)[0.25ρ_r'(1+z/1+z_mr)^4
       ==+0.21ρ_b'(1+z/1+z_mr)^3]^-1/4.

The accretion of dark matter onto a PBH will be horizon-limited and should be quite small if the spatially-averaged cosmological value for ρ_χ,0 is assumed. The dark matter accretion rate is

    Ṁ_χ   = 4π R_S^2cβ_χρ_χ,

where β_χ is defined in Eq. (<ref>). The mass density of dark matter evolves according to ρ_χ = ρ_χ,0(1+z)^3 such that Eq. (<ref>) becomes

    Ṁ_χ   = 16π G^2/c^3(Θ_χ^2+2Θ_χ)^1/2/1+Θ_χρ_χ,0(1+z)^3M^2.

If there is an enhancement of the DM density term ρ_χ due to the formation of a DM halo there will be an appropriate enhancement of the DM accretion rate. Thus Eq. (<ref>) represents a lower limit on the DM accretion rate. For a treatment of accretion from an enhanced DM halo see <cit.>. Since the temperature T_χ of dark matter decoupled from the radiation temperature at z_fr∼2.1×10^13, the dimensionless quantity Θ_χ in this redshift regime is quite small. The expansion of β_χ for Θ_χ≪ 1 is β_χ≃ (2Θ_χ)^1/2. Thus Eq. (<ref>) becomes

    Ṁ_χ   = 16π G^2/c^3[3k_BT_0/m_χ c^2(1+z_fr)]^1/2ρ_χ,0(1+z)^4M^2
       = (3.4×10^-58 gs^-1)(1+z)^4M_15^2,

which is about ten orders of magnitude smaller than the accretion rate due to the baryon-radiation coupled fluid. In this regime the accretion rate of dark matter onto a PBH becomes similar to the Hawking evaporation rate when

    M_ch,3c   = (2.3×10^28 g)(1+z)^-1.

The constraints on the accretion rates further strengthens the argument that accretion onto a critical mass PBH is unimportant and most if not all of the lifetime of such a PBH is dominated by the Hawking evaporation.






 §.§ Post-DM freeze-out accretion

In the post-DM freeze-out (z_mr<z<z_fr) era the universe is dominated by radiation. The dark matter, if it comprised of WIMPs, will be non-relativistic until redshifts higher than z_fr <cit.> and will accrete at a horizon-limited rate. The accretion of baryonic matter and radiation is horizon-limited as before. It is convenient to apply ρ_b≪ρ_r and therefore ignore the baryonic matter terms and allow c_s∼ c/√(3). Also in this redshift regime, the effective number of relativistic degrees of freedom g_⋆ begins to increase at higher redshift so it is important to express the radiation term as in Eq. (<ref>). The accretion rate is therefore

    Ṁ_b+r   = 4π R_S^2c_s(ρ_r+ρ_b)
       ≃8π^3 G^2k_B^4T_0^4/15√(3)c^8ħ^3g_⋆(1+z)^4M^2
       = (2.0×10^7 gs^-1)(g_⋆/86.25)(1+z/1+z_fr)^4M_15^2.

For a near critical mass PBH this is a large accretion rate compared to the magnitude of the Hawking rate. Thus the mass of a PBH in this redshift regime where these two rates balance is

    M_ch,4   = (1.5×10^12 g)(g_⋆/86.25)^-1/4(1+z/1+z_fr)^-1.

The period of enhanced accretion in the early universe is quite short due to the strong redshift dependence, i.e. Ṁ∝ (1+z)^4 so no significant accretion is expected for critical mass PBHs. This is consistent with the findings from previous studies on PBH accretion, i.e. <cit.>. At high redshift a critical mass PBH will not accrete significantly, but massive PBHs can grow by about an order of magnitude by z_mr.

In this redshift regime the accretion of DM onto the PBH is small. It is increasingly important at higher redshift but is never larger than the radiation accretion rate in Eq. (<ref>). At the DM freeze-out redshift the DM particles are somewhat relativistic, i.e. Θ_χ∼ 0.075, such that the accretion rate is of the same form as Eq. (<ref>) to a good approximation. In this regime Eq. (<ref>) also remains valid.






 §.§ Pre-DM freeze-out accretion

In this redshift regime the universe undergoes many changes as g_⋆ increases and all particles become relativistic. At high enough redshifts all particles have the same temperature and follow T=T_0(1+z). The accretion rate at these high redshifts is therefore the same as Eq. (<ref>). The PBH will not accrete radiation in the early universe if T_BH>T_r, which corresponds to z<4.5×10^10 if M = 10^15 g. The radiation accretion at these high redshifts is highly dependent on the particle physics model. This study employs the Standard Model with all the latest particle masses from <cit.>. The equivalent mass density in radiation changes dramatically in the early universe because of the change in g_⋆ as shown in Table <ref>.

Table <ref> summarizes the relevant properties of the universe with reference to the equations they are first noted. Table <ref> summarizes the relevant evaporation and accretion rates of PBHs in the relevant redshift regimes with reference to the equations or sections they are first noted. 









§  RESULTS

From the evaporation and accretion expressions in <ref> it is possible to construct a rough accretion or evaporation history for any PBH with mass M_i forming at redshift z_i. The critical mass holes with M_i = M_cr∼ 5.1×10^14 g will suffer no significant accretion in their entire lifetime if located in a suitably `average' volume of the universe. They will assume the evaporation timescale in Eq. (<ref>) and evaporate according to Fig. <ref>.


If the same PBHs of Fig. <ref> happened to form later, say at redshift z_i = 10^8, it would not affect their history due to the small timescales in the early universe. The accretion rate of the cosmological fluid onto larger PBHs at high redshift will be more important.

The analysis in <ref> can be summarized in a look-up plot of M_f against M_i. The regime important for M_i∼ 10^-4M_⊙ holes is shown in Fig. <ref> and the entire mass regime is shown in Fig. <ref>. Note the agreement of Fig. <ref> in the near-critical mass regime to Figure 2 of <cit.>. The holes evaporating at higher redshift must have initial masses slightly lower than M_cr. Note that no significant accretion occurs across the intermediate mass regime between M_cr and ∼10^36 g due to the low accretion rates for BHs of this mass. Isolated PBHs in this mass regime accreting the spatially-averaged cosmological fluid do not grow much. This does not account for enhancement of the accretion rates due to structure formation and thus represents a first approximation. If the accretion rate is enhanced via Ṁ_b,B→ f_bṀ_b,B where f_b = ρ_enh/ρ_b is an enhancement factor and ρ_enh is the enhanced baryonic matter density, then a PBH of given initial mass can reach a higher mass for a given final redshift. This is reflected in the dotted lines of Fig. <ref>, which show the final mass of a PBH growing from 1.0 s after the big bang to z=0.1 given an enhancement factor f_b=10^1,10^2,10^3. Even a small enhancement of the baryonic matter density leads to a large increase in the possible final mass of the accreting PBH. Since the Bondi accretion rate is proportional to M^2, higher mass PBHs will accrete more than lower mass PBHs and this increase in the accretion rate is indicated by the increasing M_f in Fig. <ref> around M∼10^36 g.

In the first few seconds of the universe (z≳10^9), PBHs approaching the formation mass limit around 10^38 g have a large accretion rate (see Eq. <ref>). This large accretion rate, though short-lived, can increase the mass of the PBH by about an order of magnitude by z=10^9. This effect is absent in lower mass PBHs and thus is visible in Fig. <ref> as a small increase beginning above M_i∼10^38 g.










§  CONCLUSIONS AND DISCUSSION

A comprehensive view of the evolution of PBHs throughout cosmic history was presented. The accretion and evaporation histories of PBHs with masses in the approximate range 10^14 g < M < 2×10^38 g were calculated. PBHs with lower masses will have evaporated prior to the current era and are not considered and PBHs with higher masses are not allowed due to the Hubble mass constraint of Eq. (<ref>). The accreted fluids were assumed to have spatially averaged cosmological densities and the details of structure formation were not included. The important quantities for accretion are the mass densities of the various cosmological fluids, the sound speed in those fluids, and the details of their behavior at all relevant redshifts. All of these details were calculated precisely for the ΛCDM concordance cosmology.

The important findings of this study are the following:


  * A PBH with initial mass near M_cr=5.1×10^14 g will not accrete radiation or matter in any significant quantity and will thus evaporate according to the timescale given in Eq. (<ref>). A PBH with initial mass less than M_cr will evaporate prior to the current era.

  * A PBH with initial mass in the approximate range 10^15 g < M_i < 10^35 g neither evaporates nor accretes significantly over a Hubble time. Such a PBH would have to grow by other means, i.e. merging with other BHs or accreting while in an overdense region of the universe. Since the Hawking evaporation rate is so small for PBHs in this mass regime, the lower limit on the final (observed) mass of such PBHs is thus simply M_f = M_i.

  * A PBH with initial mass M < 10^38 g will not grow significantly in the early universe, i.e. within the first few minutes after the Big Bang. This finding is consistent with other PBH accretion studies, e.g. <cit.>. The small increase for BHs with M_i∼10^38 g seen in Fig. <ref> results from the large accretion rate for high-mass holes in Eq. (<ref>). It represents a growth of approximately one order of magnitude in the early universe, consistent with previous studies. There is negligible growth of critical mass PBHs in the radiation-dominated era.

  * A PBH with initial mass in the approximate range 10^35 g < M_i < 10^37 g can accrete significantly during its lifetime. In the redshift regime z_th<z<z_rec, a PBH with M < (8.7×10^40 g)(1+z)^-3/2 accretes at the Bondi rate and is Eddington-limited above that. In the redshift regime z≤ z_th, a PBH with M < 3.8×10^37 g accretes at the Bondi rate and is Eddington-limited above that. A PBH with such a mass that grows at the Bondi rate for its whole lifetime can thus grow by one or two orders of magnitude.

  * When a PBH grows enough to have its baryonic matter accrete at an Eddington-limited rate, the hole can increase in mass by many orders of magnitude if evolving into the late universe z_f∼ 0. Since the PBH will grow by accreting the spatially averaged cosmological gas, this growth represents how an `average' PBH accretes at the Eddington limit. The true accretion history of course will be complicated by feedback effects which were not modeled here. The curves in Fig. <ref> thus represent an `average' growth. A true astrophysical hole of this mass may grow at either a higher or a lower rate.

  * A PBH with initial mass in the approximate range 4×10^37 g < M_i < 10^38 g will accrete at an Eddington limited rate after z_rec and the final mass of such a hole depends on its observed redshift. At z_f=30, the hole can only grow to M_f=1.01M_i. The hole can grow to M_f=7.7M_i if z_f=6. The hole can grow to M_f=6.4×10^11M_i if z_f=0.1 and to M_f=1.1×10^13M_i if z_f=10^-4. See Fig. <ref> for more details.

The PBH mass histories discussed in this study represent a first approximation of their cosmic behavior. Several astrophysical applications may be discussed in the context of the above results:


  * It is impossible to explain the large BHs with M∼10^10M_⊙ observed <cit.> at z>6 via PBHs with Eddington-limited accretion of the `average' baryonic gas, even with M_i∼ 10^5M_⊙. These holes must be explained through multiple massive PBH mergers, mergers with BH seeds from the first generation of stars, or PBHs in overdense regions accreting at super-Eddington rates.

  * PBHs do not easily grow to 30M_⊙ by z∼0.1 through Bondi accretion of the `average' cosmological fluid. These PBHs cannot easily explain the binary BH mergers observed by LIGO as the gravitational wave events GW150914 <cit.> and GW151226 <cit.> unless they experience an enhancement of the Bondi rate through various channels. One such channel is a baryonic matter density enhancement leading to Ṁ_b,B→ f_bṀ_b,B as discussed in 4. Small enhancement factors allow a lower mass PBH to reach 30M_⊙ compared to those PBHs accreting the average cosmological baryonic matter. Another possibility is the LIGO BHs were PBHs that formed with an initial mass M_i=M_f, where M_f is their mass at the merger time. According to <cit.>, the event rate for PBH mergers would be high enough to explain the GW events if the PBHs constitute a large enough fraction of the dark matter. However, PBHs in the appropriate mass range to explain these LIGO events are unlikely to be a large enough fraction of the DM as constrained from CMB measurements discussed in <cit.>, <cit.>, and <cit.>. Either LIGO has chanced upon two relatively rare PBH mergers or there is a common stellar evolution channel that produces BHs of these masses. Both explanations are interesting and more data are needed to distinguish these two possibilities.

  * Searches for PBH bursts <cit.> are ongoing. Although there are candidates for such events, no confirmed PBH burst event has been detected. The spectral properties of such bursts should be distinguishable from the `normal' GRBs. The non-detection of such an event has a few explanations. First, the fraction of PBHs that make up the dark matter must be quite low for PBHs of the relevant mass scale (see Fig. 9 of <cit.>). Thus it is plausible that not enough of these PBHs exist to subsequently evaporate and trigger gamma-ray detectors. Second, it might be possible for the critical mass holes to accrete enough to no longer evaporate in the current era. However, the accretion rate is too small and this would not explain the non-detection of PBH bursts. Even if the accretion rate onto small PBHs happened to be large enough, there would be smaller PBHs that would accrete enough to reach M_cr anyway, filling the void of critical mass holes.







§ ACKNOWLEDGMENTS

This work was supported by the Nevada NASA Space Grant College and Fellowship Training Program Cooperative Agreement #NNX10AN23H. We thank the anonymous referee for multiple suggestions that significantly improved the quality of the manuscript. J.R.R. thanks the UNLV High Energy Astrophysics Group for fruitful discussions. We extend our gratitude to Dr. Qing-Guo Huang of the Institute of Theoretical Physics, Chinese Academy of Sciences who alerted us to an error in an equation that is now fixed.



§  MAXIMUM PBH FORMATION MASS

As discussed in <ref> the maximum formation mass of a PBH will be the Hubble mass, i.e. the mass contained within the Hubble volume at a given time. The Hubble volume is

    V_H    = 4/3π R_H^3,

where R_H = 2ct is the Hubble radius in the early radiation dominated universe. The critical density is the Hubble mass in a Hubble volume and thus

    M_H(t)    = ρ_crV_H
       = 3H^2/8π G·4/3π(2ct)^3
       = 4H^2c^3t^3/G,

but in a radiation dominated universe H=1/(2t) so that

    M_H    = c^3t/G,

which recovers Eq. (<ref>). See <cit.> for more detailed discussions of the PBH mass function.



§  SOUND SPEED IN THE COSMOLOGICAL FLUID

In the late universe at redshifts lower than z_th, the temperature of the baryonic matter decouples from the CMB photon temperature. Thus the sound speed in the baryonic fluid is given by Eq. (<ref>)

    c_s,b   = (1.5×10^3 cms^-1)(1+z),

where the increase due to reionization around z∼9 is not taken into account. In the redshift regime z_th<z<z_rec the redshift dependence changes due to the temperature coupling between the baryonic matter and the CMB radiation. Thus the sound speed evolves as Eq. (<ref>)

    c_s,b   = (1.9×10^4 cms^-1)(1+z)^1/2.

In the above equations it is assumed that the baryonic matter is composed entirely of hydrogen; corrections due to the helium and metal content of the baryonic matter need to be made for a more realistic calculation.

In the early universe at redshifts higher than the recombination redshift z_rec∼1090, the baryonic matter is coupled to the CMB radiation. The sound speed in such a coupled fluid can be found by calculating

    c_s^2    = (∂ P/∂ρ)_s,

where the subscript s on the right hand side indicates taking the derivative at constant entropy. The dominant pressure term is the radiation pressure and the density is a sum of radiation and baryonic terms ρ = ρ_r + ρ_b. The dark matter does not contribute to the pressure or density terms but has an early influence when it is relativistic at redshifts greater than z_fr∼2.1×10^13.

Rewriting the partial derivatives of Eq. (<ref>) in terms of temperature gives

    c_s^2    = (∂ P_r/∂ T)_s/(∂ρ_r/∂ T)_s+(∂ρ_b/∂ T)_s.

Recalling Eq. (<ref>) and P_r = ρ_r c^2/3 the numerator of Eq. (<ref>) is

    (∂ P_r/∂ T)_s    = 4π^2/90g_⋆(T)k_B^4T^3/c^3ħ^3 = 4ρ_rc^2/3T, ignoring the small ∂ g_⋆/∂ T terms. Similarly, the first term in the denominator of Eq. (<ref>) is(∂ρ_r/∂ T)_s    = 4π^2/30g_⋆(T)k_B^4T^3/c^5ħ^3 = 4ρ_r/T.

Recalling at high redshift the radiation and baryonic gas temperatures are coupled, i.e. T_r = T_b = T and using T = T_0(1+z), the second term in the denominator of Eq. (<ref>) is

    (∂ρ_b/∂ T)_s    = ∂/∂ T(ρ_b,0T^3/T_0^3) = 3ρ_b/T.

Combining Eq. (<ref>–<ref>) into Eq. (<ref>) gives

    c_s^2    = c^2/34ρ_r/4ρ_r+3ρ_b.

It is clear that at redshifts higher than z_mr∼3400 the sound speed calculated using Eq. (<ref>) asymptotes to c_s∼ c/√(3)∼ 0.577c. The behavior of the sound speed across all relevant redshifts is plotted in Fig. <ref>.




§ REFERENCES



