Threshold and infrared singularities: time evolution, asymptotic state and entanglement entropy

Threshold and infrared divergences are studied as possible mechanisms of particle production and compared to the usual decay process in a model quantum field theory from which generalizations are obtained. A spectral representation of the propagator of the decaying particle suggests that decay, threshold and infrared singularities while seemingly different phenomena are qualitatively related. We implement a non-perturbative dynamical resummation method to study the time evolution of an initial state. It is manifestly unitary and yields the asymptotic state and the distribution function of produced particles. Whereas the survival probability in a decay process falls off as $e^{-\Gamma t}$, for threshold and infrared divergent cases falls off instead as $e^{-\sqrt{t/t^*}}$ and $t^{-\Delta}$ respectively, with $\Gamma, \Delta \propto (coupling)^2$ whereas $1/t^* \propto (coupling)^4$. Despite the different decay dynamics, the asymptotic state is qualitatively similar: a kinematically entangled state of the daughter particles with a distribution function which fulfills the unitarity condition and is strongly peaked at energy conserving transitions but broadened by the"lifetime"$1/\Gamma~;~ t^*$ for usual decay and threshold singularity, whereas it scales with the anomalous dimension $\Delta$ for the infrared singular case. Threshold and infrared instabilities are production mechanisms just as efficient as particle decay. If one of the particles is in a dark sector and not observed, the loss of information yields an entanglement entropy determined by the distribution functions and increases upon unitary time evolution.

generalizations are obtained. A spectral representation of the propagator of the decaying particle suggests that decay, threshold and infrared singularities while seemingly different phenomena are qualitatively related. We implement a non-perturbative dynamical resummation method to study the time evolution of an initial state. It is manifestly unitary and yields the asymptotic state and the distribution function of produced particles. Whereas the survival probability in a decay process falls off as e −Γt , for threshold and infrared divergent cases it falls off instead as e − √ t/t * and t −∆ respectively, with Γ, ∆ ∝ (coupling) 2 whereas 1/t * ∝ (coupling) 4 . Despite the different decay dynamics, the asymptotic state is qualitatively similar: a kinematically entangled state of the daughter particles with a distribution function which fulfills the unitarity condition and is strongly peaked at energy conserving transitions but broadened by the "lifetime" 1/Γ ; t * for usual decay and threshold singularity, whereas it scales with the anomalous dimension ∆ for the infrared singular case. Threshold and infrared instabilities are production mechanisms just as efficient as particle decay. If one of the particles is in a dark sector and not observed, the loss of information yields an entanglement entropy determined by the distribution functions and increases upon unitary time evolution.

I. INTRODUCTION
Most particles in the standard model decay, quarks and gluons are confined, and charged particles interacting with gauge fields are dressed by a cloud of soft massless gauge fields. Therefore, of all the particles in the standard model perhaps only neutrinos and photons appear as asymptotic single particle states in the S-matrix. The dressing of charged particles by massless gauge bosons results in infrared divergences in radiative corrections as a consequence of the emission and absorption of the soft gauge quanta. Understanding these infrared phenomena and the infrared finiteness of the S-matrix has been [1][2][3][4][5][6][7][8][9] and continues [10][11][12][12][13][14][15][16][17] to be the focus of a substantial body of work motivated by precision calculations of physical observables for collider experiments [18][19][20].
Infrared phenomena also plays a fundamental role in quantum aspects of gravity as a consequence of emission and absorption of gravitons [21,22].
Prior to the discovery of the Higgs boson, early work [23,24] recognized that the S-matrix approach to describing particle decay breaks down when the mass of the particle approaches the multiparticle threshold [23][24][25][26]. In particular refs. [23][24][25][26] recognized a singularity in the self-energy of the particle as its mass approaches threshold from below, and as a consequence the particle no longer appears as an asymptotic state in the S-matrix.
Notably this situation is similar to the case of infrared singularities in gauge theories that arise because the mass of the charged particle coincides with the multiparticle threshold suggesting that, perhaps, threshold and infrared singularities, although quantitatively different, are manifestations of similar phenomena suggesting a generalized decay of the particle.

Motivations and objectives:
Extensions beyond the standard model posit the existence of new particles as possible explanations of the origin of dark matter in cosmology. Some of these extensions introduce light or ultralight particles [27][28][29][30][31][32][33], and an important question in these models is to identify and assess the production mechanism for these dark matter candidates. A recent study [34] revealed certain universality of infrared phenomena in the sense that infrared divergences associated with emission and absorption of massless quanta feature similar dynamics and asymptotic states in bosonic, fermionic and (abelian) gauge theories. This study also revealed that the infrared divergences could be an effective production mechanism of soft massless particles, and was extrapolated to the realm of production of light dark matter or dark radiation during a radiation dominated cosmology [35]. In this article we extend the study of ref. [34] to compare and contrast the dynamics of decay, threshold and infrared divergences to identify hitherto unexplored production mechanisms that could be relevant in early Universe cosmology and also, perhaps, of some phenomenological interest in particle physics.
As windows beyond the standard model open to explore possible explanations of dark matter and or dark radiation, our study is motivated by its possible impact in identifying and assessing alternative production mechanisms available in the dark sector, but also to explore fundamental aspects of the dynamics of particle decay, threshold and infrared divergences that could be of a more overarching phenomenological and theoretical interest.
Objectives: Our objectives in this study are the following: i) to compare and contrast the dynamical aspects of particle decay and threshold and infrared divergences within a model quantum field theory and draw more general conclusions on the time evolution of initial towards asymptotic states, ii) to understand threshold and infrared singularities as possible production mechanisms and to explore a qualitative similarity between these seemingly different phenomena, iii) to understand the time evolution that leads from the initial to the final asymptotic state and to characterize the properties of the latter, iv) for threshold and infrared divergences the usual decay rates vanish, therefore understanding the time evolution of initial states will clarify the dynamics of relaxation towards equilibrium in these cases.
Our study does not address the important issues of the infrared finiteness of the S-matrix, a far broader subject of much current interest [17][18][19][20]. It is much more narrowly focused on understanding the time evolution of states and the emerging asymptotic states in the case of threshold and infrared divergences. A reassessment [36] of the Lehmann, Symanzik and Zimmermann reduction formula for asymptotic states beginning with a finite time analysis and extending it to the infinite time limit has highlighted the subtleties of this limit.
Our study in this article may provide complementary further insights into asymptotic theory in cases in which threshold and infrared divergences substantially modify the asymptotic long time dynamics, and may contribute to the fundamental understanding of the asymptotic states emerging from these processes.
Brief summary of results: We study decay, threshold and infrared phenomena within a simple model of a real scalar field Φ coupled to two other scalar fields of different masses that effectively captures the different phenomena by varying the various masses. The Kallen-Lehmann representation of the propagator of the Φ field including radiative corrections illustrates how decay, threshold and infrared phenomena, although seemingly disparate are qualitatively related.
Furthermore, it clearly shows the breakdown of a Breit-Wigner approximation as the mass of the particle approaches threshold.
A dynamical resummation method [34,37] is implemented to study the time evolution of an initial single particle state of the Φ field towards the final asymptotic state in all cases. This method is manifestly unitary and complementary to the dynamical renormalization group [38,39].
It yields not only the time evolution of the initial state, but also describes the emergence of the asymptotic state during the evolution and its properties.
We find that whereas the time evolution of the survival probability of a single particle state in a typical decay process is e −Γt , in the cases of threshold and infrared singularities the usual decay rate vanishes, however we find that the survival probability of the initial state indeed decays: in the case of threshold divergence it evolves as e − √ t/t * and for infrared divergences as ∝ t −∆ . Whereas Γ and the anomalous dimension ∆ are of O(g 2 ) with g the coupling, the relaxation time scale t * ∝ 1/g 4 as a consequence of the threshold singularity.
We find that despite the different time evolution, the asymptotic state is qualitatively similar: a kinematically entangled state of the daughter particles with pair correlations. We obtain the probabilities of these pairs, show that they satisfy the unitarity condition and identify them as the distribution function of the produced particles which are obtained in each case. Although these are peaked at energy conserving transitions, are much narrower in the case of threshold divergences as a consequence of a longer "lifetime" of the initial state and feature a scaling behavior with the anomalous dimension ∆ in the case of infrared divergences.
A corollary of this result is that threshold and infrared singularities are just as efficient production mechanisms as decay.
These asymptotic states are very different from those postulated in quantum electrodynamics [3,5,8,16] as solutions to the infrared problem, but are unambiguously obtained from the unitary time evolution of an initial state. We argue that the pair correlations in the asymptotic state, in other words the entanglement of the daughter particles, implies the same distribution function for each, which we obtain from the time evolution in all cases. If either one of the daughter particles is not measured for example in the "invisible decay" into a dark matter particle, the information loss leads to an entanglement entropy, which is shown to grow during the time evolution from the initial to the asymptotic state.

II. KALLEN-LEHMANN SPECTRAL REPRESENTATION:
We consider a model of a massive real scalar field Φ coupled to two other real scalar fields χ 1 , χ 2 , to illustrate the main phenomena within a simpler setting, with the objective of drawing more general conclusions. Such model has previously been investigated within the context of threshold singularities in refs. [24,26]. It is described by the following Lagrangian density This Lagrangian density provides a simple arena to study the main aspects of our focus in this article: i) if M < (m 1 + m 2 ), a single Φ particle is stable, ii) when M > (m 1 + m 2 ) a Φ particle is unstable and decays into a pair of χ 1 , χ 2 particles, iii) when M = (m 1 + m 2 ) the mass of the Φ particle is exactly at threshold and this case is a manifestation of the threshold singularity, studied originally in ref. [24], iv) infrared singularity when M = m 1 , m 2 = 0 arising from the emission and absorption of massless quanta. In this case again the mass of the particle Φ coincides with the multiparticle threshold. This latter case features the same infrared singularities as that of a charged field coupled to a massless field studied in ref. [34] within a bosonic model with Lagrangian In this model the infrared singularity emerges in the self-energy of the Φ field as a consequence of the emission and absorption of massless quanta Φ ↔ Φχ. In ref. [34] it is shown that the infrared behavior of this model is similar to that of a Dirac fermion Yukawa coupled to a massless scalar (a renormalizable theory), and in turn is similar to the infrared divergence of the fermionic self-energy in quantum electrodynamics. Hence, the Lagrangian (II.1) furnishes a simple quantum field theory that allows to study all four cases: i) stable, ii) unstable, iii) threshold singularity and iv) infrared divergence within the same model by adjusting the masses appropriately. The Lagrangian density (II.1) describes a super renormalizable theory, however, because we are interested in infrared and long time phenomena which we expect to be insensitive to the ultraviolet behavior of the theory, this model is expected to capture the long time dynamics reliably. This expectation is confirmed by the study of ref. [34] where infrared phenomena and long time dynamics were shown to be the same for a superrenormalizable and a renormalizable model. Furthermore, in ref. [37] it has been shown that ultraviolet divergences contribute to very early transients that do not affect the long time dynamics and can be safely absorbed into a renormalization of the initial amplitude. This is a consequence of the wide separation of time scales between the ultraviolet early transients and the long time infrared phenomena. Taken together the results of these previous studies serve as anchors that allow us to draw more general conclusions on the long time dynamics from the simple model described by eqn. (II.1).
We begin by studying the Kallen-Lehmann spectral representation [40] of the single Φ particle propagator including a Dyson resummation of the one loop self-energy shown in fig. (2). The propagator is given by the self-energy is calculated in dimensional regularization in dimension D = 4 − ε, and introducing a renormalization scale µ we find with γ E the Euler-Mascheroni constant and Separating the real and imaginary parts of the self-energy, we find Subtracting the real part of the self-energy at P 2 = M 2 p at which the real part of the inverse propagator vanishes, namely and to leading order in the coupling replacing M → M p in the expression for the real and imaginary parts of the self-energy (II.8), it follows that The Kallen-Lehmann spectral function is given by [40] σ( it contains the information on the asymptotic properties of the quanta of the real scalar field Φ and obeys the sum rule σ(P 2 ) dP 2 = 1 . (II.14) A single particle pole below threshold, namely for P 2 < (m 1 + m 2 ) 2 , for which Σ I (P 2 ) = 0, (II. 16) The wave function renormalization constant yields the amplitude of the single particle pole and determines the overlap between the bare single particle state and the asymptotic renormalized state of a stable particle that has been dressed by quantum fluctuations. Therefore, when the single particle pole is below threshold, the particle is stable and where σ c (P 2 ) is the contribution from the multiparticle continuum above threshold. In this case the sum rule (II.14) yields whereas if the particle decays and the single particle pole is embedded in the continuum, there is no single particle pole below threshold and the sum rule (II.14) yields namely it is saturated by the continuum "background" and the single particle quanta of the Φ field are not asymptotic states.

III. DECAY, THRESHOLD AND INFRARED SINGULARITIES:
A. Decay and threshold singularities In order to discuss both cases of decay and threshold singularities we consider the simpler case of equal masses m 1 = m 2 ≡ m when the two particle threshold is at P 2 = 4m 2 , furthermore, it is convenient to introduce the dimensionless variables In terms of these variables the spectral density becomes We study the cases r > 1 (stable particle) and r < 1 (unstable decaying particle) separately to highlight both differences and similarities.
Case I: Stable particle 4m 2 > M 2 p (r > 1) In this case the propagator features an isolated single particle pole at s = 1 below the two particle threshold at r > 1 and for s < r The full spectral density in this case when the particle pole is below the two particle threshold is given by For r > 1 the particle is present as an asymptotic state with probability Z(r) < 1. However, we note that as M 2 p → 4m 2 , namely as the position of the single particle pole approaches the threshold from below, or r → 1 from above, the residue at the isolated pole below threshold vanishes as with a square root singularity, and very sharply in weak coupling, obviously this behavior is strongly non-perturbative. Furthermore, we find that while the continuum contribution to the spectral density vanishes at threshold, it becomes sharply peaked near threshold as r → 1 from above (or M 2 p → 4m 2 from below). Fig. (3) displays the spectral density for r > 1, namely the case of a stable particle described by an isolated pole below threshold.
Defining the contribution from the two particle continuum above threshold as we have confirmed numerically that the sum rule (II.18) is fulfilled. As r → 1 from above, the residue at the pole vanishes, but the continuum contribution saturates the sum rule. Fig. (4) shows Z(r) and C(r), it clearly displays that Z(r) vanishes sharply and C(r) rises sharply as M 2 p → 4m 2 from below (r → 1 + ), in agreement with the sum rule (III.8) which can be confirmed from the figure.
Precisely at M 2 p = 4m 2 when the mass shell coincides with the multiparticle threshold there is a singularity in the sense that the amplitude of the single particle pole vanishes and the spectral density at P 2 = 4m 2 diverges in such a way as to maintain the sum rule. This behaviour has been described as a threshold singularity [24]. What is clear in the case when M 2 p = 4m 2 is that the single particle "dissolves" into the continuum and is not an asymptotic state since its residue, namely the overlap of the bare and asymptotic state vanishes. However, the particle does not "decay" in the usual manner because the imaginary part of the self-energy vanishes at P 2 = M 2 p = 4m 2 , hence the "decay rate" Γ ∝ Σ I (P 2 = M 2 p )/M p vanishes identically when M 2 p = 4m 2 .
Case II: Unstable particle M 2 p > 4m 2 (r < 1) In this case the particle "pole" moves off the physical sheet into the second (or higher) Riemann sheet becoming a decaying resonant state which is not an asymptotic state in the S-matrix. The spectral density only has support above the two particle threshold where D(s, r) is given by eqn. ( moderately large coupling to exhibit the behavior as the position of the resonance approaches the threshold from above as compared with the cases where it is far above threshold. for an unstable, decaying particle with M 2 P > 4m 2 (r < 1) for r = 0.3, 0.6, 0.96; g = 0.05.
In this case the sum rule (III.9) is saturated by the contribution above threshold since there is no support below threshold, and we have confirmed numerically in all cases that C(r) = 1 with C(r) given by eqn. (III.8).
When the distance between threshold and the position of the resonance ("pole") is much larger than the width Γ the propagator and the spectral density may be very well approximated by a Breit-Wigner Lorentzian function in the narrow width approximation, where the wave function renormalization Z bw is given by (II.16) but now above threshold, with M 2 P > 4m 2 and given by , (III. 13) we note that in contrast to the case when M 2 P < 4m 2 , in this case as M 2 P → 4m 2 from above it is straightforward to confirm that Z bw remains finite in agreement with the conclusion in ref. [24]. However, when the particle is unstable Z bw does not have the interpretation of the amplitude of the renormalized single particle state in the asymptotic state. It is clear from eqn. (III.11) that the Breit-Wigner approximation of the spectral density is only reliable for very weak coupling as it does not obey the sum rule C(r) = 1 since Z bw = 1.
When the position of the resonance ("pole") is far away from threshold and for a narrow width, the propagator may be approximated by a Breit-Wigner distribution which in the narrow width approximation becomes with Z bw given by eqn. (III.13) and to leading order in the weak coupling g, we can set Z bw = 1 and recognize Γ as the decay rate at rest obtained from the lowest order S-matrix approach, namely The long time dynamics of the retarded propagator is obtained from the Fourier transform of the Breit-Wigner propagator (III.14), namely from which we interpret Z bw not as the amplitude of the single particle in the asymptotic state but as the weight of the resonance contribution to the spectral density as evident from eqn. (III.11).
As fig.(5) clearly shows, when M 2 p → 4m 2 from above (r → 1 from below), although Γ → 0, the spectral density can no longer be described as a narrow width Breit-Wigner Lorentzian and the resonance cannot be described as a complex pole in the second (or higher) Riemann sheet. In this limit we find that as s → 1 + and for weak coupling displaying a square root singularity at threshold in agreement with fig. (5). In this case the narrow width Breit-Wigner approximation is neither valid nor useful to describe the resonance near threshold as the spectral density diverges as the threshold is approached. However, this singularity notwithstanding, we confirmed numerically the sum rule C(1) = 1, but the interpretation of a finite Z bw as a wave function renormalization associated with the resonance is no longer useful as a description of the asymptotic state.
As M 2 p → 4m 2 from below the single particle state is no longer an asymptotic state, however its amplitude does not decay in time with the usual exponential decay law because the decay rate Γ → 0 when the "pole" coincides with the two particle threshold. Furthermore, as is clear from fig. (5) and from eqn. (III.19) the Breit-Wigner approximation breaks down as M p → 4m 2 from above and the time evolution of the resonant state is not an exponential as in the case (III.18) but a more complicated function determined by the square root branch cut beginning at threshold. This time evolution will be studied in detail in the next section.

B. Infrared singularity:
The infrared singularity is associated with the emission and absorption of a massless particle by a massive one. Such is the case, for example, in quantum electrodynamics where the one loop fermion self energy features an infared divergence on the fermion mass shell. Whereas in gauge theories care must be taken to maintain gauge invariance and satisfy Ward identities, the simpler model of a charged scalar field in interaction with a neutral massless field, described by the Lagrangian density (II.2) features the same infrared divergence [34]. In turn, we can study the infrared divergence within the framework of the model described by (II.1) by taking m 1 = M, m 2 = 0. However, in order to display the emergence of the infrared singularity more clearly, let us consider the case m 2 = 0 ; m 1 = m and we will explore the limit M p → m where the infrared divergence becomes manifest.
In this case, in terms of the variables s; g introduced in eqn. (III.1) along with the ratio the spectral density is given by For R > 1 the Φ particle is stable and we find and M 2 p σ c (s) is the contribution from the two particle continuum above threshold, given by eqn. (III.21) for s > R. Equation (III.24) clearly shows that as R → 1 + the wave function renormalization vanishes, namely there is no longer an isolated single particle pole, again the particle "dissolves" into the continuum as its (renormalized) mass approaches the threshold from below.
This behavior is displayed in Fig. (6) which shows Z ir (R) vs. R for g = 0.01.
Again we have confirmed numerically the validity of the sum rule Z ir (R) + ∞ R M 2 p σ c (s)ds = 1, therefore when the single particle pole approaches the threshold from below, the sum rule is saturated from the continuum contribution M 2 p σ c (s) which is displayed in fig. (7) showing a sharp rise near threshold as it absorbs the normalization of the single particle pole when it merges with threshold.
We conclude that in the infrared limit R = 1 when the single particle pole merges with the threshold the single particle "dissolves" into the continuum and is no longer an asymptotic state.
However, as in the case of threshold singularity, the particle does not "decay" in the usual sense because the decay rate vanishes when the pole mass coincides with the two particle threshold. The infrared singularity for the case R = 1 is the same as that studied in the model of a charged scalar field coupled to a massless real scalar field [34], this previous study also revealed an emerging universality of infrared phenomena and showed that the amplitude of an initial single particle state decays with a power law with anomalous dimension. The infrared singularities in this model field theory are the same as in the general Lagrangian density (II.1), replacing χ 1 → Φ ; χ 2 → χ and The lesson that we draw from this analysis based on the Kallen-Lehmann representation is that threshold and infrared divergences result in that the probability of the single particle state vanishes, transferring the normalization to the multiparticle continuum. This "flow" of probability from single particle to multiparticle states is a manifestation of particle production, we refer to these cases as generalized decay, because, indeed the single Φ particle does decay into the multiparticle continuum despite the fact that the S-matrix decay rate Γ formally vanishes, because the imaginary part of the self energy vanishes at threshold. Furthermore, although there are quantitative differences between threshold and infrared divergences, for example in the manner that Z vanishes as the threshold is approached and the sharp rise of the continuum contribution near threshold, qualitatively the two phenomena are rather similar as evidenced by the figures displaying Z and σ c in both cases.
In the next section we study this "flow" or generalized decay from the point of view of the time evolution of an initial single particle state towards an asymptotic state.

IV. TIME EVOLUTION: DYNAMICAL RESUMMATION METHOD :
We now obtain the asymptotic state by following the time evolution of an initial single Φ particle state. For this purpose we now introduce a method that implements a dynamical resummation directly in time [34,37] and is complementary to the dynamical renormalization group [38,39]. We briefly revisit here the main aspects of this method for coherence and completeness of presentation, referring the reader to previous studies [34,37] for more details.
Consider a system whose Hamiltonian is H = H 0 + H I with H I a perturbation. The time evolution of states in the interaction picture of H 0 is given by where the interaction Hamiltonian in the interaction picture is The Schroedinger eqn. (IV.1) has the formal solution and the time evolution operator in the interaction picture U (t, t 0 ) obeys Now we can expand the time evolved state as where |n are eigenstates of the unperturbed Hamiltonian, H 0 |n = E n |n , and form a complete set of orthonormal many particle states. From eq.(IV.1) one finds the exact equation of motion for the coefficients C n (t), namelyĊ Although this equation is exact, it generates an infinite hierarchy of simultaneous equations when the Hilbert space of states spanned by {|n } is infinite dimensional. However, this hierarchy can be truncated by considering the transition between states connected by the interaction Hamiltonian at a given order in H I .
Specifically, for the model under consideration here, consider the situation depicted in figure (8) where the single particle state, |1 Φ k , couples to the two particle state |1 χ 1 p ; 1 χ 2 q , which couples back to |1 Φ k via the interaction Hamiltonian where the fields are in the interaction picture.
Consider that at the initial time t = 0 a single Φ particle state with momentum k is prepared, upon time evolution the interaction Hamiltonian connects this state with a two particle state of the χ 1 , χ 2 fields, therefore the time evolved state is given by where the dots stand for multiparticle states that connect to |1 Φ k in higher order in H I , and we have explicitly used momentum conservation which is justified by the matrix elements obtained in appendix (A). In what follows we use q = k − p to simplify notation.
The hierarchy of equations (IV.6) lead to the following coupled equations for the amplitudeṡ The initial value problem in which at time t = 0 the initial state is a single Φ particle state, namely |Ψ(t = 0) = |1 Φ k and the vacuum for the other fields corresponds to the initial conditions We solve eq.(IV.10) with these initial conditions and input the solution into eq.(IV.9) to find and we used eqn. (IV.2). It is convenient to write Σ Φ (t, t ′ ) in a spectral representation, namely where we have introduced the spectral density (t) from which we can compute the time dependent probability to populate the two particle The hermiticity of the interaction Hamiltonian H I , and the equations (IV.9,IV.10) yield which together with the initial conditions in eqs.(IV.11) yields the unitarity relation This is the statement that the time evolution operator U (t, 0) is unitary, namely The integro-differential equation (IV.13) can be solved exactly via Laplace transform [37], however finding the time evolution from the inverse transform involves a technically difficult integral with branch cut singularities. Instead, recognizing that for weak coupling there is a separation of time scales we invoke the dynamical resummation method introduced in ref. [34] which hinges on a separation of time scales warranted for weak coupling, and provides a non-perturbative resummation directly in real time equivalent to the dynamical renormalization group [38,39].
The time evolution of C Φ k (t) determined by eq.(IV.13) is slow in the sense that the time scale is determined by a weak coupling kernel Σ which is second order in the coupling. This allows us to use an approximation in terms of a consistent expansion in time derivatives of C Φ k (t). Let us Integrating by parts in eq.(IV.13) we obtain The second term on the right hand side is formally of fourth order in H I suggesting how a systematic approximation scheme can be developed. Setting and integrating by parts again, we find This process can be implemented systematically resulting in higher order differential equations.
where we used the initial condition C Φ and With this solution we find the time evolution of the coefficients of the multiparticle states from eqn. (IV.12) from which we obtain the probability of the multiparticle states in the time evolved state, and in particular the asymptotic state as t → ∞.
The survival probability of the initial state is given by (IV.32) In the long time limit we find where P stands for the principal part, yielding a renormalization of the bare frequency of the state , whereas the long time limit of γ(t) yields the decay law of the initial state.
It is illuminating to write the energy renormalization using the explicit form of the spectral density given by eqn. (IV.16), namely where the principal part in eqn.(IV.33) removes the region in momenta when the denominator vanishes denoted by the superscript prime in the sum. This is the usual quantum mechanical result for the second order energy shift.

A. Stable particles:
Before we analyze the time evolution of the coefficients C Φ k (t), we can understand their asymptotic behavior for the case of stable particles. The spectral density ρ Φ (k 0 ) vanishes for k 0 < k 0T where k 0T = k 2 + (m 1 + m 2 ) 2 corresponds to the two particle threshold (see eqn. (A.3)). In the case of a stable particle E Φ k < k 0T the denominator in γ(t), eqn. (IV.30) never vanishes and the cosine term averages out in the long time limit. Therefore in the case of a stable particle with energy below the two particle threshold energy it follows that Hence, in the case of a stable particle for which the single particle energy is below the two particle threshold (neglecting renormalization to lowest order) the time evolution of the initial single Φ particle amplitude yields the asymptotic result namely, the probability of finding the initial (bare) single particle state in the asymptotic state For the case m 1 = m 2 = m, it is straightforward to find that as M 2 → 4m 2 , namely as displaying the threshold divergence that results in the vanishing of the overlap between the asymptotic state and the initial single particle state, namely Z → 0.
As we discuss below this is a consequence of taking the infinite time limit too soon in the limit when the single particle energy approaches the two particle threshold.
It is illuminating to understand the time scale over which the integral (IV.30) approaches its asymptotic limit (IV.35). In the case of a stable particle, with E k < k 0T as mentioned above, the denominator in (IV.30) does not vanish in the domain of integration k 0 ≥ k 0T , therefore we can separate the time dependent cosine term from the expression for γ(t). Hence, consider the integral in the long time limit the cosine averages out and this integral vanishes by dephasing, on a time scale t dp ≃ 1 since k 0T − E k is the smallest frequency contributing to the integral, which, in turn dominates the long time limit. Therefore as the single particle energy approaches the threshold from below the dephasing time scale t dp diverges, and as discussed above the overlap Z vanishes. This is the case of threshold divergences, the dynamics of which will be studied in detail below.

B. Decaying particle:
In this case the (renormalized) single particle energy is above the two particle threshold, E Φ k > k 0T , and the denominator in (IV.30) vanishes within the domain of integration, therefore the cosine term cannot be separated. Let us consider the case of equal masses m 1 = m 2 = m in which case we find from eqns. (IV.30, A.3) that in terms of which where for a decaying state with M 2 > 4m 2 ⇒ E Φ k > k 0T , and (IV.44) In the long time limit M t → ∞; X(t) → ∞ we find is the correct decay rate (III.16) including the time dilation factor, and is identified with the usual result from Fermi's golden rule, and where ρ e (ξ) = (ρ(ξ) + ρ(−ξ))/2. The details of the derivation of this result are given in appendix (B). We find the long time behavior where δE Φ ∞ is a renormalization of the single particle energy.

C. Threshold singularity:
The expression for γ(t), eqn. (IV.43) in terms of ρ given by eqn. (IV.44) makes explicit the modification of the decay in the case of threshold singularity, namely M 2 = 4m 2 , in this case, we note that in this case the decay rate (IV.46) Γ k = 0, and the spectral density vanishes with a square root at threshold.
In the limit t → ∞ it follows that ρ → and the survival probability decays as with an effective lifetime where E Φ k /M is the usual time dilation factor. Namely the decay law changes from e −Γ k t → e − √ t/t * k , as the mass of the particle approaches the threshold. The square root behavior is a consequence of fact that the spectral density vanishes as a square root near threshold.
Furthermore, whereas in the case of decay there is a constant, time independent contribution in the asymptotic long time limit of γ(t) which defines the wave function renormalization, no such term arises in the case of threshold singularity. Therefore, at threshold when M 2 = 4m 2 , even when the usual decay rate (III.16) vanishes, the amplitude of the initial single particle state decays not as e −Γt but as e − √ t/t * with an effective lifetime t * given by eqn. (IV.52), reflecting the square root divergence in both in the spectral density approaching threshold from above and the wave function renormalization approaching the threshold from below. This new decay law is in qualitative agreement with a result found in ref. [26] and implies that the single particle state is not an asymptotic state in agreement with a vanishing wave function renormalization from below, and the fact that the continuum contribution of the spectral density saturates the sum rule.
In order to understand the asymptotic behavior in more detail it proves convenient to study the case k = 0 and to introduce the dimensionless combinations r = 4m 2 /M 2 and τ = M t, in terms of which, for k = 0, Since the factor (1 − cos(x))/x 2 is localized within a region of width ≃ 2π around the origin, for r < 1 the function J(r, τ ) approaches its asymptotic limit J(r, ∞) = π (1 − r) 1/2 within a τ scale ≃ 2π/(1 − r). As the threshold is approached from above, namely r → 1 from below, the asymptotic value becomes smaller and smaller taking a longer and longer time scale to reach, and for r = 1 the function J(1, τ ) ∝ 1/ √ τ for large τ . This behavior is clearly displayed in fig. (9) for r = 0.8, 0.9, 0.98, 1.
vanishing linearly as k 0 approaches the threshold k 0T = E Φ k . This situation must be contrasted with the case of threshold divergence where the spectral density vanishes as a square root as k 0 approaches threshold. However, in both cases the usual decay rate Γ given by eqn. (IV.46) vanishes.
It is convenient to introduce the dimensionless combinations in terms of which we find in this case We note that this integral features a logarithimic divergence in the region of small η. Following ref. [34] we write the above integral as and (IV.60) In the long time limit T ≫ 1 the first integral in I 1 (T ) features an infrared logarithmic divergence, whereas in the second integral and in I 2 (T ) the cosine terms average out yielding in this limit and γ E is the Euler-Mascheroni constant, z ir is infrared and ultraviolet finite. Therefore, for the infrared case we find namely the probability of the initial single particle state decays in time as a power law with anomalous dimension 2g 2 and is not an asymptotic state in S-matrix amplitudes.
In summary, we find the following asymptotic long time limits for the unstable cases in which the "mass shell" of the particle is above or at threshold, In all cases of "generalized decay" as described by eqn. (IV.64) the asymptotic state is The probabilities must obey the sum rule which is the statement of unitarity (IV.18) in the asymptotic long time limit when the amplitude of the initial state vanishes.
The question that we address is how this sum rule is fulfilled being that the probabilities |C χ p; k (∞)| 2 are formally of order g 2 . In appendix (C) we show that for all cases, and up to O(H 4 I ) the asymptotic probabilities are given by where E Φ k in this expression is the renormalized single particle energy (see eqn. (C.8)). Introducing the spectral representation (IV. 16 Armed with this general expression we can now study the individual cases by considering the different forms of γ(t) and spectral densities.
We can confirm the unitarity relation (V.2) to leading (zeroth) order at this stage by setting Z d = 1 and writing In the narrow width limit where we used the result (IV.46). To prove unitarity up to O(H 4 I ) requires a somewhat deeper analysis, which we now undertake.
Following similar steps as in appendix (B) we find In the narrow width limit ε ≪ 1 the first integral is straightforward yielding π/ε − 2/ξ + O(ε), and we can set ε → 0 in the second and third integrals 1 , yielding where we used the results (B.12,B.13). Therefore, with Z d = e −2z d ≃ 1 − 2z d + · · · we indeed find  Infrared divergence: The fulfillment of the unitarity condition (V.2) in the case of infrared divergence has been confirmed up to O(H 4 I ) in ref. [34] to which the reader is referred for further technical details. However, for the sake of completeness we here summarize the main steps to leading (zeroth) order to compare with the previous cases. In this case, the spectral density is .

B. The asymptotic state:
In the cases of decay, threshold and infrared singularities discussed above, the asymptotic state after the amplitude of the initial state has become negligible, features a common form, namely 23) or in the case of the infrared singularity for the model given by the Lagrangian density (II.2) of a charged scalar field interacting with a massless scalar, obtained by identifying χ 1 ≡ Φ ; χ 2 ≡ χ and χ a massless field, namely In the analysis below we will consider the asymptotic state (V.23) describing all cases with the implicit understanding that the case of infrared divergence is obtained by the replacement In all the cases studied in this article, namely particle decay and those that feature threshold and infrared singularities, the initial single particle state "decays" either exponentially or with a power law and is not an asymptotic state. The asymptotic states that result from the time evolution of these processes are given by (V.23,V.24), these are correlated kinematically entangled states of the daughter particles. We highlight this noteworthy point: particle decay and the processes that feature threshold and infrared divergences, while quantitatively different in the details of the dynamical evolution of the amplitudes, are asymptotically qualitatively similar and determined by the asymptotic states (V.23,V.24) which characterize the production of the daughter particles, with the total production probability fulfilling the unitarity relation, namely p |C χ p; k (∞)| 2 = 1.
Hence, in conclusion, threshold and infrared singularities result in the production of the daughter particles, much in the same manner as the usual decay process.
Out of the pure state (V.23) (or (V.24)), we can construct the (pure state) density matrix  26) or similarly, of operators that act solely on the Hilbert space of the field χ 2 . In these cases the trace over the "unobserved" fields yields a reduced density matrix, namely ̺ χ 1 = Tr χ 2 ̺ ; ̺ χ 2 = Tr χ 1 ̺ , (V. 27) and the unitarity condition obviously yields From the asymptotic state (V.23) we find These reduced density matrices describe mixed states, and are diagonal in momentum and particle number. In particular we identify the distribution function of the produced particles in terms of the expectation value of the number operator for each field N χ ( p) as therefore, as a consequence of entanglement, both daughter particles share the same distribution function. Furthermore, as a consequence of unitarity we find a result with the clear interpretation that there are in total one χ 1 and one χ 2 particles in the asymptotic state, a physically correct outcome of the generalized decay of a single Φ particle into one χ 1 and one χ 2 particles.
which can be written as the replacement (V.8) in the narrow width limit yields a sharp energy conserving delta function, however, a small but finite width introduces an energy uncertainty in the distribution of daughter particles as a consequence of the lifetime 1/Γ k of the initial state with a concomitant broadening of the distribution function.
For the case of infrared singularity, namely m 1 = M, m 2 = 0 we find This distribution function does not feature any particular scale, although it is peaked at Ω = namely energy conserving transitions, it falls off as a power law Ω −2(1−g 2 ) consistently with the scale invariance and anomalous dimension associated with infrared phenomena found in ref. [34] .
For the case of threshold singularity we can write the distribution function as although there is an analytic expression for this integral in terms of Fresnel integral functions, a graphical representation is more illuminating and is displayed in fig. (10). The distribution function is sharply peaked at Ω = E χ 1 p + E χ 2 q − E Φ k ≃ 0 with a width of the order of 1/t * consistent with the "lifetime" of the initial state. Note that this distribution is narrower than the case of decay because the lifetime t * ∝ 1/g 4 is longer as compared with 1/Γ ∝ 1/g 2 . The entanglement entropy has also been discussed in the case of decay in ref. [41], and for infrared singularity within the context of quantum electrodynamics in refs. [16,34,35]. While the entanglement entropy is a corollary of the pair correlation in the asymptotic state, it is just beginning to receive attention within particle physics [42].

VI. DISCUSSION
General lessons: common aspects of decay, threshold and infrared divergences: Although we have focused our study on a simple quantum field theory, the results obtained in the previous section suggest some universality in the asymptotic state arising from decay, threshold or infrared singularities in that in all these cases the asymptotic state is a kinematically entangled multiparticle state with a probability that saturates the unitarity constraint. While this is obviously a consequence of unitary time evolution, the corollary is that threshold and infrared divergences are just as efficient mechanisms of particle production as the process of decay. The asymptotic distribution functions |C χ (∞)| 2 are peaked at Ω = E χ 1 p + E χ 2 q − E Φ k = 0, namely energy conserving transitions, but broadened. In the case of decay the width of this distribution is O(Γ) ∝ g 2 consistent with a broadening by the lifetime ∝ 1/Γ, for threshold singularity the distribution is narrower, within a width of O(1/t * ) ∝ g 4 again consistent with a much longer "lifetime", and in the case of infrared singularity the distribution function features a scaling behavior with anomalous dimension Ω −2(1−g 2 ) , as a consequence of the scale invariance and anomalous dimension associated with infrared phenomena [34].
The detailed analysis of the different cases yield the following set of criteria on the spectral density ρ(k 0 ) and the mass of the particle that determines the time evolution of the survival probabilities: a: If the spectral density does not vanish at k 0 = E k where E k is the single particle energy of the "decaying field", the survival probability decays as usual ∝ e −Γt with Γ = 2πρ(E k ). This is simply the statement of Fermi's golden rule and is the S-matrix result for the decay width at leading order in the coupling.
b: If E k coincides with the multiparticle threshold and the spectral density vanishes at threshold as a square root ∝ |k 0 − k 0T | this case corresponds to a threshold singularity. The usual decay rate vanishes but the survival probability decays as e − √ t/t * . This result cannot be obtained within the S-matrix approach, since the transition probability per unit time in the infinite time limit, namely the usual decay rate calculated via S-matrix vanishes.
c: If E k coincides with the multiparticle threshold and the spectral density vanishes linearly at threshold ∝ |k 0 −k 0T | this case corresponds to an infrared singularity. The usual decay rate vanishes but the survival probability decays algebraically with an anomalous dimension t −∆ . Again, this result cannot be obtained via the usual S-matrix calculation for the transition probability per unit time in the infinite time limit, again such decay rate vanishes.
More generally, if the spectral density vanishes at threshold as |k 0 −k 0T | β the survival probability decays as e −C t 1−β with C a coupling dependent constant, β = 1 is the infrared singular case and yields a logarithmic behavior.
Thus, threshold and infrared singularities differ only on how the spectral density vanishes at threshold: if as a square root, then the decay is e − √ t/t * yielding a distribution function with a breadth ∝ 1/t * , if linearly, the decay is ∝ t −∆ yielding a distribution function with scaling dimension 2 − ∆.
Infrared and threshold divergences as production mechanism of ultralight particles: Although we have studied the dynamics associated with infrared divergences for the case in which the χ 2 particle is massless, the results apply to the case of ultralight particles proposed to be dark matter candidates, from "fuzzy" dark matter with a mass ≃ 10 −22 eV [31][32][33], to axions with a mass ≃ 10 −6 eV [27,28]. Consider that such particles are coupled to heavier one, with a mass 100MeV, the departure from threshold is 10 −14 of the value of the threshold position, this means that although the threshold is just above the single particle pole, the wave function renormalization is vanishingly small (see figures (4,6)) thus transferring the normalization to the continuum, and the time evolution -either as a square root or logarithmic-lasts for a very long time thus populating the asymptotic state with the ultralight degree of freedom. Thus infrared or threshold divergences are an efficient mechanism for production of ultralight dark matter candidates as proposed recently in ref. [35]. has also been recognized in ref. [24].
Relaxation and thermalization: Both for threshold and infrared divergences the usual decay rate vanishes, however the survival probability decays either as e − √ t/t * or as t −2g 2 , in either case the decay law cannnot be described by Fermi's golden rule or the S-matrix approach. This observation leads to the question of how Φ particles would thermalize with a bath of χ particles under the conditions of threshold or infrared divergence. In the usual Boltzmann equation, the thermalization rate is directly proportional to the decay rate obtained from Fermi's golden rule modified by spontaneous emission/absorption factors. This question is of relevance in cosmology and requires a treatment different from the Boltzmann equation which directly inputs the transition probabilities per unit time from S-matrix theory, these are precisely the relaxation rates from Fermi's golden rule which vanish for threshold or infrared divergences. A related question is how detailed balance emerges between decay and inverse decay processes, since in the usual formulation detailed balance is a consequence of explicit energy conservation and the energy conserving constraint is not exactly satisfied for threshold and infrared divergences, this is the reason that the usual decay rate vanishes in these cases. Work on these aspects is in progress and will be reported elsewhere [43].

Entanglement entropy, correlations and thermalization
We have discussed the emergence of the entanglement entropy upon tracing an "unobserved" degree of freedom out of the pure asymptotic state density matrix. Such tracing, or coarse graining, yields a mixed state, and a concomitant von Neumann entropy as a consequence of the loss of information in the coarse graining process. The pair correlations in the pure entangled state entail that the reduced density matrices ̺ χ 1 , ̺ χ 2 feature the same probabilities, (see eqns.V.29,V.30) which are identified as the distribution function of the produced particles, hence the same entanglement entropy.
A remarkable experiment reported in ref. [44]) shows that the entanglement entropy as a result of correlations in a closed system heralds thermalization. It is therefore an intriguing possibility that in the early Universe, indeed a closed system, the entanglement entropy associated with cosmological particle production from threshold or infrared divergences [35] may also herald the onset of a thermal state.
Entanglement plays a fundamental role in the determination of time reversal and CP violation in neutral meson systems [45], therefore it is a tantalizing possibility that correlations of particles produced via threshold or infrared divergences may prove to be also relevant in experimental particle physics. The potential relevance of the concept of entanglement entropy associated with information loss in the asymptotic final state, in particular if some of the decay products belong to a dark sector beyond the standard model, both in cosmology and in particle physics merits further study.
Phenomenological consequences of the lifetime for threshold divergences: The generalized decay as a consequence of threshold divergences with a survival probability that decays as e − √ t/t * implies that even when the usual decay rate vanishes (infinite lifetime), there is an intrinsic finite lifetime t * ∝ 1/g 4 . This result may have potentially relevant phenomenological implications, as the decay products of this process may feature displaced vertices with a very long but finite decay length.
Radiative corrections: moving away from threshold. The condition for threshold divergence, namely that the mass of the particle coincides exactly with the value of the lowest multiparticle threshold, will most likely not survive radiative corrections. However, such corrections will be proportional to a power of a small coupling, thus while not exactly at threshold, the departure from threshold is perturbatively small. Let us consider that upon radiative corrections the mass of the particle moves perturbatively below threshold so that (4m 2 − M 2 )/M 2 ∝ α with α a small coupling. In this case the particle has been rendered stable by radiative corrections, however, asymptotically its probability in the final state is Z ∝ e hence featuring an essential singularity in the coupling and for all intent and purpose the particle does not appear as an asymptotic state. If, on the other hand, the radiative correction moves the mass above threshold, the particle is unstable, decaying as e − √ t/t * during a time ∝ 1/α until it begins decaying as e −Γt , and is not an asymptotic state. In conclusion, radiative corrections while capable of moving the position of the mass shell away from threshold perturbatively, the probability of the particle to be present in the asymptotic state practically vanishes.
For the case of infrared divergence, for example for the model defined by the Lagrangian density (II.2) in which a massive charged particle emits and absorbs a massless χ particle, unless the mass of this particle is protected by some symmetry radiative corrections will induce a non-vanishing mass thus modifying the conclusions. However if such modification is perturbatively small, the mass shell of the charged particle will be very close to threshold and the near-threshold behavior will ensue as discussed above.

VII. CONCLUSIONS
Motivated by the possibility that a dark sector beyond the standard model could feature ultralight particles as dark matter candidates, in this article we study threshold and infrared divergences as hitherto unexplored possible production mechanisms that could be relevant in cosmology. In the case of threshold and infrared divergences the usual decay rates vanish, therefore understanding the time evolution in these cases will pave the way towards understanding the process of thermalization beyond the usual Boltzmann approach which inputs the transition rates per unit time in the infinite time limit. Our main objectives are to compare the usual decay process to the time evolution and particle production associated with threshold and infrared divergences and to understand the nature and characteristics of the asymptotic state. We study these different mechanisms in a model field theory that provides a simple arena to explore these phenomena within the same setting by varying the masses of the various fields yet allows to extract more general lessons. An analysis based on the Kallen-Lehmann representation of the particle propagator suggests that decay, threshold and infrared singularities, while seemingly widely different phenomena are qualitatively related, and also highlights the breakdown of a Breit Wigner approximation to propagators in the cases of threshold and infrared divergences. A dynamical resummation method complementary to the dynamical renormalization group is introduced to study the time evolution of initially prepared single particle states. This method is manifestly unitary and yields the asymptotic state, from which we obtain the distribution function of the produced particles. We find that whereas in a typical decay process the survival probability of the initial single particle state decays as e −Γt , in the case of threshold divergence it decays as e − √ t/t * and for the case of infrared divergence t −∆ , where Γ and ∆ are ∝ (coupling) 2 while t * ∝ 1/(coupling) 4 . Although the decay laws are strikingly different, the asymptotic state is more "universal" in the sense that it is a kinematically entangled state of the daughter particles. The probability of the asymptotic state is shown in each case to satisfy the unitarity condition. The distribution function of the particles in the asymptotic state are strongly peaked at energy conserving transitions, but in the case of the usual decay and of threshold singularity they are broadened by the lifetime of the decaying state 1/Γ, t * respectively, whereas in the case of the infrared divergence the distribution function falls off with a scaling behavior with scaling dimension 2 − ∆.
Therefore the results of this study indicate that threshold and infrared divergences are production mechanisms just as efficient as the usual particle decay. If either one of the particles in the final state is not observed as perhaps in an "invisible decay" into a dark matter particle, the information loss leads to an entanglement entropy which grows as a consequence of unitary time evolution. These alternative mechanisms may be relevant for production of particles in the dark sector in cosmology with possible phenomenological consequences in invisible decays with displaced vertices and long decay lengths, and also to novel thermalization dynamics, a possibility that merits further study. and passing to the continuum limit with 1 V p → d 3 p (2π) 3 we recognize that (IV.16) is given by which is the Lorentz invariant two body phase space, finally yielding To leading order in the coupling we replaced in I 1 (t) the integration interval is symmetric and the function (1−cos(x))/x 2 is even in x, therefore only the symmetric combination ρ e (x/M t) = (ρ(x/M t)+ρ(−x/M t))/2 contributes to this integral.