Production of massive bosons from the decay of a massless particle beam

Taking a two interacting scalar toy model with interaction term $g\phi\chi^2$, we study the production of $\chi$ particles coming from the decay of an asymptotic and highly occupied beam of $\phi$ particles. We perform a nonperturbative analysis coming from parametric resonant instabilities and investigate the possibility that massive $\chi$ particles are produced from decays of massless $\phi$ particles from the beam. Although this process is not present in a perturbative analysis, our nonperturbative approach allows it to happen under certain conditions. For a momentum $p$ of the beam particles and a mass $m_\chi$ of the produced ones, we find that the decay is allowed if the energy density of the beam exceeds the instability threshold $p^2 m_\chi^4/(2g^2)$. We also provide an analytical expression for the spontaneous decay rate at the earliest time.

Taking a two interacting scalar toy model with interaction term gφχ 2 , we study the production of χ particles coming from the decay of an asymptotic and highly occupied beam of φ particles. We perform a nonperturbative analysis coming from parametric resonant instabilities and investigate the possibility that massive χ particles are produced from decays of massless φ particles from the beam. Although this process is not present in a perturbative analysis, our nonperturbative approach allows it to happen under certain conditions. For a momentum p of the beam particles and a mass mχ of the produced ones, we find that the decay is allowed if the energy density of the beam exceeds the instability threshold p 2 m 4 χ /(2g 2 ). We also provide an analytical expression for the spontaneous decay rate at the earliest time.

I. INTRODUCTION
The decay of a particle into other species is one of the simplest and most relevant effects in relativistic field theories. From the theoretical point of view, the decay rate for a process φ i → φ j + φ k can be defined in general terms as [1] where p j,k are the three-momentum of φ j,k , V the volume where the theory is defined, t the time, and |S f i | 2 the transition probability of the process. When the initial and final states are asymptotic, i.e., free of interactions, the matrix element S f i ≡ f | S |i is calculated perturbatively. It gives where M f i is the Feynman amplitude and ω i,j,k the energy of the particles. In this formula, the delta function imposes that the momentum and energy conservation between initial and final asymptotic states must be exactly fulfilled in the transition. In a physical situation, the decay processes originate from beams or clumps of the decaying particles. To compute the produced particle flux the analysis is usually done for the decay of a single particle of the beam or clump and the decay rate is calculated in perturbation theory using Eq. (2), where asymptotic states of the produced particles is assumed. Finally it is integrated over the contributions of all the particles of the beam. This reasoning is valid as long as the assumption of asymptotic final states is correct or at least when it is a good approximation. Following this recipe, it is easy to find that massless particles of a beam can not decay into other * ariel.arza@gmail.com massive species because the energy-momentum conservation for asymptotic states is never fulfilled 1 . More general, this recipe based on the asymptotic state assumption forbids any decay process for which the total mass of the produced particles exceeds the mass of the decaying particles.
However, if the beam or clump carries an extremely big energy density, this assumption is no longer valid, the produced particles are not in asymptotic states because they continue to interact with the decaying field and then nonperturbative effects start to be relevant. As a consequence, Eq. (2) is meaningless and the energymomentum conservation as well as the decay rate should be accounted in a different way. This motivates us to explore the possibility that processes (not only decays) forbidden by the asymptotic final state assumption are now allowed.
In this article we study these nonperturbative effects for a simple scalar toy model. We focus on a beam of massless particles and their decay into massive species. We find that this is indeed possible when the beam energy density exceeds some threshold which is explicitly calculated. As far as we know this possibility has not been discussed in the literature and, despite our study being limited to a theoretical point of view, we expect that it could be relevant in astrophysics, collider physics, and cosmology since in these topics, it is usual to ignore processes that, under the asymptotic state assumption, seem to be impossible. As very intense fields play an important role in these subjects, we suggest that these effects should be something to look at.

II. NON-PERTURBATIVE AXION DECAY AS AN EXAMPLE
To start, let us introduce, as an example, the interesting case of the axion decay into two photons. This inspired us to write this article and certainly serves as a good example to understand what we mean with nonperturbative effects in the previous paragraphs. Assuming asymptotic photon states for the decay of a single axion, the scattering matrix of the process can be computed perturbatively from Eq. (2), then the decay rate obtained from Eq. (1) is given by in the axion rest frame. Here g aγγ is the axion to two photons coupling and m a the axion mass. For an axion clump composed by N a particles (all of them at rest), and assuming, again, asymptotically free photon states after decays, the number of photons per unit time emitted by the clump is 2N a Γ a→2γ . This straightforward result ignores the fact that once photons are produced, they continue to interact with the clump. Indeed, if the axion clump is highly occupied in a particular mode, the presence of the produced photons stimulate the decay of the remaining axions, leading to an exponential growth of the photon occupancy number [6][7][8][9][10][11][12]. This process originates from Bose statistics and is known as Bose enhancement. The effect can also be interpreted as a parametric resonance on the electromagnetic field. Parametric resonance is a well known effect with many applications in cosmology, for instance in the physics of inflation [13][14][15][16][17][18] (see Ref. [19] for a review). Let us discuss briefly this parametric resonance in more detail for the axion-photon case. As it is a nonperturbative calculation, it is better understood for time scales where the clump has not been significantly depleted. Under this assumption one can neglect backreactions on the axion field and the resulting differential equations for the electromagnetic field can be linearized. For a dense axion clump with energy density ρ a pairs of photons emitted with momenta k and − k have a time dependent occupancy number (in the linear regime) given by Here σ a = g aγγ /2 ρ a /2, k = 2k − m a and s k = The most important feature of this result is that each mode k that satisfies 2 k < 4σ 2 a , i.e., experiences exponential growth in its occupancy number. We can compute the photon number density using a saddle point approximation, for σ a t 1 we find From Eqs. (6) and (5) we see that the photon number density grows at the rate affecting a bandwidth of size centered at k = m a /2. This non perturbative axion decay features two interesting properties that we can not find in the perturbative case. First, notice that from Eq. (7) we see that the photons are emitted at a rate different from what we observe in Eq. (3). The most important discrepancy is that while the perturbative result scales as g 2 aγγ , the nonperturbative one scales as g aγγ . References [20,21] give a very intuitive explanation about what is happening. Looking at the Boltzmann equation, one findṡ n a = −Γ a→2γ 1 + 2f γ,ma/2 n a .
We can see that the decay rate is corrected by a factor 1 + 2f γ,ma/2 . When f γ,ma/2 > 1/2, it is enhanced dramatically. This stimulated decay effect is known as Bose enhancement. To probe that f γ,ma/2 > 1/2 can be reached easily, let us analyze Eq. (9) at the beginning, when f γ,ma/2 ∼ 0. Sinceṅ γ = −2ṅ a , after a small period of time δt we get n γ = 2Γ a→2γ n a δt. From the uncertainty principle the photon field is spread in the energy bandwidth δk ∼ 1/δt, then Now, using Eq. (3) and taking Eq. (8) for δk, we find It means that the effect coming from Bose enhancement becomes important instantly. Moreover, using Eq. (10), Eq. (9) can be written aṡ which is consistent with Eq. (6), at least in the order of magnitude. The other property, maybe the most important concerning our hypothesis, comes from Eq. (5). While in the perturbative calculation the energy momentum conservation for asymptotic states implies that all the produced photons must have a momentum k = m a /2, the nonperturbative analysis allows the produced photon to have momenta within the bandwidth defined in Eq. (5). As this window for the photon momenta is proportional to √ ρ a , we would recover the perturbative result in the limit ρ a → 0. This limit is indeed consistent with the asymptotic state assumption since for a small axion energy density the interaction of the produced photons with the axion field is negligible. This window for the momenta of the produced particles, photons for the axion case, is what leads us to wonder about the possibility of decays that are forbidden in the perturbative calculation, where the asymptotic state assumption is an unbroken rule. For a more general decay process of the form φ i → φ j + φ k , if the decaying field φ i is highly occupied, the linearized equations of motion for the fields φ j and φ k look quite similar to the ones of the axion-photon case. While the analog of σ a has in general a momentum dependence (σ a → σ k ) and a different form (it also depends on the form of the interaction term), it always scales as the squared root of the energy density ρ φi of the decaying field φ i . On the other hand, and also very important, k has always the same form; i.e., if the decaying field φ i is occupied in a state with momentum p, where ω and m l the corresponding mass of each field φ l . The parametric resonance and depletion of the field φ i takes place when Notice that in the limit ρ φi → 0, suitable for the single particle decay, the condition k = 0 is required. This condition is exactly the energy-momentum conservation for asymptotic states. We are interested in the case where φ i is massless 2 and φ j , φ k are massive (or at least one of them). In this case k can never be zero and, therefore, a single massless particle can not decay by itself. However, when φ i is a highly populated beam in some momentum mode, the condition for the decay is actually Eq. (14). Particles from the massless beam start to decay into φ j and φ k with momentum modes satisfying Eq. (14) even if k = 0 has no solutions for k.
When non-perturbative effects take place, of course the energy conservation can not be accounted by assuming asymptotic states. In our approach we do it by calculating the energy expectation value of the produced field 2 In models of massless preheating [22] the inflaton field suffers selfinteractions due to a λφ 4 potential. It makes the zero momentum modes of the inflaton to oscillate harmonically. It finally drives the decay of the field into SM particles by parametric resonance. Our interest is rather directed at a beam of massless particles and their possibility to decay into massive species.
directly from the solutions of the equations of motion. The details will be discussed later. At this point we have taken the parametric resonance properties to explain roughly how a beam of massless particles can be depleted by the decay into massive particles. In the following we will discuss this mechanism in detail using a simple toy model involving two scalar fields. Our analysis includes solutions of the equations of motion in the rotating wave approximation, study of the parametric resonance instabilities, computation of the energy density threshold for the decays and consistency check of energy conservation. Some technical details and validity check of the approximation are left in the Appendix.

III. A TWO SCALAR TOY MODEL
We consider the following interaction hamiltonian density where φ represents the decaying particles, χ the produced particles and g the coupling constant. For now we consider the general case where both fields are massive. Working in the Heisenberg picture, the equations of motion of the system are where m φ and m χ are the masses of φ and χ, respectively. We are interested in the case where φ is a beam of particles, highly occupied in a single momentum mode p. In this limit, φ can be perfectly considered as a classical field. On the other hand, as in the axion-photon case, we consider timescales where φ is not significantly depleted. It allows us to neglect backreactions over φ, so we make the rhs of Eq. (16) equal to 0 and Eq. (17) becomes linear. These considerations lead us to write φ as a monochromatic classical plane wave as where ρ φ is the time averaged energy density of the beam and ω p = p 2 + m 2 φ . We write the quantum field χ in terms of creation and annihilation operators as where Ω k = k 2 + m 2 χ and the operators χ k and χ † k satisfy the commutation relations [χ k , χ k ] = 0 and [χ k , χ † k ] = (2π) 3 δ 3 ( k − k ). Inserting Eqs. (18) and (19) into Eq. (17), we get the following set of equations for the new operators where we have defined Let us study the resonant solutions of Eq. (20) analytically using the rotating wave approximation (RWA) 3 . To do so we write χ k (t) = a k (t) e −iΩ k t . The RWA is a method that allows to transform a second order system of differential equations to a first order one, in our context it can be done when a k varies slowly respect to χ k . For more details of the method, see Appendix A.
We are interested in the case where at t = 0 there is no χ particles. Therefore, at the very beginning the only relevant process 4 to account is φ → 2χ. Defining and neglecting non-relevant terms, Eq. (20) becomeṡ where The general solution of Eq. (23) is where we have defined For each momentum mode k, its occupancy number is Figure 5 (see the Appendix) shows a perfect agreement between numerical solutions of Eq. (20) and solutions in the rotating wave approximation. 4 The decay φ → 2χ is the main process only at the beginning.
When χ is populated, processes of the form φ + χ → χ are activated. See Appendix A and B for a more detailed discussion. and the total density of produced particles is given by The modes that exhibit parametric resonance are the ones that satisfy Now let us talk about energy conservation (momentum conservation is satisfied since we obtained Eq. (20)). After N φ decays of φ particles, the energy expectation value of χ must be E χ = ω p N φ in order to conserve energy during the processes. We compute E χ directly by E χ = 0| H χ |0 , where the free hamiltonian H χ is given by By using Eqs. (19) and (25) we get (see Appendix C) where we have ignored the vacuum energy contribution. Defining q = p − k we find where we have used Eq. (22), the property f χ, k = f χ, q and the fact that the average number of produced χ particles N χ = V n χ is twice N φ . Since d 3 k = d 3 q, the integral in Eq. (32) cancels, getting as expected.
Notice that this result as well as the ones in the previous paragraph do not depend on the masses of the particles. The only requirement is a parametric instability. We will show later that even the case of massless φ and massive χ is unstable for some region in momentum parameter space if ρ φ is large enough. In other words, we will show that, after some threshold of ρ φ , the production of massive scalar from decays of massless particles is possible without energy-momentum violation.
Before continuing, we would like to discuss the validity of our calculation. Our results were found in a RWA where the amplitudes a k vary slowly respect to χ k , so Eq. (25) is valid under the condition Ω k s k . Assuming Ω k ∼ Ω p− k ∼ ω p we would be safe for Even though condition Eq. (34) will be useful in most of the remaining discussion, we have to pay special attention when either Ω k or Ω p− k approaches to zero because s k could blow up. For instance, when m χ = 0, Ω p− k = 0 if k = p.

IV. INSTABILITY CONDITION AND DECAY OF A MASSLESS BEAM
Now we are going to study the instability conditions in two opposite cases: (1) massive φ and massless χ, and (2) massless φ and massive χ. Case (1) has been strongly discussed in the literature, so we are using it as a matter of comparison with case (2), which is what concerns us. Let us define the dimensionless quantities κ = k/ω p , µ = m χ /ω p , v = p/ω p , and β κ = κ 2 + µ 2 . Let us also define η κ = s κ /ω p , which in terms of the defined dimensionless parameters is given by We also define θ as the angle formed by k and p.
For case (1) we can always choose a reference frame where the φ particles are at rest, then the beam becomes a clump of particles with speed v = 0. Notice that this is very similar to the axion-photon case discussed before 5 . Eq. (35) simplifies and the instabilities (η 2 κ > 0) take place for 1/2 − α < κ < 1/2 + α. For case (2) Eq. (35) must be evaluated with v = 1. In this case the analysis is more complicated, so before finding analytical expressions we will first identify some general properties by showing some plots. Figure 1 shows η κ as a function of µ for α = 10 −3 , θ = 0, and different values of κ. Of course we avoid the extremes κ = 1 and κ = 0 where our calculation breaks down. We can see instabilities only for µ 1. Figure 2 shows η κ as a function of θ for α = 10 −3 , µ = √ α/2, and different values of κ. We observe that the instabilities occur for θ smaller than 1. To have a better look of the instability regions, see Fig. 3 where η κ is plotted as a function of κ for a fixed 5 It is also analogous to the inflaton decay where the inflaton field is considered as an homogeneous classical field [20]. µ and different values of θ. We see that in terms of κ the bandwidth is of order 1 for θ = 0, much bigger than case (1) where it is of the order of α. However, we also see that the instability region becomes smaller as θ increases, which is the main advantage of case (1) where the instability window keeps its size for any angle.
We also see in Fig. 3 that, close to the extremes, κ ≈ 1 and κ ≈ 0, η κ is enhanced. It may contribute to the efficiency of the process, although we have to be always alert that s k /Ω k 1 is satisfied, otherwise we would enter in regions where our RWA breaks down.
These plots gave us some generals features of the instability conditions, in terms of physical parameters we found that m χ ω p and θ 1. Now we will use these features to find analytical expressions. The frequencies can be approximated as Ω k ≈ k + m 2 χ /(2k) and Ω p− k ≈ p − k + (pkθ 2 + m 2 χ )/(2(p − k)). When computing k , the main contributions of Ω k and Ω p− k cancel with ω p , getting In the above result we have assumed that we are far enough from k = 0 and k = p.
To find the instability condition we solve the equation η 2 κ = 0 to get the instability limits for κ. As θ = 0 gives the biggest instability window for κ, we evaluate at this angle obtaining the limits The size of the window is therefore the instability is triggered when µ < √ α or when in terms of the energy density.
It is not too difficult to check that the instability condition Eq. (39) is Lorentz invariant. It is, of course, what we would expect. Now two questions arise; why the instabilities were found for p m χ ?. What happens for p ∼ m χ or p m χ ?. We answer the first question as follows. First, the instability occurs when µ < √ α. Second, we are working in a RWA where the validity of our analysis is limited by α 1. These two facts imply that the instability must satisfy p m χ as a consequence. In other words; in our rotating wave approximation, the instability takes place only for reference frames where p m χ . To answer the second question, we take Eq. (20) in the limit p m χ . Taking into account only the modes k and k − p we get the system For solutions of the form A k = c k e γt e −iω p t/2 and A k− p = d k e γt e iω p t/2 , we find that γ is given by It is clear that the instability is also activated when Eq. (39) is fulfilled. We also expect the same result for p ∼ m χ , but it is not clear how to proceed analytically.

V. SPONTANEOUS DECAY RATE
So far we have explained that the decays of massless into massive particles are allowed due to a nonperturbative effect caused when a huge energy density of the beam prevents asymptotic states for the produced particles. However, we still do not clarify how the first χ pair is spontaneously released. To do so, we focus on the decay of a single φ-particle of the beam but, of course, taking into account the effect that the beam energy density causes on it. Working now in the interaction picture, the only non-zero element of the S-matrix expansion is where |i is the initial state that contains one φ particle with momentum p and |f the final state composed by two χ particles, one with momentum k and another with momentum q. In Eq. (43), T is the time during which the system transits from the state |i to the state |f . As the transitions that are typically studied involve processes where the final state is free of interactions, one could just make T → ∞. However, in some cases this approximation breaks down due to finite time effects [23][24][25] or simply if the transition time is not big enough to be considered as infinite [26]. In our case, this is indeed what happens. From the nonperturbative analysis we found that if the instability condition is fulfilled, there is no asymptotic states and that the time scales for transitions into unstable modes are of the order (2s k ) −1 (see Eq. (27)).
Without saying anything about T we have (44) Replacing Eq. (44) into Eq. (1) and applying an extra factor 1/2 coming from the fact the final states are identical particles, we get Let us first consider the scenario where ρ φ → 0. In this case the system transits to a free of interaction state after the decay, so we can take T → ∞. In this limit sin( k T )/ k → πδ( k ), then the usual energy-momentum conservation relation, k = Ω k +Ω p− k −ω p = 0, involving asymptotic initial and final states, holds. For case (1) we, of course, obtain the usual formula For case (2) we get Γ φ→2χ = 0 because k can never be zero. It is the usual result that prohibits the decay of massless into massive particles due to the lack of final asymptotic states that satisfy energy-momentum conservation. Now we consider the case where ρ φ is big enough to trigger the instabilities discussed before, i.e. it satisfies Eq. (39). It is clear that there is no asymptotic states for the χ particles produced from the decay, the evolution of the χ modes affected by parametric resonance are rather described by Eq. (25). In parametric resonance the number of the produced particles grows exponentially as Eq. (27), at least at the beginning, therefore the transition time of the process is T = (2s k ) −1 . Since the main contributions for the growth are provided by modes that fulfill s k k , we can use the approximation sin( k T ) = sin( k /2s k ) ≈ k /2s k with great accuracy. Equation (45) becomes For case (1) we can make σ k ≈ m φ α to preserve only terms of order g 2 . We find As we get the same as Eq. (46), the physical reason of this result should be investigated further. For case (2) the integration is not straightforward as in case (1). Assuming θ 1, which is true on behalf of previous discussions, we get where (51) We recall that κ ± (µ/ √ α) are defined in Eq. (37). In terms of the quantity ω p , which is m φ for case (1) and p for case (2), the ratio between both decay rates is simply A plot of F (µ/ √ α) is shown in Fig. 4. As we already discussed, we can see that for case (2) the decay is clearly possible if µ < √ α, but it becomes impossible for µ ≥ √ α. For µ √ α the decay rate approaches its maximum which is a half of the decay rate for case (1).
In our toy model we can also connect the g 2 dependence of the earliest spontaneous decay rate and the rate proportional to g at which χ particles are produced by parametric resonance. As for the axion-photon discussion in Sec. II, we take the Boltzmann equatioṅ n φ = −Γ φ→2χ (1 + 2f χ )n φ , where f χ ∼ (2π) 3 n χ /d 3 k and Γ φ→2χ ∼ g 2 /(8πω p ) for both case (1) and case (2). It is not difficult to find that also for both cases d 3 k ∼ πg 2ρ φ . Combining all of this and the fact thaṫ n χ = −2ṅ φ , we find from the Boltzmann equation thaṫ (1) as well as for case (2).

VI. DEPLETION OF THE MASSLESS BEAM
To finish we will briefly discuss the conditions that allow a significant depletion of the φ field. We already said that condition Eq. (39) triggers the instability, however it is a necessary but not sufficient condition for a substantial decay of the beam. The missing condition is that the time t res during which the instability takes place is enough to reach the exponential regime. It requires α ω p t res > 1. We need the produced particles to stay inside the beam extension for the parametric resonance to develop, therefore t res is the time that χ particles take to leave the beam region. As the χ particles have momenta pointing mainly to the same direction as the beam propagates, t res is given roughly by d/(1 − v χ ), where v χ is the velocity of the produced particles and d the beam length. This velocity is approximately 1 − µ 2 /2, which leads to the following depletion condition Notice that, given the resonance condition µ 2 < α, Eq. (53) is satisfied necessarily by taking d > 1/(2p).

VII. CONCLUSION
In this article we have taken a simple scalar toy model, with interaction term gφχ 2 , to study the process φ → 2χ in the particular case of massless φ and massive χ. Although the one particle analysis prohibits the decay by the requirement of energy-momentum conservation of asymptotic states, we found that when φ is highly occupied, nonperturbative effects (parametric resonance) allow the decay as long as the energy density of the decaying particle beam exceeds the threshold defined in Eq.
(39). We demonstrated that energy-momentum is indeed conserved during the particle production and, finally, we got an analytical expression for the spontaneous decay rate at the earliest time. This result suggests a deeper inspection of nonperturbative effects in particle processes, especially the ones that are relevant in cosmology, astrophysics and collider physics.
For case (2) (1) k is big, of the order of Ω k , but . It means that if the process φ p → χ k +χ p− k is efficient, once χ k and χ p− k are abundant, the processes χ k + φ p → χ k+ p and χ p− k + φ p → χ 2 p− k will become efficient too. Following the same reasoning, a cascade is produced and every process of the form χ k+n p + φ p → χ k+(n+1) p and χ n p− k + φ p → χ (n+1) p− k (n = 0, 1, 2, ...) will be eventually important after some period of time. To prove this last claim, we compute the classical version of Eq. (20) numerically assuming that only the modes k and p − k are initially occupied. Figure 5 shows |A k | 2 as a function of time for the modes k, k + p, k + 2 p and k + 3 p. We do not plot the modes p − k, 2 p − k, 3 p − k and 4 p − k because, according to the initial conditions, the time evolution of their square amplitudes are identical as the ones for k, k + p, k + 2 p and k + 3 p, respectively. We evaluate at k = ω p /2 and we assume that at the beginning the modes k and p − k have initial conditions |A k (0)| 2 = |A p− k (0)| 2 = 1. All the other modes have a null initial amplitude. We see that the squared amplitudes of modes different from k and p − k begin to grow gradually. In Fig. 5 we also compute numerically in the rotating wave approximation (see dotted lines). We basically take Eq. (A1) and solve the system neglecting fast oscillating terms. It is clear that this approximation perfectly agrees with the full solutions of the second order differential equations.