Cosmological decay of Higgs-like scalars into a fermion channel

We study the decay of a Higgs-like scalar Yukawa coupled to massless fermions in post-inflationary cosmology, combining a non-perturbative method with an adiabatic expansion. The renormalized survival probability $\mathcal{P}_\Phi(t)$ of a (quasi) particle ``born'' at time $t_b$ and decaying at rest in the comoving frame, $\mathcal{P}_\Phi(t) = \Big[\frac{t}{t_b}\Big]^{-\frac{Y^2}{8\pi^2}}~ e^{ \frac{Y^2}{4\pi^2}\,\big(t/t_b\big)^{1/4} } \,e^{-\Gamma_0\,(t-t_b)}~ \mathcal{P}_\Phi(t_b) $, with $\Gamma_0$ the decay rate at rest in Minkowski space-time. For an ultrarelativistic particle we find $\mathcal{P}_\Phi(t) = e^{-\frac{2}{3}\Gamma_0\,t_{nr}\,(t/t_{nr})^{3/2}}~ \mathcal{P}_\Phi(t_b)$ before it becomes non-relativistic at a time $t_{nr}$ as a consequence of the cosmological redshift. For $t\gg t_{nr}$ we find $\mathcal{P}_\Phi(t) = \Big[\frac{t}{t_{nr}}\Big]^{-\frac{Y^2}{8\pi^2}}~ e^{ \frac{Y^2}{4\pi^2}\,\big(t/t_{nr}\big)^{1/4} }~\Big[\frac{t}{t_{nr}}\Big]^{\Gamma_0 t_{nr}/2} \,e^{-\Gamma_0\,(t-t_{nr})}~ \mathcal{P}_\Phi(t_{nr})$. The extra power is a consequence of the memory on the past history of the decay process. We compare these results to an S-matrix inspired phenomenological Minkowski-like decay law modified by an instantaneous Lorentz factor to account for cosmological redshift. Such phenomenological description \emph{under estimates the lifetime of the particle}. For very long lived, very weakly coupled particles, we obtain an \emph{upper bound} for the survival probability as a function of redshift $z$ valid throughout the expansion history $\mathcal{P}_\Phi(z) \gtrsim e^{-\frac{\Gamma_0}{H_0}\,\Upsilon(z,z_b)}\,\mathcal{P}_\Phi(z_b)$, where $\Upsilon(z,z_b)$ only depends on cosmological parameters and $t_{nr}$.

A dynamical renormalization is introduced to describe the formation of a renormalized (quasiparticle) state. The renormalized survival probability P Φ (t) is ultraviolet finite, independent of the cutoff and decays on much longer time scales. During radiation domination, for a (quasi) particle "born" at time t b and decaying at rest in the comoving frame, , with Γ 0 the decay rate at rest in Minkowski space-time. The power with "anomalous dimension" and stretched exponential are remnants of the formation of the quasiparticle and consequence of the cosmological redshift. For an ultrarelativistic particle we find P Φ (t) = e − 2 3 Γ0 tnr (t/tnr ) 3/2 P Φ (t b ) before it becomes nonrelativistic at a time t nr as a consequence of the cosmological redshift. For t ≫ t nr we

I. INTRODUCTION
The decay and scattering of particles are some of the most fundamental processes in particle physics, within and beyond the Standard Model, with profound impact in cosmology. These processes are ultimately responsible for establishing a state of local thermodynamic and chemical equilibrium and are fundamental ingredients in kinetic processes in the early universe [1][2][3]. Particle decay is not only ubiquitous, but it plays an important role in big bang nucleosynthesis (BBN) [1,[4][5][6][7][8][9], and the generation of the baryon and lepton asymmetries [10][11][12][13]. The decay of long-lived dark matter particles is constrained by various cosmological and astrophysical probes [14][15][16][17][18], and recently it has been suggested that the two body decay of a long lived dark matter particle may relieve the tension between distance ladder and cosmic microwave background measurements of the Hubble constant [19].
Most treatments of particle decay (and/or inverse decay) in cosmology implement the S-matrix quantum field theory approach as in Minkowski space-time. In this framework, the unstable decaying state is prepared at a time far in the past (t → −∞), and one obtains the transition amplitude to a given final state far in the future (t → ∞). Taking the infinite time limit in the transition amplitude yields a total energy conserving delta function. Squaring this delta function to obtain the transition probability yields a total energy conservation delta function multiplied by the total elapsed time. Dividing by this large time and summing over all the final states for a given decay channel gives the total transition probability per unit time, namely a decay rate. Energy conservation, a consequence of the infinite time limit, yields kinematic constraints (thresholds) for decay and scattering processes.
In an expanding cosmology such an approach is at best approximate and at worst unreliable when the Hubble expansion rate is large even during a post-inflationary early stage of a radiation dominated cosmology. In a spatially flat Friedmann-Robertson-Walker (FRW) cosmology there are three space-like Killing vectors associated with spatial translational invariance and spatial momentum conservation, however, as a consequence of cosmological expansion there is no global time-like Killing vector, therefore particle energy is not manifestly conserved in scattering or decay processes.
A consistent formulation of dynamic processes in an expanding cosmology requires implementing methods of quantum field theory in curved space time [20][21][22][23][24][25][26][27][28]. Early studies revealed a wealth of novel phenomena such as particle production [20,23,24] and processes that are forbidden in Minkowski space time as consequence of strict energy conservation.
S-matrix theory was extended to simple cosmological space times to study the decay of a massive particle into two massless particles conformally coupled to gravity in ref. [29]. In references [30,31] these methods were adapted to calculate the decay of a massive bosonic particle at rest into two massless bosonic particles conformally coupled to gravity and into massless fermions Yukawa coupled to a scalar.
More recently [32] the decay of bosonic particles into two other bosonic degrees of freedom during a radiation dominated era was studied by implementing a non-perturbative method. This method was adapted to quantum field theory from the study of linewidths in quantum optics [33,34], combined with a physically motivated adiabatic expansion. While the results of this reference agreed with those obtained in ref. [30] for a particle decaying at rest in the comoving frame in the long time limit, they revealed new phenomena for highly relativistic decaying particles as a consequence of the cosmological redshift, and the relaxation of kinematic thresholds as a consequence of energy uncertainties determined by the Hubble scale.
Our study in this article is a natural extension of that in ref. [32] focusing on decay of a heavy bosonic particle into fermions, a more relevant case for standard model physics (and probably beyond) since most of the fermionic degrees of freedom in the standard model (with the possible exception of neutrinos) are Yukawa coupled to the Higgs boson.

Brief summary:
The study of fermionic degrees of freedom as decay products introduces several conceptually important distinctions with the bosonic case studied in refs. [30,32] that results in novel aspects of cosmological decay. First, fermionic degrees of freedom couple to the background gravitational field via the spin connection [20,25,28,[35][36][37][38][39][40][41][42][43][44][45]. Secondly, fermions Yukawa coupled to a bosonic degree of freedom yield a renormalizable theory. Recently the decay of a bosonic particle Yukawa coupled to fermions was studied within a non-perturbative real time framework in Minkowski space-time [46]. This study revealed novel transient dynamics associated with the dressing of the decaying particle by fermion-antifermion pairs into a quasiparticle state, which decays on a longer time scale. Such "dressing" leads to the necessity of an ultraviolet divergent renormalization of the decaying state and a detailed understanding of the various time scales to separate the many-particle dynamics of renormalization and dressing from that of the actual decay of the quasiparticle. Such dynamical effects cannot be addressed within an S-matrix framework since these effects are not secular in time and their contribution vanishes when the transition probability is divided by the total time in the infinite time limit. The dynamics of dressing and quasiparticle formation have been recently addressed in ref. [47] for a consistent interpretation of the reduction formula in asymptotic quantum field theory.
We introduce a dynamical renormalization that absorbs the ultraviolet divergences associated with fermion pairs into a renormalized survival probability at a renormalization time scale t b . The survival probability obeys a dynamical renormalization group equation with respect to t b . The cosmological redshift encodes the memory of the transient dynamics of quasiparticle formation in the decay law not seen in Minkowski space-time. If the decaying particle is ultrarelativistic, the decay dynamics depends crucially on t nr , the time scale at which it becomes non-relativistic as a consequence of the cosmological redshift. An S-matrix inspired, phenomenologically motivated, Minkowski-like decay law is shown to under estimate the lifetime of the decaying state.
Section (II) introduces the model and the adiabatic approximation, section (III) summarizes the non-perturbative framework to obtain the time evolution of the survival probability. In section (IV) we obtain the decay function for massless fermions during radiation domination, section (V) describes the dynamical renormalization method, section (VI) analyzes the decay dynamics of the renormalized survival probability during radiation domination, compares the results to an S-matrix inspired decay function, and introduces an upper bound to the decay function for very long lived, very weakly coupled particles valid all throughout the expansion history. Section (VII) discusses the various results analyzing their regime of validity and highlighting several implications. Section (VIII) presents our conclusions summarizing the main results. Various appendices contain technical details, in particular appendix (B) derives the decay law in Minkowski space time, highlighting the renormalization aspects to compare to the curved space-time case.

II. THE MODEL:
We consider a Higgs-like scalar field Yukawa coupled to one Dirac fermion in a spatially flat Friedmann-Robertson-Walker (FRW) cosmology with scale factor a(t) in comoving time. Generalizing to include Majorana fermions and/or more fermionic species is straightforward.
In comoving coordinates, the action is given by is the Ricci scalar, and ξ is the coupling to gravity, with ξ = 0, 1/6 corresponding to minimal or conformal coupling, respectively. Introducing the vierbein field e µ a (x) defined as where η ab is the Minkowski space-time metric, the curved space time Dirac gamma-matrices γ µ (x) and the fermionic covariant derivative D µ are given by [25,[35][36][37] where the γ a are the Minkowski space time Dirac matrices, chosen to be in the standard Dirac representation, and the covariant derivative D µ is given in terms of the spin connection by where Γ λ µν are the usual Christoffel symbols. For an (FRW) in conformal time dη = dt/a(t), the metric becomes where η µν = diag(1, −1, −1, −1) is the flat Minkowski space-time metric, and the vierbeins e µ a are given by The fermionic part of the action in conformal coordinates now becomes The Dirac Lagrangian density in conformal time simplifies to where i ∂ = γ a ∂ a is the usual Dirac differential operator in Minkowski space-time in terms of flat space time γ a matrices. Introducing the conformally rescaled fields and neglecting surface terms, the action becomes The effective time dependent masses are given by and In the non-interacting case, Y = 0, the Heisenberg equations of motion for the spatial Fourier modes with comoving wavevector k for the conformally rescaled scalar field are The Heisenberg fields are quantized in a comoving volume V , the real scalar field χ is expanded as where the mode functions g k (η) obey The mode functions are chosen to obey the Wronskian condition (II.19) and a, a † obey the usual canonical commutation relations.
For Dirac fermions the field ψ( x, η) is expanded as where the spinor mode functions U, V obey the Dirac equations [38][39][40][41][42][43][44][45] These equations become simpler by writing with U λ ; V λ being constant spinors [44,45] obeying Multiplying the Dirac equations on the left by γ 0 , it is straightforward to confirm that We choose the normalizations so that the operators b, b † , d, d † obey the canonical anticommutation relations. Furthermore, we will choose particle-antiparticle boundary conditions so that h k (η) = f * k (η) (see below). We note that for m f = 0 the conformally rescaled fermi fields obey the same equations as in Minkowski spacetime but in terms of conformal time, whereas this only occurs for bosons if they are conformally coupled to gravity, namely with ξ = 1/6, or for a radiation dominated cosmology (see below). The equivalence of massless fermions to those in Minkowski space-time will allow a direct comparison with the case of decay in flat space time studied in ref. [46] and summarized in appendix (B), and to interpret the differences with the curved space-time case.
A. Adiabatic approximation in post-inflationary cosmology: The standard (post-inflation) cosmology is described by radiation (RD), matter (MD) and dark energy (DE) dominated stages, we take the latter to be described by a cosmological constant.
Friedmann's equation in comoving time is where the scale factor is normalized to a 0 = a(t 0 ) = 1 today. We take as representative the following values of the parameters [48][49][50]: Passing to conformal time η with dη = dt/a(t), where the metric is given by (II.5) and C(η) ≡ a(t(η)), it follows that a eq is the scale factor at matter-radiation equality.
Hence the different stages of cosmological evolution, namely (RD), (MD), and (DE), are characterized by We will begin the study the dynamics of particle decay during the (RD) dominated era, generalizing afterwards to the case of a very long lived, very weakly coupled particle. During (RD) and (MD) we find, 36) and conformal time in terms of the scale factor is given by During the (RD) stage 38) and the relation between conformal and comoving time is given by 39) a result that will prove useful in the study of the decay law during this stage.

Bosonic fields:
Solving the mode equations (II. 18 We recognize that is the local energy measured by a comoving observer, and k p (t) is the physical wavevector redshifting with the cosmological expansion.
Writing the solution of (II.40) in the WKB form [23,[25][26][27][28] 43) and inserting this ansatz into (II.40) it follows that W k (η) must be a solution of the equation [25] (II. 44) This equation can be solved in an adiabatic expansion (II.45) 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 this is most easily recognized in comoving time t in terms of the comoving local energy (II. 41,II.42) and the Hubble expansion rate (II. 47) In terms of these variables, the first order adiabatic ratio (II.46) becomes [32] ω ′ k (η) .
is the local Lorentz factor.
The adiabatic approximation relies on the smallness of the (time dependent) adiabatic ratio corresponding to the physical wavelength ∝ 1/k p (t) and/or the Compton wavelength of the particle 1/m φ being much smaller than the size of the particle horizon d H (t) ∝ 1/H(t) at a given time.
During (RD) the particle horizon grows as a 2 (t) and during (MD) it grows as a 3/2 (t) whereas the physical wavelength grows as a(t). Therefore, if at a given initial time the adiabatic approximation is valid and H(t) ≪ E k (t) the reliability of the adiabatic expansion improves with the cosmological expansion.
To understand the origin of this approximation consider that the decaying particle is produced in the (RD) stage during which where g eff 100 is the number of ultrarelativistic degrees of freedom. Therefore, An upper bound on this ratio is obtained by considering that the decaying particle is produced at the scale of grand unification with T ≃ 10 15 GeV, assuming that this scale describes the onset of the (RD) era. Taking a typical comoving energy E k (t) ≃ T (t), one finds that H(t)/E k (t) 10 −3 and diminishes with cosmological expansion and diminishing temperature. This argument suggests that for typical particle physics processes the adiabatic ratio H(t)/E k (t) ≪ 1 throughout the postinflation thermal history.
In terms of this adiabatic ratio, we find where R(t) is the Ricci scalar (II.2). Furthermore, it is straightforward to find that therefore, this ratio is of second adiabatic order and can be safely neglected to the leading adiabatic order pursued in this study, justifying the simplification of the mode equations to (II.40) even for non-conformal coupling to gravity.
In this study we consider the zeroth-adiabatic order with the mode functions given by Since the decay function is ∝ Y 2 , keeping the zeroth adiabatic order yields the leading contribution to the decay law. Furthermore, as shown in detail in ref. [32] particle production as a consequence of cosmological expansion is an effect of higher order in the adiabatic expansion, thus it can be safely neglected to leading order.
The phase of the mode function has an immediate interpretation in terms of comoving time and the local comoving energy (II.41,II.42), namely which is a natural generalization of the phase of positive frequency particle states in Minkowski space-time.
During the (RD) era with C(η) given by (II.38) we find that the criterion (II.50) for the validity of the adiabatic approximation implies Fermi fields: The adiabatic expansion is straightforwardly applied to the fermionic case and has been discussed in the literature [39][40][41][42][43]. Beginning with the mode equations (II. 26,II.27) with , (II.58) therefore the purely imaginary term in these mode equations are of first adiabatic order and will be neglected to leading (zeroth) adiabatic order. Hence, to leading order we find In what follows we will refer to ω 2 k (η) = k 2 + M 2 (η) for both bosonic and fermionic degrees of freedom with M 2 (η) = m 2 C 2 (η) for either case. To leading (zeroth) order in the adiabatic expansion the Dirac spinor solutions in the standard Dirac representation and with the normalization conditions (II.29) are found to be To leading adiabatic order these spinors satisfy the completeness relations where the projector operators at different times Λ + k (η, η ′ ) ; Λ − k (η ′ , η) and their properties are given in appendix (A).

III. NON-PERTURBATIVE APPROACH TO THE DECAY LAW:
In Minkowski space-time, the decay rate of a particle is typically computed via S-matrix theory by obtaining the transition probability per unit time from an in-state prepared in the infinite past to an out-state in the infinite future. Obviously, such an approach -taking the infinite time limit-is not suitable in a time dependent gravitational background. An alternative approach in Minkowski space-time considers the Dyson-resummed propagator in frequency space that includes radiative corrections through the self-energy. The imaginary part of the self-energy evaluated on the mass shell in frequency space is identified with the decay rate, and a Breit-Wigner approximation to the full propagator, namely approximating the self-energy near the (complex) pole yields the exponential decay law. Again such an approach is not available in an expanding cosmological background where the lack of a time-like Killing vector prevents Fourier transforms in time-frequency, and makes the self-energy explicitly dependent on two time arguments, not only on the difference.
Instead we implement a quantum field theory method that complements and extends the Wigner-Weisskopf theory of atomic linewidths, that is particularly suited to study time evolution in time dependent situations. This method is manifestly unitary and yields a non-perturbative description of transition amplitudes and probabilities directly in real time. We summarize below the main aspects of the method as it applies to this study, referring the reader to [32][33][34] for details.
The total Hamiltonian in conformal time is given by H 0 +H I where H 0 is the free field Hamiltonian and is the interaction Hamiltonian in the interaction picture. Passing to the interaction picture wherein a given state is expanded in the Fock states associated with the creation and annihilation operators a, a † , b, d, etc. of the free theory, namely |Φ(η) I = n C n (η)|n , the amplitudes obey the coupled This is an infinite hierarchy of integro-differential equations for the coefficients C n (η). Consider that initially the state is |Φ so that C Φ (η i ) = 1 ; C κ (η i ) = 0 for |κ = |Φ , and consider a first order transition process |Φ → |κ to intermediate multiparticle states |κ with transition matrix elements κ|H I (η)|Φ . Obviously the state |κ will be connected to other multiparticle states |κ ′ different from |Φ via H I (η). Hence for example up to second order in the interaction, the state |Φ → |κ → |κ ′ . Restricting the hierarchy to first order transitions from the initial state |Φ ↔ |κ , and neglecting the contribution from vacuum diagrams which just yield a re-definition of the vacuum state 1 (see discussion in ref. [32]) results in the following coupled equations and inserting this solution into equation (III.3) we find We study the decay of a single particle bosonic state into a fermion-anti-fermion pair to leading order in the Yukawa coupling and the adiabatic approximation. Therefore the initial state is a single particle bosonic state with momentum k, namely |Φ ≡ |1 χ k . The set of states |κ with a non-vanishing matrix element of H I with this single particle state are |κ ≡ |1 f p,λ , 1 f q,λ ′ where λ, λ ′ are the polarization of the fermion and antifermion states. The matrix elements entering in the evolution of the amplitudes are and the self-energy (III.7) to leading order in the adiabatic expansion becomes where q = k − p. This is the fermionic one-loop self energy in curved space time to leading order in the adiabatic expansion.
Obviously the differential equation (III.6) cannot be solved exactly with the above self-energy.

In Minkowski space time the self-energy is a function of the time difference allowing a solution via
Laplace transform [33,34]. However, in a time dependent expanding cosmology such an approach is not available. This is a consequence of the lack of a global time-like Killing vector. Instead for weak coupling we resort to a Markov approximation [32]. While details are available in ref. [32][33][34] to which the reader is referred, we summarize here the main aspects of this approximation.
We begin by introducing with the condition Then (III.6) can be written as which can be integrated by parts to yield the first term on the right hand side of (III.14) is of order Y 2 , whereas the This expression clearly highlights the non-perturbative nature of the Wigner-Weisskopf approximation. The imaginary part of the self energy Σ Φ yields a renormalization of the adiabatic frequencies and will not be addressed here [33,34], whereas the real part determines the decay law for the survival probability directly exhibits the non-perturbative nature of the method. The self-energy is given by (III.9) to leading order in Yukawa coupling.

IV. MASSLESS FERMIONS:
Our goal in this article is to study the decay of a heavy Higgs-like scalar field into much lighter fermions, neglecting the fermion masses. This is a suitable scenario for the standard model where the Higgs scalar can decay into all the charged leptons and quarks but for the top, and the quark and lepton masses may be safely neglected. Such scenario also includes the possibility of decay into neutrinos in the case that neutrino masses originate in Yukawa couplings to a Higgslike scalar beyond the standard model. We postpone the study of decay into heavier fermionic degrees of freedom to a companion article. Focusing on the case of massless fermions allows a direct comparison with results in Minkowski space time, which are summarized in appendix (B).
Furthermore, understanding this simpler case provides a pathway towards the more general case of massive fermions to be studied elsewhere.
For massless fermions ω ψ k (η) = k, in this case the projector operators Λ ± in (III.9) are given by eqn. (A.12) in appendix (A), and the self-energy (III.9) can be written in dispersive form as where the spectral density is given by We carry out the k 0 integral in (IV.1) by introducing an upper (comoving) ultraviolet cutoff Λ and a short time convergence factor η − η ′ → η − η ′ − iǫ with ǫ → 0 + and replacing k 2 yielding the final result for the self-energy In our analysis we will keep Λ fixed but large and take the limit ǫ → 0 + first, clearly this is the correct limit when the theory is considered as an effective field theory valid below a cutoff Λ. We note that the flat space time limit is obtained by replacing η → t, and the frequency ω φ k to be time independent (see appendix (B)).
It remains to perform the time integrals to obtain Γ Φ (η) and The total time derivative in (IV.4) is integrated by parts and consistently with keeping the leading order in the adiabatic expansion, terms of the form ω ′ /ω 2 are neglected since these yield higher order adiabatic corrections. In the limit ǫ → 0 + for fixed Λ we find the decay function In obvious notation the contributions I 2b (k, η), I 3b (k, η) are the Λ independent terms in I 2,3 respectively. These three contributions are studied separately below, analyzing their cutoff dependent and independent terms extracting the different physics of each term.
A. Analysis of I 1,2,3 : In the following analysis we will take the cutoff Λ to be the largest of all scales, in particular Λ ≫ ω k (η) at all times.
I 1 : I 1 vanishes identically as η → η i and the oscillatory terms become negligibly small for Λ(η − η i ) ≫ 1, therefore I 1 grows to its asymptotic value very rapidly, on a time scale η − η i ≃ 1/Λ. This divergent contribution corresponds to a renormalization of the amplitude and is similar to a linearly divergent renormalization in Minkowski space time [46] (see appendix (B)).
The technical details of the analysis of I 2 are relegated to appendix (D). The main result is where γ E = 0.577 · · · is Euler's constant and I 2b (k, η) is given by equation where this contribution is analyzed in detail. We discuss this contribution in further detail in sections (V,VI) below.
With Λ ≫ ω k the argument of the sine function in the first term in eqn. (IV.8), namely in , (IV.12) therefore in this limit we find /m 2 being the Lorentz factor whose time dependence is a consequence of the cosmological redshift.
In appendix (E) we provide the analysis for I 3b , gathering both terms we find that where S(k, η ′ ) is given by (E.25) with asymptotic limit S(k, η ′ ) → 1 for large η ′ . Therefore I 3 = I 3a + I 3b does not depend on Λ in the limit Λ → ∞. This is similar to the case in Minkowski ).

V. RENORMALIZATION: DYNAMICS OF "DRESSING"
The final result for the decay function in (IV.5), I(Λ, k, η) is given by , where I f in (k, η) is independent of the cutoff Λ in the limit Λ → ∞, and for (η − η i ) ≫ 1/Λ it is given by The linear and logarithmic dependence on the cutoff Λ are exactly the same as in Minkowski space time [46], as obtained in the appendix (B). This similarity is expected as the cutoff dependence arises from the short distance behavior of the self-energy correction which should be insensitive to the curvature of space time. As discussed in ref. [46] the origin of this divergence is the "dressing" of the bare single particle state by a cloud of fermion-anti-fermion pairs into a renormalized quasiparticle state. In a renormalizable theory the growth of the density of states at high energy implies that this cloud of excitations contains high energy states. The dynamical build-up of the cloud of excitations occurs on a time scale η − η i ≃ 1/Λ at which the divergent contributions to I 1,2 saturate, see eqn.
(IV.6) and the discussion in appendix (D).
The "dressing" of the bare into the physical renormalized quasiparticle state is accounted for by the wave-function renormalization of the amplitude [46]. For large cutoff scale Λ and for a weakly coupled theory with Y 2 ≪ 1 there is a wide separation betweeen the time scales of formation of , which for weak coupling is the longest scale. Therefore, we can evolve the initial state in time up to an intermediate time scale , so that the initial state had enough time to be "dressed" by fermion-antifermion pairs into the renormalized quasiparticle state, but did not have time to decay. For example, taking η b − η i = 1/ω φ k (η i ) fulfills the conditions of time scale separation because ω φ k ≪ Λ, and because for Y 2 ≪ 1 there will be many oscillations of the field before it decays. Taking this renormalization scale is tantamount to an "on-shell" renormalization scheme. We identify η b as the time of formation -or "birth" -of the "dressed" or quasiparticle state [46], which after formation decays on a much longer time scale.
The time evolution of the "bare" single particle state until it is renormalized or "dressed" is implemented by the following procedure. Writing where, taking (η b − η i ) ≫ 1/Λ, the subtracted quantity The subtracted contributions and are obtained explicitly in appendices (D,E) respectively. During (RD) we find (see appendix (D) for definitions and eqn. (D.12)) J(ξ ′ ) is given by eqn. (D.9) in appendix (D), and where S(η) is given by eqn. (E.25) in appendix (E). The contribution from I(Λ, k, η b ) is absorbed into wave-function renormalization Z as follows. Writing equation (III.17) as where the renormalized probability is given by The exponent in the wave function renormalization Z(η b ) is given by yielding an ultraviolet divergent wave function renormalization. The renormalized probability obeys The decay function that describes the time evolution of the renormalized survival probability is given by it is finite and independent of Λ in the large cutoff limit. The time scale η b acts as a renormalization scale, obviously the survival probability P Φ,r (η) is independent of this renormalization scale, hence it obeys a dynamical renormalization group equation, namely 14) The solution of this equation is, obviously 2 , describes the probability of the renormalized quasiparticle state. This "dressed" state decays with the finite and cutoff independent decay function η η b Γ Φ (η ′ )dη ′ on time scales much longer than the "dressing" or renormalization scale η b .
In the following analysis we will drop the subscript r from P Φ,r to simplify notation since we will be strictly dealing with the renormalized survival probability.
The decay function (V.13) depends explicitly on the initial time η i (see explicit expressions in appendix (D)). However, P Φ,r (η b ) is defined at the renormalization scale η b and it is taken to be the initial probability of the fully renormalized state after all the short time transient dynamics that result in the "dressing" of the bare into the renormalized quasiparticle state have subsided.
Therefore, the dependence of the contributions (V.6,V.8) on η i must be traded for a dependence on η b .

Let us write
with β ≫ 1 so that the Λ dependent terms in I 1,2 reached their asymptotic behavior. For example, the "on-shell" renormalization scheme corresponds to β ≡ Λ/ω k (η i ). Therefore, in terms of the Hubble rate and the physical cutoff Λ ph (η i ) = Λ/C(η i ) at the initial time H(η i ) we find in (RD) Since the cutoff scale Λ is taken to be much larger than any of the energy scales and the adi- the second term in the bracket is at most of first adiabatic order, this is the case for the "on-shell" renormalization scheme for which β ω k (η b )/Λ = 1. Hence, to leading adiabatic order we can safely replace ω k (η i ) → ω k (η b ) in the expressions. Using the results of appendix (D) we find that similar arguments justify the replacement γ k (η i ) → γ k (η b ) along with η i → η b in all the quantities that enter in the decay function. In the limit of large cutoff Λ the trade-off between the variables at the initial time and those at the renormalization scale η b does not depend on the cutoff as it must be for a consistent effective field theory description well below the cutoff scale. We note that the adiabatic approximation plays an important role in this separation and is a necessary ingredient because the frequencies depend on time unlike in Minkowski space-time. In particular for the "on-shell" renormalization scheme because the adiabatic condition (during (RD)) corresponds to ω k (η i ) η i ≫ 1 (see equation (II.57)).

VI. DYNAMICS OF DECAY.
Once we have absorbed the ultraviolet divergences into a renormalization of the amplitude, we now proceed to analyze the main physical aspects of the decay dynamics leveraging the adiabatic approximation.
A. Decay during radiation domination: We assume that the decaying particle has been produced early during the (RD) stage by some (unspecified) particle physics process at a high energy/temperature scale, focusing first on the dynamics of decay during this era. The subtracted decay function I S (k, η, η b ) (V.4) can be written in a compact manner amenable to a numerical study as and where ξ, W [ξ] are defined in eqn. (V.7), and we have introduced the following functions (see appendices D,E) where J[ξ] is defined in eqn. (D.9) in appendix (D). To leading adiabatic order (see appendix (E)) where Si[x] is the sine-integral function (see equation (E.26) and discussion in appendix (E)).
We highlight that the contribution I R S is a distinct feature of the renormalizable Yukawa interaction and of the fermionic density of states, whereas I 3S in (V.8) is very similar to the decay function found in the scalar case studied in ref. [32].
As discussed above, to leading adiabatic order we set η b = η i in I 3S and obtain (see appendix where I 3S (k, η, η b ) must be obtained numerically.
However, before we engage in a numerical study we analyze the different contributions to extract a physical picture of which terms dominate at different time scales. In order to analyze the behavior in the different regimes, we write the Lorentz factor both in terms of the variable ξ = η η i − 1 (see appendix (C)) as well as in terms of comoving time with the equivalence 1 + ξ ≡ t/t i (see also appendix (C)), where t nr is the comoving time scale at which the decaying particle becomes non-relativistic, given by (VI.10) Whence the limits Let us focus first on the contribution I R S (k, η, η b ) given by (VI.2). In Minkowski space-time the frequencies are time independent, therefore W [ξ ′ ] = 1 and J(ξ ′ ) = (ω k − k)η i ξ ′ . The analysis of appendix (B) shows that in Minkowski space time for ξ ≫ 1 the second term in (VI.2), namely ii) ξ m < ξ γ i : during this interval the function F 1 [ξ, ξ b ] continues to rise monotonically a behavior summarized by figure (7)  (approximately) cancelling the logarithm from the first term in I R S , whereas F 2 ≃ 2 ln[ξ m /ξ b ] remains constant, yielding a plateau in I R S . This approximate cancellation is effective during a time interval that increases for γ i ≫ 1 (see discussion in appendix (D)). According with eqn. (VI.9) and the limits (VI.11) during this interval, wherein I R S is approximately constant, the decaying particle is in the ultrarelativistic regime. In this stage the constancy of I R S is expected because in the ultrarelativistic regime the frequencies are nearly time independent since ω k (η) ≃ k ≃ ω i . Therefore W [ξ] ≃ 1 yielding F 1 ≃ 2 ln[ξ/ξ b ] thereby cancelling the logarithmic time dependence of the first term in (VI.2), similarly to Minkowski space-time.
If γ i ≫ 1 the decaying particle is "born" ultrarelativistically and there is a (long) time window thereby approximately cancelling the first term in I R S whereas F 2 [ξ, ξ b ] remains nearly constant. Therefore for γ i ≫ 1 it follows that I R S (k, η, η b ) rises rapidly on a time scale ≃ ξ m reaching a maximum and remaining nearly constant iii) ξ ≫ γ i : The cosmological redshift eventually makes the decaying particle to become nonrelativistic when ξ ≫ γ i ≫ 1. During this stage the particle is non-relativistic as a consequence of the cosmological redshift. The time dependence of the frequency now yields W [ξ ′ ]+1/ W [ξ ′ ] ≫ 2, hence F 1 > 2 ln [ξ]. In this stage it follows that W [ξ] ≈ ξ/γ i , therefore for ξ ≫ γ i ≫ 1 we find by splitting it into the stages ξ b ≤ ξ ≤ γ i and ξ > γ i . The first stage yields 2 ln[γ i /ξ b ] since during this (ultrarelativistic) stage W [ξ ′ ] ≃ 1, and the second yields (approximately) 2 ξ/γ i since during In summary, for a particle that is "born" ultrarelativistically, namely with γ i ≫ 1, the contribution I R S rises rapidly up to a value ≃ 2 ln[ξ m /ξ b ] on a time scale ξ m ≪ γ i given by (D.16), remains nearly constant up to a time scale ξ ≃ γ i at which the particle becomes non-relativistic, and begins to fall-off as −2 ξ/γ i for ξ ≫ γ i .
In the opposite limit when γ i ≃ 1 the decaying particle is non-relativistic already at the initial time and ω k (η) ≃ m φ C(η). In this case F 2 [η, η b ] saturates rapidly, on a scale ξ m ≃ π/ω i η i ≪ 1, and F 1 [η, η b ] grows faster than logarithmically, hence F 1 − F 2 becomes larger than the logarithm in the first term of I R S and negative. This behavior leads to an early suppression of decay. This analysis is approximately summarized during the ultrarelativistic (UR) and non-relativistic (NR) regimes, by (see eqn. (D. 19) in appendix (D)), (VI.12) The main aspects of this analysis are confirmed by a numerical study summarized in figures (1) and (2) for γ i = 2, 10 respectively. Notice the different scales in the figures highlighting the emergence of the plateau and the crossover to a diminishing (negative) square root behavior at a scale ξ ≃ γ i .
Decay at rest: for a very massive particle "born" and decaying at rest in the comoving frame, namely for γ i = 1, and ω i η i ≫ 1 we can provide an analytic form of the decay function for time scales ξ ≫ ξ b ≃ 1/ω i η i for on-shell renormalization. As discussed in appendices (D,E), (VI.14) Setting η i = η b to leading adiabatic order, we find for γ i = 1 and t ≫ t b leading to the survival probability for t ≫ t b Decay of particles with γ i ≫ 1: these are particles that are "born" ultrarelativistically. For γ i ≫ 1 the contribution from I R S (k, η, η b ) has been summarized by eqn (VI.12) and is displayed in fig. (2): a rapid rise on a time scale ξ m ≪ γ i given by (D.16) up to I R S ≃ 2 ln[ξ m /ξ b ] followed by a near plateau during the stage while ξ γ i . This contribution falls off slowly as − ξ/γ i during the non-relativistic stage, ξ ≥ γ i (see eqn. (VI.12)). While a quantitative analysis of I 3S requires a numerical study, we can obtain a fairly accurate estimate as follows. The contribution from S to I 3S (see equation (VI.3)) is discussed in appendix (E), and can be approximately summarized as: S ≈ 0 for ξ < ξ s and S(η) ≈ 1 for ξ > ξ s with ξ s given by (E.28,E.29).
With γ i ≫ 1, the ultrarelativistic stage corresponds to γ i ≫ ξ, during the stage γ i ≫ ξ s ≫ ξ it follows that S ≈ 0, using 1 + ξ = t/t i , and eqns. (VI.10,VI.14), during this stage I 3S is given in also describes the decay function in the case of a scalar field decaying into two massless scalars [32].
During this stage for t ≪ t nr we find For γ i ≫ ξ ≫ ξ s it follows that S ≃ 1, therefore the above result is multiplied by a factor 2.
Hence, during the ultrarelativistic stage with γ(t) ≫ 1, or t ≪ t nr , and S = 1 in (VI.3), it follows that which when combined with the result (VI.12) yields in this ultrarelativistic regime, for Neglecting the perturbatively small non-secular constant in the decay function from the first term in (VI.21) 3 , we find in the time interval for t nr ≫ t ≫ t b during which the decaying particle is ultrarelativistic and the transient dynamics of quasiparticle formation has saturated We can now use the property (V.15) and write for t > t nr (VI.24) After the decaying particle becomes non-relativistic for ξ ≫ γ i or t ≫ t nr when γ(t) ≃ 1, the contribution S ≃ 1 and I 3S (ξ) − I 3S (ξ nr ) becomes the dots in the above expression stand for terms of higher order in the ratio t nr /t.
Finally, combining with the result given by eqn. (VI.12), the total decay function after the particle has become non relativistic ξ ≫ γ i (or t ≫ t nr ≫ t b ) is given in comoving time by where we have neglected a perturbatively small constant term and approximated t i γ 2 i ≃ t nr for γ i ≫ 1. Hence for t ≫ t nr ≫ t b we find It would be expected that after t nr when the particle has become non-relativistic as a consequence of the cosmological redshift, the time evolution of the survival probability would be similar to that of a particle born and decaying at rest. However, the result (VI.27) features an extra power law with exponent Γ 0 t nr /2 as compared to the decay function for the particle born at rest, eqn. (VI.16). This difference reflects the memory of the past evolution in the form of the integral (VI.24).
We can provide a measure of the impact of curved space time effects on the decay function by comparing the results above to a phenomenological, S-matrix inspired Minkowski decay law allowing for a local time dilation factor to account for the cosmological redshift, namely is the decay width at rest in Minkowski space time, and γ(t) the local Lorentz factor (VI.9). The comparison to the cutoff independent subtracted decay function (VI.1) is facilitated by introducing so that the Minkowski-like decay function is given by where a factor is included in (VI.30) to ensure that I M (t i = t b ) = 0 consistently with the subtraction definining (VI.1). For t ≫ t i this phenomenological decay function is interpreted as that of Minkowski space-time but with the instantaneous Lorentz time dilation factor. For t ≫ t i it provides a "benchmark" to compare the results obtained above for the decay function to an S-matrix inspired instantaneous Minkowski decay law.
Before we engage in a numerical comparison, it is illuminating to analyze the cases discussed above.
Non This can be understood from the following argument.
As discussed above and in appendix (D), for γ i ≫ 1 the contribution I R S (see eqn. (VI.2)) rises on a time scale ξ m ≃ (3πγ 2 i /ω i η i ) 1/3 up to a maximum ≃ 2 ln(ξ m /ξ b ) after which it remains nearly constant up to ξ ≃ γ i yielding the logarithmic term in (VI.21). For example, for γ i ≃ 200 , ω i η i ≃ 100 and "on-shell" renormalization with ξ b = 1/ω i η i , the contribution from I R S rises up to a value 7 on a comoving time scale t m /t i ≈ 240. Since the Hubble time scale 1/H(t) = 2t during (RD), it follows that I R S rises up to the plateau over ≃ 240 Hubble times, with the possibility that during this time I S (t) > I M (t). However, after the particle has become nonrelativistic, namely for t ≫ t nr , the cosmological decay function I S (t) is given by (VI.26) whereas showing that I S (t) ≪ I M (t) for t ≫ t nr . This suggests a crossover behavior for very large values of γ i : there is an early time window during the ultrarelativistic stage wherein the cosmological decay function may be larger than the Minkowski one, however as the decaying particle eventually becomes non-relativistic the latter will ultimately dominate. This behavior is borne out by a detailed numerical study. However, for γ i = 200 fig. ( 5) shows that the contribution from I R S dominates at early time, rising on a time scale t/t i ≃ 100. In this case I M is smaller than I S during a substantial time window, ≈ 500 Hubble times from the "birth" of the quasiparticle, before crossing over to becoming the largest decay function.
Therefore we conclude that in the ultrarelativistic case, for very large values of γ i , the decay function is larger than the phenomenological Minkowski one within a substantial time interval but eventually becomes smaller at a time scale that depends on the various parameters. In either case, at long time the decaying particle lives longer than predicted by a Minkowski decay law The discussion above focused on decay during the radiation dominated era that lasts until C(η) = a eq ≃ 10 −4 corresponding to an ambient temperature T ≃ eV at a time t eq ≈ 10 12 secs. If the decaying particle is very long lived as would befit a dark matter candidate, it would continue to decay during the matter and perhaps dark energy dominated eras. This case corresponds to an extremely small Yukawa coupling, which allows to safely neglect early transient effects that . This is because both terms saturate on short time scales therefore they yield perturbatively small corrections to the decay function for very weak Yukawa coupling as compared to the terms that continue to grow in time. Hence, neglecting these perturbatively small transient contributions for very weak Yukawa couplings, the decay function simplifies to where ξ = (η − η i )/η i and F 1 is given by (VI.4) and the dots in (VI.32) stand for constant terms For general scale factor W [ξ] is given by With "on-shell" renormalization (ξ b = 1/ω i η i ), we find Taking the initial time to correspond to an initial temperature 10 15 GeV yields C(η i ) ≃ 10 −28 therefore the logarithm contribution to the decay function for η ≃ √ a eq /H R yields Obtaining the contribution from F 1 over the whole history from early (RD) into (MD) can be done numerically, although this is a rather challenging task because of the enormous dynamic range with the scale factor varying over twenty four orders of magnitude. However, we can provide a simple estimate of the remaining two terms of the decay function at long time during the (MD) era and/or beyond. If the particle remains ultrarelativistic, then as discussed in the previous sections the contribution from F 1 cancels the logarithmic time dependence of the first term, hence the combination of the first two terms saturates (this is the plateau in fig. (2)) and yields a perturbatively small time independent contribution to the decay function. Hence during this ultrarelativistic stage the last term in (VI.32) dominates the decay function.
After the particle has become non-relativistic then W [ξ] ≃ C(η)/γ i C(η i ) ≫ 1 and during (MD) using eqn. (II.35) and taking as an upper bound η ≃ √ a eq /H R we find Finally, we can estimate the last term in (VI.32) during the stage when the particle is nonrelativistic and (MD) dominated, taking γ(η ′ ) ≃ 1 and taking η ≃ √ a eq /H R , we find Since during the ultrarelativistic stage the time dependence of the first and second term cancel out and the last term dominates the decay dynamics, we conclude that the last term in (VI. 32) dominates the decay dynamics of a very long-lived particle with very weak Yukawa coupling, all throughout the time evolution. Since the first (logarithmic) term is always subdominant, and the second term is negative, larger in magnitude than the logarithmic term but also subdominant at ; a nr ≡ k/m φ , we find that the upper bound to the decay function at redshift z is given by where Γ 0 = Y 2 m φ /8π is the decay rate at rest in Minkowski space time, and depends solely on the cosmological parameters and a nr = k/m φ the scale factor at which the decaying particle transitions from ultrarelativistic to non-relativistic, and we have taken z b ≫ 1.
The redshift evolution of the survival probability all throughout the expansion history is summarized concisely as The inequality in eqn. (VI.41) reflects that eqn. (VI.39) yields an upper bound to the decay function. For a nr = 0, namely when the decaying particle is "born" at rest, it follows that Decay at rest in the comoving frame (γ i = 1): The time evolution of the survival probability is given by where Γ 0 = Y 2 8π m φ is the decay width of a particle at rest in Minkowski space-time. The power law and stretched exponentials are both a remnant of the renormalization, or "dressing" of the bare into the quasiparticle state and a distinct consequence of the cosmological redshift. Indeed, in Minkowski space-time the terms that give rise to these contribution become time independent after the transient dynamics, whereas, in curved space-time, the origin of these contributions is the time dependence of the frequencies via the cosmological redshift.
The methods that we implemented in this study, a non-perturbative formulation combined with a physically motivated adiabatic expansion including a consistent treatment of renormalization, are very different from those implemented in ref. [31]. The decay law of a particle at rest (VII.1) is also very different from that reported in ref. [31]. The origin of the discrepancy is not clear to us.
However, since the power law and stretched exponentials originate precisely from the contributions to the renormalization of the survival probability, we suspect that the discrepancy originates in the treatment of the ultraviolet divergences. These are of the same form as in Minkowski space-time (see appendix (B) and ref. [46]) as expected since these are short distance divergences, but have not been discussed or addressed in ref. [31]. As explained above, the time dependence of the frequency yields an unexpected contribution to the decay law on longer time scales that originates in the dynamics of quasiparticle formation.
Born ultrarelativistically: if the particle is "born" or produced ultrarelativistically, namely , which determines when the particle transitions from being ultrarelativistic (γ(t) ≫ 1 or t ≪ t nr ) to non-relativistic (γ(t) ≃ 1 or t ≫ t nr ) as a consequence of the cosmological redshift. The dynamical evolution of the survival probability is different in these stages. a) ultrarelativistic stage: Although for t ≫ t nr the particle has become non-relativistic because of the cosmological redshift, as compared to the case of decay at rest (VII.1), this decay law features a new power with exponent Γ 0 t nr /2. Its origin is the memory of the decay function manifest in the form of the integral of the cosmological redshift in (V.8) over the whole history of the decay process. Therefore, even well after the decaying particle has become non-relativistic, the survival probability features an enhancement factor that "knows" about the past history when the particle was ultrarelativistic. The dynamics during the transition from the ultrarelativistic to the non-relativistic behavior must be studied numerically, and the previous section shows such study for several values of the parameters.
Massless fermions vs massless bosons: Ref. [32] studied the decay of a scalar into two massless scalars, therefore we can now compare the results of that study to those obtained here for the case of scalar decay into massless fermions. The main difference is in the contribution I R S in eqn.(VI.1) which is given by eqn. (VI.2). The contribution from I 3S to the decay function is the same for fermions and bosons, for example the function G[x] is the same that enters in scalar decay [32]. The extra contribution, namely I R S has the same origin as the ultraviolet divergent contributions that are absorbed in wave function renormalization. This is also the case in Minkowski space-time [46] as shown in appendix (B). Whereas in Minkowski space time this contribution becomes time independent after a short time transient and is absorbed into wave function renormalization, in a FRW cosmology, it is time dependent as a consequence of the cosmological redshift and becomes important for non-relativistic particles. Namely, I R S is a remnant of the physical process of quasiparticle formation. There is no such contribution in the case of decay into two scalars because the theory in this case is superrenormalizable, hence there is no equivalent of the I R S term. This contribution suppresses the decay function at long time, thereby enhancing the lifetime of the decaying particle. This behavior is yet another source of discrepancy with the results of ref. [31], which finds a larger rate in the fermionic case. The source of this discrepancy are precisely the "anomalous" power and stretched exponential which are a consequence of the quasiparticle formation and wave function renormalization. Although the decay probability requires an ultraviolet divergent wave function renormalization even in Minkowski space-time, this seems to be an aspect missing in the treatment of ref. [31]. The cumulative effect of these differences results in that a meaningful comparison to our study has eluded us.
For decay at rest γ(t) = 1, this decay law misses the power with anomalous dimension and the stretched exponential, whose combination is negative. Therefore the Minkowski-like decay law overestimates the suppression of the survival probability in the case of decay at rest. For a particle that is produced ultrarelativistically, during the stage wherein γ(t) ≫ 1, namely t ≪ t nr one finds which is smaller than (VII.2). For t ≫ t nr when the decaying particle has become non-relativistic This analysis also suggests that the corrections to the decay law are more important for particles "born" very early during (RD) and very long lived a situation that befits most descriptions of a dark matter candidate.
Caveats: We have focused on studying scalar decay into massless fermion pairs, a situation that approximates most of the fermionic decay channels of a Higgs scalar in the standard model.
An important aspect of this decay process is that it does not feature thresholds. Including the mass for the decay products introduces kinematic thresholds, a consequence of strict energy-momentum conservation. In ref. [32] it was argued that the Hubble rate of expansion introduces a natural energy uncertainty leading to a relaxation of the kinematic thresholds, thereby allowing processes that are forbidden in Minkowski space-time by energy conservation. Furthermore, ref. [46] has shown that energy uncertainties associated with transient non-equilibrium aspects of the decay allow decay into heavier particles during a time interval. In an expanding cosmology these effects may combine with the energy uncertainty from Hubble expansion to enhance the decay by opening up novel channels that would be otherwise forbidden by strict energy conservation. These aspects associated with the masses of the decay products will be the subject of further study.
The inclusion of masses for the decay products becomes a more pressing issue In this study we have neglected finite temperature corrections to decay vertices and masses, their inclusion requires studying the time evolution of an initial density matrix. Furthermore, if the decay products thermalize with the medium, their population build-up will lead to Pauli blocking factors thereby suppressing the decay of the parent particle. These effects remain to be studied but are beyond the scope and goals of this article.
Possible implications: The time dependence of the decay function reveals non-equilibrium aspects that have not been previously recognized, not only from the transient build-up of the quasiparticle but also the memory effects that yield the unexpected power laws and stretched exponentials. These novel non-equilibrium effects may lead to interesting and perhaps important dynamics relevant to baryogenesis and leptogenesis. In particular, we envisage corrections to quantum kinetic processes for particle production and their inverse processes. Typically quantum kinetics inputs transition rates perhaps with finite temperature contributions but ultimately obtained from S-matrix theory. Namely, such transition rates are obtained in the infinite time limit and the forward and backward probabilities input strict energy conservation, and as a consequence, obey detailed balance . The richer time dependence of the decay function revealed by this study, with the hitherto unexplored novel non-equilibrium aspects, suggests that similar dynamical processes may enter in a modified quantum kinetic description in the early universe. We expect to report on these and other related issues in future studies.

VIII. SUMMARY, CONCLUSIONS AND FURTHER QUESTIONS.
In this article we studied the decay of a bosonic particle into massless fermions via a Yukawa coupling in post-inflation cosmology. The approximation of massless fermions is warranted for a heavy Higgs-like scalar within or beyond the standard model decaying into most charged leptons or quarks (but for the top) of the standard model. We implemented a non-perturbative method that yields the time evolution of the survival probability P Φ (t) combined with a physically motivated adiabatic expansion. This expansion is justified when H(t)/E k (t) ≪ 1 where H(t) is the Hubble rate and E k (t) the local energy of the particle as measured by a comoving observer. We have argued that this approximation is valid for typical particle physics processes during the radiation dominated era and beyond. In a standard cosmology the reliability of this approximation improves with the cosmological expansion, therefore if the adiabatic condition is fulfilled at the initial time when the decaying particle is produced, its reliability improves along the expansion history. We carried out a detailed analytic and numerical study of the decay function during the radiation dominated era. The dynamics of decay depends crucially on whether the particle is non-relativistic or relativistic. For a particle that is "born" at rest in the comoving frame during (RD) we find that after short time transients, the survival probability is given by where Y is the Yukawa coupling and Γ 0 is the decay rate at rest in Minkowski space-time. For the case in which the decaying particle is "born" ultrarelativistically the time evolution over the whole history during (RD) must be obtained numerically. Different regimes emerge depending on whether the particle is ultrarelativistic for t ≪ t nr or non-relativistic for t ≫ t nr where is the time scale at which the decaying particle of mass m φ that is born ultrarelativistically with comoving momentum k transitions to being non-relativistic as a consequence of the cosmological redshift. During the ultrarelativistic regime (t ≪ t nr ) we find for t ≫ t b that the decay function is a stretched exponential whereas for t ≫ t b and after the particle has become non-relativistic (t ≫ t nr ) we find The extra power of t/t nr as compared to the case when the particle is born at rest (see eqn. (VIII.1)) is a consequence of the memory of the decay function on the past history during the ultrarelativistic stage.
The cosmological decay law is compared to a phenomenological Minkowski-like, S-matrix inspired decay law with an instantaneous Lorentz time dilation factor P (M ) we found that this phenomenological law describes at long times a much faster decay thereby under estimating the lifetime of the decaying particle.
The decay dynamics revealed by this study during (RD) allows us to extrapolate to the case of very long lived, i.e. very weakly coupled particles. We obtain a decay function that yields an upper bound to the survival probability all throughout the expansion history under the assumption of two body decay into a massless fermions, it is given by where Υ(z, z b ) is given by (VI.39) and depends only on the cosmological parameters and the scale factor at which the particle transitions from ultrarelativistic to non-relativistic.
One important conclusion from these results is that using a decay rate as measure of the decay dynamics is not a useful concept and misses the correct dynamical evolution. An S-matrix calculation of transition amplitudes or probabilities, where the time interval is taken to infinity not only it does not capture the various different dynamical scales and temporal behaviour of the survival probability, but substantially under estimates the lifetime of the decaying state.
An important corollary of this study is that the S-matrix approach to describe quantum decay in the cosmological setting is in general inadequate, while it may yield a good approximation for processes of decay at rest for weakly coupled particles late in the cosmological history, it misses important non-equilibrium dynamics. The non-equilibrium effects revealed by our study, from the transient dynamics of the formation to the quasiparticle, to the memory of the decay function on the past history of the decaying particle could be relevant in the quantum kinetics of processes in the very early universe. These could have potential impact in CP-violating non-equilibrium dynamics, baryogenesis and leptogenesis and merit further study.
where I is the 2 × 2 identity matrix. These expressions can be written more compactly introducing the following functions (suppressing the momentum and conformal time arguments), Two relevant cases: 1:) Equal time, η = η ′ , 2:) Massless fermions, Appendix C: Useful identities: In this appendix we gather some useful identities valid during the radiation dominated stage (see also appendix (D)).
The local Lorentz factor in conformal time is given by yielding the identity The relationship with comoving time t is obtained via eqn. (II.39), namely The conformal and comoving time scales η nr , t nr respectively, determine the scale at which the decaying particle transitions from relativistic with γ(η) ≫ 1 for η ≪ η nr or t ≪ t nr to nonrelativistic with γ(η) ≃ 1 for η ≥ η nr or t ≥ t nr . In terms of η, η i we find, Consider the first term in I 2 eqn. (IV.7), For Λ ≫ k, m φ the argument of the cosine becomes simply Λ (η ′ − η i ). Define x = Λ (η ′ − η i ), and change integration variable to x, with x f = Λ (η − η i ), yielding In the limit Λ → ∞ we find where γ E = 0.577 · · · is Euler's constant and for x f ≫ 1 we find Ci(x f ) = sin(x f )/x f + · · · . We confirmed the result (D.3) numerically. Therefore for Λ ≫ k, m φ , 1/(η − η i ) we find Let us now consider Using the identities obtained in appendix (C) for a radiation dominated cosmology, we write where we introduced the definitions and γ i is the local Lorentz factor at time η i .
In terms of these variables we find We note that the fulfillment of the adiabatic condition at all times implies that For ξ ′ ≪ 1 it is straightforward to find that J[ξ ′ ] features the expansion In terms of these variables we find that the subtracted integral I 2b (k, η, η b ) defined by eqn. (V.4) is given by Consider the two contributions to this function, During the time scale when J(ξ ′ ) ≪ 1 the term cos[J(ξ ′ )] ≃ 1 therefore F 2 (ξ) ≃ F 1 (ξ) and I 2b ≃ 0. Figures (6, 7) display F 1,2 (ξ) for ω i η i = 100 and γ i = 2, 10 respectively for ξ b = 1/ω i η i . Although in general the value of ξ m must be obtained numerically, for ω i η i ≫ 1 there are two limits that afford an analytic estimate: a:) for ω i η i ≫ 1 and γ i ≃ 1, we assume, self consistently that ξ m ≪ 1, therefore from eqn. (D.11) we obtain for γ i ≃ 1 , (D. 15) this expression confirms the assumption that ξ m ≪ 1 for γ i ≃ 1. b:) for γ i ≫ 1, it is convenient to carry out the integral (D.9) by expanding ω φ k (η 1 ) ≃ k + m 2 φ C 2 (η 1 )/k + · · · and keeping the leading order term, we find This latter expression is fairly accurate even for γ i ≃ 2, 3. We have confirmed numerically the validity of these approximate values of the maxima of F 2 (ξ) (see figures (6,7)). In both cases we find that for ω i η i ≫ 1 the value at the maxima fulfill ξ m /γ i ≪ 1. In summary, we find that during the time interval ξ b < ξ < ξ m F 1 (ξ) ≃ F 2 (ξ) ≃ 2 ln[ξ/ξ b ] and I 2b [k, η, η b ] ≃ 0. For ξ > ξ m the contribution F 2 (ξ) ≃ F 2 (ξ m ) ≃ 2 ln[ξ m /ξ b ] remaining constant while F 1 (ξ) increases monotonically. The above analysis shows that for ω i η i ≫ 1 it follows that ξ m ≪ γ i in the whole range of γ i , therefore during the interval ξ m < ξ < γ i and W [ξ ′ ] ≃ 1, hence The behavior of F 1,2 in the ultrarelativistic case γ i ≫ 1 is summarized as follows: Let us now consider the following integral in I 3b , namely the second contribution to eqn. (IV.8): For η ′ ≫ η i and z(η ′ ) ≫ 1 the integral in (E.10) has an adiabatic expansion, for τ ≪ z we find δ k (τ ; η ′ ) = − τ 2 2γ 2 z + · · · , (E.13) therefore δ k is of adiabatic order one and higher, furthermore, R[τ ; z] = 1 − τ z γ 2 + · · · (E.14) and to leading (zeroth) adiabatic order we can replace P[τ ; η ′ ] = 1. The τ integral in (E.10) is dominated by the region τ ≃ 0, and the region for which τ ≃ z yields a contribution ∝ 1/z, hence of first adiabatic order or smaller. Therefore to leading (zeroth) adiabatic order we neglect δ k in (E.6) and replace P → 1 in (E.10).
Although the variables (E.2,E.3) allow an explicit identification of the nature of the adiabatic expansion, the most suitable variables to merge the results for I 3b with those of the contributions from I 2b are those introduced in appendices (C) and (D). We now recast the results for I 3b in terms of these variables. Introducing 15) in terms of which we find using (C.1) which fulfill the identity We have argued above in this appendix that for ω φ k (η ′ ) η ′ ≫ 1 the term δ ≡ ∆ is higher order in the adiabatic approximation and can be neglected, and that to leading order in this approximation we can set P ≡ P → 1. We now test this assertion numerically in terms of the new variables. Since in the new variables the product ω φ k (η ′ ) η ′ = ω i η i γ i y f (y) it follows that ω φ k (η ′ ) η ′ ≫ 1 at all times implies that ω i η i ≫ 1 which is precisely the statement of the validity of the adiabatic approximation at the initial time. beyond which S(ξ) ≃ 1.
In particular for γ i = 1 (the particle decaying at rest) and ω i η i ≫ 1, S(ξ ′ ) reaches a maximum at ξ ′ = ξ s ≃ π/ω i η i + O(1/(ω i η i ) 2 ) with S(ξ s ) ≃ 1. In the opposite limit, for an ultrarelativistic particle with γ i ≫ 1 and ω i η i ≫ 1 we find self-consistently that S(ξ ′ ) reaches a maximum at ξ s with α(ξ s ) = π, where ξ s is a solution of For 2π γ 2 i /ω i η i ≪ 1 we find and for 2π γ 2 i /ω i η i ≫ 1 (E. 30) In both cases we find that ξs γ i ≪ 1 whenever γ i ≫ 1. Using the relations derived in appendix (C) along with the identities C(η ′ ) = C(η i ) (η ′ /η i ) = C(η i ) (1 + ξ ′ ) and m φ C(η i ) = ω i /γ i , it follows that I 3 (k, η) given by (E.11) can be written in terms of the same variables as I 2 , namely ξ = η/η i − 1 and η b = η b /η i − 1. We find The contribution from the term with S in the integrand must be done numerically, however, the first term can be done analytically, yielding where we have set ξ b = 0 to leading adiabatic order.