Real-time dynamics of axion particle production due to spontaneous decay of a coherent axion field

We show that the coherent axion field spontaneously decays by emitting axion particles via quantum tunneling if axion potential has more than one minimum. By employing the $\hbar$-expansion and mean field approximation, we develop a formalism to trace the real-time dynamics of the axion particle production including its backreaction to the coherent axion field. We also present numerical results for the time evolution of various physical quantities including the strength of the coherent axion field, phase-space density of produced axion particles, energy density, and pressure. Phenomenological implications of our results are also discussed.


I. INTRODUCTION
The Peccei-Quinn mechanism is one of the most compelling resolutions to the strong CP problem [1,2]. The idea is that the θ-angle of quantum chromodynamics (QCD) is promoted to be a dynamical field θ →φ by introducing an axial U(1) symmetry (Peccei-Quinn symmetry). Low energy dynamics of the fieldφ is determined by the topological structure of the QCD vacuum, because of which an effective potential for the field has a stable minimum at φ = 0 [3]. Therefore, the vacuum expectation value of the field φ dynamically evolves toward the minimum φ = 0 no matter what the initial expectation value is, and the strong CP problem is resolved.
In the present paper, we propose a novel axion particle production mechanism by the spontaneous decay of coherent axion field induced by quantum tunneling. Namely, we consider a situation in which coherent axion field is excited φ = 0 at some instant. Such a coherent excitation may be seeded by, for example, the misalignment mechanism [31][32][33], QCD sphalerons (e.g. the cosmological QCD phase transition, heavy-ion collisions), and classical electromagnetic fields with nonvanishing E · B (e.g. solar flares). As the QCD vacuum has 2π-periodicity in terms of the θ-angle, the axion potential would have infinitely degenerate minima in terms of the fieldφ. This implies that if one considers quantum fluctuations on top of the coherent field, ϕ ≡φ − φ , the coherent field initially peaked at some value φ = 0 would be delocalized because of quantum tunneling between different minima and the quantum fluctuation |φ| 2 grows. This is equivalent to saying that the coherent axion field is unstable against quantum fluctuations because of quantum tunneling, and it decays spontaneously by emitting axion particles. Notice that the existence of many minima is important here. For a single-minimum potential, quantum fluctuations cannot grow and the coherent axion field is stable because that quantum tunneling does not take place. The physical mechanism is essentially the same as the particle production mechanism in the presence of classical field such as Hawking radiation [34,35], Schwinger mechanism [36][37][38], and Landau-Zener transition [39][40][41][42], in which classical field becomes unstable and decays spontaneously by emitting particles because of quantum tunneling.
We shall discuss the aforementioned axion production mechanism based on quantum field theory (QFT). Field theoretical study of real-time dynamics of particle production in the presence of classical field has been developed significantly over the past decades. In particular, in the context of the Schwinger mechanism, not only the dynamical evolution of the production number, but also the backreaction by the produced particles was elegantly formulated based on QFT within mean field approximation by using a proper ultraviolet (UV) regularization scheme [43][44][45][46][47]. Inclusion of backreaction is essential for the energy conservation of the system and for describing the decay dynamics of coherent axion field. The mean field approximation neglects scattering between produced particles (which is responsible for mode-couplings). Therefore, the use of the mean field approximation is justified when coupling constants are sufficiently small [48], which is actually the case for axion. More precisely, the product between the phase-space density f and the coupling constant g determines whether the scattering effects beyond the mean field approximation are important. This is intuitively because scattering takes place more frequently if there are more particles. For the axion self-coupling constant g ∼ χ/f 4 a 10 −40 , where χ ∼ Λ 4 QCD is the topological susceptibility of QCD, huge phase-space density f 1/g ∼ 10 40 is required for the scattering effects to be important. For such huge phase-space density, the axion dynamics would be analyzed classically within, for example, the lattice simulation technique [49]. In the present study, however, we shall explicitly show that the phase-space density of produced axions via the present axion production mechanism is far from dense f / 1/g, and therefore our mean field treatment is justified.
The present paper is organized as follows: In Sec. II, we formulate the axion production mechanism based on QFT within the -expansion and mean field approximation. In Sec. III, we present numerical results based on our formulation and discuss the real-time dynamics of the axion particle production mechanism quantitatively. Section IV is devoted to summary and discussion.
Throughout this paper, we use the mostly minus metric g µν = diag(+1, −1, −1, −1), and work in the Heisenberg picture. We adopt the natural units, and is not written explicitly in the following as long as it is not confusing.

II. FORMALISM
In this section, we explain the theoretical basis of this study. In Sec. II A, in order to describe the axion particle production mechanism from coherent axion field, we employ the -expansion and mean field approximation to derive the equations of motion for the coherent axion field and dynamical axion field on top of the coherent field. In Sec. II B, we explain how to formulate axion particle production with the derived equations of motion. In this study, we adopt the adiabatic particle picture, and canonically quantize the dynamical axion field at each instant of time. We, then, directly evaluate the vacuum expectation value of the axion number operator, and see that it becomes nonzero if axion potential has more than one minimum. In Sec. II C, we discuss the energy balance of the system to check the consistency of our formulation and to see the importance of the backreaction for the energy conservation and describing the decay dynamics of the coherent axion field.

A. Equations of motion
We consider an axion Lagrangian given by whereφ is the axion field and V is the axion potential.
In this paper, we assume for simplicity that the system is not expanding and is homogeneous in space.
Since we are interested in axion particle production from coherent axion field, it is convenient to separate the total axion fieldφ into dynamical axion fieldφ and coherent axion fieldφ aŝ φ =φ +φ. ( Here, the coherent fieldφ is defined as a classical expectation value of the total field operatorφ for a given in-state |in asφ The coherent axion field is a classical quantity and we assume that it is of orderφ = O(1). On the other hand, the dynamical axion field is a quantum fluctuation, so we assumeφ = O( 1/2 ). With this power counting rule, one may expand the Lagrangian (1) in terms of (orφ) as From Eq. (4), one can derive the equation of motion up to O( 3/2 ) as In the following, we ignore the O( 3/2 ) term 2 . By taking the in-in expectation value of Eq. (5) (mean field approximation), one obtains equations of motion for the coherent fieldφ and the dynamical fieldφ as In computing the third term, one has to be careful about how to regularize the UV-divergence in the twopoint function in|φ 2 |in → in| :φ 2 : |in . Note that we will give a precise definition of the regularization :• : later in Sec. II C 1 (see Eq. (25)) based on the adiabatic regularization scheme [50], which can be understood as an analog of the Pauli-Villars regularization scheme in usual QFT in equilibrium. For the moment, we just use :• : to mean a regularized value of • . We also note thatφ 2 − in| : ϕ 2 : |in is higher order in the -expansion compared toφ 2 and in| : ϕ 2 : |in . This is because the expectation value of the difference in| :φ 2 − in| : ϕ 2 : |in : |in is vanishing, and it can give nonvanishing contribution after, at least, it is squared in| : (φ 2 − in| : ϕ 2 : |in ) 2 : |in = 0, which is higher in . Therefore, as long as the -expansion is justified, one can safely neglect the termφ 2 − in| : ϕ 2 : |in in Eq. (6b). An important point in Eqs. (6a) and (6b) is that the dynamical axion field and the coherent field are coupled with each other when the axion potential V is higher than quadratic inφ, i.e., the axion potential V has more than one minimum. This point is essential for the axion particle production as we shall discuss in the next section II B.
The third term in Eq. (6a) represents backreaction of axion particle production to the coherent field. This term is higher in the -expansion than the other terms, however, one should not dismiss this term because it is crucially important for the energy conservation of the system. Indeed, we shall explicitly see later that without this term axion particles are endlessly created so that the energy eventually diverges (see Sec. III B).
Before passing, it is instructive to discuss more about the validity of the -expansion in terms of the size of the phase-space density f and the axion self-coupling con- where m is the axion mass and f a is the axion decay constant. As will be shown later (see Eq. (48)), the two-point function in| :φ 2 : |in 2 For a single-minimum potential with higher order nonlinear terms, higher order -corrections become important. For example, if we consider a quartic potential V ∝ φ 4 , the 2 -correction effectively gives a mass term ∼ φ 2 φ 2 for the Lagrangian of the dynamical field. This term gives a strong negative feedback to the dynamical field to grow, i.e., axion particle production is strongly suppressed. Physically, this is because a singleminimum potential does not permit quantum tunneling and the dynamical field cannot go far away from the minimum. and the total number of produced axions N are related with each other as in| :φ 2 : |in ∼ N/(mV ). Then, by assuming that the typical momentum scale of produced axions is given by the axion mass scale m as d 3 p ∼ (4π/3)m 3 ∼ m 3 (which is actually consistent with our numerical simulations; see Figs. 4 and 9), one may estimate the magnitude ofφ in terms of the phasespace density f ≡ (2π) 3 d 6 N/dx 3 dp 3 ∼ 10 2 × N/(m 3 V ) as By noting thatφ n is always accompanied by V (n) ∼ gf −n a in the -expansion in Eq. (4), one may understand that the -expansion is equivalent to an expansion in terms of gf . Therefore, the -expansion (4) is justified as long as the condition, is satisfied. In other words, the -expansion becomes invalid if the phase-space density becomes extremely huge as f 100/g 10 42 . We shall show explicitly that this is not the case for our axion production mechanism (see Figs. 8 and 9).

B. Axion particle production
We shall solve the coupled equations (6a) and (6b) as an initial value problem to trace the dynamical evolution ofφ andφ. In this subsection, we first explain how to obtain information of (or define) axion particles from the field operatorφ on the basis of the adiabatic particle picture [50]. We, then, show that dynamical axion particles are spontaneously produced from the coherent axion field φ = 0 if the axion potential has more than one minimum.

Adiabatic particle picture
In order to explain the adiabatic particle picture, let us first recall how to define "particle" in a system in equilibrium with the standard canonical quantization procedure. First, we expand the dynamical field operatorφ in terms of a mode function ϕ where we normalized the mode function ϕ (eq) p as 1 = +iϕ (eq) * where In general, the choice of ϕ (eq) p is not unique. That is, one can choose arbitrary ϕ (eq) p as long as it satisfies the equation of motion. However, an equilibrium system should be invariant with respect to time translation. Therefore, it is natural to respect this invariance, and so one can uniquely identify the mode function ϕ (eq) p as an eigen-function of the invariance (i.e., the plane wave) as where axion mass m is defined as the value of √ V in equilibrium and is a positive constant; otherwise the system cannot be equilibrated. After the mode expansion, we impose the canonical commutation relation ontoφ, which in turn quantizesâ (eq) p to define particles. Because of the normalization condition (10), the annihilation operator satisfies Also, the annihilation operatorâ (eq) p can be expressed in terms of the mode function ϕ It is evident from this expression that a different mode function defines a different annihilation operator. We stress that, in equilibrium, one can uniquely define the annihilation operatorâ (eq) p because one can uniquely identify the mode function ϕ (eq) p thanks to the timetranslational invariance. In nonequilibrium, however, the time-translational invariance is explicitly broken so that one has to adopt (or assume) another guiding principle to identify a mode function to define an annihilation operator.
In the adiabatic particle picture, we assume that a given nonequilibrium system is sufficiently adiabatic 3 . 3 The "adiabaticity" is guaranteed if times time derivatives of ωp divided by appropriate powers of ωp (i.e., n ω (n) p /ω n+1 p (n ∈ N)) is small because a time derivative appears with in the equation of motion. This condition is necessary to justify the use of the lowest order WKB solution (14) as the adiabatic mode function ϕ  becomes very large such that it compensates the smallness of . If this is the case, it could be more appropriate to use a higher order WKB solution as the adiabatic mode function; although this is not the case in the present problem because we are always assuming the smallness of throughout our formulation and ω (n) p remains small as long as the typical energy scale of axion condensate is characterized by the QCD energy scale, which is a reasonable assumption in axion phenomenology. A different adiabatic mode function would not only affect the intermediate particle number, but also would change the time evolution of two-point functions because we adopt the adiabatic regulariza-Then, it may be natural to assume that the "true" mode function should not deviate significantly from that in equilibrium (i.e., the plane wave ϕ (eq) p ). The (lowest order) adiabatic particle picture adopts this assumption as a guiding principle to define a mode function in nonequilibrium. Namely, the adiabatic particle picture approximates the "true" mode function by the lowest order Wentzel-Kramers-Brillouin (WKB [53][54][55]) solution ϕ (ad) p of the equation of motion (6b), where Here, we normalized the adiabatic mode function ϕ (ad) p in the same way as ϕ (eq) p (Eq. (10)). Notice that ω p depends on time t through the time dependence ofφ if V is higher than quadratic. ϕ (ad) p is a natural generalization of ϕ (eq) p as they coincide with each other when the system is in equilibrium, in which ω p (or V ) is merely a constant. In the adiabatic particle picture, we, then, expand the field operatorφ in terms of the adiabatic mode function ϕ (ad) p , and employ the canonical quantization procedure to obtain an adiabatic annihilation operatorâ  is not an exact solution of the equation of motion (6b). This difference results in spontaneous axion particle production as we see below.
tion scheme to eliminate UV-divergence in two-point functions, whose subtraction does depend on the choice of the adiabatic mode function (see Sec. II C). The dependence of the truncation order of the WKB expansion for the adiabatic mode function to the intermediate particle number was previously discussed, for example, in Refs. [51,52]. In particular, Ref. [51] considered an analytically solvable model of the Schwinger mechanism, and proposed that there exists an optimal truncation order, around which effects of the truncation order become less significant and universal particle number close to the exact one can be obtained. Therefore, it could be desirable for large and/or ω (n) p to consider the optimal truncation order for the adiabatic mode function. However, it is, in general, impossible to determine the optimal order without knowing the exact solution. Also, it is known that higher WKB solutions can lead to unphysical time evolution of a system (cf. backreaction problem in the Schwinger mechanism [47]). Furthermore, higher WKB solutions contain negative powers of ωp. This leads to severe infrared divergences for massless (or very light) particles, which is actually the case for axion.

Axion particle production
We shall show that axion particles are spontaneously produced from the coherent axion field if the axion potential has more than one minimum by explicitly computing the in-vacuum expectation value of the number operator.
To this end, let us for simplicity assume that the system is in equilibrium V = const. > 0 initially t ≤ t in . Then, the in-vacuum state |vac; t in can be defined as a state such that it is annihilated byâ (eq) p as 0 =â (eq) p |vac; t in for any p.
On the other hand, by using Eq. (13), one can show that the time evolution of the adiabatic annihilation operator a (ad) p can be expressed in terms of a Bogoliubov transformation asâ (ad) where the Bogoliubov coefficients α p , β p are given by Here, ϕ p is the exact solution of the equation of motion (6b) in the momentum space, with initial conditions Notice that β p = 0 holds only when ϕ p = ϕ (ad) p . In other words, β p = 0 must hold if the axion potential has more than one minimum. Indeed, ϕ p couples to the coherent axion field because of the existence of many minima as explained below Eqs. (6a) and (6b), and hence ϕ p should dynamically evolve to deviate from ϕ (ad) p because ϕ (ad) p is no longer a solution of the equation of motion. Note also that the Bogoliubov coefficients always satisfy which is a consequence of the unitarity 4 . 4 A proof is the following: Let us define S-matrix as a matrix such It is evident from Eq. (17) that the in-vacuum state |vac; t in is no longer annihilated by the annihilation operatorâ (ad) p for t > t in . Thus, the in-vacuum expectation value of the number operator vac; t in |â (ad) † pâ (ad) p |vac; t in is nonvanishing, which implies that axion particles are spontaneously produced from the coherent axion field. Indeed, one can explicitly compute the expectation value as 5 where δ 3 (p = 0) = V /(2π) 3 was used.

C. Energy conservation
Finally, let us discuss the energy balance of the system to check the consistency of our formulation and also the axion particle production from the viewpoint of energy conservation. To this end, we first explain how to regularize the UV-divergence of two-point functions. In this paper, we adopt the adiabatic regularization scheme [50], which is compatible with the adiabatic particle picture. Then, we compute the expectation value of the energymomentum tensor, and show how the total energy of the system is conserved. We shall explicitly see that the energy of the coherent axion field is converted into axion particles, and the backreaction term in Eq. (6a) is crucial for the energy conservation.

Adiabatic regularization
The UV-divergence of in-in expectation values of twopoint functions originates from the uninterested vacuum contribution and must be subtracted by a regularization procedure. Usually, the vacuum is a stable state (i.e., does not decay spontaneously), so that the divergent vacuum contribution is constant in time. However, in the presence of a time-depending external field, the vacuum state is, in general, no longer stable and it decays spontaneously because of the particle production. This implies that the vacuum contribution cannot be constant in time (cf. the same problem appears in the real-time dynamics of the Schwinger mechanism [43][44][45][46][47]). Indeed, from the time dependence of the adiabatic annihilation operator (17), one may define the vacuum at time t, which we write |vac; t , as a state such that 6 0 =â (ad) p (t) |vac; t for any p.
Using |vac; t , one may identify the vacuum contribution at time t of some two-point functionφΓφ as which obviously depends on time. Note that it is impossible to express the vacuum contribution at time t (24) as an expectation value in terms of the in-vacuum |vac; in , which does not know the vacuum structure at time t.
In the adiabatic regularization scheme, we subtract the vacuum contribution (24) from the bare invacuum expectation value vac; t in |φΓφ|vac; t in , and define the regularized in-vacuum expectation value vac; t in | :φΓφ : |vac; t in as We remark that it is inevitable in the adiabatic regularization scheme to evaluate an expectation value in terms of the adiabatic vacuum state |vac; t in order to compute an in-in expectation value in a well-defined/consistent manner.
Let us see explicitly how Eq. (25) eliminates the UVdivergence. For this purpose, let us consider Γ = 1, i.e., ϕΓφ =φ 2 as an example. We first identify the divergent structure of the bare value vac; t in |φ 2 |vac; t in by solving the equation of motion (19) perturbatively in terms of 1/ω p . To do this, it is convenient to make a WKB-type ansatz for ϕ p as Then, we substitute this expression into the equation of motion (19) to find 6 Note that |vac; t is a Heisenberg state, so that it is time independent, i.e., |vac; t (t 1 ) = |vac; t (t 2 ) for any t 1 , t 2 . Nevertheless, as |vac; t is defined as an eigen-state of the timedepending adiabatic number operatorâ . In this sense, one could say |vac; t looks time-dependent although it is a Heisenberg state. This equation can be solved perturbatively in terms of 1/ω p 7 as Notice that the first term coincides with the lowest order WKB solution, which is nothing but the adiabatic mode function ϕ (ad) p by definition. From Eq. (28), one can identify the divergent structure of the bare expectation value vac; t in |φ 2 |vac; t in as It is evident that only the first term in the brackets is divergent. In other words, only the contribution from the lowest order WKB solution is UV-divergent, and is always identical to the vacuum contribution 8 Therefore, we have in| : Thus, the regularized expectation value in| :φ 2 : |in defined by the adiabatic regularization scheme (25) is always finite and well-defined. Note that in| :φ 2 : |in is free from logarithmic divergences. This is because we dropped higher order terms in because of the weakness of the axion coupling constant (see Sec. II A), i.e., interactions among dynamical axions are neglected (i.e., no loop corrections to mass and/or coupling constants). If one includes higher order terms in , which are responsible for interactions, a logarithmic divergence would appear (cf. V = λφ 4 case was discussed in Ref. [60]), which can be eliminated by a renormalization procedure [43][44][45][46]60]. Also, remark that the adiabatic regularization scheme does not violate the energy conservation law as we will explicitly check below. The subtraction in nonequilibrium QFT should be time-dependent, so that it is nontrivial whether the subtraction is consistent with conservation laws. This is not a specific issue in the adiabatic regularization scheme, but is equally nontrivial in other regularization schemes. Hence, it is important to check conservation laws explicitly so as to ensure that the regularization scheme is not unphysical.

Energy and pressure
Now, we evaluate the expectation value of the energymomentum tensor T µν within the adiabatic regularization scheme. The (symmetric) energy-momentum tensor T µν is defined aŝ where higher order quantum corrections O( 3/2 ) are neglected in the last equality so as to be consistent with Eq. (5).
The energy density operatorÊ and the pressure oper-atorP are defined as the diagonal components ofT µν . Their expectation value can be evaluated within the adiabatic regularization scheme as ≡ vac; tin|:Êφ:|vac; tin Note that ∂φ = 0 andP =T xx =T yy =T zz because we assumed spatial homogeneity. Let us check the finiteness of the expressions (33) and (34). By using Eq. (28), one finds vac; t in | :Ê : and vac; t in | :P : Therefore, the energy density vac; t in | :Ê : |vac; t in and the pressure vac; t in | :P : |vac; t in are indeed finite. Note that there is no logarithmic divergence because there are no interactions among dynamical axions within our order of the -expansion as explained in Sec. II C 1.

Energy conservation
Let us explicitly check the total energy conservation of the system. The time derivative of the energy density of axion particles can be evaluated as where the use is made of the equation of motion for ϕ p (19). On the other hand, the time derivative of the energy density of the coherent axion field is given by Here, we used the equation of motion (6a) to reach the second line. Therefore, by combining Eqs. (37) and (38), we find that the total energy is strictly conserved as 0 = ∂ t vac; t in | :Ê : |vac; t in .
It is now evident from Eq. (38) that if the energy of axion particles increases because of the spontaneous axion particle production, the same amount of energy of the coherent axion field must decrease. That is, the energy of the coherent axion field is the source of axion particle production. Notice that the right-hand side of Eq. (38) vanishes if one neglects the backreaction term (which is higher in ) in the equation of motion for the coherent field (6a). In this case, axion particles are endlessly created and the energy eventually diverges because the coherent axion field supplies energy to the axion particle production endlessly without any energy loss. This is clearly unphysical. Therefore, it is important to maintain the higher order 2 -term in the equation of motion (6a) when discussing quantum corrections to the coherent field.

III. NUMERICAL RESULTS
We numerically solve the equations of motion (6a) and (6b), and discuss quantitatively the real-time dynamics of the axion particle production from coherent axion field.

A. Setup
In order to solve the equations of motion (6a) and (6b), we need to specify the axion potential V and initial conditions forφ andφ.
For the potential V , we consider a periodic potential with infinitely degenerated minima described by Here, χ is the topological susceptibility of QCD, f a is the axion decay constant, and θ is the intrinsic θ-angle.
Notice that one of the stable minima of the potential is φ = θf a , at which holds.
That is, we consider a situation such that the coherent axion fieldφ sitting at one of the minima is perturbed suddenly at time t in . Such a perturbation may be seeded by, for example, the misalignment mechanism [31][32][33] 9 . QCD sphaleron transitions (e.g. the cosmological QCD phase transition, heavy-ion collisions), and classical electromagnetic fields with nonvanishing E · B (e.g. solar flares). The magnitude of ξ determines the initial (kinetic) energy of the coherent field. In this paper, for the sake of simplicity, ξ is assumed to be a constant which does not depend on x, i.e., the system is spatially homogeneous. We also assume ξ is smaller than the QCD energy scale ξ χ 1/2 ∼ Λ 2 QCD because the axion potential (39) is an effective potential, which is valid only below the QCD energy scale.
The initial condition for dynamical axion fieldφ is already specified in Sec. II B 2. We impose that the field operatorφ is expressed in terms of the plane wave at time t = t in as where ϕ (eq)

B. Without backreaction
In this subsection, we artificially neglect the backreaction term in the equation of motion for the coherent axion field (6a). Although this treatment is unphysical because it violates the energy conservation law as explained in Sec. II, it is instructive to understand some basic features of the axion particle production mechanism. In particular, we numerically show that axion particles are endlessly created, no matter how weak the initial energy ξ is. This implies that a coherent axion field configuration φ = 0 is unstable with respect to quantum fluctuations, 9 Precisely, the misalignment mechanism is initiated with finitē φ − faθ = 0 and vanishingφ = 0. For this case, the coherent axion field initially has potential energy, instead of kinetic energy. This difference is essential neither in our production mechanism nor in the evolution ofφ. The result is basically the same as long as the initial total energy is the same. and thus inclusion of quantum -corrections is inevitable in understanding the real-time dynamics of axion. Also, an analytical calculation based on the standard perturbation theory is doable in the absence of the backreaction term (see Appendix A). We compare the numerical results with the analytical perturbative calculation and see that the present axion particle production mechanism is genuinely nonperturbative which cannot be described within the perturbation theory because it is caused by the nonperturbative quantum tunneling. Figure 1 shows the time evolution of the coherent axion fieldφ. The coherent fieldφ oscillates in time with (almost) constant frequency and amplitude. This is simply because we assumed that the initial kinetic energy of the coherent field is smaller than the potential height ξ 2 χ. Therefore, it is classically forbidden for the coherent field to climb up the potential hill to move to a next minimum, so that it oscillates around the initial minimum just like a harmonic oscillator. Mathematically, one may expand Eq. (6a) in terms of (φ −φ| t=tin )/φ| t=tin 1 to get

Coherent axion field
One immediately gets a solution to this equation as The constant amplitude implies that the coherent fieldφ never loses energy and endlessly supplies energy to the dynamical fieldφ because of the absence of backreaction. Note that, in expanding systems, there appear additional terms in Eq. (44), because of whichφ decays in time.  Figure 2 shows the time evolution of the variance of the axion field (or the magnitude of the dynamical axion field) vac; t in | :φ 2 : The variance vac; t in | :φ 2 : |vac; t in diverges exponentially for any values of ξ. That is, although the coherent axion fieldφ oscillates just around the initial minimum because of the classical restriction, the quantum fluctuationφ can go far away from the minimum because of quantum tunneling. The fluctuation eventually extends to infinity because the axion potential (39) has infinite number of minima, which are distributed infinitely in the field space. We stress that quantum tunneling plays an essential role here. Indeed, the fluctuation cannot grow if there exists only one minimum, which is in contrast to the axion particle production via parametric resonance due to oscillating coherent axion field [23][24][25][26].
This quantity vac; t in | :φ 2 : |vac; t in is directly related to the total number of produced axions In other words, if vac; t in | :φ 2 : |vac; t in increases because of quantum tunneling, N/V should increase, i.e., axion particles are produced. Indeed, by using one can re-express vac; t in | :φ 2 : |vac; t in in terms of the Bogoliubov coefficient β p as Note that the second term has an oscillating factor in time, whose frequency is ∼ m 0 . This is why the curves in Fig. 2 look like a band. The first term is a classical contribution because |β p | 2 equals to the phase-space Time evolution of the phase-space density (2π) 3 d 6 N/dp 3 dx 3 = |βp| 2 of axion particles without backreaction. The initial energy ξ is fixed as ξ/ √ χ = 0.5. density (2π) 3 d 6 N/dp 3 dx 3 (22). The second term arises because of a quantum interference between the positive and negative frequency modes, and gives a quantum correction to the classical value. Usually, the second term gives only a small contribution compared to the first one after the momentum integration because the integrand contains severely oscillating factor in momentum ∼ exp −i ω p dt . As we shall see soon, produced axions are typically soft |p| m 0 (see Fig. 4), for which energy cost to produce becomes smaller. Thus, it is good to approximate ω p in Eq. (47) as ω p ∼ m 0 . Then, after neglecting the quantum correction, one obtains vac; t in | :φ 2 : |vac; 3. Axion particle production Figure 3 shows the time evolution of the total number of produced axion particles N/V = d 3 p|β p | 2 /(2π) 3 . N/V increases exponentially for any value of ξ. This implies that the coherent axion field is unstable against quantum fluctuations and it spontaneously decays by producing axion particles no matter how small or large the initial fluctuation ξ is.
The time evolution of the phase-space density (2π) 3 d 6 N/dp 3 dx 3 of produced axion particles is shown in Fig. 4. Low momentum |p| m 0 axion particles dominate the production. This is a natural result because it is energetically costful to create hard particles.
It is instructive to compare these results with a perturbative calculation to clarify the nonperturbative nature of our production mechanism: if the initial energy is sufficiently small compared to the QCD energy scale ξ √ χ, the coherent axion field does not deviate significantly from the initial minimum δφ ≡ (φ −φ t=tin )/φ t=tin 1 (see Fig. 1). Then, one may naively expect that the axion particle production rate can be evaluated perturbatively in terms of δφ (or equivalently ξ). As the coherent field φ decouples from the dynamical fieldφ in the absence of the backreaction term, one can carry out the perturbative calculation analytically as (see Appendix A for the detail) and An important point here is that the perturbative production number N (pert) is proportional to t. Apparently, this is inconsistent with the exponential growth observed in Figs. 3 and 4, and the perturbation theory agrees with our results only at very early times. This disagreement is a natural result because our production mechanism is a consequence of quantum tunneling (as we discussed in Sec. III B 2), in which the existence of many minima is important. By the perturbative expansion in terms of δφ, only the local structure of the axion potential around one minimum is taken into account, and the global information of the axion potential (i.e., the existence of many minima) is lost. Figure 5 shows the time evolution of the total energy density vac; t in | :Ê : |vac; t in and total pressure vac; t in | :P : |vac; t in . As expected, the energy and the pressure diverge because of the divergent contribution from axion particlesφ. Therefore, one cannot dismiss the backreaction term to describe real-time dynamics of axion.

Energy and pressure
The energy density is dominated by the mass energy of produced axion particles because they are soft (see Fig. 4). To see this analytically and to find out relationship to other observables such as N and the variance, let us rewrite vac; t in | :Ê ϕ : |vac; t in in terms of the Bogoliubov coefficient. By using Eq. (33), one obtains vac; t in | :Ê ϕ : Then, by approximating ω p ∼ m 0 , we find On the other hand, quantum corrections dominate the pressure.
To show this, let us rewrite vac; t in | :P ϕ : |vac; t in (see Eq. (34)) in terms of the Bogoliubov coefficient as As we know that produced axion particles are soft (see Fig. 4), it is good to approximate ω p ∼ m 0 . Then, only the second term in Eq. (53), which is nothing but quantum corrections, survives as in| : Equation (54) oscillates in time, which is the reason why we see an oscillating behavior in Fig. 5.

C. With backreaction
Let us recover the backreaction term, and study the role of the energy conservation. Below, we explicitly see that because of the energy conservation, (i) dynamics of the coherent fieldφ is no longer independent of the dynamical fieldφ; (ii) the variance cannot extend to infinity, i.e., vac; t in | :φ 2 : |vac; t in < ∞; (iii) the number of produced axion particles stays finite; and (iv) not only the axion particle production process, but also an axion particle pair annihilation process takes place. Figure 6 shows the time evolution of the coherent field φ. Since the energy of the coherent field is consumed in the axion particle production, the amplitude ofφ decreases from the initial value. However, the amplitude never goes to zero even at late times. In other words, the coherent field never falls into a minimum of the potential (i.e., no equilibration), but it oscillates around a minimum forever if the system is non-expanding. This is because axion particle pair annihilation occurs. Indeed, because of the backreaction term, not only the axion particle production process, but also its inverse process, i.e., pair annihilation process (2φ →φ) can occur. Due to the pair annihilation process, the time evolution of the coherent field is modified as follows: For the first moments, the coherent field loses its energy (i.e., amplitude decreases) and the phase-space density of axion particles increases because of the axion particle production process. As the phase-space density increases, produced axion particles begin to overlap with each other and the pair annihilation process becomes more efficient than the production process. Once the pair annihilation process dominates, the phase-space density decreases and the coherent field would obtain energy from axion particles (i.e., amplitude increases) because of the energy conservation; see Eq. (38). After some time, the phasespace density becomes dilute and the pair annihilation process should be suppressed. Then, the production process would again be more efficient. In this way, the dominant process changes as time goes, and hence the coherent field never gets equilibrated. Notice that we can see small oscillations in the amplitude at late times in Fig. 6. This is a manifestation of the above mechanism. Also, in Sec. III C 3, we shall discuss the time evolution of the phase-space density. There, we can explicitly see that the annihilation process and the production process, indeed, alternatively take place. Note that the same oscillation mechanism was already discussed in the real-time dynamics of the Schwinger mechanism [45][46][47], in which classical electric field shows an oscillating behavior (plasma oscillation). Note also that in order for the coherent field to completely fall into a minimum, the energy of the total axion fieldφ must be dissipated by, for example, coupling to other particles and attaching to heat bath.

Coherent axion field
The frequency of the oscillation deviates slightly from m 0 because of the backreaction term. Indeed, by assuming that vac; t in | :φ 2 : |vac; t in is sufficiently adiabatic, one can expand the equation of motion (6a) in terms of (φ −φ| t=tin )/φ| t=tin 1 as Thus, the frequency is given, roughly, by m 0 × 1 − vac; t in | :φ 2 : |vac; t in /2f 2 a < m 0 . Intuitively, produced axion particles act as friction to the coherent field, which lowers the frequency of the oscillation.  the number of produced axion particles is small and the backreaction effect can be neglected, the variance grows exponentially as we found in Fig. 2. However, after sufficient number of axion particles is produced, the backreaction effect becomes important and the variance stops growing. Thus, as expected, the variance never diverges, i.e., the quantum fluctuationφ cannot go far away from the initial minimum because of the energy conservation. Note that the variance shows oscillating behaviors at late times because of the pair annihilation process as explained in Sec. III C 1. Figure 8 shows the total number N/V = d 3 p|β p | 2 /(2π) 3 of produced axion particles. It is evident that the total number stays finite because of the energy conservation. The time dependence is basically the same as the variance vac; t in | :φ 2 : |vac; t in (see log(d 6 N/dp 3  Time evolution of the phase-space density (2π) 3 d 6 N/dp 3 dx 3 = |βp| 2 of axion particles with backreaction. The initial energy ξ is fixed as ξ/ √ χ = 0.5. Fig. 7) because they are closely related with each other as shown by Eq. (48). The only difference is that we cannot see the "band structure" in this figure because the variance receives quantum corrections (the second term in Eq. (47)), which oscillate severely in time. Figure 9 shows the time evolution of the phase-space density (2π) 3 d 6 N/dp 3 dx 3 of produced axion particles. The phase-space density is significantly modified by the backreaction (see Fig. 4). In particular, the pair annihilation process plays an important role: The phase-space density is softened by pair annihilation of hard axion particles. Also, we can see a striped pattern in the figure, which is because the annihilation process and the production process alternatively take place as explained in Sec. III C 1.

Axion particle production
With our parameter choice ξ/ √ χ = O(0.1), i.e., the initial energy of the coherent axion field is at most the QCD energy scale, the axion particle production is completed, roughly, within t prod ∼ O(1000) × m −1 0 . The current axion mass bound is 10 −5 eV m 0 10 −2 eV, so that the time-scale can be estimated as 10 5 eV −1 t prod 10 8 eV −1 (or 10 −10 sec t prod 10 −7 sec). Figure 10 shows the time evolution of the energy density vac t in | :Ê : |vac t in and the pressure vac t in | :P : |vac t in . The energy density is strictly conserved because of the backreaction. With our parameter choice ξ/ √ χ = O(0.1), about one-half of the initial energy of the coherent axion field is eventually converted to axion particles. On the other hand, the coherent fieldφ dominates the total pressure. The produced axion particlesφ only give a relatively small contribution because axion particles are soft.

IV. SUMMARY AND DISCUSSION
We discussed the axion particle production from coherent axion field. We found that the coherent axion field is unstable against quantum fluctuations and it spontaneously decays by emitting axion particles when the axion potential has more than one minimum, for which quantum tunneling between different minima takes place.
In order to study the above axion particle production mechanism, we developed a formalism which describes the real-time dynamics of the axion particle production. Namely, we derived equations of motion for coherent axion field and dynamical axion field on top of the coherent field by employing the -expansion and mean field approximation. The equations are coupled with each other when axion potential is higher than quadratic, i.e., axion potential has more than one minimum. Because of this coupling, the dynamical axion field obtains energy from the coherent field, which eventually results in the axion particle production. We included not only 1 -terms, but also the next leading order 2 -term in the equations of motion in order to take into account backreaction by the axion particle production to the coherent axion field. The backreaction is essential for the energy conservation of the system and for describing the decaying dynamics of the coherent field.
We solved the equations of motion numerically to trace the real-time dynamics of the axion particle production quantitatively. By adopting the adiabatic particle picture and canonically quantizing the dynamical axion field at each instant of time, we evaluated the phase-space density of the produced axion particles as a function of time. We also computed the energy density and the pressure of the system, and followed their time evolution. As a result, we found, in particular, that (i) soft axion particles are produced from the coherent axion field, whose number grows exponentially until the production stops because of the energy conservation; (ii) the production number is much more abundant than usual perturbative axion production mechanisms, which are suppressed by the large decay constant f a , and, roughly, one half of the initial energy of the coherent axion field is eventually converted into axion particles; (iii) not only axion particle production, but also an axion particle pair annihilation process occurs due to the backreaction, which leaves a characteristic striped pattern in the phase-space density; (iv) the coherent axion field loses its energy through the axion particle production, but does not completely fall into a minimum because of the pair annihilation process; and (v) inclusion of the backreaction term is important, without which the system evolves unphysically and various physical quantities including the energy density diverge. In addition, we also compared our results with the perturbation theory. We found that our axion particle production mechanism is genuinely nonperturbative, which cannot be described within the perturbation theory even if the initial energy of the coherent field is very small. This is because the quantum tunneling is the phys-  for various values of initial energy ξ with backreaction. The left panel shows the total value of the system, and the center and the right panels show the contribution from the coherent axion field and the produced axion particles, respectively.
ical origin of our production mechanism. There are several interesting applications of the present study. One example is an application to cosmology, in particular, to axion dark matter scenario. In the standard axion dark matter scenario, the misalignment mechanism [31][32][33] is expected to produce coherent axion field in the early Universe, whose energy density can be comparable to that of dark matter in the present Universe. Since coherent fields have vanishing momentum, it would be a good candidate for cold dark matter. According to our results, coherent axion field (or axion condensate) is unstable against quantum fluctuations and it decays quickly by producing "warm" axion particles, whose momentum is not so large |p| m 0 but nonvanishing. In addition, the momentum distribution has a characteristic striped pattern due to the annihilation process. These features should affect the small-scale structure of the present Universe, which might be tested by future experiments/observations. Note that we considered a nonexpanding system in this paper for simplicity. Inclusion of expansion is important to make quantitative prediction for cosmology because the production process should be affected significantly by expansion (e.g. red-shift in momentum distribution).
Ultra-relativistic heavy-ion collisions are another example. This is because heavy-ion collisions can produce sizable number of QCD sphalerons. QCD sphalerons would supply energy of the order of the QCD energy scale Λ QCD to coherent axion field, which in turn decays into axion particles through our production mechanism. As an example, let us consider Pb-Pb collisions at the LHC. In each heavy-ion collision, a quark-gluon plasma (QGP) with temperature T ∼ 200 MeV, vol-ume V ∼ 10 4 fm 3 , and lifetime τ ∼ 10 fm is produced. As the QCD sphaleron rate Γ is given by Γ ∼ 30α 4 s T 4 [61], we may estimate the number of sphalerons in a QGP created in a heavy-ion collision as ΓV τ ∼ 10 4 for α s ∼ 0.3. The rate of collisions is of the order 1 kHz, and thus, roughly, 10 7 sphalerons would be produced per second in heavy-ion collisions. On the other hand, we found that about one-half of the energy of coherent axion field is eventually converted into axion particles (see Fig. 10). Therefore, one QCD sphaleron may produce 0.5 × Λ QCD /m 0 ∼ 10 +12∼+15 very soft axions. Combining the sphaleron production number, we estimate that N ∼ 10 +19∼+22 axion particles may be produced in a second in heavy-ion collisions. This is a very large number compared to usual perturbative axion production processes (e.g., solar axion flux on the Earth is the order of 10 9 cm 2 s −1 [62]), which is always suppressed strongly by the large decay constant f a . Therefore, heavy-ion collisions may be another good source of axion particles. Note, however, that spatial homogeneity is assumed in the present paper. Since a QCD sphaleron is strongly localized in space with the QCD energy scale, it is very important to consider spatial inhomogeneity to make a more realistic estimation.
It is also interesting to study the emission of gravitational waves induced by the present axion particle production mechanism. Emission of gravitational waves by production of axion particles via parametric resonance was previously discussed in Refs. [24][25][26]. We expect a similar gravitational wave emission occurs from our production mechanism as well. If this is the case, it may leave a characteristic signature in gravitational waves due to the oscillating pattern in the momentum spectrum.
This could serve as a novel way to detect (although indirect) axion particles through observation of gravitational waves. We leave these topics for a future work.
In this appendix, we analytically discuss the axion particle production from coherent axion field within the lowest order perturbation theory by neglecting the backreaction effect.
To this end, we consider a Lagrangian for dynamical axion fieldφ given by where O( 3/2 ) terms are neglected (see Eq. (6b)). Here, the coherent fieldφ obeys the equation of motion (6a) without the backreaction term, i.e., Now, let us assume that the coherent fieldφ does not deviate significantly from a minimum of the axion potential V , which we writeφ 0 . Physically, this assumption is equivalent to assume, as we shall see explicitly later, that the initial energy of the coherent field is sufficiently smaller than the height of the potential. We, then, perturbatively expand the Lagrangian (A1) in terms of Here, In the standard S-matrix formalism, the S-matrix is defined as a time-ordered product of the interaction Lagrangian L int as Then, the in-and out-state annihilation operators are related with each other by the S-matrix aŝ a (out) Therefore, in the lowest order in L int , we have For the periodic axion potential (39), the interaction Lagrangian L int starts from n = 2, i.e., where we used Eq. (45) to get the second line. By plugging this expression into Eq. (A8), one obtains The total number N (pert) can be obtained by integrating this expression over p as where we used ∞ 0 d|p|δ(|p|) = 1/2.