Derivative couplings in gravitational production in the early universe

Gravitational particle production in the early universe is due to the coupling of matter fields to curvature. This coupling may include derivative terms that modify the kinetic term. The most general first order action contains derivative couplings to the curvature scalar and to the traceless Ricci tensor, which can be dominant in the case of (pseudo-)Nambu-Goldstone bosons or disformal scalars, such as branons. In the presence of these derivative couplings, the density of produced particles for the adiabatic regime in the de Sitter phase (which mimics inflation) is constant in time and decays with the inverse effective mass (which in turn depends on the coupling to the curvature scalar). In the reheating phase following inflation, the presence of derivative couplings to the background curvature modifies in a nontrivial way the gravitational production even in the perturbative regime. We also show that the two couplings -- to the curvature scalar and to the traceless Ricci tensor -- are drastically different, specially for large masses. In this regime, the production becomes highly sensitive to the former coupling while it becomes independent of the latter.


Introduction
Our current understanding of cosmology is based on the ΛCDM model, which introduces a dark matter sector to explain the dynamics of the Universe. This dark matter sector is yet poorly understood as the fundamental nature of its constituents remains still a mystery. Due to the seemingly negligible interactions between dark matter and the Standard Model (SM) sector [1,2], many authors have proposed in the last decades that dark matter may have been produced gravitationally during the early stages of the Universe [3][4][5][6][7][8]. Gravitational production has been studied since the pioneer work of Parker [9] in which the effect of spontaneous particle production by the expansion of the Universe was first studied. In general, we know now that particle creation is a general feature of quantum field theory in curved spacetimes which is the paradigm we use to study the early epochs of the Universe.
Inflation [10,11] provides the seeds for structure formation giving a nearly scaleinvariant power spectrum of fluctuations due to the quantum fluctuations of the inflaton field [12,13]. The bridge between inflation and the usual hot big bang scenario is obtained by an intermediate phase known as reheating in which the energy stored in the inflaton field is transmitted to the Standard Model degrees of freedom. Particle production due to the expansion of the Universe is specially relevant during these two stages of the Universe and the transition between them. The density of created particles during these primordial cosmological epochs could potentially explain the dark matter abundance that is currently observed.
With this goal in mind, we will consider a relativistic massive scalar field nontrivially coupled to gravity. In particular, we will include terms involving the field ϕ, its derivatives ϕ ,µ , the Ricci tensor R µν , and the scalar curvature R in the action. The Rϕ 2 term has been largely studied in the previous literature [3][4][5][6][7][8], but terms containing derivatives of the field such as ϕ ,µ ϕ ,ν R µν and ϕ ,µ ϕ ,µ R have not been studied in this gravitational production context as far as we know. The main goal of this work is to study particle production in the early universe in the presence of these derivative couplings. The existence of nonsuppressed nonminimal couplings of this kind appear in the context of (pseudo-)Nambu-Goldstone bosons associated with the breaking of global symmetries [14,15].They are also generic for disformal scalar fields [16], that are coupled to the Einstein tensor. An example of this type of fields is provided by branons in flexible brane-world models [17][18][19][20][21]. The phenomenology of such couplings have been extensively studied in colliders [19,[22][23][24][24][25][26], with precision observables [25,26], and with astrophysical observations [24,[27][28][29].
This manuscript is organised as follows. In section 2 we introduce the model of a massive scalar field with a linear coupling to the Ricci scalar and linear derivative couplings to both the Ricci tensor and scalar, neglecting the back-reaction of the field on the background and the interaction with the SM particles. We then obtain the equation of motion for this field in a Friedmann-Lemaître-Robertson-Walker background. We quantise it and briefly describe the formalism of particle production in such curved spacetimes. In section 3 we solve this equation for de Sitter geometry as a special case to see how the derivative couplings affect particle production. As we will see, in this case the addition of derivative couplings does not alter the type of equation for the modes of the field, simply introducing an effective mass which cannot be distinguished experimentally from the physical mass. In section 4 we investigate the gravitational effects in a simplified model for the reheating background, allowing us to delve into the effects of the derivative couplings. The investigation in the reheating epoch is carried out numerically and we use an interpolation regime for the transition between the de Sitter phase and the reheating one. We conclude in section 5 with a summary of the main results obtained in this work.

Particle creation in curved spacetimes
Let us consider a massive scalar field with several direct linear couplings to the background curvature. The dynamics of this field in a given spacetime is encoded in the action where g is the determinant of the spacetime metric g µν , R is the Ricci curvature scalar, ξ is a dimensionless constant, m is the mass of the field ϕ, and S der will be defined shortly. Apart from this latter term this is the well-known action for a nonminimally coupled scalar field without modifying the kinetic term. The derivative-coupling contribution to the action where R µν is the Ricci tensor and γ and σ are parameters with units of inverse energy squared, has been written in this particular form in order to separate the part that describes the coupling to the curvature scalar from the part that contains the coupling to the traceless Ricci tensor. In fact, as we will see, these two contributions play an entirely different role in the particle production.
In cosmological scenarios the background geometry is well described, in average, by the Friedmann-Lemaître-Robertson-Walker metric. For simplicity and in agreement with the current cosmological observations [30], we are considering flat spatial sections. Moreover, for later convenience, we will write the FLRW metric using conformal time: For this specific geometry, the Ricci tensor is diagonal and its spatial components, due to homogeneity and isotropy, are all equal. It can be expressed in terms of the scale factor a(η), the Hubble parameter H = a /a 2 , and the curvature scalar R = 6a /a 3 as where the prime denotes derivative with respect to the conformal time η. From the action (2.1) we obtain the equation of motion for the field ϕ: where, for convenience, we have defined the functions An appropriate rescaling of the field variable allows us to eliminate the friction term. Indeed, in terms of the rescaled field χ(x, η) = f (η)ϕ(x, η) with the equation of motion for χ reads As a general comment, from (2.6) and (2.9) we see that some terms of the field equation vanish if σ = −4γ, including a term proportional to the second derivative of the curvature scalar R . In this case the field couples only to the Einstein tensor, which is also the case for disformal scalar fields and branons [16][17][18][19][20][21]. Moreover, additional significant simplifications may appear for specific functional forms of the scale factor. For instance, in a radiationdominated phase (a ∝ η), the curvature scalar vanishes; and in a de Sitter phase, which mimics the dynamics of inflation (a ∝ η −1 ), all the curvature terms are proportional to the (constant) Hubble parameter in this phase.
Exploiting the fact that the spatial sections are flat, we can expand the field χ in Fourier amplitudes which satisfy harmonic oscillator-like equations of motion where k = |k|, with time-dependent frequency We can write any solution as a linear combination of two complex conjugate solutions {v k (η), v * k (η)}. Due to the isotropy of the background geometry, we can choose bases of solutions which only depend on k. Then the Fourier amplitude χ k can be expressed as a linear combination of the mode function v k and its complex conjugate, the coefficients being the creation and annihilation variables a k and a * −k : In order to preserve the standard Poisson bracket structure for the annihilation and creation variables, the mode functions must be normalised. We will follow the convention of Ref. [13]: (2.14) The Fock quantisation procedure (see e.g. Ref. [31]) can be carried out by promoting the annihilation and creation variables a k and a * k to creation and annihilation operatorŝ a k andâ † k satisfying the canonical commutation relations The Hilbert space for the states of the field is spanned in the usual way by the action of finite linear combinations of products of these operators on the vacuum state, defined as the one which is annihilated by all the annihilation operators, i.e., Generally speaking, in curved spacetimes, the definition of the vacuum state for the field is not unique but suffers from different ambiguities that can lead to potentially inequivalent quantum theories. Different choices of mode functions in (2.13) will define a different set of creation and annihilation operators and hence, different vacuum states. As each basis of mode functions is complete, different bases must be related by a linear transformation, i.e., if we consider a different set of mode functions {u k , u * k }, we can write v k (η) = α k u k (η) + β k u * k (η), (2.17) with α k , β k complex coefficients satisfying appropriate normalisation conditions that ensure the normalisation of the new modes. This leads to new creation and annihilation operators related to the previous ones by the so-called Bogoliubov transformations, The vacuum states corresponding to each quantisation |0 a and |0 b are related by provided that both quantisations are unitarily equivalent [32][33][34], that is to say, if and only if k |β k | 2 < ∞. The Bogoliubov coefficients α k and β k can be straightforwardly computed from (2.17) imposing the Wronskian condition (2.14): If we introduce the particle number operatorN (a) = kâ † kâ k (and likewise for the quantisation (b)), it is easy to see that the total number of particles a 0|N (b) |0 a corresponding to the quantisation (b) in the vacuum |0 a , is divergent because the spatial volume is infinite. Introducing the particle-number densityn (a) to avoid this problem, it can be shown that (2.21) In cosmology, the gravitational production occurs because of the time evolution of the background geometry. In these scenarios, the vacuum state for a quantum field is not stationary. Particles are produced because the evolved vacuum state at a given instant of time is in general an excited state with respect to the instantaneous vacuum defined at that instant of time.
We can define the vacuum state of the scalar field through the initial conditions given in the asymptotic past. These initial conditions define a set of modes v k and v * k to be discussed in the following sections. On the other hand, a comoving observer will have a different instantaneous notion of vacuum. We can use the so-called adiabatic prescription [35] to define this instantaneous vacuum. This prescription is appropriate and well defined if the geometry evolves slowly enough as compared with the characteristic time scales of the field, i.e., whenever the frequency ω k (η) is a slowly changing function during the time ∆η = 1/ω k (adiabatic condition): The corresponding modes are then given by the WBK approximation to first order [35] at a given η 0 , which satisfy the following initial conditions:

Particle creation in de Sitter spacetime
In order to explain the current cosmological observations with the standard cosmological model, people usually introduce an inflationary epoch prior to the usual big bang model. For our purposes, we can consider this inflationary epoch as a pure de Sitter geometry, which can be described with a scale factor a(η) = −1/(H 0 η). Here the Hubble parameter H 0 is constant and the conformal time η lies on the interval (−∞, 0). In more realistic scenarios, the background geometry is governed by a slowly varying Hubble parameter in a quasi-de Sitter expansion. As the de Sitter geometry is a maximally symmetric solution, the Ricci scalar is proportional to the metric, R µν ∝ Rg µν . Therefore, in this background, the coupling to the derivatives of the field depends only on the single parameter γ. This is not surprising because for the de Sitter geometry the traceless Ricci tensor is identically zero. Computing the curvature related quantities for this scale factor and substituting them in (2.4) and (2.12), we obtain the time-dependent frequency for the field during the inflationary phase: where we have defined Note that µ behaves as a (constant) effective mass. Thus we see that the effect of the different considered couplings to the curvature is to modify the effective mass of the field during the inflationary phase. Inspecting the possible values of the coupling constants, we note that: 1. The effective mass µ becomes imaginary for certain combinations of ξ and γ. In fact, the potential in these cases is not even bounded from below. These scenarios correspond to tachyonic states for the field.
2. For m = 0, ξ = 1/6, γ = 0 the effective mass vanishes and we recover the well-known result that there is no production of purely conformal massless particles.
Effectively, we have only a single parameter to characterise the scalar field, meaning that we cannot discriminate between different sets of the parameters m, ξ and γ leading to the same µ in pure de Sitter spacetime. The adiabatic condition (2.22) in this scenario is fulfilled either for large effective masses (µ 1) or whenever the combination k|η| 1. From (3.2) we see that a positive value of the ξ parameter tends to move the system towards the adiabatic regime, whereas nonvanishing values of the derivative coupling parameter γ will move it away from this regime.
The solution to the mode equation (2.11) can be written in terms of Hankel functions [36]: The so-called Bunch-Davies vacuum is defined by mode functions that in the asymptotic past η → −∞ behave as plain waves [37] in analogy with the definition of vacuum in Minkowski spacetime. Demanding that in the asymptotic past the vacuum state only has positive-frequency modes, and taking into account the asymptotic behaviour of the Hankel functions for |η| → ∞ [36], the mode function v k must have the form (3.3) with In order to simplify the computations and for later convenience, we will perform a change of variables from the set (k, η) to the new set (y, η), where the new variable is defined as y = k|η|/|ν|. In this new set of variables, the modes defining the Bunch-Davies vacuum read: It is now easy to see, after some algebra, that in terms of the independent variables (y, η) the Bogoliubov coefficients β k (η) := β(y) depend only on the variable y and not on y and η independently.
Once the β(y) Bogoliubov coefficient is computed, the created number density of particles is obtained by integrating its modulus squared over all the possible momenta and dividing by the fiducial volume a(η) 3 : Note that this expression is independent of the conformal time η. This is a remarkable result: the density of produced particles in cosmological de Sitter spacetime is constant over time for a comoving observer in the adiabatic regime.
In figure 1, we show the spectral particle production for small and large masses as compared with H 0 at the time η = −1/H 0 . In the left panel we see that as the mass increases within the sub-Hubble regime, the amplitude of the production peak decreases while the resulting spectrum broadens, in such a way that the total production grows as the mass increases. This result can be easily understood because the particle creation is a consequence of the conformal symmetry breaking, and the parameter controlling the degree of breaking is precisely the mass. On the other hand, for masses beyond the Hubble parameter H 0 (right panel), the peak amplitude also decreases as the mass increases but the spectrum broadens much less as compared with the left panel (not enough to overcome the dampening of the amplitude). The global result is quite different: the net production decreases with increasing masses. We can understand this effect if we think of the Hubble parameter as an indicator of the energy stored in the gravitational field. Then, creating particles more massive than the energy of the background is more difficult.
We want to analyse the asymptotic behaviour of the number density (3.6) for the regime of large effective masses µ 1. In this regime, ν is pure imaginary, i.e. ν = i|ν|. Hence, we will perform a series expansion of the mode solutions in 1/|ν|. One could be tempted to ignore the contribution 1/4 to ν but it turns out to be crucial. Then, the asymptotic series for the Hankel function reads [36] H (1) When we compute the coefficient β(y), the first and second terms of the series expansion cancel out, so the leading contribution to the particle density comes from the third order terms in 3.7. Then, the leading order for |β(y)| 2 is It should be noted that for masses larger than a few times H 0 the agreement with the exact spectrum is almost perfect. Finally, integrating (3.6) using the residue theorem, we find that the behaviour of the particle density in the large effective mass regime is 4 Particle creation during reheating

Background evolution
Our current understanding of the early universe includes two phases before the usual big bang cosmology: inflation and reheating. The dynamics of both phases can be described as a FLRW model sourced by a scalar field which has a slow-roll regime, giving rise to the inflationary phase; it also has a minimum in its potential where it oscillates provoking its decay into the particles that permeate the universe nowadays in the process known as reheating. As we are interested in the effects of the derivative couplings discussed above in the gravitational production during the early universe, we need to take into account not only the inflationary phase, which is well described for our purposes with a de Sitter phase, but also we need to include the reheating epoch. We are mimicking the background spacetime dynamics of both phases using a de Sitter phase for inflation and a matterdominated universe as the reheating phase. This choice of a matter-dominated solution comes from the mean behaviour of the background during the oscillations of the inflaton field φ in the massive chaotic inflationary scenario with a potential ∝ m 2 φ φ 2 [38]. We can use the parameterisation for the scale factor introduced in [39] which provides a continuous and differentiable function: where H 0 = √ 12πm φ is the value of the Hubble parameter at the onset of the single-field inflationary epoch for the model we are considering and we have set the end of inflation at η = 0 (with m φ being the inflaton mass). From now on we will use H 0 as the basic quantity for units. This parameterisation is good enough to get a continuous Hubble parameter. However, it fails to provide continuous and differentiable curvature scalar and tensor. In order to have a satisfactory model, we need the frequency of our field to be sufficiently smooth. Hence, we need to introduce an interpolation regime for the curvature quantities between the end of inflation (i.e., the end of the de Sitter phase) and a moment η * in the reheating epoch. This time η * is chosen by equating the Hubble parameter from the above parameterisation and the one obtained through the dynamics of the inflaton field, which for any single-inflation model can be written during the oscillatory phase of the field around the minimum of its potential as (4.2) where the inflaton field has been expanded at first consistent WKB order. In these expressions, M P = 1.22 · 10 19 GeV is the Planck mass, m φ = 10 13 GeV is the inflaton mass, and Φ 0 = M P /(2 √ π) is the amplitude of the field at the onset of the oscillations. These values for the inflaton mass and amplitude are constrained by the CMB observations [30].
Equating (4.2) and the Hubble parameter from the parameterisation of the scalar factor, we determine η * = 2.4/H 0 .

Particle production
In a model driven by an inflaton field, all the background quantities would be sufficiently smooth functions. Hence, we want to interpolate R and R µν so they reproduce this differentiability to sufficiently high order. In order to do so, we introduce polynomials of degree 24 for the value of the functions as well as for each of their derivatives independently, so each of them is a C 12 function. The degree of the polynomials has been chosen so it is the minimum possible one so the results are independent of the interpolation order. Note that without the interpolation, the frequency would exhibit a discontinuity when the de Sitter epoch finishes. This discontinuity would abruptly excite modes of arbitrarily high momenta on the field and hence the expected production of particles would be greater than the one obtained considering the background evolution due to an inflaton field.
The equation of motion for the modes during the post inflationary epoch, after the junction time η * , acquires the form where the frequency of the modes has a rather complicated functional form in terms of the background quantities: where A, B, and F are given in (2.6) particularised for the background described above.
The expression of the frequency as a function of the conformal time η is obtained using the Hubble parameter, the Ricci tensor, and the curvature scalar coming from the averaged behaviour of the inflaton field during the reheating epoch. We will solve the equation for the modes of the field numerically, using as initial conditions the value of the field and its first derivative at the end of the de Sitter phase, given by (3.3) with the constants given by (3.4). Despite the cumbersome expression for the frequency, a direct evaluation of the adiabatic parameter |ω k /ω 2 k | yields that we safely lie in an adiabatic regime on the reheating phase for all the region of the parameter space we will consider. Therefore the definition of the adiabatic vacuum is well suited to define a proper instantaneous vacuum to compare with the evolved initial state. Now we are ready to compute the gravitational particle production during the early universe using the model we have discussed, delving into the influence of the derivative couplings in this process. Once we have computed the numerical mode functions which define the initial vacuum state, we can compute the Bogoliubov coefficients relating this initial vacuum state and the adiabatic ones during the post inflationary epoch via equation (2.20).
We will consider the derivative couplings as perturbations, meaning that their contribution is subdominant with respect to the one coming from the mass term and the coupling to the Ricci scalar. Furthermore, even for perturbative values of these derivative couplings there is a strong limitation coming from the de Sitter phase, as each case in which their contributions lead to imaginary effective masses, i.e., µ 2 < 0, the vacuum state of the field would become unstable [35]. Hence, we have an additional constraint on the derivative couplings given by the inequality:

Perturbative analysis
For sufficiently small values of the derivative couplings γ and σ, we can perform an expansion of the particle production. As we argue in the following this requires two different regimes of approximation: a perturbative expansion of the action that can still induce nonlinear effects due to the nontrivial kinetic term, and an additional expansion for small derivative couplings as compared with the mass (in appropriate units). In this subsection, we will highlight the main aspects of these perturbative and linear expansions for the particle production.
There are reasons that suggest that we should restrict the values of the derivative couplings to the perturbative regime. For instance, these derivative couplings affect the kinetic term of the field, modifying the propagator in a highly nontrivial way. In particular, the derivative couplings of the scalar field to the Riemann tensor at first order are not renormalisable. If such terms are present, higher order terms are expected. This means that if the first order contributions are not small, neglecting higher terms may not be well-motivated.
In view of the action (2.1), we can qualitatively advance that, in the perturbative regime, the coupling σ multiplied by the largest value of the traceless Ricci tensor along the whole period relevant to particle production must be much smaller than 1. Similarly, γ multiplied by the largest value of the curvature scalar must be much smaller than 1. More explicitly, it is easy to see that for metrics of the form (2.3) the factor multiplying to ϕ 2 is proportional to 1 + A and the factor multiplying to ∂ i ϕ∂ i ϕ is proportional to 1 + B, where A and B are given by in (2.6). The perturbative regime is then achieved when |A| and |B| are much smaller than 1 during the whole evolution. Considering the maximum values of the background quantities multiplying γ and σ, we can write the perturbative conditions in terms of conditions on the maximum absolute value of the derivative couplings: (4.6) Here and from now on we will use the rescaled dimensionless derivative couplingsγ = 12H 2 0 γ andσ = H 2 0 σ, and the dimensionless mass parameterm = m/H 0 . Even in this perturbative regime, the way in which the quantity A enters the definition of the quantisation field χ makes the final result nonlinear in the parametersγ andσ. In particular, the derivative couplingγ contributes to the effective mass of the field, even during the de Sitter phase, whileσ does not (see e.g. Eq. (3.2)). In the following we focus on the case ξ = 1/6 both for simplicity and to isolate the effects of the derivative couplings.
Let us check the conditions on the derivative couplings that are necessary for the changes on the particle production to be small. To begin with, we expand the frequency to first order in the derivative couplings, If we impose the condition that the linear term Ω k is small compared with the unperturbed frequency squared ω 2 0k , once we evaluate the background geometrical quantities accompanying the derivative coupling at their maximum value during the whole evolution, we can extract a necessary condition on the derivative parameters: These conditions are necessary but not sufficient to have small changes in the particle production |β k | 2 as we discuss in the following. Let us start will the zeroth order |β 0k | 2 . This is calculated using the zeroth order term u 0k of the adiabatic modes u k (i.e., (2.23) with ω 0k instead of ω k ) and the zeroth order v 0k of the vacuum modes v k . These zeroth order modes are solutions to the equation v 0k + ω 2 0k v 0k = 0. The contribution to |β k | 2 linear in the derivative couplings can also be easily found, although the calculations are quite tedious and not very illuminating. Let us describe the ingredients involved in the final expression other than the zeroth order quantities already mentioned, namely, ω 0k , u 0k , v 0k , and their derivatives. These extra ingredients are the already calculated linear order contribution to the frequency Ω k and its derivatives, the corresponding perturbations to the adiabatic modes, and the perturbation ζ k to the modes v k , defined as v k = v 0k (1 + ζ k ). This is the only nontrivial quantity that we have to calculate. Taking into account the equation that it must satisfy (namely, (2.11) with ω k expanded to first order), it is easy to see that ζ k is given by the quadrature (4.9) In view of this expression, we can assume that ζ k (η) is small throughout the evolution provided thatγ andσ satisfy some extra conditions. For illustrative purposes let us study the large mass regime. Taking into account the form of Ω k and a quick numerical estimation of the values of the time integrals involved, we obtain the following rough estimation of the maximum value throughout the evolution for ζ k (η): which apply for large masses. The mass-dependent factors come from the double time integral. The huge difference in their values can be traced back to the fact that the dominant evolution quantity accompanyingγ is proportional to H(η) 2m2 and the one accompanyingσ is proportional to H (η)m 2 : the former has a nonvanishing value during the de Sitter phase and decays afterwards while the latter vanishes in the de Sitter phase and grows in absolute value from there to vanish again. The double integration with the highly oscillatory function v 0k (η) 2 (with a mass dependent frequency) then reduces the amplitude of theγ term by a factor 1/m, while that ofσ has as an additional factor 1/m 7/2 . The linear modification of |β k | 2 also depends on the derivative of the mode but nothing new comes from there. We learn two things from this upper bound to ζ k . First, the contribution of the term linear inγ is small in the perturbative regime. Second, in the same regime, the contribution ofσ is negligible. It turns out that, for large masses, although the linear contributions are small in the perturbative regime as we have argued, this does not mean that the quadratic terms in the derivative couplings are negligible because the mass-dependent factors accompanying them can be very large and in fact dominate. As we will see in the next section, this is indeed the case for theγ contribution. On the other hand, we will also see that the quadratic terms inσ are small in the large mass regime, provided that (4.8) are satisfied, becoming absolutely negligible. For small masses, the linear analysis is less conclusive. However it seems to indicate that small values of σ can lead to significant changes in the production.
Performing the second order expansion in order to check the conditions for the linear regime to be valid is calculationally very expensive, more so taking into account that we can provide numerical results to all orders and directly see the role of each of the couplings to the particle production. In fact, as already suggested above, the most interesting results are obtained within the perturbative regime but for nonlinear contributions to the production. This is done in what follows.

Results
We will start our analysis discussing the effects of the different couplings on the spectral particle distribution for the field. We will cover the perturbative range for the derivative couplings limited by the restriction (4.5) that ensures that there are no tachyons in the theory but allowing for nontrivial behaviours of the particle production. Figure 2 shows the spectral particle production for small masses (we usem = 0.1 as the working value). The left panel shows the production for various values of the derivative couplingγ keepingσ = 0, while the right panel shows the production for various values ofσ keepingγ = 0, so we can study the effect of each coupling separately. Both panels are very different as far as the considered range of the derivative couplings is concerned. In the left panel, the main graphics correspond to values of theγ coupling (not much) smaller than the boundaries (4.8) of the regime in which we expect small variations in the production. The inset displays the production for couplings beyond this regime but still well inside the perturbative limit (4.6). In the right panel, however, the main graphics correspond to values of theσ coupling close to 1. These values are inside both the perturbative regime and the range of (4.8). The inset displays the production for couplings also inside this regime but with significant nonlinear changes. This choice of values for the derivative couplingσ (associated with the traceless Ricci tensor) is motivated by the fact thatσ has a much smaller effect on the particle production (as compared withγ). This is suggested by the arguments and bounds discussed in the previous subsection and confirmed by the numerical results. Indeed, for values ofσ comparable to those ofγ in the left panel, the result would be just that all curves are superimposed. Even more, as we can see in figure 2 the change in the spectral production for values ofσ close to 1 are of the same order as for values ofγ inside the regime (4.8). The specific ranges that we have used for these couplings lie in the intervalsσ ∈ [−1/2, 1/2] in the main graphics. For the inset we have considered the values σ = ±2, still inside both the perturbative limit and the regime (4.8), but for which we obtain nontrivial contributions. We have consideredγ ∈ [−m 2 /2,m 2 /2] ⊂ (−1, 1) (with the regime in which we expect no significant changes being |γ| m 2 as deduced from (4.8)). Note that we cannot go beyondγ = −m 2 due to the restriction (4.5) even being within the perturbative regime. For the curves on the inset we have considered the values γ = (2m 2 , 3m 2 ) for which we expect nontrivial contributions to particle production, though still well inside the perturbative limits.

Spectral particle production
Let us start with the main graphics on the left panel (varyingγ). It is not surprising that within the considered range forγ the changes in the amplitude of the spectra are linear in the value of the coupling. We see that for any value ofγ the spectral production has the same dominant peak as in the conformal case and that its amplitude grows for negative values ofγ and diminishes for positive values. There are two a priori factors behind these changes in the peak amplitude. First, the curvature scalar (coupled throughγ) contributes to the value of the effective mass of the field during the de Sitter phase diminishing and broadening the spectral production for larger effective masses (see left panel in figure 1). Second,γ also contributes through an additional pure curvature effect, present only after the de Sitter phase, when the curvature scalar evolves with time. Let us briefly discuss this mechanism in more detail. If we consider the evolution after the de Sitter phase with only a mass term equivalent to the effective mass (including the effect ofγ on the mass, i.e., µ(m,γ), but neglecting any other effect during the evolution) we see that the behaviour of the amplitude of the spectra withγ is inverted with respect to that in figure 1, with all the peaks centred around the same momentum and with similar width. This is precisely the opposite behaviour to that in the left panel of figure 2. For instance, note that the purple curve (corresponding toγ = −m 2 /2) gives an effective conformal behaviour (µ = 0) in the de Sitter epoch, hence no particle production occurs during that epoch, and still the spectral production is significantly augmented during the transition phase.
The inset in the left panel shows the spectra for situations beyond the regime (4.8) but still for perturbative values of the couplingγ. We can see that the effect of the coupling to the curvature scalar induces significant changes to the spectral production as a secondary peak emerges at higher momenta. This secondary peak dominates the production for values of the coupling beyond the linear regime.
The main graphics on the right panel (varyingσ close to 1 and withγ = 0) also show a dominant production peak for momenta that were sub-Hubble during inflation. In this case, there is no contribution of the derivative couplings to the effective mass µ and the whole effect is due to pure tensor (i.e. traceless) Ricci curvature enhancements. We also see that the production for large momenta is significantly enhanced by the very same mechanism over a broad region of the spectra, augmenting the production even further. This effect is specially important on the red and purple curves which correspond to the larger values of σ within the considered range. In the inset we can see that the spectral production grows wildly, despiteσ being inside the regime (4.8), as discussed in the previous subsection. For positive values ofσ there are two contributions: the primary peak corresponding to the de Sitter production enhanced by the subsequent evolution and a broad band on larger momenta, while for the negative values all the production is concentrated on the latter broad band.
As a final remark, we can conclude that the coupling to the traceless Ricci tensor contributes significantly less than the coupling to the curvature scalar as we have already inferred from the perturbative analysis of the contributions to the modes.
In figure 3 we show the effects of the derivative couplings on the spectral particle production for large masses of the field (we have consideredm = 10 as a representative value). The left panel shows the production for various values of the derivative couplingγ forσ = 0, while the right panel shows the production for different values ofσ keepinḡ γ = 0. In both cases, the main graphics correspond to values of the considered derivative coupling smaller than the perturbative boundaries (4.6). In particular, the parameters lie in the intervalγ ∈ [−0.01/2, 0.01/2], andσ ∈ [−1/2, 1/2], respectively (the same as in figure 2). Note that these values are inside the regime (4.8), where we expect small linear The figure on the left panel depicts the dependence of the spectral particle distribution on different values of the derivative couplingγ, within the perturbative regime, forσ = 0. We see that for negative values of the coupling constantγ, the contribution of the curvature scalar is to enhance the particle distribution obtained at the end of the de Sitter phase, as the spectral production is dominated by the same single peak. On the other hand, for positive values of the couplingγ to the curvature scalar we see the same effect as for couplings beyond the linear regime for small masses (inset of figure 2-left), meaning that the de Sitter peak is damped and that there appears a broad band of production on larger momenta, becoming the dominant contribution to the spectral production. As happened for small masses, there are two different factors to explain these contributions: the changes in the de Sitter effective mass and the pure curvature effects during the transition phase. It is not difficult to see, by direct comparison with the right panel of figure 1, that the wild variation on the amplitude of the peaks has its major contribution from the effect of the curvature variation during the transition between inflation and reheating. It is important to note that even though we are considering values of the coupling constant well inside the perturbative regime and much smaller than the mass of the field, the resulting contributions to the particle distribution are nonlinear in the strength of the couplingγ, as advanced in subsection 4.3.
On the right panel, the main figure shows that the effect of the traceless Ricci tensor mediated by theσ coupling on the spectral production is negligible for values of the coupling within the perturbative limit inasmuch as all the curves lie one on top of the other. This is precisely the behaviour expected from (4.10), where theσ linear contribution to the mode perturbation is very much suppressed by high powers of the mass.
Comparing the results for large masses with our prior analysis of the small mass regime we can see two important differences: on the one hand, the importance of the derivative couplingσ diminishes, going from giving important contributions to the spectral density (figure 2-right) to becoming negligible in the large mass scenario ( figure 3-right). On the other hand, we see that for the same values of the couplingγ to the curvature scalar, the contributions of this term go from being linear for small masses to become nonlinear in the large mass regime. Hence, we can conclude from the analysis of the spectral distribution that, for large masses, the produced density becomes very sensitive to the value ofγ, while it is not affected by changes inσ.

Total density of produced particles
Once we have discussed the effect of the derivative couplings on the spectral distribution of produced particles, we are ready to focus on the total produced density at a time η > η after reheating We take η = 50/H 0 as a working value. The subsequent evolution ofN (η) is constant for η > η . Therefore, the evolution of the density n(η) will be solely due to the expansion of the background as the gravitational production is already stabilised. Let us start discussing the effects of the derivative couplings in the small mass regime for the field. In figure 4 we show the number density of particles computed numerically at η for the small mass case (m = 0.1, ξ = 1/6) as a function ofσ andγ.
The graphics in the left panel depict the density of produced particles as a function of the derivative couplingσ for several different constant values ofγ. As we can see and have already mentioned in the previous analysis, we are in the nonlinear regime of this parameter except for values around |σ| 0.1. One of the main features of the dependence withσ is that the density attains a minimum value for some negative value of the coupling (actually the specific value varies withγ) and grows significantly for positive and larger negative values, which is not surprising as we are already near the limit (4.8). Although these curves may seem parabolic, they actually have higher order contributions. Another important observation at the light of this figure is that the same density is produced for different pairs (σ,γ). In particular, if we take the production in the case with no derivative couplings, we can see that the same total production can be obtained for different sets of nonvanishing values of the derivative couplings, meaning that their respective contributions compensate each other. Although the production may be same, the spectral distribution is certainly different and this provides a mechanism for distinguishing between the case with vanishing derivative couplings and any other with the same total production.
The right panel of figure 4 shows the density for small masses of the field as a function ofγ for several different constant values ofσ. As the values of the coupling we are considering lie well inside the regime (4.8), it is not surprising that we obtain a linear behaviour of the density with the strength of the couplingγ. The intensity of the linear effect of theγ term increases with increasing absolute value of the coupling constantσ. It is worth noting that for positive values ofγ the production diminishes, which is the opposite effect of just changing the effective mass in the de Sitter phase, giving a strong indication that the curvature effects during the transition are of vital importance for understanding the derivative couplings. For the large mass regime, the contribution of theσ coupling is negligible, as can be deduced from the spectral distribution (figure 3). Hence in figure 5, we have focused only on the nonlinear dependence of the total production (form = 10) on theγ coupling. The production is very sensitive to the value ofγ. It attains a minimum value for a certain strength of the coupling. Again we see that the total production may have the same value for differentγ, although the spectral production is very different. The total particle production is a quadratic function of the derivative couplingγ to very good approximation with mass dependent coefficients.
Finally, let us see how the particle production depends on the mass of the fieldm for given values of the derivative couplings, as shown in figure 6. We have chosen these values so they are representative of the regimes we have studied above. More explicitly, we consider the limiting cases used in the discussions and figures of the spectral production both inγ and inσ. Smaller values ofσ do not provide additional information other than inducing small perturbative changes in the total production. We can see from figure 6 that the two regimes we have discussed in the previous subsections are indeed representative of the behaviour for large and small masses, so we will analyse the figure in these terms.  Let us start with the small mass regime. The effect of negativeσ couplings close to 1 in absolute value on the total production is less prominent than the effect of positive values. Indeed the contribution for positive values ofγ grows rapidly as the mass decreases, giving larger total production for small masses. In the rest of cases the maximum production in the small mass regime is attained for masses of orderm = 1. In this region of small mass, perturbative values ofγ induce linear changes in the production (as also would small values ofσ, although not shown in the figure) that grows in importance as the mass decreases.
The large mass regime presents some nontrivial characteristic as we have already mentioned. To begin with, the purple and blue curves describing the production for perturbative values ofγ quickly grow for large masses. The contribution of positive values ofγ changes from diminishing the total production for intermediate masses to wildly enhancing it after a threshold mass that depend on the specific value of the coupling. We can understand this effect using the previous analysis of the spectral production: as the mass becomes larger, nonlinear effects enter the game. The linear effect for positive values of this coupling is to damp the amplitude of the de Sitter peak. The nonlinearity introduces a broad band of excited modes that grows with the strength of the coupling and, as we have already mentioned, becomes the principal contribution to the production. These effects are generic for any value ofγ. Note that the growth in importance of the γ coupling with larger mass is already suggested by the linear analysis of the mode as the linear contribution grows with the mass (4.10). A similar behaviour could be expected for theσ coupling for the same reason. However, we see a total production that is entirely unaffected by this coupling to the traceless Ricci curvature tensor. This negligible behaviour continues for masses greater than the ones represented. This is in agreement with the study performed in subsection 4.3, where we concluded that theσ contribution to the perturbed mode is very much suppressed (see (4.10)). This is also in consonance with the spectral analysis above. The production is dominated by the de Sitter phase and it is not surprising that the production decays asN ∝ m −1 , given the adiabatic evolution of the background geometry.
The behaviour of the total production in the intermediate mass region is simply an interpolation of the two regimes (small and large mass) already described.
In conclusion, we have seen that including derivative couplings to the background curvature changes in a highly nontrivial way the gravitational production in the early universe for the scalar field even in the perturbative regime.

Conclusions
We have studied the role of derivative couplings to the background curvature in the gravitational production for a quantum scalar field. We have considered a coupling to the curvature scalar mediated through a coupling constant γ and a term proportional to the traceless Ricci tensor through a constant σ.
During the inflationary epoch of the Universe, mimicked by the de Sitter solution, only the term proportional to γ is nonvanishing, so there is no contribution to the production due to different values of σ. We have analytically calculated both the spectral and the total production for the large mass regime, obtaining that the total production decreases with the inverse of the effective mass, which depends on the mass itself and the couplings to the curvature scalar ξ and γ.
After inflation we have considered a reheating phase dominated by the oscillations of the inflaton field around the minimum of its potential. We have interpolated the behaviour of all the background quantities from their values at the end of the de Sitter phase to their values at the onset of the inflaton oscillations using sufficiently smooth functions so the results are independent of the considered order of the interpolation functions. The production in the transition phase is strongly dependent on the value of the couplings γ and σ. To focus only on the relevance of the derivative couplings we have set ξ = 1/6. The gravitational production quickly stabilises during the reheating phase (we are neglecting the oscillations of the background quantities during reheating).
The perturbative regime for the derivative couplings comes from requiring that the derivative terms in the action are perturbative in the sense that they are much smaller than one. However, small values for the coupling constants can induce nonlinear effects because of nontrivial contributions as compared with the field mass. In addition, the validity of the linear approximation for the spectral production requires additional conditions that depend on the mass. The summary is that the perturbative regime requires that the (dimensionless) couplings be much smaller than one, while the linear regime for the production requires conditions that depend on the mass. In this sense, it is not surprising to have nonlinear effects in the production even though the couplings are perturbative, as happens for sufficiently large masses.
The effects of both derivative terms are of different relative importance for a given mass of the field. Moreover, their contributions to the gravitational production depend appreciably on the mass of the field. We have seen that for small masses, values of γ lying within the perturbative regime give contributions to the particle production that are comparable to the ones due to significantly higher values of σ. In this regime, the gravitational production can change notably for perturbative values of the coupling constants so their effects cannot be neglected.
The behaviour of the production drastically changes for large masses. For values of γ still well within the perturbative regime, the effects on the production become nonlinear while the effects of the σ term are negligible. Hence, for large masses the production becomes very sensitive to the value of the derivative coupling γ. This can be understood in terms of the different conditions that γ must satisfy to remain in the linear regime and which are violated for sufficiently large masses. On the other hand note that the σ coupling becomes irrelevant for large masses.
To conclude, we have seen that the presence of derivative couplings to the Riemann tensor can significantly affect the gravitational production. Therefore, they can change the current picture of gravitational production of dark matter [4][5][6][7] as the production changes in a highly nontrivial way even when the coupling constants are in the perturbation regime.
The derivative couplings are often neglected in all the discussions on gravitational production but, for certain theoretical field models, they can be important, and hence their gravitational production will be dominated by the effects we have discussed. This situation will be general for (pseudo-)Nambu-Goldstone bosons, whose disformal couplings dominate their phenomenology. In this context, the presence of disformal scalars, such as branons, are well-motivated from a quantum approach, and its gravitational production will differ from a traditional scalar field with nonderivative couplings.