Particle Creation in Bouncing Cosmologies

# Particle Creation in Bouncing Cosmologies

Diogo C. F. Celani    Nelson Pinto-Neto    Sandro D. P. Vitenti Centro Brasileiro de Pesquisas Físicas, Rua Dr. Xavier Sigaud 150
22290-180, Rio de Janeiro – RJ, Brasil
July 19, 2019
###### Abstract

We investigate scalar particle creation in a set of bouncing models where the bounce occurs due to quantum cosmological effects described by the Wheeler-DeWitt equation. The scalar field can be either conformally or minimally coupled to gravity, and it can be massive or massless, without self interaction. The analysis is made for models containing a single radiation fluid, and for the more realistic case of models containing the usual observed radiation and dust fluids, which can fit most of the observed features of our Universe, including an almost scale invariant power spectrum of scalar cosmological perturbations. In the conformal coupling case, the particle production is negligible. In the minimal coupling case, for massive particles, the results point to the same physical conclusion within observational constraints: particle production is most important at the bounce energy scale, and it is not sensitive neither to its mass nor whether there is dust in the background model. The only caveat is the case where the particle mass is larger than the bounce energy scale. On the other hand, the energy density of produced massive particles depend on their masses and the energy scale of the bounce. For very large masses and deep bounces, this energy density may overcome that of the background. In the case of massless particles, the energy density of produced particles can become comparable to the background energy density only for bounces occurring at energy scales comparable to the Planck scale or above, which lies beyond the scope of this paper: we expect that the simple Wheeler-DeWitt approach we are using should be valid only at scales some few orders of magnitude below the Planck energy. Nevertheless, in the case in which dust is present, there is an infrared divergence, which becomes important only for scales much larger than today’s Hubble radius.

###### pacs:
04.62.+v, 98.80.-k, 98.80.Jk

## I Introduction

In cosmological bouncing scenarios, the initial singularity present in the standard cosmological model is absent: the Universe contracts from a very large size, bounces when it becomes sufficiently small, and expand afterwards. The bounce can occur once or many times, depending on the model. This fact (the absence of singularities) is, per se, already very important, and a sufficient reason to investigate these models more deeply. Furthermore, other puzzles of the standard cosmological model are absent in these scenarios (as the horizon and flatness problems), and they can also supply a mechanism to generate primordial cosmological perturbations from quantum vacuum fluctuations, with almost scale invariant spectrum Peter et al. (2007); Peter and Pinto-Neto (2008), as in inflationary models. Hence, they can also be viewed as alternatives to inflation, although they are not necessarily contradictory to it.

There are nowadays many mechanisms to generate the bounce, and there are also many open questions and issues to be investigated concerning these models, for reviews see Refs. Novello and Bergliaffa (2008); Battefeld and Peter (2015). One of them concerns the phenomenon of quantum particle creation around the bounce: is it large enough to modify the background and/or induce some sort of reheating, making the model asymmetric around the bounce, or it is always negligible?

The aim of this paper is to investigate particle creation in the set of bouncing models that have been investigated by some of us in the last decades. In these models, the bounce occurs due to quantum cosmological effects when the curvature of space-time becomes very large: the related Wheeler-DeWitt equation is solved, and interpreted using the de Brogli-Bohm quantum theory (the usual Copenhagen point of view cannot be used in quantum cosmology, see Ref. Pinto-Neto and Fabris (2013) for a review on this subject). The trajectories describing the scale factor evolution are calculated and they are usually non-singular, presenting a bounce due to quantum effects at small scales, and turning to a classical standard evolution when the scale factor becomes sufficiently large Acacio de Barros et al. (1998); Peter et al. (2007); Alvarenga et al. (2002); Pinto-Neto and Fabris (2013). These models usually contain one single hydrodynamical fluid, or two fluids Pinto-Neto et al. (2005): the usual observed radiation and dust contents which are present in our Universe. We expect that the simple Wheeler-DeWitt approach we are using should be valid only at scales some few orders of magnitude below the Planck energy, being a limit of some more involved theory of quantum gravity suitable for energy scales close or above the Planck scale (see Refs. Ashtekar (2013, 2014) as examples of approaches in this direction).

The paper is organized in the following way: in the next section (Sec. II), we present the essential properties of the bouncing models we are considering, and the relevant quantities for the particle creation computation. In Sec. III, we calculate the production of particles in bouncing models with a single radiation fluid, for the conformal and minimal coupling case, and for massive and massless scalar fields. In Sec. IV, we will perform the same calculations for bouncing models with two fluids, dust and radiation. We conclude in Sec. V with a summary of our results, and a discussion concerning its perspectives. All numerical calculations performed in this work were done using the Action Angle (AA) variables method as described in Appendix A. Moreover, we also used the Wentzel-Kramers-Brillouin (WKB) approach, described in Appendix B. Both methods gave the same results within the required precision.

## Ii Particle Production in Bouncing Models

We will consider bouncing models where the physical effect which avoid the usual classical cosmological singularity are quantum corrections to the classical Friedmann equations. These quantum effects are calculated from a Wheeler-DeWitt quantization of the background and interpreting the solution using the de Broglie-Bohm quantum theory, see Ref. Pinto-Neto and Fabris (2013) and references therein.

For a single hydrodynamical fluid with , the Wheeler-DeWitt equation reads

 i∂Ψ(0)(a,T)∂T=14∂2Ψ(0)(a,T)∂χ2, (1)

where we have defined

 χ=23(1−λ)−1a3(1−λ)/2,

is the scale factor, and is the degree of freedom associated to the fluid, which plays the role of time.

This is just the time reversed Schrödinger equation for a one dimensional free particle constrained to the positive axis. As and are positive, solutions which have unitary evolution must satisfy the condition

 Ψ⋆(0)∂Ψ(0)∂χ∣∣∣χ=0−Ψ(0)∂Ψ⋆(0)∂χ∣∣∣χ=0=0 (2)

(see Ref. Alvarenga et al. (2002) for details). We will choose the initial normalized wave function

 Ψ(init)(0)(χ)=(8Tbπ)1/4exp(−χ2Tb), (3)

where is an arbitrary constant. The Gaussian satisfies condition (2).

Using the propagator procedure of Refs. Acacio de Barros et al. (1998); Alvarenga et al. (2002), we obtain the wave solution for all times in terms of :

 Ψ(0)(a,T)=[8Tbπ(T2+T2b)]1/4exp[−4Tba3(1−λ)9(T2+T2b)(1−λ)2]exp{−i[4Ta3(1−λ)9(T2+T2b)(1−λ)2+12arctan(TbT)−π4]}. (4)

Due to the chosen factor ordering, the probability density has a non trivial measure, and it is given by . Its continuity equation coming from Eq. (1) reads

 ∂ρ∂T−∂∂a[a(3λ−1)2∂S∂aρ]=0, (5)

which implies in the de Broglie–Bohm quantum theory that

in accordance with the classical relations and .

Inserting the phase of (4) into Eq. (6), we obtain the Bohmian quantum trajectory for the scale factor:

 a(T)=ab[1+(TTb)2]13(1−λ). (7)

Note that this solution has no singularities, and tends to the classical solution when . We are in the time gauge , thus is related to conformal time through

Although this solution describe a bounce for any choice of barotropic fluid, in this work we are interested in a bouncing scenario without any exotic component. Hence, the matter component with largest equation of state that will dominate the evolution for small will be radiation. Thus, choosing , we obtain that , while is simply the conformal time, i.e., .

The solution (7) can be obtained for other initial wave functions (see Ref. Alvarenga et al. (2002)). In Ref. Peter et al. (2016) it was shown that defining the wave function at an arbitrary initial time also leads to bouncing solutions. For the symmetric bounce given in Eq. (7), the quantum effect which causes the bounce is equivalent to adding to the Friedmann equation a term corresponding to a negative energy density of a stiff matter fluid. This same behavior has been obtained in other quantum cosmological bounce scenarios, as in Bergeron et al. (2014), and in Refs. Ashtekar et al. (2006); Taveras (2008) when the background fluid has a dust-like equation of state.

A more elaborated and detailed model containing two fluids, dust and radiation, can be found in Ref. Pinto-Neto et al. (2005). The model parameters can be chosen such that the radiation fluid dominates during the bounce, and the dust fluid dominates far from the bounce scale. The Bohmian solution for the scale factor coming from the phase of the wave solution of the corresponding Wheeler-DeWitt equation reads

 a(η)=ae⎡⎣(ηη∗)2+2ηbη∗√1+(ηηb)2⎤⎦, (9)

where is the scale factor at matter-radiation equality, and parameters and are related to the wave-function parameters [similarly to the case of a single fluid where the spread of the initial Gaussian distribution ends up being the bouncing time scale in Eq. (7)].

It is, nonetheless, more convenient to reparametrize the bounce solution using observable related quantities. In this section, all quantities calculated at a time when the scale factor has the same value as today will be denoted by the subscript . Expanding Eq. (9) for large we obtain the Hubble function

 H2≈4aeη2∗(1a3+aea4), (10)

from where we can readily identify the dimensionless density parameters today and as the coefficients of and , respectively

 Ωm0=aea04R2Hη2∗,Ωr0=(aea0)24R2Hη2∗, (11)

where is the critical density today, is the co-moving Hubble radius, and the energy densities of matter and radiation. Next, expanding the Hubble function for large , i.e., considering the dust domination only in the far past, we get

 H2≈H20Ωr0(x4−η2bx6R2H),

where we introduced the redshift-like variable . From the expression above, we see that near the quantum bounce the Hubble function evolves as a classical Hubble function in the presence of a radiation fluid with density parameter , and a stiff matter fluid with negative density parameter given by

 Ωq0=−Ωr0x2b,xb≡RHηb√Ωr0. (12)

Hence, the quantum effect we have calculated, which stops the contraction and realizes the bounce, is dynamically equivalent to a bounce caused by the presence of an additional stiff matter fluid with negative energy, besides the usual matter and radiation fluids, in a classical cosmological scenario obeying the Friedmann equation. Note, however, that this equivalence is valid only at the background level.

Using these new parameters, we obtain

 H2≈H20Ωr0x4[1−(xxb)2]. (13)

Consequently, provides the scale factor where the bounce takes place (apart from a small correction coming from the dust matter density). Finally, we can invert the expression above to obtain all wave-function parameters in terms of the new observable related ones,

 ae=a0Ωr0Ωm0,η∗=2RH√Ωr0Ωm0,ηb=RHxb√Ωr0. (14)

The curvature scale at the bounce can be calculated as

 Lb =1√R∣∣∣η=0=√a3(η)6a′′(η)∣∣ ∣∣η=0 =abηb√6(2γb+1)=1√(2γb+1)a0RHx2b√6Ωr0,

where is the four dimensional Ricci scalar and

 γb≡Ωm0(4xbΩr0), (15)

is the ratio of the dust and radiation matter density at the bounce (where the factor of 4 was included for later convenience). Imposing that the bounce scale is larger than the Planck scale, , we can obtain an upper bound on . This bound is relevant since we expect that the Wheeler-DeWitt equation should be a valid approximation for any fundamental quantum gravity theory only at scales not so close to the Planck length. Using we obtain

 a0RHLp≈8×1060,

and consequently,

 xb≲√81030(6Ωr0)1/4≈2×1031.

In the above calculation, we assumed because the bounce must also happen at energy scales higher than the start of the nucleosynthesis (around ), which implies (using cosmic microwave background radiation temperature ), and we are assuming that should not be much smaller than its usual value . Hence we get

 1011≪xb<2×1031. (16)

Using the parameters above, we make the definitions

 ¯η=ηηb,¯k=kηb,andrb=mabηb,

which are the natural parameters appearing in the equations we will solve, as we will see in the following sections. With this definition, it is easy to see that

 rb=abηbLC≈LbLC,LC≡1m, (17)

where is the Compton wavelength of the massive particle. Note that usually because the curvature scale at the bounce is much smaller than the Compton wavelength, or the mass of the particle is much smaller than the mass-energy scale at the bounce.

In terms of the new parameters we write the Bohmian trajectory as

 a(η)=ab(¯η2γb+√1+¯η2), (18)

Note that, the dust and radiation terms have equal weight at , which is the same result one would obtain substituting in the equation above. In the case of pure radiation ( and, therefore, ), the scale factor reduces to

 a(η)=ab√1+¯η2, (19)

which is exactly the trajectory given in Eq. (7).

The Klein-Gordon equation for the modes of a free massive scalar field in a Friedmann background reads Birrell and Davies (1982)

where . Using the scale factor from Eq. (18) and the natural parameters, the mode equation reads

 ϕ′′k+(ν2−V)ϕk=0, (21) V≡(1−6ξ)[(1+¯η2)−3/2+2γb]√1+¯η2+¯η2γb, (22) ν2≡¯k2+r2ba2a2b, (23)

where the prime denotes derivatives with respect to , i.e., , represent the potential felt by the modes , and the frequency of the mode . It should be pointed out that, in the presence of the dust fluid, the potential changes twice: as one can see from Eq. (22), the radiation fluid becomes dominant in the denominator at a different time than in the numerator (at , and , respectively, see also Fig. 1). In the Appendix A we also realize that the function also plays the role of a potential for . Using the definitions above we get

 ν′ν=r2b¯η[(1+¯η2)−1/2+2γb](¯η2γb+√1+¯η2)¯k2+r2b(¯η2γb+√1+¯η2)2. (24)

Given a complete set of solutions for the mode equation Eq. (21), a set of creation and annihilation operators is defined, and consequently a vacuum state Wald (1994); Birrell and Davies (1982); Chung et al. (2003). The ambiguity in defining the vacuum state stems from the fact that we do not have a general procedure to define a unique set of modes when space-time does not possess a global time-like Killing vector. One special and suitable choice is the so called adiabatic vacuum Birrell and Davies (1982); Parker and Toms (2009); Fulling (1989). One of the main physical properties of this vacuum state choice is that its vacuum expectation value of the number operator varies minimally when the expansion rate of the universe becomes arbitrarily slow (see also Ref. Chung et al. (2003) for a good review on that). As discussed in Ref. Chung et al. (2003), for a given mode at a time , the adiabatic vacuum can be defined up to a maximum order .111This is a consequence of the fact that the adiabatic expansion is asymptotic. However, in some special cases the series is convergent, and the vacuum can be defined up to an arbitrary function that decreases faster than any finite power of . The maximum order is a monotonically increasing function of . Therefore, large ’s have less ambiguity in their vacuum definition than small ’s.

The adiabatic vacuum state has two points relevant to our problem. First, it may depend on the time chosen to define it. If we impose the adiabatic vacuum condition at a time and evolve the modes through Eq. (21) until , we may obtain a different set of mode functions we would otherwise get by imposing the adiabatic vacuum condition at . The other point about this procedure, which is a consequence of the existence of the maximum adiabatic order discussed above, is that it cannot be applied for all modes , as it depends on the behavior of the mode functions, and for a given time only a subset of modes behave in an adiabatic manner.222Alternatively, one can impose a vacuum state by choosing a boundary condition. This asymptotic state is sometimes called “Bunch-Davies vacuum”. Nonetheless, this choice is not free from ambiguities, and coincides with the adiabatic vacuum up to its approximation order Chung et al. (2003).

The first point can be laid down as follows: given a mode at , we impose that the mode initial conditions for are given by the adiabatic approximation up to the maximum order . Using this as initial condition on Eq. (21), we obtain the solution . Repeating the process at a time and comparing both solutions at , we have

 ϕ(i)k(ηf)−ϕ(f)k(ηf)≲O(k−Nk,[ηi,ηf]), (25)

where is the maximum adiabatic order attainable in the interval .333See for example Eq. (33) of Chung et al. (2003). To measure the difference between vacua, we introduce the norm squared of the Bogoliubov coefficients given by Birrell and Davies (1982); Parker and Toms (2009)

 ∣∣β(i,f)k∣∣2=∣∣ϕ(i)kϕ(f)′k−ϕ(f)kϕ(i)′k∣∣2≡n(i,f)k, (26)

The quantity is the number density of particles with mode measured by observers in the adiabatic vacuum defined at if the initial state was the adiabatic vacuum defined at .

Suppose that the adiabatic approximation is valid through the whole interval . It means that there is some , and consequently

 n(i,f)k≲O(k−Nk,[ηi,ηf]). (27)

In other words, in the ultraviolet limit (UV) , the frequency present in Eq. (21) dominates over all other terms. Also, it can be shown that the maximum adiabatic order increases to infinity in this limit. This is equivalent to say that in the ultraviolet limit there is a strong suppression in the number of particles created: as the adiabatic order goes to infinity, any ambiguity in the vacuum definition must fall faster than any finite power of , as an exponential decay. Hence, there is no divergence in the UV limit, and the particle production is finite (unless some infrared divergence is present).

Inspecting Fig. 1, we note that the potential has a maximum at the bounce, . Thus, any mode with will be in the adiabatic regime during its whole evolution, including through the bounce itself, and hence particle production of such modes will be exponentially suppressed. This was verified numerically, as we will see. In other words, the maximum of the potential provides a natural cutoff: any mode larger than it will be strongly suppressed in the particle creation. On the other hand, for modes smaller than the potential maximum, there is a time interval where the adiabatic approximation fails, and particles will be produced.

An important comment has to be made now: in Ref. Quintin et al. (2014), all calculations are done for modes much less than the Hubble radius at all epochs of their cosmological scenario, which is physically the ultraviolet limit (). Hence, as we discussed above, there is an exponential cutoff for these modes, and particle production is heavily suppressed. Another way to phrase it can be: for modes which never cross the potential, the adiabatic vacuum solution, which matches the boundary condition in the far past before the bounce, is always a good approximation at any time. Hence, it coincides with the solution obtained through the adiabatic boundary condition prescribed in the far future after the bounce. Consequently, one must have for these modes. The particle production which is obtained in Ref. Quintin et al. (2014) comes from the matching conditions they impose, which does not capture the precise quantitative evolution of mode function for these large modes. In other words, they artificially introduce a background discontinuity through the matching approximation. Note that, in Ref. Haro and Elizalde (2015), a discontinuity in the background between the different phases is assumed from the beginning, and this gives rise to particle production which depends explicitly on the assumed discontinuity. In practice, the discontinuity creates an infinite potential, invalidating the adiabatic approximation at that point.

As for the second point, we want to calculate the amount of particle creation using the above adiabatic vacuum prescription and for all modes . Thus, if the adiabatic vacua are defined at and , the first order adiabatic vacuum state, given by the zeroth order WKB solution of equation Eq. (21), coincides with the infinite order adiabatic vacua, and hence they precisely define state solutions for all modes . In practice, we will consider that the scalar field is initially in the adiabatic vacuum state in the far past , compute the evolution of such modes until the expansion era far from the bounce (), and compare it to the scalar field solution at with boundary condition of being the adiabatic vacuum at . That is why we numerically prescribe our initial conditions far from the bounce, in the past and in the future, through first order adiabatic approximated solutions, and we verify that the initial approximations coincide with the numerical solutions obtained with such boundary conditions for a long interval of time before the adiabatic approximation looses its validity.

We are also interested in the energy density of the created particles. Using a vacuum definition at as our system state, the expected value of the energy density at any time with respect to this vacuum is

 ⟨ρ⟩(i)=1a4ηb∫d3k(2π)312[ ∣∣ϕ(i)′k∣∣2+(ν2+H2)∣∣ϕ(i)k∣∣2 −H(ϕ(i)∗kϕ(i)k)′], (28)

where .

As it is well known, the energy density of a scalar field in curved space time is divergent. Not only the usual divergence obtained in Minkowsky but new ones must be taken care of. Nonetheless, in this work we want the study the amount of energy resulting from the particle creation during the bounce phase. For this reason we introduce the expected value of the energy density for the adiabatic vacuum defined at , i.e., . Then, the difference

 Δρ=⟨ρ⟩(i)−⟨ρ⟩(f), (29)

given us a finite quantity (see Birrell and Davies (1982)) representing the amount of energy created by the bounce phase. We can relate the mode functions associated to both states as

 ϕ(f)k(¯η)=α(i,f)kϕ(i)k(¯η)+β(i,f)kϕ(i)∗k(¯η). (30)

The constants and can be readily calculated using Eqs. (74) and (75). In the far future, when the modes are deep in the adiabatic phase, the energy difference is given by

 Δρ=12π2η4ba4∫∞0d¯k¯k2n(i,f)kν, (31)

and the number density of created particles is

 n=12π2η3ba3∫∞0d¯k¯k2n(i,f)k. (32)

In the massive case, this provides the energy density when is large enough and consequently for modes relevant to the integral above, i.e., when modes satisfying are dominant for the integral. As we will see below, the particle number density has an exponential cutoff in the ultraviolet limit, therefore, for the massive case and large enough, we have

 Δρ≈mn, (33)

In the massless case, , where the frequency term is , the energy density yields the usual result for a relativistic fluid,

 Δρ=12π2η4ba4∫∞0d¯k¯k3n(i,f)k. (34)

## Iii Bouncing with a Radiation Fluid

In this section, we consider a Bohmian solution of the Wheeler-DeWitt equation obtained in Peter et al. (2007) for the case of a universe dominated by a radiation perfect fluid only. In this case the scale factor is described by Eq. (19).

The Klein-Gordon equation Eq. (21) simplifies to

 ϕ′′k+⎡⎣¯k2+r2b(1+¯η2)−(1−6ξ)(1+¯η2)2⎤⎦ϕk=0. (35)

Note that the dependent term in (35) goes to zero at both past and future infinity, , for whatever value of . Thus, the vacuum solutions at these asymptotic limits do not depend on . The term [Eq. (24)] simplifies to

 ν′ν=r2b¯η¯k2+r2b(1+¯η2). (36)

Comparing to this term also goes to zero at but at a slower rate, i.e., and .

Figure 1 shows the gravitational potential and the mass term for many different cases. The potential goes to zero in two different ways, depending on the presence of dust matter. Thus, in the massless case, the total modification to the frequency vanishes asymptotically, giving simple plane-waves as asymptotic solutions of (35) for any value of . In the massive case, the mass term dominates for , while the dependent term vanishes in this limit. Hence, the asymptotic solutions are mass dependent, but still do not depend on .

### iii.1 Conformal Coupling

Let us start with the simpler conformally coupled case (), which is well-known in the literature Audretsch and Schäfer (1978); Audretsch and Schafer (1978); Birrell and Davies (1982). The time dependent term of the frequency is just the mass term, exemplified in Fig. 1 and Eq. (35), reduces to

 ϕ′′k+[¯k2+r2b(1+¯η2)]ϕk=0. (37)

This equation has exact solutions in terms of parabolic cylinder functions Gradshteyn and Ryzhik (2014). The normalized solutions that match the adiabatic vacuum solution are Birrell and Davies (1982)

 ϕ(i)k(¯η) =e−π8λ(2rb)1/4Diλ−12[(i−1)√rb¯η], ¯η<0, (38) ϕoutk(¯η) =ϕink(−¯η)∗, ¯η>0, (39)

where

 λ=rb(1+¯k2r2b). (40)

To calculate the Bogoliubov coefficients, we use the identity Gradshteyn and Ryzhik (2014)

 Dν(z)=eiπνDν(−z)+√2πΓ(−ν)eiπ2(ν+1)D−ν−1(−iz), (41)

to show that, for ,

 ϕ(i)k(¯η)=√2πeiπ/4Γ(1−iλ2)e−πλ/4ϕ(f)k(¯η)−ie−πλ/2ϕ(f)∗k(¯η). (42)

It follows that

 n(i,f)k=e−πλ=exp[−πrb(1+¯k2r2b)], (43)

which falls off faster than any inverse finite power of or . As it was explained in Ref. Audretsch and Schäfer (1978), this is the spectrum of a non-relativistic thermal gas of particles with physical momentum , chemical potential, and temperature, respectively given by

 μ=−rb2abηba2ba2,T=12πabηbκba2ba2,

where is the Boltzmann constant.

From Eq. (32), the number density of created particles is

 n=12π2a3η3b∫d¯k¯k2n(i,f)k=(√rb2aηb)3exp(−πrb),=(√Ωr0a0RHLC)3/2x38e−πrb (44)

The quantity is usually very small. For instance, assuming (which gives roughly ) and the Higgs particle, one of the most massive scalar particles in the standard model, one gets . Hence, we can neglect the exponential in Eq. (44) due to the smallness of . In other words, the exponential provides a very large mass cutoff

 mcutoff=(xb1030)24.7×1015GeV.

Then, using Eq. (33), one can write the energy density of created particles for very large scale factors as

 Δρ=mn≈m[x(Ωr)1/42(RHLc)1/2]3,≈(mmH)5/2x3×10−44g/cm3, (45)

where is the Higgs mass (), yielding

 Ωϕ≡Δρρcrit0≈(mmH)5/2x3×10−15. (46)

The particle density does not depend on the bounce depth, at least not for the cases considered, for masses closer or larger to the cutoff the particle number density is further reduced by the exponential factor.

We have also calculated Eq. (43) numerically, using the technique described in Appendix A, in order to test our algorithm. Later we will follow analogous steps to calculate the Bogoliubov coefficients for the minimal coupling case. To do so, we have used the definition of the Bogoliubov given in Eq. (26) [and in terms of AA variables by Eq. (77)]. Examining Eqs. (58) and (59), we see that controls the evolution. If is always smaller than one, then the amplitude of the forced oscillations on will also be smaller than one. In that case the angle will evolve accordingly to the frequency . When this happens, the adiabatic approximation will be valid during the bounce and, as discussed above, the particle production will be heavily suppressed. In Fig. 2 we see two examples of modes, top panels, where never grows larger than one and the adiabatic approximation is valid everywhere. On the other hand, for modes having for some time interval, the evolution of both the angles and the adiabatic invariants are dominated by this term stopping the oscillations and increasing the value of . Hence, when drops bellow one again, the further evolution of will only add a small negligible oscillatory term. This kind of evolution can be seem in the lower panels of Fig. 2. To compute the numerical evolution of these modes, we used the approximated solutions (Eqs. 7982) using the condition [Eq. (84)] with , i.e., . After this point we compute the numerical evolution.

It is worth noting that, from the discussion above, if never grows larger than one, then the particle production must be strongly suppressed for all modes where that happens. It is easy to check that the maximum (in absolute value) of reads [see Eq. (36)]

 ν′ν2∣∣∣¯ηm=±2rb3√3(¯k2+r2b),¯ηm=± ⎷¯k2+r2b2r2b, (47)

This maximum attains its largest value at . If this value is smaller than one, i.e.,

 rb>23√3≈0.385,

then the particle production is severely weakened. Note from Eq. (44) that the exponential cutoff is active when

 rb>1π≈0.318. (48)

Hence, we can see that the particle creation cutoff can be approximately deduced by examining the term alone.

In contrast, if the mass appearing in satisfies

 rb<23√3≈0.385,

then all modes in the interval

 ¯k<√rb√23√3−rb,

will have a period in their evolution where the term becomes important, and particle creation takes place. For example, in Fig. 2, a Higgs mass was used (), so particles with (lower panels) are produced, while particle creation with (upper panels) is suppressed, because the adiabatic approximation is valid everywhere.

For the modes where particle creation occur, away from the cutoff, we can also infer the amplitude. In this case, the term controls most of the behavior of the AA variables. Equation (86) describes the evolution of during this period. In our case, the function is an even function of , hence, if the adiabatic evolution stops at (here the subindex denotes an arbitrary initial time), then it will restart at . In this case, Eq. (86) shows that, the adiabatic invariant will return to its original value at the point where ceases to dominate its evolution, i.e., , independently of . Also, this point is defined as , for any . Hence, looking at Eq. (79), we see that at this point all modes evolve back to the same amplitude , and we conclude that the amplitude of does not depend on for modes away from the cutoff. Note that this is consistent with Eq. (43), i.e., for small , .

Figure 4 compares our numerical solution using Eq. (77) with the exact one given by (43). The comparison shows that our numerical calculation predict correctly the particle production in the conformal coupling case.

### iii.2 Minimal Coupling

Next, we calculate numerically the Bogoliubov coefficients for Eq. (35) in the minimal coupling case, , the generalization for any constant value of being simple enough. As it can be seen from Eq. (35), the solutions do not depend on at a sufficient time distance from the bounce, . However, in this case, we have to distinguish the massive case from the massless case: the massless case is trivial only in the conformal coupling case, where there is no particle production [Eq. (35) reduces to a collection of free harmonic oscillators in this situation].

#### iii.2.1 Massless Case

In the massless minimally coupled case, the Eq. (21) for the modes reduces to

 ϕ′′k+⎡⎣¯k2−1(1+¯η2)2⎤⎦ϕk=0. (49)

It can be seen from the equation above that the minimal coupling term goes to zero at infinity, and the equation reduces to that of a simple harmonic oscillator. The vacuum mode solution is then given by

 ϕk(¯η)=(2¯k)−1/2.exp(−i¯k¯η). (50)

In practice, from (49), the plane-wave vacuum solution can be consistently used for . The same conclusions arise from the inspection of condition (84) and the approximated solutions Eqs. (7982).

Following the discussion for the conformal coupling case, here, the controlling term is simply

 Vν2=1¯k2(1+¯η2)2.

Its maximum with respect to the time variable is simply . Therefore, any will result in particle creation, i.e., there is an interval where , and the adiabatic approximation no longer holds. On the other hand, the modes are adiabatic during the whole evolution, and any particle production is suppressed. Thus, in this case, we expect an exponential cutoff at .

The adiabatic regime ends at approximately

 η0=−√1¯k−1,

where we are considering only modes that produces particles, i.e., . Them, using the solution for the dominated given in Eq. (88), and evolving from to results in

 I=I0[1+(¯k|η0|+tan−1|η0|)2],≈I0(π4¯k2−2π3√¯k+1). (51)

In this case is also equal for all modes. The end of the adiabatic regime is determined by , and the amplitude of during this regime is also determined by this same quantity. Thus, we expect from the computation above that for . In Fig. 3 we can see clearly that the frozen value of and increases when gets smaller, consistently with the discussion above.

We solved the full system using the AA variables and Eqs. 58, (59) and Eqs. (61), (62), (Eq. (49)) numerically for as described in Appendix A, using to obtain the initial integration time. As in the conformal coupling case, for greater than the maximal height of the potential, , the adiabatic approximation is valid everywhere and is negligible. In fact, as we discussed before, in this regime decays exponentially with . For , one also could use the approximation methods developed in Refs. Mukhanov et al. (1992); Celani (2011) to show our result above that . We have verified numerically these assertions. This means physically that only those modes with wave number near the bouncing energy scale contribute effectively to particle creation.

To calculate the Bogoliubov coefficients, we use Eq. (77) together with our numerical results. Figure 5 shows the Bogoliubov coefficients. Integrating numerically in , we obtain the particle number density

 n≈6.7×10−2a3η3b≈3×10−2x3(xb1030)3[cm−3]. (52)

This result is the same for the massless and massive cases, within the numerical precision. For massless particles, the energy density calculated through Eq. (34) gives

 Δρ≈3×10−2a4η4b,Ωϕ=2.5×10−10x4(xb1030)4. (53)

Hence, only for one should have a significant amount of massless scalar particles, which could exceed the radiation density in the universe. Note, however, that such would imply a , which goes beyond the scope of the calculations done in this paper.

#### iii.2.2 Massive Case

In the limit , the coupling term depending on is negligible because it falls off with while the mass related term goes with . Thus, in this region, we can use the exact solutions (38) and (39), being understood that they are only approximately valid away from the bounce. In the numerical calculation we used both the first order adiabatic approximation and the exact solution to provide the initial conditions for the numerical system (again, as described at Appendix A). The results are the same, within the required precision. In order to test the consistency of the choice of initial time, we evaluated the initial exact solutions along with our numerical solutions. We did this for a number of values of the parameters, with and , and found no inconsistency: the numerical evolution was indistinguishable from the exact solution (38) for sufficient early times.

Figure 5 shows the Bogoliubov coefficients in the minimal coupling case with many different masses. It can be seen that the Bogoliubov coefficients are many orders of magnitude greater than in the conformal case, see Fig. 4. In Fig. 6, we can see examples of the power-laws present in the Bogoliubov coefficients. For small , the spectrum has its shape similar to the massive conformal coupling case, but the amplitude is increased by a factor of (this factor was found by fitting the small interval with a power-law). This can be understood looking at the lower panels in Fig. 6. In these plots, it is evident that the term dominates earlier, and controls the way the field leaves the adiabatic regime. Nevertheless, around the bounce phase, and differently from the conformal coupling case, the term dominates and increases substantially the amplitude of the adiabatic invariants. For larger , the maximum of is much smaller than that of [see Eqs. (47)]. Therefore, in these cases it is the the one which halts the adiabatic regime and, consequently, the spectrum ends up similar to the massless case.

Finally, as discussed in Sec. III.2.1, the total particle number density for the massive and massless are equal, within the required precision, and they are both given by Eq. (52). In the ultra-relativistic limit, the energy density is given by Eq. (53). In the non-relativistic limit, the energy density is described by a dust like fluid [Eq. (33)], namely

 Ωϕ=mnρcrit0=x3(mmH)(xb1030)32.1×106. (54)

Hence, for the Higgs particle, any bouncing model of this type with will produce a dust-like fluid of Higgs particles with energy density larger than the critical density today.

## Iv Bouncing with Two Fluids: Radiation and Dust

In this section, we consider a model of the universe where its energy content is dominated by two fluids, dust and radiation, such that the radiation fluid dominates during the bounce and the dust fluid dominates in the far past before the bounce, and in the far future after the bounce. This is a more realistic bouncing model, not only because it takes into account the observed dark matter distribution in the universe, but also because the matter domination epoch can account for the almost scale invariant spectrum of cosmological perturbations indicated by recent observations Planck Collaboration et al. (2014, 2015).

The solution for the scale factor is given by Eq. (18). Using the interval for defined at Eq. (16), we get the following interval for ,

 3.7×10−29<γb≪7.5×10−9, (55)

where we used the values for and discussed in Sec. II. To understand the effect of adding a dust-like fluid to the background model, we first study the potential term , which dictates the mode cutoff on . This controlling term is a monotonically decreasing function of [see Eq. (22)], and has its maximum at ,

 Vν2∣∣∣¯η=0=1+2γb¯k2+r2b. (56)

Within the allowed interval (55), we see that the presence of dust does not change significantly the maximum, which will be approximately the same for this whole interval. Therefore, the cutoff does not change significantly in the presence of dust.

Away from the cutoff, the two controlling terms are and , which have different intervals of importance for different values of , and . For this reason, have different forms as a function of , as appearing in Figs. 7 and 8. For instance, in the massless case, the infrared increase of begins for modes when , implying that the potential term be larger than one, and the adiabatic regime terminates during this interval. We can see this behavior in Fig. 7, where changes its shape around , and in Fig. 8 where it changes at .

Thus, apart from the infrared divergence in the massless case, the amount of particle creation not only does not depend on the mass of the particles (except for the conformal case), but also does not depend on the radiation-matter equality constant , for values within the constraint (55). The Figs. 5, 7 and