There are two popular phenomenological explanations of the observed accelerated universe expansion. It is either prescribed to the existence of the so called dark energy (DE) with the equation of state

$$ P\approx - \rho, $$
(1)

where P and ρ are respectively pressure and energy densities. Another way to describe cosmological acceleration is a modification of the classical action of general relativity (GR) by an addition of a non-linear function of curvature, F(R):

$$ A =-\frac{m_{\mathrm {Pl}}^2}{16\pi}\int d^4x\sqrt{-g} \bigl[R+F(R) \bigr] , $$
(2)

where the 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

$$ 3\mathcal{D}^2 F'_R-R+RF'_R-2F= 8\pi T^\mu_\mu/m_{\mathrm {Pl}}^2, $$
(3)

where \(\mathcal{D}^{2}\equiv \mathcal{D}_{\mu} \mathcal{D}^{\mu}\) is the covariant D’Alambertian operator, \(F'_{R}\equiv d F/ d R\), and T μν is the energy–momentum tensor of matter.

To describe the astronomical data regarding cosmological acceleration the solution of the equation

$$ R_c-R_c F'_R (R_c) + 2 F(R_c) = 0 $$
(4)

should be equal to \(R_{c} = - {32\pi \varOmega_{\lambda}\rho_{c}}/{m_{\mathrm{Pl}}^{2}}\), where Ω λ ≈0.75 is the fraction of the observed vacuum-like energy density and ρ c ≈10−29 g/cm3 is the total cosmological energy density.

The pioneering suggestion [1, 2] of gravity modification with F(R)∼μ 4/R suffered from strong instabilities in celestial bodies [3]. To cure this instability further modifications of gravity have been suggested [4–6]. These and some other versions are reviewed e.g. in Refs. [7–9]. The suggested modifications, however, may lead to infinite-R singularities in the past cosmological history [7] and in the future in astronomical systems with rising energy/matter density [10–13]. 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 [14–17]. Possibly in astrophysical systems a 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 [6]:

$$ { F(R) = -\lambda R_0 \biggl[1- \biggl(1+ \frac{R^2}{R_0^2} \biggr)^{-n} \biggr]-\frac{R^2}{6m^2}, } $$
(5)

where n is an integer, λ>0, and \(| R_{0} | \sim 1/t_{U}^{2}\), where t U ≈13 Gyr is the universe age. The parameter m is bounded from below, m≳105 GeV, to preserve successful predictions of BBN, see e.g. [20].

Other models [4, 5] 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 work [17–19]. There was renewed interest to this problem [20, 21] 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 the equation of motion (3) contains non-linear derivative terms:

$$ \bigl[(\partial_t R)^2 - (\partial_j R)^2 \bigr]/R, $$
(6)

which makes it very inconvenient for qualitative analysis. To avoid that, we introduce, instead of R, the new function, proportional to F′(R):

$$ \xi\equiv \frac{1}{2\lambda n} \biggl(\frac{T_0}{R_0} \biggr)^{2n+1}F'_R, $$
(7)

where \(T_{0} = 8\pi \rho_{0}/m_{\mathrm{Pl}}^{2} \). It is also convenient to introduce dimensionless quantities z=ρ(t)/ρ 0, y=−R/T 0, g=1/[6λn(mt U )2](ρ m0/ρ c )2n+2, and dimensionless time \(\tau = m t \sqrt{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:

$$ \xi''+ dU/d\xi =0 , $$
(8)

where a prime denotes derivative with respect to τ and dU/dξ=z−y(ξ). Unfortunately y cannot be expressed through ξ analytically. We have only the inverse relation

$$ \xi = {1}/{y^{2n+1}} -gy , $$
(9)

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 agreement of the analytical results with numerical ones where the latter can be done allows to conclude that analytical estimates can be valid for high frequency and small g as well.

The minimum of the potential U(ξ) is located at y(ξ)=z(τ), so it moves with time according to

$$ \xi_{\min} (\tau) ={z(\tau)^{-(2n+1)}}-gz(\tau). $$
(10)

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

$$ \varOmega^2 = \frac{\partial^2 U}{\partial \xi^2}\biggm{|}_{y=z} = \biggl(\frac{2n+1}{z^{2n+2}}+g \biggr)^{-1} . $$
(11)

Note that the physical frequency is \(\omega = \varOmega m \sqrt{g}\).

For small oscillations ξ=ξ min+ξ 1 takes the form

$$ \xi(\tau)\approx \xi_{\min}(\tau)+\alpha(\tau)\sin \bigl[F(\tau)+\delta \bigr] $$
(12)

where δ is a phase determined by the initial conditions and

$$ F (\tau) \equiv \int^\tau_{\tau_0}d\tau'\, \varOmega\bigl(\tau'\bigr). $$
(13)

For positive and negative ξ the potential can be approximated as

$$ U(\xi) = U_+(\xi)\varTheta(\xi) + U_-(\xi) \varTheta(-\xi) , $$
(14)

where

(15)

There is a kind of the conservation law for the energy of field ξ:

$$ \frac{1}{2} \xi'^2 + U(\xi) - \int^\tau_{\tau_1}d\tau' \frac{\partial z}{\partial\tau'} \xi\bigl(\tau'\bigr)=\mathrm{const}, $$
(16)

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 the assumption made above z linearly grows with time as z(τ)=1+κτ, where

$$ \kappa = (m t_{\mathrm{contr}} \sqrt{g})^{-1}. $$
(17)

This simple law may be not accurate when t/t contr>1, but probably the results obtained are not too far from the realistic case.

Equation of motion (8) for small oscillations ξ 1 can be rewritten as

$$ \xi_1''+ \varOmega^2\xi_1= -\xi_{\min}'' . $$
(18)

Term \(\xi_{\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

(19)

Here and in what follows sub-0 means that the corresponding quantity is taken at initial moment τ=τ 0.

We impose the following initial conditions:

$$ \begin{cases} y(\tau=\tau_0)=z(\tau=\tau_0)=1 ,\\ y'(\tau=\tau_0)=y'_0 , \end{cases} $$
(20)

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

$$ { \alpha_0 = \bigl(\kappa - y'_0\bigr) (g+ 2n+1)^{3/2} }. $$
(21)

To keep ξ 1 small the initial value of the derivative \(y'_{0}\) should be also small. In this case the numerical results, shown in Fig. 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.

Fig. 1
figure 1

Oscillations of ξ 1(τ), in the case n=2, κ=0.01, g=0.01 and initial conditions y 0=1, \(y'_{0}=\kappa/2\). The amplitude of the oscillations is in very good agreement with analytical estimate (12)

Full size image

Our primary goal is to determine the amplitude and shape of the oscillations of y. In contrast to ξ the oscillations of y are strongly anharmonic 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 Fig. 2.

Fig. 2
figure 2

“Spikes” in the solutions. The results presented are for n=2, g=0.001, κ=0.4, and \(y'_{0}=\kappa/2\). Note the asymmetry of the oscillations of y around y=z and their anharmonicity

Full size image

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

(22)

The value of z(τ 1)≡z 1 is found from the condition α=ξ min i.e.

$$ z_1^{-(2n+1)} - g z_1 = \big|\kappa - y_0' \big| ( g +2n +1 )^{5/4} \biggl( \frac{2n+1}{z_1^{2n+2}} + g \biggr)^{1/4} . $$
(23)

In the limit of small g, \(g< 1/z_{1}^{2n+2}\), asymptotically |ξ max|=(g/n)1/2 z(τ 1)−n. In the same limit Eq. (23) gives

$$ z_1 = \bigl[ \bigl( \kappa - y_0' \bigr)^2 (2n+1 )^3 \bigr]^{-1/(3n+1)}. $$
(24)

This finally determines

$$ y_{\max} = (ng)^{-1/2} z_1^{-n}. $$
(25)

This result is much larger than the naive estimate \(g y_{\max}^{2n+2} \sim 1\) [13]. For negative ξ potential (15) behaves as U≈ξ 2/(2g), so the characteristic frequency of the oscillations in the region of negative ξ is about \(\varOmega \sim 1/\sqrt{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 [22].

According to calculations of Ref. [20], harmonic oscillations of curvature with frequency ω and amplitude R max transfer energy to massless particles with the rate (per unit time and volume):

$$ \dot{\rho}_{PP}\simeq {R^2_{\max} \omega}/{(1152\pi)} , $$
(26)

The life-time of such oscillations is \(\tau_{R} = {48 m_{\mathrm {Pl}}^{2}}/{\omega^{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 a slowly varying amplitude B(t), superimposed on smooth background A(t):

$$ R(t) = A(t)+B(t)\sum_{j=1}^{N} \exp \biggl[-\frac{(t-jt_1)^2}{2\sigma^2} \biggr] . $$
(27)

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 [22]). Finally we find

$$ \bigl \vert \tilde{\mathcal{R}}(\omega)\bigr \vert ^2\simeq \frac{4\pi^2 B^2 \sigma^2e^{-\omega^2\sigma^2}\Delta t}{t_1^2}\sum_j\delta \biggl(\omega- \frac{2\pi j}{t_1} \biggr) . $$
(28)

Identifying B with R max=y max T 0 and integrating over frequencies we obtain

$$ \dot{\rho}_{PP} \simeq {y_{\max}^2 T_0^2}/({576\pi t_1}) . $$
(29)

Time interval t 1 is approximately equal to \(2\pi/\omega_{\mathrm{slow}}= 2\pi/(\varOmega_{\mathrm{slow}}m \sqrt{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

$$ \dot{\rho}_{PP} = C_n \frac{\rho_0^2 (m t_U)^2 z^{n+1}}{m_{\mathrm{Pl}}^4 t_{\mathrm{contr}}} \biggl( \frac{t_U}{t_{\mathrm{contr}} } \biggr)^{\frac{n-1}{3n+1}} \biggl( \frac{\rho_c}{\rho_0} \biggr)^{\frac{(n+1)(7n+1)}{3n +1}}, $$
(30)

where

$$ C_n = (2n+1)^{\frac{9n-1}{2(3n+1)}} ( 6\lambda n )^{\frac{7n+1}{2(3n +1)}}/ (18 n). $$
(31)

It is convenient to present numerical values: \(\rho_{c}/ m_{\mathrm{Pl}}^{4} \approx 2\times 10^{-123}\) and \((mt_{U})^{2} \approx 3.6\times 10^{93} m_{5}^{2}\), where ρ c ≈10−29g/cm3 and m 5=m/105 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:

$$ \rho_{CR} \approx 10^{-24} m_5^2 z^{n+1} \ \mathrm{GeV}\,\mathrm{s}^{-1}\,\mathrm{cm}^{-2}. $$
(32)

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 the assumptions R>R 0 and ρ 0>ρ c . However, if we take ρ 0 somewhat larger than ρ c the results would be approximately correct.

For m=1010 GeV the predicted flux of cosmic rays with energy around 1019 eV is close to or even higher than the observed flux.