Fermionic decay of a massive scalar in the early Universe

We derive a curved space generalization of a scalar to fermion decay rate with a Yukawa coupling in expanding Friedmann-Robertson-Walker universes. This is done using the full theory of quantum fields in curved spacetime and the added-up transition probability method. It is found that in an expanding universe the usual Minkowskian decay rates are considerably modified for early times. For conformally coupled scalars the decay rate is modified by a positive additive term proportional to the inverse of mass and related to the expansion rate of the Universe. We compare and contrast our results with previous studies on scalar to scalar decay and find that in general the decay channel into fermions is the dominant channel of decay in the very early Universe.


I. INTRODUCTION
Particle decay processes in the early Universe have deep and profound implications in cosmology, from baryogenesis [1,2] to Big Bang nucleosynthesis [3] and reheating scenarios after inflation [4]. When considering these early Universe processes in depth, the influence of spacetime curvature cannot be neglected anymore. In this regime the Minkowskian quantum field theory is ultimately only an approximation and of limited applicability. Instead, when studying particle processes in the early Universe, or in any kind of situation where the effects of gravitation cannot be neglected, quantum field theory in curved spacetime must be used. As a result of spacetime curvature and the related violation of conservation laws, the particle decay rates, cross sections and lifetimes are modified compared to usual flat space results [5,6]. New particle processes, forbidden in Minkowski space, are to be considered leading to new Feynman diagrams even at first-order [5,7,8]. The lack of energy conservation, e.g., in Friedmann-Robertson-Walker (FRW) universes, leads to processes where a particle can even decay into quanta of its own field [9,10]. All these suggest that considerable care must be exercised when studying mutually interacting fields in curved spacetime.
The groundwork for interacting quantum fields in expanding universes was laid in the 1980s expanding the familiar S matrix formulation of in-out states to simple cosmological models. Studies focused on renormalization of self-interacting fields [11] but also on massive scalar particles decaying into two conformally coupled scalar particles [5,[12][13][14]. The method known as added-up probability, introduced in [5], was used to study photon decay in radiation dominated universe [15] and more recently in investigating scalar particle decay in expanding FRW universes [6,16,17]. These investigations into scalar decay revealed that the Minkowskian decay rates are significantly modified at early times, but in the longtime limit approach the Minkowskian results asymptot-ically. Lately, a nonperturbative method, the so called Wigner-Weisskopf method [18] familiar from quantum optics, was adopted to study particle decay in the context of inflationary cosmology [19] as well as post inflationary cosmology [20] and to study decay in renormalizable field theories [21].
Although these studies have to a great degree increased our knowledge of particle decay in curved spacetime, they have mainly been concerned with a scalar channel decay, see e.g., [5,6,12,16,17]. The negligence of the fermionic channel may be attributed to the more complex nature of the problem involving spinors which may lead to some mathematical difficulties encountered in the technical calculations. Clearly, the study of the fermionic decay channel would further increase and complete our knowledge of particle decay in curved space. This serves as the motivation for this paper where we will consider a scalar field interacting with a Dirac spinor through a Yukawa interaction in curved space. By applying the added-up method, we calculate the decay rate for a fermionic decay of a massive scalar in an expanding spatially flat FRW universe. The scalar field is taken with an arbitrary coupling ξ to gravity hence avoiding the restrictions commonly encountered in literature where ξ is usually restricted to be that of minimally or conformally coupled. The fermions are assumed to be massless in order to to avoid interpretational issues of the decay rate when using the added-up method.
This article is structured in the following way. In Sec. II we introduce the necessary theoretical background and give a brief review of the method of added-up probability. In Sec. III we obtain the positive mode solutions for both the spinor field and the scalar field. These modes are used in Sec. IV to derive a transition amplitude in curved space for the scalar to decay into a fermion-antifermion pair. In this section we will also show that the transition amplitude has the correct Minkowskian limit. After that we will compare the decay into fermions with the recently obtained decay into massless scalars in Sec. V where it is shown that the fermionic decay channel is dominant to the scalar channel in the early Universe. Finally, we will discuss the implications and reasons behind our results in Sec. VI and end with conclusions in Sec. VII. arXiv:1904.05084v2 [gr-qc] 23 Apr 2019

II. PRELIMINARIES
As was stated in the introduction, there are at least two distinct ways to calculate the decay probability in curved spacetime. The difference between these two methods lies in the quantity which is eventually being calculated. The added-up method is fundamentally a global expression for the decay rate, while the Wigner-Weisskopf method is concerned more with the differential decay rate. In this article we will use the added-up method introduced in [5] and the purpose of this section is to summarize this method. We begin by introducing the necessary theoretical background and then give a brief review of the method of added-up probability. Natural units = c = 1 are used and the metric is chosen with a positive time component.

A. Theoretical background
We will consider a model described by a fourdimensional spatially flat Friedmann-Robertson-Walker spacetime with the metric given in standard coordinate time t with a dimensionless scale factor a(t). On this classical curved background we consider propagating quantum fields mutually interacting with each other. In particular we consider a Yukawa type interaction where a massive scalar particle decays into two massless fermions. The Lagrangian density L for this theory is given by which consists of a Lagrangian density for the complex scalar field L φ , a fermion Lagrangian L ψ and an interaction term L I . For a complex scalar field in curved space the Lagrangian density is straightforwardly given as where g is the determinant of the metric, R the Ricci scalar and the coupling of the scalar field to gravity is controlled by the dimensionless parameter ξ. In four dimensions, the value ξ = 1/6 is known as conformal coupling, while the minimal coupling is given by ξ = 0. Denoting the covariant d'Alembert operator by , the Klein-Gordon equation in curved spacetime for the scalar field φ is given by The fermion part of the Lagrangian requires a little more thought in curved space. Because the spinor does not transform like a tensor under Lorentz transformation, one cannot just replace the derivatives with their covariant counterparts. To incorporate spinors into general relativity, one can introduce a set of four covariant vector fields e a µ known as a tetrad. We adopt a convention where the latin indices refer to local inertial coordinates while greek indices refer to general coordinates. With this formalism, the Lagrangian for a massless spinor field in curved spacetime is given by where ψ denotes the Dirac conjugate spinor ψ = γ 0 ψ † . The curved space gamma matrices are defined via the tetrad as γ µ = e a µ γ a , where γ a denotes the usual flat space gamma matrix. The curved space gamma matrices satisfy the usual anticommutation relations {γ µ , γ ν } = 2g µν (6) and the covariant derivative is defined with the help of a spin-connection Γ µ as where By varying the action one obtains the Dirac equation for a massless ψ-particle in curved space The choice of a massless spinor field is motivated by its necessity in the added-up method, but this choice can also be thought from a perspective of symmetry. If one considers a complex scalar field with non-zero charge, a chiral spinor field ψ L with opposite charge and a zerocharge field ψ R , in a U (1)-symmetric theory a mass term in the fermion Lagrangian density violates this symmetry while the scalar Lagrangian density is left unchanged. It should also be noted that in this case no mass appears in the renormalization procedure either. Therefore, if one wants to consider a (globally) U (1)-symmetric theory with chiral fields, then the fermionic Lagrangian density should be taken as massless. Moreover, if the field φ is considered as an SU (2) doublet, by defining the (global) U (1)-charges of the theory in a suitable way, the fermionic fields may be thought of as fermions of the Standard Model.
Because a complex scalar field can always be decomposed into two real scalar fields as φ = (φ 1 + iφ 2 )/ √ 2, from now on we will consider only real scalar fields. This is merely to make the technical calculations and theory more manageable. In what follows, the only effect of this onto the calculations is that the total transition amplitude is to be multiplied by a factor of two in the complex case. For a Yukawa type interaction the interaction part of the Lagrangian when considering real scalar fields is of the form where h = 0 is a dimensionless coupling constant chosen to be real. To describe the interaction, we will apply the S matrix scheme where the S matrix is given as whereT denotes the time-ordering operator. The exponential factor e −βt acts as a switch-off for the interaction for large times with β being a positive constant and called the switch-off parameter. The perturbative expansion of the S matrix for this interaction gives where Furthermore, we consider only tree level processes for which the transition amplitude is defined as B. The added-up probability Recall a result concerning free fields in curved spacetime: conformally coupled massless particles are not created as a result of spacetime expansion [22,23]. This result, which applies both for spin-0 as well as spin-1/2 particles, is at the heart of the added-up probability method introduced in [5]. Because of this result, a detection of a massless conformally coupled particle in the asymptotic region of spacetime indicates that this particle must have been solely created or influenced by a decay process.
In curved spacetime there is still one other complication arising from the lack of energy conservation in curved space; the creation of all three particles with a priori unrestricted momenta. This lack of kinematic thresholds was also considered in [18] using the Wigner-Weisskopf method. Because of this, it is necessary in the added-up method to further restrict the out momenta to fulfill the three-momentum conservation law p = k 1 + k 2 . Only after that one is left with what resembles closest a decay process in curved space. This process can be described at tree level with two Feynman diagrams (Fig.1). In the original treatment of Audretsch and Spangehl only decay into scalar particles was concerned. Nevertheless, the formalism used in the added-up method applies equally well for other types of decay product particles as long as they are massless and conformally coupled to facilitate a meaningful decay rate interpretation in curved space. Especially this also holds for the case where a particle decays into two conformally coupled massless fermions. For a more exhaustive treatment of the addedup method, we refer the reader to the original paper [5].
The added-up transition amplitude corresponding to the two Feynman diagrams in Fig. 1 reads as where k 1 = k and k 2 = p − k and |out refers to the state in the out region, i.e., the asymptotic future. In the formalism one could use both the in or out regions for the one-particle and vacuum states but not the usual in and out states used in flat space. This is because such a matrix element suffers from physical interpretational issues in curved spacetime [5]. With this in mind we choose the asymptotic out region for the simple reason that the massive field modes can be recognized in the out region. The total transition probability w tot is obtained by summing over all the k modes, With this brief review of the added-up method we turn our attention to finding the positive field modes needed in the calculations.

III. THE POSITIVE MODE SOLUTIONS
One of the most important issues in quantum field theory in curved spacetime is finding analytic solutions to the field equations. This is a task of great complexity because in general one is lead to solving highly nontrivial differential equations. Even when the solution is obtained, it is not all that certain that the positive mode solutions can be obtained. In any case, one must fix the scale factor in order to even begin finding the field modes.
The literature is filled with numerous studies for finding scalar field modes for different types of universes, see e.g., [24,25] and references therein. Many models studied are non-physical with the purpose of merely illustrating the effects of particle creation from the vacuum. Even if a positive field mode can be found for a more general realistic model, it usually suffers from quite strict restrictions [6]. The fermionic case presents no simplifications into finding the positive modes. Some exact solutions have been reported [26][27][28] and also more recently by Moradi [29,30]. The goal of this section is to construct the positive field modes for the spinor and scalar fields used in the calculation of the transition probability in the next section.

A. Spinor modes
As noted, finding a solution to the curved space Dirac equation is a highly technical task and generally one needs to fix the scale factor explicitly in order to obtain an exact solution. This is the case even when considering a massless Dirac equation like Eq. (9). In a spatially flat FRW universe, the most general scale factor we can choose in this paper is a typical power-law expansion where b is a positive constant controlling the expansion rate of the Universe. We take n ∈ [0, 1) which covers most of the interesting cosmological cases like radiation and matter dominated universes. The de Sitter solution is excluded in this paper because a positive mode solution for the massive scalar field cannot be found for de Sitter space [6]. Let us begin by expanding the curved space Dirac equation (9) in flat FRW spacetime. With the substitution of ψ = a −3/2ψ we can further simplify the ensuing equation to obtain Here k i denotes the vector component of the momentum vector k and γ a are given in chiral representation as where σ 0 is the identity matrix and σ i are the Pauli matrices. The advantage of using chiral representation is that Eq. (18) reduces upon solving e.g., the spinor components ψ 2 , ψ 4 into sets of two second-order differential equationsψ for the components ψ 1 and ψ 3 . This differential equation can readily be solved and upon substituting the answers to the remaining equations for ψ 2 and ψ 4 obtained from Eq. (18) yields the unnormalized solutions for the spinor components. From the solutions we obtain a set of mode solutions for the massless fermion in chiral representation, where the positive and negative energy spinors are to be normalized according to Applying this normalization condition the normalized modes are given as where we have defined k := |k| and k ± := k 1 ± ik 2 . The fields can thus be expanded as

B. Scalar modes
For the scalar field modes, one needs to solve the curved space Klein-Gordon equation (4) and identify the positive modes. The scalar field φ can be expanded in the usual way as and because of the homogeneity of the spatial sections, the positive mode solutions are separable into time and space parts where p := |p|. The positive mode is found by solving the ensuing differential equation for χ p (t) and recognizing the positive frequency modes in the usual way [24]. The normalized solutions χ p (t) were obtained in [6] for a massive scalar field mode with an arbitrary coupling to gravity in the rest frame of the particle p = 0. The rest frame was used in order to get an exact solution for the decay rate in [6] as well as in previous studies [16,17]. In what follows, it turns out that in order to find an exact solution for the transition amplitude, we will also in this paper need to restrict to the rest frame of the decaying particle. Hence the solution obtained in the rest frame is sufficient for our purposes. The solution for the massive scalar mode from [6] reads as where η denotes the conformal time and index α is defined as At this point we have to reconcile the issue that different scale factors with different time variables were used for solving the Dirac equation above and what was obtained for the solution of the scalar particle. Since we are using a scale factor a(t) = bt n and in [6] the scale factor was given in conformal time as a(η) = b η n /2 , we must find the relation between the primed and unprimed variables. The coordinate time can be transformed into coordinate time using the usual relation dt = a(η)dη.
Comparing the scale factors in coordinate time gives the relations between the primed and unprimed variables, In terms of the standard coordinate time t and the unprimed variables, the scalar field mode reads as and the index α as This completes the necessary background formalism.

IV. TRANSITION AMPLITUDE IN CURVED SPACE
Having established the proper formalism, we turn to the calculation of the total transition amplitude and the decay rate. In this section we will explicitly calculate the transition amplitudes corresponding to the two Feynman diagrams of Fig. 1, verify the transition amplitude has a correct form in the Minkowskian limit and discuss about convergence issues near the spacetime singularity. We will end this section with a calculation of the exact decay rate for conformally coupled massive scalars.

A. Total transition amplitude
The matrix element corresponding to the diagram (a) can be written with the positive field modes as For diagram (b), the matrix element is obtained by substituting the scalar mode with its complex conjugate w p (t, x) → w p (t, x) * . Upon inserting the positive field modes given in Sec. II, this integral can further be simplified by taking the limit inside the integral by introducing a time cutoff at t = T . Moreover, the spatial integration can readily be done resulting in taking out the delta function δ(p − k 1 − k 2 ) for diagram (a) or δ(p + k 1 + k 2 ) in case of diagram (b). These correspond to the usual conservation of three-momentum. Hence, in the end, after also performing the k 2 integration, one is left with the amplitude squared corresponding to diagram (a) as and similarly for diagram (b) with the substitution of To proceed ,we switch to the rest frame of the decaying particle where, p = 0. The reasons for this are twofold. First, we only know the massive scalar field mode for a general power-law expansion in its rest frame. Second, by taking p = 0 the integrand simplifies to a form where we can use the properties of distributions to obtain an exact result. We begin by calculating the spin sum in the rest frame giving where k = |k|. Because the spin sum is proportional to the square of the momentum vector length, we can combine the k-integration of both diagrams into by a change of variables k → −k in amplitude squared of diagram (b). Upon inserting the scalar field mode (32) we obtain the following integral The above integral can be written as a three dimensional integral which allows the use of the distribution identity on the k-integration, which reduces Eq. (38) into a two-dimensional integral. The remaining integrations can be performed in a symmetric fashion using derivative identities of distributions 1 . Using also the fact that δ (y − x) = −δ (x − y) we can switch the order of integration to obtain an integral which has two derivative terms. In the end, by performing the remaining delta integration, we arrive at the exact form for the total transition probability as where we have changed to a dimensionless variable n) . Because of the absolute value squared, the realness of Eq. (40) is manifest. Although derived using the real scalar field modes, the total transition probability using the complex field is just Eq. (40) multiplied by a factor of two.
The equation (40) is our main result of this paper. It is an exact formula for the transition probability from which the decay rate may be obtained. The integration can be performed exactly resulting in an extremely long and complicated expression involving hypergeometric and Bessel functions. This expression is of no illuminating value and we do not present it here.
In its current form, Eq. (40) gives the mean decay rate when divided by the time t. Of interest besides this is also on the form of a differential decay rate which can readily be obtained as the integrand of the transition probability. If we make this interpretation the differential decay rate is given as (41)

B. Minkowskian limit
As with any theory in physics, a generalized theory should reproduce the known results in the more specialized situation. Therefore a decay rate in an expanding spacetime should yield the flat spacetime results when the appropriate limit is taken. This limit is given when the spacetime is static, i.e., the scale factor a(t) = 1. This is achieved for the value n = 0 and, although not needed in the final equation (40), one may take b = 1 since it is only a scaling factor. With the value n = 0, the parameter α = 1/2 and the transition probability reads as Because H 1/2 (u) = i 2/(πu)e −iu , we see immediately that the integrand in (42) is equal to 2/π and the transition amplitude reduces to From this, the decay rate is obtained in the usual way by dividing the transiton probability by the infinite time t. Thereby dividing equation (43) by ∆t := t − t 0 , we obtain the Minkowskian decay rate which corresponds to the decay rate calculated in the added-up method using only Minkowskian plane waves. This consistency serves as a validation of our procedure.

C. Convergence near spacetime singularity
Equation (40) gives the transition probability beginning from some time t 0 up to time t, assuming that the total transition probability w tot 1. Although e.g., in post-inflationary scenarios the time t 0 truly differs from zero, in some instances it may be necessary to take the limit t 0 → 0. Whether for mathematical reasons to simplify the integration or to investigate the decay process all the way from singularity. It is therefore illuminative to study also the zero time limit of the transition probability.
The convergence of the integral of Eq. (40) in the limit where t 0 → 0 can be investigated in a simple fashion. First, the Hankel function of the second kind has the following limiting form when the argument approaches zero: H (2) α (u) ∼ iΓ(α)(u/2) −α /π. Second, requiring that the behavior of the power of the variable u is greater than −1, to avoid singular behaviour, we can restrict the values of α for which we obtain a finite result. The condition for which the transition probability is finite is given by the requirement α < 0. This implies that the integral diverges for all values of n and ξ when t 0 → 0 is taken.
There is one value, namely α = (1 − n)/2, however for which the limit is finite. The reason for this is that for this value, the power of the limiting form of the Hankel function exactly cancels out the power of the prefactor leaving only u n to be integrated. Hence, the transition probability Eq. (40) converges near the spacetime singularity when the following condition holds, This condition is valid for Minkowski limit n = 0 as well as for radiation dominated universe n = 1/2 for any value of the coupling. But the most peculiar feature of it is that the equality holds for conformal and only conformal coupling ξ = 1/6 and for any n ∈ (0, 1). Now that we know the integral of Eq. (40) is convergent for conformal coupling at the spacetime singularity, we can calculate it explicitly.

D. Conformally coupled massive particles
As noted, the transition probability (40) is finite for all n ∈ [0, 1) in the limit t 0 → 0 only for conformally cou-pled massive particles. For a conformally coupled massive scalar particle, the decay probability reads as where, to get to the second line, we have used the derivative identity [u α H α−1 . This integral can be integrated exactly to yield where the last constant term arises from the lower limit of the integration. For further insight, we take a look at the asymptotic behaviour of this transition amplitude when the time is large. Given by the leading terms in the asymptotic expansion, the transition probability is given by We see that it is linear in time t, which would prompt us to divide the transition probability by the time t to obtain the decay rate in the familiar fashion of Minkowski space field theory. A complication however arises from the constant term reflecting the fact that in curved space one needs to be more careful concerning usual Minkowski space procedures. To deal with the constant term we adopt the same procedure as was introduced in [5] where this constant term is divided by a time t grav , known as time of gravitational influence. This can be defined as t grav := t f −t i , where t i indicates the time when the gravitational field begins its influence and t f its end. With this definition the asymptotic decay rate is given by It is now seen that asymptotically the decay rate obtains a gravitational correction term to the Minkowskian decay rate. This term is positive implying, in the case of conformal coupled particles, that the effect of the gravitational field is to enhance the decay rate. With Eq. (49) we are now in a position to compare this to the process where a massive scalar decays into two massless scalars.

V. COMPARISON WITH SCALAR TO SCALAR DECAY
In this section we will mainly be concerned in comparing the results of this paper into studies where scalar to scalar decay was concerned using also the added-up method [6,16,17]. This allows a direct comparison due to the same methods of calculations being used.
The asymptotic decay rate for a massive scalar φ decaying into two conformally coupled massless scalars χ was given in [6] as where λ is the coupling constant. This expression was also derived in the limit of t 0 → 0 making the comparison valid. Comparing the expressions (49) and (50) a major difference is immediately seen corresponding to the sign of the relative correction term. For a decay into scalars, the sign is negative, implying that the effect of a gravitational field is to diminish the rate of decay. On the contrary, in fermionic decay the sign is positive and the effect of the gravitational field is to enhance the decay. Both relative correction terms are increasing functions of the parameter n on the interval n ∈ [0, 1), so the gravitational correction term increases in magnitude as n increases. This means that the decay into fermions increases and decay into scalars diminishes as the parameter n increases.
Although the ordering of these decay rates is most prominently seen in the asymptotic decay rates, it cannot be inferred from these that the decay rates are so ordered for all time. To investigate the small time behaviour, one needs to take a look at the exact decay rates. This is done in Fig. 2 for three cosmologically interesting situations, namely universes dominated by stiff matter, radiation and ordinary matter. From Fig. 2 we see that this ordering of the decay rates is the same for these universes for all time and not just asymptotically. Finally, it should be noted that in practical situations t grav is usually much longer than inverse of mass, when these corrections are in practice small. These relative correction terms might not, however, be neglegted altogether, because for t ∼ m the exact equation (42) divided by t should be used instead.
As was noted, the integral in Eq. (40) diverges when t 0 approaches the spacetime singularity. When looking at the behavior of the differential decay rate (41), interpreted as the integrand of (42), we notice that it behaves like t −2α−1 in the vicinity of the spacetime singularity. On the other hand, the scalar to scalar differential rate obtained in [6], behaves like t −2α+1 near the singularity. These observations together imply that near the singularity, as long as 2α + 1 > 2α − 1, the decay channel into fermions is dominant over the decay channel into scalars. However, this inequality is always true implying that in the very early Universe, the decay channel into fermions is always the dominant one to the scalar channel no matter what is the scalar coupling or universe matter content. It must be stressed, however, that this holds only near the singularity and the exact time when it is true is not evident from this analysis because the decay rates are differently proportional to the couplings and masses in the prefactors which affect the decay rate.

VI. DISCUSSION
Our study of scalar particle decay in the early Universe has presented us with novel features concerning the effect of gravitation on the decay process. Not only does the curved background modify the Minkowskian results, but the effect is in some cases opposite for fermionic and scalar channels.
We wish to emphasize three key elements of this paper: (a) We have derived a curved space generalization for a massive scalar to decay into two massless fermions.
(b) For conformally coupled scalars, the decay into fermionic channel is enhanced and decay into scalar channel diminished by the curved background.
(c) It is found that in the very early Universe, the fermionic decay channel is dominant to the scalar channel.
These two last points bear some further elaboration. For conformally coupled scalars, the integration from the singularity introduces an additive constant to the decay rate. In the paper [16] a proposition was put forward in the case of a scalar to scalar decay that the faster the Universe is expanding, the smaller is the decay rate into scalars. We may elaborate this point by speculating within the same framework, i.e., by taking a look at the Hubble parameter H =ȧ/a = n/t. The Hubble parameter shows us that as n increases, the Universe expands relatively faster. In the case of a fermionic decay channel, the faster the Universe is expanding the larger is the decay rate. We may speculate the reason behind this from a statistical point of view. As the Universe expands faster, more and more states are becoming available for the fermions to occupy. For bosons, on the other hand, this same expansion reduces the Bose enhancement thereby diminishing the decay. This would statistically explain the observed phenomena, but further investigations into this feature are surely needed.
As was noted, independent of the matter content of the Universe or the coupling of the massive scalar, the fermionic decay channel is the dominant one in the very early Universe. How close this time should be to the singularity for the fermionic channel to be the dominant one cannot be inferred from the analysis made in Sec. IV. For times farther away from the singularity, it cannot be said that for any fixed n the decay into fermions is more probable because the prefactors on the decay rates differ; the scalar to scalar is proportional to m −1 while scalar to fermion is proportional to the mass m with different couplings. Furthermore our study has been restricted to massless fermions in the fermionic decay channel. The main reason for this has been the fact that in order to apply the added-up method and to obtain a meaningful decay rate interpretation in curved space this needs to be done. For fermions this restriction is, however, less severe when considering times before electroweak phase transition because then fermions are usually massless as in the Standard Model. If the time when the fermion channel is dominant over the scalar channel is before this phase transition, then it could be inferred that the decay into Standard Model fermions is more preferred than massless scalars. Which brings up another point.
The analysis for scalar to scalar decay in [16] was done using massless scalars as decay products. The inference that the fermion decay channel is dominant to the scalar channel only applies when the product scalar particles are massless, although one could argue that for very light scalars this is still valid. The case where the decay products are also massive cannot be straightforwardly calculated using the added-up method and it cannot be said if the fermionic channel is indeed the dominant one in all cases. The decay into massive scalars has been considered in [20] using the methods introduced in [18,19] for post inflationary cosmology where the authors obtained similar type of results indicating slower decay for scalars. It would be a worthwhile study to see whether similar types of results are obtained for a fermionic channel as well.

VII. CONCLUSIONS
Particle decay rates in the early Universe play an ever increasing role in modern cosmology. In this paper we have used quantum field theory in curved spacetime and the added-up probability method to calculate the decay rate of a φ → ψψ process in curved spacetime. We have found that the Minkowskian decay rates are modified in the early Universe by the presence of a gravitational field and that the fermionic decay channel of a scalar is dominant to a scalar decay channel in the very early Universe.
This modification of decay rates should be taken into account when considering with more precision early universe processes involving scalar decay scenarios. For fermionic decay of a scalar, possible scenarios involve e.g., fermionic preheating scenarios. Also, as cosmological data and measurements become increasingly more accurate, it may be necessary in the future to include the effects of curved space also in particle decay rates.