Curvature oscillations in modified gravity and high energy cosmic rays

It is shown that F(R)-modified gravitational theories lead to curvature oscillations in astrophysical systems with rising energy density. The frequency and the amplitude of such oscillations could be very high and would lead to noticeable production of energetic cosmic ray particles.

There are two popular phenomenological explanations of the observed accelerated universe expansion.It is either prescribed to existence of the so called dark energy (DE) with the equation of state: where P and ̺ are respectively pressure and energy densities.Another way to describe cosmological acceleration is a modification of the classical action of the general relativity (GR) by an addition of a non-linear function of curvature, F (R): where function F (R) is chosen in such a way that the modified GR equations have a solution R = const in absence of any matter source.Taking the trace g µν δA/δg µν = 0, we find: where D 2 ≡ D µ D µ is the covariant D'Alambertian operator, F ′ R ≡ dF/dR, and T µν is the energymomentum tensor of matter.
To describe the astronomical data about cosmological acceleration the solution of the equation should be equal to R c = −32πΩ λ ̺ c /m 2 P l , where Ω λ ≈ 0.75 is the fraction of the observed vacuumlike energy density and ̺ c ≈ 10 −29 g/cm 3 is the total cosmological energy density.
The pioneering suggestion [1] of gravity modification with F (R) ∼ µ 4 /R suffered from strong instabilities in celestial bodies [2].To cure these instability further modifications of gravity have been suggested [3][4][5].These and some other versions are reviewed e.g. in refs.[6,7].The suggested modifications, however, may lead to infinite-R singularities in the past cosmological history [6] and in the future in astronomical systems with rising energy/matter density [8,9].These singularities can be successfully eliminated by an addition of the R 2 -term into the action.Such a contribution naturally appears as a result of quantum corrections due to matter loops in curved space-time [10,11].Possibly in astrophysical systems singularity may be avoided even without R 2 if one takes into account distortion of the background flat Minkowsky space-time by an impact of rising R.However, it would take place if R very much differs from its GR value and such a deviation from normal gravity would be observable.
In what follows we consider the version of the modified gravity suggested by Starobinsky [5]: where n is an integer, λ > 0, and |R 0 | ∼ 1/t 2 U , where t U ≈ 13 Gyr is the universe age.Parameter m is bounded from below, m 10 5 GeV, to preserve successful predictions of BBN, see e.g.[14].
Other models [3,4] demonstrate similar behavior.We study the evolution of R in a contracting astrophysical system and show that R oscillates around its GR value with quite a large amplitude, efficiently producing elementary particles with masses below the oscillation frequency.The particle production damps the oscillations and may suppress or even eliminate the singularity, which would appear if R 2 term was absent.
Cosmology with R 2 term and in particular the particle production by the oscillating curvature was studied in the early works [11][12][13].There was renewed interest to this problem [14,15] stimulated by the study of the interplay of the infrared and ultraviolet terms in eq.(5).
In this paper we study solutions of eq.(3) in the system of non-relativistic particles with rising energy density which we approximate as ̺ = ̺ 0 (1 + t/t contr ), where t contr is the effective time of contraction.We assume that the gravity of matter is not strong and thus the background metric can be considered as flat.Written in terms of R equation of motion (3) contains non-linear derivative terms: which makes it very inconvenient for qualitative analysis.To avoid that, we introduce, instead of R, the new function, proportional to F ′ (R): where , and dimensionless time τ = mt √ g.In terms of these quantities and in the limit of R ≫ R 0 and sufficiently homogeneous ξ (both conditions are naturally fulfilled) the equation of motion for ξ takes a very simple form: where prime denotes derivative with respect to τ and dU /dξ = z − y(ξ).Unfortunately y cannot be expressed through ξ analytically.We have only inverse relation so we have to use different approximate expressions in different ranges of ξ.
In what follows we make analytical estimates of the amplitude of the curvature oscillations and compare them with numerical calculations for some values of the parameters.The agreement is generally quite good.However, we cannot do numerical calculations for realistically tiny values of g or huge frequencies of oscillations due to numerical instability and/or too long computational time.So an agreement of the analytical results with numerical ones where the later can be done allows to conclude that analytical estimates can be valid for high frequency and small g as well.
The minimum of potential U (ξ) is located at y(ξ) = z(τ ), so it moves with time according to It is intuitively clear that even if initially ξ takes its GR value ξ = ξ min it would not catch the motion of the minimum and as a result it starts to oscillate around it.Dimensionless frequency of small oscillations, Ω, is determined by: Note that physical frequency is ω = Ω m √ g.For small oscillations ξ = ξ min + ξ 1 takes the form: where δ is a phase determined by the initial conditions and For positive and negative ξ the potential can be approximated as: where There is a kind of the conservation law for the energy of field ξ: where τ 1 is an arbitrary fixed time moment.The last term appears because U explicitly depends on time through z.If ∂z/∂τ is positive, which is the case for a contracting body, the value of U (ξ) would in general grow with time.According to assumption made above z linearly grows with time as z(τ ) = 1 + κτ , where This simple law may be not accurate when t/t contr > 1, but probably the results obtained are not too far from realistic case.
Equation of motion (8) for small oscillations ξ 1 can be rewritten as Term ξ ′′ min is proportional to κ 2 , which is usually assumed small, so in first approximation it can be neglected, though an analytic solution for constant Ω or in the limit of large Ω can be obtained with this term as well.Using expansion (12) and neglecting α ′′ , we obtain Here and in what follows sub-0 means that the corresponding quantity is taken at initial moment τ = τ 0 .We impose the following initial conditions the first of which corresponds to GR solution at the initial moment.The initial value of the amplitude derivative α 0 can be expressed through y ′ 0 , which we keep as a free parameter, as: To keep ξ 1 small the initial value of the derivative y ′ 0 should be also small.In this case the numerical results, shown in figure 1, are in remarkable agreement with eq. ( 12).For y ′ = κ the oscillations seem not to be excited.However, this is an artifact of the approximation used.For example, the "source" term in the r.h.s. of eq. ( 18) creates oscillations and hence deviations from GR with any initial conditions.
Our primary goal is to determine the amplitude and shape of the oscillations of y.In contrast to ξ the oscillations of y are strongly unharmonic and for negative and even very small |ξ| the amplitude of y may be very large because y ≈ −ξ/g, according to eq. ( 9).This feature is well demonstrated by the results of numerical calculations shown in figure 2.
We can estimate the amplitude of the spikes using the energy evolution law (16).Let us use it at the moment when the maximum value of ξ reaches zero.So at this moment ξ ′ = 0. Since U (0) = 0, the constant in the r.h.s. of eq. ( 16) turns to zero.Now let us go to larger time and neglect the oscillating part of ξ under the integral.The maximum absolute value of negative ξ is determined by the equation: In the limit of small g, g < 1/z 2n+2 1 , asymptotically |ξ max | = (g/n) 1/2 z(τ 1 ) −n .In the same limit eq. ( 23) gives: This finally determines This result is much larger than the naive estimate gy 2n+2 max ∼ 1 [9].For negative ξ potential (15) behaves as U ≈ ξ 2 /(2g), so the characteristic frequency of the oscillations in the region of negative ξ is about Ω ∼ 1/ √ g in dimensionless time which corresponds to ω ≈ m in physical time.Evidently frequency of oscillations of y in this region is the same.It can be shown that the calculated amplitude of y (25) becomes noticeably larger with rising z [16].
According to calculations of ref. [14], harmonic oscillations of curvature with frequency ω and amplitude R max transfer energy to massless particles with the rate (per unit time and volume): The life-time of such oscillations is τ R = 48 m 2 P l /ω 3 .In our case the oscillations are far from harmonic and we have to make Fourier expansion of the spiky function y(τ ).To this end we approximate the "spike-like" solution as a sum of gaussians of amplitude modulated by slowly varying amplitude B(t), superimposed on smooth background A(t): We assume that σ ≪ t 1 , that is the spacing between spikes is much larger than their width.The Fourier transform of expression ( 27) is straightforward but rather tedious (details can be found in our work [16]).Finally we find: Identifying B with R max = y max T 0 and integrating over frequencies we obtain Time interval t 1 is approximately equal to 2π/ω slow = 2π/(Ω slow m √ g), where Ω slow is given by eq. ( 11).Amplitude y max is defined in eq. ( 25); parameter κ is determined by eq. ( 17) and "coupling constant" g is defined below eq.( 7).Taking all the factors together we finally obtain: where , where ̺ c ≈ 10 −29 g/cm 3 and m 5 = m/10 5 GeV.Now assuming that the particle production lasts during t ≈ t contr and taking ̺ 0 = ̺ c , we find the energy flux of cosmic rays produced by oscillating curvature: This result is a lower limit of the flux of the produced particles.With larger z when the minimum of the potential shifts deep into the negative ξ region the production probability significantly rises.Analytical evaluation in this case is more difficult but simple qualitative arguments indicate to it and numerical calculations supports this assertion.Recall that our derivation has been done under assumption R > R 0 and respectively ̺ 0 > ̺ c .However if we take ̺ 0 somewhat larger than ̺ c the results would be approximately correct.
For m = 10 10 GeV the predicted flux of cosmic rays with energy around 10 19 eV is close or even higher than the observed flux.