Superheavy dark matter in $R+R^2$ cosmology with conformal anomaly

Cosmological evolution and particle creation in $R^2$-modified gravity is considered for the case of the dominant decay of the scalaron into a pair of gauge bosons due to conformal anomaly. It is shown that in the process of thermalization superheavy dark matter with the coupling strength typical for the GUT SUSY can be created. Such dark matter would have the proper cosmological density if the particle mass is close to $10^{12}$ GeV.


Introduction
The most popular and natural hypothesis that dark matter consists of the lightest supersymmetric particles was eliminated by the LHC lower bounds on the masses of the supersymmetric partners [1]. These bounds not only excluded low energy SUSY with the energy scale around 10 2 − 10 3 GeV but also excluded supersymmetry at arbitrary high scale, if one assumes the standard cosmological model and the stability of the lightest supersymmetric particle (LSP). Indeed, the energy density of LSP, according to the conventional calculations, is where M LSP is the mass of LSP and ρ (obs) DM ≈ 1 keV/cm 3 is the observed value of the cosmological density of dark matter.
The equation (1.1) and our results obtained in Ref. [2] and in the present paper do not demand full supersymmetry and are valid for any massive stable particle with the coupling strength typical to that in supersymmetry. So in what follows we will not use the abbreviation LSP for these particles but instead call them X-particles.
In our recent work [2], we have shown that in modified R + R 2 cosmology the relative density of LSP can be considerably smaller than that predicted in the standard scenario. This opens the window for the lightest supersymmetric particle with the mass about 1000 TeV to be a viable dark matter candidate. The frozen number density of massive relics is calculated in terms of the present day density of photons of the cosmic microwave background radiation, see e.g. [3]. The relative decrease of the LSP density in R 2 -cosmology is related to an efficient particle production by the oscillating curvature scalar after freezing of the LSP production and annihilation, as it is shown in our paper [4]. Consequently, the number density of CMB photons rises and the ratio of the frozen number density of LSP to the number density of photons drops down.
Production of dark matter particles by the scalaron decay in a different aspect was considered also in Ref. [5].
According to our work [4], the cosmological evolution in R + R 2 theory is considerably different from the standard one, based on the classical General Relativity (GR). In the modified gravity, the cosmological evolution can be categorized into the following four epochs. At first, there was the exponential expansion (Starobinsky inflation [6]), when the curvature scalar R(t) (called scalaron) was very large and was slowly decreasing down to zero. The next epoch began when R(t) dropped down to zero and started to oscillate, periodically changing sign. The oscillations of R led to particle production and this epoch can be called Big Bang. Next, there was the transition period from the scalaron domination to the relativistic matter domination. Finally, after scalaron had decayed completely, we arrived to the standard cosmology which is governed by General Relativity.
The frozen density of massive species strongly depends upon the probability of particle production by R(t). In our previous papers [2,4], we considered the decays into minimally coupled massless scalar particles and into massive fermions or conformally coupled scalars. However, as it is argued in Ref. [7], the production of massless gauge bosons due to conformal anomaly may be significant. We avoided this problem assuming a version of supersymmetric model, where conformal anomaly is absent. Here we clear out this restriction and consider freezing of massive species in the theory where the particle production by oscillating curvature predominantly proceeds through anomalous coupling to gauge bosons.
The paper is organized as follows. In Sec. 2 we summarize our results on cosmological evolution in R 2 -gravity and present the known theoretical estimates of the variation of the coupling constants with changing momentum transfer. In Sec. 3 an estimate of the cosmological number density of X-particles created by direct decay of the scalaron is presented. It is shown that the energy density of X-particles would be close to the energy density of DM if the Xparticle mass is rather small, M X ∼ 10 7 GeV. However, as shown in Sec. 4, in this case the X-particle production by thermal processes in plasma, in turn, would be unacceptably strong.
To avoid this crunch we assume that X-particles are Majorana fermions because in this case their direct production by the scalaron is forbidden. According to the calculations in Sec. 4 the cosmological density of X-particles would be equal to the observed density of DM if M X ∼ 10 12 GeV. In Conclusion, the results are discussed and compared to the other cases studied earlier.

Cosmological evolution in R gravity
This section contains a condensed summary of the main results of our works [2,4]. The action of the theory has the form: where m P l = 1.2·10 19 GeV is the Planck mass, S m is a matter action. Here g is the determinant of the metric tensor g µν taken with the signature convention (+, −, −, −). The Riemann tensor describing the curvature of space-time is determined according to R α µβν = ∂ β Γ α µν + · · · , R µν = R α µαν , and R = g µν R µν . We use here the natural system of units = c = k B = 1. As we see in what follows, M R is the mass of the scalaron field. The spectrum of the temperature fluctuations of the cosmic microwave background radiation (CMB) demands [5,8]: We consider homogeneous and isotropic matter distribution with the linear equation of state: where w is usually a constant parameter. For non-relativistic matter w = 0, for relativistic matter w = 1/3, and for the vacuum-like state w = −1.
Equation of motion for the curvature which follows from action (2.1) has the form: (2.4) where H =ȧ/a is the Hubble parameter and Γ is the total scalaron decay rate, which is determined by the dominant decay channel. See discussion and the list of references in our works [2,4]. Note that our definition of Γ in the present paper differs by factor 2 from that used in our earlier works. We assume that the scalaron field is homogeneous R = R(t), neglecting small perturbations generated in the course of inflation. After inflation is over, the scalaron field starts to oscillate as where θ is a constant phase determined by initial conditions. This and the subsequent equations are valid in the limit M R t 1, but Γt 1. The Hubble parameter is similar to that at the matter dominated stage (MD), but with fast oscillations around the MD value: The cosmological energy density of matter at this period depends upon the decay width of the scalaron, which in turn depends upon the dominant decay channel.
If there exists scalar particle minimally coupled to gravity, the decay width of scalaron into massless scalars would be: In this case, the energy density of predominantly relativistic matter is equal to: If there are several species of massless scalars, the expressions (2.7) and (2.8) should be multiplied by g S , where g S denotes the number of species. For massive scalar with the mass m s the width of two-body decay would be somewhat suppressed due to the phase space factor proportional to 1 − 4m 2 s /M 2 R . If scalaron predominantly decays into fermions or conformally coupled scalars the decay width vanishes in the limit of massless final state particles and is equal to [5]: where m f is the mass of fermion or conformally coupled scalar. The width is dominated by the heaviest final particle. The corresponding matter density is: (2.10) -3 -Now let us turn to the scalaron decay induced by the conformal anomaly. Production of massless gauge bosons by conformally flat gravitational field was first studied in Ref. [9,10] and applied to the problem of heating in R 2 -inflation in Ref. [7]. 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 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 Fig. 1, 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. 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: (2.12) Correspondingly the energy density of matter created by the decay into this channel would be: It is instructive to compare the energy density of 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. To do that we need to define canonically normalized scalar field Φ linearly connected with R with the kinetic term in the Lagrangian equal to (∂Φ) 2 /2: (2.14) Hence, the energy density of the scalaron field is equal to: where expression (2.5) has been used. This result coincides with the expression for the total cosmological energy density in spatially flat matter dominated universe.
The presented equations are valid if the energy density of matter remains smaller than the energy density of the scalaron until it decays. Comparing Eqs. (2.8), (2.10), and (2.13) with (2.15) we find that in all the cases t cr Γ = 5/6, 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. 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 (GR). 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 two-body reactions between relativistic particles. The typical value of this cross-section at high energies, E m, is [12]: 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 Γ rel ≡ṅ n = σ rel vn rel (2.19) 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 [13][14][15]. Correspondingly we arrive to the following thermal equilibrium condition: Using Eq. (2.17), we find that equilibrium is established at the temperatures below Here we took α R = 0.025 and α = 0.033.
The time corresponding to this temperature is equal to where C is defined in Eq. (2.17). 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/Γ an : (2.24) 3 X-particle production through direct scalaron decay There are two possible channels to produce massive stable X-particles: first, directly through the scalaron decay into a pair XX and another by inverse annihilation of relativistic particles in plasma. Firstly, let us consider the scalaron decay. The probability of the scalaron decay into a pair of fermions is determined by decay width (2.9) 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 Eq. (3.1), n R = ρ R /M R , and ρ R is defined in Eq. (2.15). So Eq. (3.3) turns intoṅ It 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 (2.24) and 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: Consequently, the energy density of X-particles in the present day universe would be: The last approximate equality in the r.h.s. is the condition that the energy density of X-particles is 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. And it does exist. If X-particles are 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: where σ ann v is the thermally averaged annihilation cross-section of X-particles and n eq is their equilibrium number density.
This equation was originally derived by Zeldovich in 1965 [16], and in 1977 it was applied to freezing of massive stable neutrinos in the papers [17,18]. After that it was unjustly named as Lee-Weinberg equation.
The thermally averaged annihilation cross-section of non-relativistic X-particles, which enters Eq. (4.1), for our case can be taken as where the last factor came from thermal averaging of the velocity of X-particles squared, equal to v 2 = T /M X , which appears because the annihilation of Majorana fermions proceeds in Pwave. 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. The expression (4.2) is an order of magnitude estimation. The exact form of the cross-section depends upon particle spins, the type 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 Eq. (2.21). Equation (4.1) 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 Eq. (4.1) numerically. In the limit n X n eq Eq. where the subindex "0" means that the solution is valid for n X n eq , y = M X /T and we have used Eq. (2.17) and the expression for C 0 below this equation.
For the initial temperature we take T in = 6 · 10 −3 M R , according to Eq. (2.21), and T f in = T h = 5.1 · 10 −6 M R (2.24). Correspondingly y in = 1.7 · 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 (4.3): where we have taken g s = 2, β ann = 10 and lastly, according to the line below Eq. (2.21), α = 0.033 and α R = 0.025. The ratio F 2 (y) is depicted in Figs. 2, 3 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 (4.4) is a good approximation to the exact solution. In the opposite case, when y f in > y max , we have to solve Eq. (4.1) numerically.  To solve the equation (4.1) it is convenient to introduce the new function according to: where a(t) is the cosmological scale 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 (4.1) is reduced to: Next, let us change the variables from t to y = M X /T . Evidentlyẏ = −y(Ṫ /T ). Using time-temperature relation (2.17), we find Keeping in mind that Numerical solutions 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 CM B at the preset time, where n CM B = 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 Eq. (4.6). 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 Eq. (4.4), 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 Eq. (4.10). However, the solution demonstrates surprising feature: its derivative changes sign at y 10, when n X n eq , as it is seen from the value of F 2 presented in Figs. 2 and 3. Probably this evidently incorrect result for z(y) originated from a very small coefficient in front of the brackets in Eq. (4.10).
The problem can be avoided if we introduce the new function u(y) according to: . (4.14) 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 Eq.(4.14) and is depicted in Figs. 4 and 5 for different values of µ.   The logarithm of the energy density of X-particles (4.15) with respect to the observed energy density of dark matter as a function of M X is presented in Fig. 6. If M X ≈ 5 · 10 12 GeV, X-particles may be viable candidates for the carriers of the cosmological dark matter.

Conclusion and discussion
There is general agreement that the conventional Friedmann cosmology is incompatible with the existence of stable particles having interaction strength typical for supersymmetry and heavier than several TeV. A possible way to save 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 reduction presents popular now Starobinsky inflationary model [6]. 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 production of radiation from the scalaron decay.
The concrete range of the allowed mass values depends upon the dominant decay node 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 Xparticles, 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. However, so heavy X-particles would decay though formation of virtual black holes, according to Zeldovich mechanism [19]. The corresponding decay life-time would be To make τ X,BH larger than the universe age t U ≈ 5 · 10 17 sec we either need M X < 10 7 GeV, or to apply the conjecture of Ref. [20] which leads to a strong suppression of the decay through virtual black holes for spinning or electrically charged X-particles. If X-scalars are conformally coupled to curvature or X-particles are fermions, then the probability of the scalaron decay is proportional to M 2 X . For sufficiently small M X the production of X-particles would be quite weak, so that their cosmological energy density would be close to the observed density of dark matter if M X ∼ 10 6 GeV [4].
There is another complication due to conformal anomaly, which leads to efficient decay of scalaron into massless or light gauge bosons. There are some versions of supersymmetric theories where conformal anomaly is absent, which were considered in Ref. [4]. In the present work we have not impose this restriction and studied a model with full strength conformal anomaly. In this case the thermalization of the cosmological plasma started from the creation of gauge bosons and the reactions between them created all other particle species.
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. To restrict the density of X-particles produced by the direct decay the observed value M X should be below 10 7 GeV. But 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 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. It opens the possibility for X-particles to make proper amount of dark matter, if their mass is about 5 · 10 12 GeV.
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 chance to discover these particles in accelerator experiments but they may be observable in low background experiments or through the products of their annihilation in galaxies.