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

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.

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).

In the models we will analyze, the Belinsky-Khalatnikov-Lifshitz (BKL) instability Belinskii et al. (1970) is not addressed. Of course it would be better to work out models without such instability. However, to insert an extra ekpyrotic scalar field Khoury et al. (2001); Buchbinder et al. (2007); Lehners et al. (2007) with suitable but ad-hoc potential seems excessive to solve this problem, turning the advantage of this type of models quite subjective. Note that any cosmological model, either inflationary, bouncing, or any other approach to primordial cosmology, has a much more serious problem to deal with: the large degree of initial homogeneity necessary to turn all these models compatible with observations (or at least to turn them computationally feasible). This is infinitely more serious than the BKL problem, as long as one has to turn identical the infinite many possible functions of time per space point characterizing a general inhomogeneous gravitational field in order to obtain a homogeneous geometry. If some unknown physical mechanism or theory of initial conditions is capable to justify such extreme fine tuned state, then it would not be a big surprise that it could also make identical the three remaining time functions characterizing the three directions of space. Once one has initially assumed a homogeneous and isotropic universe, one can show for the models we are considering that the shear perturbation will never overcome the background degrees of freedom, even growing as fast as in the contracting phase. This is because the shear perturbation in such models is multiplied by a very small number, a factor , where is the Planck length and is the Hubble radius today, and if the bounce is not very deep, the shear will always remain sufficiently small (see Ref. Vitenti and Pinto-Neto (2012); Vitenti et al. (2014) for details and Battefeld and Peter (2015) for a discussion in a broader context). Note that, besides the Ekpyrotic solution, there are also other possibilities to isotropize the universe, e.g., by adding non-linear effects in the matter content at high energy scales Bozza and Bruni (2009) or through quantum effects Pinto-Neto et al. (2000). However, we do not think it is necessary at this stage to add such effects to our models. Without them, and in a model already containing a radiation fluid from the beginning, there is no need to have some sort of reheating after the bounce in our models. Hence the present work differentiates substantially from previous work on this subject Quintin et al. (2014); Haro and Elizalde (2015); Hipolito-Ricaldi et al. (2016), where the bounce mechanism is completely different, and an ekpyrotic phase is present. It is worth mentioning that such approach, where one assumes an initial contracting phase already homogeneous and isotropic containing no exotic components, were previously study in many papers Peter et al. (2003); Pinto-Neto et al. (2005); Peter et al. (2005, 2006, 2007); Peter and Pinto-Neto (2008); Bessada et al. (2012).

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


where we have defined

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


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


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 :


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


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:


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


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


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


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

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


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


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,


The curvature scale at the bounce can be calculated as

where is the four dimensional Ricci scalar and


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

and consequently,

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

Figure 1: The gravitational potential and the mass term for the parameters discussed in Sec. II. In the initial phase, the potential grows as a power law. If there is dust, then the power-law changes during the dust radiation transition. The potential attains its maximum near the bounce. In the minimally coupled case (), the potential maximum is . These features of the potential are shown in the continuum and dashed lines. In the massive case, the mass term dominates at early times, larger the mass longer the time interval it dominates the mode evolution. This is shown by the dotted and dot-dashed lines in the figure. Note that the mass term dominates the gravitational potential up to the radiation dominated epoch, unless the particle mass is very small (). Hence the presence of dust does not affect much the particle production. For the same reason, the solutions at past and future infinity do not depend on , unless it becomes very close to the conformal coupling.

Using the parameters above, we make the definitions

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


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


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


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


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


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


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)


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


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


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


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


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


and the number density of created particles is


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


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

Figure 2: Evolution of the adiabatic invariants for the conformal coupling case using the Higgs mass . The four figures show the evolution for different values of . The first one, in the top left panel, shows a mode high enough such that the angle evolution never gets dominated by the term. The oscillations, in this case, never cease. The gap appearing in this panel is just the result of using a logarithm scale in , otherwise, we would get regular oscillations in . The top right panel shows the evolution for a larger . In this case, it becomes more evident the two oscillatory phases appearing in Eqs. (7980) which also show that the amplitude of their evolutions is . Consequently, any change in the evolution of this factor can be seen in and . The lower left panel shows unequivocally that when the term becomes larger than one, it dominates the evolution of the adiabatic invariants, and stops their oscillatory behaviors. In this case, the adiabatic invariants no longer oscillate after the bounce and instead freeze out to constant values. We can also see that during the dominated time interval, the adiabatic invariants change very little, growing in the contracting phase and shrinking back to one in the expanding phase, when the potential gets smaller than one. Finally, in the lower right panel, we can see the same behavior for a smaller value of .

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


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


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


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)




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


to show that, for ,


It follows that


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

where is the Boltzmann constant.

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


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

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


where is the Higgs mass (), yielding


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.

Figure 3: Evolution of the adiabatic invariants for the minimal coupling and massless case without dust. The four figures show the evolution for different values of . The description follows closely the discussion of Fig. 2, the main difference being that in this case the cutoff happens at . Hence, in the sequence we see the adiabatic regime valid during all times for and the amplitude amplification of the adiabatic invariants as the value of decreases.

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)]


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

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


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

then all modes in the interval

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


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


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

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

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


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.

Figure 4: The Bogoliubov coefficients in the conformal coupling case with the Higgs mass (), mass () and mass (). We compare our numerical solutions with the exact solution given by Eq. (43).

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


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


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.

Figure 5: The Bogoliubov coefficients in the minimal coupling case with the Higgs mass (), mass (), mass () and the massless case. In this figure we also plotted the first order approximation of the amplitude obtained in Eq. (51), showing an excellent agreement with the full numerical calculation.

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 6: Evolution of the adiabatic invariants for the minimal coupling, massive case (Higgs mass), without dust, for different values of . In this situation one must consider the effect of together with . In the top panels, the values of are such that is less important than , hence, the behavior is similar to the massless case showed in Fig. 3. The lower panels, show the situation for small . In this case, the contribution of the mass term is important and overcome the potential term before and after the bounce, yielding the conformal coupling spectrum for these modes. The dominance of near the bounce, absent in the conformal coupling case, makes the amplitude increases substantially with respect to the former case.

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


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 ,


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 ,


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 .

Figure 7: The Bogoliubov coefficients in the minimal coupling case with the Higgs mass (), mass (), mass () and the massless case. In contrast with Fig. 5 here the background includes a dust like fluid (). The result is similar to the radiation only, with the notable exception of the massless case. Note that, for very large wavelengths there is a growing slope that creates a infrared divergence in the total particle number.

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