Cosmic ray production in modified gravity

This paper is a reply to the criticism of our work on particle production in modified gravity by D. Gorbunov and A. Tokareva. We show that their arguments against efficient particle production are invalid. $F(R)$ theories can lead to an efficient generation of high energy cosmic rays in contracting systems.

in the early universe. It is shown that the former is much more efficient. Section V is dedicated to an important effects of non-harmonicity and deviation from the adiabaticity for stimulation of particle production. In Sec. VI a possible impact of the discreteness of matter is considered and it is found there that such effects do not lead to suppression of the R-oscillations. In Section VII we summarize our counterarguments and conclude.

II. THE MODEL AND BASIC EQUATIONS
In this section we briefly present the main features and notations of the papers [6,7] to make the presentation self-contained. For more details we address the reader to the quoted works.
A. Modified gravity without R 2 term We consider the model with the action: where m P l = G −1/2 N ≃ 1.22 × 10 19 GeV is the Planck mass, R is the scalar curvature, and S m is the matter action. The first term in Eq. (1) is the usual General Relativity (Einstein-Hilbert) action and F (R) is a non-linear function of the scalar curvature. In this subsection, following Refs. [6,7], we take F (R) equal to: Here n is an integer, λ = const. > 0 (in what follows we take λ = 1), and |R c | is of the order of 8π̺ c /m 2 P l , where ̺ c is the present day value of the total cosmological energy density, so R c ∼ 1/t 2 U , where t U ≈ 4 × 10 17 s is the universe age.
This model is not entirely realistic because it requires a past-time singularity [5] and, which is even worse, in systems with rising energy density it can lead to a singularity in a finite time in the future [1][2][3]. Still, investigating this model is instructive, because it is technically much simpler than a more realistic model where the singularities are eliminated by the R 2 term added to the action, see below Eq. (13). Many essential features of the realistic scenario, particularly at small curvatures, can be understood with this simpler model. In what follows we consider objects with |R| ≫ |R c |.
The equations of motion which follow from the action (1) have the form: where F ′ R = dF/dR, D µ is the covariant derivative, and T (m) µν is the energy-momentum tensor of matter. Taking the trace of Eq. (3) leads to the equation: This is a closed equation for R except for metric depending terms in the covariant derivative and in T µ µ . If the metric slightly differs from the flat Minkowski one, equation (4) would contain only one unknown scalar function which completely determines the evolution of R. In this limit the equation can be reduced to the simple oscillator form: for the function where the potential U (w) is defined as: with ν = 1/(2n + 1), q = 2nλ, U ′ (w) = dU/dw, and T = 8πT µ µ /m 2 P l . Moreover, we are considering the case |R| ≫ R c , since in realistic astrophysical systems T ≫ R c . Their ratio is about T /R c ∼ ̺ m /̺ c ≫ 1 and hence w ≪ 1. Thus the first term in square brackets in Eq. (7) dominates. The potential U would depend on time if the mass density of the object under scrutiny also changes with time, T = T (t).
If only the dominant terms are kept in equations (5) and (7) and if the space derivatives are neglected, the equation (5) simplifies to:ẅ It is convenient to introduce the dimensionless quantities: where β and γ are so chosen (see below) that the equation for ζ becomes particularly simple: Here a prime denotes differentiation with respect to τ and the trace of the energy-momentum tensor of matter is parametrised as: The constants γ and β are equal to It was shown in Ref. [1] that in systems with rising mass density the position of the minimum of the potential ζ min = (1 + κτ ) −1/ν moves towards zero and so does ζ(τ ) itself, practically independently on the initial conditions, i.e. on ζ(0) and ζ ′ (0). The function ζ(τ ) oscillates around ζ min (τ ) and at some moment it passes beyond the minimum of U and reaches the point ζ = 0. This value of ζ however corresponds to the singularity R = ∞. The system arrives to the singularity in a finite time, while the external energy density still remains finite. This conclusion is supported by the numerical calculations, shown in Fig. 1. For these particular figures we took ζ 0 = 1 (i.e. the GR value) and ζ ′ 0 = 0. A more rigorous choice of initial conditions is discussed in section III. In what follows we shift the initial time moment in such a way that τ = 0 corresponds to the onset of the rise of density perturbations, keeping the same notation τ for the shifted time.

B. Regularization by R 2
To eliminate the past and future singularities the function (2) was modified by adding a term proportional to the curvature squared [10]: where some parameters are defined below Eq. (2) and m has dimension of energy. The additional term is relevant only at very large curvatures, because m 10 5 GeV is necessary in order to preserve the successful predictions of the standard BBN [11]. In our works [6,7] we have used the model based on Eq. (13). A large m implies that the stabilisation takes place at very high R. Though R does not become infinite, it can reach huge values in systems with rising mass density. This rise normally originates after the onset of structure formation at z ∼ 10 4 or any time later.
We are particularly interested in the regime |R c | ≪ |R| ≪ m 2 , in which F (R) can be approximated by We consider a nearly-homogeneous distribution of pressureless matter, with energy/mass density rising with time but still relatively low (e.g. a gas cloud in the process of galaxy or star formation). In such a case the space derivatives can be neglected and, if the object is far from forming a black hole, the space-time metric is approximately Minkowski. Then Eq. (4) takes the form Let us introduce the dimensionless quantities 1 where ̺ c ≈ 10 −29 g/cm 3 is the cosmological energy density at the present time, ̺ m0 is the initial value of the mass/energy density of the system under scrutiny, and T 0 = 8π̺ m0 /m 2 P l . Next let us introduce the new scalar field: In terms of this field Eq. (15) can be rewritten in the simple oscillator form: where a prime denotes derivative with respect to τ . The potential of the oscillator is defined as: The substitution (17) is analogous to what is done in [1] but now y cannot be analytically expressed through ξ and we have to use two approximate expression for positive and negative ξ, see Eqs. (45a, 45b) below. It is clear that (18) describes oscillations around y = z (the "bottom" of the potential), which corresponds to the usual GR solution R + T = 0. So we can separate the solutions into the average and the oscillatory part. For small deviations from the minimum of the potential, the solution takes the form: where and the dimensionless frequency Ω is defined as taken at y = z. From (18), we find that it is equal to The conversion into the physical frequency ω is given by It is assumed that at the initial moment τ = τ 0 the function ξ(τ ) sits at the minimum of the potential, otherwise we would need to add a cosine term in (17). If initially ξ(τ 0 ) was shifted from the minimum, the oscillations would generally be stronger and the effect of particle production would be more pronounced.

C. Scalaron Potential and Evolution
One cannot analytically invert Eq. (17) to find the exact expression for U (ξ). However, we can find approximate expressions for gy 2n+2 ≪ 1 (ξ > 0) and gy 2n+2 ≫ 1 (ξ < 0). The value ξ = 0 separates two distinct regimes, in each of which Ω has a very simple expression [see Eq. (23)] and ξ is dominated by either one of the two terms in the r.h.s. of Eq. (17). Hence, in those limits the relation ξ = ξ(y) can be inverted giving an explicit expression for y = y(ξ), and therefore the following form for the potential: where By construction U and ∂U/∂ξ are continuous at ξ = 0. The shape of this potential is shown in Fig. 2. The bottom of the potential, as it is obvious from Eq. (19), corresponds to the GR solution R = − T , or y(ξ) = z, and its depth (for gz 2n+2 < 1) is We will use a very simple form for the external energy density z, namely Here, κ −1 and t contr are respectively the dimensionless and physical timescales of the contraction of the system; analogously, τ 0 and t 0 are respectively the dimensionless and physical initial times, which for simplicity and without loss of generality will be taken equal to 0, though in Sec. III we shift the time. It is also useful to express the physical parameters such as m, the initial energy density ̺ m0 , etc., in terms of their respective "typical" values. Let us define where ̺ c = 10 −29 g cm −3 is the present (critical) energy density of the Universe. In terms of these quantities, we can rewrite g and κ as The equation of motion (18) for small oscillations ξ 1 , which is defined in Eq. (20), can be rewritten as The term ξ ′′ min is proportional to κ 2 , which is usually assumed to be 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 an account of this term as well. Using expansion (20) 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 scalaron amplitude α 0 can be expressed through y ′ 0 , which we keep as a free parameter, as: In our works [6,7] we have chosen the initial conditions: |y ′ 0 | κ and |κ − y ′ 0 | ∼ κ. Since, as we know, the evolved collapsing cloud of matter gradually deviates from GR, the closeness of the derivative y ′ 0 to GR at the moment when the oscillating y(τ ) crosses the GR value is quite unnatural. So we disregarded the case |κ−y ′ 0 | ≪ κ. Under these assumptions ξ 1 remains small, ξ 1 ≪ 1, but non-negligible. In this case the numerical results, shown in figure 3, are in remarkable agreement with Eq. (20). From Eq. (33), one is led to think that for y ′ 0 = κ the oscillations would not be excited at all. This is only true at first order in |y ′ 0 − κ|, so their amplitude would be of order κ 2 . It is one of the points of the criticism raised by GT. They claim that the system under scrutiny must start from the state determined by General Relativity (GR). It means that y 0 = 1 and y ′ 0 = κ. According to the GT statement, the curvature scalar in this case would be much smaller than it is found in our works. However, as we shall see below, such initial conditions could be realized only at the very onset of structure formation when the cosmological decrease of density stopped and the contraction started. Indeed, systems with decreasing energy density approach General Relativity and thus the GR initial conditions must be true. Evidently near the minimal density its evolution is quadratic in time and density behaves as δ̺ ∼ (κτ ) 2 . However, it can be shown (see the next section) that the contracting system quickly deviates from GR approaching a state in which density grows approximately linearly with time and ultimately realizing, at some later time, the conditions which we took in our works as the initial conditions. In realistic situation the density rise is closer to the exponential one in accordance with the Jeans law. It leads to somewhat faster excitation of the curvature oscillations and in this sense is more favorable for the discussed effect.
Our primary goal is to determine the amplitude and shape of the oscillations of y. Expanding y as it is easy to show [7] that during the initial oscillating phase, when |β| < z, we have Let us stress that, in contrast to ξ, the oscillations of y quickly become strongly anharmonic and even for slightly negative ξ the amplitude of y may be very large because y ≈ −ξ/g, according to Eq. (17). This feature is well demonstrated by the results of numerical calculations shown in figure 4.
We can estimate the amplitude of the spikes analytically using the energy evolution law: 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 κ is given by Eq. (27). 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. Let us use this law when the minimum value of ξ reaches zero. Evidently in the minimum ξ ′ = 0. Since U (0) = 0, the constant in the r.h.s. of Eq. (36) turns to zero. Now let us go to larger time and neglect the oscillating part of ξ under the integral. The minimum value of ξ (maximum absolute value) ξ min of negative ξ is determined by the equation: The value of z(τ 1 ) ≡ z 1 is found from the condition α = ξ max i.e.
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. (38) gives: This finally determines Let us estimate now the maximum value of R or, what is essentially the same, y max (40). In our works we did it both ways, analytically and numerically. We present analytical asymptotic estimates for a small value of parameter g given by Eq. (16). The numerical solutions has been found in our works [5,6]. The agreement between numerical and analytical results is excellent. However, numerical results can be done only for non-realistic values of parameters, e.g. small g but not so tiny as its real value, but the smaller is g, the better must be asymptotic solution. So we can trust both numerical and analytical calculations.
Analytical results, presented in this work, are copied from our papers [6,7]. The maximum value of R is equal to R max = R GR y max , where y max is given by eq. (40) Here R GR is the value which the curvature would have in the GR limit, i.e. R GR = ̺/m 2 P l , where ̺ is the matter density in the contracting cloud, ̺ can be written as where ̺ c is the cosmological energy density and Q is some enhancement factor of the rise of mass density in the contracting cloud. It is reasonable to expect Q ∼ 10 5 . Here the numerical factors of order ten are omitted. They are not important.
Next we use expression (16) for g, Eq. (27) for κ and taking z 1 ≈ (κ) −2/(3n+1) in Eq. (39) we find for the maximum amplitude at the spike R max = R GR y max : where σ = 2n/(3n + 1). The last factor is so small, that for any reasonable values of parameters R << m 2 and thus the theory remains below its ultraviolet cut-off.

III. INITIAL CONDITIONS
As mentioned above, the disagreement with GT is about the initial curvature velocity y ′ 0 , which in turn determines the initial amplitude of the scalaron oscillations, α 0 . We assumed that α 0 ∼ O(κ), whereas GT argued that, because the GR solution is y ′ 0 = κ, α 0 should be of order κ 2 . Note that (31) vanishes for y ′ 0 = κ due to Eq. (33). Evidently, the source term in the r.h.s. of Eq. (30) produces oscillations and hence deviations from GR of order κ 2 with any initial conditions, including y ′ 0 = κ. Using Eqs.(23), (31), and (35), we observe that which grows for increasing z regardless of the value of α 0 , whether ∼ κ or ∼ κ 2 . Similarly, deviations in the first time derivative y ′ from z ′ = κ increase, and at some moment we will inevitably have |y ′ 0 − κ| ∼ κ. Redefining the initial time as such moment, we conclude that the estimates of our works [6,7] are indeed reliable.
Let us study now the evolution of the system from the moment when the expansion of the cloud changed into contraction, so at t = t in = 0 the density reaches a minimum value ̺ = ̺ 0 . As is well known, the Tolman solution for an over-density, which at some moment (which we fix at t = 0) decouples from the cosmological expansion and starts contracting, is: which is clearly different from the linear dependence (11). However, we can show that this difference does not undermine the validity of our results.
To avoid misunderstanding, expressed by one of the referees, let us stress again the following. The difference between our way of fixing the initial conditions and that of Refs. [8,9] is that we made an order of magnitude guess for the value of, say, y ′ at a moment t 1 , while GT suggested that a physically justified way is to fix y ′ at a previous time moment t 2 . We agree with the GT choice, but using it and numerically proceeding from t 2 to t 1 , we have found that the value of y ′ at t 1 is roughly the same as it was guessed in our work. So our guessed initial condition, though taken without rigorous justification, happens to be correct. Here and below the notation t is used instead of τ , though y is a dimensionless curvature which depends on dimensionless time.
We agree that the initial conditions should be rigorously fixed at the minimum value of ̺. But using this prescription we have observed even more efficient generation of the curvature oscillations than in our more naive approach. So the production of cosmic rays by the mechanism of our papers [6,7] is by no means suppressed in contrast to the claim of papers [8,9]. For negative t, the density ̺ drops down with t → 0. At this stage y(t) tends to its GR value, while y ′ (t) tends to zero independently on the initial conditions taken at negative t. It is well known that R(t) in the systems with decreasing ̺(t) approaches the GR value with decreasing speed. So a study of the cosmological history implies the following initial conditions for Eq. (18): On the other hand, in systems with rising ̺ the curvature R runs away from GR and oscillates around the GR value with rising amplitude. As mentioned above, R even tends to infinity in finite time if the R 2 /m 2 term is absent. Hence at that stage the derivative y ′ can naturally reach values of order κ (or even higher). For a qualitative understanding of the solution we first consider the case when the R 2 term is absent and the evolution of R(t) or, to be more precise, of ζ(τ ) is governed by equation (10) with the unique analytic expression for the potential at positive ζ. We remark that when ζ reaches zero the curvature becomes infinite. The treatment of this case is technically much simpler than that with R 2 present, but this simpler case well illustrates the evolution of R(t) in the more realistic model determined by the F (R) of Eq. (13).
We solve numerically equation (10) taking the initial value of ζ at the minimum of the potential, that is ζ 0 = ζ min = (1 + κτ ) −(n+1) with some initial τ in not necessarily equal to zero. Note that ζ = ζ in corresponds to the GR value of curvature. We took several interesting initial values of ζ ′ 0 , i.e. of the "velocity" of ζ with respect to the bottom of the potential. The least favorable case for the excitation of the curvature oscillations is ζ ′ 0 = −(2n + 1)κ, when initially ζ runs with the same speed as the bottom of the potential. The results are presented in figs. 5. We see that in all the cases ζ reaches zero roughly speaking at κτ ∼ 1. It means that the In all the cases ζ initially sits in the minimum of the potential, ζ0 = ζmin, i.e it takes the GR value, while the initial "velocities" ζ ′ 0 relative to motion of the minimum of the potential are different by magnitude and direction. The least favorable case is when initially ζ moves together with the minimum of the potential, i.e. ζ ′ 0 = −(2n + 1)κ. Still in all the cases the singularity ζ = 0 is reached after a time such that κτ ∼ O(1). singularity R = ∞ or, better to say, high spikes of R are quickly reached. Now let us turn to the physically justified initial conditions taken at the moment when the energy density given by the expression (43) acquired the minimal value, in other words, we choose the initial conditions (44) taken at the initial time t in = 0. We calculate how ζ evolves starting from these initial conditions and check to which values we arrive going to the time when the density rises according to the law (11). In this way we find that at that time the "velocity" ζ ′ 0 is generically of order κ but not close to the special least favorable value ζ ′ 0 = −(2n + 1)κ. But we have seen above that even in this case the formation of spikes with large amplitude is sufficiently fast to ensure efficient particle production.
We found no essential difference between the evolution of ζ starting from many different initial values for ζ ′ , in particular from those determined by the Tolman solution, recommended in Refs. [8,9]. Thus we conclude that the correct choice of the initial conditions at the onset of the density increase does not inhibit the efficiency of the curvature rise in contrast to the claim of Refs. [8,9]. In other words, the choice of the initial conditions based on the Tolmann model is quite good, even more favorable for creation of large amplitude curvature oscillations than the initial condition ζ ′ ∼ κ, as it was taken in our works [6,7].
We have done similar calculations taking ̺(t) evolving according to the Tolmann solution in the theory with the R 2 -regulator. (Notice that we have different notations for the unknown function: ζ if R 2 is absent and ξ if R 2 is present.) The equation of motion in presence of the R 2 -term becomes much more complicated because we have to use different expressions for the potential at positive and negative ξ.
The equation of motion still has the simple form (18,19) but y has different expressions through ξ at positive and negative ξ. In the limit of small g we have: y = ξ −1/(2n+1) , if ξ > 0, and ξ = −gy, if ξ < 0, so there are the following two equations in the two regions: The first equation is solved numerically, leading to a result quite close to the solution of Eq. (10) for the case without R 2 . The initial conditions for this equation are determined from the Tolman solution in the same way as is done above, i.e. the initial values are taken according to Eq. (44). The evolution of ξ(t) does not differ much from the previous case without the R 2 -term as long as ξ remains positive. After a while at the moment t = t 1 the almost harmonically oscillating ξ crosses zero and we have to use equation (45b). The initial values for the solution of this equation are taken as ξ(t 1 ) = 0, by definition, and the initial value ξ ′ (t 1 ) is determined from the numerical solution of Eq. (45a). For κ ≈ 0.003 the initial value of the derivative is calculated to be ξ ′ (t 1 ) ≈ −0.1. For smaller κ e.g. for κ = 0.0003 the derivative ξ ′ (t 1 ) is about −0.03. In all the cases the spikes of R at negative ξ are not suppressed in comparison with our simple estimates in papers [6,7]. (We do not distinguish here between physical time t and the dimensionless one, τ . Hopefully it will not lead to confusion.) To avoid possible misunderstandings let us stress again that the dynamical initial time is t = 0, while t 1 corresponds to the initial time taken in our works [6,7]. We have shown that the order of magnitude guess of the initial conditions taken in Refs. [6,7] is very well supported by the calculated evolution from the rigorously found initial conditions at t = 0 determined according to the Tolman solution.
If g is much smaller than unity, Eq. (45b) can be quite accurately solved analytically. To this end we rewrite it in the form: where ξ 1 = ξ + g(̺ m /̺ 0 ). Since the duration of the spike is very short, i.e. about g 1/2 , we can neglect the time variation of (̺ m /̺ 0 ) during this period and the solution can be easily found. For example if we take the initial conditions as discussed above, namely ξ(t 1 ) = 0 and ξ ′ (t 1 ) = −0.1 the solution has the form: In this way we reconstructed the solution found in our papers [6,7] in the limit of very small g obtaining a spike with very high amplitude. We want to stress again that in our original works we fixed the initial conditions when the perturbations have already evolved and are growing roughly linearly with time. We agree that the initial conditions can be rigorously fixed at the minimum of ̺. Nevertheless using this prescription we have found a similar or even more efficient generation of the curvature oscillations than in our more naive approach. So the production of cosmic rays by the mechanism of our papers [6,7] is by no means suppressed in contrast to the statement of the papers [8,9].
Thus the rigorous fixing of the initial conditions based on the Tolman-type solution leads to essentially the same huge amplitude of curvature oscillations as it has been found in our works [6,7] with an intuitive choice of the initial conditions. We conclude that the criticism of [8,9] related to our choice of the initial conditions is irrelevant. We can summarize our counter-argument on this issue as follows: • indeed the exact GR initial conditions lead to a smaller initial amplitude of the oscillations of both y and ξ, but oscillations of y grow at the same relative rate. Therefore, at some later moment we will have a more sizeable |y ′ 0 − κ| ∼ κ, and our results are recovered. • the fact that the oscillations of y grow in time is independent of the assumed form z(τ ), as long as z ′ > 0 (contraction), as shown numerically for a different functional form (43) in Fig. 6.

IV. GENERATION OF CLASSICAL CURVATURE FIELD VERSUS QUANTUM PRODUCTION OF SCALARONS
The second point of GT criticism refers to the probability of the scalaron production. According to their estimates [9] "one cannot expect such a large contribution of scalarons because less than one particle inside the horizon may be created in the present (or recent) Universe, see [8]. Scalarons created in the very early Universe were very heavy (with mass m) and hence decayed to the SM particles. So the initial conditions in [7] seem to be irrelevant." Out of that GT concluded that the initial number density of scalarons in the present day universe must be very small. Evidently it is not so. The production of heavy quanta in the very early universe has nothing to do with light particle production at the present epoch. In Ref. [8] GT also considered particle production at the present day universe in the objects with rising energy density due to gravitational (Jeans) instability. We comment on this at the end of this Section.
In our works [6,7] we studied the generation of the classical field R(t) in the sufficiently late or present day universe after the onset of the structure formation. We have done explicit calculations of generation of the classical curvature field R(t), which are not related to the remnants of the early produced heavy scalaron quanta. In quantum language, which is perfectly compatible with the classical approach, "our" quanta of scalarons, produced at the present epoch, are light. Their number can easily be huge in the same way as low-energy radio waves correspond to a huge number of massless photons.
Clearly the heavy scalaron particles produced in the early universe would not survive to the present day, though they may be the source of all matter in the universe e.g. in R 2 -inflation. On the other hand, our mechanism of classical generation of R(t) field may result in a noticeable contribution to energetic cosmic rays. The classical generation of R(t), of course, agrees with quantum production rate of light quanta, though technically the calculations are more complicated. For example, the standard way to treat the universe heating after inflation is to calculate particle production by a classical oscillating scalar field (inflaton). It is equivalent to the consideration of the decay of the inflaton condensate.
These facts in no way undermine the validity of our results, which apply to regions of high and increasing density. As we have shown explicitly in the previous section, even starting from purely GR initial conditions, the increasing external density acts as a source term for curvature oscillations and moves solutions away from GR. In such dense, contracting systems, where ̺ m /̺ c ≫ 1, certainly a contribution even much higher than 10 −4 ̺ c 2 would not affect other physical phenomena significantly, because these increase is only local, but gives much weaker contribution to the average cosmological density.
In Ref. [8] GT have considered particle production in the Einstein frame, while our calculations have been done in the Jordan frame. They managed to simplify the equation of motion for the scalaron mode function to a simple linear equation, Eq. (17) of their work. However, in this approximation all essential effects of our approach are missing. The high production probability advocated in our papers is induced by the non-linear term in the equation of motion and by the strong non-harmonicity of the curvature oscillations.

V. EFFECTS OF NON-HARMONIC OSCILLATIONS OF CURVATURE
The authors of Refs. [8,9] write: "... production of such energetic particles in a slow and smooth process of contraction looks very surprising from any point of view..." Another quotation from GT: "This result [of the efficient particle production [6,7]] is rather unexpected, because high frequency oscillations (i.e. heavy particles) are produced by a slow process (the structure formation) very inefficiently. Moreover, cosmological evolution naturally gives zero initial amplitude for such oscillations, which can be associated only with the scalar mode (heavy scalaron)." In our opinion this statement is grossly incorrect. Perhaps the picture, that Gorbunov and Tokareva have in mind, originated from the assumption of harmonic oscillations of the scalaron field whose frequency and amplitude vary adiabatically. In [7], we have indeed computed the predicted particle production in the initial, harmonic oscillatory phase, and found it to be negligible. However, the behaviour of R(t) in the potential (25) is neither necessarily "slow" nor "smooth" at all times. The fact that the statement of GT is incorrect is especially striking if the R 2 -term is absent and R tends to a singularity. The presence of R 2 prevents the formation of a curvature singularity but does not make the process smooth and slow. The temporal evolution of R(t) is drastically different from the adiabatically slow changing quasi-harmonic oscillations.
Our approach to particle production is very much different from that by GT. We considered the particle creation by the classical R(t)-field, while GT calculated the decay of the scalaron quanta, bearing in mind that the number of these quanta is negligibly small. Our results strongly contradict the GT statement: "In all the realistic situations the scalaron production is very inefficient. Always less than one particle inside the corresponding Jeans volume is produced. Subsequent scalaron decay contribution to the cosmic ray flux is found to be in infinitesimal." Clearly the oscillating curvature scalar which is the solution of our equations (45a, 45b) corresponds to a very large number of the produced scalarons at rest. Even with the initial conditions which according to GT lead to a strong suppression of our results, the produced number density of scalarons is much larger than GT evaluated by their quantum production mechanism.
According to GT the "Scalarons created in the very early Universe were very heavy (with mass m) and hence decayed to the SM particles. So the initial conditions in [6,7] seem to be irrelevant." (we repeat the quotation here to stress a different point of disagreement). Indeed the decay rate of scalarons in the early universe, as calculated in Refs. [12][13][14] and later confirmed in [11,15], is equal to: With such a decay rate the expected density of scalarons should quickly drop down, if they were created in the very early universe. The similar conclusion might be valid for the classical curvature oscillations. However, this is not so in our scenario. As explained above, GT use arguments valid in a cosmological scenario but apply them to a different physical situation, namely contracting clouds with ̺ ≫ ̺ c . In our picture the curvature oscillations were induced at a late cosmological stage when the frequency of the oscillations in the quasi-harmonic regime (i.e. when ξ is positive) is much lower than m. So the field is relatively light, and the corresponding life-time is long. In the spike regime the frequency is large but the duration of the spike is too short to create a relevant damping. Moreover, the expression for Γ R (48) is true only for harmonic oscillations [11]. Furthermore, a very important fact is that there is a constant energy supply for the oscillations of R(t) coming from the cloud of the collapsing matter. This slow process induces the low frequency oscillations which transform their energy into high frequency modes (spikes) due to the non-linearity of the equation. Because of that the damping of the oscillations due to particle production is negligible. The validity of these statements is supported by numerical calculations.

VI. EFFECTS OF DISCRETENESS OF MATTER
As was noticed by GT there is another mechanism which might lead to the suppression of particle production by high frequency oscillations of R (spikes), namely, the inhomogeneous distribution of the background matter consisting of separate particles. As written in Ref. [9]: "The approximation that is certainly doubted is the homogeneous energy-momentum tensor for the background matter. If one considers a set of discrete particles instead of continuous medium one obtain that they produce the scalaron field like point sources. Every particle produces the field ξ ∼ exp(−ωr)/r like the usual massive scalar with mass ω. If the average distance d between particles is large, d ≫ 1/ω then the scalaron field is strongly inhomogeneous and approximation used for the energy-momentum tensor is not valid. In this case we expect that the energetic particle production may still take place only in a small region around the point sources, i.e. in a very small spherical regions of radius 1/m. Therefore, the total number of produced particles is suppressed by a very small number n/m 3 where n is the matter number density." We have already explained that the process is not smooth and slow but the problem of discrete "constituents" of the cosmological matter is certainly of interest. Before addressing this problem we note that this problem is not inherent to our mechanism. Indeed, the averaging problem and the matching between the small scale picture of point-like particles to a continuous perfect fluid (CDM) at larger scales is a problem of great interest even in GR. Needless to say, this problem is well beyond the scope of our original works. Therefore, as in all previous works on f (R) gravity before us, we assume that a matter distribution homogeneous on average is for all practical purposes equivalent to a truly homogeneous density distribution.
Let us stress again that similar or exactly such problem exists in the standard cosmology. The universe is described in the ideal liquid approximation with the average density distribution. This is true in the case of linear equations, when the averaging of the equation leads to the equation for the average values. However, this is not so for the nonlinear GR equations. Nevertheless no other work found/assumed any difference between "strictly homogeneous" and "homogeneous on average" in the models of modified gravity (and there are very good reasons that the equality holds rather well in GR, probably thanks essentially to Gauss's Law). So, unless GT show explicitly that this equivalence is wrong, we see no reason why we should not use a fluid description for Dark Matter.
Moreover, the average distance between the matter particles is model dependent and might be sufficiently short or they may be even overlapping in the quantum mechanical sense. This allows to avoid the problem indicated by GT. As a possible counter-example we can mention the case of Bose-condensate dark matter, composed of light (e.g. axions) or heavy scalars, where the particle constituting dark matter forms a continuous, nearly homogeneous condensate.
We can also consider normal dark matter particles or neutral atoms of the usual baryonic matter and recall that they are not point-like particles but quantum mechanical wave packets. The wave packet size is determined by the coherence length of the particle propagation in cosmic gas/plasma so the packets would easily overlap and the quantum mechanical density matrix of the cosmological medium would be sufficiently homogeneous.
The size of the wave packet of a particle produced by the decay of a pure state of some parent particle is equal to the parent particle lifetime. For example the size of the wave packet of a muon created in the pion decay π → µν is macroscopically large, L = τ = 2.6 · 10 −8 s ≈ 70 m. However, it is realistically much shorter because of the breaking of coherence of the pion wave packet, see e.g. [16]. For example, if the pion is in a thermal bath, then the coherence is maintained for no longer than a typical time δt ∼ v/l f ree , where l f ree is the mean-free path of the source particle and v is the relative (Möeller) velocity of the colliding particles. Bearing in mind applications to the cosmological plasma dominated by relativistic matter we take v = 1. So it seems natural to suggest that the wave packet size in thermal gas is determined by the particle mean free path in the gas, where this particle is created and propagates (see e.g. [17,18] and references therein). So the packet size is: where σ tot and n are respectively the total cross-section of the particle interactions with the cosmic plasma and the plasma density. If we take some heavy matter particles, their decoupling from the rest of the plasma, or better to say, gas takes place when: where H ∼ T 2 /m P l is the Hubble parameter, T is the plasma temperature and m P l is the Planck mass. Assuming that the plasma is dominated by relativistic matter we take v = 1. The plasma density is n ∼ T 3 and the interaction cross-section is σ tot ∼ α 2 /m 2 DM , where m DM is the mass of the dark matter particle and α is the coupling constant, which is typically small, α 10 −2 . Correspondingly we find for the decoupling temperature: Hence the wave packet size at the decoupling would be: The wave packet would be stretched out by the cosmological expansion and by the quantum mechanical packet spreading. The cosmological expansion gives the factor z = T d /T 0 = m 2 DM /(α 2 T 0 m P l ), where T 0 is the CMB temperature at the present day, so the wave packet size would become: The cosmological number density of dark matter particles is equal to n DM ≈ (z + 1) 3 (keV/m DM )/cm 3 , where z is the redshift with respect to the present time. Correspondingly the average distance between the DM-particles is d DM = (m DM /keV) 1/3 (z + 1) −1 cm. So Evidently this ratio is very small for any reasonable m DM and thus the dark matter particles in cosmology are not well separated rigid "stones" but overlapping quantum states with very homogeneous density matrix. However, if the dark matter particles are primordial black holes, then they would not overlap. Now let us show that the GT claim that the distance between particles d is much larger than the interaction scale in the harmonic regime, d ≫ 1/ω (see the beginning of this section), is grossly overestimated. Indeed, according to Ref. [19]: Yukawa profile outside the particles. The first issue is indeed very interesting and relevant even in GR, but we opted to follow all literature before us in considering that one can indeed approximate a distribution homogeneous on average as truly homogeneous. As for the second assumption, it is strictly true only in the quasi-static approximation which is clearly violated in our quickly oscillating, spiky solutions, so we do not believe it to be a crucial point against our results.
Moreover, the particles constituting the system might not be classical but quantum particles whose wave packet size is much larger than ℓ. There are also scenarios with Bose-Einstein condensate dark matter for which the density distribution is continuous and the GT arguments are not applicable.
3. There is an essential difference between our works and the critical scenario by GT. We solved the classical equation of motion governing the evolution of practically massless field R(t) and calculated particle production by this oscillating field. Naturally in a quantum language the number of the R-quanta in this classical state can be large. Such situation is realised for instance in the process of emission of an electromagnetic wave, which carries an almost infinite number of photons. The solution of the classical equations is straightforward and simple, and the number of the produced quanta is huge for any initial conditions, ours and those taken by GT.
4. One more critical argument by GT is the statement that an adiabatic slow process cannot lead to an efficient creation of the classical field. However, this reasoning might be applied only to slow harmonic oscillations, while in our case the oscillations are strongly non-harmonic and non-adiabatic.