Particle decay in post inflationary cosmology

We study scalar particle decay during the radiation and matter dominated epochs of a standard cosmological model. An adiabatic approximation is introduced that is valid for degrees of freedom with typical wavelengths much smaller than the particle horizon ($\propto$~Hubble radius) at a given time. We implement a non-perturbative method that includes the cosmological expansion and obtain a cosmological Fermi's Golden Rule that enables one to compute the decay law of a parent particle of mass $m_1$, along with the build up of the population of daughter particles of mass $m_2$. The survival probability of the decaying particle is $P(t)=e^{-\widetilde{\Gamma}_k(t)\,t}$ with $\widetilde{\Gamma}_k(t)$ being an \emph{effective momentum and time dependent decay rate}. It features a transition time scale $t_{nr}$ between the relativistic and non-relativistic regimes and for $k \neq 0$ is always smaller than the analogous rate in Minkowski spacetime, as a consequence of (local) time dilation and the cosmological redshift. For $t \ll t_{nr}$ the decay law is a"stretched exponential"$P(t) = e^{-(t/t^*)^{3/2}}$, whereas for the non-relativistic stage with $t \gg t_{nr}$, we find $P(t) = e^{-\Gamma_0 t}\,(t/t_{nr})^{\Gamma_0\,t_{nr}/2}$. The Hubble time scale $\propto 1/H(t)$ introduces an energy uncertainty $\Delta E \sim H(t)$ which relaxes the constraints of kinematic thresholds. This opens new decay channels into heavier particles for $2\pi E_k(t) H(t) \gg 4m^2_2-m^2_1$, with $E_k(t)$ the (local) comoving energy of the decaying particle. As the expansion proceeds this channel closes and the usual two particle thresholds restrict the decay kinematics.

cosmology use decay rates obtained from S-matrix theory in Minkowski spacetime. In this formulation, the decay rate is obtained from the total transition probability from a state prepared in the infinite past (in) to final states in the infinite future (out). Dividing this probability by the total time elapsed enables one to extract a transition probability per unit time. Energy conservation emerging in the infinite time limit yields kinematic constraints (thresholds) for decay processes.
The decay rate so defined, and calculated, is an input in analyses of cosmological processes.
In an expanding cosmology with a time-dependent gravitational background, there is no global time-like Killing vector; therefore, particle energy is not manifestly conserved, even in spatially flat Friedmann-Robertson-Walker (FRW) cosmologies, which do supply spatial momentum conservation. Early studies of quantum field theory in curved space-time revealed a wealth of unexpected novel phenomena, such as particle production from cosmological expansion [20][21][22][23][24][25][26][27][28][29] and the possibility of processes that are forbidden in Minkowski space time as a consequence of energy/momentum conservation. Pioneering investigations of interacting quantum fields in expanding cosmologies generalized the S-matrix formulation for in-out states in Minkowski spacetimes for model expansion histories. Self-interacting quantized fields were studied with a focus on renormalization aspects and contributions from pair production to the energy momentum tensor [23,24]. The decay of a massive particle into two massless particles conformally coupled to gravity was studied in Ref. [30] within the context of in-out S-matrix for simple cosmological space times. Particle decay in inflationary cosmology (near de Sitter space-time) was studied in Refs. [31,32], revealing surprising phenomena, such as a quantum of a massive field decaying into two (or more) quanta of the same field. The lack of a global, time-like Killing vector, and the concomitant absence of energy conservation, enables such remarkable processes that are forbidden in Minkowski spacetime. More recently, the methods introduced in Ref. [30] were adapted to study the decay of a massive particle into two conformally massless particles in radiation and "stiff" matter dominated cosmology, focusing on extracting a decay rate for zero momentum [33]. The results of Ref. [33] approach those of Minkowski spacetime asymptotically in the long-time limit.
Motivation, goals and summary of results. The importance and wide range of phenomenological consequences of particle decay in cosmology motivate us to study this process within the realm of the standard post inflationary cosmology, during the radiation and matter dominated eras.
Our goal is to obtain and implement a quantum field theory framework that includes consistently the cosmological expansion and that can be applied to the various interactions and fields of the standard model and beyond.
Brief summary of results: We combine a physically motivated adiabatic expansion with a non-perturbative method that is the quantum field theoretical version of the Wigner-Weisskopf theory of atomic line-widths [34] ubiquitous in quantum optics [35]. This method is manifestly unitary, and has been implemented in both Minkowski spacetime and inflationary cosmology [36,37], and provides a systematic framework to obtain the decay law of the parent along with the production probability of the daughter particles. One of our main results, to leading order in this adiabatic expansion, is a cosmological Fermi's Golden Rule wherein the particle horizon (proportional to the Hubble time) determines an uncertainty in the (local) comoving energy. We find that the parent survival probability may be written in terms of an effective time-dependent decay rate which includes the effects of (local) time dilation and cosmological redshift, resulting in a delayed decay. This effective rate depends crucially on a transition time, t nr , between the relativistic and non-relativistic regimes of the parent particle, and is always smaller than that in Minkowski spacetime, becoming equal only in the limit of a parent particle always at rest in the comoving frame.
An unexpected consequence of the cosmological expansion is that the uncertainty implied by the particle horizon opens new decay channels to particles heavier than the parent. As the expansion proceeds this channel closes and the usual kinematic thresholds constrain the phase space for the decay process. While in this study we focus on the radiation dominated (RD) era, our results can be simply extended to the subsequent matter dominated (MD) and dark energy dominated eras.
In appendix (A) we implement the Wigner-Weisskopf method in Minkowski spacetime to provide a basis of comparison which will enable us to highlight the major differences with the cosmological setting.

II. THE STANDARD POST-INFLATIONARY COSMOLOGY
We focus on the decay of particles in the post-inflationary universe, described by a spatially flat (FRW) cosmology with the metric in comoving coordinates given by g µν = diag(1, −a 2 , −a 2 , −a 2 ) . (II.1) The standard cosmology post-inflation is described by three distinct stages: radiation (RD), matter (MD) and dark energy (DE) domination; we model the latter by a cosmological constant.

(II.3)
It is convenient to pass from "comoving time," t, to conformal time η with dη = dt/a(t), in terms of which the metric becomes (a ≡ a(η)) g µν = diag(a 2 , −a 2 , −a 2 , −a 2 ) . In the standard cosmological picture and the majority of the most well-studied variants, most of the interesting particle physics processes occur during the RD era and so we focus most of our attention on this epoch; however, we also contemplate the possibility of long-lived dark matter particles that would decay on time scales of the order of 1/H 0 . The RD and MD epochs cover approximately half of the age of the Universe and during these stages the evolution of the scale factor can be written as which facilitates the explicit analytical study of the decay laws. In turn, the conformal time at a given scale factor a is given by During the (RD) stage the relation between conformal and comoving time is given by 10) a result that will prove useful in the study of the decay law during this stage.

III. THE MODEL:
We consider two interacting scalar fields φ 1 , φ 2 in the FRW cosmology determined by the metric (II.1), with action given by is the Ricci scalar, and ξ 1,2 are couplings to gravity, with ξ = 0, 1/6 corresponding to minimal or conformal coupling, respectively. We identify φ 1 as the field associated with the decaying ("parent") particle, and φ 2 as that of the decay product ("daughter") particles.
Expressing the action of Eq. (III.1) in terms of the comoving spatial coordinates and the conformal time, while rescaling the fields as neglecting surface terms as usual, where For the standard cosmology, using (II.5) Quantization: We begin with the quantization of free fields [23,[25][26][27][28] (λ = 0) as a prelude to the interacting theory. The Heisenberg equations of motion for the conformally rescaled fields in conformal time are It is convenient to consider the spatial Fourier transform in a comoving volume V , namely, for either field, respectively.
Although solutions of (III.9) can be found for separate stages or model expansion histories [33], solving for the exact mode functions for the standard cosmology with the different stages, even when neglecting the term with a ′′ /a, is not feasible. Instead we focus on obtaining approximate solutions in an adiabatic expansion [23, 25-28, 41, 42] that relies on a separation of time scales between those of the particle physics process and that of the cosmological expansion. As an example, let us consider a physically motivated setting wherein the decaying particle has been produced ("born") early during the radiation dominated stage by the decay of heavier particle states at the Grand Unification (GUT) scale ≃ 10 15 GeV. Assuming that the mass of the (DM) particle is much smaller than this scale, the production process will endow the (DM) particle with a physical momentum k p (η) = k/a(η) ≃ 10 15 GeV with k being the comoving momentum. If the environmental temperature of the plasma is T ≃ T GUT ≃ 10 15 GeV and neglecting the processes that reheat the photon bath by entropy injection from massive degrees of freedom, then T GUT ≃ T CMB /a(η i ) implying that the scale factor at the GUT scale a(η i ) ≃ 10 −28 . In turn this estimate implies that the comoving wavevector k with which the (DM) is produced is k ≃ 10 −13 GeV.
The result (III.6) suggests that when considering initial conditions at the GUT scale (or below) corresponding to a(η i ) ≥ 10 −28 the term a ′′ /a in (III.9) can be neglected for ω k (η i ) ≫ 10 −30 GeV for scalar fields minimally coupled to gravity (or for any |ξ j | ≃ O(1)), since ω 2 . This condition is most certainly realized for particles produced from processes at the GUT scale, since as argued above such processes would yield comoving wavectors k ≃ 10 −13 GeV, hence ω k (η i ) ≥ 10 −13 GeV for (DM) particles (or daughters) with masses below the GUT scale. Therefore under these conditions we can safely ignore the term with a ′′ /a in (III.9). Below (see eqn. (III.26) and following comments) we show explicitly that this term is of second order in the adiabatic expansion and can be ignored to leading order. The mode equations (III.9) now become Field quantization is achieved by writing where the mode functions g k (η) obey the equation of motion with the Wronskian condition so that the annihilation a k and creation a † k operators are time independent and obey the canonical commutation relations [a k , a † k ′ ] = δ k, k ′ . Writing the solution of this equation in the WKB form [23,[25][26][27][28] and inserting this ansatz into (III.10) it follows that W k (η) must be a solution of the equation [25] (III. 15) This equation can be solved in an adiabatic expansion We refer to terms that feature n-derivatives of ω k (η) as of n-th adiabatic order. The nature and reliability of the adiabatic expansion is revealed by considering the term of first adiabatic order for generic mass m: this is most easily recognized in comoving time t, introducing the local energy E k (t) and Lorentz factor γ k (t) measured by a comoving observer in terms of the physical momentum k p (t) = k/a(t) and the Hubble expansion rate H(t) =ȧ (t) a(t) = a ′ /a 2 . In terms of these variables, the first order adiabatic ratio (III.17) becomes .
(III. 20) In similar fashion the higher order terms in the adiabatic expansion can be constructed as well: where R(t) is the Ricci scalar (III.2). Consequently, (III.16) takes the form: Consider that the decaying (parent) particle is produced during the radiation dominated stage at the GUT scale with T ≃ 10 15 GeV, with m ≪ T and k p ≃ T corresponding to E k (t) ≃ T and γ k ≫ 1 (ultrarelativistic). With the number of ultrarelativistic degrees of freedom g eff ≃ 100 the expansion rate is and it follows that This analysis clarifies the separation of scales: the Hubble expansion rate H(t) ≪ E k (t), namely there are many oscillations of the field during a Hubble time and the ratio is further suppressed by large local Lorentz factors. This ratio becomes smaller as the scale factor grows and the Hubble rate slows, thereby improving the reliability of the adiabatic expansion. For example, today H(t 0 ) ≃ H 0 ≃ 10 −42 GeV, which is much smaller than the typical particle physics scales even for very light axion-like (DM) candidates.
Therefore we adopt the ratio as the small, dimensionless adiabatic expansion parameter. The physical interpretation of this (small) ratio is clear: typical particle physics degrees of freedom feature wavelengths that are much smaller than the particle horizon proportional to the Hubble radius at any given time (see discussion section below for caveats).
Consequently, when considering the term a ′′ /a in the equation of motion (III.9), we find that Therefore the ratio a ′′ /ω 2 k a is of second adiabatic order and can be safely neglected to the leading adiabatic order which we will pursue in this study, justifying the simplification of the mode equations to (III.10).
In this article we consider the zeroth-adiabatic order with the mode functions given by postponing to future study higher adiabatic corrections (see discussion section below). The phase of the mode function has an immediate interpretation in terms of comoving time and the local comoving energy (III.18), namely which is a natural and straighforward generalization of the phase of positive frequency particle states in Minkowski space-time.

IV. PARTICLE INTERPRETATION: ADIABATIC HAMILTONIAN
Unlike in Minkowski space-time where the full Lorentz group unambiguously leads to a description of particle states associated with Fock states that transform irreducibly and are characterized by mass and spin, the definition of particle states in an expanding cosmology without a global time-like Killing vector is more subtle [20,23,[25][26][27][28].
Our goal is to study particle decay implementing the adiabatic approximation described above, focusing on the leading, zeroth order contribution with the mode functions (III.27). Field quantization in terms of these modes entail that the creation and annihilation operators of the adiabatic particle states depend on time so that the quantum field obeys the (free field) Heisenberg equations of motion. Passing to the interaction picture to obtain the transition amplitudes and probabilities, we would need the explicit time dependence of the creation and annihilation operators. In this section we show explicitly that to leading adiabatic order the operators that create and annihilate the adiabatic states are time independent. This is an important simplification that allows the calculation of matrix elements in a straightforward manner.
In order to establish a clear identification of the zeroth order adiabatic modes with particles we analyze the free-field Hamiltonian, which in terms of the conformally rescaled field operators is given by Writing the field operators in terms of their Fourier expansions, we have Integrating over d 3 x, gathering terms and neglecting the term a ′′ /a in (III.9) as discussed above, we find We can now expand these coefficients Ω k (η) and ∆ k (η) in terms of the functions W k (η) by using the explicit form of the mode functions and using the relation (III.15) the frequencies Ω k (η); ∆ k (η) can be written as It is convenient to introduce which allows us to rewrite the Hamiltonian as This Hamiltonian can be diagonalized by a time-dependent Bogoliubov transformation. We do this in two steps. First we writeã and choose θ k (η) to absorb the phase of ∆ k , i.e., tan θ k (η) = W ′ k (η)/α k (η). Then where We introduce the Bogoliubov transformation to a new set of creation and annihilation operatorŝ noting that u k , v k are real functions of η and | k| only. For theb k ,b † k to obey the canonical commutation relations, it follows that u 2 k − v 2 k = 1. Then the Hamiltonian can be written 15) and the u k and v k chosen to make off-diagonal terms vanish. Then writing u k = cosh φ k and v k = sinh φ k , we find In the second step we absorb the fast phases into the redefinition in terms of which the Hamiltonian can be written as This is a remarkable result: the new operators b † k , b k define a Fock Hilbert space of adiabatic eigenstates, the exact frequencies of which are the zeroth order adiabatic frequencies ω k (η) = k 2 + m 2 a 2 (η). We emphasize that b † k (η), b k (η) depend explicitly on time because the Bogoliubov coefficients u k (η), v k (η) depend on time, while the original operators a k , a † k are time independent in the Heisenberg picture. This is also clear by inverting the relations (IV.13), and using (IV.10) the redefinition (IV.17) along with u 2 (IV.20) Using (III.15) and the adiabatic expansion (III.16) it is straightforward to find that Hence, to zeroth order in the adiabatic expansion b k = a k and the annihilation and creation operators of adiabatic particle states are independent of time. Time dependence of the operators b k , b † k emerges at second order in the adiabatic expansion.
Therefore, the study in this section justifies our identification of particle states to leading (zeroth) order in the adiabatic expansion, namely the time independent operators a † , a create and annihilate zeroth order adiabatic particle states of time dependent frequency ω k (η). This is important because below we cast the interaction picture in terms of these operators and the mode functions g k (η). The analysis above explicitly shows the consistency of this approach to leading order in the adiabatic approximation. In higher order the time evolution of the operators b, b † entail particle production [20, 23, 25-28, 41, 42], an important aspect that will be relegated to future study (see discussion section below). In the analysis that follows we will consider the leading (zeroth) order adiabatic modes.

V. THE INTERACTION PICTURE IN COSMOLOGY
The creation and annihilation operators a k , a † k for each respective field define Fock states, with a vacuum state |0 defined by a k |0 = 0. Since to leading order in the adiabatic approximation a, a † coincide with b, b † associated with single particle adiabatic states, it follows that a † k |0 are identified (to this order) with the single particle states associated with the mode functions(III.27).
In the Schrödinger picture, quantum states obey where in general the Hamiltonian carries explicit η dependence. The solution of (V.1) is given in Consider a Hamiltonian that can be written as is the free theory Hamiltonian and H i (η) the interaction Hamiltonian. In the absence of interactions with H i = 0, the time evolution operator of the free theory U 0 (η, η 0 ) obeys It is convenient to pass to the interaction picture, where the operators evolve with the free field Hamiltonian and the states carry the time dependence from the interaction, namely and their time evolution is given by The unitary time evolution operator in the interaction picture U I (η, η 0 ) obeys For the conformal action (III.4) it follows that where the fields are given by the free field expansion (III.11) with the mode functions solutions of (III.12,III.13) and time independent creation and annihilation operators for the respective fields.
The perturbative solution of eqn. (V.6) is Amplitudes and probabilities in perturbation theory.
Before we consider the non-perturbative Wigner-Weisskopf method, we study the transition amplitudes and probabilities in perturbation theory as this will yield a clear interpretation of the results.
Let us consider the amplitude for the decay process χ 1 → 2 χ 2 given by p , a = 1, 2 are the single particle states associated with the respective fields. With the expansion (V.8) we find to lowest order in perturbation theory, The total transition probability is given by and taking the continuum limit yields Noting the property and introducing the identity Θ(η 2 − η 1 ) + Θ(η 1 − η 2 ) = 1, relabelling terms and using the property (V.14), we find We define the transition rate We emphasize to the reader that in typical S-matrix calculations in Minkowski spacetime, the presence of a time-like Killing vector (and the implementation of the infinite time limit) lead to a time independent transition rate and subsequently a standard exponential decay law. In FRW spacetime, this approach is in general invalid. Rather, the transition rate introduced above will define the decay law obtained within the non-perturbative Wigner-Weisskopf framework described below.

VI. WIGNER-WEISSKOPF THEORY IN COSMOLOGY
The quantum field theoretical Wigner-Weisskopf method has been introduced in refs. [36,37], where the reader is referred to for more details. As discussed in these references, this method is manifestly unitary and leads to a non-perturbative systematic description of transition amplitudes and probabilities directly in real time. Here we describe the main aspects of its implementation within the cosmological setting. Consider an interaction picture state |Ψ(η) I = n C n (η)|n , expanded in the Hilbert space of the free theory; these are the Fock states associated with the annihilation and creation operators a k , a † k of the free field expansion (IV.2) for each field. Inserting into (V.6) yields an exact set of coupled equations for the coefficients In principle this is an infinite hierarchy of integro-differential equations for the coefficients C n (η); progress can be made, however, by considering states connected by the interaction Hamiltonian to a given order in the interaction. Consider that initially the state is |A so that C A (η i ) = 1 ; C κ (η i ) = 0 for |κ = |A , and consider a first order transition process |A → |κ to intermediate multiparticle states |κ with transition matrix elements κ|H I (η)|A . Obviously the state |κ will be connected to other multiparticle states |κ ′ different from |A via H I (η). Hence for example up to second order in the interaction, the state |A → |κ → |κ ′ . Restricting the hierarchy to first order transitions from the initial state |A ↔ |κ results in a coupled set of equations These processes are depicted in fig. (1).

|A |κ |κ
|A κ|H I |A A|H I |κ and inserting this solution into equation (VI.2) we find This integro-differential equation with memory yields a non-perturbative solution for the time evolution of the amplitudes and probabilities. In Minkowski space-time and in frequency space, this is recognized as a Dyson resummation of self-energy diagrams, which upon Fourier transforming back to real time, yields the usual exponential decay law [36,37]. Introducing the solution for C A (η) back into (VI.3) yields the build-up of the population of "daughter" particles.
The equation (VI.5) is in general very difficult to solve, but progress can be made under the weak coupling assumption by invoking the Markovian approximation. A systematic implementation of this approximation begins by introducing with the condition Then (VI.5) can be written as which can be integrated by parts to yield Since E A ∝ O(H 2 I ) the first term on the right hand side is of order H 2 I , whereas the second is Therefore to leading order in the interaction, the evolution equation for the amplitude becomes This expression clearly highlights the non-perturbative nature of the Wigner-Weisskopf approximation. The imaginary part of the self energy Σ A yields a renormalization of the frequencies which we will not pursue here [36,37], whereas the real part gives the decay rate, with (VI.14) Finally, the amplitude for the decay product state |κ is obtained by inserting the amplitude (VI.13) into (VI.4), and the population of the daughter particles is |C κ (η)| 2 .
In our study the state |A is a single particle state of momentum k of the decaying parent particle.

A. Disconnected vacuum diagrams
Before we set out to obtain the self-energy and decay law for a single particle state of the field χ 1 into two particles of the field χ 2 we must clarify the nature of the vacuum diagrams. Consider the initial single particle state |A = |1 (1) k → |1 (2) p ; 1 in which the initial state is annihilated and the two particle state produced, and b): a four particle state in which the initial state evolves unperturbed and a three particle state |1 (1) − p− q is created out of the vacuum by the perturbation. These contributions are depicted in fig. (2). (1) k to first order in H I . Solid lines single particle states of the field χ 1 , dashed lines are single particle states of the field χ 2 .
These processes yield two different contributions to κ 1 The second disconnected diagram (b) corresponds to the "dressing" of the vacuum. This is clearly understood by considering the initial state to be the non-interacting vacuum state |0 ; it is straightforward to repeat the analysis above to obtain the closed set of leading order equations that describe the build-up of the full interacting vacuum state. One finds that diagram (b) without the non-interacting single particle state precisely describes the "dressing" of the vacuum state.
Clearly, similar disconnected diagrams enter the evolution of any initial state. In the case under consideration, namely the decay of single particle states, the disconnected diagram (b) does not contribute to the decay but to the definition of a single particle state obtained out of the full vacuum state. In S-matrix theory these disconnected diagrams are cancelled by dividing all transition  matrix elements by 0|S|0 . Within the Wigner-Weisskopf framework, we write the amplitude for the single particle state |A = |1 where C 0 (η) is the amplitude for the interacting vacuum state obeying (3). With the total self energy given by the sum of the decay (a) and vacuum (b) diagrams as in figure (3), it follows that the amplitudẽ A is determined only by the connected (decay) self energy diagram (a). This is precisely the same as dividing by the vacuum matrix element in S-matrix theory where similar diagrams arise in Minkowski space time with a similar interpretation [36,37]. This is tantamount to redefining the single particle states as built from the full vacuum state. Therefore we neglect diagram (b). We emphasize that this is in contrast with the method proposed in ref. [33] wherein following ref. [30] the disconnected diagram (b) is kept in the calculation of the decay process. Now we are able to calculate the general form of the decay law by considering the decay process χ 1 → 2χ 2 in the interacting theory with H I (η) given by (V.7) to leading order in λ and keeping only the connected diagrams.
The initial state corresponds to single particle state of the χ 1 field |A = |1 (1) k , and the decay process corresponds to a transition to the state |κ = |1 Taking the continuum limit, summing over intermediate states, and accounting for the Bose symmetry factor in the final states yields This is precisely the self-energy (V.13) obtained in the perturbative description of the transition probability and amplitude, equation (V.12), which enters in the evolution equation (VI.5) for the single (parent) particle. Following the steps of the Markovian approximation leading up to the final result (VI.14), we find This expression for the probability makes manifest the non-perturbative nature of the Wigner-Weisskopf method.

VII. DECAY LAW IN LEADING ADIABATIC ORDER.
In this article we study the decay law in the theory described above to leading adiabatic order, namely zeroth order. The study of higher adiabatic order effects, primarily associated with the production of particles by the cosmological expansion, is relegated to a future article (see discussion section below).
In the leading (zeroth) order adiabatic approximation the mode functions are given by , ω k (η ′ ) = k 2 + m 2 a 2 (η ′ ) , (VII.1) and Σ k takes the following form where q = | k− p|. Obviously even to this order both the time and momentum integrals are daunting.
However, progress is made by first considering the case of a massive parent particle decaying into two massless daughter particles. This study will reveal a path forward to the more general case of all massive particles.
A. Massive parent, massless daughters in RD cosmology: Setting m 2 = 0, the self energy becomes The momentum integral in (VII.3) is carried out exactly by introducing a convergence factor after which it becomes Changing integration variables from d(cos(θ)) to q = | k − p| this integral becomes The decay width Γ k (η) is obtained from eqn. (VI.21), and is given by where a factor of 1 2 originates from the integration of the delta function in η (the "prompt" term), and where we set η i = 0 and introduce The expression for S can be simplified substantially, revealing a hierarchy of time scales associated with the adiabatic expansion in radiation domination, during which First we address the integral To begin with we introduce the dimensionless variable (in what follows we suppress the η dependence of z to simplify notation) where E k (t) = k 2 p (t) + m 2 is the physical energy measured locally by a comoving observer with k p (t) = k/a(η) the physical momentum, and H(t) = a ′ (η)/a 2 (η) = 1/(η a(η)) during radiation domination, while H(t) = 2/(η a(η)) during matter domination. The dimensionless ratio (VII. 13) is the inverse of the adiabatic ratio H(t)/E k (t) (we have suppressed the momentum and η dependence in z to simplify notation). The inequality in (VII.13) is a consequence of the adiabatic approximation wherein the physical wavelengths are much smaller than the Hubble radius (∝ the particle horizon). Next, we write η ′′ = η 1 − (η − η ′′ )/η and introduce (VII.14) in terms of which This relation allows us to write where we introduced with the local Lorentz factor given by .
(VII. 18) During (RD) the Lorentz factor can be written as the conformal time η nr determines the time scale at which the parent particle transitions from relativistic η ≪ η nr to non-relativistic η ≫ η nr . In the following analysis we suppress the ηdependence of γ k , z for simplicity.
We emphasize that the relations (VII.15,VII.16) are exact in a radiation dominated cosmology.
Changing integration variables from η ′′ to x given by (VII.14) and using the above variables we find that the integral (VII.12) simplifies to the following expression where δ k (τ ) is of adiabatic order ≥ 1 and given by where we recall that both z and γ k depend explicitly on η. Inserting these results into (VII.8,VII.9,VII.10), and using the new variables z, τ given by eqns. (VII.13,VII.14) we find where the term in the bracket follows from k/ω (1) The expression (VII.23) is amenable to straightforward numerical analysis. However, before we pursue such study, it is important to recognize several features that will yield to a simplification in the general case of massive daughters. The various factors above display a hierarchy of (dimensionless) time scales widely separated by 1/z ≪ 1 in the adiabatic approximation: the "fast" scale τ , the "slow" scale τ /z etc. It is straightforward to find that confirming that δ k is of first and higher adiabatic order. This has a simple, yet illuminating interpretation in terms of an adiabatic expansion of the integral (VII.12). If the frequencies ω (1) k were independent of time, this integral would simply be J k (η, η ′ ) = ω (1) k (η − η ′ ) ≡ τ . Therefore an expansion of J k [η, η ′ ] around η ′ = η would necessarily entail derivatives of ω (1) k , namely terms of higher adiabatic order. Consider such an expansion: In terms of τ = ω (1) where we used (III.20) and (VII.13). The second term is precisely the leading contribution to δ k (VII.26). This analysis makes explicit that an expansion of the integral (VII.12) in powers of η − η ′ is precisely an adiabatic expansion in terms of derivatives of the frequencies. Since the n-th power of η − η ′ in such expansion is multiplied by the n − 1 derivative of the frequencies, and when (η − η ′ ) is replaced by τ /ω (1) k (η) the n − 1 derivative of the frequencies is divided by (ω (1) k (η)) n yielding a dimensionless ratio of adiabatic order n − 1.
Let us now consider the full integral expression for S(η) given by (VII.23) with the corresponding expressions for P [τ ] and δ k (τ ). For z ≫ 1 the terms of the form τ /z, τ 2 /z 2 will be negligible in most of the integration region but for the region of τ ≈ z where these terms become of O(1).
However, because of the factor τ in the denominator of the integrand in (VII.23), a consequence of the momentum integration, this region is suppressed by a factor 1/z ≪ 1 yielding effectively a contribution of first (and higher) adiabatic order. Therefore the contribution from adiabatic corrections, proportional to powers of τ /z are, in fact, subleading. This argument suggests that the zeroth order adiabatic approximation to S(η), namely should be a very good approximation to the full function S(η) for z ≫ 1 with closed form expression where Si[x] is the sine-integral function with asymptotic behavior Si[x] → π/2 − cos(x)/x + · · · as x → ∞. This function rises and begins to oscillate around its asymptotic value at x ≃ π. This behavior implies that the rise-time of Si[A 0 (z; η)] to its asymptotic value in the ultrarelativistic case γ k ≫ 1 increases ∝ γ 2 k . In fact one finds that the full function S(η) and its first order adiabatic approximation S 0 (η) vanish as γ k → ∞. Namely, the contribution from S 0 (and similarly from S) is negligible while the particle is ultrarelativistic. This expectation is verified numerically. Figures (4,5) display S(z) and S(z) − S 0 (z) vs. z for the non-relativistic limit γ k = 1 and for the relativistic regime γ k = 10. In both cases these figures confirm that the zeroth adiabatic approximation S 0 (z) is excellent for z ≫ 1. They also confirm the slow rise of this contribution in the ultrarelativistic case, note the scale on the horizontal axis for the case γ k = 10 compared to that for γ k = 1. For γ k > 1 we have displayed the results for a fixed value of γ k to illustrate the main behavior for the non-relativistic and relativistic limits and highlight that the relativistic case features a larger rise-time. Obviously a detailed numerical study including the η dependence of γ k will depend on the particular values of k, m 1 .
Replacing the function S(η) by the zeroth order approximation S 0 (η) is also consistent with our main approximation of keeping only the zeroth order adiabatic contribution in the mode functions.
Therefore consistently with the zeroth adiabatic order, we find that the decay rate for the case of a massive parent decaying into two massless daughters is given by The decay law of the probability, given by (VI.21) requires the integral of the rate (VII.31). It is now convenient to pass to comoving time in terms of which we find (again setting η i = 0) where Γ 0 is the decay rate of a particle at rest in Minkowski space-time and γ k (t) the time dilation factor, which depends explicitly on time as a consequence of the cosmological redshift of the physical momentum.
Up to the factor F(t ′ ), the decay rate in comoving time has a simple interpretation: namely a decay width at rest divided by the time dilation factor. During (RD) it follows that where t nr (k) is the transition time scale between the ultrarelativistic (t ≪ t nr ) and non-relativistic (t ≫ t nr ) regimes, assuming that the transition occurs during the (RD) era, which is a suitable assumption for masses larger than a few eV.
In the (RD) era we find (using VII.13, VII.18, VII.19, and VII.31) (VII. 37) In Minkowski space time, the calculation of the decay rate in S-matrix theory takes the initial and final states at t = ∓∞ respectively, in which case the Si function attains its asymptotic value and F = 1. The derivation of the decay rate in Minkowski space-time but in real time implementing the Wigner-Weisskopf method is described in detail in appendix (A) and offers a direct comparison between the flat and curved space time results.
In general the integral in (VII.32) must be obtained numerically. However, in order to understand the main differences resulting from the cosmological expansion we first focus on the non-relativistic and the ultra-relativistic limits respectively.
Non-relativistic limit: In this limit we set k = 0 with γ k (t) = 1 and for simplicity we take the Si function to have reached its asymptotic value, therefore replacing F(t ′ ) ≃ 1 inside the integrand yielding 1 This is precisely the decay law in Minkowski space time and coincides with the results obtained in ref. [33]. However this is the case only if the parent particle is "born" at rest in the comoving frame, otherwise the time dilation factor modifies (substantially, see below) the decay rate and law.
Ultra-relativistic limit: In this limit we set m 1 = 0 corresponding to γ k → ∞ in the argument of the Si function, in which case its contribution vanishes identically, yielding F(t ′ ) = 1/2 and with physical wavevector k p (t) = k/a(η(t)). During (RD) this result yields the following decay law of the probability This decay law is a stretched exponential, it is a distinct consequence of time dilation combined with the cosmological redshift of the physical momentum.
Although obtaining the decay law throughout the full range of time entails an intense numerical effort and depends in detail on the various parameters k, m 1 , H R etc. We can obtain an approximate but more clear understanding of the transition between the ultrarelativistic and nonrelativistic regimes by focusing solely on the time integral of the inverse Lorentz factor, because the contribution from F is bound 1/2 ≤ F ≤ 1. Therefore, setting F = 1 and during (RD) we find For the ultrarelativistic regime t ≪ t nr we find the result (VII.40) up to a factor 1/2 because we have set F = 1, whereas in the non-relativistic regime, for t ≫ t nr , we obtain the transition again, the extra power of time is a consequence of the cosmological redshift in the time dilation factor. For k = 0, namely t nr = 0, we obtain the non-relativistic result (VII.38).
The function G k (t) interpolates between the ultrarelativistic case ∝ t 3/2 for t ≪ t nr and the non-relativistic case ∝ t for t ≫ t nr and encodes the time dependence of the time dilation factor through the cosmological redshift.
In Minkowski space time the result of the integral in (VII.41) is simply Γ 0 t which is conveniently written as as Γ 0 t nr (t/t nr ). Because G k is a function of t/t nr , a measure of the delay in the cosmological decay compared to that of Minkowski space time is given by the ratio G k (x)/x with x ≡ t/t nr . This ratio is displayed in fig. (6), it vanishes as x → 0 as x 1/2 and G k (x)/x → 1 as This analysis suggests that the effect of the cosmological expansion can be formally included by defining a time dependent effective decay rate, and t nr depends on k (see (VII.41)), so that the decay law for the survival probability of the parent particle becomes This effective decay rate is always smaller than the Minkowski rate for k = 0 as a consequence of time dilation and its time dependence through the cosmological redshift, coinciding with the Minkowski rate at rest only for k = 0, namely particles born and decaying at rest in the comoving frame.
Finally, the effect of the function F must be studied numerically for a given set of parameters k, m 1 . However, we can obtain an estimate during the (RD) era from the expression (VII.37) for the argument of the Si-function. Writing as x → 0 and approaches π/2 for x ≃ π the large pre-factor in (VII.45) for typical values of k, m 1 entails that the transition between these regimes is fairly sharp, therefore we can approximate the function F(t ′ ) as

B. Massive parent and daughters
We now consider the self energy (VII.2) for the case of massive daughters. Unlike the case of massless daughters, in this case neither the time nor the momentum integrals can be done analytically. However, the study of massless daughters revealed that the adiabatic approximation in the time integrals is excellent when the adiabatic conditions H(t)/E k (t) ≪ 1 are fulfilled for all species. The analysis of the previous section has shown that inside the time integrals we can replace a(η ′ ) → a(η) ; ω k (η ′ ) → ω k (η) since the difference is at least first order (and higher) in the adiabatic approximation (see the results for the factor P (τ ) in eqn. (VII.23)). Furthermore, carrying an adiabatic expansion of the time integrals of the frequencies is tantamount to expanding these in powers of η − η ′ , with the first term, proportional to η − η ′ yielding the zeroth adiabatic order and the higher powers of η − η ′ being of higher adiabatic order. Replacing η − η ′ = τ /ω (1) k (η) associates the higher powers of τ with higher orders in the adiabatic expansion as discussed above. However, this argument by itself does not guarantee the reliability of the adiabatic expansion because for τ ≃ z = E k /H each term in this expansion becomes of the same order. What guarantees the reliability of the adiabatic expansion is the momentum integral that suppresses the large η − η ′ regions. This is manifest in the 1/τ suppression of the integrand in the case of massless daughters (see eqn. (VII.23)). This can be understood from a simple observation. Consider the momentum integral in the full expression (VII.2), setting η = η ′ in the exponent yields a linearly divergent momentum integral. This is the origin of the singularity as η → η ′ in (VII.5). The contributions from regions with large η − η ′ oscillate very fast and are suppressed. Therefore the momentum integral is dominated by the region of small η − η ′ . In appendix (B) we provide an analysis of the first adiabatic correction and confirm both analytically and numerically that it is indeed suppressed by the momentum integration also in the case of massive daughters.
Therefore we consider the leading adiabatic order that yields (VII.47) It is convenient to recast this expression in terms of the physical (comoving) energy and momenta: ω k (η) = a(η)E k (t) absorbing the three powers of a(η) in the denominator in the momentum integral in (VII.47) into the measure d 3 p → d 3 p ph where p ph (η) ≡ p/a(η) is the physical momentum, keeping the same notation for the integration variables (dropping the subscript "ph" from the momenta) to simplify notation, we obtain (VII.48) The variable , (VII.49) corresponds to the physical particle horizon, proportional to the Hubble time. Obviously the momentum integrals cannot be done in closed form, however (VII.48) becomes more familiar with a dispersive representation, namely (VII.50) with the spectral density ρ(k 0 , k; η) = λ 2 a(η) , (VII.51) we refer to (VII.50) the cosmological Fermi's Golden Rule. In the formal limit T → ∞ The density of states (VII.51) is the familiar two body decay phase space in Minkowski space-time for a particle of energy k 0 into two particles of equal mass. It is given by (see appendix (A)), where k ≡ k ph (η) is the the physical momentum, and the theta function describes the reaction threshold.
Before we proceed to the study of Γ k (η) for m 2 = 0, we establish a direct connection with the results of the previous section for m 2 = 0, where the momentum integration was carried out first.

Threshold relaxation:
However, before taking the infinite time limit we recognize important physical consequences of the rate (VII.50). The Hubble time T introduces an uncertainty in energy ∆E ≃ 1/ T ≡ H which allows physical processes that violate local energy conservation on the scale of this uncertainty.
In particular, this uncertainty allows a particle of mass m 1 to decay into heavier particles, as a consequence of the relaxation of the threshold condition via the uncertainty. This remarkable feature can be understood as follows. The sine function in (VII.50) features a maximum at k 0 = E (1) k (η) with half-width (in the variable k 0 ) ≈ π H, narrowing as T increases. The spectral density ρ(k 0 , k; η) has support above the threshold at k * 0 = k 2 + 4m 2 2 , it is convenient to write this threshold as k * 0 = (E (1) k (η)) 2 + (4m 2 2 − m 2 1 ). For 4m 2 2 − m 2 1 < 0 the position of the peak of the sine function, at k 0 = E (1) k (η) lies above the threshold, but for 4m 2 2 − m 2 1 > 0 it lies below it. In this latter case, if the condition is fulfilled, the "wings" of the sine function beyond the peak feature a large overlap with the region of support of the spectral density. This is displayed in figs. (7,8 ). This phenomenon entails the opening of unexpected new channels for a particle to decay into two heavier particles as a consequence of the energy uncertainty determined by the Hubble time. In the adiabatic approximation with E (1) which shows that this condition becomes more easily fulfilled for a relativistic parent. This is clearly displayed in fig. (8).
To gain better understanding of this condition, let us consider the specific case of an ultrarelativistic parent with mass m 1 ≃ 100 GeV with a GUT-scale comoving energy E k ≃ 10 15 GeV decaying into two daughters with mass m 2 ≃ 1 TeV for illustration. We can then replace E k ≃ k/a(η) with k ≃ 10 −13 GeV being the comoving momentum that yields a physical momentum k ph ≃  The range of E, T are chosen to comply with the validity of the adiabatic condition ET ≫ 1.
This figure shows that the uncertainty "opens" the threshold to decaying into heavier particles, the example in the figure corresponds to m 2 = 2m 1 . We have numerically confirmed that as T increases R(E) diminishes as a consequence of a smaller overlap. As the scale factor increases these new decay channels close, allowing only the two body decay for m 1 > 2m 2 and the decay rate is given by the long time limit (VII.55) where Γ 0 is the usual decay rate at rest in Minkowski space time. Following the analysis of the previous section, one now finds a decay law similar to that in eqn. (VII.41) but with Γ 0 now given by (VII.60).

Daughters pair probability:
With the solution for the amplitude of the single particle state, we can now address the amplitude for the decay products from the result (VI.4) with |κ = |1 (2) p , 1 (2) q and |A = |1 (1) k . The decay product is a correlated pair of daughter particles. The corresponding matrix element is given by (VI.19) in terms of the zeroth order adiabatic mode functions (VII.1). Writing the solution for the decaying amplitude where Re E (1) k (η) = Γ k (η)/2, and neglecting the contribution from the imaginary part which amounts to a renormalization of the frequencies [36,37], we find (using VI.4) The time integral is extremely challenging and can only be studied numerically. We can make progress by implementing the same approximations discussed above. Since Γ k depends on the slowly varying frequency, it itself varies slowly, therefore we will consider an interval in η so that the decay rate remains nearly constant, replacing the exponentials by their lowest order expansion in η ′ − η i . During this interval we find the following approximate form of the daughter pair probability, where we set η i = 0. This expression is only valid in restricted time interval, its main merit is that it agrees with the result in Minkowski space time (see appendix A) and describes the early build up of the daughters population from the decay of the parent particle. The occupation number of daughter particles is obtained by calculating the expectation value of the number operators a † q a q ; a † p a p in the time evolved state, it is straightforward to find a † q a q = a † p a p = |C p, k (η)| 2 , (VII.64) the fact that these occupation numbers are the same is a consequence of the pair correlation.
A more detailed assessment of the population build up and asymptotic behavior requires a full numerical study for a range of parameters.

VIII. DISCUSSION
There are several aspects and results of this study that merit further discussion. Medium corrections: In this study we focus on the corrections to the decay law arising solely from the cosmological expansion as a prelude to a more complete treatment of kinetic processes in the early Universe. In this preliminary study we have not included the effect of medium corrections to the interaction vertices or masses. Finite temperature effects, and in particular in the early radiation dominated stage, modify the effective couplings and masses, for example a Yukawa coupling to fermions or a bosonic quartic self interaction would yield finite temperature corrections to the masses ∝ T 2 . These modifications may yield important corrections to the spectral densities and may also modify threshold kinematics. However, the dynamical effects such as threshold relaxation, consequences of uncertainty and delayed decay (relaxation) as a consequence of cosmological redshift of time dilation are robust phenomena that do not depend on these aspects. Our formulation applies to the time evolution of (pure) states. In order to study the time evolution of distribution functions it must be extrapolated to the time evolution of a density matrix, from which one can extract the quantum kinetic equations including the effects of cosmological expansion described here. This program merits a deeper study beyond the scope of this article. We are currently pursuing several of these aspects.
Cosmological particle production: Our study has focused on the zeroth adiabatic order as a prelude to a more comprehensive program. We have argued that at the level of the Hamiltonian, the creation and annihilation operators introduced in the quantization procedure create and destroy particles as identified at leading adiabatic order and diagonalize the Hamiltonian at leading (zero) order. Beyond the leading order, there emerge contributions that describe the creation (and annihilation) of pairs via the cosmological expansion. We have argued that these processes are of higher order in the adiabatic expansion, therefore can be consistently neglected to leading order.
For weak coupling, including these higher order processes of cosmological particle production (and annihilation) in the calculation of the decay rate (and decay law) will result in higher order corrections to the rate of the form λ 2 × (higher order adiabatic). However, once these processes are included at tree level, namely at the level of free field particle production, they may actually compete with the decay process. It is possible that for weak coupling, cosmological particle production (and annihilation) competes on similar time scales with decay, thereby perhaps "replenishing" the population of the decaying particle. The study of these competing effects requires the equivalent of a quantum kinetic description including the gain from particle production and the loss from decay (and absorption of particles into the vacuum). Such study will be the focus of a future report.
Validity of the adiabatic approximation: The adiabatic approximation relies on the ratio . In a radiation dominated cosmology the Hubble radius (H −1 (t)) grows as a 2 (t) and during matter domination it grows as a 3/2 (t) whereas physical wavelengths grow as a(t), with a(t) the scale factor. During these cosmological eras, physical wavelengths become deeper inside the Hubble radius and the ratio H(t)/E k (t) diminishes fast. Therefore if the condition as, sterile neutrinos for example, whose decay may inject energy into the plasma with potential implications for BBN. Such a possibility has been raised in refs. [7]- [14] with regard to the abundance of 7 Li. The decay law of these other species of particles (such as sterile neutrinos beyond the standard model) could be modified and their efficiency for energy injection and potential impact on BBN may be affected by these modifications. Such possibility remains to be studied.
Wave packets: We have studied the decay dynamics from an initial state corresponding to a single particle state with a given comoving wavector. However, it is possible that the decaying parent particle is not created ("born") as a single particle eigenstate of momentum, but in a wave packet superposition. Taking into account this possibility is straightforward within the Wigner-Weisskopf method, and it has been considered in Minkowski space time in ref. [36]. Consider an initial wave packet as a linear superposition of single particle states of the parent field, namely (1) k , where C Caveats: The main approximation invoked in this study, the adiabatic approximation, relies on the physical wavelength of the particle to be deep inside the physical particle horizon at any given time, namely, much smaller than the Hubble radius. If the decaying parent particle is produced ("born") satisfying this condition, this approximation becomes more reliable with cosmological expansion as the Hubble radius grows faster than a physical wavelength during an (RD) or (MD) cosmology. However, it is possible that such particle has been produced during the inflationary, near de Sitter stage, in which case the Hubble radius remains nearly constant and the physical wavelength is stretched beyond it. In this situation, the adiabatic approximation as implemented in this study breaks down. While the physical wavelength remains outside the particle horizon, the evolution must be obtained by solving the equations of motion for the mode function. During the post inflationary evolution well after the physical wavelength of the parent particle re-enters the Hubble radius the adiabatic approximation becomes reliable. However, it is possible that while the physical wavelength is outside the particle horizon during (RD) (or (MD)) the parent particle has decayed substantially with the ensuing growth of the daughter population. The framework developed in this study would need to be modified to include this possibility, again a task beyond the scope and goals of this article.

IX. CONCLUSIONS AND FURTHER QUESTIONS
Motivated by the phenomenological importance of particle decay in cosmology for physics within and beyond the standard model, in this article we initiate a program to provide a systematic framework to obtain the decay law in the standard post inflationary cosmology. Most of the treatments of phenomenological consequences of particle decay in cosmology describe these processes in terms of a decay rate obtained via usual S-matrix theory in Minkowski space time. Instead, recognizing that rapid cosmological expansion may modify this approach with potentially important phenomenological consequences, we study particle decay by combining a physically motivated adiabatic expansion and a non-perturbative quantum field theory method which is an extension of the ubiquitous Wigner-Weisskopf theory of atomic line widths in quantum optics [35]. The adiabatic expansion relies on a wide separation of scales: the typical wavelength of a particle is much smaller than the particle horizon (proportional to the Hubble radius) at any given time. Hence we introduce the adiabatic ratio H(t)/E k (t) where H(t) is the Hubble rate and E k (t) the (local) energy measured by a comoving observer. The validity of the adiabatic approximation relies on H(t)/E k (t) ≪ 1 and is fulfilled under most general circumstances of particle physics processes in cosmology.
The Wigner-Weisskopf framework allows to obtain the survival probability and decay law of a parent particle along with the probability of population build-up for the daughter particles (decay products). We implement this framework within a model quantum field theory to study the generic aspects of particle decay in an expanding cosmology, and compare the results of the cosmological setting with that of Minkowski space time.
One of our main results is a cosmological Fermi's Golden Rule which features an energy uncertainty determined by the particle horizon (∝ 1/H(t)) and yields the time dependent decay rate.
In this study we obtain two main results: i) During the (RD) stage, the survival probability of the decaying (single particle) state may be written in terms of an effective time dependent rate Γ k (t) as P (t) = e − Γ k (t) t . The effective rate is characterized by a time scale t nr (VII.41) at which the particle transitions from the relativistic regime (t ≪ t nr ) when P (t) = e −(t/t * ) 3/2 to the nonrelativistic regime (t ≫ t nr ) when P (t) = e −Γ 0 t t tnr Γ 0 tnr /2 where Γ 0 is the Minkowski space-time decay width at rest. Generically the decay is slower in an expanding cosmology than in Minkowski space time. Only for a particle that has been produced ("born") at rest in the comoving frame is the decay law asymptotically the same as in Minkowski space-time. Physically the reason for the delayed decay is that for non-vanishing momentum the decay rate features the (local) time dilation factor, and in an expanding cosmology the (local) Lorentz factor depends on time through the cosmological redshift. Therefore lighter particles that are produced with a large Lorentz factor decay with an effective longer lifetime. ii) The second, unexpected result of our study is a relaxation of thresholds as a consequence of the energy uncertainty determined by the particle horizon.
A distinct consequence of this uncertainty is the opening of new decay channels to decay products that are heavier than the parent particle. Under the validity of the adiabatic approximation, this possibility is available when 2πE k (t)H(t) ≫ 4m 2 2 − m 2 1 where m 1 , m 2 are the masses of the parent, daughter particles respectively. As the expansion proceeds this channel closes and the usual kinematic threshold constrains the phase space available for decay. Both these results may have important phenomenological consequences in baryogenesis, leptogenesis, and dark matter abundance and constraints which remain to be studied further.

Further questions:
We have focused our study on a simple quantum field theory model that is not directly related to the standard model of particle physics or beyond. Yet, the results have a compelling and simple physical interpretation that is likely to transcend the particular model. However, the analysis of this study must be applied to other fields in particular fermionic degrees of freedom and vector bosons. Both present new and different technical challenges primarily from their couplings to gravity which will determine not only the scale factor dependence of vertices but also the nature of the mode functions (spinors in particular). As mentioned above, cosmological particle production is not included to leading order in the adiabatic approximation but must be consistently included beyond leading order. The results of this study point to interesting avenues to pursue further: in particular the relaxation of kinematic thresholds from the cosmological uncertainty opens the possibility for unexpected phenomena and possible modifications to processes, such as inverse decays, the dynamics of thermalization and detailed balance. These are all issues that merit a deeper study, and we expect to report on some of them currently in progress. therefore the approach to asymptotia and to the full width takes a much longer time for an ultrarelativistic particle with t asy ∝ 2k/m 2 , whereas it is much shorter in the non-relativistic case t asy ∝ 1/m. In S-matrix theory in Minkowski space time one takes t → ∞, and obviously in this limit the Si− function reaches its asymptotic value, therefore the time dependence of the rate cannot be gleaned.
Integrating in time first: massive particles and Fermi's Golden rule.
Let us consider now the full dispersion relations for the daughter particles, calling E k that of the parent decaying particle and ω p = p 2 + m 2 2 that of the daughter. From (VI.7) and (VI.21), we need (A.7) We find ; q = | k − p| , (A.8) the asymptotic long time limit this is simply Fermi's Golden rule which yields the standard result for the decay rate Although E 2 k − k 2 = m 2 1 we have left the result in the form shown to make use of it in the cosmological case and to highlight the threshold.
Before taking the limit t → ∞ the real time rate (A.8) can be conveniently written in a dispersive form, namely with the spectral density as a function of z in this case, remaining perturbatively small when compared to I 0 . Therefore this study confirms that the first order adiabatic corrections are indeed subleading as compared to the leading (zeroth) order contribution for large z = E k (t)/H(t).  Figure 11: The integrals I 1 (z), I 2 (z) vs. z, for m 2 /m 1 = 0.25 , k = 0.