The QCD Axion: A Unique Player in the Axiverse with Mixings

In an axiverse with numerous axions, the cosmological moduli problem poses a significant challenge because the abundance of axions can easily exceed that of dark matter. The well-established stochastic axion scenario offers a simple solution, relying on relatively low-scale inflation. However, axions are typically subject to mixing due to mass and kinetic terms, which can influence the solution using stochastic dynamics. Focusing on the fact that the QCD axion has a temperature-dependent mass, unlike other axions, we investigate the dynamics of the QCD axion and another axion with mixing. We find that the QCD axion abundance is significantly enhanced and becomes larger than that of the other axion for a certain range of parameters. This enhancement widens the parameter regions accounting for dark matter. In addition, we also find a parameter region in which both axions have enhanced abundances of the same order, which result in multi-component dark matter.

In the axiverse scenarios, it is known that the abundance of axions produced by the misalignment mechanism [41][42][43] will be too large unless the initial angles are fine-tuned. This is nothing but the cosmological moduli problem [44,45]. In particular, the larger the mass of the axion, the greater the abundance. For example, the QCD axion is known to have excessive abundance when the decay constant is larger than 10 12 GeV, which sets the upper end of the so-called axion window.
The cosmological moduli problem can be mitigated if the Hubble parameter during inflation is sufficiently small and the inflation lasts long enough. For example, eternal inflation with a low energy scale fits this scenario [46]. This is because the initial misalignment angle θ i becomes much smaller than O(1) if the inflation lasts long enough, allowing the axion field to follow the so-called Bunch-Davies distribution, which balances quantum fluctuations and classical motion. It has been found that the QCD axion with a decay constant f a ≳ 10 12 GeV, which exceeds the upper limit of the axion window, is allowed without overclosing the universe, as long as the Hubble parameter is below the QCD scale during inflation [47,48]. In Ref. [49], this stochastic approach was applied for the first time to the cosmological moduli problem posed by numerous axions appearing in the axiverse, not just the QCD axion. In this case, the axion abundance is greater for the lighter axion, because although the energy density at the onset of oscillations is the same, the lighter axion starts oscillating later. Their results show that there is an upper bound on the Hubble parameter during inflation, H inf < keV -MeV, depending on the typical decay constant of axions, and that the cosmological moduli problem in the axiverse is solved when this bound is satisfied. This was also subsequently confirmed in Ref. [50].
Note that the stochastic dynamics can be altered, for example, when a Hubble-induced mass is present [51], when the potential deviates from the quadratic potential [52,53], or when the QCD gauge coupling is strong during inflation [54]. In such cases, the axion abundances are known to change. In particular, an important implicit assumption in stochastic axion scenarios is that the axion minima do not change during and after inflation [48]. This assumption is violated, e.g., if the inflaton is an axion that mixes with other light axions. In such cases, the probability distribution of axions could shift from the potential minimum to near the maximum, leading to a significant increase in the abundance of axions. Consequently, it becomes possible to explain all dark matter with axions with the decay constants as small as the astrophysical lower bound [53,55,56].
In addition, most analyses of axions in the axiverse so far have not considered mixing between axions, and cosmological and astrophysical effects have been studied for individual axions because the axion masses are hierarchical. On the other hand, it has been pointed out in Refs. [57][58][59][60][61][62][63][64] that the mixing of many axions through mass and kinetic terms can have significant cosmological consequences, and a typical example is the so-called clockwork/alignment mechanism. Another interesting phenomenon that is characteristic when the QCD axion is composed of multiple axions is resonance phenomena similar to the MSW effect in neutrino oscillations [65][66][67]. Through the resonance, the QCD axion can be converted to lighter axion-like particles and vice versa. In some cases, the axion starts to run along a lighter flat direction, going over potential hills and troughs [68]. When multiple axions are present and mixed with each other in this way, the dynamics can lead to complex and interesting phenomena.
In this paper, we study for the first time the mixing effect between string axions and QCD axion under stochastic dynamics. Normally, in stochastic axion scenarios, the initial field values are determined by the equilibrium distributions during inflation. However, in the presence of mixing effects and temperature dependence of the mass, we find that the initial values set during inflation can be significantly modified by the post-inflationary axion dynamics. Specifically, when the axion potential is generated by QCD instanton effects and another non-perturbative effect, the axions mix through the mass term, causing the mass eigenstates to vary in time due to the temperature dependence of the QCD potential. We find that if each of the two potentials has a mass of the same order at the onset of field oscillations, the total energy density of the axions can be significantly enhanced compared to the case without mixing. In particular, while the lighter axion tends to have a larger abundance in the stochastic scenario, the QCD axion, which is heavier than the mixing partner in the vacuum, can dominate the abundance due to this enhancement. This enhancement breaks the one-to-one correspondence between the axion mass and the abundance, thus broadening the parameter range that explains the dark matter. This effect is analogous to the aforementioned shift of the potential minimum due to mixing with the inflaton, but in our scenario, it involves only simple dynamics of axions. Also, unlike resonance phenomena, there need be no adiabatic invariants and therefore no large mass hierarchy. This is therefore an example of how mixing between axions and temperature dependence can be very important in axiverse scenarios, especially with stochastic axions.
The rest of this paper is organized as follows. In Sec. 2, we show the model and explain the dynamics of the axions during and after inflation. In Sec. 3, we show the results of the numerical simulations of the axion dynamics and demonstrate the enhancement of the axion energy density. Finally, we summarize and discuss the results in Sec. 4.

Stochastic axions with mixings 2.1 Set-up
In the string axiverse, there exist numerous axions at low energies, which acquire potentials from non-perturbative effects such as strong dynamics in hidden gauge sectors. Additionally, compactification of extra dimensions typically induces small instanton effects that generate potentials for axions. As a result of these potential terms, axions generically get mixed with each other. To solve the strong CP problem using these axions, we need axions coupled to gluons, and at least one of them must be extremely light when we switch off the nonperturbative QCD effects. The requirement for such light axions coupled to gluons is nothing more than the quality problem of the Peccei-Quinn symmetry. This problem can be solved naturally in the axiverse [5], where there are many axions whose mass is very light and spans a very wide range.
Although many axions may exist over a wide range of scales, for our interest in the mixing between the QCD axion and other axions, it is sufficient to consider two axions whose masses are not too far apart. We introduce two axions, a and ϕ, and identify their linear combination as the QCD axion, as described below. We consider a low-energy effective Lagrangian for a and ϕ given by with the potential of 1 Here χ(T ) denotes the topological susceptibility of QCD, where we adopt n = −8.16 [71], χ 0 = (75.6 MeV) 4 , and T QCD = 153 MeV. We have neglected the higher-order QCD contributions from the a-meson mixings in the potential since the axion field evolution will turn out to be around the vicinity of the CP-conserving minimum, where the higher-order terms are irrelevant. Note that we have used a field redefinition without loss of generality such that a is the combination that couples to gluons. Then, by taking m ϕ → 0, a becomes the QCD axion, while it is a component of the QCD axion with m ϕ ̸ = 0. However, since we are mainly interested in the case where m ϕ is smaller than the QCD axion mass in the vacuum, then a is the main component of the QCD axion. Therefore, we often refer to a as the QCD axion below. Also, the constant phase in each potential is absorbed into a and ϕ. Thus, the minimum of V QCD , a = 0, is the strong CP conserving point. Throughout the paper, we concentrate on the possibility since we consider that they are both the string axions. We define the temperature-dependent axion mass m a (T ) from χ(T ) = m 2 a (T )f 2 a , which gives the zero-temperature mass Here and hereafter, we denote quantities at the present time by subscript 0. For later use, we define Φ and A as Then, V aϕ (a, ϕ) is a function of Φ only, and it is flat in the direction of A.
Around the origin, a = ϕ = 0, the potential V (a, ϕ) can be approximated by the quadratic terms as The mass matrix, M (T ), is diagonalized by an orthogonal matrix U as where α is the mixing angle, and we assume m H0 > m L0 > 0 without loss of generality. Note that both α and U continuously depend on T . The heavy and light mass eigenstates, s H and s L , are related to a and ϕ as Note that U , m H,L , and s H,L are all temperature-dependent quantities. We show the dependences of m H and m L on m a /m ϕ for f a = f ϕ and |N | = 1 in Fig. 1. For m a ≫ m ϕ , the mass eigenstates are largely determined by V QCD as s H ≃ a and s L ≃ ϕ. Then, the mass eigenvalues become m H ≃ m a and m L ≃ m ϕ . On the other hand, for m a ≪ m ϕ , V QCD is negligible, and Φ and A correspond to s H and s L , respectively. We also see that m H > m L holds for all m a /m ϕ . Thus, we refer to s H and s L as the heavier and lighter modes regardless of temperature, respectively. Around m a = m ϕ , the heavier mode transitions from Φ to a as m a /m ϕ increases. If the adiabatic condition is satisfied during the transition, the resonant conversion between the two axions can take place. In the following, however, we do not need the adiabatic condition, and we will see that the axion dynamics is more complicated. If a has couplings to the SM particles other than gluons, such as photons, then s H and s L will also be coupled to them through the mixing. From the relation we define the effective decay constant related to a as  Using these quantities, we can interpret that s H and s L are coupled to the SM particles with the effective decay constant f eff,H and f eff,L , respectively. We show the dependence of f eff,H on m ϕ for |N | = 1, m a = m a0 , and f a = f ϕ in Fig. 2

Stochastic initial conditions set during inflation
Here, we discuss the mass eigenstates and their typical field values during inflation. During inflation, we assume that the Gibbons-Hawking temperature [72], T GH ≡ H inf /2π, is much lower than T QCD . This assumption implies an upper bound on H inf as In this case, the QCD axion acquires a potential during inflation, and the topological susceptibility is given by χ = χ 0 . Then, the axion fields are around the potential minimum at a = ϕ = 0. Thus, the mass eigenstates are s H0 and s L0 . If the axion masses, m H0 and m L0 , are much smaller than H inf , the axion fields diffuse around the origin due to quantum fluctuations. If the duration of inflation is sufficiently long, the axion field values follow the Bunch-Davies distribution. In the Bunch-Davies distribution, the variances of the mass eigenstates are given by [47,48] The field values at the end of inflation play a role of the initial condition for the field dynamics after inflation. In the following, we parameterize s H0 and s L0 at the end of inflation by where c H and c L are typically of O(1). Then, the initial values of a and ϕ are given by Let us consider two limiting cases: m a0 ≫ m ϕ and m a0 ≪ m ϕ with no hierarchy between f a and N f ϕ . First, in the limit of m a0 ≫ m ϕ , the potential is dominated by V QCD , and the mass eigenvalues become during inflation. The matrix U 0 is given by where one can see that the mixing angle, α, is very small with f a ∼ f ϕ . In this case, the initial condition is approximated by Note that we have |a init | ≪ |ϕ init | for c H , c L = O(1) in this case. Next, we consider the limit of m a0 ≪ m ϕ . In this limit, the potential is dominated by V aϕ , and the mass eigenvalues become The matrix U is given by Thus, the mass eigenstates are given by and the initial conditions become In this case, Φ approximately remains the mass eigenstate during and after inflation, and the dynamics of the two mass eigenstates always decouple. So in the following we will focus on the case of m a0 ≳ m ϕ to see how the two fields evolve via the mixing. Moreover, if m ϕ ≪ H(T QCD ), with H(T ) being the Hubble parameter at the cosmic temperature T in the radiation dominated Universe, the QCD axion begins to oscillate first, and the other axion begins to oscillate after m a (T ) becomes constant. In this case, the field dynamics are again independent for each of the two axions. Thus, we focus on m ϕ in the mass range of in the following.
Here, we make two comments on the assumptions of our analysis. First, to obtain the Bunch-Davies distributions for the axions, the duration of inflation should be sufficiently long as with N inf being the e-folding number of the inflationary period and N i rela = H 2 inf /m 2 i being the e-folding number required for relaxation of the i-th axion with the mass eigenvalue m i . With multiple axions, the lightest one requires the longest relaxation time, and thus we evaluate the e-folding number for the lightest mass eigenvalue in our setup. For example, this condition becomes N inf > 10 32 for H inf = 5 MeV and m L0 = 5 × 10 −10 eV. Such long inflation does not necessarily require eternal inflation 2 since the upper bound to have non- [47,75]. Second, the fluctuations of the axion fields can spatially modulate the Hubble parameter during inflation. This can be checked by using the aforementioned relaxation time scale, N i rela . Without any mixing effect and assuming the quadratic potential, we can consider the constraint for each axion potential, which contributes to the expansion via 6. This is consistent with the one from Ref. [47].
In the presence of mixings or generic potential shapes, we need to solve the Fokker-Planck equation [76,77] with the volume effect [78][79][80][81]. Here, we analytically derive a conservative bound. With multiple axions, the inflationary Hubble parameter is altered at most by where Λ 4 tot is the total potential height for the multiple axions. If ∆H inf /H inf ×max i N i rela ≪ 1 is satisfied, we need not worry about the back reaction. For our system, we obtain Λ 4 tot = 2(χ 0 + m 2 ϕ f 2 ϕ ) ≃ 2χ 0 , and ∆H inf ≈ 7 × 10 −38 H inf . Thus, ∆H inf /H inf × max i N i rela ∼ 10 −5 and we certainly have a parameter region that we can safely neglect the backreaction. In particular, we can neglect the backreaction effect in the whole range of parameters shown in Fig. 9.

Post-inflationary dynamics before QCD phase transition
Next, we consider the field dynamics after inflation. After the end of the inflation, the inflaton decays into the SM particles, which form a hot thermal plasma. As the temperature of the universe increases and becomes higher than the QCD scale, the axion potential from non-perturbative QCD effects, V QCD , disappears, and we have only the potential as which is only a function of Φ. Depending on the relative size of m Φ (or m ϕ ) and m a , the evolution of the system is different. When m Φ ≃ 3H(T ) > 3H(T QCD ), Φ starts to evolve while A remains constant. After that, the field dynamics depends on when V QCD becomes relevant.
For example, let us consider the limit of m Φ ≫ H(T QCD ) and m a0 ≫ m ϕ . In this limit, Φ damps due to oscillations well before T = T QCD . Thus, the axions settle down at the potential minimum of V aϕ , (a m , ϕ m ), determined by − As a result, we obtain where we have used Eqs. (18) and (19) in the second equalities. Note that, since we are assuming m a0 ≫ m ϕ , we have a init ≪ ϕ init . Thus, for f a ∼ |N |f ϕ , 3 the field values of the axions are modified by the post-inflationary dynamics as We note that a m is enhanced compared to a init , which is crucial for the evaluation of the axion abundances, as we will see in the next section. Since a m /a init is proportional to m −1 ϕ , the enhancement will be more significant for smaller m ϕ as long as Φ starts to oscillate before V QCD becomes relevant. Thus, we expect that the enhancement is most significant for m ϕ such that V QCD and V aϕ becomes relevant at about the same time. In other words, we expect the enhancement for m ϕ ∼ m enh ≡ 3H(T a,osc ) with T a,osc satisfying m a (T a,osc ) = 3H(T a,osc ). Considering H(T ) ∝ T 2 and m a (T ) ∝ T −4.08 , we obtain a m /a init ∝ f −2.04/3.04 a for m ϕ ∼ m enh , which leads to the maximal enhancement of the QCD axion abundance with a factor proportional to f . Note that here we neglected the temperature dependence of the effective degrees of freedom of radiation, g * , in the Hubble parameter.

Enhancement of the QCD axion abundance
In the previous section, we have seen that the amplitude of the QCD axion becomes larger than the initial value due to post-inflationary dynamics caused by mixing. Numerical calculations are needed to determine when and to what extent the QCD axion abundance indeed increases.

Setup for numerical calculations
Here, we perform the numerical calculation of the axion dynamics and show how the enhancement of the axion abundance depends on the model parameters.
The equations of motion for a and ϕ are given bÿ where the dots represent derivatives with respect to the physical time t, and the Hubble parameter, H, is given by the cosmic temperature T through the Friedmann equation in the radiation-dominated era: The time evolution of the temperature is determined by the conservation of the entropy in the physical volume ∝ R 3 with the scale factor R: which leads toṪ = − π 2 10 g * s (T ) We use the temperature dependence of the effective degrees of freedom of radiation for energy density and entropy density, g * (T ) and g * s (T ), given in Ref. [71]. We set the initial conditions with the initial temperature T init satisfying which corresponds to a time well before the onset of oscillations of ϕ. The final time of the simulations is set to be well after the energy densities of s H and s L come to follow ∝ R −3 . The energy density can be expressed in terms of the current density parameter as where ρ mix is the sum of the energy density of a and ϕ, and ρ c is the critical density. Since the energy density of oscillating scalars scales proportionally to the entropy density, we evaluate ρ mix /s at the end of the numerical calculations and obtain where ρ c /s 0 ≃ 3.6 × 10 −9 h 2 GeV with the reduced Hubble constant h ≃ 0.67. The input parameters of the numerical calculations are For simplicity, we fix Since we are interested in the low-scale inflation and the Bunch-Davies distribution whose width is much smaller than the decay constant, the potential can be approximated by the quadratic terms. In this case, |N | is degenerate with f ϕ , and the sign of N can also be absorbed into the definition of ϕ. The dependence of the axion abundances on the initial conditions, c H and c L , will be discussed later. Now, the remaining parameters are f , m ϕ , and H inf . As long as the quadratic approximation is valid, the field values always scale as ∝ H 2 inf , and thus the choice of H inf does not affect the axion dynamics qualitatively. To see the non-trivial dynamics of axions due to mixing effects, we focus mainly on 3H(T QCD ) ≲ m ϕ ≲ m a0 as mentioned above. From the Friedmann equation, we obtain Thus, we will investigate the mass range of 10 −11 eV ≲ m ϕ ≲ m a0 .

Axion dynamics
In the following, we show the numerical results for which corresponds to m a0 = 5.7 × 10 −9 eV .
We consider the following three values of m ϕ , m ϕ = 10 −10.5 eV , 10 −9.5 eV , and 10 −8 eV , as examples for the dynamics with enhancement (m ϕ = 10 −9.5 eV) and without enhancement (m ϕ = 10 −10.5 eV and 10 −8 eV). First, we show the result for m ϕ = 10 −10.5 eV in Fig. 3. The top panel shows the trajectory of a/f a and ϕ/f ϕ . Initially, a/f a is smaller than ϕ/f ϕ because of m a0 ≫ m ϕ . At the very beginning, the axion fields slowly roll down the potential in the Φ-direction. In the figure, it first moves to the right. Then, V QCD grows and the axion field starts to oscillate rapidly in the a-direction. After that, the fields also start to oscillate in the ϕ-direction. The bottom panel shows the time evolution of the energy density. ρ H and ρ L are the energy density of the heavier and lighter modes and ρ mix is their sum. As a comparison, we also consider the  (Bottom panel) Energy densities of the two fields, heavier mode, and lighter mode with mixing and a and ϕ without mixing. ρ mix almost overlaps ρ ϕ,no-mixing . The vertical line represents m ϕ t = 1, which almost corresponds to T = T QCD . As is usually the case with stochastic axions, the lighter axion has a larger abundance. case with N = 0, where a and ϕ decouple from each other. By solving the dynamics of each field with the initial Bunch-Davies distribution, we obtain the energy densities of a and ϕ, ρ a,no-mixing and ρ ϕ,no-mixing . After V QCD arises, the heavier mode is approximately a. Since a grows due to the slow roll in the Φ-direction before oscillations, ρ H is enhanced compared with ρ a,no-mixing . On the other hand, ρ L is almost the same as ρ ϕ,no-mixing since the slow roll of Φ or the oscillation of a has little effect on the time evolution of ϕ. As a result, the total energy density is hardly enhanced.
Next, we show the result for m ϕ = 10 −9.5 eV in Fig. 4. As before, a/f a is initially smaller than ϕ/f ϕ because of m a0 > m ϕ . For T ≫ T QCD , the potential is dominated by V aϕ , and Φ starts to roll down the potential while A remains constant. As the temperature decreases, V QCD becomes relevant and then dominant. Thus, the field motion changes its direction, and a starts to oscillate rapidly. In this process, a acquires a larger field value than the initial condition, and the energy density of the two fields is enhanced compared with the case where the two fields evolve independently. In particular, ρ H is significantly enhanced compared with ρ a,no-mixing while ρ L is not so different from ρ ϕ,no-mixing . As a result, the total energy density is also enhanced due to the interplay of the two fields.
Finally, we show the result for m ϕ = 10 −8 eV in Fig. 5. In this case, the heavier mode (≃ Φ) starts to oscillate and damps well before the emergence of V QCD . Then, ρ L becomes dominant later. As a result, the total energy density is different from the case with N = 0 only by a factor of ≃ 1.6. This is because the lighter mode is approximately equal to A, which has an effective decay constant f eff,A = √ 2f a (see Eq. (11)). Considering the result of the standard misalignment mechanism for the QCD axion, Ω a ∝ f 1.17 a [82][83][84], we expect ρ L /ρ a,no-mixing ≃ 2 1.17/2 ≃ 1.5.

Enhancement in QCD axion abundance
Next, we look at the m ϕ dependence of the enhancement factor. We show the ratio of the abundances of the heavier mode with mixing, Ω H , and the QCD axion without mixing, Ω a,no-mixing , for f = 10 14 , 10 15 , and 10 16 GeV in Fig. 6. Note that the heavier mode s H almost corresponds to the QCD axion a at low temperatures for m a0 ≫ m ϕ . We see that Ω H is enhanced around m ϕ = m enh as expected. On the other hand, the ratio becomes less than unity for larger m ϕ . This is because s H starts to oscillate due to V aϕ earlier than a without mixing for m ϕ > m enh . The enhancement is most significant for f = 10 14 GeV, with which the ratio is Ω H /Ω a,no-mixing = O(10 3 ) for m ϕ ≃ m enh . The maximum ratio for f = 10 14 GeV is larger than that for f = 10 16 GeV by a factor of 333 ≃ 100 1.26 , which validates the relation obtained in Sec. 2.3, Ω H /Ω a,no-mixing ∝ f −1.34 as a rough estimate.     Fig. 3 but for m ϕ = 10 −8 eV. The lighter axion (mostly the QCD axion) has a larger abundance. The slight enhancement over the unmixed case is due to the mixing effect on the decay constant, not the dynamics (see Fig. 2). We also show the ratio of the total abundances of the two fields between the cases with N = −1 and N = 0 against m ϕ for f = 10 14 , 10 15 , and 10 16 GeV in Fig. 7. The enhancement is most significant for f = 10 14 GeV, with which the ratio is Ω mix /Ω no-mixing = O(100) at the peak. The wiggle in the right side of the peak corresponds to the oscillation phase of Φ when the QCD potential becomes relevant. For smaller m ϕ , we obtain Ω mix ≃ Ω no-mixing . On the other hand, for larger m ϕ , we obtain Ω mix /Ω no-mixing ≃ 2 1.17/2 ≃ 1.5 as explained above.
We have seen that ρ H becomes dominant when the enhancement is significant in Fig. 4 in contrast to the case without enhancement in Figs. 3 and 5. To visualize this trend, we show the contributions of the heavier and lighter modes to the enhancement in Fig. 8. Here, we choose f = 10 14 GeV, with which the enhancement is most significant in Figs. 6 and 7. We see that the heavier mode dominates the energy density for the mass region where the enhancement is significant. In this mass region, the heavier mode corresponds to the QCD axion a with m a0 ≃ 5.7 × 10 −8 eV for f = 10 14 GeV. For m ϕ ≃ 4.0 × 10 −9 eV, both the heavier and lighter modes have the same order of energy densities larger than Ω no-mixing . For m ϕ ≲ 3H(T QCD ), the lighter mode ≃ ϕ is dominant, and, for m ϕ ≳ m a0 , the lighter mode ≃ A is dominant.

Viable parameter space for dark matter
If the axion potential can be approximated by mass terms, then the squares of the oscillation amplitudes are proportional to H 4 inf at the end of the inflation, and so are the energy densities at any epoch after inflation with the other parameters fixed. Using this approximation, the energy density derived for H inf = 5 MeV can be converted to H inf such that the axion field explains all dark matter. To see the validity of the mass approximation, we define the typical amplitude of the Bunch-Davies distribution in the ϕ-direction: Note that the typical misalignment angle in the a-direction is much smaller for the parameters of interest:θ These typical angles,θ ϕ andθ a , correspond to the initial values of ϕ/f ϕ and a/f a with c H = c L = 1 in the limit of m a0 ≫ m ϕ . We show H inf explaining all dark matter for f = 10 14 GeV and c H = c L = 1 in Fig. 9. The solid lines represent H inf with which the two fields with mixing (thick gray), a without mixing (red), and ϕ without mixing (blue) explain all dark matter, respectively. Forθ ϕ > π (the gray shaded region), the assumption that a and ϕ oscillate around the origin a = ϕ = 0 is invalid and our analysis cannot be applied directly. In this region, however, the typical initial amplitude of ϕ is of order f ϕ , and the axion abundance is larger than the observed dark matter abundance. For π/2 <θ ϕ < π (between the two gray dashed lines), the mass approximation of the potential becomes inaccurate, and the result requires some correction. We can see that the QCD axion can explain all dark matter with H inf much smaller than the case without mixing effects.

Initial condition dependence
So far, we have taken c H = c L = 1. Here, we discuss the importance of the initial conditions focusing on two parameter sets. First, we consider the peak of enhancement: f = 10 14 GeV , m ϕ = 7.9 × 10 −10 eV , 10 -10 10 -9 10 -8 10 -7 5 10 50 100 Figure 9: Hubble parameter during inflation with which the two fields with mixing (thick gray), a without mixing (red), and ϕ without mixing (blue) explain all dark matter for f = 10 14 GeV and c H = c L = 1. The gray dashed lines representθ ϕ = π/2 and π. and investigate how the enhancement depends on the initial condition. As long as the mass approximation of the potential is valid, the field dynamics is linear and the energy density is proportional to the square of the field amplitudes. Thus, the enhancement factor Ω mix /Ω no-mixing depends on the initial condition through c H /c L . We parameterize the initial condition by 0 ≤ β < π as c H = cos β , c L = sin β , and perform the numerical simulations for both N = −1 and N = 0. Note that V (a, ϕ) is invariant with (a, ϕ) → (−a, −ϕ) and that π ≤ β < 2π leads to the same energy densities as 0 ≤ β < π. We show the dependence of the enhancement factor on β in Fig. 10. The enhancement is suppressed around β = 0 and π, where c L is small compared with |c H |. This behavior can be understood by the observation that the enhancement is due to the conversion of ϕ to a in the early stage of field oscillations (see Fig. 4). For typical initial conditions, ϕ init /f ϕ is larger than a init /f a , and the motion in the Φ-direction enhances a/f a . On the other hand, for c L ≪ |c H |, ϕ init /f ϕ ≲ a init /f a and the enhancement does not occur. Next, we consider m ϕ with which Ω H ≃ Ω L in Fig. 8:  We show the energy ratio of the heavier mode, Ω H /Ω mix , in Fig. 11. We see that both the heavier and lighter modes have non-negligible energy density except for two regions near β = 0 and π. These two regions lead to Φ init ≃ 0 and A init ≃ 0, respectively. Then, one of ρ H and ρ L is highly suppressed as the initial condition because (Φ, A) correspond to (s H , s L ) at high temperatures T ≫ T QCD . Since the fields start oscillations before the emergence of V QCD for the parameters in Eq. (54), the energy densities of the heavier and lighter modes are not transferred to each other due to the resonant conversion [65][66][67]. As a result, one of the heavier and lighter modes becomes dominant for the two regions near β = 0 and π.

Conclusions and discussions
The stochastic axion scenario fits well with the string axiverse, where there are many axions with masses spread over a wide parameter range. As long as the inflationary scale is kept relatively small, all axions will stay near the potential minimum, thus avoiding the notorious cosmological moduli problem in the string axiverse. In this paper, we have shown that among the axions in the axiverse, the QCD axion is special in the context of the stochastic axion scenario because it necessarily has a temperature-dependent potential. Even if the axions have suppressed initial misalignment angles in the stochastic scenario, the mixing between the QCD axion and another axion and the time dependence of the QCD axion potential make the field trajectory after inflation quite non-trivial. Especially when the mixing is nonresonant, the two axion exhibits a highly complicated behavior. Through this dynamics, the QCD axion abundance can be enhanced by many orders of magnitude compared to the case where mixing is neglected. Let us see if the QCD axion can be a dominant component of dark matter in the axiverse with stochastic axions. In this scenario, the lighter axion tends to have a greater abundance for the same decay constant. This is because while the energy density at the onset of oscillations is of the order of H 4 inf , the lighter axion starts oscillating later. On the other hand, the initial amplitude cannot exceed the decay constant, and therefore there is a lower bound on both H inf and the axion mass m ϕ to account for all dark matter by the (lightest) axion in this scenario, as shown in Ref. [49]. For example, for f ϕ = 10 16 (10 14 ) GeV, the lower bound is H inf ≳ 10 keV (10 MeV) and m ϕ ≳ 10 −18 (10 −10 ) eV. Now, from Fig. 9, one can see that the QCD axion can explain all the dark matter for H inf ≃ 7 MeV and m ϕ ≃ 10 −9 eV with mixing. Thus, if the decay constant for other axions is universal and equal to 10 14 GeV, the contributions of the other axions are subdominant, i.e., the QCD axion is the dominant component of dark matter in the string axiverse scenario. This decay constant is somewhat lower than those conventionally adopted in the string axiverse, but could potentially be realized through a large volume compactification scenario [39]. We also note that the decay constant can be larger than 10 14 GeV, if there is no axion near the lower bound on the mass.
Inflation with H inf of order MeV is a low-scale inflation, but it is high enough for successful cosmology. This is because the corresponding energy scale of the inflaton potential is about 10 7-8 GeV, and the reheating temperature can be significantly higher than the weak scale. Consequently, we could use the active sphaleron reaction along with several potential baryogenesis mechanisms, such as leptogenesis and electroweak baryogenesis.
Interestingly, the QCD axion dark matter with the decay constant of O(10 14 ) GeV can be searched for through e.g., lumped element experiments [85][86][87][88]. If the non-trivial dynamics was caused by the mixing between the QCD axion and another axion, there should be another axion in the mass range of 10 −11 eV to the mass of the QCD axion. The existence of such an axion with a mass close to that of the QCD axion has been discussed in various contexts [11,[64][65][66][67], and if we can find both of them, such an axiverse scenario with the stochastic axions would be one of the plausible possibilities.