Axion Fragmentation

We investigate the production of axion quanta during the early universe evolution of an axion-like field rolling down a wiggly potential. We compute the growth of quantum fluctuations and their back-reaction on the homogeneous zero-mode. We evaluate the transfer of kinetic energy from the zero mode to the quantum fluctuations and the conditions to decelerate the axion zero-mode as a function of the Hubble rate, the slope of the potential, the size of the barriers and the initial field velocity. We discuss how these effects impact the relaxion mechanism.


Introduction
Axion-like particles (ALPs) are ubiquitous in a large number of high-energy completions of the Standard Model (SM). They are for instance generic predictions in the low-energy spectrum of string compactifications [1,2]. ALPs denote any particles which are a pseudo-Nambu Goldstone Bosons enjoying a discrete symmetry. They are light and have weak self-interactions. Interestingly, they can be used to solve open problems of the SM as they are well-motivated candidates to explain e.g. dark matter [3][4][5], inflation [6][7][8][9], and baryogenesis [10,11]. They were originally introduced to solve the Strong CP problem [12,13] and still remain the most popular explanation to this puzzle. Independently of their virtues in providing solutions to some major open problems in the SM, they are also interesting purely from their specific phenomenological properties, notably in cosmology. They can leave imprints in many different ways such as in the Cosmic Microwave Background, in Large Scale Structures, in Gravitational Waves, in stellar astrophysics and they can be searched for directly in dedicated laboratory experiments such as haloscopes (see e.g [14,15] for recent reviews). The role of axions in cosmology has been the subject of a large number of studies. Specifically, zero-mode axion oscillations around the minimum of the axion potential can provide a large component of the energy density of the universe and mimic the effect of dark matter. Surprisingly, the effect of axion quantum fluctuations in the early universe has been mostly overlooked so far in the literature where only the zero mode homogeneous mode has been considered. In this paper, we investigate in details the effect of axion particle production during the evolution of the homogeneous zero-mode.
The production of axion quanta when the zero-mode field oscillates around one minimum of the potential is generally suppressed unless the initial position of the field is extremely close to the maximum of the potential [16]. 1 In this work, we are instead investigating the exponential production of axion quanta when the axion is rolling down its potential with a large velocity and the axion is crossing a large number of minima/maxima during its evolution. Such situation was somehow rarely investigated before, although it is quite generic. We are interested in the resulting friction force on the zero-mode. Our findings do not require any tuning of initial conditions. Our only assumption is the existence of a small slope together with some wiggles, and an initial kinetic energy larger than the barriers potential energy. In fact, this is precisely the type of potential introduced in the relaxion mechanism proposed in [20] to resolve the electroweak scale hierarchy problem. This type of axion potential, with a linear term plus a cosine, was first considered in the context of string cosmology, in models of axion monodromy [21,22] In this paper, we show that the production of axion particles generates a friction that decelerates the rolling of the field. This occurs if the field velocity is large enough to overcome the sinusoidal term and the field goes over a large number of wiggles. As we will see, the equation of motion for the axion fluctuation can be described by the Mathieu equation and parametric resonance gives an exponential production of the particles with specific wave numbers. We denote this phenomenon as axion fragmentation. This provides a novel mechanism to stop axion rolling. We focus our attention on the most dramatic effect of fragmentation, i.e., the regime in which the axion stops its motion by transferring all of its kinetic energy to the fluctuations. Other effects can arise from fragmentation. For example, in situation in which a field rolls down a steep potential, fragmentation can be a way of generating a slow-roll regime not sustained by the sole Hubble friction. At the same time, the effect discussed in this paper is mathematically very similar to the amplification of fluctuations in oscillons [25][26][27][28][29] 3 .
The phenomenon of scalar field fragmentation has been studied in the context of preheating [30][31][32][33][34]. The excitation of gauge field quanta from an axion has been studied extensively in different contexts, mainly in axion inflation models where the axion has a coupling to gauge fields, see e.g. [6,8]. We provide a model-independent detailed analytical treatment of axion fragmentation that can be applied to various setups. We discuss the precise conditions for fragmentation to stop the field (even far away from the global minimum), which we check over a numerical analysis. We discuss the implications of these findings for the relaxation mechanisms of the electroweak scale. In this case, the relaxion scans the Higgs mass-squared term in the early universe, and dynamically realizes an electroweak scale which is suppressed compared to the cutoff scale. A key ingredient of this scenario is the friction that slows down the relaxion field rolling. In the original paper [20], Hubble friction is responsible for the slow roll of the relaxion. Alternatively, a coupling to SM gauge bosons φF µνF µν can provide the necessary friction, through tachyonic particle production [35][36][37][38][39][40]. In [41] the necessary friction is provided by parametric resonance 2 In Ref. [21], particle production was not alluded to. In the follow-up paper [22], the axion is the inflaton and axion particle production was mentioned, but in the case of small wiggles, which does not lead to exponential particle production. Friction from backreaction on the zero mode was therefore not discussed. Later, Ref. [23] studied a quadratic potential with wiggles. The particle production, from the zero-mode oscillations around the global minimum with large amplitude, was discussed at the linearized level. However, the backreaction on the zero mode was not. In a follow-up study [24], the particle production was discussed by non-perturbative numerical analysis. The aim is not to exploit this effect as a stopping condition (the field eventually stops at a global minimum) but as a dark matter production mechanism. It is there suggested that axions quanta produced during the rolling stage could be (warm) dark matter candidates.
3 A recent analysis [29], although not related to axions, suggests that quanta produced by oscillons could be the dark matter candidates. of the Higgs zero mode. In [42] the field is stopped thanks to a potential instability. Finally, in [43] the relaxion is slowed down via the production of dark fermions. So far, the excitation of axion particle themselves, was not considered, although they are present in the most minimal models where the relaxion has no extra couplings to gauge fields. We analyse in details how this can be used as a stopping mechanism for the relaxion in our companion paper [44], the results of which are summarised here.
The paper is organized as follows. In section 2, we describe our setup and derive the equation of motion of the fluctuations at the leading order. We also give an intuitive discussion on the particle production effect. Then, we derive semi-analytic formula in section 3, and show numerical results in section 4. We discuss the analysis beyond the leading order in section 5. In Sec. 6 we discuss how our results apply to the relaxion mechanism. Our conclusions are drawn in Sec. 7. In the appendices A and B we provide further details on the derivation of our results.

Axion fragmentation in a nutshell
In this section we discuss how the axion field φ evolution is affected by the axion fragmentation phenomenon. The dynamics of axion quantum fluctuation is described by a Mathieu equation with time varying coefficients, from which we can estimate how the fluctuations back-react on the zero mode. Our goal is to study the dynamics of axion particle production when the axion field is rolling down a wiggly potential. For concreteness, we consider the following potential and we assume the height of the barrier Λ b as a constant for simplicity (see Fig. 1 for a sketch of the potential). For simplicity, we do not include any cosmological constant term in Eq. (2.1), and we assume that the vacuum energy is cancelled in the late universe by some other mechanism, about which we remain agnostic. We are interested in the case in which the barriers are large, i.e.
which corresponds to say that the potential has local minima. Additionally, we assume that the kinetic energy of φ is large enough to overcome the barriers,φ > Λ 2 b . The equation of motion of the axion φ(x, t) is given bÿ Figure 1. Sketch of the axion potential given in Eq. (2.1). Once axion fragmentation starts, the field φ takes a time ∆t frag and moves a distance ∆φ frag until it gets trapped in one of the wiggles.
where a is the scale factor of the Friedmann-Lemaître-Robertson-Walker metric and H =ȧ/a is the Hubble expansion rate. Let us decompose φ(x, t) into a classical homogeneous mode φ(t) and small fluctuations δφ (with no risk of confusion, we will denote the homogeneous mode as φ(t)): with a k , a † k being respectively the annihilation and creation operators which satisfy 5) and the initial condition of the mode function u k at t → −∞ is given by In the analysis of this work, we treat δφ as small perturbation. We discuss the validity of this approximation in Sec. 5. By using this approximation, we expand the last term of the LHS of Eq.
The third term of this expansion gives the dominant source of the backreaction to the zero mode from the particle production. The equations of motion of φ and δφ are given byφ where . . . indicates a quantum expectation value. Equations (2.7) and (2.8) can be rewritten in terms of the mode functions u k as Let us denote withφ 0 the initial velocity of the field φ. We assume thatφ 0 is large enough to overcome the barrier, i.e.,φ otherwise φ is trapped in the first valley. This marks a crucial difference with the wellstudied case of parametric resonance due to a scalar field which oscillates coherently at the minimum of its potential [16,[25][26][27][28][29][30][31][32][33][34]. Our case of study is sketched in Fig. 1. The field φ rolls over many wiggles until it gets trapped.
To get insights about the axion particle production, let us estimate the effect of friction. In the limit of constantφ and H = 0, Eq. (2.10) is simplified to the Mathieu equation [45]:ü Solutions to the Mathieu equation have instability and grow exponentially if the parameters are in some specific regions. 4 Forφ 2 0 > Λ 4 b , the instability region presents a band structure, as it is shown in Fig. 2. In the limitφ 2 0 Λ 4 b , the solution has an instability if the momentum k is close to nφ 0 /2f with integer n ≥ 1. For n ≥ 2, the speed of the growth is slow and the size of the instability band is small. Hence, the instability band with n = 1 gives the most important source of the friction to decelerate the axion rolling, which is given bẏ Figure 2. In grey, the instability bands of Eq. (2.11). Inside those bands, the solution grows exponentially. For the calculation of the boundaries of the instability bands, see e.g., Ref. [46].
Equivalently, one can write the instability band as |k −k cr | < δk cr forφ 2 Λ 4 b , where k cr and δk cr are defined as where, initially,φ =φ 0 . Inside the instability band, the asymptotic behavior of u k at large t is given by . (2.14) Let us now estimate the energy of the growing modes. The number of modes which exponentially grow per unit volume is ∼ (k 0 cr ) 2 δk 0 cr . The energy density of the fluctuations is As long asφ is constant, this grows as The homogeneous mode gradually looses its kinetic energy because of back-reaction, and the instability band moves towards the region of small values of k (see Fig. 3). The exponential growth of the modes with the wave number k 0 cr stops when this mode goes out from the instability band. At that time, the critical wave number has changed by δk 0 cr . Using the definition of k cr , the kinetic energy of the zero mode The timescale that the mode with wave number k 0 cr spends inside the instability band can be estimated combining Eqs. (2.19) and (2.16): (2.20) The time evolution of kinetic energy dρ/dt is roughly given as ∼ δK/δt amp : where we have dropped the subscript 0 because this equation of motion is now valid at any velocity. Eq. (2.21) can be integrated exactly fromφ 0 to 0, giving the time ∆t frag and the field excursion ∆φ frag from the beginning of particle production until the field stops: Hereφ 0 is the velocity of φ at the beginning of the particle production. In the equations above we neglected O(1) factors as the calculation above was only approximate.
In the next section, we will derive these expressions using a more precise treatment.

Analytical discussion
In what follows we discuss in detail the axion fragmentation dynamics introduced in the previous section. We establish the conditions to decelerate an axion field uniquely due to particle production friction from the axion field itself. The approximate analytical formulae derived here will be compared with the numerical solutions in the next section.
In the intuitive discussion in Sec. 2 we considered the limit in which the Hubble expansion is negligible, i.e., H = 0. Before deriving the conditions to make the field decelerate, let us consider the effect of the Hubble friction in the equation of motion for the fluctuation in Eq. (2.10). We can anticipate two additional effects once we consider the cosmic expansion. Most importantly, the growth of the modes is suppressed by the friction term 3Hu k in Eq. (2.10). In addition, since the instability band moves towards smaller values of momentum when the zero-mode decelerates, Hubble expansion makes a given mode to spend more time inside the instability band due to the red-shift of the physical momentum k/a(t).
Let us assume for the moment H to be constant. By definingũ k ≡ e 3Ht/2 u k , Eq. (2.10) can be rewritten as For H = 0, the last equation simply reproduces Eq. (2.11). According to Eq. (2.14), the exponential growth of u k is at most exp(Λ 4 b t/2fφ). Thus, in order forũ k to grow, the Hubble expansion rate should be bounded by Equivalently, the last equation can be rewritten aṡ is the slow-roll velocity of the field in the linear potential −µ 3 φ for a constant Hubble rate H. With a slight abuse of notation, in the following we will use this definition also in the case in which H is not constant or the potential is not linear, and this quantity will turn out to be useful, even without representing a proper slow-roll velocity. In addition, to have particle production active, the field should go over the barriers, H <φ/2πf . 5 However, this condition is trivially satisfied when both of Eq. (3.2) andφ > Λ 2 b are satisfied. Hence, we assumė The first of this assumption is valid until the field keeps rolling. Asφ ≈ Λ 2 b , the field stops. Contrary to the former, the second becomes easier to satisfy as the velocity decreases, therefore it is enough to assume that it is valid for the initial conditions. The assumptions above allow us to simplify the analysis due to the following three reasons. First, we can regardφ as a smooth function of the time t. The numerical solution ofφ has a smooth component and a rapidly oscillating component with frequencyφ/2πf . This oscillating component is caused by the wiggles and its relative size compared to the smooth component is ∼ Λ 4 b /φ 2 . We then neglect this oscillating term in this section. Second, we can assume thatφ is constant during the amplification time δt amp , which we defined before as the time it takes for the mode k cr to exit the instability band. This can be calculated, using Eq. (2.18), as Using the result in (3.6), we impose that | ... φ |δt amp |φ|, which we will justify later in Eq. (3.24). This condition can be rewritten as which shows that ifφ depends onφ polynomially, this condition is satisfied if Λ 4 b φ 2 . As a third simplification, we can drop the friction term 3Hu k in Eq. (2.10) without changing the physical momentum by O(1) fraction during the amplification. Now that we have presented our simplifying assumptions, let us go back to discuss the equation of motion. For a given velocityφ, once cosmic expansion is taken into account, the critical mode k cr and the width δk cr are obtained dividing the left hand side of Eq. (2.13) by the scale factor a: . (3.8) The condition H <φ/(2πf ), which was discussed above Eq. (3.5), ensures that k cr > a H. In other words, the amplification process takes place and ends when the modes are well inside the horizon. This marks an important difference with the case in which quantum fluctuations grow until they exit the horizon, as it happens, for example, with the gauge fields generated at the end of axion inflation [6]. For our discussion, we do not need to specify when the axion fragmentation dynamics takes place. The latter can be embedded in the cosmological history of the universe at different epochs. The usual dynamics of the modes crossing the horizon during inflation and then re-entering after the Big Bang does not affect our discussion. By taking the initial condition Eq. (2.6) and assuming constantφ, the asymptotic behaviour of u k with k/a = k cr /a =φ/2f after the amplification is Equation (3.9) is derived in Appendix A in the following way. The equation of motion Eq. (3.1) is solved by means of a WKB approximation in three separate time intervals, before the mode k cr enters the instability band, when the mode is well inside the instability band and after it has left it. In the two transition regions, when the mode enters and exit the instability band, the solution is found in terms of Airy functions. Finally, the five intervals are matched by imposing the continuity of the solution. From Eq. (3.9) we see that for exponential particle production occurs, one needs: where the factor of 2 difference compared to Eq. (3.9) depends on the fact that the energy density scales with |u k | 2 . We discuss the validity of this assumption later, around Eq. (3.20). Using Eq. (2.4) and Eq. (3.9), we can estimate the energy density per volume in comoving momentum space right after the end of the amplification as Using the initial condition for u k in Eq. (2.6), the energy density before the amplification is Then, the energy of the fluctuation is amplified by a factor of When this factor is much larger than 1, particle production is efficient to provide the friction for the homogeneous mode. The increment of the energy density of the fluctuations because of the particle production is where dk cr /dt is the velocity of the instability band in the momentum space. For non-zero H, k cr is determined as Therefore, the velocity of the instability band is given as dk cr /dt = a(φ + Hφ)/(2f ). Thus, (3.14) The kinetic energy of the homogeneous mode isφ 2 /2 and its time derivative isφφ. As a result, from conservation of energy, we obtain the following equation: Compared to Eqs. (2.7) and (2.9), Eq.(3.15) describes the effect of fragmentation after averaging over many oscillations of the sinusoidal potential, and is more suitable for our analysis. Using this equation we will determine the general conditions to stop the field due to axion fragmentation, which are obtained in the next sections. The first and second term in the right-handed side of Eq. (3.15) are the effect of Hubble friction and acceleration by the slope, respectively. This equation can be regarded as a consistency condition forφ during the fragmentation phase. By solving the above equation,φ can be calculated as a function ofφ, Λ b , f , µ 3 , and H.

General condition to stop the axion
Let us discuss conditions to stop the axion field. Here we summarize the results, while the details of the derivation are given in Appendix B. φ < 0 must hold from the initial time until the field has come to a complete stop. As detailed in Appendix B.3, this is realized if (and only if) the following condition holds for the initial velocity: Here W n (x) is the product logarithm function (also known as Lambert W function), 6 whose real branches are plotted in Fig. 4. Eq. (3.16) expresses an equilibrium between the slope and the Hubble expansion which allows for efficient fragmentation: if the slope increases, in order to avoid the field acceleration, the friction due to cosmic expansion should compensate this effect. Alternatively, forφ 0 < µ 3 /(2H) = 3/2φ SR , one can see Eq. (3.16) as a lower bound on Λ b that expresses, for given µ 3 and H, the necessary amount of fragmentation needed in order to slow down the field. If Eq. (3.16) is satisfied, Eq. (3.15) has only one solution, with negativeφ, which 6 The product logarithm function is the inverse function of W e W = x. In general, there exist infinite number of solutions for this equation, and W 0 and W −1 are the two real ones. In Fig. 4. This function is available ProductLog in Mathematica or special.lambertw in SciPy. See e.g., http://mathworld.wolfram.com/LambertW-Function.html. is given by 7φ Here a and b are dimensionless parameters which are defined as Let us discuss the validity of the assumption Eq. (3.10). The effect of the axion fragmentaion in Eq. (3.15) has a exponential factor with a exponent N ≡ πΛ 8 b /2fφ 2 |φ + Hφ|. As we have discussed, for exponential particle production to occur, N should be larger than 1. By using Eq. (3.17), we obtain As we can see in Fig. 5, for fixed a, N is an monotonously increasing function of 1/b. For a 1, as we will always assume, N becomes small at 1/b 0, and is well 7 Before fragmentation is active, the field is only subject to its potential and to Hubble friction, and its equation of motion is simplyφ = µ 3 − 3Hφ, where we neglected a small oscillating term. This cannot be obtained as the Λ b → 0 limit of Eq. (3.17). The reason is that the equation of motion Eq. (3.15) was derived assuming fragmentation is active. In particular, we assumedφ Λ 4 b /f . The acceleration is initiallyφ = µ 3 − 3Hφ, and asymptotes to Eq. (3.17) as fragmentation starts. approximated as N −b. Thus, by requiring N > 1, we obtain is positive and the inequality is satisfied. Thus, the above condition can be rewritten aṡ As long as Eq. (3.20) (or equivalently Eq. (3.21)) is satisfied, we can safely useφ given in Eq. (3.17). Note that 2Hφ − πΛ 8 b /(2fφ 2 ) is a monotonously increasing function ofφ. Thus, we can see that Eq. (3.21) is satisfied for anyφ <φ 0 if (and only if) If this condition is satisfied, Eq. (3.17) for the accelerationφ can be used to describe the fragmentation process from its begininng until the end of that. In phenomenologically interesting applicationsφ ≤φ SR , and Eq. (3.22) is always satisfied. To summarize, in order for fragmentation to stop the rolling of the axion field we need to simultaneously impose the conditions of Eq. (3.3) and Eq. (3.16) to stop the axion. We are interested in the case in which the potential has local minima (i.e. Λ 4 b /f > µ 3 ), and we will assumeφ 0 ≤φ SR . In this case, Eq. (3.3) is trivially satisfied, and the only condition that must hold is Eq. (3.16). Figure 6 and Fig. 7 show examples of the parameter space in which Eq. (3.16) is satisfied. In the case ofφ 0 φ SR , there is an upperbound on µ 3 because the acceleration effect by the slope should be weaker than the particle production effect.
On the other hand, the bound on µ 3 is very weak forφ φ SR and it becomes trivial . We show this velocity in red in Fig. 6 and Fig. 7. As long asφ SR f 2 , the modification to the slow roll velocity is small, i.e., φ SR(frag) φ SR .

Stopping conditions in several limits
In the previous section, we described the generic negative solution ofφ and the stopping conditions. Here we discuss several cases in which the conditions can be simplified by taking some limits of the parameters.
modified slow roll velocity H < H cr No exponential particle production  No exponential particle production This can be regarded as a refinement of Eq. (2.21). In this case, the amplification factor of the fluctuation energy in Eq. (3.12) is given by In the derivation of Eq. (3.9) we assumed thatφ f 2 , which is necessary for the validity of the low-energy EFT of the axion field, thus the amplification factor is much larger than 1, enhancing the efficiency of fragmentation. The time scale and the field excursion during the axion particle production are where we dropped subleading terms. The number of wiggles which the axion travels until it stops is For example, forφ 0 /Λ 2 b = 10 2 and log (32π 2 f 4 /φ 2 0 ) = 10, this number is ∼ 3 × 10 7 .  3.2.2 H 0 and µ 3 = 0 Next, let us discuss the case in which µ 3 = 0 and H is small enough to be neglected. Forφ 0 φ SR , Eq. (3.16) can be simplified as This corresponds to the left part of Figs. 6 and 7, which show how H is irrelevant ifφ 0 φ SR(frag) ≈φ SR . In the limit of H = 0, the negativeφ solution Eq. (3.17) is simplified asφ where a, b are defined in Eq. (3.18). In the limit µ 3 = 0 we recover Eq. (3.24). The acceleration |φ| monotonically decreases as a function of µ 3 . In Fig. 8, we show the relative variation ofφ when the slope µ 3 goes from 0 to µ max defined in Eq. (3.29).
As long asφ < f , the decrements ofφ is at most ∼ 15 %. Hence, the estimates for 8 Eq. (3.24) can be obtained by separately taking the limits µ 3 − 2Hφ → 0 + and µ 3 − 2Hφ → 0 − of the two branches of Eq. (3.17), respectively, and using the expansions the fragmentation time and the total field excursion of Eqs. (3.26, 3.27) still provide a reliable approximation.

3.2.3φ =φ SR
Let us discuss the case in which the initial velocity is equal to the slow roll velocity, i.e.,φ 0 =φ SR . This is the case if the dynamics takes place during inflation, since the velocity is exponentially driven to the attractor slow-roll velocity irrespectively of the initial conditions. In this case, the condition Eq. (3.16) can be simplified to the following equation: (3.31) By replacingφ SR = µ 3 /(3H), the above condition Eq. (3.31) can be equivalently rewritten as a condition for µ 3 and H for given Λ b and f : This condition tells us that the Hubble friction is required to prevent acceleration in order to work the fragmentation mechanism if the slope is steeper than the threshold value µ 3 th . Alternatively, one can think of Eq. (3.31) as a lower bound on Λ b to have enough fragmentation to stop the field, for given µ 3 and H.

Numerical analysis of the equations of motion
In this section, we test the analytical understanding developed in Sec. 3 against a numerical solution of the equations of motion for the homogeneous mode φ(t), Eq. (2.9), and for the fluctuations δφ(x, t), Eq. (2.10). We still limit ourselves to a linear level analysis, the validity of which will be further discussed in Sec. 5. For this calculation, we discretize the integral in Eq. (2.9) and take 10000 modes whose  momentum is between 10 −4 k 0 cr and k 0 cr + 10δk 0 cr . The momenta are evenly spaced in logarithmic scale. The differential equations are solved numerically by the fourth order Runge-Kutta method.
In Fig. 10   The phase diagrams of the axion particle production are shown in Figs. 14 and 15 for general values of µ 3 and H. In these figures, we take φ = −πf /2 as the initial condition so that Λ 4 b cos φ/f = 0 at the beginning. In Fig. 14, we take H = 0 and show the parameter region in which the particle production is efficient inφ-µ 3 plane. The figure shows that the condition Eq. (3.29) successfully reproduces the numerical result for the maximal slope µ 3 that allows stopping due to fragmentation. In Fig. 15, we take nonzero H and show the parameter region in which the particle production is efficient in theφ SR − µ 3 plane. The figure shows the excellent agreement between Eq. (3.31) and the numerical results.
In Fig. 16, we show the time evolution of the energy spectrum of the fluctuations. This quantity can be easily estimated as follows. For H = 0 and µ 3 = 0, by using   Then, we can calculate the energy spectrum after fragmentation as where we dropped the oscillating terms. Fig. 16 shows that this estimation agrees with the result of the numerical calculation. From Eq. (4.2), we see that the peak frequency coincides with the initial position of the instability band, k 0 cr =φ 0 /(2f ). Sinceφ 0 > Λ 2 b , the emitted particles are typically relativistic. We expect that nonperturbative effects will broaden this spectrum (see the discussion in the next section).

Beyond the perturbative analysis
Let us comment on the validity of the leading order expansion of δφ. For H = 0 and µ 3 = 0, the size of fluctuation which is generated by the axion particle production can be estimated by using Eq. (4.1) as where we averaged out the oscillating term in the integral. The time evolution of the ratio between δφ 2 and f is shown in Fig. 17. We can see that δφ 2 at late times becomes large and we need to use non-perturbative methods for a concrete analysis in this regime. However, Fig. 17 shows that δφ 2 /f O(1) is satisfied during most of fragmentation process. Thus, we can expect that the estimation on the time scale Eq. (3.26) and field excursion Eq. (3.27) during fragmentation do not considerably vary from the ones obtained using a non-perturbative analysis, unless non-linear effects lead to a significant suppression of the growth of the field perturbations. This suppression is not expected to happen in our scenario due to the periodic potential, as we elaborate in the next paragraph.
A non-linear analysis would require a dedicated lattice study, on which we are currently working [47]. Preliminary results indicate that our main quantitative result, namely the estimate of the total field excursion before the axion is stopped Eq. (3.27), is correct up to a factor of order 1. The main difference that we expect to emerge from a non-linear analysis is the spectrum of excited fluctuations, which will be broadened by interactions among different modes. Moreover, the final field configuration will be inhomogeneous, with the scalar field populating more than one minima of the potential. This observation does not change significantly the picture we described so far, since the spread δφ 2 will always be much smaller compared to the total field excursion ∆φ frag . Interestingly, formation of domain walls can occur in this scenario, depending on the non-linear evolution of the system and on other model dependent inputs such as the axion lifetime. We postpone this discussion to a future publication [47].
A thorough comparison with other resonant systems is non-trivial, and it can hardly provide any insight on the non-linear behaviour of our model. Many models which were studied in the context of preheating feature a strong suppression of the growth at the non-linear level, generically due to the appearance of large effective mass terms. As an example, in Ref. [48], the resonant production of scalar particles χ from the oscillations of the inflaton ϕ is studied, and it is shown that both the perturbations of the inflaton field δϕ (which are generated non-linearly) and the inclusion of the quartic coupling χ 4 , suppress the further growth of χ modes through a mass term ( δφ 2 + χ 2 )χ 2 (with appropriate coupling constants). Moreover, as fluctuations grow and drain energy from the zero-mode, the amplitude of the latter decreases, and hence the force driving the growth also progressively decreases. The case under study here has two peculiarities with respect to the one above. First, since the field traverses many periods of the periodic potential, the oscillating term that stimulates the growth of fluctuations has effectively a constant amplitude. Secondly, because of the approximate shift symmetry, no effective mass scale is generated in our case. Instead, all corrections enter the equation of motions only through the cosine potential, and thus there is no reason to expect any suppression. A setup similar to ours is discussed in [24], in which a monodromy potential m 2 φ 2 + Λ 4 b cos φ/f is studied. The paper shows that the evolution of the zero-mode stops shortly after the fluctuations have entered the non-linear regime, in accordance with our expectations.

Consequences: Relaxation of the electroweak scale
The axion fragmentation dynamics explored in this work should be taken into account in the evolution of any axion field which rolls down a wiggly potential. This phenomenon can fundamentally impact on a broad range of models, such as axion monodromy constructions and relaxion scenarios. In this section, we consider the effects of axion fragmentation on the relaxation mechanisms of the electroweak scale.
The relaxion mechanism is a solution to the electroweak hierarchy problem in which the Higgs mass term is controlled by the evolution of an axion-like field, the relaxion [20]. This field evolves classically in the early universe until it stops close to a critical point, defined as the field value at which the Higgs VEV is zero. A key ingredient in this picture is a potential that features periodic wiggles, similar to the one discussed in this work. Relaxion fragmentation affects this construction in a substantial way [44], as we detail below in the two main implementations of the relaxion idea which have been discussed in the literature so far.

Higgs dependent barriers
In the original proposal [20], the cosine term in the relaxion potential Eq. (2.1) has an amplitude dependent on the Higgs VEV, Λ 4 b ∝ h n , with n = 1, 2. For n = 1 (QCD relaxion), Λ 4 b ∼ m q Λ 3 QCD where m q is the light quark masses. For the case with n = 2, the scale Λ b cannot be far from the electroweak scale, satisfying Λ b TeV. The potential contains an interaction between the Higgs and the relaxion: Initially, the Higgs mass term µ 2 h ≡ Λ 2 − g Λφ is positive, and the VEV is zero. As soon as φ > Λ/g , µ 2 h turns negative, a VEV develops, and the cosine term grows. In Ref. [20], it is assumed that the entire evolution takes place during a long period of inflation, and that Hubble friction is strong enough to stop the field as soon as the wiggles become larger than the average slope and the potential develop local minima, i.e. for gΛ 3 ≈ Λ 4 b /f . In particular, this happens when the time that it takes to roll over one period of the cosine term is longer than one Hubble time, i.e. for Relaxion fragmentation offers an additional source of friction for the relaxion rolling. As discussed in Ref. [44], this opens up two possibilities: on the one hand, the relaxion can be stopped by fragmentation even when Eq. (6.2) is not satisfied. On the other hand, it is possible to stop the relaxion field with a much shorter period of inflation or even in the absence of an inflationary background, with a negligible Hubble friction. This opens new possibilities for relaxion model building, independent from constraints on the inflationary sector. If relaxation takes place after inflation, it is possible, at least in principle, to conceive a model in which this phase has observable features, most probably in gravitational waves.This study could open the way to observable relaxion models.
Higgs independent barriers An alternative relaxion construction was proposed in [35], in which the amplitude of the cosine term is Higgs-independent, and the friction is mainly provided by gauge boson particle production. The relaxion couples to the Chern-Simons term of the massive SM Z boson, through a term In the presence of this coupling, the equation of motion for the transverse polarization of the Z has a tachyonic instability for small mass m Z Contrarily to the case discussed above, initially the Higgs has a large VEV, and the SM particles are heavy. As the relaxion approaches the critical point and the gauge bosons become lighter, the tachyonic instability is triggered and the relaxion kinetic energy is dissipated through the production of Z bosons. Fragmentation poses a serious threat to this model [44]. Since the amplitude of the cosine term is constant, fragmentation is always active, and the relaxion can be slowed down and stopped when the Higgs mass is large and close to the cut-off Λ, thus spoiling the successful relaxation of the Higgs VEV to its current value. In particular, • if relaxation takes place after inflation [39], the parameter space is restricted by the condition of avoiding excessive fragmentation. Moreover, once cosmological constraints are taken into account, the mechanism is excluded at least for a cutoff larger than few TeV.
• If relaxation happens during inflation, the constraints from fragmentation reduce the available parameter space but do not exclude the model. The dark matter scenario discussed in [40], in particular, is not affected.

Summary and outlook
In this paper, we discussed the production of quantum fluctuations during the evolution of an axion-like field rolling down a potential featuring wiggles, as given in Eq. (2.1). We refer to this effect as axion fragmentation. While the production of quanta is suppressed when an axion oscillates around the minimum of its potential, unless the initial amplitude is very large and the initial position of the field is tuned close to the maximum of the sinusoidal potential, the effect is very large in the case where the axion field crosses many of its maxima. We studied in detail under which conditions axion fragmentation can efficiently stop the evolution of the field. We computed the time scale needed for stopping and the corresponding field excursion. The wavefunction of the fluctuations obeys the Mathieu equation and the energy of the modes within the instability band Eq. (2.12) grows exponentially. If both the slope of the potential and the Hubble expansion rate are sufficiently small, this particle production effect decelerates the homogeneous mode. The condition is given as Eq. (3.3) and Eq. (3.16) in terms of the initial field velocity, the linear slope of the potential, the Hubble rate, the size of the barriers and the periodicity of the sinusoidal potential. The corresponding accelerationφ is given in Eq. (3.17).
Axion fragmentation is a generic effect which can have interesting phenomenological implications. It is particularly relevant for the mechanism of cosmological relaxation of the electroweak scale. We dedicate a separate paper to study in details these implications in Ref. [44], where we conclude that new regions of parameter space open and novel directions for relaxion model building are offered by this effect.
In the present work we study the regime in which the potential has local minima, see Eq. (2.2). It would be interesting to investigate the effect of fragmentation in the case in which the oscillating term in potential does not generate local minima. This may have implications for some relaxion models where loop effects generate small Higgs-independent barriers, so that there are constant wiggles with small amplitude during the whole scanning of the Higgs mass parameter.
Another promising direction will be to explore the impact of axion quanta on the cosmological history of the universe. As discussed in Ref. [44], depending on the equation of the state of the universe during axion rolling, the produced quanta may represent a significant fraction of the energy density of the universe. Whether they can be viable dark matter candidates, depends on the time of fragmentation. Such quanta may in turn induce gravitational waves. They may be diluted or leave observable imprints. These effects deserve detailed studies which we postpone for future work.

A Approximate solution of Eq. (2.10)
Let us discuss the evolution of the wave functions u k . The boundary condition at t → −∞ is given by Eq. (2.6). For the duration of the amplification process, we can neglect the Hubble friction term 3Hu k in the equation of motion. This is justified for k φ /2f if H Λ 4 b /φf is satisfied (see the discussion around Eq. (3.2)). The mode function u k satisfies the following equation of motion: Note that the Hubble expansion still has its effect via a −2 k 2 term although we have dropped 3Hu k .
A.1φ = 0 and H = 0 (Mathieu equation) First, let us consider the case of constantφ and H = 0 as in Sec. 2. Taking the scale factor a = 1, Eq. (A.1) is simplified as This is the Mathieu equation, and it is known that u k grows exponentially when k satisfiesφ Let us see this behavior explicitly. We define δ and for convenience: Then, Eq. (A.2) is rewritten as We assume 1 and we will expand perturbatively in . On the other hand, we assume δ = O(1), as we are interested in the first instability band. In the limit = 0, the solution is of the form u k = a cos(φt/2) + b sin(φt/2) for constant a and b. This motivates the following ansatz for u k : Here cos and sin encode the rapid oscillations of the solution, while A n 's and B n 's are slowly varying coefficients. The terms with n ≥ 2 are introduced to maintain a consistency with Eq. (A.5). We define a dimensionless time τ : . (A.7) By plugging Eqs. (A.6, A.7) into Eq. (A.2), we obtain the following differential equations for A n 's and B n 's: for n = 1 and, for n ≥ 2, The differential equations above indicate a simple hierarchy between the coefficients At the leading order in , we can neglect A n and B n with n ≥ 2, and also the second derivative in Eqs. (A.8), (A.9). We obtain For −1 < δ < 1, the asymptotic behavior of A 1 and B 1 are then .
For |δ| < 1 the solution is unstable. On the contrary, for δ 2 > 1 it oscillates as sin(

A.2 Small non-zeroφ and H
Let us now introduce a small, constant accelerationφ: φ =φ(0) +φt, (A. 16) whereφ(0) =φ(t = 0). The assumption of neglecting higher derivatives is justified in the main text, see the discussion around Eq. (3.7). Let us consider a given mode k. By a simple time shift, here we define t = 0 as the time at which k is at the center of the instability band, which is now defined as As the velocity decreases, the instability band moves to lower modes. For a given k, we will solve the equations of motion from a time slightly before it enters the instability band, until slightly after it exits, and we will see how the initial oscillatory behaviour is then amplified inside the instability band until the mode exits. We assume for u k an ansatz similar to Eq. (A.6), but this time the cosine will depend on φ(t): Again, we only consider the case in which the velocity is larger than the wiggles, In order to keep the evolution under perturbative control, we assume that the time spend in the instability band is short, and that the velocity changes only slightly during this time. Furthermore, we assume that the effect of Hubble friction is small: In Eq. (A.18) we make explicit the a −1 dependence of the mode functions on the scale factor. When taking the derivativesu k ,ü k we will keep this factor as a constant, consistently with the assumption that the amplification time for any mode k is much shorter than a Hubble time. By plugging the above u k into the equation of motion Eq. (A.1), we obtain the equations for the coefficients of the sine and cosine terms. Similarly to before, the n ≥ 2 equations show an hierarchy A n , B n ∼ (A 1 , B 1 )× n−1 . This can be extracted from the equations of motion or by analogy with the constant velocity case, remembering that the acceleration gives only a small correction to the evolution according to Eq. (A.19). Thus, we can neglect the contribution from the n ≥ 2 terms. For n = 1 we obtain the equation for the coefficients of cos φ/2f and sin φ/2f : By using the results in the previous section, we can estimate that d 2 A 1 /dt 2 , d 2 B 1 /dt 2 ∼ (Λ 4 b /φf )dA 1 /dt, (Λ 4 b /φf )dB 1 /dt. Thus, by assuming Λ 4 b φ 2 , we can drop d 2 A 1 /dt 2 and d 2 B 1 /dt 2 in the above equations. We require the time evolution ofφ is sufficiently slower than that of A 1 and B 1 , i.e., |φ/φ| |A −1 1 dA 1 /dt|, |B −1 1 dB 1 /dt|. This condition can be understood as follows. The time evolution of the energy of zero mode is given by dρ 0 /dt =φφ = ρ 0 (φ/φ). The energy stored in the fluctuations inside the instability band is Thus, dropping second time derivatives, the growth of δρ fluc is given by The growth in energy of the modes inside the instability band correspond to the slow-down of the zero mode. Thus, dρ 0 /dt ∼ dδρ fluc /dt, and we obtainφ/φ ∼ (Ȧ 1 + B 1 )/(A 1 + B 1 ) δρ fluc /ρ 0 . Due to the smallness of the band width δk, δρ fluc /ρ 0 1, and henceφ/φ (Ȧ 1 +Ḃ 1 )/(A 1 + B 1 ). The last condition can be rewritten as which impliesφ φ 2 /f . Thus, we can drop the termsφA 1 /2f andφB 1 /2f . Finally, we define the auxiliary quantities where k/a =φ(0)/2f . Furthermore, since we assume the Hubble expansion rate is sufficiently small as explained in Eq. (A. 19), we expand the scale factor as a −2 (t) a −2 (0)(1 − 2Ht). In terms of τ and z, we obtain For |zτ | < 1, (A.39) For zτ > 1, where c A and c B are constants. A i can be obtained using Eq. (A.42): Here Ai (x) = dAi/dx and Bi (x) = dBi/dx. Similarly, we obtain the generic solution for τ 1/z: where c A and c B are constants. To connect each region of τ , we will need the asymptotic form of Airy functions: and We will also need the expansion of the phases of the WKB solutions close to the critical points zτ = ±1: A.2.2 Matching with the initial condition Eq. (2.6) The solutions described have free coefficients that can be obtained from matching at the intersection of their regime of validity and with the initial condition where Eq. (A.54) is specified up to a phase, which we choose in order to make the phases of A 1 and B 1 real and positive for τ > −1/z. The initial conditions for A 1 (τ ) and B 1 (τ ) are given as where θ is a phase which we do not specify here. Matching Eqs. (A.36, A.37) with the initial condition Eq. (A.55) at τ < −1/z, we have Here we dropped exponentially suppressed term exp . For τ ∼ 1/z, B 1 can be written as a linear combination of Airy functions. Matching onto Eq. (A.60, A.61) gives Therefore, after amplification has ended, the asymptotic behavior of u k is in Eq. (3.15). For z < 0, the initial conditions for A 1 and B 1 is The coefficients A 1 , B 1 for τ < −1/|z| are given as A 1 and B 1 for τ > 1/|z| are given as Therefore, after amplification has ended, the asymptotic behavior of u k is By using the definition of z given in Eq. (A.25), we obtain Eq. (3.9).

B Detailed analysis on Eq. (3.15)
In this section, we discuss Eq. (3.15) in detail.

B.1 The solutions of Eq. (3.15)
Let us discuss the solutions of Eq. (3.15), which we report here for simplicity: By defining Eq. (3.15) can be rewritten as where a and b are dimensionless parameters which are defined as Eq. (B.3) can be rewritten as which has the following solutions: where W 0 , W −1 are the two main branches of the Lambert function (product logarithm). The first two solutions have x < 0, the other two have x > 0. Formally, the first solution only exists for a < 1. By using the consistency condition of EFT, we can assumeφ 2 < 32π 2 f 4 , i.e., 0 < a < 1. In terms of the accelerationφ, the analytic φ 1 −Hφ 0φ 2φminφ3 0 Figure 19. A typical behavior of f (φ,φ) as a function ofφ.

B.2 The condition not to have positiveφ solution
In order to stop the axion rolling, the accelerationφ should always be negative. Before fragmentation starts, the field is only subject to the slope and Hubble friction, φ = µ 3 − 3Hφ. Let us start by deriving a condition not to have a solution 0 <φ < Eq. (3.15) is equivalent to f (φ,φ) = 0. We can easily see the following property of f (φ,φ): The function f (φ,φ) monotonously increases forφ < −Hφ because ofφ 2 /32π 2 f 4 < 1.
There exists one local minimum atφ =φ min > −Hφ. A typical behavior of f (φ,φ) is shown in Fig. 19. The solutions of f (φ,φ) = 0 are classified depending on the sign ofφ min (φ), f (φ,φ min (φ)), and f (φ, 0), which are given as    Table 1. Solutions of f (φ,φ) = 0 for givenφ and their classification. Note that lation is not reliable. Still, in this regime, we can assume thatφ varies continuously. When fragmentation turns on and Eq. (B.1) becomes valid,φ will smoothly decrease from its initial value until it reaches one of the solutions in Eqs. (B.7-B.9), and it will stick to it for the rest of the evolution. In order to understand the behaviour during this phase, it is also useful to notice that, if f (φ,φ) > 0,φ will decrease. Oppositely, if f (φ,φ) < 0,φ will increase. Hence,φ 1 is a stable solution, as well asφ 3 (if it exists). On the contrary,φ 2 , if it exists, is an unstable solution.
Let us now discuss the necessary and sufficient conditions in order to guarantee that the field slows-down until it stops. The conditions in Eqs. (B.20, B.21), or equivalently in Eq. (B.26), give a bound on the slope µ 3 that depends on the velocitẏ φ. If this condition is satisfied forφ =φ 0 , initially the solutions to Eq. (B.1) havë φ < 0: .

(B.29)
This is a necessary but not sufficient condition to stop the rolling of the field φ. Indeed,φ should be negative until the axion will be stopped. In the case µ 3 < µ 3 1 (φ 0 ), the only solution of f (φ 0 ,φ) = 0 isφ 1 . Thus, at the beginning,φ becomesφ 1 (φ 0 ) andφ starts to decrease. Asφ decreases,φ changes continuously andφ =φ 1 (φ) is always satisfied even ifφ 2 orφ 3 solution appear at some velocityφ <φ 0 . As we have seen in App. B.1,φ 1 is always negative and the field will be stopped in the end.

B.4 Modified slow roll velocity
In this section, we want to determine the existence conditions and the value of a constant-velocity solution of the equation of motion, under the effect of fragmentation. Fragmentation acts as an additional source friction, thus it is clear that this velocity, which we callφ SR(frag) , will be smaller than the usual slow-roll velocitẏ φ SR ≡ µ 3 /3H. As we discussed in App. B.3, if such a solution exists the field will decelerate until it reachesφ SR(frag) , and its evolution cannot be stopped. First, let us discuss the existence condition of the solution of f (φ, 0) = 0. Let us define ζ ≡φ/φ SR > 0. Then, we can rewrite f (φ, 0) as    where, after simplifying a factor ζ, we expanded separately the right-hand side and the exponent of the left-hand side for ζ 1. This procedure is justified by the observation that the exponent is very large and thus the LHS changes much more rapidly than the RHS with ζ. By solving f (φ, 0) = 0 and (∂f /∂φ)(φ, 0) = 0 simultaneously, we obtaiṅ

(B.33)
For H > H cr , Eq. (B.31) admits two solutions. The smallest of the two has dφ/dφ > 0 and is unstable. The largest solution is instead stable, and can be regarded the modified slow roll velocity, given bẏ (B.34)