Particle production from oscillating scalar field and consistency of Boltzmann equation

Boltzmann equation plays important roles in particle cosmology in studying the evolution of distribution functions (also called as occupation numbers) of various particles. For the case of the decay of a scalar condensation ϕ into a pair of scalar particles (called χ), we point out that the system may not be well described by the Boltzmann equation when the occupation number of χ becomes large even in the so-called narrow resonance regime. We study the particle production including the possible enhancement due to a large occupation number of the final state particle, known as the stimulated emission or the parametric resonance. Based on the quantum field theory (QFT), we derive a set of equations which directly govern the evolution of the distribution function of χ. Comparing the results of the QFT calculation and those from the Boltzmann equation, we find non-agreements in some cases. In particular, in the expanding Universe, the occupation number of χ based on the QFT may differ by many orders of magnitude from that from the Boltzmann equation. We also discuss a possible relation between the evolution equations based on the QFT and the Boltzmann equation.

(Hereafter, such a scalar condensation is called φ.) If φ can decay into a pair of bosonic particles (called χ) as φ → χχ, the effects of the stimulated emission can be important since the daughter particles may be enormously populated. In particular, recently, a new mechanism of producing bosonic dark matter from the inflaton decay has been proposed [17], in which a stimulated emission of the bosonic dark matter plays an important role.
Such systems have been studied by employing the Boltzmann equation or in the context of the parametric resonance [2,[18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36]. In particular, particle production from an oscillating scalar field has been intensively investigated, particularly using the Mathieu equation [37,38]. In the previous studies, the evolution equation of the expectation value of the field operator (which we call wave function) for the final-state particle is converted to the Mathieu equation, based on which the occupation number of the final state particle has been estimated. Then, it has been shown that there exist resonance bands and that the occupation numbers in the resonance bands grow exponentially. In the broad resonance region, it has been known that the particle production is non-perturbative and that the JHEP03(2021)296 process cannot be described by the Boltzmann equation. On the contrary, sometimes it has been argued that the Boltzmann equation can provide a proper description in the narrow resonance regime.
In this paper, we study the particle production from an oscillating scalar field in the narrow resonance regime. We pay particular attention to the relation between the results based on the Boltzmann equation and those from the quantum field theory (QFT). Even though the effects of the stimulated emission may be included in the Boltzmann equation [1,2], it is unclear whether the Boltzmann equation can properly describe the particle production from the oscillating scalar field from the QFT point of view, in particular, the parametric resonance. Based on the QFT, we derive a set of equations which directly govern the evolution of the occupation number of the final state particle χ with an oscillating φ background. We will see that the results of our evolution equations are in good agreement with those from the conventional approach with the Mathieu equation. Then, we compare the occupation number obtained from the QFT calculation with that from the Boltzmann equation. We will see that the results from the QFT and the Boltzmann equation may differ even in the narrow resonance regime when the occupation number becomes larger than ∼ 1. In particular, in the expanding Universe, growth factors derived from two approaches may differ many orders of magnitude. We also discuss how our evolution equations can be related to the conventional Boltzmann equation and consider a possible explanation of the discrepancy. This paper is organized as follows. In section 2, the setup of our analysis is given. In section 3, we derive the evolution equations and discuss particle production in the flat spacetime. In section 4, we study particle production taking into account the cosmic expansion. Section 5 is devoted to the summary of this paper.

Setup
Here, we study particle production from an oscillating real scalar field φ, whose mass is m φ . We concentrate on the case that φ is spatially homogeneous and the scalar field φ depends only on time t. We consider the timescale much shorter than the dissipation of the motion of φ and the back reaction. Then, assuming that the potential of φ is well approximated by a parabolic one, the motion of φ is described as whereφ is the amplitude. The number density of φ is given by The scalar field φ is assumed to interact with a real scalar field χ. Although our discussion holds as far as the interaction Hamiltonian has the form of eq. (3.5) given below, we consider the following interaction Hamiltonian to make our discussion concrete: 3)

JHEP03(2021)296
where A is a dimension 1 coupling constant. With the above interaction Hamiltonian, the decay rate for the process of φ → χχ is obtained as where p χ is the three-momentum of χ produced by the decay: with m χ being the mass of χ. In eq. (2.4), the superscript "(0)" indicates that Γ (0) φ→χχ is the perturbative decay rate in the vacuum.
The production of χ in the system introduced above has been also studied in the context of parametric (or tachyonic) resonance [19,23,24,29]. In particle production through the parametric resonance, the modes in resonance bands are effectively produced. In the present case, the widths of the bands are determined by the following parameter: As we will see below, the band width is ∼ qm φ . Hereafter, for the comparison with the Boltzmann equation, we concentrate on the case that the band widths are narrow, i.e., In the narrow resonance regime, χ is produced almost "on-resonance," i.e., the momentum of χ produced from the φ oscillation is | k| p χ . In the perturbative description, effects which are higher order in A are suppressed when q is small [25].

Evolution equations from QFT
In studying particle production in the QFT, an approach by employing the Mathieu equation has been often used in the context of parametric resonance. In our analysis, however, we adopt a different approach, which should be equivalent to the one with the Mathieu equation, by deriving evolution equations for the distribution function. An advantage of our approach is that the equations directly describe the evolution of the distribution function. Because of this, we can consider a possible relation between the evolution equation of the distribution function in the QFT and the Boltzmann equation.
We first discuss particle production in the flat (i.e., Minkowski) spacetime. We decompose χ by using the creation and annihilation operators, denoted as a k and a † k , respectively:

JHEP03(2021)296
The creation and annihilation operators satisfy the following commutation relation: With the creation and annihilation operators, the free part of the Hamiltonian is given by where In addition, the interaction Hamiltonian is expressed as In order to study the production of χ, we introduce the distribution function of χ: where V is the spatial volume, and the expectation value of the operator O is defined by using a density matrix ρ as Using the distribution function, the number density of χ is given by In the following, we work in the interaction picture. The time dependence of the operator O is given byȮ where the "dot" denotes the derivative with respect to time, and hencė In addition, the density matrix evolves aṡ Consequently, the evolution of the expectation value of the operator O is governed by the following differential equation:

JHEP03(2021)296
In order to study the evolution of the distribution function, we also introduce Then, the time derivatives of f k and g k are given bẏ 14) We consider isotropic solutions, and the momentum dependences of f k and g k are only via k ≡ | k|; thus, in the following, these functions are denoted as f k and g k , respectively. We decompose the complex function g k as where ξ k and η k are real functions. Then, the evolution equations becomė We note that the evolution equations above also apply to the case of q 1, although we do not consider such a case in this paper. Before solving these differential equations, we comment on the initial condition. Because we study the χ production from the φ oscillation, we consider the case that the χ sector is initially at the ground state (vacuum), denoted as |0 ; thus, We numerically solve eqs. (3.17)-(3.19) with the initial condition given by eq. (3.20). The evolutions for ω k = 1 2 m φ and m φ are shown in figures 1 and 2, respectively, with taking q = 10 −2 . As one can see, f k monotonically increases when t m φ irrespective of ω k (as far as ω k m φ ). The behavior at t → ∞ strongly depends on ω k . For ω k = 1 2 m φ , the functions f k and η k grow exponentially at t → ∞, while ξ k is oscillating. On the contrary, for ω = m φ , f k does not show the growing behavior at t → ∞.
The behaviors of f k , ξ k , and η k at t m φ can be obtained by expanding these functions as a power series of t. With t being small enough, the effects of the oscillation of JHEP03(2021)296

JHEP03(2021)296
φ is unimportant and we obtain We can also understand the exponential growths of the modes with ω k ∼ 1 2 m φ at t → ∞. We first adopt the following ansatz: where α k and λ k are positive constants. Substituting the above ansatz, eq. (3.15) becomes where we neglect terms which are not exponentially growing, and For the study of the modes with ω k ∼ 1 2 m φ , we neglect the second term of the right-hand side of the above equation because it is sub-dominant. Combining eqs. (3.14) and (3.26), we obtaiṅ (3.27) We neglect terms which are oscillating with the timescale of ∼ m −1 φ because we are not interested in terms with rapid oscillations. We define the time average of the function F (t) as where the timescale t * is taken to be m −1 Using the relation cos 2 m φ t t 1 2 and sin 2m φ t t 0, we can see that eqs. (3.14) and (3.15) are satisfied with the ansatz given in eq. (3.24) (neglecting rapidly oscillating terms and terms which do not grow exponentially) if  Table 1.
, based on the numerical calculation (Numerical) and the approximation adopting the ansatz of the exponential increase (eq. (3.24)) with the growth rate given in eq. (3.29). We take q = 10 −2 , and several values of ω k (parameterized The growth rate λ k should be real, which determines the resonance band. We can see that the width of the resonance band is of O(qm φ ). Then, keeping the leading term in q, we may use 30) and the resonance band is given by Notice that the resonance band given above corresponds to the first resonance band in the study of the parametric resonance, and that the growth rate given in eq. (3.30) is consistent with that given in [21,23,24]. In order to check the validity of the ansatz of the exponential increase given in eq. (3.24), we compare the value of f k (t = 100(qm φ ) −1 )/f k (t = 10(qm φ ) −1 ) obtained from the numerical calculation and that predicted by eq. (3.24) (with the growth rate given in eq. (3.29)). The results with taking q = 10 −2 are shown in table 1 for several values of ω k in the resonance band. Notice that f k (t = 10(qm φ ) −1 ) is of O(10) for the frequencies considered in table 1, and hence the effect of stimulated emission is already important when t = 10(qm φ ) −1 . We can see that the occupation number in the resonance band is well described by the exponential increase with the growth rate given above once the effect of the stimulated emission becomes effective.

QFT and ordinary Boltzmann equation
The relation between the analysis so far and that based on the conventional Boltzmann equation is interesting. Sometimes particle production from the oscillating scalar field φ is studied by using Boltzmann equation. φ can be regarded as a coherent state of a nonrelativistic scalar field and the production of χ may be regarded as the decay of φ. Then, the collision term in the Boltzmann equation is given by (see, for example, [1,2])

JHEP03(2021)296
where the first and the second terms in the square bracket describe the effects of the decay and the inverse decay of φ, respectively. In the following, we consider how the collision term can be related to the argument based on the QFT. We treat the cases of t (qm φ ) −1 and t (qm φ ) −1 separately. When t (qm φ ) −1 , f k is smaller than 1 and the solution of eq. (3.15) is approximately given by where c is a constant. Here, we have neglected the last term of the right-hand side of eq. (3.15) because it is of O(q) and is sub-dominant relative to the first term. 2 Because g k (0) = 0 and also because g k is non-singular at ω k → 1 2 m φ , c = 1 and henceḟ k becomeṡ In the above expression, we neglect the terms which are rapidly oscillating with the timescale of ∼ m −1 φ (even in the limit of ω k → 1 2 m φ ). Notice that the function γ(ω; t) has the following property: ∞ 0 dωγ(ω; t) = 1 + cos 2m φ t. (3.36) We consider the time averaged value ofḟ k and neglect terms which are rapidly oscillating: Using eq. (3.36), we obtain The homogeneous solution of eq. (3.15), i.e., the solution with taking (1 + 2f k ) → 0, is given by with c being a constant. In eq. (3.33), the second term in the bracket is neglected because it gives a higher order contribution in terms of q.

JHEP03(2021)296
In addition, Thus, assuming that t * can be chosen to be large enough, we may approximate and obtain When t (qm φ ) −1 , f k in the resonance band shows the exponential growth and is much larger than 1, while g k behaves as eq. (3.25). In such a case, the conventional Boltzmann equation is obtained if one performs the following replacement in eq. (3.25): Here, the relation 1 x−i0 = P ( 1 x ) + iπδ(x) (with P being the principal value), as well as the smallness of λ k relative to m φ , are used. However, we should note that this replacement is allowed if λ k does not depend much on ω k and also if f k (t) is insensitive to k around the pole. These requirements may not be satisfied in the case of our interest (see eq. (3.30)), which may result in the non-agreements between the result of the QFT calculation and that of the Boltzmann equation. Substituting eq. (3.25) (with using the above replacement) into eq. (3.14), we obtaiṅ Neglecting the rapidly oscillating term with taking the time average, Combining the results for t (qm φ ) −1 and t (qm φ ) −1 , which are given in eqs. (3.42) and (3.45), respectively, and estimating the collision term aṡ we find that the collision term given in (3.32) well describes the evolution of f k if the assumptions and approximations adopted in the above argument are valid. In order to see the validity of the collision term from the QFT point of view, we calculate the number density of χ with two different methods, taking m χ = 0 (and hence ω k = k) and q = 10 −2 .
• First, we use the QFT approach; the number density based on the QFT is denoted as n Notice that, in the above expression, the distribution function in the parenthesis is that for k = p χ and is obtained by solving eqs.  with larger t requires a computational cost.) For t (qm φ ) −1 , on the contrary, the peak width of f k is smaller than the width of the resonance band (see below). For the replacement given in eq. (3.43), it is required that the function f k is well approximated by f pχ for |k − 1 2 m φ | λ k and that the ω k dependence of λ k is unimportant. These cannot be the case in particular when t (qm φ ) −1 . Consequently, n (Boltzmann) χ becomes larger than n (QFT) χ , as shown in figure 3. We can understand the asymptotic behaviors of n (QFT) χ and n (Boltzmann) χ for t (qm φ ) −1 as follows. At t (qm φ ) −1 , as we see below, f k is sharply peaked at ω k = 1 2 m φ . We expand λ k around the peak and obtain (3.49) Because of the second exponential factor, the width of f k (t) gets smaller as t increases. By using the above expression with neglecting m χ , n We can see that the above expression is in a good agreement with the numerical result. In addition, Thus, when t (qm φ ) −1 , n (Boltzmann) χ (t) becomes larger than n (QFT) χ (t). One may wonder if we can introduce an "averaged" occupation number (or growth rate) in the resonance band to make two approaches consistent. However, as we will see in the next section, the discrepancy in the case with the cosmic expansion cannot be solved with such a prescription. (In addition, in appendix A, we give a consideration about the relation between the QFT and Boltzmann equation.) Before closing this section, we comment on the back reaction. In our analysis, the effects of the back reaction are neglected; we assume that the amplitude of the φ oscillation does not depend on time and that the motion of φ is well described by eq. (2.1). This is the case when the energy density transferred to the χ sector is smaller than the initial energy density in the φ sector. Let us denote the typical value of the occupation number of the modes in the resonance band as f (res) ; conservatively, we may take f (res) ∼ f k=pχ . Then, by using the fact that the width of the resonance band is ∼ qm φ , the energy density in the χ sector is estimated as ρ χ ∼ f (res) qm 4 φ . Requiring that ρ χ is smaller than the initial energy density of φ, we obtain and, equivalently, f (res) A φ . The above constraint can be satisfied for any value of q by taking large enoughφ (and small enough A). Our results are applicable to the parameter JHEP03(2021)296 region consistent with the above constraint. We may also define the effective decay rate Γ (eff) φ→χχ as the inverse of the timescale with which a single φ becomes a pair of χ; we can estimate Γ (eff) φ→χχ . With the constraint (3.52) being satisfied, Γ (eff) φ→χχ is always smaller than qm φ as far as q 1. Another back reaction may be due to the scattering process like χφ → χφ (with φ in the initial state being that in the coherent oscillation). Such a scattering process may remove χ from the resonance band. Effects of the scattering processes are not taken into account in our analysis because φ is treated as a classical field.
One can estimate the interaction rate as Γ (scat) ∼ q 4 m 3 φ φ 2 which can be neglected compared with the growth rate which is of O(qm φ ).

Particle production with cosmic expansion
In this section, we study particle production taking into account the effects of the cosmic expansion. We show that the analysis based on the Boltzmann equation may result in a significant overestimation of the occupation number of χ compared to the analysis based on the QFT.

Particle production with cosmic expansion from QFT
With the cosmic expansion, the momentum of each mode redshifts. Thus, if we consider a mode which has a frequency larger than ω + in an early epoch, it enters the resonance band as the Universe expands, then the frequency becomes smaller than ω − and the mode exits the resonance band. The occupation number may exponentially increase when the mode is in the resonance band, as we see below.
Here, we are particularly interested in the behavior of the occupation number when the mode is around the resonance band. The timescale of the mode to go through the resonance band is ∼ qH −1 , where H is the expansion rate: with a being the scale factor. The timescale of the evolution of the occupation number is much shorter than the timescale of the cosmic expansion because q 1 and hence the change of H is unimportant. In the following analysis, we neglect the time dependence of H and take with a 0 being a constant. Due to the same reason, again we treatφ as constant.
Because of the hierarchy between the timescales of the cosmic expansion and the evolutions of f k and g k , for the timescale shorter than H −1 , evolutions of f k and g k are expected to be the same as those in the flat spacetime. Then, f k and g k obeẏ

JHEP03(2021)296
where the first terms of the right-hand sides of the above equations describe the effect of redshift. We can simplify solve the above equations by introducing functions for a fixed comoving momentum: where K(t) is the physical momentum for the given comoving momentumk (which is independent of time): We decompose the function g using real functions ξ and η as where with Then, eqs. We numerically solve eqs. (4.10)-(4.12). We take m χ = 0 for simplicity, and i.e., Ω(0) = 1 2 m φ . We choose t i and t f such that Ω(t i ) ≥ ω + and Ω(t f ) ≤ ω − , and follow the evolutions of f and g from t = t i to t = t f ; here, we take and hence The cosmic time at which the mode enters (exits) the resonance band is denoted as t + (t − ):  Table 2. The ratio f (t − )/f (t + ), based on the numerical calculation (Numerical) and eq. (4.20), taking q = 10 −2 . We take several values of q and H (parameterized as H = hq 2 m φ ). For comparison, we also show the value of the right-hand side of eq. (4.23).
i.e., Ω(t ± ) = ω ± . We impose the following initial condition at t i : We show the behavior of f for q = 10 −2 in figures 4 and 5, in which we take H = 10 −2 q 2 m φ and H = 5 × 10 −3 q 2 m φ , respectively. We can observe exponential growths of f when ω − Ω(t) ω + , while there is no significant increase of the occupation number when the mode is outside of the resonance band.
For the case that f (t − ) 1, the total amount of the increase of f can be estimated by using the growth rate λ k given in the previous section. Because the enhancement of f occurs when the mode is in the resonance band, we concentrate on the epoch of t + ≤ t ≤ t − ; during such an epoch, the occupation number is expected to evolve aṡ f k kH(∂f k /∂k) + λ k f k . Then, adopting the growth rate given in eq. (3.30),ḟ is given byḟ which results in For the case that the mass of χ is negligible, the above ratio is estimated as [20,22,23] In order for a significant enhancement of the occupation number, H should be smaller than q 2 m φ .
To show the validity of the above estimation, we calculate the ratio f (t − )/f (t + ) by using eq. (4.20) and also by numerically solving eqs. (4.10)-(4.12); for the numerical calculation, we adopt the setup given by eqs.  JHEP03(2021)296

Breakdown of Boltzmann equation
Finally, let us come to our main point. We compare the above results with that obtained by using the Boltzmann equation. Adopting the collision term given in eq. (3.32), the Boltzmann equation with the cosmic expansion iṡ 21) or equivalently,ḟ = 2Γ One can easily solve this equation and obtain We can see that, when the effect of the stimulated emission is effective, the enhancement factor suggested from the Boltzmann equation is exponentially larger than that from the argument based on the QFT (see eq. (4.20)); the exponent is doubled in the result from the Boltzmann equation. (For comparison, we also show the right-hand side of (4.23) in table 2.) Notice that the discrepancy in the present case is rather serious than that in the flat spacetime. This is because, here, we consider the evolution of the mode with a fixed comoving momentum which goes through the resonance band; thus the discrepancy cannot be solved even if we consider a prescription to average the growth rate in the resonance band. We may naïvely set an ansatz for a Boltzmann equation with including the quantum effect. We can find that the result based on the QFT is well described by eq. (4.21) with replacing (1 + 2f k ) → (1 + f k ) in the right-hand side, i.e., (4.24) The solution of the above equation is given by which has the same growing behavior as that in the QFT (see eq. (4.20)). Interestingly, eq. (4.25) well describes the numerical result for any value of f (t → ∞). The theoretical justification of this equation will be considered in the future [39].

Summary
In this paper, we have studied particle production from an oscillating scalar field φ, assuming that the final state particle χ is very weakly interacting. We have paid particular attention to the consistency of the results from the Boltzmann equation and those from the QFT calculation. We have concentrated on the case that the production of χ is via the process φ → χχ.

JHEP03(2021)296
First, we have considered particle production in the flat spacetime. In such a case, we have discussed the evolution of the occupation number of each mode (i.e., mode with a fixed momentum k) separately in the narrow resonance regime. We have derived the evolution equations for the occupation number of each mode based on the QFT. A resonance band shows up at ω k close to 1 2 m φ , which corresponds to the lowest resonance band in the context of the parametric resonance. The modes within the resonance band can be effectively produced. For the timescale much longer than (qm φ ) −1 , the occupation numbers of the modes in the resonance band exponentially grow; the growth rate obtained in our analysis is consistent with that given by the study of the parametric resonance using the Mathieu equation. Then, comparing the occupation number obtained from the QFT calculation with that from the Boltzmann equation, we have found that they do not agree well when the occupation number is larger than ∼ 1. On the contrary, when f k 1, we have found a good agreement of two results. We have also argued how our evolution equation based on the QFT could be related to the ordinary Boltzmann equation. When the occupation number is larger than ∼ 1, some of the approximation and assumption necessary for such an argument cannot be justified, which, we expect, causes the disagreement.
Then, we have studied particle production taking into account the effects of cosmic expansion. With the cosmic expansion, the physical momentum redshifts. The momentum of each mode stays in the resonance band for a finite amount of time and then exits the resonance band. The exponential growth of the occupation number occurs only in the resonance band. The growth factor has been studied numerically and analytically, adopting the evolution equations based on the QFT. The agreement between numerical and analytical results is excellent. We have also analyzed the system by using the conventional Boltzmann equation and found that the growth rate obtained by solving the Boltzmann equation is a factor of 2 larger than that based on the QFT. Thus, the occupation number from the Boltzmann equation may become exponentially larger than that from the QFT, and a naïve use of the conventional Boltzmann equation may result in a significant overestimation of the number density of χ. 3 In this paper, we have considered the production of a bosonic particle, concentrating on the lowest resonance band of the parametric resonance. Consideration of the production of fermionic particles and the study of the higher resonance bands, as well as the use of the evolution equations based on the QFT to other phenomena, are left as future works [39].