Cosmological Perturbations of Axion with a Dynamical Decay Constant

A QCD axion with a time-dependent decay constant has been known to be able to accommodate high-scale inflation without producing topological defects or too large isocurvature perturbations on CMB scales. We point out that a dynamical decay constant also has the effect of enhancing the small-scale axion isocurvature perturbations. The enhanced axion perturbations can even exceed the periodicity of the axion potential, and thus lead to the formation of axionic domain walls. Unlike the well-studied axionic walls, the walls produced from the enhanced perturbations are not bounded by cosmic strings, and thus would overclose the universe independently of the number of degenerate vacua along the axion potential.


Introduction
The QCD axion is a Nambu-Goldstone boson which arises in association with spontaneous breakdown of Peccei-Quinn (PQ) symmetry, and it dynamically solves the strong CP problem [1,2,3]. The strength of interactions with the standard model particles as well as the periodicity of the axion potential are determined by the axion decay constant, which is currently constrained as 10 9 GeV f a 10 12 GeV. (1.1) Here the lower bound is set by astrophysical arguments of star cooling [4], while the upper bound comes from the requirement that the abundance of the axion be less than or equal to that of cold dark matter (CDM) without fine-tuning the initial misalignment angle [5,6,7]. Cosmological considerations further set bounds on a combination of the decay constant f a and the inflationary Hubble scale H inf . Firstly, f a should be larger than H inf , otherwise the PQ symmetry breaking would happen after inflation and thus would produce domain walls which overclose the universe. An exception to this statement is when the number of degenerate vacua along the bottom of the PQ scalar's Mexican hat potential (the so-called domain wall number) is N = 1; then the walls are bounded by cosmic strings without being connected to other walls, and such walls of finite size can annihilate soon after the QCD phase transition [8,9,10,11]. Furthermore, if the axion constitutes the CDM, then even tighter bounds on f a and H inf are obtained from constraints on CDM isocurvature perturbations [12,13]. For instance, with a decay constant of f a = 10 12 GeV, current constraints on CDM isocurvature from measurements of the cosmic microwave background (CMB) [14] require H inf 10 7 GeV. The bound on H inf is particularly strong for f a 10 11 GeV, due to the nonquadratic form of the axion potential [15], dubbed the anharmonic effect. The anharmonic enhancement of the axion isocurvature perturbations becomes so strong for small values of f a that a decay constant of f a 10 9 GeV is excluded [16]; thus the lower bound of (1.1) can also be derived from cosmology when the dark matter consists of axions. These observations indicate that, any detection of primordial gravitational waves from inflation in the near future would exclude the QCD axion as a dark matter candidate.
Here it should be noted that the above arguments on f a and H inf from domain walls and isocurvature perturbations are modified if the decay constant takes different values during the inflation epoch and today. The works [10,17] pointed out that if f a was larger during inflation, then the axion isocurvature perturbations would actually be smaller than that inferred from the present value of f a , allowing high-scale inflation without contradicting observations. (See also Refs. [18,19,20,21] for recent works along this line. A similar effect can also be induced by a non-minimal coupling of the axion kinetic term to gravity [22].) While this is certainly an intriguing possibility, however, the introduction of a dynamical decay constant not only leads to smaller-than-expected isocurvature perturbations on the CMB scales, but also dynamically enhances the isocurvature perturbations on smaller length scales. This is understood by viewing the decay constant f a and the axion field θ respectively as the radial and angular directions of the PQ scalar; if the radial component f a is forced to quickly shrink while the angular velocityθ is nonzero, then the motion along the angular direction would speed up, leading to the growth of the fluctuation δθ among different regions of space. Moreover, this effect not only works on sub-horizon scales, but further stretches out to scales much larger than the Hubble radius at the time when the decay constant varies. 1 In this paper we examine the evolution of the full axion perturbation spectrum under a timedependent decay constant. We show that, while the dynamical decay constant can help to evade the isocurvature constraints on CMB scales, it can also strongly enhance the axion fluctuations on smaller scales, which may even exceed the periodicity of the axion potential and thus lead to the formation of axionic domain walls. These domain walls are produced due to the enhanced axion fluctuations, although the PQ symmetry is broken before inflation, and thus they are different from the well-studied axionic walls. In particular, the axionic walls from the enhanced fluctuations are not bounded by cosmic strings, and thus would overclose the universe even in the case of N = 1 [24]. The strong enhancement of the axion fluctuations may also lead to the recovery of the PQ symmetry in the post-inflation epoch, which would upset the cosmological expansion history. 2 On the other hand, if the axion perturbations are only mildly enhanced, then such small-scale isocurvature could cause a variety of observable consequences, as discussed in e.g. [30,31].
Before moving on, we also comment on another possibility to suppress the axionic isocurvature 1 A similar effect was invoked in [23] for the generation of large-scale magnetic fields in the post-inflationary universe. 2 In this paper we investigate the enhancement of the sub-and super-horizon scale axion fluctuations due to the shrinking of the decay constant. We also note that if the radial component of the PQ scalar oscillates about its minimum, both the axion and the radial component can be produced by parametric resonance. Such sub-horizon fluctuations can lead to a nonthermal symmetry restoration, which has been studied in [25,26,27,28,29,21].
perturbations. If the PQ symmetry is badly broken during inflation, the axion mass may become heavier than or comparable to the Hubble parameter. Then, it does not acquire sizable super-horizon fluctuations, suppressing the isocurvature perturbations [32,18,33,34,35,36,37]. In this case, one has to make sure that the explicit PQ symmetry breaking is sufficiently suppressed in the present universe to be consistent with the neutron electric dipole moment experiments. The layout of this paper is as follows: We describe the setup and notation in Section 2. In Section 3 we give general discussions on the evolution of the axion fluctuations under a dynamical decay constant, then in Section 4 we compute the fluctuations in the post-inflationary universe. In Section 5, we discuss the cosmological implications of the enhanced fluctuations, and also carry out case studies. We conclude in Section 6.

QCD Axion with a Dynamical Decay Constant
Throughout this paper we denote the PQ symmetry breaking scale by f , i.e., the field value of the PQ scalar at the minimum of its Mexican hat potential being |Φ| min = f / √ 2. Writing the PQ scalar in terms of two real scalar fields r and θ as then the angular direction θ along the potential minimum serves as the axion field. We refer to f and θ as the axion decay constant and the (dimensionless) axion field, respectively. We couple Φ to N heavy PQ quark pairs (Q i ,Q i ) in the fundamental and anti-fundamental representation of SU(3) c , where the flavor index i runs from 1 to N . Then, the axion couples to gluons through the color anomaly of the PQ symmetry. Integrating out the heavy PQ quarks and setting the radial component to the potential minimum, the action of the axion is given by This gives an effective potential for θ with periodicity 2π/N , Here, the temperature-dependent axion mass m QCD (T ) emerges when the cosmic temperature cools down to around Λ QCD ≈ 200 MeV. On the other hand at T Λ QCD , the axion is a free field with only the kinetic term. We also note that the discussions in Section 1 can be rewritten for f by the replacement f a → f /N ; in particular, (1.1) is read as a bound on f /N . Now, as we have discussed in the introduction, the coefficient of the kinetic term can be timedependent when the decay constant f is actually a dynamical field [10,17]. (See also e.g. [21,26,38] for explicit models of axions with dynamical decay constants.) Let us now study the axion dynamics with a dynamical f , at very early times when the axion potential is absent. We consider a flat FRW background ds 2 = a(τ ) 2 (−dτ 2 + dx 2 ), (2.4) and assume the decay constant to be homogeneous, i.e. f = f (τ ). Then the equation of motion of the axion reads (a 2 f 2 θ ) where a prime denotes a derivative in terms of the conformal time τ . Going to the Fourier space, the equation of motion is rewritten as We stress that this is the exact equation of motion of the axion before the potential emerges, and in particular that there is no mode-mode coupling. When quantizing the theory by promoting θ to an operator, one can check that the commutation relations of θ and its conjugate momentum are equivalent to those of the annihilation and creation operators when the mode functions are independent of the direction of k, i.e., with k = |k|, and further obey the normalization condition (2.9) By choosing the vacuum |0 to be annihilated by all the annihilation operators, the two-point correlation function of the axion field is computed as where the power spectrum is (2.11)

Dynamics of Axion Perturbations
Although the axion potential is absent in the very early universe, the axion velocity, i.e. the angular velocity of the PQ scalar, is not exactly zero due to quantum fluctuations, which are stretched to cosmic scales and become classical during inflation. The velocity fluctuations can get strongly enhanced even on super-horizon scales when the decay constant f , corresponding to the radial component of the PQ scalar, quickly shrinks. The enhanced velocity fluctuations would then lead to the growth of the axion field fluctuations. The basic picture can be seen by looking at the homogeneous mode k = 0, for which the equation of motion (2.7) gives the scaling of This clearly shows that the phase velocity, given an initially nonzero value, increases when f decreases quickly enough. However we should also note that for a nonzero k mode, even when it is well outside the horizon, the scaling behavior of the velocity fluctuation is not necessarily the same as for the homogeneous mode (3.1).
Let us now look into the inhomogeneous modes, i.e. k = 0, which represent the axion fluctuations. To make the discussions concrete, we focus on an FRW background with a constant equation of state parameter w that is not −1/3, i.e., Then for the Hubble rate H = a /a 2 , its time derivative obeys The condition of w = −1/3 guarantees that aH does not stay constant, which turns out to be convenient when we rewrite the equation of motion below. We further assume that the decay constant varies in time as a power-law of the scale factor, with a constant n. Then, introducing the equation of motion (2.7) can be rewritten in the form The general solution of this equation is ν ) is a pair of Hankel or Bessel functions of the first and second kinds, and κ, λ are constants. The phase velocity can be expressed in terms of lower order functions as where we have used a relation for the Hankel/Bessel functions, We also note that the variable u varies in time as u ∝ a (1+3w)/2 in terms of the scale factor, and u → 0 (u → ∞) describes the modes being well outside (inside) the Hubble horizon.

Outside the Horizon
Let us study the behavior of the solution (3.7) when the modes are well outside the horizon. Here we discuss the general solution using Bessel functions, i.e., The limiting forms of the Bessel functions for u → 0 are [39], (3.11) Here, since the leading order expressions of J ν−1 and Y ν−1 scale differently in term of u, they do not continually cancel each other in the expression (3.8) for the phase velocity. 3 Thus on super-horizon scales k aH, the κ and λ terms in (3.8), together with the prefactor u ν+1 , scale as ∝ a ν(1+3w) and ∝ a 1+3w , respectively; whichever grows faster eventually dominates the super-horizon phase velocity.
To summarize, unless one of the constants κ, λ are set to zero, the phase velocity on modes well outside the horizon k aH scales eventually as The latter case (3.13) is not captured by the analysis of the homogeneous mode (3.1). One can already see in the equation of motion (3.6) that the super-horizon limit of the inhomogeneous modes may behave differently from an exactly homogeneous mode; there k = 0 corresponds to u = 0, and thus the dynamics can be essentially different between homogeneous and inhomogeneous modes. We also note that, although the parameter ranges for (w, ν) shown in (3.12) and (3.13) do not cover all the possibilities, they will be sufficient for the cases we discuss in the following sections.
Expressing the scaling behaviors of the two cases (3.12) and (3.13) collectively as with p being either ν(1 + 3w) or (1 + 3w), then one sees that the phase velocity corresponds to a dynamical mode for the phase fluctuations which scales as The total super-horizon phase fluctuation contains this dynamical mode as well as the usual constant mode.

Inside the Horizon
The behavior of the modes when deep inside the horizon is seen from the asymptotic forms of the Bessel functions for u → ∞ (recall that u is a positive variable), Therefore, θ k inside the horizon is an oscillating solution with an oscillation amplitudeθ k that scales in time asθ In particular when f is a constant, i.e. n = 0, the fluctuation simply damps as ∝ a −1 due to the expansion of the universe.

Evolution of Axion Perturbations after Inflation
Let us now apply the above discussions to study the evolution of axion perturbations under a dynamical f in the post-inflation epoch. We suppose the PQ symmetry to be broken before inflation, and we particularly study the case where f takes a constant value f inf until some time after inflation, 4 then varies with time as ∝ a −n until it reaches the present-day value f 0 , i.e., Here a i and a f denote the scale factors in the post-inflation epoch when the time-evolution of f starts and terminates, respectively. Quantities measured at a = a i(f) will be represented by the subscripts "i(f)". We further assume that the universe becomes (effectively) matter-dominated, i.e. w = 0, right after inflation, and that f evolves in time during this phase. The dynamics of the axion fluctuations can be analyzed by connecting the solutions for θ k in each epoch. However for clarity, instead of connecting the exact solutions (3.7), (3.8), we will mostly use the asymptotic scaling solutions we have obtained for the super-horizon modes (3.12), (3.13), and the sub-horizon modes (3.17). This will provide useful approximations, as we will see later when we compare with the exact solutions.

Inflation Epoch
The axion fluctuations during inflation with constant Hubble rate H = H inf and decay constant f = f inf is given by (3.7) with w = −1 and ν = 3/2. The coefficients κ and λ are fixed from the normalization condition (2.9) and the requirement that the mode function approaches a positive frequency solution in the asymptotic past, i.e., requiring a Bunch-Davies vacuum. For this purpose it is convenient to use the Hankel function of the first kind, with which the solution is expressed as, up to an unimportant phase factor, In the super-horizon limit k aH inf , these solutions are approximated by (4.5) One clearly sees from (4.4) that θ k is constant at leading order (we denote this value by θ const k hereafter), and that the dynamical mode described by the phase velocity (4.5) is subdominant on super-horizon scales.

Post-Inflation Epoch with Constant f
The (4.9) The phase fluctuation is a sum of the constant mode and the growing mode obtained by integrating the phase velocity, , (4.10) which shows that the growing mode dominates over the constant mode on wave numbers k larger than . (4.11) One sees that k grow initially coincides with the horizon scale, i.e. k grow (τ i ) ∼ a i H i , then k grow becomes smaller than aH as f shrinks. Thus the wave number k grow (τ f ) corresponds to the smallest k-mode (largest length scale) where the axion fluctuation is enhanced due to the shrinking of f . After f becomes constant, i.e. a > a f , the fluctuations outside the horizon becomes constant once again. 5 Let us focus on the wave number that re-enters the horizon when f approaches its present value f 0 , i.e., The fluctuation on k f when a = a f is computed by extrapolating the super-horizon growing mode in (4.10) to the time of horizon entry; in terms of the power spectrum (2.11) we find where we have ignored the factor (2n − 3 2 ) −1 in (4.10). Note that this expression does not depend directly on H f itself; the enhanced fluctuation amplitude is independent of when f varies.
After f approaches f 0 , i.e. a > a f , the wave mode k f is inside the horizon and thus the fluctuation is damped as One can also check that in this epoch, until the mode k grow (τ f ) enters the horizon, the fluctuation spectrum has a plateau in the range aH k k f , as we will see in Figure 1. Hence, denoting the time when k grow (τ f ) enters the horizon as τ = τ grow , the fluctuation amplitudes of the modes coming into the horizon k = aH are the same as (4.14), When naively connecting the super-horizon scalings (3.12) and (3.13) across a = a f , the super-horizon fluctuations might seem to grow even after f has stopped its evolution. However it should be kept in mind that (3.12) and (3.13) only describe the asymptotic behaviors in each epoch, and so does not necessarily provide good approximations right after the transition. To obtain the correct behavior of θ k ∼ const., one needs to take into account both terms in (3.11) and/or the sub-leading terms that are dropped there.

Backreaction
Before closing this section, let us comment on the backreaction from the phase fluctuations on the expanding universe. By computing the energy momentum tensor sourced by the axion kinetic term in the action (2.3), the energy density of the phase fluctuations is obtained as In the above discussions we have treated the decay constant as a classical background that, while decreasing in time, injects energy into the phase fluctuations. Here, since the energy density of the field(s) that sets the time evolution of the decay constant cannot exceed the total energy density of the universe, ρ θ 3M 2 p H 2 should always be satisfied. In other words, the phase fluctuations cannot be enhanced to become ρ θ ∼ 3M 2 p H 2 without significantly backreacting on the field(s) driving the f -dynamics, and also on the background universe. The coupling between f and θ in the axion kinetic term further sources direct backreaction from the enhanced axion fluctuations. We also note that if ρ θ is enhanced to become comparable to ∼ f 4 , then it can backreact on the radial component of the PQ scalar. On the one hand, the backreaction may force the PQ scalar to climb up the Mexican hat potential, while on the other hand the enhanced angular velocity may prevent the PQ scalar from approaching the origin. It would be interesting to study whether the enhanced axion fluctuations assist or prevent the restoration of the PQ symmetry in the post-inflation universe.

Cosmological Implications of Axion Perturbations
Let us now discuss the cosmological implications of the enhanced axion fluctuations. We will show that while a shrinking f may relax constraints from large-scale isocurvature perturbations, the strongly enhanced axion fluctuations on small scales can lead to disastrous formation of axionic domain walls.

Domain Walls
In the previous section we have seen that, after the decay constant approaches its present-day value, the enhanced axion fluctuation P 1/2 θ holds a plateau over the wave modes aH k k f . These enhanced fluctuations that are scale-invariant up to the Hubble radius damps as ∝ a −1 as the universe expands. Eventually the universe undergoes reheating, and the axion's effective potential (2.3) emerges when the universe cools down to T ∼ Λ QCD ; if the fluctuation P 1/2 θ then is still larger than the period of the potential ∆θ = 2π N , (5.1) then the axion settles down to different potential minima in different patches of the universe, leading to the formation of cosmic domain walls.
Here a few things are worth noting. Since the PQ symmetry breaking happens before inflation in our case of interest, unlike the standard axionic domain walls [8,9,10,11], the walls produced here are not connected to cosmic strings. (To be more precise, some walls are bounded by strings produced before inflation, but the number density of the strings is extremely tiny because of the dilution during inflation.) The axion potential (2.3) can be viewed as a one-dimensional potential with an infinite number of degenerate vacua, and we stress that the formation of domain walls is due to the axion fluctuations being spread over multiple vacua, 6 with the number of different vacua being ∼ P 1/2 θ /∆θ. This is why domain walls can form without having wall junctions, even with numerous vacua. It should also be noted that in the one-dimensional axion potential (2.3), each vacuum only has two vacua on its sides. As a consequence, adjacent domains in the universe that are separated by a single wall can only hold vacua that are adjacent also in the axion potential.
The important question for cosmology is whether infinitely large domain walls form. If they do, they would overclose the universe and spoil the cosmological expansion history. When P 1/2 θ /∆θ 1, with the vast number of vacua, naively, it would seem difficult for a single wall to extend throughout the universe. However, this is not necessarily the case due to the specific features mentioned above. In order to see this, let us refer to the unperturbed background field value of the axion before the potential emerges byθ, and the field value that becomes the nearby potential minimum by θ min s . In other words,θ and θ min s satisfy θ min s − ∆θ/2 <θ < θ min s + ∆θ/2. We also label the other vacua along the one-dimensional axion potential as θ min s±1 = θ min s ± ∆θ, θ min s±2 = θ min s ± 2∆θ, · · · , (5.2) and further divide the vacua into two groups, depending on which side of θ min s they are located on, (Here we included θ min s in Θ min− , but one may choose to include it in the other group as well.) Now, let us crudely model the universe as a collection of cells with size of the Hubble radius, each of which holding a vacuum belonging to either group with an independent probability. If P 1/2 θ /∆θ 1 when the potential emerges, then since the axion fluctuation spreads equally on both sides ofθ, the probability for each cell to fall in Θ min− or Θ min+ would be roughly the same, p ∼ 0.5. Here, in order to tell whether the collection of cells with Θ min− or Θ min+ forms an infinitely large connected region, useful insights can be obtained from studies of percolation theory: Since the percolation threshold p c in three dimensions is typically smaller than 0.5 (e.g. p c ∼ 0.3 for a cubic lattice [44]), an infinitely large region with vacua belonging to Θ min− exists, and likewise for Θ min+ . This means that an infinitely large wall between Θ min− and Θ min+ appears in the universe. Because of the specific vacuum distribution discussed above, this indicates that there exists an infinitely large domain wall between the θ min s and θ min s+1 vacua. Such infinitely large domain walls can also form between other pairs of vacua as well. Over time the domain walls can annihilate each other, however, as each patch of the universe randomly holds Θ min− or Θ min+ at T ∼ Λ QCD , we expect there to be always at least one domain wall that extends throughout the observable universe. The domain walls without junctions will likely obey the usual scaling law [45,46,47] soon after formation.
Therefore we conclude that axionic domain walls would form and overclose the universe, if is satisfied on the mode k f , and thus also for the horizon size mode, when the axion potential emerges. A full treatment of domain walls from enhanced axion fluctuations is beyond the scope of this paper, but it would be interesting to carry out detailed analyses (which presumably requires numerical simulations) of the formation and evolution of the wall network.

Isocurvature Perturbations
When the axion constitutes a non-negligible fraction of CDM, the large-scale fluctuations of the axion serve as CDM isocurvature perturbations, which are severely constrained by CMB measurements.
Here we briefly review how such isocurvature constraints on the axion are relaxed with a dynamical decay constant [10,17]. Let us focus on the vicinity of one of the minima θ min of the axion potential (2.3), where the potential is approximately quadratic, Then the energy density (isocurvature) fluctuation of the axion is estimated as whereθ denotes the axion's unperturbed background field value before the axion acquires its potential. 7 Hence the power spectrum of the CDM isocurvature perturbations, which we denote by P S , is expressed in terms of the field fluctuation spectrum (2.11) as Here Ω θ(CDM) is the axion (CDM) abundance today, and upon moving to the far right hand side we have used the relation between the axion abundance and the phase [48] (note that, as we are considering the axion potential to emerge after the decay constant has approached its present-day value, the dynamics of the unperturbed axion is basically unmodified unless axionic walls are formed): A scale-invariant and uncorrelated CDM isocurvature is constrained on CMB scales by Planck [14] as P S (k * ) 0.040 × P ζ (k * ) (95% C.L., TT, TE, EE + lowP) (5.10) where the pivot scale is k * /a 0 = 0.05 Mpc −1 , and the adiabatic power is P ζ (k * ) ≈ 2.2 × 10 −9 . Supposing k * k grow (τ f ) so that the axion fluctuations on the CMB scales are not affected by the evolution of f , then from (4.4), This combined with (5.8) and (5.10) gives an upper bound on the inflation scale, For instance if f 0 /N = 10 12 GeV and Ω θ = Ω CDM (which is realized with N |θ − θ min | ∼ 1, cf. (5.9)), without the evolution of f (i.e. f inf = f 0 ), the CMB bound on isocurvature perturbations requires the inflation scale to be as low as H inf 10 7 GeV. On the other hand when the decay constant is allowed to shrink, then the upper bound on H inf is relaxed by a factor of f inf /f 0 . However with a dynamical f , the axion perturbations can get enhanced on smaller scales, as we have seen in the previous sections. We will study this effect in more detail in the following subsection.
Before moving on, we should also remark that the above estimation of the isocurvature perturbations breaks down for N |θ − θ min | 1, i.e., when the axion is initially sitting away from the potential minimum. In particular for cases with Ω θ = Ω CDM , such anharmonic effects become important for f 0 /N 10 11 GeV; there the expression (5.9) for the axion abundance is modified, and most importantly, the isocurvature perturbations become much larger than in (5.8). Consequently, the upper bound on the inflation scale is much stronger than in (5.12). See [15,16] for detailed discussions on the anharmonic enhancement of the axion isocurvature perturbations.

Enhancement of Small-Scale Perturbations
We now study two example cases and see how the small-scale axion fluctuations are actually enhanced, possibly leading to formation of domain walls. In both examples we assume the time evolution of the decay constant as was discussed in Section 4, and further suppose the axion to account for all the dark matter, Ω θ = Ω CDM . (5.13)

Case 1
We first study the case with f 0 N = 10 11 GeV, H inf = 10 14 GeV, (5.14) where the inflation scale is set to be roughly the current upper limit [14]. Were it not for a dynamical f , with these parameters the PQ symmetry breaking would happen after inflation and then domain walls would overclose the universe unless N = 1. However this can be avoided with a timedependent f ; one sees from (5.12) that even the isocurvature constraint can be satisfied if f shrinks by i.e. f inf /N 10 18 GeV. 8 Let us now look into the small-scale fluctuations. The axion fluctuation on the mode k f when the evolution of f ceases is estimated from (4.13) as , (5.16) in terms of the potential period (5.1). In particular when n > 5/2, the constraint (5.15) gives a lower bound on the fluctuation. For instance, for n = 4, P After τ = τ f , the universe eventually undergoes reheating (as we have been assuming that f evolves during the matter-dominated era), and the axion potential emerges when T ∼ Λ QCD . Since the fluctuations on the mode k f as well as the horizon size mode are damped as ∝ a −1 from the value (5.17), one sees that domain walls would form if the axion potential emerges before the universe expands by a/a f ∼ 10 4 since when f has approached its present-day value.
We also note that when n is as large as n 12, then the energy density of the fluctuations (4.16) becomes at least comparable to the total density of the universe at τ = τ f , which indicates the breakdown of our treatment of the decay constant and the cosmological expansion as homogeneous backgrounds. With such a large n, the backreaction from the axion fluctuations becomes nonnegligible before the time τ = τ f . On the other hand, whether the fluctuation energy density exceeds f 4 depends also on when f varies. For instance if f shrinks quickly soon after inflation, then f 4 can become smaller than the total density of the universe. In such cases one will have to worry about the recovery of the PQ symmetry due to the gravitational background or the enhanced axion fluctuations.
The fluctuation P 1/2 θ /∆θ is shown in Figure 1, for the case of n = 4 and f inf /N = 10 18 GeV. Here we have plotted the exact solution for θ k , obtained by starting from the initial conditions (4.2) and (4.3) set during inflation, and then connecting to the general solutions (3.7) and (3.8) in each epoch of the f evolution (4.1) after inflation. Here it should be noted that the fluctuation amplitude at τ = τ f (4.13) depends on the ratio f inf /f 0 , but not directly on when f evolves; hence for the wave modes shown in the figures, the shape of the fluctuation spectrum is basically independent of parameters such as H f , as long as H inf , f inf /N , f 0 /N , and n are fixed. The left figure displays the fluctuation spectrum as a function of k when a = a i (purple), a f (blue), 10 4 a f (red). One sees that the initially flat super-horizon spectrum is enhanced while f evolves on wave numbers k k grow (τ f ). After f ceases to evolve, the spectrum is flattened as the fluctuations enter the horizon and are damped. Hence, until the mode k grow (τ f ) enters the horizon, the fluctuations coming into the horizon have roughly the same amplitude as that on the mode k f inside the horizon, cf. (4.15).
The right figure shows the fluctuation as a function of a, focusing on the wave number k f . One sees that the fluctuation's dynamical mode dominates over the constant mode some time after f aH| a=10 4 a f k P 1/2 This estimate is seen to agree well with the oscillation amplitude of the exact solution. In order to avoid domain walls from forming in this case, the axion potential should not emerge at least until a ∼ 10 4 a f , before which P 1/2 θ (k = aH) is larger than ∆θ.

Case 2
Now let us consider the parameter set f 0 N = 10 10 GeV, H inf = 10 14 GeV. (5.18) For this f 0 /N , the unperturbed axion field value needs to be initially located close to the potential hilltop in order for the axion to account for all the dark matter, as then the onset of the axion oscillation is delayed and the axion abundance is enhanced compared to that shown in (5.9). However, the anharmonic effects also strongly enhance the axion isocurvature fluctuations, resulting in a much stronger upper bound on the inflation scale than shown in (5.12). For f 0 /N = 10 10 GeV and Ω θ = Ω CDM , were it not for a dynamical f , the upper bound on H inf is actually ∼ 10 2 GeV; see Figure 3 of [16]. On the other hand when f is allowed to vary, the isocurvature bound constrains the rescaled inflation scale (f 0 /f inf )H inf , 9 and thus we have (f 0 /f inf )H inf 10 2 GeV. Hence for the parameters (5.18), the isocurvature constraint is satisfied if f shrinks as i.e. f inf /N 10 23 GeV. Although such a large decay constant may be difficult to obtain in a controlled setting [49,50], here let us proceed with this bound.
The large hierarchy between f inf and f 0 significantly enhances the small-scale fluctuations. For the mode k f , cf. (4.13), the fluctuation is . (5.20) Hence, e.g., for n = 4, P and domain walls form if the axion potential emerges before the universe expands by a/a f ∼ 10 7 after the evolution of f ceases. For n 5, the k f mode contribution to the fluctuation density (4.16) becomes at least comparable to the total density of the universe at τ = τ f , and thus the backreaction from the fluctuations needs to be taken into account.

Conclusions
While a dynamical decay constant of the QCD axion can help accommodate high-scale inflation, we have shown that it also enhances the small-scale axion isocurvature perturbations. This effect stretches out to super-horizon modes, and thus the enhanced perturbation modes will continue to enter the horizon for some time after the decay constant stops its evolution. The axion fluctuation can become much larger than the potential period, thus unless the fluctuation redshifts away by the time the universe cools down to T ∼ Λ QCD , it would lead to the formation of axionic domain walls. These walls produced by the enhanced axion fluctuations are not bounded by cosmic strings, and thus would dominate the universe even in the case of N = 1.
In our analyses we have assumed the decay constant f to vary as a power-law of the scale factor. However in actual cases the dynamics of f could be more complicated. If, for instance, f instead oscillates around its final value f 0 , then the axion perturbations may get a "kick" during the first oscillation, and be enhanced much stronger than in the cases analyzed in this paper. The oscillations of f could further induce resonant amplifications of the sub-horizon axion perturbations [25,26,27,28,29,21]. We also note that, for the evolution of f to take place not too far from the QCD phase transition, a relatively light degree of freedom is required. In a simple case, the mass of the field that sets the axion decay constant f needs to be much lighter than f itself. Such a hierarchy in scales is realized in a supersymmetric axion model, where the saxion (corresponding to the radial component r, and thus setting f ) acquires a mass only from supersymmetry breaking effects. When the fluctuations of the axion are enhanced by the dynamical f , the saxion may also acquire large fluctuations, which results in large spatial fluctuations of f . It would be important to study the evolution of perturbations in explicit models of axions with dynamical decay constants.
While we have mainly focused on the QCD axion in this paper, similar discussions hold also for axion-like fields with dynamical decay constants, or more generally, for fields that have kinetic terms with time-dependent coefficients. If such a field possesses a periodic potential at energy scales higher than the QCD scale, then even if the evolution of the decay constant happens in the very early times, the enhanced perturbations could still be transformed into domain walls before redshifting away. (However in such cases one may have to take into account the effect of the potential, and in particular the induced mode-mode coupling, during the evolution of the decay constant.) We also note that, depending on when the decay constant evolves, the enhanced axion(-like) isocurvature modes may leave observable signals for upcoming CMB and large scale structure surveys.