Tracing Primordial Black Holes in Nonsingular Bouncing Cosmology

We in this paper investigate the formation and evolution of primordial black holes (PBHs) in nonsingular bouncing cosmologies. We discuss the formation of PBH in the contracting phase and calculate the PBH abundance as a function of the sound speed and Hubble parameter. Afterwards, by taking into account the subsequent PBH evolution during the bouncing phase, we derive the density of PBHs and their Hawking radiation. Our analysis shows that nonsingular bounce models can be constrained from the backreaction of PBHs.


Introduction
The matter bounce scenario [1][2][3] is one type of nonsingular bounce cosmology [4][5][6][7][8][9], which is often viewed as an important alternative to the standard inflationary paradigm [10][11][12][13]. By suggesting that the universe was initially in a contracting phase dominated by dust-like fluid (with a vanishing equation-of-state parameter w = 0), then experienced a phase of nonsingular bounce, and afterwards entered a regular phase of thermal expansion. The matter bounce cosmology can solve the horizon problem as successfully as inflation and match with the observed hot big bang history smoothly. Based on primordial fluctuations generated during matter contracting and their evolution through the nonsingular bounce, one can obtain a scale invariant power spectrum of cosmological perturbations. Unlike inflation, the matter bounce scenario does not need a strong constrain on the flatness of the potential of the primordial scalar field that drives the evolution of the background spacetime [14,15]. Also, this scenario can avoid the initial singularity problem and the trans-Planckian problem, which exists in inflationary and hot big bang cosmologies [16,17].
The aforementioned scenario has been extensively studied in the literature, such as the quintom bounce [18,19], the Lee-Wick bounce [20], the Horava-Lifshitz gravity bounce [21][22][23], the f (T ) teleparallel bounce [24][25][26], the ghost condensate bounce [27], the Galileon bounce [28,29], the matter-ekpyrotic bounce [30][31][32], the fermionic bounce [33,34], etc. (see, e.g. Refs. [1,35] for recent re-* Corresponding author. views). In general, it was demonstrated that on length scales larger than the time scale of the bouncing phase, both the amplitude and the shape of the power spectrum of primordial curvature perturbations can remain unchanged through the bouncing point due to a no-go theorem [36,37]. A challenge that the matter bounce cosmology has to address is how to obtain a slightly red tilt on the nearly scale invariant primordial power spectrum. To address this issue, a generalized matter bounce scenario, which is dubbed as the -Cold-Dark-Matter ( CDM) bounce, was proposed in [38] and predicted an observational signature of a positive running of the scalar spectral index [39,40].
As a candidate describing the very early universe, the matter bounce scenario is expected to be consistent with current cosmological observations and to be distinguishable from the experimental predictions of cosmic inflation as well as other paradigms [7,41]. Meanwhile, a possible probe of primordial black holes (PBHs) may offer a promising observational approach to distinguish various paradigms of the very early universe [42,43]. PBHs could form at very early times of the universe, where a large amplitude of density perturbations would have obtained. Correspondingly, the formation process and the abundance of PBHs strongly depend on those early universe models, in which fluctuations of matter fields are responsible for such large amplitudes of density perturbations [44].
In the literature, most of attentions were paid on the computation of PBH predictions from the inflationary paradigm (for instance see [45][46][47][48][49]), while so far, only a few works addressed the PBH formation in a bouncing scenario [50,51]. Furthermore, those studies of PBHs in a bouncing scenario have not yet been discussed in detail, for specific cosmological paradigms or been applied to http falsify various early universe cosmologies, especially the matter bounce scenario. In the context of matter bounce cosmology, there are several differences on the computation of the PBH abundance comparing with that in an expanding universe. First, comparing with inflation where the primordial fluctuations become frozen at the moment of the Hubble exit, those primordial fluctuations on matter fields in bounce cosmology would continue to increase after the Hubble exit during the contracting phase until the universe arrive at the bouncing phase [9,20], and the contracting phase would yield a different initial condition for the PBH formation and evolution. Second, once these PBHs have formed, the contraction of spacetime could also compress and enlarge the primordial matter density, thus change the PBH horizon radius which then can lead to effects on their evolution.
In this paper, we perform a detailed survey on the PBH formation and evolution in the background of the matter bounce cosmology. In Section 2, we briefly introduce the matter bounce scenario and describe the formation of the power spectrum of primordial curvature perturbation in an almost model-independent framework. In Section 3, a physical picture of the PBH formation in the contracting background is presented. After a process of detailed calculations, the threshold for forming PBHs and the corresponding mass fraction are provided. In Section 4, we discuss the evolution of PBHs in the bouncing phase by taking into account the effects arisen from the contraction of the background and the Hawking radiation. In Section 5, we summarize our results and discuss on some outlook of the PBH physics within the nonsingular bouncing cosmology.

Nonsingular bounce
Nonsingular bounce can be achieved in various theoretical models, namely, to modify the gravitational sector beyond Einstein, to utilize matter fields violating the Null Energy Condition (NEC), or in the background of non-flat geometries (see e.g. [52,53]). It is interesting to notice that, in general, on length scales larger than the time scale of the nonsingular bouncing phase, primordial cosmological perturbations remain almost unchanged throughout the bounce [36,37]. In this regard, one expects that the effective field theory approach should be efficient to describe the information of a nonsingular bounce model at background and perturbation level. Recently, it was found in [30] that a nonsingular bounce model can be achieved under the help of scalar field with a Horndeski-type, non-standard kinetic term and a negative exponential potential. Within this model construction, the matter contracting phase can be obtained directly by including the dust-like fluid or involving a second matter field [31]. Note that, in the present study we assume that the effective field approach of bouncing cosmology is valid through the whole evolution without modifications to General Relativity.

The model
It turns out that, under the description of the effective field theory approach, the background dynamics of the nonsingular bouncing cosmology can be roughly separated into three phases: the matter-dominated contraction, the non-singular bounce, and the thermal expansion. We consider a simple model starting with a matter contracting phase (t < t − ) from an initial time t initial → −∞, and then entering into a nonsingular bouncing phase at t − , which lasts till t + . After the bounce ends at t + , the universe begins the hot big bang expansion, which is in accordance to the current observations.
The evolution of the matter bounce cosmology in each stage can be approximately described as follows.
(i) In the matter contracting phase, the scale factor of the universe shrinks as where a − is the scale factor at time t − , and t − is related to the Hubble parameter at t − via the relation t − −t − = 2 3H − . The equation-of-state parameter during this phase is w = 0, which can be realized in many ways, such as by cold dust, by massive field or by the gravity sector involving non-minimal couplings. We parameterize these different mechanisms by introducing the sound speed c s , which can affect the propagation of primordial perturbations in the gradient terms.
(ii) In the nonsingular bouncing phase, the scale factor of the universe can be approximately described as [30] a(t) = a B e  into the Friedmann equation, one can see that the null energy condition ρ + p > 0 is violated around the bouncing point. It is the negative pressure that avoids the singularity and drives the universe to evolve from a contracting phase to an expanding phase.
(iii) In the era of radiation-dominated expansion, we have In present analysis we have adopted the assumption that the heating process happens instantly after the bounce (see [54,55] for relevant analyses).
From such parameterizations, the matter bounce cosmology can be approximately described by model parameters H − , H + , ϒ , and also c s if perturbations are taken into account. In Fig. 1 we depict the evolution of comoving Hubble length |H −1 | = |aH| −1 , and one can read that one Fourier mode of cosmological perturbation in the matter bounce cosmology could exit the Hubble radius during the contracting phase, and then enter and re-exit the Hubble radius during the bouncing phase, and eventually re-enter the Hubble radius again in the classical Big Bang era. Note that the change of the Hubble radius in the vicinity of the bouncing point can be very large. In the literature, observational constraints upon the bounce cosmology can be derived from various cosmological experiments such as the cosmic microwave background [41] and primordial magnetic fields [56]. In this work we will provide the independent constraints on the model parameters from PBHs.

Curvature perturbation during matter contracting phase
During the matter contracting phase, equation of motion for the curvature perturbation can be expressed as where v k = zR k is the Mukhanov-Sasaki variable with R k being the comoving curvature perturbation, z = a √ ρ+p c s H relies on the detailed evolution of the background dynamics, and the prime represents the derivative with respect to the comoving time η.
Assuming that primordial perturbations originated from vacuum fluctuations at initial times, 1 one can derive the solution of Eq.
where H 3/2 is the 3 2 -order Hankel function of the first kind. From Eq. (2.6), the power spectrum of R k is then given by during the matter contracting phase. At large scales c s k is almost scale independent, and R k This result is different from that in an expanding spacetime, in which R k is time independent at large scales.
Moreover, at small scales c s k |H|, curvature perturbation is of

PBH formation during matter contracting phase
PBHs originate from the collapsed over-dense regions seeded by cosmological perturbations in early universe. When an over-dense region starts to form a black hole (BH), its radius R is expected to be larger than the Jeans radius R J ≡ c s 2 π ρ = 2 3 πc s |H −1 | and smaller than the Hubble length |H −1 | [42,44,48,58,59], which implies, 2 3 πc s |H −1 | ≤ R ≤ |H −1 | . (3.1) For scales smaller than R J , the pressure gradient force would prevent the collapsing of the over-dense region. It is clear that R J ≥ 1 Note that the vacuum state is not necessarily the only choice of the initial condition for primordial curvature perturbations, namely, they may also arise from fluctuations of a thermal ensemble [20,57].
|H −1 | when c s ≥ 0.39, and hence, there would be no PBH in this case.
We note that the PBH formation in the contracting phase is different from that in a purely expanding universe. In the expanding universe, the PBH formation would begin only when the size of a region reenters the Hubble scale. However, in the contracting phase, when the size of the region exits the Jeans scale, the PBH formation would have started if this region is dense enough. Thus, we would like to focus on the fluctuations within the Jeans scale instead of the Hubble scale in the contracting universe.

PBH mass fraction
The background dust-like fluid inside the over-dense regions may collapse into PBHs, and the abundance of PBHs is depicted by the mass fraction, which in the matter contracting phase is (3.2) where ρ PBH denotes the density of PBHs and ρ bg the density of the remnant dust-like fluid. In general, PBHs and remnant dustlike fluid together drive the evolution of background ρ PBH + ρ bg = 3M 2 p H 2 . Note that this relation only holds when the subsequent evolutions, including the accretion and Hawking radiation, of the formed PBHs are not considered.
The mass fraction can be obtained from the Press-Schechter theory [60] (see [42,44,48] also): where δ ≡ δρ ρ is the fractional density fluctuations, σ ≡ δ 2 is the mean mass fluctuation at the Jeans scale, which can be derived from the power spectrum of curvature perturbation 2 R , δ c is the threshold of the PBH formation, and the upper limit δ m ensures that the PBHs are no larger than the Hubble scale [42,44,48]. Note that Eq. (3.3) is based on the assumption that δ obeys the Gaussian The value of threshold δ c is determined as follows. Considering a spherical over-dense region with a radius R, the space inside the region satisfies Friedmann equation [61,62] where H is the Hubble parameter of the background, a is the scale factor inside the over-dense region, δ H and δ K are the perturbed Hubble parameter and curvature respectively, and δ ≡ δ − δ K a 2 H 2 . It is noticed that H < 0 for a contracting phase, and δ H < 0 in the over dense region. The outer region (> R) is thought to be unperturbed for simplicity. We assume that when the region collapses to a BH, its surface has an additional physical speed v = 1 with respect to the conformally static background, i.e. collapsing at the speed of light. In the comoving slicing, it is written as δ H · R = −1 .  The mass of BH with a radius R, in the contracting background, can be determined by Note that, both the density and curvature fluctuations could contribute to the BH mass, as shown in [48]. We mention that the above description of BH is a rough estimate, and hence, Eq. (3.7) slightly differs from the Misner-Sharp mass M = R 2G . The accurate description of the BH in the matter contracting phase should be based on the Tolman-Bondi-Lemaitre metric [63][64][65], which will be addressed in our follow-up study. In present analysis, however, we note that at small scale R |H −1 | the above model naturally recovers the solution of a Schwarzschild BH M = R 2G in a flat spacetime, which indicates that this estimate remains reliable.
To derive the values of δ c and σ , one needs to know the relations among the variables δ, δ K and R k . It is convenient to take the comoving gauge, in which where (x, t) is the Bardeen potential and R(x, t) is the curvature perturbation [45,66,67]. The relation between their Fourier components k and R k is given by [45,67] − , (3.9) where the dot denotes the derivative with respect to the cosmic time t. The solution during the matter contracting era R k ∝ a −3/2 differing from k = − 3 5 R k in the matter dominant expanding era.
One also obtains (3.11) where δ k is the Fourier component of the density perturbation δ. From (3.6) and (3.8), one also has δ = 3δ K 2a 2 H 2 = 3δ. Therefore, the threshold is δ c = 3 where the first term on the left side is the mass of the overdense region and the second term is the mass outside the R J but inside the Hubble radius. As a result, we obtain δ m = 3δ m = (3.13) Accordingly, σ is determined by 2 In this case, the PBH is formed by an over-dense core with radius R J , and a shell of the background fluid within the radius from R J to |H −1 |.  (3.14) where W kR J a is the window function. For simplicity, we take for c s < 0.39; and β(t) = 0 for c s ≥ 0.39. It can be seen that the value of β increases as a(t) becomes smaller along with the contracting and reaches its maximum value at the end of the matter contracting phase t − . From the above expression, one can read that the maximal value of β is unity, which corresponds to the case that the universe is dominated by black holes. Accordingly, we expect that β 1 so that the universe is dominated by the background dust-like fluid. By numerically solving Eq. matter-radiation equality (see for example [68] for related analyses). Therefore, if one considers the case of β 0.1 and c s = 10 −4 as a specific example, then Eq. (3.16) yields |H − | 10 −1 M p or |H − | 10 4 M p . Note that, |H − | 10 −1 M p can be easily satisfied, but |H − | 10 4 M p is disfavored from model construction.

Evolution of PBH in bouncing phase
It is important to point out that the subsequent evolutions of PBHs after their formation are not considered in the above section. Associated with the evolutions, including the accretion and the evaporation due to Hawking radiation, one can obtain constraints upon bouncing cosmology. This is the main subject in this section. For simplicity, we only consider the PBH evolutions during the bouncing phase.

The growth of PBH mass
At first, we calculate the growth of PBH by accreting the mass around in the contracting background. For simplicity, we still assume that the space outside the BH horizon is unperturbed Friedmann universe and the matter which flows into the horizon is due to the cosmic contracting. For a BH with mass M and radius R, the mass increased equals to that flowing into the horizon, with the speed |H|R, We shall investigate the evolution in the contracting era of the bouncing phase t − < t < 0. During this stage, the Hubble length diverges quickly (see Fig. 1), and all PBHs can be treated as Schwarzschild BHs and R = 2G M, according to the arguments in the last section. Therefore, Eq. (4.1) reduces tȯ (4.2) and the solution is . (

4.3)
It is obvious that, the BH radius increases with t approaching to the bouncing point. Eq. (4.3) also gives a constraint as follows, For the models violating the constraints above, the mass of PBH will grow to infinity before the bouncing point. According to (3.1), one can read that R(t − ) ≥ 2 3 πc s |H −1 − |, and hence, the above constraint yields explicitly that 3 for c s ≤ 0.39. Eq. (4.5) implies that the bouncing phase is expected to proceed rapidly so that the mass of PBHs would not become divergent and the universe is safe from collapsing into the black hole completely.

Backreaction and theoretical constraints
The bouncing phase is driven by the background fluid which violates NEC in the frame of general relativity. From the Friedmann equations not considering PBHs, the density and pressure of the background fluid at the bouncing point are ρ bg = 0 and p bg = −2M 2 p ϒ respectively. To include the contribution of backreaction from PBHs, the Friedmann equation at the bouncing point can be written as (4.6) where ρ BR = ρ PBH + ρ γ + 3p γ is the effective density of the backreaction, ρ PBH is the density of PBH, ρ γ and p γ are the density and pressure of Hawking radiation respectively. For simplicity, we assume all Hawking radiation particles are ultra-relativistic, which satisfy the relation p γ = ρ γ /3. If the back reactions neutralize the negative pressure of background ρ BR − 6M 2 p ϒ > 0, the universe would not expand any more, as mentioned in Section 2. In the following, we will constrain the model through PBHs.
Recall that the energy density of PBHs is which depends on the mean PBH mass M or the mean PBH radius R , and the mean PBH separation L. At the beginning of the bouncing phase t − , the PBH density is known PBHs. Here we assume that all PBHs have the same mass, which is a simplification generally used [49], and all PBHs are of the Jeans At the present investigation, we have neglected the newly formed PBHs during the bouncing phase. Therefore, the mean separation follows L(t) ∝ a(t), and accordingly, Note that the PBH evaporation due to Hawking radiation has been studied clearly in [43]. Following this, one can estimate the lifetime of PBH to be τ 13.7 Gyr  , the Hawking radiation can be ignored and ρ γ (0) = 0. In the limit that PBHs finish their evaporation at t − , one has ρ γ (t − ) = ρ PBH (t − ) and then ρ γ ∝ a −4 . In the limit that the effective evaporation happens at the bouncing point, ρ γ (0) = ρ PBH (0), and ρ PBH ∝ a −3 before its evaporation. Therefore, one has (4.12) where ρ PBH (t − ) is given by Eq.  the asymptote for all isolines of ρ BR /|p|. Therefore, the constraint ρ BR /|p| < 1 reduces to Eq. (4.5) at the low energy.

Estimates of observational constraints
After the bouncing phase, PBHs will evolve in the regular expanding universe, and hence can be constrained by today's observations. In this subsection, we briefly estimate the constraints upon the bouncing models though the observational limits of PBHs.
For massive PBHs with an initial mass M > 10 15 g, they could survive until today and their density scales as a −3 approximately in the expanding universe [43]. As a result, the current PBH (M >  (4.14) For the light PBHs with M < 10 15 g, they have been completely evaporated today. Therefore, the density parameter of the Hawking radiation HR is approximately estimated as HR ρ PBH (0) where a r is the scale factor when PBHs complete the evaporation, and a b < a r < a 0 . HR is not allowed to be larger than the total density of today's radiation r ∼ 10 −4 , measured by the CMB experiments. As a consequence, one gets In Fig. 4, we provide the constraints on ρ PBH (0) along with a b /a 0 from the bounds (4.14) and (4.16), respectively. The PBH mass M and density ρ PBH (0) are functions of c s , H − and ϒ , which are provided in the previous subsection. As a result, we find that the nonsingular bounce models can be constrained by observations due to (4.14) and (4.16). We note that these constraints are quite loose. More stringent constraints can be obtained by taking the observational bounds of PBH and HR from the CMB, gamma-ray burst and BBN etc. [43].

Conclusion and outlook
In this paper, we have investigated the formation and evolution of PBHs in the matter bounce cosmology. Firstly, we described the general matter bounce models by some parameters like c s , H − and ϒ . The comoving curvature perturbation R k is also calculated during the matter contracting phase, which seeds the PBH formation. Then we had a discussion about the condition of PBH formation in the contracting background, which is different from that in the expanding universe. By taking a simple collapsing model, the threshold of the density fluctuation for forming a PBH is derived. Furthermore, in the comoving gauge, the density fluctuation and its threshold are all related to the curvature perturbation R k . Therefore, we can calculate the mass fraction of PBHs β in the contracting phase from the Press-Schechter theory, and constrain the bouncing models from PBHs. PBH formation depends on the model parameters c s and H − . For instance, PBHs can form in the contracting phase only if c s ≤ 0.39, since the BHs are not allowed to be larger than the cosmic apparent horizon. When c s ≤ 0.39, β is small and the model is safe for |H − | M p or |H − | M p , but only the low energy regime is favored from the model construction.
The subsequent evolution of PBHs in the bouncing phase is also investigated, with the PBH accretion and Hawking radiation considered. The growth behavior of PBH yields a constraint to the model ϒ ≥ c 2 s π 2 H 2 − , in case that the mass of PBHs grow to infinity before the bouncing point. Moreover, the back reaction of PBH and its Hawking radiation is calculated to constrain models, in order not to neutralize the negative pressure p bg = −2M 2 p ϒ , since in such case the universe cannot expand again. The constraint reduces to ϒ ≥ c 2 s π 2 H 2 − at the low energy scale H 2 − 10 9 c 5 s M 2 p . Afterwards, a rough constraint of bouncing model though the PBH observations is given, and the constraint is stringent only when ϒ is slightly larger than c 2 s π 2 H 2 − . A more precise analysis taken into account both the PBH growth and the Hawking evaporation as well as the detailed constraints from cosmological observations will be addressed in our follow-up works.
We note that while our paper was being prepared, an independent work was being carried out by another group [59], which explores similar features of BH formation in bouncing cosmology.