ournal of C osmology and A stroparticle hysics Gravitational waves from axion monodromy

. Large ﬁeld inﬂation is arguably the simplest and most natural variant of slow-roll inﬂation. Axion monodromy may be the most promising framework for realising this scenario. As one of its deﬁning features, the long-range polynomial potential possesses short-range, instantonic modulations. These can give rise to a series of local minima in the post-inﬂationary region of the potential. We show that for certain parameter choices the inﬂaton populates more than one of these vacua inside a single Hubble patch. This corresponds to a dynamical phase decomposition, analogously to what happens in the course of thermal ﬁrst-order phase transitions. In the subsequent process of bubble wall collisions, the lowest-lying axionic minimum eventually takes over all space. Our main result is that this violent process sources gravitational waves, very much like in the case of a ﬁrst-order phase transition. We compute the energy density and peak frequency of the signal, which can lie anywhere in the mHz-GHz range, possibly within reach of next-generation interferometers. We also note that this “dynamical phase decomposition” phenomenon and its gravitational wave signal are more general and may apply to other inﬂationary or reheating scenarios with axions and modulated potentials.


Introduction
The central predictions allowing us to discriminate between inflationary models are the slow roll parameters, most prominently the tilt of the scalar power spectrum n s 0.96 and the tensor-to-scalar ratio r 0.08 [1,2]. However, even with growing precision many models are expected to remain consistent with this limited set of data. Hence it is of great importance to identify additional predictions characterizing specific classes of models.

JCAP11(2016)003
In short, our message is the following: due to the typical instantonic modulations of the potential, first order phase-transition-like, violent dynamics may occur after the end of inflation and before reheating. This leads to additional gravitational waves, which are of course very different in frequency from those studied in the CMB. Thus, monodromy models may (in addition to or independently of their prediction of r) be established by future ground-or space-based interferometers. The final word would thus come from the new field of gravitational-wave astronomy [27] (for a recent review see e.g. [28]).
Of course, gravitational waves are a well-known signature of cosmic strings, including axionic strings. For recent work on the non-trival late-time dynamics of axionic models and gravitational waves which is closer in spirit to our proposal see, e.g. [29][30][31][32][33][34]. Emission of gravitational waves from axionic couplings in inflationary setups has recently been considered in [35,36] (see also [37,38] for constraints). Independently of the gravitational wave signal, related dynamical phenomena may also occur in the dark matter context [39,40].
Let us now explain the physics underlying our scenario in more detail: the basic building block of monodromy inflation is an axion with sub-planckian decay constant f . Nonperturbative effects induce the familiar cos(φ/f )-type potential. When the monodromy effect is included, this cosine-potential shows up in the form of modulations of the long-range, polynomial term. Even if the relative size of these modulations is small at large field values, where slow-roll inflation is realized, they can become dominant near the minimum of the polynomial potential. After inflation, the field oscillates with decreasing amplitude such that its motion eventually becomes confined to the vicinity of one of these local minima. However, due to field fluctuations, different local minima may be chosen in different regions of the same Hubble patch. In other words, the Universe is decomposed into phases.
Two comments are in order. First, the field fluctuations inducing the above phenomenon can have different origin. On the one hand, there are inflationary super-horizon fluctuations which have become classical at the time when they re-enter the horizon. On the other hand, the field is subject to an intrinsic quantum uncertainty at any given time. As we will see, this second type of uncertainty is, after parametric amplification [40], more likely to source the desired phase decomposition.
Second, while we for definiteness identify the monodromic axion with the inflaton, the phenomenon can occur also in different contexts. In particular, similar physics may arise in any non-inflationary axion model with a monodromy (the relaxion being one recent popular example [41], see also [42]). Moreover, even within the inflationary context, our proposal of 'dynamical phase decomposition' is not restricted to monodromy models. Indeed, models of 'aligned' or 'winding' inflation [43][44][45] can naturally exhibit short-range modulations on top of a long-range periodic potential. The Weak Gravity Conjecture [46] for instantons may in fact demand such modulations [47,48] and the simplest string constructions naturally provide them [49]. Furthermore, the Weak Gravity Conjecture for domain walls constrains models based on monodromy [42,50]. In particular, the size of the wiggles is bounded [50].
Let us now complete the discussion of the cosmological dynamics: after a phase decomposition has occurred, the regions with the lowest-lying populated minimum will expand. This is very similar to the way in which a strong first-order phase transition is completed through the collision of cosmic bubbles of true vacuum. However, in our case the transition occurs before reheating and very far from thermal equilibrium. It may thus be more appropriate to talk about 'dynamical phase decomposition' rather than about a phase transition in the usual sense. Nevertheless, the concept of bubble formation and collision is still appropriate in our setting.

JCAP11(2016)003
Thermodynamic cosmological phase transitions have been widely explored in various contexts (see [51][52][53] for early seminal work and [54] and refs. therein for the case of the electroweak phase transition). In particular, it is well known that they source gravitational radiation (see [55][56][57] for recent reviews and e.g. [58] for radiation from other cosmological sources such as cosmic strings and preheating). We will rely on these results. This paper is structured as follows: section 2 provides the basic setup. More specifically, section 2.1 introduces our axion monodromy setting with dominant quadratic potential and a series of local minima, section 2.2 briefly discusses reheating, and section 2.3 explains the post-inflationary dynamics and the phase decomposition. In section 3 we estimate the probability that a phase decomposition occurs as a consequence of field fluctuations. In particular, we treat inflationary fluctuations in section 3.1 and the intrinsic quantum uncertainty in section 3.2. In section 4 we address the crucial issue of the possible enhancements of fluctuations due to background-field-oscillations in the modulated potential. Ultimately, in section 5 we estimate the spectrum and abundance of the gravitational radiation produced during the phase transition before we conclude in section 6. Additionally, we devote appendix A to a more detailed discussion of scalar field fluctuations after inflation. Appendix B provides some details of an inflection point model of inflation in which our gravitational wave signal may also arise.

Phases from axion monodromy
In this section we introduce a string-motivated scenario of the inflationary universe. It is based on the framework of Axion Monodromy Inflation [3,4]. We begin by explaining the basic features of the axion potential and the possibility of having coexisting populated axionic vacua.

Local minima in axion monodromy
The inflaton potential of a model of axion monodromy inflation contains an oscillatory term, which respects a discrete shift symmetry, and a polynomial term, which breaks it explicitly. In this work we will take the polynomial term to be quadratic, as this will be sufficient to demonstrate the effect we wish to study: The potential (2.1), plotted in figure 1, can have local minima, whose existence depends on the values of the prefactor Λ 4 , the so-called axion decay constant f and the inflaton mass m. Although we have in mind the specific case in which φ is the inflaton, many of our considerations apply to the case of a generic axion-like field with potential (2.1). 1 We begin with an analysis of the classical evolution of φ. The equation of motion of φ reads:φ

JCAP11(2016)003
Therefore, the amplitude of φ decreases. This conclusion remains valid even for time-varying H. In fact, since H(t) ∼ ρ 1/2 φ /M p , the Hubble rate decreases as the amplitude of φ falls. Once the amplitude is sufficiently small, the cosine oscillations in (2.1) cannot be neglected any longer. Eventually, the field is caught in one of the cosine wells. One can observe that the field is more likely to get trapped in one of the lowest-lying minima. This can be understood as follows: the friction term in (2.2) becomes less and less relevant as H(t) falls. This implies that the fractional energy loss per oscillation also decreases. Therefore, even though the field can in principle get stuck at any time, it is more likely to do so late in its evolution, when it is oscillating near the bottom of the well containing the lowest-lying minima. 2 For this reason, we will mostly focus on the last two wells.
The existence of different local minima implies that the universe is potentially decomposed into several phases. This happens if the field settles in different minima in different parts of the universe. Due to fluctuations, the scalar φ can end up in one or the other minimum in different regions of the same Hubble patch. In this paper, we are interested in studying the conditions under which such a phase decomposition can take place.
If such a phenomenon occurs, eventually the field will settle in the state of lower energy as a consequence of the expansion of bubbles containing the true vacuum. This corresponds to a phase transition, which can in principle have strong cosmological signatures, above all the radiation of gravitational waves. We would like to provide an estimate for the spectrum of gravitational waves produced in such an event (see also [29][30][31][32][33] for related work, but in different contexts).
Before moving on to study the details of this scenario, let us determine the condition on the parameters Λ, f, m for the potential to exhibit local minima. In order to have local minima, the equation V = 0 must have non-vanishing solutions. Let us, for a moment, simplify by setting γ = 0. We then have: Graphically, it is clear that this equation has non-vanishing solutions only if (2.5) Here we have used γ = 0, but the equation remains parametrically valid even for γ = 0. Under this condition the potential has the form represented in figure 1. Practically however, as we have already remarked, we will focus only on the two lowest local minima, which are in general non-degenerate for γ = 0 (see figure 2). 3 We are now ready to move on to a detailed discussion of phase decomposition in our scenario. However, before doing so, a few comments about reheating are in order.

Reheating
Typically, after inflation the Universe undergoes a so-called reheating phase, where the energy density in the inflaton sector is transferred to standard model degrees of freedom, and possibly to dark sectors. Our analysis of the dynamics of the inflaton after inflation needs to take 2 The actual argument to show that wells at the bottom of the potential are more likely to host the inflaton requires more care. In section 4 we provide numerical examples that support this statement. 3 The choice γ = 0 may lead to stable domain walls, which are generically problematic. this into account. We therefore first focus on the following question: does reheating happen before or after the inflaton field is caught in one of the cosine wells?
In order to answer this question, we need to specify the interactions of the inflaton with matter and/or radiation. Here let us consider the example of a Planck-suppressed moduluslike coupling to a scalar field χ. The largest decay rate, barring the possibility of parametric resonances, is obtained if χ enjoys a coupling (1/M p )χ 2 φ, leading to This can arise if χ is a Higgs boson (see e.g. [59,60]). Decay rates to gauge bosons may also be of this type. The inflaton may also couple to scalar matter as (1/M p )φχ χ + c.c., but this leads to even smaller decay rates. The strategy is now as follows: we assume that the field is trapped in one of the cosine wells, and compare the perturbative decay rate (2.6) with the Hubble rate H 2 ∼ Λ 4 /M 2 p , as we are assuming that the field is oscillating in the last wells. For the same reason, we should take m 2 φ = V | min Λ 4 /f 2 in (2.6). As usual in cosmology, the condition for the decay to be efficient is Γ H. If we find Γ < H, then we will be consistent with our assumption that  V φ is as in (2.1) and V λ = 0. The initial condition is φ(t 0 ) = M p . Furthermore we have chosen κ = 60, f = 10 −2 M p . The scalar field oscillates around φ = 0 over a wide field range crossing several wells, before getting caught in one of the local minima at t ≈ 11/m. the field is first caught in one of the wells and decays only later. Now, The inequality (2.7) is easily satisfied in our setup, as in quadratic inflation m ∼ 10 −5 M p and f may be only slightly smaller than M p , while κ O(1). Therefore the field generically decays perturbatively only after getting caught in one of the cosine wells. In what follows, we will therefore not consider the decay of φ any longer.

Field oscillations and damping
We first focus on how the Hubble parameter changes after inflation. The energy density of the universe is a sum of three terms: ρ = V φ + T φ − V λ , where V λ is the energy density due to the cosmological constant. We absorb V λ in V φ and define V φ such that the global minimum has vanishing potential energy. The evolution of the Hubble parameter and of the scalar field φ is dictated by the Friedmann equation, together with the equation of motion of φ: Since the energy density is decreasing due to friction, it is clear that also H(t) will decrease with t. However, wheneverφ = 0, φ is undamped and the energy density is constant. In consequence H(t) is stationary as well. The typical behaviour of φ(t) and H(t) is shown in figure 3.
Let us now examine the energy density in the inflaton field in more detail. Its amplitude decreases after each oscillation due to Hubble friction. At some point the energy density ρ φ is comparable to the height of the last cosine wells. The field is then caught in one of the local minima, depending on the initial conditions. In the absence of spatial field inhomogeneities, the inflaton will populate only one of these two minima. This situation is shown in figure 4.   Figure 5. Same as figure 4, but now only the last two wells are shown. The width of the shaded bands represent the uncertainty of the scalar field energy density δρ, due to field fluctuations. The distance between the two bands corresponds to the energy loss ∆ρ due to friction.
The conclusions can radically change in the presence of field fluctuations δφ(t, x). In this case, the field may end up in one or the other minimum in different regions belonging to the same Hubble patch. The lines drawn in figure 4 become bands of a certain width, corresponding to the uncertainty δρ in the energy density induced by δφ (see figure 5). Now suppose that, as a consequence of friction, the energy density has decreased to a value close to the height of the barrier separating the two last minima in figure 4. During the next oscillation, the field will start rolling inside one of the two wells, say the one on the left of figure 5. At the end of the oscillation, since its energy density is smaller than the height of the barrier separating the two minima, the field is very likely to remain confined in the left well. However, due to the field fluctuations, there is a non-vanishing probability that the field actually reaches the other well on the right hand side and remains there in some regions of the Universe. If this is the case, at different point in the same Hubble patch the field lives in different minima. Therefore, the Universe decomposes into two phases. This is precisely the situation we are interested in. Quantitatively, let us compute the energy density lost during one oscillation, as a consequence of Hubble friction. In general this is not an easy task given the complicated shape of the potential in figure 4. However, it is greatly simplified by focusing only on the last cosine wells. In this case, the loss of energy during a half oscillation inside a single well can be estimated by a quadratic approximation of the potential: with M 2 = |V | min , i.e. the curvature of the potential (2.1) around the minimum of the well, where cos(φ/f ) ≈ −1. We obtain since we are interested in the regime κ > 1, we can take M 2 ≈ Λ 4 /f 2 . Immediately after inflation the Hubble rate evolves approximately as during matter domination, i.e.
so that the relative decrease in energy density in one half period ∆t ∼ M −1 /2 is given by where ∆ρ ≡ |ρ f − ρ i |. We now focus on the last two wells, such that ρ ∼ Λ 4 . Using Friedmann's equation, (2.14) We can now quantitatively discuss the probability of having a phase decomposition. To this end we have to compare the decrease in energy density due to friction ∆ρ and the uncertainty due to field fluctuations δρ. If we have δρ ∼ ∆ρ, the field will populate more than one vacuum with O(1) probability, as should be clear from figure 5. The term probability here refers to the exact choice of model parameters which, at that level of precision, appears arbitrary to the low-energy effective field theorist. The task of the next section is therefore to present two possible origins of the fluctuations δρ. These are respectively classical inflationary inhomogeneities and quantum uncertainties of the scalar field φ.

Fluctuations and phase decomposition
In this section we will analyse two sources of fluctuations of the inflaton field. In section 3.1 we focus on (classical) inflationary fluctuations. Those originate as sub-horizon size quantum fluctuations, but are then stretched to super-horizon size and become classical. After inflation, they re-enter the horizon and may lead to the phase decomposition described above. In section 3.2 we focus instead on the intrinsic quantum uncertainty which characterizes a quantum field at any given time. Independently of any inflationary pre-history of our field, this effect is present directly during the oscillatory stage and may also lead to phase decomposition. In fact, the estimates that we obtain in this section imply a small probability of phase decomposition. However, in section 4 we will see that the probability can be much larger because fluctuations may be enhanced after inflation.

Inflationary fluctuations
The evolution of inflationary fluctuations on sub-and super-horizon scales is well-known (see e.g. [61,62]). For us the only crucial point is that, once a certain inflationary mode re-enters the horizon it behaves like a dark matter fluctuation, i.e. it is a decaying oscillation (see [63,64] for a detailed study of scalar field fluctuations after inflation). From now on we focus on the amplitude of such an oscillation, which we denote by δφ inf k . The background is denoted by φ 0 and satisfies the equation of motion (2.2).
The initial conditions on δφ inf k are determined by matching with the power spectrum of the gauge-invariant curvature perturbation R. This quantity is conserved on superhorizon scales and is given by where the right hand side is evaluated at horizon exit, i.e. for k ∼ H. In the slow-roll regime the quantity in (3.1) coincides with the familiar expression = 1 To determine the probability of phase decomposition we need to understand how field fluctuations δφ inf k , once they re-enter the horizon, give rise to density fluctuations δρ inf k . More specifically, we wish to determine δρ inf k at time t Λ . This is the time at which the amplitude of the background φ 0 has decreased to values comparable to the width of the last wells, and it is given by Let us now consider a mode which exits the horizon just before the end of inflation, i.e. when H ∼ m. Such a mode re-enters the horizon shortly afterwards, at time t 0 ∼ m −1 < t Λ and has therefore momentum k(t 0 ) ∼ m. At t 0 we can determine the size of δφ inf k by matching with the curvature perturbation. However, the field fluctuation δφ inf k thus obtained will be out of phase with the oscillation of the background φ 0 . It is maximal whenφ 0 is maximal and vanishes when φ 0 is at a turning point. If this remained the case for the subsequent evolution until t Λ , the fluctuation could not give rise to a phase decomposition. This can only occur if we have a sizable fluctuation at a turning point of φ 0 .
However, we expect decoherence between δφ inf k and φ 0 after only a few oscillations. The reason is that φ 0 oscillates with frequency m while a mode with k ∼ m will oscillate with frequency √ k 2 + m 2 ∼ √ 2m. Thus, at some time after t 0 the field fluctuation δφ inf k will be an admixture of out-of-phase but also in-phase-oscillations w.r.t. to φ 0 . It is exactly the in-phase-oscillations which do not vanish at turning points and it is these fluctuations which give rise to density perturbations δρ inf k . In the following we will thus assume that once a mode enters the horizon, while initially out of phase with φ 0 , it will give an O(1) contribution to an in-phase oscillation with corresponding δρ inf k after only a few periods. 4 To estimate the size of fluctuation at time t Λ given a fluctuation at t 0 we need to take the expansion of the universe into account. In particular, the energy density scales as Thus, for a density fluctuation with k ∼ m there is a dilution in the time span between t 0 and t Λ .

JCAP11(2016)003
Let us now determine the probability of phase decomposition due to a mode with k ∼ m at the time t 0 of horizon re-entry. As argued before, the field fluctuation will quickly give rise to a density fluctuation. Instead of taking the intermediate step via field fluctuations, let us match the density fluctuations directly to the curvature fluctuations at re-entry: Now, using where the second equation follows from H ∼ m, ∼ 1 at the end of inflation, we obtain The above probability was derived for modes with k ∼ m at horizon re-entry. However, we are interested in the situations when the above probability is largest. Thus let us consider how this result is modified if we consider modes that re-enter the horizon later. Such modes have k < m at re-entry and they spend less time inside the horizon before the moment t Λ . Hence they would in principle give rise to a larger probability than (3.6). However, they cannot enter too late as they need to have enough time to decohere w.r.t. φ 0 . A more detailed analysis would be needed to determine how late a mode can enter the horizon and nevertheless give rise to a sizable density fluctuation. Thus, (3.6) should be seen as a reasonable estimate. We shall comment more on the size of this probability at the end of section 3.2.

Quantum fluctuations
There is another potentially relevant source of fluctuations of the field φ. This is the intrinsic uncertainty due to the quantum nature of our scalar. It can be simply written as In order to estimate the maximal effect of these fluctuations, let us consider the following setting: consider a scalar field φ with fluctuations δφ q k approaching the local maximum of some potential. This is basically as in figure 2, with the field approaching the maximum from the left side. We then ask the following question: what is the field distance from the local maximum which the background field has to reach such that its fluctuations can lift it over the potential barrier? The relevant energy scale around the maximum is O(M ), where M 2 is the curvature of the potential at the maximum. One can convince oneself that only fluctuations of the order δφ ∼ M are relevant for overcoming the barrier: modes with k M have an effective non-tachyonic mass and are insensitive to the instability. In contrast, modes with k M are sensitive to the tachyonic instability, but their fluctuations are smaller than those of modes with k ∼ M .
Alternatively, we can understand this point by considering tunnelling. Hence we take a homogeneous scalar field approaching the maximum and study the conditions under which quantum tunnelling to the other side of the barrier becomes efficient. In order to answer this question, let us use the following standard tunnelling formulae: the tunnelling rate is given by -10 -JCAP11(2016)003 e −S 0 , where S 0 is the action of critical bubble formation. When the thin wall approximation is applicable, this reads (see e.g. [65]) where σ ∼ M δφ 2 is the bubble wall tension and ∆V can be estimated with a quadratic approximation, ∆V ∼ M 2 δφ 2 . However, in most of the cases the thin-wall calculation is not appropriate. Nevertheless, we still expect S 0 ∼ δφ 2 /M 2 , with a different prefactor. 5 According to (3.8) the tunnelling rate is unsuppressed when δφ 10 −1 M , which is the same condition we found with the previous approach up to a numerical prefactor. Given the uncertainty in the prefactor, for the time being we use the parametric dependence δφ ∼ M .
Following these two arguments the uncertainty in the energy density induced by such field fluctuations is If quantum fluctuations induce an energy gain which is larger than the loss due friction one expects phase decomposition. Therefore, we compare (3.9) with (2.14), using also M 2 (3.10) Comparing (3.6) and (3.10) we conclude that the probability of a phase decomposition due to inflationary fluctuations is larger than the one due to quantum ones by at least a factor Let us now discuss the size of the probabilities (3.6), (3.10). In quadratic inflation, m ∼ 10 −5 M p . For κ ∼ O(10) and f ∼ 10 −13/5 M p the probability (3.6) is of the order of 0.1. This implies that a phase decomposition is rather likely for these values of parameters. For the same choices, the probability (3.10) is only slightly smaller, i.e. P q ∼ 10 −6/5 . However, further numerical suppression is expected in (3.10). Therefore we conclude that in the regime κ ∼ O(10), f 0.3 · 10 −2 M p phase decomposition is likely to happen as a consequence of inflationary fluctuations. If f grows above 10 −2 phase decomposition quickly becomes improbable. Whether such small values of f are natural depends on the details of the model leading to (2.1). We will comment more on the size of f at the end of section 4.
Further there exist arguments [50] based on the Weak Gravity Conjecture (WGC) [46] which constrain the size of modulations of the potential in axion monodromy inflation (for an earlier somewhat different perspective see [42]). It is thus important to check whether the region of parameter space considered in this work is consistent with these bounds. Given a domain wall with tension T and charge e the electric WGC demands T eM p . In our case we have T = Λ 2 f and e = 2πmf (see [50] for details). Using our definition As a result, while κ is bounded from above, our preferred parameter region for phase decomposition (κ ∼ O(10), f 10 −2 M p ) is consistent with the WGC. Interestingly, let us notice that the region of parameter space which is ruled out by the WGC is also constrained JCAP11(2016)003 by current observations, which require Λ 4 /(m 2 φ 2 ) (10 −3 − 10 −2 ) during inflation [14,66], where φ M p .
To close this section let us make two important remarks. First, note that phase decomposition can still occur even if the probability is small. In this case we only expect very few bubbles per Hubble patch. The second comment concerns once again the distinction between classical and quantum fluctuations. These two sources can also be distinguished based on the length scale R at which their effect is strongest. As we explained, classical inhomo- Since f < M p , the size R of the latter inhomogeneities is parametrically smaller than that of the inflationary ones.

Enhancement of fluctuations
Until now we have assumed that fluctuations, whether they are of classical or quantum nature, remain small during the evolution of the universe after inflation. The aim of this section is to discuss a possible enhancement of the fluctuations δφ k due to the functional form of the potential (2.1). Before going into details, let us summarise the main result. In this section we are mainly interested in fluctuations with k ∼ m at time t Λ . When the field oscillates at the bottom of the potential containing only the last few well, these fluctuations can be enhanced for certain values of f and κ. Crucially, these modes never exit the horizon during inflation, because at t < t Λ their wavelength is smaller than the Hubble radius. Therefore, these modes are never classicalised. Nevertheless, in this section we study their enhancement treating them as classical. We expect that our analysis will still capture the main effect. We provide a more detailed discussion at the end of this section.
A large growth of fluctuations can severely affect our conclusions. On the one hand, large fluctuations of the inflaton field may be desirable to some degree in our setup: the larger δφ k , the easier it is to cross the barrier between two local minima. Furthermore, if a mode with k m has a large amplitude, it may induce a phase decomposition independently of the modes with k ∼ m that we have studied in the previous section. If this is the case (3.6) underestimates the probability for phase decomposition.
On the other hand, large fluctuations with k m can lead to short-range violent dynamics rather than to the formation of well defined bubbles (which need a length scale 1/M ). As we will describe in more detail in section 5, this can negatively affect the strength of the gravitational wave signal related to our setting.
For these reasons, it is crucial to assess if any large growth of fluctuations occurs in our setting. As we have already mentioned, we focus on the enhancement of classical fluctuations and comment later on the applicability of the results for quantum modes. Such an analysis involves solving the coupled equations of motion for the background field and for the fluctuation in an expanding background and in the presence of a non-zero gravitational field. In appendix A.3 we present a first step towards this goal by providing the relevant equations of motion.
Here we will perform a somewhat different, simplified analysis, assuming that the gravitational field is negligible for the following reason. As we describe in appendix A, the addition of a gravitational field leads to a dark matter-like growth of fluctuations, which is negligible compared to the exponential growth that we are seeking in this section. Let us rewrite the equations of motion in a way that is more suitable for a numerical analysis. Namely, we where 'prime' now denotes a derivative w.r.t. t and where we have used a a = 2 3t during matter domination.
In the absence of friction, the background solution ϕ 0 is periodic and the equation of motion for δφ k is a Hill's equation. Solutions to such an equation exhibit a resonant behaviour for certain values of k [40,67]. This is similar to the resonances encountered in the context of preheating (for a review see [68]). However, the inclusion of a time-dependent background as well as friction, may affect the growth of the solution. Generically, one expects that modes with k m may experience exponential growth for certain values of f and κ. We investigated the behaviour of φ 0 and δφ numerically, for certain values of the parameters. In figure 6 the background φ 0 is plotted as a function of t for f = 10 −2 M p , κ = 60. For this parameter choice the field is caught near one of the local minima after only one oscillation.
For the same values of f and κ we study the evolution of the mode δφ 5m . In particular, we distinguish two regimes. Firstly, we focus on the behaviour of δφ 5m until the time at which the field is caught in one of the local minima. This is the interesting regime for the mechanism of phase decomposition that we have presented in section 2.3. As shown in figure 7, the mode does not grow in this time interval. Note that our equations are homogeneous in δφ k . Hence our numerical determination of the enhancement is not affected by the initial value of δφ k .
Secondly, we show the behaviour of the fluctuations after the field is caught around one of the local minima in figure 8. We observe that the mode grows after the background settles in one of the wells and decays again at late times. In particular, the field is caught around t 1 ≈ 8/m and by the time t 2 ≈ 13/m the size of its fluctuation δφ 5m has grown by roughly one order of magnitude with respect to its initial value δφ| initial ≈ 10 −5 M p . 6 Meanwhile, the amplitude of the background field φ keeps decreasing between t 1 and t 2 , as shown in figure 6: 6 This growth is nevertheless small compared to the growth by several orders of magnitude that we find below for the case where the field performs several oscillations before being caught in one of the minima.  in particular |φ(t 1 )| − |φ(t 2 )| ≈ 1.5f . The mild growth of δφ described in this particular example is therefore not likely to induce a late-time phase decomposition. Other modes may nevertheless experience a stronger growth after the field settles in one of the local minimum. However, as our previous example shows, a fluctuation δφ f is generically needed to overcome the decay of the background amplitude. This possibility clearly cannot be studied by means of linearised equations of motion for δφ. Therefore, in what follows we will focus purely on the growth of fluctuations before the field is caught in one of the cosine wells. We leave the interesting possibility of a late-time phase decomposition before reheating for future study.
Focusing on the first regime, the situation can be radically different from the previous example, depending on the values of the parameters appearing in the potential. In figure 9 we plot the background scalar field for f = M p /300, κ = 20. The field gets stuck in one of the cosine wells after six oscillations. As a consequence of the longer time that the field spends oscillating across several cosine wells, fluctuations can now grow significantly. In figure 10, we plot the logarithm of the absolute value of δφ k , again for k = 5m, f = M p /300, κ = 20. We see that the amplitude of this mode grows by three orders of magnitude before the background  field settles in one of the local minima. In fact, such strong growth takes us out of the regime of validity of the linearized equation of motion (as δφ k becomes comparable to f ).
Recently, the growth of fluctuations in a potential with cosine modulations was studied in [40]. The authors argue that, for a relatively small Hubble scale and neglecting gravity, the fluctuations for k ∼ m can grow as e N k , where N k is roughly given by Here F (κ) is a function of the order up to a few whose value depends on the initial amplitude of φ. In our case, since H k∼m ∼ Λ 2 /M p ∼ κ 1/2 m(f /M p ), we conclude that fluctuations should grow with exponent: where we have dropped any κ-dependence due to our ignorance regarding F (κ). Our numerical examples, which take Hubble expansion into account explicitly, confirm that for small f such an exponential growth does happen for most values of κ 10. In contrast, we do not observe enhancement if the value of f is chosen too large. Apart from this qualitative discussion, we are unfortunately unable to provide an analytical understanding of the dependence of the enhancement on the parameters κ and f . This is partly due to the fact that the phenomenon strongly depends on the precise minimum the field ultimately settles in. The latter question can only be addressed in a probabilistic approach, i.e. we can only say where the field is more likely to get trapped. Therefore, we do not have a precisely monotonic dependence of the enhancement in terms of f and κ.
In the absence of an analytical treatment, we performed a numerical search for enhanced fluctuations, focusing on modes with k m at the time when H ∼ O(Λ 2 /M p ). The results are reported in figure 11 in the form of a grid of points. Each point corresponds to a value of κ and f . We observe that in the region of interest fluctuations tend to be enhanced whenever f M p /200. Here, we define 'enhancement' as follows: we will refer to a mode to be enhanced if its original amplitude has grown by roughly two orders of magnitude before getting caught. The enhanced δφ k is comparable to f and we are therefore at the boundary between the linear and non-linear regime. Interestingly, this boundary corresponds to parameter values such that the probability (3.6) is of the order 10 −1 . However, note that the probability in (3.6) was determined for modes which will exhibit k < m at t Λ . While such fluctuations may experience growth, the generic expectation is that enhancement does not occur for modes with k < m.
Let us now briefly discuss enhancement of quantum fluctuations. This is in principle a complicated issue: we cannot use the classical equations of motions to analyse the behaviour of the quantum system. However, it is important to notice that, if quantum fluctuations are initially enhanced, they quickly become classical. Here by "classical" we mean that their occupation number becomes large, such that (4.2) can be used to study their evolution. Parametric resonance in quantum mechanics and quantum field theory has been studied analytically and numerically: the conclusion is that quantum modes do experience exponential growth (see e.g. [69,70]). We expect that the same effect will occur also in the system analysed here. Therefore, our study of classical fluctuations should extend, at least partially, to quantum fluctuations. If more enhancement occurs in the quantum case, then phase decomposition is even more likely. We leave a more detailed study of this effect for future work.

JCAP11(2016)003
Finally, let us summarise our findings concerning the probability of phase decomposition: 1. For f 0.5 · 10 −2 M p enhancement is generically not observed. The probabilities (3.6) and (3.10) are small, so that a phase decomposition is unlikely. Nevertheless, it is still possible that very few bubbles of tiny size containing the state of lower energy are nucleated. As we describe in the next section, observational signatures from such a situation may be quite strong. Classical inflationary fluctuations are the dominant cause of phase decomposition in this regime.
2. For f 0.5 · 10 −2 M p we have the following situation. On the one hand, fluctuations with k ∼ m at t = t Λ are generically enhanced. These modes are genuinely quantum modes, since they never exited the horizon. In this region, the enhancement may be just large enough to give rise to a probability of phase decomposition of order O(1). Furthermore, we observe numerically that, at fixed f and κ, modes with k m do not experience the same exponential growth. This will turn out to be a useful observation when examining the gravitational wave signal from the associated phase transition.
On the other hand, according to (3.6) and (3.10), classical and quantum modes with k m can lead to a phase decomposition, even if they are not enhanced. Therefore we conclude that a phase decomposition is very likely to be induced. Assuming that our analysis of enhanced classical fluctuations extends to the quantum ones, the dominant cause of phase decomposition are quantum modes with k ∼ m at t Λ .
3. For f 10 −2 M p fluctuations with k ∼ m are strongly enhanced. In this region phase decomposition is very likely to occur. However, it is hard to provide any description of such a highly non-linear regime. Classical and quantum fluctuations with k m are generically not enhanced, but would also lead to phase decomposition according to (3.6) and (3.10).
One more comment is in order before moving on to the phenomenological signatures of our setup. Phase decomposition happens generically for rather small axion decay constants. One may question whether such values of f are plausible in the spirit of axion monodromy. The answer depends very much on the framework in which monodromy is implemented. In a stringy setup it is a question of moduli stabilisation: e.g. in the Large Volume Scenario (LVS) [71,72] decay constants are generically suppressed by the volume of the compactification manifold and are therefore naturally small.

Gravitational radiation from phase transitions
In the previous sections, we have described a mechanism which can potentially lead to phase decomposition in the early universe after inflation. In this section, we will assume that such phase decomposition indeed occurs. In the presence of different populated vacua, bubbles containing the state of lowest energy can form and expand. Their collisions are a very interesting and well-known source of gravitational radiation. This has been studied in detail in the literature in various contexts and regimes (see [56] and references therein).
The aim of this section is twofold. First of all, we would like to elucidate the peculiarities of our setup concerning the energy released into gravitational waves during the collision of bubbles. Rather than focusing on precise calculations, we will give a qualitative discussion and provide formulae analogous to the more familiar case of bubbles colliding in a relativistic -17 -

JCAP11(2016)003
plasma. In this case, there are three possible sources of gravitational radiation: the collision of bubble walls, sound waves in the plasma and its turbulent motion. The second goal of this section is to give estimates of the relic density and frequency of the gravitational wave signal which can be obtained in our setup.

Gravitational waves from bubble collision
The focus of this subsection is the collision of bubbles and the possible shocks in the fluid surrounding them. These phenomena are usually studied in the so-called envelope approximation [53]. The latter consists in assuming that the energy liberated in gravitational waves resides only in the bubble walls before the collision. Furthermore, it is assumed that only the uncollided region of those walls contributes to the production of gravitational waves, i.e. the interacting region is neglected. Such an approximation has been initially applied to the case of vacuum-to-vacuum transitions [52] and later to collisions in a radiation bath.
In a thermal phase transition, the energy released into gravitational waves depends on four parameters. First of all, there is the time scale of the phase transition δ −1 or, equivalently, the initial separation between two bubbles d ∼ δ −1 . Secondly, there is the ratio η of the vacuum energy density released in the transition to that of the thermal bath, i.e.
where specifies that the quantity is evaluated at the time of completion of the phase transition. Thirdly, the efficiency factor λ characterizes the fraction of the energy density which is converted into the motion of the colliding walls. Finally, the bubble velocity v b is not necessarily luminal, as the walls have to first displace the fluid around them. The energy released into gravitational waves of peak frequency is then given by [54]: where ρ tot is the background energy density at completion of the phase transition. The parameters v b and λ are actually expected to be functions of η, in such a way that for In our case, bubbles collide before reheating, therefore there is no radiation bath around them. However, as we describe in appendix A, an oscillating scalar field corresponds to the presence of a matter fluid. Crucially, the time scale of the field oscillations is set by m −1 , and may be smaller than the time of collision. Therefore, oscillations of the scalar field cannot be generically neglected. Unfortunately, we do not have specific formulae for this case. Since we are interested only in an order of magnitude estimate for the spectrum of gravitational waves, it seems reasonable to extend (5.2) to our setup, with the obvious modification Furthermore, we shall hide our ignorance about the dependence of λ and v b on η by defining θ 0 = θλ 2 v 3 b and leave the determination of these parameters for future work. Therefore, we base our estimates on the following formula for the energy released in gravitational waves from the collision of bubbles and shocks in the matter fluid: where in our case ρ tot ∼ Λ 4 . In addition, the peak frequency of gravitational waves in the envelope approximation is given by ω peak σδ, (5.5) where σ O(0.1) should be fixed numerically and includes effects due to subluminal bubble walls velocity. The next task is to estimate δ and η. Let us start with the ratio H /δ. The energy density at the time of the phase transition corresponds to the height of the barrier separating the two minima, therefore H ∼ ρ 1/2 /M p ∼ Λ 2 /M p . We expect the typical frequency of the phase transition to be set by the momentum k of the spatial inhomogeneities of the field φ. The phase transition can be induced by any mode which is present at t = t Λ . The largest frequency that one can take is set by k ∼ M , as we have already discussed in section 3.2. This corresponds to a scenario where phase decomposition is likely. However, bubble collisions are most violent when the field makes it over the barrier separating the two minima only very rarely. In this case there are only few bubbles per Hubble patch. This latter scenario gives the strongest signal as it corresponds roughly to In order to understand how strong can the signal be in our setup, we assume (5.6) in what follows, but one should keep in mind that this is optimistic. The estimate of η is less straightforward, at least conceptually. If we adopt the envelope approximation, then we need to compute the vacuum energy density released in the phase transition. This is simply the difference between the energy density of the two minima in figure 12. Using a quadratic approximation, we find where ∆φ ∼ f is the approximate field separation between the two minima. The energy in the matter fluid ρ matter is roughly given by the height of the deepest well, i.e. ρ matter Λ 4 . This is because at completion of the phase transition the oscillations of the scalar field span almost the whole well. Therefore, in the envelope approximation we obtain As we have mentioned, deviations from this simple picture may arise in our case. On the one hand, a certain fraction of the energy of the walls might for example be dissipated into the matter fluid. In this case, only a fraction of would lead to production of gravitational waves. This effect might be captured by the efficiency prefactor λ. 7 On the other hand, the energy released into the fluid while the bubbles expand and collide might also contribute to the production of gravitational waves. Namely, this energy might be converted into bulk motion of the fluid. In this case, the energy released in gravitational radiation should be larger than , and could possibly be as large as Λ 4 . This effect is captured by studying the fluid as a source of gravitational waves. We comment very briefly on this topic in the next subsection. Figure 12. Two-well potential. In the picture is the energy difference between the two minima, while Λ 4 is the value of the potential at the local maximum.

Gravitational waves from the matter fluid
In analogy with the case of radiation there are at least two effects which can further contribute to the total energy released in gravitational radiation during the phase transition. Here we just provide the formulae given in [56] for the thermal case, keeping in mind that they may not straightforwardly extend to our setup: 1. Sound waves in the fluid : this arises because a certain fraction λ v of the energy of the walls is converted after the collision into motion of the fluid (and is only later dissipated). In the case of radiation this gives a contribution The prefactor θ sw is expected to be smaller than θ in (5.2).
2. Turbulence in the fluid : one expects further contributions as a certain fraction λ turb of the energy of the walls is converted into turbulence. In the case of radiation one obtains: The prefactor is expected to be larger than θ in (5.2). Note, that for these two mechanism the dependence on H/δ is only linear.
In certain regimes, these two effects may be larger than the one due to bubble collisions and shocks in the fluid. However, they are not fully understood, even in the case of radiation. Therefore, in the next subsection we will neglect them, and obtain only a lower bound on the relic abundance of gravitational waves. This should still be useful to understand the approximate size and frequency of the signal. Nevertheless, the reader should keep in mind that there are other possible contributions even beyond the ones mentioned in this subsection (see e.g. [73] for recent progress).

Frequency and signal strength of gravitational waves
In order to compute the relic abundance and frequency of gravitational waves emitted during the phase transition, we need to know the equation of state of the background energy density -20 -

JCAP11(2016)003
from the end of the phase transition to today. Assuming standard evolution after reheating the behaviour of the scale factor until today is essentially fixed. 8 Furthermore, the inflaton field generically behaves as non-relativistic matter after inflation. It remains to be addressed whether deviations from the equation of state w ≈ 0 might occur immediately after the phase transition, before reheating. Due to the very large release of energy during the collision of bubble walls, it is conceivable that the fluid initially behaves relativistically. This would correspond to an early phase of radiation domination, i.e. w = 1/3, in a similar fashion to some preheating scenarios. Eventually, the fluid cools down and its non-relativistic behaviour is restored. Depending on the reheating temperature, this may or may not happen before the inflaton decays. If the system were in a thermal ensemble, the fluid would behave non-relativistically after T ∼ m φ ∼ M . For the time being, we allow for a general equation of state parameter w after the phase transition and before T ∼ M .
Therefore the background energy density at reheating is given by where from now on the subscript N R denotes that a certain quantity is evaluated at the time when the fluid becomes non-relativistic. Let us define the prefactors The prefactor ν w quantifies the duration of a period of matter domination before reheating, while ν nr quantifies the duration of an early epoch of matter domination. Obviously 0 < ν w , ν nr ≤ 1. The energy density in gravitational waves scales as a −4 . According to (5.4), at reheating we have 14) where in the last step we have used (5.11). The relic energy density of gravitational waves is then where t 0 is the current age of the Universe. The ratio a RH /a 0 can be determined by imposing entropy conservation. Furthermore, ρ RH can be computed using the standard formula for the energy density of radiation It is characterized by the effective number of degrees of freedom.

JCAP11(2016)003
Finally, the density parameter today Ω GW (t 0 ) ≡ today is given by where h ≡ H 0 /(100 km · s −1 · Mpc −1 ). Let us now estimate the peak frequency of the emitted radiation. For this quantity the only relevant parameter is δ. Frequencies scale as ∼ a −1 , therefore we have: By combining (5.18), (5.5), H 2 ∼ ρ/M 2 p and (5.16) one obtains We have therefore determined the relevant parameters of the emitted radiation as function of δ and η. Now we can plug in (5.6) and (5.8), to obtain final formulae. Then we have: In the envelope approximation, it is also possible to compute the full spectrum of the gravitational radiation emitted from the collision of the bubble walls. This reads [56]: with S env (ω) = 3.8(ω/ω 0 ) 2.8 1 + 2.8(ω/ω 0 ) 3.8 .
In order to estimate the maximal possible size of our signal, let us now specify to the case in which the energy present after the bubble collisions is converted into radiation. This corresponds to setting w = 1/3 in (5.22). At T ∼ M the inflaton goes back to a nonrelativistic behaviour. The prefactors ν w and ν nr can now be explicitly computed, using . from (5.11), the background energy density scales as matter from completion of the phase transition until reheating.
In figure 13 we plot the spectrum (5.22) for three different choices of parameters, using w = 1/3. We fix m 10 −5 M p as required by observations. We have chosen parameters in such a way as to maximize the overlap with sensitivity regions of current and future space-and ground-based detectors, which are bounded by dashed lines in the plot. We have also fixed θ 0 = 10 −2 , σ = 10 −1 . Our plot provides examples of the wide range of frequencies that can be obtained in our setting, simply varying the axion decay constant, the reheating temperature and the number of local minima. Interestingly, for reasonable choices of parameters the signals are in the reach of future detectors.

Conclusions
In this paper we investigated the production of gravitational waves from post-inflationary dynamics in models of Axion Monodromy inflation. We expect such phenomena to also occur for generic axionic fields with potentials with sizable modulations, albeit the numerical results may differ significantly in these cases.
The main observation is that in models of axion monodromy inflation, the inflaton potential consists of a monotonic polynomial with superimposed cosine-modulations. While these modulations have to be small to allow for successful inflation, they tend to dominate near the bottom of the potential. In fact, these cosine 'wiggles' can be large enough such that the potential exhibits a series of local minima.
After inflation ends, the inflaton is exploring this wiggly part of the potential. As a consequence of Hubble friction, the rolling axion can get stuck in one of the wells. This -23 -may happen well before the inflaton reheats the standard model degrees of freedom. Since the energy density in the axion field decreases together with H, the field is more likely to get caught in one of the last minima. We therefore focused on a two-well setting, which is obtained by "zooming" into the full monodromy potential.
Taking fluctuations of the inflaton field into account, the inflaton field does not necessarily get caught in one unique local minimum in the entire Hubble patch. If field fluctuations are sufficiently large, a phase decomposition occurs such that at least two different vacua are populated after inflation. The probability of this occurring is given by P ∼ δρ/∆ρ, where δρ is the uncertainty in the axion energy density induced by the field fluctuations and ∆ρ is the frictional loss of energy in one oscillation.
We can distinguish two sources of field fluctuations which may lead to a phase decomposition. Firstly, there are the classical inhomogeneities naturally inherited from inflation. Secondly, there are the intrinsic quantum uncertainties characterising any quantum field. These two sources are essentially indistinguishable at very early and late times. However they are in principle of different size in the intermediate regime that we are interested in. In particular, we have found that inflationary fluctuations are more likely to induce a phase decomposition: the probability that they do so is P ∼ κ −1/3 (m/M p )(M p /f ) 5/3 . Here m is the mass of the inflaton-axion, f the axion decay constant which defines the axion periodicity 2πf and κ/π roughly counts the number of local minima of the potential. Therefore, we observe that a phase decomposition is likely for κ ∼ O(10) and f 0.3 · 10 −2 M p . The probability that quantum fluctuations induce a phase decomposition is smaller by a factor (f /(M p κ)) −4/3 m/M p .
Furthermore, due to the oscillatory term in the axion potential, fluctuations can experience exponential growth. This effect arises generically for modes with k ∼ m at the time when the field is rolling over the last wells. These fluctuations never exited the horizon, therefore they are effectively quantum modes. We extended the study of enhancement in [40] to the case with varying H, but still neglecting the gravitational field in the equations of motion. Numerically and using a classical approximation, we observe the existence of a region of parameter space where a phase decomposition is likely to occur. This happens roughly in the same regime where inflationary classical fluctuations are also likely to induce a phase decomposition. However, given the exponential enhancement of quantum modes, the latter are more likely to be the dominant cause of phase decomposition. Enhanced quantum modes quickly become classical, therefore we expect our main results to hold even after a more detailed analysis, which we leave for future work. For larger values of f modes are not enhanced. Phase decomposition, although unlikely, might still occur as a result of quantum or classical fluctuations. For smaller decay constants fluctuations are very strongly enhanced. Phase decomposition occurs but it is hard to understand the physics in such a highly non-linear regime.
If a phase decomposition occurs, bubbles containing minima of lowest energy expand. Collisions of these bubbles source gravitational waves. We estimated the energy density and frequency of the emitted radiation in terms of the axion parameters, in the envelope approximation. Furthermore, we note that the matter fluid associated to the oscillating inflaton may also radiate gravitational waves. This is similar to the case of a thermal phase transition. The spectrum of the emitted radiation can peak in a wide range of frequencies (from mHz to GHz), depending on the reheating temperature and on the time of the phase transition. In this sense, our source is similar to other post-inflationary phenomena, such as preheating and cosmic strings. However, it is interesting to observe that the frequency may be lowered -24 -JCAP11(2016)003 in our case since a matter dominated phase can follow the phase transition. The spectrum is at least partially in the ballpark of future space-and ground-based detectors. Thus, we can hope that axion monodromy may one day be investigated by means of gravitational wave astronomy.

A Scalar field fluctuations after inflation
In this appendix, we provide a detailed analysis of the evolution of scalar field fluctuations after inflation. We begin by deriving the Klein-Gordon equations for a scalar field and its fluctuations with potential (2.1), including the gravitational field. We then focus on the simple case of a purely quadratic potential. This gives us the opportunity to review why scalar field fluctuations behave like dark matter perturbations. We then provide equations to study the fluctuations in the full potential containing the 'wiggles'.

A.1 Equations of motion
Let us begin with the equations of motion of a scalar field in the post-inflationary universe. The starting point is the Klein-Gordon equation: The metric appearing in (A.1) is the perturbed FRW metric (here we follow [61]): where B i = ∂ i B +B i and the hat denotes divergenceless vectors and traceless tensors. In particular, we can consistently focus only on scalar modes, since the vectors can be gauged away and the tensors are not sourced by δφ. We therefore have and (A.2) are not gauge invariant. However, the following quantities are gauge invariant: In what follows we will perform our computation in the Newtonian gauge, defined by:

JCAP11(2016)003
In the latter, the perturbed metric reads: The components Φ and Ψ are related by the perturbed Einstein equations. In the absence of off-diagonal terms in the spatial components of the perturbed stress-energy tensor the Einstein equations impose Φ = Ψ. With this constraint, the metric (A.5) provides the newtonian limit of general relativity.
We are now in a position to write down the Klein-Gordon equation (A.1), expanding the scalar field as φ = φ 0 + δφ and keeping only the leading order terms in the perturbed quantities δφ, Φ. Then the background φ 0 obeys: while the inhomogeneous δφ(t, x) satisfies: Furthermore, there are three Einstein equations relating the gravitational field Φ to the fluctuation δφ:Φ where δρ =φ 0δ φ + V (φ 0 ) is the perturbed energy density. The Einstein equations (A.8) and (A.9) lead to a Poisson equation with a gauge invariant energy density on the right hand side. In total, we have five equations for three quantities and therefore it is enough to consider only one of the Einstein equations.

A.2 Scalar field as dark matter
In order to understand why perturbations of a scalar field behave like dark matter fluctuations, we now focus on the case of a purely quadratic potential V (φ) = 1 2 m 2 φ 2 . In this case it can be shown that fluctuations obeying (A.7) behave like an ideal pressureless fluid [63,75]. The strategy is as follows [75]: by using a WKB ansatz for φ 0 and δφ one can define three fluid quantities for a scalar field: energy density, velocity potential and pressure. One then performs an expansion of (A.7) in powers of (m −1 ). The leading order in this expansion corresponds to the Euler equation for a perfect pressureless fluid. The subleading order gives the continuity equation. This establishes a dictionary between scalar field and fluid quantities. More precisely, the correspondence involves field quantities which are averaged over one period T ∼ m −1 . One can then use the standard fluid approach to prove that scalar field energy fluctuations grow like a on subhorizon scales.
We now briefly review the approach of [75]: the aim is to provide equations which can be extended to include non-linearities. The following analysis is valid in the regime m H, k/a H, k 2 /(a 2 H) m. The first two conditions are satisfied in our setup.

JCAP11(2016)003
However, this is not necessarily the case for the third one. Indeed, as explained in section 3, we are mainly interested in modes with k m at horizon re-entry, as these are most efficiently leading to a phase decomposition. However, due to the expansion of the universe k 2 /(a 2 H) decreases as a −1/2 . Therefore, even if the condition k 2 /(a 2 H) m is originally not satisfied, it rapidly becomes fulfilled.
The starting point of the analysis is a WKB ansatz for the background and for the fluctuations (see also [76,77] for related discussions): The crucial point is that u, w, B and A vary on time scales which are much longer than m As expected, we see that the effective pressure of the fluctuations vanishes at leading order. A perfect pressureless newtonian fluid in an expanding background is subject to the familiar continuity and Euler equations (see e.g. [61]): where δ ≡ δρ/ρ and v i ≡ ∂ i v. According to (A. 18), v ∼ B. We see therefore that (A.13) closely resembles the Euler equation (A.21). In fact, it can be shown that (A.13), rewritten in terms of v, is equivalent to the relativistic generalisation of (A.21).

JCAP11(2016)003 B Tuned small field inflation
In this appendix we would like to briefly address an issue related to our setup. In principle, while oscillating along the full potential (2.1), the inflaton may come close to a local maximum, without being able to actually reach it as a consequence of friction. In this case, the field motion may satisfy the slow roll conditions, therefore leading to a further inflationary phase. If the local maximum is flat enough, we might be in a scenario which closely resemble that of inflation at an inflection point (see [78] and references therein, see also [16,20,21,79] for related previous work). If inflation is viable in this setup, then current observational constraints on m may be relaxed. Indeed in this case the axion mass would not necessarily be proportional to the inflationary power spectrum. This may turn out to be particularly relevant for the signatures that we described in this paper, as the probability of having a phase decomposition is suppressed by the smallness of m. Let us also remark that the potentials that we will consider in this appendix may be relevant also for non-monodromic scenarios with more than one axion, such as alignment and winding models [43][44][45]. Let us then start by assuming that the inflaton slowly approaches a stationary inflection point φ * of the potential. We first estimate how many e-foldings would be generated as the field comes close to φ * . In this appendix we work in Planck units, i.e. we set M p ≡ 1.
Around an inflection point, the potential can be approximated as The slow roll parameters are therefore given by The number of e-foldings can be estimated as: The latter equation may represent a conflict with observations: indeed in this model is very small, i.e. 10 −2 , therefore we need η ∼ 10 −2 to satisfy the observational constraint on n s − 1 = 2 − 6η 3 · 10 −2 . From (B.3), we obtain N ∼ 10 2 . Therefore it seems that in this simple case we obtain too many e-foldings.
However it is very unlikely that the field gets caught around such an inflection point. Indeed, the latter is the point after which there are no more local minima of the potential (2.1), and therefore sits far from the minimum of the parabola. As we have already remarked, the field is more likely to get caught in one of the minima at the bottom of the potential. Therefore, we now consider a more generic stationary point and modify the potential (B.1) by including also a small quadratic contribution.
We focus on the following expansion around a stationary point φ 0 :

JCAP11(2016)003
Here we take α and β to be positive without loss of generality. We want to assess the feasibility of such an inflationary potential (see [79] for the case in which a linear term is included, instead of a quadratic term). We take α β, and focus on small field ranges, so that |φ − φ 0 | < 1 all along the inflationary trajectory.
The slow roll parameters are: The number of efoldings is given by: where φ is the field value at the beginning of inflation. It is convenient to express η in terms of N , by inverting (B.7): where in the last step we have kept only the first order in α. Similarly, we can express in terms of η , then in terms of N by means of (B.8): therefore, under the assumption α β, η e = −1 determines the end of inflation. Notice also that the ratio α/β sets precisely the inflationary range. Now we can write down the spectral index and the tensor-to-scalar ratio in terms of α, β, N : n s = 1 − 6 + 2η 1 − 1 2 + N 2 + 3 4β 2 (2 + N ) 3 − α 1 + 2 3β 2 (2 + N ) 3 (B.11) r = 16 4 β 2 (2 + N ) 3 α + 1 2(2 + N ) . (B.12) We are therefore ready to study the parameter space of the model described by (B.5). First, we fix N = 50 and plot the constraints on the parameters α and β obtained by imposing the observed values for n s and r: n s = 0.9652 ± 0.0047, r < 0.10 at 95% C.L. [1]. We then repeat the analysis for N = 60.
We present the results in figure 14 and 15. We see that a model based on (B.5) is still a viable candidate to explain the CMB observables.