Explosive Axion Production from Saxion

The dynamics of saxion in a supersymmetric axion model and its effect on the axion production is studied in detail. We find that the axion production is very efficient when the saxion oscillation amplitude is much larger than the Peccei-Quinn scale, due to a spike-like behavior of the effective axion mass. We also consider the axino production and several cosmological consequences. The possibility of detection of gravitational waves from the non-linear dynamics of the saxion and axion is discussed.


Introduction
Supersymmetric (SUSY) axion models combine the idea of SUSY and the Peccei-Quinn (PQ) mechanism [1] to solve both the hierarchy problem and the strong CP problem. SUSY axion models are also interesting from the cosmological point of view, mainly due to the existence of (relatively) light particles, the saxion and axino, which are the scalar and fermionic partners of the axion, respectively .
The cosmological saxion dynamics is highly non-trivial and its cosmological effects are still not fully understood despite its importance. In a simplified treatment, the saxion is assumed to begin a coherent oscillation at some epoch and eventually it decays into an axion pair as well as the minimal SUSY standard model (MSSM) sector. The produced axions compose dark radiation which is constrained from cosmological observations. In such a scenario, it is often assumed that the saxion perturbatively decays into the axion pair. It is still non-trivial, however, whether the axion production can be treated perturbatively especially when the saxion field value significantly deviates from the potential minimum.
In this paper we revisit the saxion dynamics in a model described by the superpotential (2) shown below to correctly estimate the axion production rate. Ref. [8] partly addressed this issue and found that there is a parametric resonant enhancement of the axion when the saxion amplitude is large enough. We analyze the dynamics in more detail and identify the origin of the resonant axion production. We find a "spike"-like feature of the effective axion mass, which originates purely from multi-field property of the dynamics. Due to the spike, the axion production rate is much more enhanced compared to the naive perturbative calculation of the decay rate. Actually the axion production is so explosive that the whole saxion-axion system becomes completely non-linear after only a few saxion oscillations. It may drastically change the conventional picture of the saxion oscillation and the resulting particle production processes.
We first discuss the saxion dynamics in Sec. 2 in terms of the canonically normalized fields. We estimate the axion and axino production rate in Sec. 3. Several cosmological implications are also mentioned. We conclude with some discussion including implications for gravitational wave observation in Sec. 4.

Saxion dynamics 2.1 Model
We consider the axion model in which the Kähler and superpotential are given by where φ andφ are PQ superfields with PQ charges +1 and −1 respectively, and X is a PQ-singlet superfield. The constants λ, f and W 0 are taken to be real and positive without loss of generality. The R-symmetry under which only X has a charge +2 ensures this type of superpotential. Including the SUSY breaking effects, the scalar potential is given as where m φ and mφ are soft SUSY breaking masses and m 3/2 = W 0 /M 2 P denotes the gravitino mass with M P being the reduced Planck scale. Below we assume f m φ , mφ, m 3/2 . The potential minimum is Below we redefine m 2 φ + λ 2 X 2 and m 2 φ + λ 2 X 2 as m 2 φ and m 2 φ , respectively. In the SUSY limit m φ , mφ → 0, there are two massless modes, called axion and saxion, corresponding to the flat direction φφ = f 2 . The saxion obtains a mass of the order of the soft SUSY breaking, but still it is much lighter than the PQ scale ∼ f . #1 Thus the saxion dynamics can have significant impact on cosmology. In particular, we are interested in the case that the saxion has an initial value as large as M P . #2 We reconsider the saxion dynamics in such a case, paying particular attention to its effect on the axion production.

Equation of motion
Let us write the PQ fields as The dynamics is almost constrained to the flat direction φφ = f 2 since the scalar degrees orthogonal to it are heavy enough. This constraint is written as The PQ scale f a is given by f 2 a = (ϕ 2 +φ 2 )/N 2 DW with N DW being the domain wall number, which depends on the PQ charge assignments on the MSSM sector.
#2 Actually, if the PQ field obtains a negative Hubble-induced mass squared through the coupling to the inflaton, the saxion gets a VEV of ∼ M P during inflation.
Using this, we can replaceφ and aφ in terms of ϕ and a φ which correspond to the saxion and axion, respectively. Thus we obtain the kinetic term for the saxion and axion as where The scalar potential is where ∆m 2 ≡ m 2 φ − m 2 φ . The equation of motion of ϕ reads as far as the backreaction from a φ is neglected. In the small amplitude limit ϕ = ϕ + δϕ (|δϕ| f ), it is simplified as The equation of motion of a φ reads It may be useful to rewrite it in terms of θ ≡ a φ /ϕ. The kinetic term for θ is given as Then the equations of motion for θ is Note that it has a shift symmetry θ → θ + const.

Canonical saxion and axion
For later convenience, we rewrite the equations of motion in terms of the canonically normalized fields. Let us define the canonical saxion ϕ and canonical axion a φ as Note that The equation of motion of ϕ is Thus the canonical saxion ϕ behaves almost as a harmonic oscillator with its initial value ϕ ini . The equation of motion of the canonical axion a φ is where Using the equation of motion of ϕ, it is rewritten as It is approximated as Thus the axion is massless only if the background is settled to the potential minimum (δϕ = 0 andφ = 0). It obtains effective mass if the background deviates from the minimum and time-dependent. Especially, when the saxion is oscillating with its amplitude much larger than f , the last term in the second line of (25) becomes important as we will see below.

Small saxion amplitude regime
First let us consider the small saxion amplitude case: ϕ amp f , to compare with the conventional calculation of the saxion decay rate. In this case, the term proportional toφ 2 is negligible in the axion mass expression (second line of (25)). Then we can extract the axion production rate from the δϕ dependence of the axion mass term. Approximating δϕ as a harmonic oscillator, we obtain the saxion decay rate into the axion pair as The first term is understood as the usual perturbative decay rate of the saxion into the axion pair around the vacuum. The second term may be regarded as annihilation of the saxion into the axion pair, which exists for finite saxion amplitude. As is known, the perturbative decay rate vanishes in the limit ∆m 2 = 0 [6]. Even in such a case, there is a finite contribution to the axion production from the saxion annihilation represented by the second term, although it cannot lead to the complete saxion decay since it decreases faster than the Hubble rate. However, the axion production is very efficient in the large saxion oscillation regime even in the case of ∆m 2 = 0 as we will show below. In the following we take ∆m 2 = 0 to derive the lower bound on the axion dark radiation abundance.

Large saxion amplitude regime
Let us consider a large saxion oscillation amplitude regime: ϕ amp f . #3 We assume m φ ϕ amp < f 2 since otherwise the saxion dynamics is not confined to the flat direction ϕφ = 2f 2 and the dynamics would be more chaotic, which may lead to the nonthermal PQ symmetry restoration. #4 The saxion oscillation is approximated as Then we have Notice that the axion temporary becomes very massive when the saxion passes through the potential minimum ϕ 0 where m 2 a ∼ m 2 φ ϕ 2 amp /f 2 m 2 φ . This is much different from the case of single PQ field models: in the single PQ field models, the canonical axion mass is always zero when the saxion passes through the potential minimum. A crucial difference between the single and multi PQ field model is that the definition of "axion" becomes timedependent in the multi-field model. Actually one easily recognizes thatã a φ for ϕ f and a aφ for ϕ f . The eigenstate suddenly changes around ϕ f and the time scale for this sudden transition is ∆t ∼ (f /φ) ϕ f ∼ f /(m φ ϕ amp ). This is the reason for the appearance of temporal large mass scale in the effective axion mass. We call this as a "spike"-like behavior of the effective axion mass due to the saxion dynamics. The left panel of Fig. 1 shows time evolution of |m 2 a (t)| for ϕ amp = 100f by numerically solving the equation of motion of the saxion (13) with ∆m 2 = 0.
The particle production rate with such a spike-like mass is calculated using the result shown in App. A. The axion energy density ρ a after one half saxion oscillation is given by where ρ φ ∼ m 2 φ ϕ 2 amp is the saxion energy density. Since we assumed m φ ϕ amp < f 2 , this ratio is smaller than one. Formally we may define a "decay rate" of saxion into the axions as It is clear from this expression that the axion production is significantly enhanced for ϕ amp f . Moreover, one should take account of the effect of parametric resonance, since the occupation number f k of the axion exceeds unity for k m sp ∼ m φ ϕ amp /f . We numerically evaluated the time evolution of the phase space density of the axion. The linearized equation of motion of a φ in the momentum space is given bÿ where m 2 a (t) is given by (24) which is evaluated by numerically solving the equation of motion of saxion (13). The initial condition is taken to be a k = 1/ √ 2ω k and˙ a k = −i ω k /2. #5 The phase space density is given by The right panel of Fig. 1 shows the phase space density of the axion f (a) k as a function of wave number k in units of m φ for 1, 3 and 5 half saxion oscillation (t = π/m φ , 3π/m φ , 5π/m φ ) for ϕ amp = 100f . Hence the spike scale is m sp ∼ 100m φ . It is seen that after one half saxion oscillation, the spectrum shows a k −2 behavior for m φ < k < m sp and cutoff for k > m sp consistent with those studied in Refs. [27,28]. At this stage, the energy density is dominated by the modes k ∼ m sp . However, after a few oscillation, modes with k ∼ m φ is most enhanced and the axion energy density will soon be comparable to the saxion energy density. Then the backreaction will become important and the linearized treatment breaks down. Another important feature is that the low frequency modes k m φ do not grow after a few oscillation despite the fact that the occupation number is much larger than one. #6 This feature is much different from the ordinary broad resonance. To understand this, one may find that Eq. (17) has a solution θ = const. in the limit k → 0, due to the shift symmetry of θ. Actually we numerically checked that a k (t) ∝ ϕ(t) roughly followed for k m φ . It is difficult to rigorously follow the non-linear evolution of the whole system, but it is reasonably expected that the axion backreacts to the saxion condensate so that the most saxion energy density is transferred to k ∼ m φ modes. Then the saxion and axion will have comparable energy density until the mean saxion amplitude decreases to ∼ f due to the Hubble expansion. After that, the axion becomes relativistic (see Eq. (25)) and the axion energy density decreases faster than the saxion. Since we are considering the initial saxion amplitude as large as the Planck scale, it is likely that the saxion dominates the universe at this stage. Eventually the saxion perturbatively decays into MSSM particles: it decays into gluons in the KSVZ model [29] and into Higgs boson or higgsinos in the DFSZ model [30]. (Note that it does not perturbatively decay into the axion pair since we are assuming ∆m 2 = 0.) We parametrize the perturbative saxion decay rate as #7 #5 This is the initial condition for the zero-point fluctuations in the Minkowski vacuum. Note that the long wave de-Sitter fluctuation generated during inflation may affect this initial condition depending on the inflation scale. However, the following result that the most enhanced mode is k ∼ m φ remains intact and the phenomenological consequence also remains the same.  10.75 where g * denotes the relativistic degrees of freedom at the saxion decay. This is a lower-bound on the axion dark radiation abundance in the sense that this amount is necessarily produced by the saxion oscillation even if ∆m 2 = 0 so that the perturbative saxion decay into the axion pair is forbidden. Although ∆N eff is typically much smaller than the upper bound from the Planck observation [31], still it can have phenomenological impacts [32][33][34][35][36]. Gravitational waves produced due to the explosive axion production may be another interesting observable as discussed in Sec. 4.

Axino
Next we study the production of the axino, or the fermionic superpartner of the axion. There are three chiral fermions in the model (2). Denoting the fermionic components of φ,φ and X by ψ,ψ and χ, χ and one linear combination of ψ andψ obtain the Dirac mass of λ |φ| 2 + |φ| 2 . On the other hand, there is a light mode, which we call the axino, defined by Its mass is given by Thus during the saxion oscillation the axino mass also oscillates and the axino is produced. Below we evaluate the axino production rate for the small and large saxion amplitude regime, respectively.

Small saxion amplitude regime
First let us consider the axino production in the small saxion amplitude case: ϕ amp f . The axino mass is expanded as Thus the perturbative decay rate of the saxion into the axino pair is given by Note that it also vanishes in the case of ∆m 2 = 0. Also such a process can be kinematically forbidden if m φ < 2m A even if ∆m 2 = 0. #8 Therefore it is possible to avoid the nonthermal axino production from the saxion dynamics as long as its amplitude is small enough. However, as we shall see below, the axino production is unavoidable in the opposite regime ϕ amp f . Below we assume ∆m 2 = 0 to derive the lower bound on the axino abundance.

Large saxion amplitude regime
Now we consider the large saxion regime: ϕ amp f . In this case, the axino mass (37) shows a sharp behavior due to the saxon dynamics. Approximately we have Fig. 2 plots |m A | during one-half saxion oscillation for ϕ amp = 100f . It is not hard to imagine that this peculiar behavior of the axino mass results in the axino production. We numerically evaluated the time evolution of the phase space density of the axino. The linearized equation of motion of the axino in the momentum space is given by (see Refs. [37][38][39][40] for more details on the fermion production) with the initial condition A k = (ω k + m A )/ω k andȦ k = −iω k A k , corresponding to the zero-point fluctuation in the Minkowski vacuum. The phase space density is given by The right panel of Fig. 2 shows the phase space density of the axino f (A) k as a function of wave number k in units of m φ for 1, 3 and 5 half saxion oscillation (t = π/m φ , 3π/m φ , 5π/m φ ) for ϕ amp = 100f . We have taken m 3/2 = m φ . We find that one-half saxion oscillation yields f The production is not as violent as the axion case studied in the previous subsection, but still there can be a significant axino production. #8 Here we omitted the δϕ 2 term in the expansion of the axino mass. Such a term induces the saxion "annihilation" into the axino pair even for ∆m 2 = 0, although it is also kinematically blocked for m φ < m A . The dominant contribution to the axino abundance may come from those created at ϕ amp ∼ f , after which the axino production will be suppressed as shown in Sec. 3.2.1. The typical momentum of the axino is ∼ m φ and the phase space density around k ∼ m φ is likely to be saturated as f This may not be negligible depending on the mass of the LSP, although it is smaller than the contribution from the thermal axino production [11,12,16,17].

Summary and discussion
We revisited the saxion dynamics in a SUSY axion model (2) with two PQ fields, and studied the axion production due to the saxion oscillation. We found a strong spike-like behavior in the canonical axion field during the course of saxion oscillation, which induces explosive production of axions if the saxion oscillation amplitude is larger than the PQ scale. The production is so efficient that the axion energy density becomes comparable to the initial saxion energy density. The typical momentum of the saxion and axion is of order the saxion mass scale m φ , and the initial saxion homogeneous condensate is expected to be transferred to the modes with k ∼ m φ . At this stage, both the saxion and axion have effective masses of order m φ hence they are semi-relativistic. Then the saxion and axion will keep to have a comparable energy density until the effective saxion amplitude becomes smaller than the PQ scale due to the Hubble expansion. After that the axion becomes relativistic while the saxion is non-relativistic and hence the saxion begins to dominate over the axion. Eventually the saxion decays in a perturbative way. This final saxion decay process is analyzed perturbatively around the vacuum and it is possible to suppress the decay into the axion pair by choosing appropriate parameters. But the non-perturbative axion production is unavoidable when the saxion oscillation amplitude is large, which translates into a lower bound on the axion dark radiation abundance (35). Some comments are in order. We focused only on the axion and axino production since they are almost model independent. Actually, however, the saxion has model-dependent couplings to the visible sector. For example, in the KSVZ model the saxion couples to vector-like (s)quarks and in the DFSZ model it couples to the Higgs bosons and higgsinos. These coupled particles are also efficiently produced due to the saxion oscillation through the preheating process. In an extreme case, it may be possible that these processes make high temperature thermal bath and the saxion and axion are thermalized much before the saxion decays perturbatively, although the complete analysis is highly model-dependent. In any case, the most important feature that the saxion-axion system becomes inhomogeneous shortly after the onset of saxion oscillation remains intact.
One of the possible phenomenological consequences is the production of the gravitational waves (GWs) due to the violent axion production and the non-linear dynamics of saxionaxion system. Assuming that the saxion and axion are non-linearly inhomogeneous and their gradient energy density with a physical wave number k ∼ m φ are comparable to the total energy density of the universe, the energy of GWs emitted from the region k −1 in one Hubble time is evaluated as ∆E GW ∼ H −1 M −2 P (d 3 Q/dt 3 ) 2 where Q ∼ ϕ 2 /k 3 is the quadrupole moment of such a region. Noting that there is (m φ /H) 3 such regions in the Hubble volume, the relative GW energy density to the total energy density is given by 4 . It is O(1) just after the onset of saxion oscillation for ϕ ini ∼ M P . Thus the GW production is efficient at the early stage after the oscillation, but it will be eventually diluted by the saxion-induced entropy production after the resonant axion production stops. Quite roughly, we estimate the present abundance of the stochastic GW background as where Ω rad 8.5 × 10 −5 denotes the present radiation density parameter and β 0.3 accounts for the deviation from ρ rad ∝ a −4 scaling law due to the change of relativistic degrees of freedom. Note that it may overestimate the GW abundance by an order of magnitude. The peak GW frequency at present in the present universe is ∼ 50 Hz C 1/6 228.75 g * (H = Γ φ ) 1/3 m φ 10 6 GeV 5/6 It may lie in the observable range in the future GW detectors such as DECIGO [41] for some reasonable parameter choices. Note that typical physical wavenumber of GWs just after the emission is always ∼ m φ , hence GWs emitted later experience less redshift and constitute higher frequency modes in the present epoch. The precise estimation of the GW spectrum requires numerical simulation, which is beyond the scope of this paper. One may also wonder whether topological defects are formed or not due to the efficient axion production. This issue is discussed in App. B where it is shown that we have no topological defects in this class of models. Thus the energy density of produced χ particles is dominated by those with momentum k ∼ m sp and it is given by In the case of axion production studied in the main text, we have It is consistent with numerical calculation shown in Fig. 1. Thus particles with f k 1 (k < m sp ) are produced after one saxion oscillation. After that, the production is even more enhanced through the parametric resonance or the Bose enhancement effects.

B Comment on domain wall formation
Here we comment on the possibility of domain wall (DW) formation due to the saxion oscillation. It is known that if the PQ field oscillates with large amplitude, the parametric resonant enhancement of the axion fluctuation can lead to the nonthermal symmetry restoration [8,[42][43][44]. In such a case, the formation of topological defects, i.e., cosmic strings after the PQ phase transition and DWs after the QCD phase transition, would be unavoidable.
However, Ref. [8] also claimed that DWs may be formed even in the case without nonthermal symmetry restoration, because the saxion oscillation in a SUSY axion model leads to the resonant enhancement of the axion fluctuation and the phase of the PQ field might be uniformly distributed in the universe. Below we show that the DW formation does not likely to happen in such a setup. Let us consider a toy model Writing φ = (f + σ + ia)/ √ 2, we set the initial condition σ = σ i and a = 0. If 0 < σ i < ( √ 2 − 1)f , the initial potential energy is lower than the potential energy at the origin φ = 0 and hence nonthermal symmetry restoration cannot occur with this initial condition. #9 Still, the axion fluctuation as well as the fluctuation of the radial mode is amplified during the σ oscillation around the vacuum σ = 0. The question is whether such an enhancement of the axion fluctuation leads to the formation of DWs or not. A naive estimation is that the parametric resonance stops when the axion fluctuation δa 2 ≡ a 2 − a 2 becomes roughly equal to f 2 , around which the backreaction becomes important. Thus we need to simulate it numerically to correctly estimate the field variance.
We performed classical lattice simulation using LATTICEEASY [46]. The initial condition is set to be σ i = 0.3f and we have taken λ = 10 −5 and f = 0.1M P with the grid size of 128 3 . The results are shown in Fig. 3. The evolution of the field variance δa 2 and δσ 2 ≡ σ 2 − σ 2 (in units of f 2 ) are shown as a function of time m σ t where m σ = √ 2λf . The Hubble expansion is not included in the left figure while a power-law Hubble expansion R(t) ∝ t 1/2 is assumed in the right figure. It is seen that, although the fluctuation is resonantly enhanced, the resonance stops before reaching δa 2 f 2 . It means that the axion field is not uniformly distributed from θ = 0 to ±π where θ ≡ arg(φ) even after the parametric resonance. Moreover, the fluctuation decreases as the universe expands. #10 Since the QCD phase transition happens far after the σ field oscillation, we expect θ is almost zero in the whole universe around the QCD phase transition epoch, and hence there is no DW formation. The story is the same for the SUSY axion model studied in the main text.