R2-Cosmology and New Windows for Superheavy Dark Matter

Evolution and heating of the universe in R2-modified gravity are considered. It is shown that the universe’s history can be separated into four different epochs: (1) inflation, (2) heating due to curvature oscillations (scalaron decay), (3) transition to matter dominated period, and (4) conventional cosmology governed by General Relativity. Cosmological density of dark matter (DM) particles for different decay channels of the scalaron is calculated. The bounds on the masses of DM particles are derived for the following dominant decay modes: to minimally coupled scalars, to massive fermions, and to gauge bosons.


Introduction
The mystery of dark matter (DM) is one of the central problem of modern cosmology. There is a large set of independent pieces of observations which unambiguously point to the existence of a new invisible form of matter disclosing itself through its gravitational action. For a review and list of references see [1][2][3].
According to multiple observations, summarized in Ref. [4], the fractional mass density of dark matter is: Here ρ crit is the critical energy density of the universe: where M Pl = 1.22 × 10 19 GeV = 2.18 × 10 −5 g is the Planck mass related to the Newton gravitational constant as G N = 1/M 2 Pl (We use the natural system of units with c = k = h = 1) and H 0 is the present day value of the Hubble parameter, for which we took H 0 = 100 h km s −1 Mpc −1 ≈ 70 km s −1 Mpc −1 . ( Thus, the observed mass density of dark matter in contemporary universe is: The existence of DM and the magnitude of its contribution into the total mass density of the universe follow from the analysis of the data on: • flat rotational curves around galaxies; • equilibrium of hot gas in rich galactic clusters; • the spectrum of the angular fluctuations of Cosmic Microwave Background (CMB) Radiation; • onset of Large Scale Structure (LSS) formation at the redshift z LSS = 10 4 prior to hydrogen recombination at z rec = 1100. Presently, possible carriers of dark matter are supposed to belong to two distinct groups: microscopically small (elementary particles) and macroscopically large. The latter may include primordial black holes (PBH) with masses starting from 10 20 g up to tens solar masses, topological or non-topological solitons, and possible macroscopic objects consisting e.g., from the so-called mirror matter. A popular form of elementary particle dark matter is called WIMPs (Weakly Interacting Massive Particles), while quite a few other particles types are also discussed in the literature. Macroscopic dark matter of stellar size is abbreviated as MACHO(Massive Astrophysical Compact Halo Object).
The Zoo of the dark matter particles includes, in particular, axions with masses 10 −5 eV or even smaller, heavy neutral leptons with the mass of a few GeV, particles of mirror matter, the so-called superheavy Wimpzillas, and many others. Until recently, one of the most popular candidates was the Lightest Supersymmetric Particle (LSP) which, according to the low energy minimal supersymmetric (SUSY) model, should have a mass of several hundred GeV. With the interaction strength typical for this model the cosmological mass density of these particles was predicted to be close to the observed density of DM. However, the absence of signal from supersymmetric partners at Large Hadron Collider (LHC), if not excluded, but considerably restrict the parameter space of the minimal supersymmetric model.
In the present review, we consider cosmological evolution in the theory where the linear in curvature term is supplemented by the quadratic one. In other words, the Lagrangian density of General Relativity (GR) is modified as follows: where R is a curvature scalar and M R is a constant parameter with dimension of mass. We show that in such a theory, cosmological evolution can be separated into four distinct epochs.

1.
Quasi-exponential (inflationary) expansion, which was realised if the initial value of R was sufficiently large. During this period, R slowly decreased to zero. 2.
Curvature (scalaron) dominated regime, when R oscillated around zero with the amplitude inversely proportional to time. At this stage oscillations of R led to particle production and to universe heating. 3.
Transition period when R dominance was changing to the dominance of (relativistic) matter.

4.
Onset of the conventional cosmology governed by General Relativity.
In the course of this review, we assume that the cosmological background is described by the spatially flat Friedmann-Lemaître-Robertson-Walker (FLRW) metric: where a(t) is the cosmological scale factor and H =ȧ/a is the Hubble parameter at an arbitrary time moment.

Freezing of LSP in Conventional Cosmology
The number density of stable dark matter particles in contemporary universe is determined by the relation of their annihilation rate to the rate of the cosmological expansion. At a certain moment, the annihilation efficiently stops and the number density of the particles tends to a constant value in the comoving volume. This process is called the freezing of species and was first studied in Refs. [5,6].
Calculations of the frozen number density of stable relic particles, X, are based on the kinetic equation, derived by Zeldovich in 1965 [5]: where n X is a number density of particles X, σ ann is XX-annihilation cross-section into light particles, v is the center-of-mass velocity. The angle brackets mean thermal averaging over the medium. We have made the standard assumption of absence of the charge asymmetry between X andX-particles, i.e., assumed equality of their number densities: n X = nX.
For annihilation of non-relativistic particles in S-wave the thermal averaging is trivially reduced to the product σ ann v in the limit of the vanishing velocity of colliding X-particles. As a guiding example for the cross-section estimate we take the well known cross-section of non-relativisitc e + e − -annihilation into photons (see e.g., [7]): where α em = 1/137 is the fine structure constant at low energies and v is the electron (or positron) center-of-mass velocity. Using this result we estimate the annihilation crosssection of heavy X-particles into all light relativistic species as where M X is the X-particle mass, α is a coupling constant, in supersymmetric theories α ∼ 10 −2 , and β ann is a numerical parameter proportional to the number of open annihilation channels, which can be of order of ten or even larger. The contribution to β ann from bosons and fermions may be somewhat different but we ignore this complication.
In the case that X-particles are the Majorana fermions, hence they should be in antisymmetric state, the annihilation proceeds in P-wave. So we take for the annihilation cross-section the expression where the last factor came from thermal averaging of the squared velocity of X-particles, equal to v 2 = 3T/M X . The equilibrium number density of X-particles, n eq , in non-relativistic limit is given by the expression n eq = g s M X T 2π where g s is the number of the spin states of X-particles. The Zeldovich Equation (7) was first applied to calculations of cosmological number density of heavy neutral leptons, which could be viable dark matter particles, in Refs. [8,9]. After publication [8] Equation (7) took the name "Lee-Weinberg equation", but the proper name, in our opinion, is the Zeldovich equation. Equation (7) can be solved analytically with quite high accuracy. The standard procedure of the solution is the following. At high temperatures, the number density of X-particles should be close to the equilibrium one, since the annihilation rate Γ ann = σ ann v n eq is typically much larger than the Hubble expansion rate. According to the Friedmann equation, the Hubble parameter in a spatially flat universe is related to the energy density of the cosmological matter as ρ = 3H 2 M 2 Pl /(8π). If the cosmic plasma is in the state of thermal equilibrium, then ρ = π 2 T 4 g * /30. Hence in the conventional Friedmann cosmology the Hubble parameter is given by the expression: where T is the temperature of the equilibrium cosmological plasma and g * is the total effective number of particle species, to which bosons and fermions contribute respectively 1 and 7/8 for each spin state. Since the number of particle species decreases with dropping temperature g * is a function of temperature and drops down from 100 down to 10 at T ∼ 1 MeV. Initially, the reaction rate is supposed to be much higher than the Hubble rate and at this stage the particle number density is close to the equilibrium one, n X = n eq (1 + δ) with δ 1. So Equation (7) can be approximately linearized in first order in δ. It is convenient to consider the dimensionless number density, f (t), in the comoving volume: where n in is the value of X-particle density at a = a in and T = T in M X , where a is the cosmological scale factor defined below Equation (6), so the X-particles can be considered to be relativistic and thus because the number density of massless particles in thermal equilibrium is 0.12g s T 3 for bosons and 0.09T 3 for fermions, we have taken the approximate value for the numerical coefficient equal 0.1 valid for both cases. As we see in what follows, the final result does not depend upon initial conditions. During the radiation dominated stage the Hubble parameter and the scale factor behave as: Using Equation (12) we find the connection between temperature and time in cosmology based on General Relativity: Correspondingly,Ṫ In terms of new function f Equation (13) and dimensionless variable x = M X /T Equation (7) takes the form: where and f eq is defined as where n eq and n in are given by Equations (11) and (14) respectively. The product a 3 T 3 in Equation (18) is approximately constant up to the corrections induced by the heating of the plasma by the massive particle annihilation when the temperature drops below these particle masses. That is why the number density of X-particles calculated in GR cosmology are normalised to the entropy density which is conserved in the comoving volume. Up to this effect the temperature drops down as the inverse scale factor, a/a in = T in /T. According to Equation (14) the coefficient (a 3 in n in )/(a 3 T 3 ) in the GR regime is equal to 0.1g s and Equation (18) turns into: where f eq is given by Equation (20) with n eq and n in defined in Equations (11) and (14) respectively, so: Due to a very large value of Q, the solution, f (x), is initially close to equilibrium: with δ 1. In this quasi-equilibrium regime Equation (21) can be algebraically solved as because d f eq /dx ≈ − f eq for large x, see Equation (22). This solution is approximately valid till δ remains smaller than unity. The condition δ(x f r ) = 1 defines the so-called freezing temperature T f r = M X /x f r , where x f r is equal to: Dimensionless number density f at the moment of freezing is given by the expression: When x becomes larger than x f r , one can neglect f 2 eq in r.h.s. of Equation (21) and this equation takes the form We need to integrate Equation (27) with the initial condition f in = f f r at x in = x f r , where f f r is defined by Equation (26).
For S-wave annihilation Q does not depend upon x and using Equations (19) and (9) we find: and Equation (27) is integrated as: The last limit transition is true if (0.1 g s Q f f r )/x f r 1. Remind that the asymptotically constant value, f as = const, is the dimensionless concentration of X-particles in the present day universe, defined according to Equation (13). It means that asymptotically the number density of X-particles drops down as 1/a 3 , or, in other words, it is conserved in the comoving volume. Let us denote as T as the temperature below which the asymptotically constant value f as is reached.
For P-wave annihilation we have to substitute into Equations (24) and (27) the factor Q P = 3Q/x instead of Q and take into account that Equations (25) and (26) are changed as: In this case, Equation (27) is transformed into: with initial conditions (30) and (31). The solution has the form We apply the obtained above results for calculations of the present day number and energy densities of relic X-particles. To this end we will use the law of entropy conservation in the comoving volume, see e.g., [10]: where T, ρ, and P are respectively the temperature, energy density and pressure of the cosmological matter, which is predominantly relativistic. We normalised entropy at the value of the scale factor, a as , at the temperature T as , when f reaches its asymptotic value f as . For relativistic matter where g * ∼ 100 at T ∼ M X . Since the asymptotic number density of X-particles conserves in the comoving volume, its ratio to the entropy stays constant: According to Equations (13) and (14): as a as a 3 f as .
Entropy density is equal to: Hence the ratio (36) is: where g * (T as ) ≈ 100 is the number of effective degrees of freedom at T = T as . Due to entropy conservation the quantity g * (T) T 3 remains constant in comoving volume for the primeval plasma in thermal equilibrium. So we can rely on this conservation down to temperatures T 1 ≈ 1 MeV, when g * (T 1 ) = 10.75. Below this temperature neutrinos decouple from the primeval plasma, while photons and e + e − -pairs remain in equilibrium with the effective number of degrees of freedom g em * (T 1 ) = 5.5. After e + e − -annihilation the effective number of degrees of freedom in electromagnetic component of the plasma becomes g em * (T m e ) = 2, where m e = 0.511 MeV is the electron mass. The issues of neutrino decoupling and the effective number of degrees of freedom are described in detail in Ref. [11].
The number density of X-particles at T = T 1 , as follows from the conservation law (36), is: where according to Equation (37) n as = n X (a as ) = 0.1 g s T 3 as f as (41) and s as is given by Equation (38) at T = T as . Below T = T 1 entropy is conserved only in the electromagnetic part of the plasma, so where s em is the entropy of electromagnetic component of the plasma consisting of photons and e + e − -pairs and sub-zero indicates that the corresponding quantities are taken in the present day universe. The entropy of photons in the contemporary universe is expressed through the photon number density, n γ0 = 0.24T 3 0 , as where n γ0 = 412 cm −3 , see e.g., astrophysics summery table in [4]. So the present day number density of X-particles is: n as f as s as s em0 = g * (T 1 ) g em * (T 1 ) 0.1g s f as g * (T as ) g em * n γ0 0.24 Let us apply the obtained above results for determination of cosmological energy density, ρ X0 = M X n X0 , of LSP in the low energy supersymmetric model. Substituting in (28) α = 10 −2 , β ann = 10, and g * (T in ) = 100, we achieve by an order of magnitude the estimate: where M TeV = M X /TeV. With the chosen values of parameters the cosmological energy density of X-particles at the present time would be: which is quite close to the observed value of the dark matter energy density (4). Typical lower bounds on masses of supersymmetric partners, according to LHC data, are in TeV region, see e.g., [4,12]. LSPs with the masses of order of a few TeV or higher would overclose the universe and so they should be excluded from the pool of possible dark matter candidates, provided the conventional Friedmann cosmology is valid. In what follows we consider modified cosmological scenario which allows lifting this constraint and to reopen the chance for superheavy particles with interaction strength typical for SUSY to be DM carriers.

Action Modification
The canonical action of General Relativity, the Einstein-Hilbert action [13][14][15], has the form: where R is the curvature scalar and g is the determinant of the metric tensor g µν with the signature convention (+, −, −, −). The Riemann tensor describing the curvature of space-time is determined according to R α µβν = ∂ β Γ α µν + · · · , R µν = R α µαν , and R = g µν R µν . This action is supplemented by the action of matter, S m , which provides the source of gravitational field, the energy-momentum tensor of matter, T µν : Taking variation of the total action, S EH + S m , over δg µν we arrive to the classical Einstein equations: In Refs. [16][17][18] corrections to vacuum Einstein action have been calculated by performing quantum averaging of T µν in external curved space-time. In particular the addition to the action contains the following terms: where the coefficient C = 1/(144πM 2 Pl ), while B depends upon the renormalisation and is unknown. According to Ref. [18] sufficiently large coefficient B leads to exponential expansion of the universe (inflation) at the early stage. The corresponding effective action is taken as: where M R is a parameter with the dimension of mass and, according to assumption, M R M Pl . The modified Einstein equations for theory (51) are the following: where Taking the trace of Equation (52) yields The General Relativity limit can be recovered when M R → ∞. In this case, we expect to obtain the usual algebraic relation between the curvature scalar and the trace of the energy-momentum tensor of matter:

Modified Evolution Equations
Let us start from an empty universe, where T µ µ = 0 and the curvature R does not depend on space coordinates, R = R(t). In FLRW-metric (6), Equation (53) takes the form: This is exactly the same equation as the Klein-Gordon equation for a homogeneous massive scalar field in curved space-time background. That is why R has the name "scalaron" with M R being the scalaron mass.
In metric (6) the curvature scalar is expressed through the Hubble parameter as: Equations (55) and (56) completely define the cosmological evolution in the absence of matter.
We assume that the cosmological matter distribution is homogeneous and isotropic with the equation of state where ρ is the energy density, P is the pressure of matter, and w is a constant parameter during a certain sufficiently long-lasting cosmological epoch e.g., for non-relativistic matter w = 0, for relativistic matter w = 1/3, and for the vacuum-like state w = −1. Correspondingly, the energy-momentum tensor of matter T µ ν has the following diagonal form: The energy-momentum tensor satisfies the covariant conservation condition D µ T µ ν = 0, which in FLRW-metric (6) has the form: In the presence of matter Equation (55) turns into: It is convenient to introduce dimensionless time variable and dimensionless functions: Equations (56), (59) and (60) now become: where prime means derivative over τ and µ = M R /M Pl .

Inflationary Stage
Here we show that with sufficiently large initial value of R, the devoid of matter universe would expand quasi exponentially [18] long enough to provide solution of flatness, horizon and homogeneity problems existing in Friedmann cosmology, for the review see e.g., textbook [19].
We perform a numerical and analytical solution of the system of Equations (55) and (56), which in terms of dimensionless quantities (61) take the form: The initial conditions should be chosen in such a way that at least 70 e-foldings during inflation are ensured: where τ in f is the moment when inflation terminated, to be determined in what follows. This can be achieved if the initial value of r is sufficiently large, practically independent on the initial value of h. Following Refs. [18,20] (see also the subsequent work [21]), we can roughly estimate the duration of inflation neglecting higher derivatives in Equations (65), so we arrive to the simplified set of equations: These equations are solved as: where r 0 is the initial value of r at τ = 0. According to Equation (67), the Hubble parameter behaves as h(τ) = ( √ −3r 0 − τ)/6. The duration of inflation is roughly determined by the condition h = 0: After h(t) reaches zero it started to oscillate with the amplitude decreasing as 2/(3t) and the exponential (inflationary) rise of the scale factor turns into a power law one.
The number of e-folding is equal to the area of the triangle below the line h(τ), thus N e ≈ |r 0 /4|. It is in excellent agreement with numerical solutions of Equation (65) depicted in Figure 1. The evolution of the dimensionless curvature scalar, r, during inflation is presented in the left panel of Figure 2. Inflation terminates when both h and r reach zero. After that, curvature starts to oscillate, as is presented in the right panel of Figure 2, efficiently producing elementary particles.

Post-Inflationary Stage Prior to the Universe Heating
In this subsection we consider the evolution of r and h solving numerically Equation (65) for the values of τ exceeding duration of inflation, τ > τ in f (70). We assume that the density of matter is negligibly small and take y = 0. Numerical solution for r(τ) immediately after the end of inflation is presented in the right panel of Figure 2. For larger τ the solutions, r(τ)τ and h(τ)τ, take very simple forms depicted in Figure 3. Both r(τ)τ and h(τ)τ oscillate with constant amplitudes, so that the curvature, τr(τ), oscillates around zero, while the Hubble parameter, τh(τ), oscillates around 2/3 almost touching zero at the minima. Simulated by the numerical solution and following Refs. [22,23] we search for the asymptotic expansion of h and r at τ 1 in the form: Here r j and h j , j = 1, 2, are some constant coefficients to be calculated from Equation (65), while the constant phases θ r, h are determined from the initial conditions and will be adjusted by the best fit of the asymptotic solution to the numerical one.
Substituting expressions (71) and (72) into the r.h.s. of second equation of system (65) we obtain: where h is found by the differentiation of (72): Comparison of the 1/τ and 1/τ 2 leads respectively to the equations: Neglecting the oscillating sine-terms, which vanish on the average, and taking the average of sin 2 (τ + θ) = 1/2, we find: Similarly, we explore the first equation of system (65) for r and to this end we need to find expressions for r , and r from Equation (71): Now, substituting the expressions for r, r , r , and h into Equation (65) we arrive at: The leading terms proportional to τ −1 neatly cancel out. The terms of the order of τ −2 leads to: Here, as usually, we have taken sin 2 (τ + θ) = 1/2. Now using Equations (75) and (77), we find: so finally: These asymptotic solutions perfectly well agree with numerical ones presented in Figure 3.
The law of the scale factor evolution corresponding to the Hubble parameter (83) was obtained in Ref. [18].
The effect of particle production on the evolution of curvature, R(t), is usually described by an addition of ΓṘ-term into the l.h.s. of Equation (55) with a constant Γ: This is exactly the same as the Klein-Gordon equation for unstable scalar particle with the decay width Γ. For Γ M R and negligible H the solution of Equation (85) looks as The value of Γ depends on the decay channel of the scalaron and is calculated in the following section for different types of the decay.
Note in conclusion that this simple description of the decay by ΓṘ-term in Equation (85) is valid only for harmonic oscillations of R(t). For arbitrary temporary evolution of R the equation becomes integro-differential, non-local in time one [22,24].

Decay into Minimally Coupled Scalars
The scalaron decay width into two massless (or very low mass) scalar bosons was calculated in refs. [18,22,25]. Here we follow our paper [22], where more rigorous approach was used based on the paper [24], which allows deriving closed equation for an arbitrary time evolution of the source field (in our case the scalaron, R(t)), while the traditional methods are valid only for the harmonic oscillations of the source.
The action of a real scalar field φ with an arbitrary non-minimal coupling to curvature has the form: Minimal coupling to gravity means that ξ = 0, this is the case which is considered in the present subsection.
In spatially flat FLRW background (6) action (87) with ξ = 0, leads to the following equation of motion:φ In turn field φ enters the equation of motion for R (53) via the trace of its energymomentum tensor: It is convenient to study particle production in terms of the conformally rescaled field and the conformal time defined according to the equations: Using these quantities we can rewrite the equations of motion as: while action (87) turns into: Above and below ∂ η denotes derivative with respect to conformal time.
Our aim is to derive a closed equation for R taking the average value of the χdependent quantum operators in the r.h.s. of Equation (91) over vacuum, in the presence of an external classical gravitational field R. Our arguments essentially repeat those of Ref. [24], where the equation was derived in one-loop approximation.
The free field χ (0) is quantized in the usual way: where x µ = (η, x), k µ = (E k , k), and k µ k µ = 0. The creation-annihilation operators satisfy the commutation relations: Equation (93) can be transformed into the following integro-differential equation convenient for perturbative solution: where the massless Green's function is The particle production effects are assumed to weakly perturb the free solution, so Equation (97) can be solved as Now we calculate the vacuum expectation values of the various terms in the r.h.s. of Equation (91), keeping only the contribution from the terms linear in χ (1) . The terms of zero order which contain only χ (0) and its derivatives have nothing to do with particle production and can only change the parameters of the theory through the renormalization procedure.
The terms relevant to the particle production are calculated on the basis of the equations: Thus we obtain: Inserting these expressions into (91), we obtain a closed integro-differential equation for R, which will be transformed into ordinary differential equation for harmonic oscillations of R neglecting the slow power law decrease of their amplitude at the scale of very fast oscillation time.
By the same reason the scale factor a varies very little during many oscillation times ω −1 = M −1 R . Thus, we expect that dη/η ∼ dt/t and that the dominant part in the integrals in Equations (102)-(104) is given by derivatives of R, since R ∼ ωR + t −1 R ωR. So the dominant contribution of particle production is given by Equation (103), and is reduced to: The equation is naturally non-local in time since the effect of particle production depends upon all the history of the system evolution.
The rigorous determination of the decay width of the scalaron is described in Ref. [22]. Here we present it in a simpler and intuitively clear way. We will look for the solution of Equation (105) in the form: where R amp is the slowly varying amplitude of R-oscillations, θ is a constant phase depending upon initial conditions, and ω and Γ is to be determined from Equation (105). The term 3HṘ is not essential in the calculations presented below and will be neglected.
Assuming that Γ is small, so the terms of order of Γ 2 are neglected and treating the r.h.s. of Equation (105) as perturbation we obtain: The first, logarithmically divergent, term in the integrand leads to mass renormalization and can be included into physical M R , while the second term is finite and can be analytically calculated at large upper limit of integration, ωt 1, according to: by closing the integration contour in the upper complex half-plane. Comparing the l.h.s. and r.h.s. of Equation (107) we can conclude that ω = M R and where sub-indices (φ, 0) indicate that we consider the decay width of scalaron into a pair of scalar bosons φ minimally coupled to curvature, i.e., ξ = 0. Our result coincides with the width of the scalaron decay into two massless minimally coupled scalars calculated in paper [25] and presented also in Ref. [26]. However, the width of the same process obtained in Ref. [27] is twice smaller, seemingly because of a numerical error.

Decay into Minimally Coupled Scalars, another Method
Now we calculate the rate of the scalaron decay into the same channel in a different way dealing with the energy loss of the scalaron into the produced particles. To this end we will use the equation of motion of the decay products and calculate the energy density of particles, φ, created by the oscillating gravitational field of the scalaron per unit time,ρ φ . Then we compare it to the energy density of the canonically normalized scalaron field [23]: In what follows we use for R the dominant part of solution (84), which in physical time has the form: where θ is a constant phase determined by the initial conditions. For this R the energy density of the scalaron is equal to: Please note that this is formally equal to the critical energy density of matter dominated universe.
Due to the energy conservationρ φ = −ρ R . So for the rate of the energy dissipation of the scalaron we find: whereρ φ is calculated along the standard lines of particle production theory in external time-dependent fields [28].
According to the Parker theorem [29,30] massless particles conformally coupled to gravity are not created by the gravitational field induced by the conformally flat FLRW metric. This is fulfilled for massless spin 1/2-fermions, massless gauge boson (up to conformal anomaly [31,32]), but is not always true for scalar bosons, as in the following lines.
The action of a real massive scalar field φ with non-minimal coupling to gravity, ξφ 2 R, has the form: S leading to the equation of motion: It is convenient to study particle production in terms of the conformally rescaled field, and the conformal time according to Equation (90). In terms of these variables, Equation (115) transforms into where R = −6 ∂ 2 η a/a 3 . The particle production by external time-dependent field V(t) was studied in many works dedicated to the universe heating, see e.g., [33,34]. In particular, the production rate was calculated perturbatively (but not only) in textbook [10] for the (inflaton) field with the harmonic dependence on time: where V 0 as well as Ω c may slowly depend on time. The field χ satisfies the equation: Since dt = adη, and Ωdt = Ω c dη, the conformal frequency is Ω c = aΩ = aM R . It was shown in Ref. [33] that the number density of particles created per unit of conformal time is: where the number density is expressed in canonical way through the transformed field Hence in physical timeṄ Now we apply this result for χ satisfying Equation (116) with ξ = 0, m φ = 0, R determined by Equation (84), and the potential: Keeping in mind that each particle χ carries energy M R /2 we find for the energy production rate per unit of physical time: And finally which coincides with the result obtained in the previous subsection.

Decay into Fermions
Here we calculate the probability of fermion-antifermion pair production by scalar field Φ: We start from the calculation of the decay width of the Φ-meson with mass M R into a pair of light fermions, assuming that the interaction Lagrangian has the form: where g is a dimensionless coupling constant and ψ andψ are the quantum field operators of the Dirac fermions. Making standard calculations we need to take matrix element of operators of the above Lagrangian between one-particle initial Φ-state with zero momentum and fermion-antifermion final state, using the plane wave decomposition of the wave-function operators. In particular the positive energy projection of state (125) where a Φ is the annihilation operator of one-particle Φ state, hopefully it is not confused with the cosmological scale factor a. Taking matrix element squared and integrating over the phase space we find for the decay width the known expression Here dτ 2 is the invariant phase space of the fermion-antifermion pair and the factor 1/(2M R ) is related to standard relativistic normalisation of the plane wave function of Φ (see any introductory book on elementary particle physics, e.g., [35]). The chosen normalisation corresponds to a single particle over all space. Now let us turn to the number density of fermions produced by the oscillating field (125). In the case of the decay the initial scalar particle is described by the plane wave exp(−iM R t) with positive energy. The number density of the initial particles is equal to: So the densities of fermions and antifermions produced per unit time can be obtained from Equation (127) by multiplication with M R Φ 2 0 /2: Now we have to go to conformal time and conformally transformed field, ψ = ψ c /a 3/2 , Φ = Φ c /a. The equation of motion for ψ c becomes as is well know the free Dirac equation with mass scaled as m ψ → m ψ a. The relation between the conformally transformed quantities and the physical ones on the l.h.s of Equation (129) has the form: Analogously the factor a 4 arises in the r.h.s. of Equation (129). The time dependence a(t) is determined by the Hubble parameter (83): Integrating this equation we obtain up to unessential constant factor: where a c (t) ∼ t 2/3 and transition to conformal time is achieved with the relation dη = dt/a c (t).
In the considered case the particle production is induced by the oscillating part of the scale factor and we have to substitute in Equation (129) instead of gΦ the product of the fermion mass, m ψ , by the amplitude of the oscillating part of the scale factor: Therefore, we obtain:Ṅ The decay width of scalaron by definition is Γ = −ρ R /ρ R . Due to energy conservatioṅ ρ R = −2ρ ψ . Produced by the oscillating curvature, fermions have the energy M R /2 each, soρ ψ = M RṄψ /2. Consequently, the width of the scalaron decay into a pair of massive fermions is: whereṄ ψ and ρ R are given by Equations (134) and (112) respectively. Finally we obtain: This result coincides with that presented in Refs. [18,26].

Generalities
In this section we consider the impact of particle production on the damping of the curvature oscillations and the effect of the created matter on the cosmological evolution. Consequently, Equations (55) and (60) acquire respectively the damping term, ΓṘ. The equation for the curvature evolution in the presence of matter with the equation of state P = wρ takes the form: Due to particle production the covariant conservation (59) of the energy density of the primeval plasma is no longer valid: For the description of the energy influx we introduced in the right hand side of this equation non-zero source term,S[R], to be identified in what follows, see Equations (143) and (144). Equation (56) for the Hubble parameter remains the same.
For dimensionless functions and time (61) we obtain the following system of equations: In what follows we solve these equations for different decay modes of the scalaron.

Minimally Coupled Scalar Mode
According to the calculations of Sections 4.1 and 4.2, the width of the scalaron decay into two minimally coupled scalars is equal to: Hence the source which is the contribution toρ φ = Γρ R according to Equations (123) and (112) is:S where R ampl = 4M R /t is the amplitude of R(t)-oscillations. In terms of the dimensionless quantities the source term can be written as: Here r 2 means amplitude squared of harmonic oscillations, r 2 ampl , of the dimensionless curvature r(τ). However, r(τ) not always oscillates harmonically. In this case we approximate r 2 as 2(r ) 2 or (−2r r). For harmonic oscillations these expressions averaged over oscillation period coincide with r 2 ampl . The function r 2 slowly changes with time.

Asymptotic Solution at t
1, but Γt ≤ 1. Let us find asymptotic behaviour of the solutions of Equations (139)-(141), assuming that the matter is relativistic (w = 1/3) and considering τ 1, but γτ ≤ 1. In such case we can neglect γ in comparison with h in Equation (140) and the system (139)-(141) takes the following simple form: The first two equations of this system coincide with Equation (65), so we can use asymptotic solutions (83) and (84) for the dimensionless Hubble parameter and the curvature scalar to find the energy density, y(τ), from Equation (147). In the r.h.s. of Equation (147) r 2 can be understood as the square of the amplitude of the harmonic oscillations of r, i.e., r 2 = 16/τ 2 . Equation (147) can be analytically integrated as: where τ 0 τ is some initial value of the dimensionless time. The asymptotic result weakly depends upon τ 0 .
Taking h(τ) from Equation (83) we can partly perform integration over dτ 1 as It is convenient to introduce new integration variables: In terms of these variables we lastly obtain: As we see below, the integral in the exponent is small, so the exponential factor in expression (151) is close to unity and thus the dominant asymptotic term is y(τ) ∼ 1/(240πτ). Higher order oscillating corrections we estimate as follows. To calculate the asymptotic behavior of the integral at large τ we present the oscillating factor as sin(τη 1 + θ) = 1 2i e i(τη 1 +θ) − e −i(τη 1 +θ) .
The integral over η 1 along the real axis from η 2 to 1 can be reduced to two integrals over ζ from 0 to ∞ along η 1 = η 2 ± iζ and η 1 = 1 ± iζ. The signs in front of iζ are chosen so that the corresponding exponent in Equation (153) vanishes at infinity. Finally we obtain: The effective value of ζ in this integral is evidently ∼ 1/τ, thus I is inversely proportional to τ in the leading order. At large τ it is much smaller than unity, so we can expand the exponential function exp(−8I/3), see Equation (151), up to the first order and obtain: where the subindex (1/3) indicates that w = 1/3 and = τ 0 /τ 1. The last integral is proportional to 1/τ 2/3 and is subdominant. We neglect it in what follows.
In Figure 4 the dimensionless energy density 240πy(τ) is presented as a result of numerical calculation of the integral (148) or equivalently (151) (solid), which is the exact solution of the differential Equation (147). It is compared with the analytic asymptotic expansion (155) of the same integral (148) (dotted). The calculations are done for mildly and very large τ. One can see that the agreement is very good. Now let us consider the asymptotic solution at τ 1, γτ ≤ 1 for non-relativistic matter with w = 0.
In this case Equations (62)-(64) take the form Since µ 1 the impact of the r.h.s. in Equation (157) is not essential and we can use the expressions (83) and (84) for h and r from the previous subsection. Moreover, the numerical solutions presented in Figure 3 strongly support this presumption. The only essential difference with the w = 1/3 case arises in Equation (158) governing the evolution of the energy density. So to calculate y(τ) we can use slightly modified obtained above results. We need to solve Equation (158). Correspondingly, there appears a coefficient (−3) in the exponent, instead of (−4), as in Equation (148): Repeating the same calculations as above we obtain instead of Equation (155):

LSP Density for the Scalaron Decay into Minimally Coupled Scalars
In this section, we find the frozen number density of X-particles assuming that the scalaron decays into two minimally coupled scalars with the decay width (142). According to the solution (155) the energy density of relativistic massless scalars, ρ s , created by the oscillating curvature, R(t), drops down with time as: This expression can be compared with the energy density of matter in the standard GR cosmology: Because of the difference between the cosmological evolution in the R 2 -theory and General Relativity, the conditions for thermal equilibrium in the primeval plasma also differs very much. Assuming that the equilibrium with temperature T is established, we estimate the particle reaction rate as where α is the coupling constant of the particle interactions, typically α ∼ 10 −2 , and β scat is the number of scattering channels, β scat ∼ 100, see discussion below Equation (9). Equilibrium is enforced if Γ scat > H ∼ 1/t or α 2 β scat Tt > 1. The energy density of relativistic matter in thermal equilibrium is expressed through the temperature as: where g * is the number of relativistic species in the plasma. We take g * ∼ 100. This number is generally accepted and includes 3 families of coloured quarks and antiquarks (36), gluons (16), intermediate bosons and photons (6), three families of leptons (12), Higgs bosons (4) and possibly supersymmetric partners. All together this is the number about 100. Using Equations (161) and (164), we find the equilibrium condition for the case of scalaron decay into a pair of massless scalars: On the other hand, for the GR-cosmology, it follows from Equations (162) and (164) that the equilibrium is established when: Let us estimate now the so-called heating temperature, i.e., the temperature of the cosmological plasma after complete decay of the scalaron. It can be estimated from the expression for the energy density of matter at the time moment equal to the inverse decay width t d = 1/Γ. For the decay into scalars, Γ is given by Equation (142) and Hence the heating temperature for the dominant decay of the scalaron into scalar particles is equal to: For M R = 3 × 10 13 GeV and g * = 100 the heating temperature is T hs ≈ 10 9 GeV. We see that the heating temperature is considerably lower than the temperature at which thermal equilibrium is established, see Equation (165).
Equating the energy density (161) to the energy density of relativistic plasma with temperature T given by Equation (164) we obtain: Comparing Equations (170) and (171) with analogous Equations (16) and (17), we see that the connection between temperature and time in R 2 -theory very much differs from that in the conventional cosmology governed by GR.
Considering relations (170) and (171) we find that the Zeldovich Equation (7) for the dimensionless function f (13) with the annihilation cross-section (9) takes the form: In R 2 theory T ∼ t −1/4 , a ∼ t 2/3 , and a 3 T 3 ∼ 1/T 5 , so the ratio (a in /a(t)) 3 can be estimated as where x = M X /T. The initial number density of X-particles can be taken as n in = 0.12g s T 3 in = 0.12g s M 3 X (the final asymptotic result does not depend upon the intial value), and so Since the density of X-particles is given by Equation (13), its ratio to the number density of photons, n γ = 0.24T 3 , is equal to: Equation (172) governing the evolution of X-particles in the scalaron dominated regime is transformed to where Q sc = 0.03 g s α 2 β ann Since the coefficient Q sc in front of ( f 2 − f 2 eq ) in Equation (176) is normally huge, then initially the solution is close to the equilibrium one, f = f eq (1 + δ), with δ equal to Please note that d f eq /dx ≈ − f eq for large x. This solution is valid till δ remains small, δ ≤ 1. According to Equations (11), (13), and (173) and according to Equations (178) and (173) we have: where Q sc is given by Equation (177).
The freezing temperature T f r = M X /x f r at which deviation from equilibrium becomes of order of unity in this case is approximately: Since Q sc 1, then x f r is also large, typically x f r ∼ (10 − 100) depending upon the interaction strength.
After x becomes larger than x f r , f 2 eq can be neglected in comparison to f 2 and Equation (176) with the initial condition f = f f r at x = x f r is simply integrated giving the asymptotic result at x → ∞: The last asymptotic limit is valid when x x f r and Q sc f f r /(4x 4 f r ) > 1. It is fulfilled because x f r ∼ ln Q sc is sufficiently large and according to Equations (180) and (179) (183) Thus, f (x) tends to a constant value, f as (182), but, according to Equation (175), the ratio of the number densities of X-particles with respect to photons drops down strongly, as 1/x 5 , in contrast to the analogous ratio in GR regime. This decrease is induced by the rise of the density of relativistic species created by the scalaron decay. This drop continues till Γt ∼ 1, when scalaron field disappears and the cosmology returns to the usual GR one. It happens at the temperature given by Equation (168). This temperature should be compared with the temperature of the establishment of thermal equilibrium T eq ≈ 2 × 10 −3 M R , as follows from Equation (165).
Our results are valid if T eq > M X > T hs . However, if the condition T eq > M X is not fulfilled, the cosmological number density of heavy massive particles with masses larger than T hs would still be suppressed and even stronger than in the case of T eq > M X .
After the complete scalaron decay, the ratio of the number densities of X-particles to that of photons remains essentially constant, as it normally happens in GR cosmology. Therefore, the present day ratio n X /n γ can be estimated as the value of this ratio at T = T hs (168), i.e., at x equal to: where x f r is determined by Equation (181).
To calculate the present day energy density of X-particles one needs to multiply the r.h.s. of Equation (185) by M X and by the present day number density of photons n γ,0 = 412/cm 3 . Taking g * ≈ 100, α ≈ 0.01, β ann ≈ 10, and M R = 3 × 10 13 GeV we obtain the estimate: This is to be compared with the observed energy density of dark matter ρ DM ≈ 1 keV/cm 3 . We see that X-particles must have huge mass, much higher than M R to make reasonable dark matter density. However, if M X > M R , the decay of the scalaron into XX-channel would be strongly suppressed and such LSP with the mass slightly larger than M R could successfully make the cosmological dark matter. We will not further pursue this possibility here but turn in the next section to LSP being a fermion.

Fermion Decay Mode
Let us assume now that the scalaron decays only to fermions. According to the calculations in Section 4.3 the width of the scalaron decay into a pair of massive fermions with mass m f , as given by Equation (136) is equal to: which agrees with Refs. [18,26]. Decay probability is dominated by the heaviest fermion. The largest contribution into the cosmological energy density at scalaron dominated regime is presented by the decay into heaviest fermion species.
According to Equation (134): This expression can be rewritten in terms of the amplitude of the curvature oscillations, Integrating Equation (147) and making transformation to physical quantities, we achieve that the energy density of the produced relativistic fermions is equal to: From (164) and (190) follows that thermal equilibrium in the case of the decay into a pair of fermions is established, when Using Equations (187) and (190), we find the temperature of the universe heating for the case of scalaron decay into fermions is: We assume that the mass of the lightest supersymmetric particle is considerably smaller than the masses of the other decay products, m X << m f , at least as m X ≤ 0.1m f . Then the direct production of X-particles by R(t) can be neglected. In such a case LSP are dominantly produced by the secondary reactions in the plasma, which was created by the scalaron production of heavier particles.
Using expression (192) for the temperature of the universe heating after the scalaron decay, we find According to Equation (191) cosmic plasma thermalised at temperatures below The time-temperature dependence, as follows from Equations (190) and (164), is: Kinetic equation for freezing of fermionic species can be solved in complete analogy with what was done in the previous section. The relative number density, f , is defined by the same relation (13) and the kinetic equation for f has the same form: where and n in = 0.09g s m 3 X is the initial number density of X-particles at the temperatures T ∼ m X . We take as in Section 5.2.2 where The equilibrium relative number density, f eq , is slightly different from the similar quantity for bosons (179) and is equal to: The freezing temperature is defined by and the so-called frozen value of f is equal to Using expression (197) for Q f we find that the frozen number density of X-particles, i.e., taken at T = T f r = m X /x f r is Some additional burning of X-particles takes place during the period, when f > f eq and Equation (196) is simplified to: The solution of this equation with the initial condition f (x f r ) = f f r for the asymptotic value at large x x f r is trivially found: This would be the asymptotic value of the relative number density of the heavy stable relics in the standard approach. However, n X does not drop down as 1/a 3 , but much faster due to the extra heating of plasma by the scalaron decay, which does not create X-particles if their coupling to the scalaron is sufficiently weak, as is assumed above. One should, however, remember that there exists a continuous production of heavier fermions, which as is mentioned in the second paragraph of this Section, is much stronger than the direct production of LSP, i.e., of the X-particles. However, the heavy fermions f are produced with huge energy E f ∼ M R /2. This energy is thermalized and is transformed into the energy of relativistic species producing M R /T relativistic particles per one X particles created by the f decays. Moreover, some heavy fermions could annihilate without creation of X particles, but this effect is rather weak, even with an account of relativistic delay of the decay.
Let us calculate now the ratio of the number densities of X-particles, to the density of the relativistic species. The latter is taken as m 3 X at the initial temperature, T in = m X , but in realistic case it can be higher by a factor of a few. E.g. the number density of photons is 0.24T 3 and the same is true for other relativistic species: leptons, quarks, gluons, and even for the electroweak bosons, so our assumption leads to some overestimate of X density. The precise value depends upon the concrete model.
After freezing, the number density of X-particles remains constant in the comoving volume, i.e., where n f r is frozen number density of X-particles, a f r is the value of the cosmological scale factor at the moment of freezing, and we used the expansion law a ∼ t 2/3 and the relation between time and temperature tT 4 = const.
where n γ ≈ 412/cm 3 and we take α = 0.01, β ann = 10, g * = 100, m f = 10 5 GeV, and m X = 10 4 GeV. For the chosen values of the parameters Q f ≈ 2.7 × 10 4 , and ln Q f ≈ 10. This energy density should be close to the energy density of the cosmological dark matter, ρ DM ≈ 1 keV/cm 3 . It can be easily achieved with m X ∼ 10 6 GeV and m f ∼ 10 7 GeV:

Gauge Bosons Mode
In this section we consider the scalaron decay induced by the conformal anomaly. Production of massless gauge bosons by conformally flat gravitational field was first studied in Refs. [31,32] and applied to the problem of heating in R 2 -inflation in Ref. [36]. The scalaron decay width for this channel is equal to: where β 1 is the first coefficient of the beta-function, N is the rank of the gauge group, and α is the gauge coupling constant. We take β 2 1 = 49, N = 8. The coefficient β 2 1 N includes also the contribution from fermionic loops into charge renormalisation. However, the fermion corrections are relatively unimportant and are not included.
The coupling constant α at very high energies depends upon the theory and is, strictly speaking, unknown. The evolution of α in the minimal standard model (MSM) is presented in Figure 5, left panel, and the same in the minimal standard supersymmetric model (MSSM) with supersymmetry at TeV scale is presented in the right panel. We can conclude that at the scalaron mass scale, Q = 3 × 10 13 GeV, α 3 ≈ 0.025 in MSM, while in MSSM it is α 3 ≈ 0.04. At Q = 10 10 GeV they are α 3 ≈ 0.033 for MSM and α 3 ≈ 0.05 for MSSM. The values of the running coupling constants are known to depend upon the set of contributing particles. In the case of MSM we take into account all known particles, while in MSSM there is some freedom depending on the explicit form of the SUSY model. However, the variation of the couplings related to this uncertainty does not lead to strong variation of our estimates of the allowed range of massess of dark matter particles.
Here we have taken into account the evolution of the coupling constants only, because their values are known at low energies and the values at high energies can be calculated knowing the underlying theory. The variation of the couplings is quite significant. On the other hand, the mass of the scalaron is fixed at high energies by the magnitude of the density perturbations generated at inflation. As for the masses of the decay products they are fixed on their mass shell and the quantum renormalization corrections are relatively weak. They can acquire some temperature corrections, which could lead to a difference of mass values over vacuum and high temperature plasma. However, we do not claim to know the exact values of masses of the possible DM particles with the interaction strength typical for supersymmetry. We have only presented an order of magnitude estimates of their masses.
Since, according to our results presented below, supersymmetry may possibly be realized at energies about 10 12 -10 13 GeV, the running of couplings according to MSM without inclusion of SUSY particles is probably correct below the SUSY scale. Recall that for particles produced at the scalaron decay Q = 3 × 10 13 GeV, while at the universe heating temperature after the complete decay of the scalaron it is near 10 10 GeV.
So numerically the decay width is: (217) The calculations of the energy density of created gauge bosons can be done in complete analogy with calculations of the energy density of massless scalars bosons performed in Section 5.2. Of course we need to take into account the difference of the corresponding decay width (142) and (216). Correspondingly, the energy density of matter created by the scalaron decay into a pair of massless gauge bosons would be: Following our paper [38] we present alternative, more rigorous arguments leading to introduction of canonically normalised field Φ (110). To this end we rewrite the first term of action (51) in the Jordan frame in the high frequency limit in terms of the cosmological scale factor a(t) in the way analogous to the derivation of the Friedmann equations performed in Ref. [39]. For high frequency oscillations and large value of M R t we have found the solutions (83) and (84), which in physical quantities have the form: Curvature scalar is related to the Hubble parameter according to: The last relation is valid in high frequency limit and for the oscillating parts of H and R which presumably give dominant contribution to the energy density.
Keeping this in mind we can rewrite action (51) as: The last equality is obtained through integration by parts.
Varying action (221) over the scalar field H we obtain the equation of motion: compare to Equation (137). This equation has the oscillating solution multiplied by a slow function of time, such as the presented above solution H ∼ sin(M R t + θ)/t. Now we need to introduce canonically normalized scalar field Φ linearly connected with H for which the kinetic term in the Lagrangian is equal to (∂Φ) 2 /2: In high frequency limit it is essentially the same field as defined above in Equation (110). According to the standard theory the energy density of the scalar field Φ is In high frequency limit Equation (220) leads to identification R = −6Ḣ ∼ M R cos(M R t) and Equation (224) can be rewritten in terms of R as where expression (219) has been used. This result coincides with the expression (112) for the total cosmological energy density in spatially flat matter dominated universe. This agreement confirms the validity of our approach.
The presented equations are valid if the energy density of matter remains smaller than the energy density of the scalaron until it decays.
It is instructive to compare the rate of the energy transferred to matter produced in three different cases of the scalaron decay into minimally coupled scalars, fermions, and gauge bosons due to conformal anomaly with the energy density of the scalaron. Comparing Equations (161), (190), and (218) with (225) we find that in all the cases t cr Γ = 5/3, where t cr is the time when the matter energy density, formally taken, is equal to the scalaron energy density. So the used above equations are not unreasonable. Let us note that the energy density of matter in the case of scalaron decay into fermion-antifermion pair is equal to doubled energy density of fermions (190).
The scalaron completely decays at t ∼ 1/Γ (up to log-correction) and the cosmology turns into the usual Friedmann one governed by the equations of General Relativity. Before that moment the universe expansion was dominated by the scalaron.
If the primeval plasma is thermalized, the following relation between the cosmological time and the temperature is valid: where subindex R at α R means that the coupling is taken at the energies equal to the scalaron mass, since the energy influx to the plasma is supplied by the scalaron decay, and g * ≈ 100 is the number of relativistic species. Consequently, Thermal equilibrium is established if the reaction rate is larger than the Hubble expansion rate H = 2/(3t). The reaction rate is determined by the cross-section of twobody reactions between relativistic particles. The typical value of this cross-section at high energies, E m, is [40]: where β ∼ 10 is the number of the open reaction channels and s = (p 1 + p 2 ) 2 = 4E 2 is the total energy of the scattering particles in their center-of-mass frame, where E is the energy of an individual particle. Hence the reaction rate is where angular brackets mean averaging over thermal bath with temperature T, n rel ≈ 0.1g * T 3 (we do not distinguish between bosons and fermions in the expression), v = 1 is the particle velocity in the center-of-mass system. We perform thermal averaging naively taking E = T in all expressions so s → 4T 2 , instead of m 2 we substitute the particle thermal mass in plasma, i.e., m 2 → 4παT 2 /3 [41][42][43]. Since thermal equilibrium is established when the reaction rate Γ exceeds the Hubble parameter H = (2/3t) we obtain the following equilibrium thermal condition: Using Equation (227), we find that equilibrium is established at the temperatures below Here we took α R = 0.025 and α = 0.033, where α R is defined below Equation (226) and we took the value of coupling constant from Figure 5. Constant α determines the scattering of relativistic particles (228), which for the energies of the order of temperature is equal to 0.033.
The time corresponding to this temperature is equal to where C is defined in Equation (227). Hence M R t eq 1, which is sufficiently long time for efficient particle production.
Another essential temperature for our consideration, is the temperature of the universe heating, when scalaron essentially decayed and the expansion regime turned to the conventional GR one. This temperature is determined by the scalaron energy density at the moment t = 1/Γ anom :

Direct X-Particle Production through the Scalaron Decay
There are two possible channels to produce massive stable X-particles: first, directly through the scalaron decay into a pair of XX and another by inverse annihilation of relativistic particles in plasma.
First, let us consider the scalaron decay. The probability of the scalaron decay into a pair of fermions is determined by decay width (187) with the substitution M X instead of m f : The branching ratio of this decay is equal to: The number density of X-particles created by the scalaron decay only, but not by inverse annihilation of relativistic particles in plasma, is governed by the equation: where Γ X is given by Equation (235), n R = ρ R /M R , and ρ R is defined in Equation (225). So Equation (237) turns intoṅ In the scalaron dominated regime non-oscillating part of the Hubble parameter is H = 2/(3t) and the equation is solved as The equations presented above are valid if the inverse decay of the scalaron can be neglected. This approximation is true if the produced particles are quickly thermalized down to the temperatures much smaller than the scalaron mass.
We are interested in the ratio of n X to the number density of relativistic species at the moment of the complete scalaron decay when the temperature dropped down to T h (234) after which the universe came to the conventional Friedmann cosmology and the ratio n X /n rel remained constant to the present time. This ratio is equal to: If dark matter totally consists of X-particles, their energy density in the present day universe should be equal to the observed energy density of dark matter: From this condition it follows that M X ≈ 10 7 GeV. For larger masses ρ (0) X would be unacceptably larger than ρ DM . On the other hand, for such a small, or smaller M X , the probability of X-particle production through the inverse annihilation would be too strong and would again lead to very large energy density of X-particles, see the following section.
A possible way out of this "catch-22" is to find a mechanism to suppress the scalaron decay into a pair of X-particles. In addition, it does exist. If X-particles are the Majorana fermions, then in this case particles and antiparticles are identical and so they must be in antisymmetric state. Thus, the decay of a scalar field, scalaron, into a pair of identical fermions is forbidden, since the scalaron can produce a pair of identical particles in symmetric state only.

Production of X-Particles in Thermal Plasma
Here we turn to the X-production through the inverse annihilation of relativistic particles in the thermal plasma. The number density n X is governed by the Zeldovich Equation (7):ṅ where σ ann v is the thermally averaged annihilation cross-section of X-particles and n eq is their equilibrium number density. The thermally averaged annihilation cross-section of non-relativistic X-particles, which enters Equation (7), for our case can be taken as (see Equation (10)): where the last factor came from thermal averaging of the velocity squared of X-particles, equal to v 2 = T/M X , which appears because the annihilation of Majorana fermions proceeds in P-wave. We take the coupling constant at the energy scale around M X equal to α = 0.033 and the number of the annihilation channels β ann = 10. This expression is only an order of magnitude estimate. The exact form depends upon particle spins, the form of the interaction, and may contain the statistical factor 1/n!, if there participate n identical particles. In what follows we neglect these subtleties. The equilibrium distribution of non-relativistic X-particles has the form: where y = M X /T and g s is the number of spin states of X-particles. The non-relativistic approximation is justified if M X > T eq ≈ M R /100 = 3 × 10 11 GeV, see Equation (231). Equation (242) will be solved with the initial condition n X (t in ) = 0. This condition is essentially different from the solution of this equation in the canonical case, when it is assumed that initially n X = n eq and in the course of the evolution n X becomes much larger than n eq , reaching the so-called frozen density. As we see in what follows, for certain values of the parameters the similar situation can be realized, when n X approaches the equilibrium value and freezes at much larger value. The other limit when n X always remains smaller than n eq is also possible.
For better insight into the problem we first make simple analytic estimates of the solution when n X n eq and after that solve exact Equation (242) numerically. In the limit n X n eq Equation (242) is trivially integrated: where the new integration variable is defined as z = y 1 /y, the subindex "0" means that the solution is valid for n X n eq , y = M X /T and we have used Equation (227) and the expression for C 0 below this equation.
For the initial temperature we take T in = T eq ≈ 7 × 10 −3 M R , according to Equation (231), and T f in = T h = 5.1 × 10 −6 M R (234). Correspondingly y in = 1.4 × 10 2 M X /M R , and y f in = 2 × 10 5 M X /M R and so y f in /y in ≈ 10 3 .
To check validity of this solution we have to compare n X0 (y) to n eq (244): where we have taken g s = 2, β ann = 10 and lastly, according to the line below Equation (231), α = 0.033 and α R = 0.025. The ratio F 2 (y) is depicted in Figures 6 and 7 as function of y for different values of y in . The ratio remains smaller than unity for sufficiently small y < y max = 50 − 150 depending upon y in . If y f in < y max , the assumption n X n eq is justified and the solution (245) is a good approximation to the exact solution. In the opposite case, when y f in > y max , we have to solve Equation (242) numerically.
To solve Equation (242) it is convenient to introduce the new function according to: where a(t) is the cosmological scale factor and a in is its initial value at some time t = t in , when X-particles became non-relativistic. In terms of z, Equation (242) is reduced to:  Next, let us change the variables from t to y = M X /T. Evidentlyẏ = −y(Ṫ/T). Using time-temperature relation (227), we find Keeping in mind that a in a we find finally: dz dy = 6π g s C 0 α 2 β ann µ 3 y 8 in y 6 where µ = M R /M X . With the chosen above values of α R and α, see the discussion after Equations (227) and (243), we find that the value of the coefficient in the r.h.s. of Equation (251) is 6πg s C 0 α 2 β ann ≈ 10 −7 .
Numerical solution of this equation indicates that z(y) tends asymptotically at large y to a constant value z asym . The energy density of X-particles is expressed through z asym as follows. We assume that below T = T h the ratio of number density of X-particles to the number density of relativistic particles remains constant and hence is equal to the ratio n X /n CMB at the preset time, where n CMB = 412/cm 3 is the contemporary number density of photons in cosmic microwave background radiation. The number density of X-particles is expressed through z according to Equation (247). Thus, the asymptotic ratio of the number densities of X to the number density of relativistic particles is We assume that y in ≈ 10 2 /µ, y f in = y h = 2 × 10 5 /µ, according to the discussion after Equation (245), and so y f in /y in ≈ 2 × 10 3 . Hence the energy density of X-particles today would be equal to: where z asym is the asymptotic value of z(y) at large y but still smaller than y h . The value of z asym can be found from the numerical solution of Equation (251). However, the solution demonstrates surprising feature: its derivative changes sign at y ≤ 10, when n X n eq , as is seen from the value of F 2 presented in Figures 6 and 7. Probably this evidently incorrect result for z(y) originated from a very small coefficient in front of the brackets in Equation (251).
The problem can be avoided if we introduce the new function u(y) according to: In terms of u(y), kinetic equation ( The numerical solution of this equation does not show any pathological features and may be trusted, so we express the contemporary energy of dark matter made of stable X-particles through the asymptotic value of u(y) as Remind that y in ≈ 100/µ and presumably µ > 1.
The asymptotic value u asym is found from the numerical solution of Equation (255) and is depicted in Figures 8 and 9 for different values of µ.  The logarithm of the energy density of X-particles (256) with respect to the observed energy density of dark matter as a function of M X is presented in Figure 10. If M X ≈ 5 × 10 12 GeV, X-particles may be viable candidates for the carriers of the cosmological dark matter.

Possible Observations
This section contains a discussion of some possibilities of observation of the products of X-particle slow decay or enhanced annihilation in ultra high energy cosmic rays. More detailed study of the phenomena considered below demands more work.
There are two possibilities to make X-particles visible: first, due to possible high density of XX-systems and, secondly, because of hypothetical instability of X-particles.
According to results of our papers [23,38,44] the mass of dark matter particles, with the interaction strength typical for supersymmetric ones, can be in the range from 10 6 to the value exceeding the scalaron mass, M R = 10 13 GeV. It is tempting to find if and how they could be observed, except for their gravitational effects on galactic and cosmological scales.
The average cosmological energy/mass density of dark matter particles in the universe is approximately 1 keV/cm 3 , while in galaxies it is about 1 GeV/cm 3 . So their number densities should respectively be: where M 6 = M X /(10 6 GeV). The characteristic annihilation time in a galaxy is: where we have taken σ ann v ≈ 10 −2 /M 2 X . The total energy flux from all annihilations in the Galaxy of the size R gal ≈ 10 kpc = 10 22 cm would be L gal = n gal ER gal /τ ann gal ≈ 10 −15 M −3 6 GeV/cm 2 /s = 3 × 10 2 M −3 6 GeV/km 2 /year, with characteristic energy of the order of E ∼ M X . The annihilation would be strongly enhanced in clusters (clumps) of dark matter [45], especially in neutralino stars [46]. Based on the latter reference, for the annihilation cross-section σ ann v = 4 × 10 −42 M −2 6 cm 2 ≈ 10 −31 M −2 6 cm 3 /s, we can conclude that the observation of XX-annihilation from neutralino stars is not unrealistic.
Due to their huge mass, relic X-particles might form gravitationally bound states and then annihilate such as positronium. Instead of fine structure constant α = 1/137 we must use the gravitational coupling constant α G = (M X /M Pl ) 2 . In complete analogy with para-positronium decay the lifetime of such bound state with respect to annihilation would be where M 13 = M X /(10 13 GeV). The flux of ultra-high energy cosmic rays (UHECR) with energies 10 21 -10 22 eV produced by the population of the bound states of XX, say, from the sphere of the radius of R = 1 Gpc would be: where f is the fraction of the gravitationally bound states of XX with respect to the total number of X-particles.
Comparing this result with the data presented in ref. [47] we can conclude that the flux of the UHECR produced in the decay of XX bound states would agree with the data if f ∼ 10 −11 .
Calculation of f is subject to many uncertainties and demand an extensive work. It will be done elsewhere.
Another chance for observation of X-particles appears if they are unstable. Heavy Xparticles would decay though formation of virtual black holes, according to the Zeldovich mechanism [48,49]. If X-particles are composite states of three fundamental constituents, as proton made of three quarks, their life-time with respect to virtual BH simulated decay would be To make the time τ X,BH larger than the universe age t U ≈ 4 × 10 17 s, we need M X < 10 7 GeV. In this case the products of the decays of X-particles with such masses could be observable in the flux of the cosmic rays with energy somewhat below 10 7 GeV.
The life-time may be further diminished if we apply the conjecture of Ref. [50] which leads to a strong suppression of the decay through virtual black holes for spinning or electrically charged X-particles. However, this suppression does not operate for spinless neutral particles. Moreover it would not be efficient enough to sufficiently suppress the decay probability of the superheavy particles of dark matter with masses of the order of 10 13 GeV. The decay rate may be strongly diminished if X-particles consist of more than three fundamental constituents. For example, if X-particles consist of six fundamental constituents, then the decay life-time would be τ X,BH ∼ This life-time is safely above the universe age t U ≈ 4 × 10 17 s.

Conclusions and Discussion
There is general agreement that the conventional Friedmann cosmology is in tension with the existence of stable particles with interaction strength typical for supersymmetry and heavier than several TeV. A possible way to save the life of such particles, we call them here X-particles, may be a modification of the standard cosmological expansion law in such a way that the density of such heavy relics would be significantly reduced. A natural way to realize such a reduction is the now popular Starobinsky inflationary model [18]. If the epoch of the domination of the curvature oscillations (the scalaron domination) lasted after freezing of massive species, their density with respect to the plasma entropy could be noticeably suppressed by the production of the relativistic species created by the scalaron decay.
The concrete range of the allowed mass values of X-particles, if they constitute the cosmological dark matter, depends upon the dominant decay mode of the scalaron. If the scalaron is minimally coupled to scalar particles X S , the decay amplitude does not depend upon the scalar particle mass and leads to too high energy density of X-particles, if M X S M R . An acceptably low density of X S can be achieved if M X S ≥ M R ≈ 3 × 10 13 GeV, since in this case X-particles should be produced in multi-scalaron processes.
If scalaron decays into fermion-antifermion pair, its decay width is proportional to m 2 f and the decay probability is dominated by the heaviest fermion. In this case X-particles could be produced in secondary reactions in plasma. For sufficiently small m f the energy density of the produced fermions would be small enough and respectively the energy density of X-particles produced either through the fermion decays or annihilation could be close to the observed density of dark matter for M X ∼ 10 6 GeV [23]. The scalaron could also decay into a pair of gauge bosons due to existence of conformal anomaly. In this case, the cosmological plasma initially would dominantly consist from the gauge bosons. In the process of thermalisation of this primeval plasma the X-particles of would-be dark matter could be created. There are two possible processes through which X-particles could be produced: direct decay of the scalaron into a pair ofXX and the thermal production of X's in plasma. If scalaron decays into a pair of XX-fermions, the energy density of X-particles would be close to the observed energy density of dark matter for the X-particle mass below 10 7 GeV. However, in this case the thermal production of X's would be too strong. We can resolve this inconsistency if the direct decay of the scalaron into X-particles is suppressed and due to the fact that a larger M X is allowed, so the thermal production would not be dangerous. The direct decay can be very strongly suppressed if X-particles are Majorana fermions, which cannot be created by a scalar field in the lowest order of perturbation theory. This opens the possibility for X-particles to make proper amount of dark matter, if their mass is about 5 × 10 12 GeV [38].
However, conformal anomaly is not necessarily present in supersymmetry inspired theories. There are some versions of SUSY theories where conformal anomaly is absent, for example N = 4 supersymmetry for which beta-function vanishes, see e.g., review ([51], Section 13.2), and references therein.
On the other hand, N = 4 super Yang-Mills theories are believed to be unrealistic because they do not allow introducing chiral fermions, even if the symmetry is broken spontaneously. Though spontaneous symmetry breaking is considered to be the most appealing way to deal with the theories with broken symmetries, it is not obligatory and the symmetry can be broken explicitly. It is possible to break the symmetry "by hand" introducing different masses to particles in the same multiplet. This would allow to construct a phenomenologically acceptable model. Since the symmetry is broken by mass, the theory would remain renormalizable. At higher energies, much larger than the particle masses, it would behave as N = 4 super Yang-Mills theory and at this energy scale the trace anomaly would vanish.
Apart from that there are phenomenologically acceptable N = 1 and N = 2 supersymmetric theories which possess the so-called conformal window, i.e., in this theories with a certain set of the multiplets trace anomaly vanishes. For a review and the list of references see [52].
Thus, a supersymmetric type of dark matter particles seems to be possible if their mass is quite high from 10 6 up to 5 × 10 12 GeV, or even higher than the scalaron mass, M R = 3 × 10 13 GeV. There is not a chance to discover these particles in accelerator experiments in foreseeable future, but they may be observable through cosmic rays from their annihilations in high density clumps of dark matter, or from annihilation in their gravitationally bound two-body states, or through the products of their decays, since they naturally should be unstable.
It is worth mentioning that with the standard inflation is impossible to relax the conventional bounds on the masses of dark matter particles. The reason for the difference between the scalaron and inflaton cosmologies, implicitly mentioned in the review, is that the scalaron decays much slower than the inflaton and due to that there is a constant and slow influx of energy to the plasma from the scalaron, which diluted the fraction of very massive DM particles. As is mentioned in the review the essential processes of the scalaron decay took place during long time with Γt ≤ 1, while for the inflaton decay the thermalization took place almost instantly. The difference between inflaton and scalaron is related to a huge difference between their couplings to matter. The typical value of the coupling constant of the inflaton decay into a pair of fermions through the interaction gΦψψ is g ∼ 10 −3 , while in the case of the scalaron it is inversely proportional to the Planck mass squared. The equivalent of g for the scalaron is m 2 ψ /M 2 Pl . Even for the case of the minimally coupled scalars g 2 ∼ (M R /M Pl ) 2 ≈ 10 −11 .
Let us note that according to Equations (168), (192), and (234) the temperature of the primordial plasma after the complete scalaron decay is higher than the temperature of BBN by several orders of magnitude. So the temperature when the universe turned to the usual FLWR cosmology is much higher than the BBN temperature equal to 1 MeV and thus BBN remained undisturbed. The energy density of stable would-be particles of dark matter is negligibly small at BBN. Let us summarise. In the conventional cosmology masses of the stable supersymmetric relics (candidates for the DM particles) should be typically below 1 TeV. This is in conflict with the LHC bounds on the low energy SUSY. On the other hand, in R 2 -gravity the masses of the stable particles with the interaction strength typical for SUSY could be much higher depending upon the dominant decay mode of the scalaron. As we showed, for the decay into minimally coupled scalars the mass of the stable relics could be as large as the scalaron mass, M R = 3 × 10 13 GeV. If the scalaron predominantly decays into a pair of massive fermions the mass of DM particles could be about 10 6 GeV. If the scalaron predominantly decays into gauge bosons due to trace anomaly, the allowed value of mass of DM particles could be about 5 × 10 12 GeV.