Exploring ALPs beyond the canonical

Axion-like particles (ALPs) are interesting dark matter candidates both from the theoretical as well as from the experimental perspective. Usually they are motivated as pseudo-Nambu-Goldstone bosons. In this case one of their most important features is that their coupling to other particles is suppressed by a large scale, the vacuum expectation value of the field breaking the symmetry that gives rise to them. This naturally endows them with very weak interactions but also restricts the maximal field value and therefore the regions where sufficient dark matter is produced. In this paper we investigate deviations from this simplest setup, where the potential and interactions are as expected for a pseudo-Nambu-Goldstone boson, but the kinetic term has singularities. This leads to a significantly increased area in parameter space where such particles can be dark matter and can be probed by current and near future experiments. We discuss cosmological limits and in the course of this give a simple derivation of a formula for isocurvature fluctuations in models with general anharmonic potentials. As an application of this formula we give an update of the isocurvature constraints for QCD axion dark matter models, using the most recent results for the QCD topological susceptibility and the newest Planck data.


Introduction
Axions and axion-like particles (ALPs) are a prediction of some of the best-motivated beyond the standard model physics scenarios (see, e.g. [1][2][3] for reviews). Many of their properties are determined by two quantities: the mass, m and the so-called decay constant, f a . An important feature that all these particles share is that they enjoy a shift symmetry, a discrete version of which is preserved at the quantum level. The existence of this symmetry protects their potential from quantum corrections that could otherwise be very large. In the framework of quantum field theory, such particles arise as pseudo Nambu-Goldstone bosons of approximate global chiral symmetries [4][5][6][7][8][9]. In other setups such as supergravity or string theory, particles with similar properties appear in the spectrum. For instance, ALPs are a general consequence of the compactification of extra dimensions and string theory [10][11][12][13][14][15][16]. In that context, there can be dozens of such particles whose potentials, kinetic terms and interactions may contain a large number of free parameters. In an attempt to accommodate all these similar particle candidates, we will talk about ALPs in the general sense of a light (pseudo-)scalar particle, and we will reserve the term "axion" to refer to ALPs that couple to the gluon field strength tensor through the QCD topological term and solve the strong CP problem.
Axion-like particles are excellent candidates to account for some or all the dark matter that we observe in the universe [17][18][19][20]. Cosmological and astrophysical observations tell us that dark matter particles should be weakly interacting, stable at cosmological scales and cold. ALPs can naturally fulfil all these requirements. First, the discrete shift symmetry constrains their possible couplings to other fields, and those that are allowed are typically suppressed by f a , which can be a large energy scale. This fact, together with their small mass which limits the possible number and type of decay products as well as the phase space, makes them extremely stable. Naively, the fact that they are very light might seem to contradict the requirement that the ALP dark matter population should be cold. However, it is easy to see that this is not necessarily the case. Because of their feeble interactions with other particles, ALPs are not produced thermally, but rather by the socalled misalignment mechanism, which yields a very non-relativistic population of ALPs that behave as cold dark matter [17,19,[21][22][23].
All in all, ALPs and axions are well motivated dark matter candidates, but their possible mass and decay constant span many orders of magnitude thereby providing a significant challenge for experimental tests. Fortunately, their properties, in particular their low mass, also provides for new opportunities for experimental searches and theoretical arguments that can be used to probe their parameter space (see [24] for a recent review).
Experimental tests are usually dependent on the coupling to Standard Model particles. One example is a coupling to two gluons, This coupling also induces a coupling to a nucleon electric dipole moment (EDM), that is particularly important for searches when φ is dark matter 1 . The coupling constants are related via [26,27] g d ≈ 2.4 · 10 −16 f a e · cm ≈ 3.4 · 10 3 GeV −2 GeV f a .  On the horizontal axis we plot the mass of the ALP, while the vertical axis gives the decay constant f a on the right and the effective coupling to nucleons g d ∝ f −1 a on the left. Canonical ALP models with a constant mass can only generate enough dark matter via the misalignment mechanism in the yellow and grey shaded areas. Accounting for the anharmonicities of the potential and allowing for a fine-tuned initial condition, this region can be enlarged to also include the orange band (we take the lowest viable Hubble scale of inflation, H I ∼ 4.5 · 10 −23 GeV). Note that the QCD axion models are restricted to lie on the magenta line. Taking the interaction to be given by Eq. (1), the region to the left of the QCD axion line is disfavoured by the unavoidable (temperature dependent) contribution to the mass from QCD effects [28] (see also §5). This region is shown in light grey. The dark blue region is excluded by the supernova limits estimated in [29]. Shaded in brown is the area where experiments looking for a static nuclear electric dipole moment (nEDM, see [30]) would have found the oscillating one, while the dotted lines represent sensitivity estimates for future oscillating EDM experiments [31,32]. In the dark green region "BBN" ALPs coupled to QCD are inconsistent with the production of the observed abundance of light elements during Big Bang Nucleosynthesis [28]. The violet and dark red lines dubbed "Earth" and "Sun" correspond to constraints from the ALP field being sourced by dense astrophysical objects [33]. The dark grey area is disfavoured by the observation of quickly rotating stellar black holes which would have been spun down in a superradiant process (from [34]). The area above the dashed black lines, plotted for different values of the Hubble scale of inflaton H I , is disfavoured due to the generation of too much power in isocurvature perturbations at the scales probed by the Planck satellite [35] (see more details in §4). Finally, f a is (softly) bounded from above by the requirement that it does not exceed the Planck scale. Figure 1 summarises the constraints that can be cast on the canonical ALP dark matter scenario from these interactions with the visible sector. In addition we show limits that arise from unavoidable gravitational interactions.
Unfortunately, some of the theoretically favoured existing models require high decay constants for the ALPs to be able to account for all the dark matter energy density that we observe in our Universe. This means that some of the better motivated combinations of (m, f a ) are not in the best position to be tested, be it through gravitational interactions or through couplings to gluons and nucleons or photons. It is therefore timely to search for models that can accommodate low enough values of f a that can be in reach of these searches, while still being able to produce the required dark matter abundance. One option is to enlarge the field range by a monodromy [36][37][38] as done in [39].
In this paper we pursue the same goal by employing a non-standard kinetic term for the ALP field. This is a possibility that has been exploited in the literature [40,41] in the context of inflationary models (though not so much for axion inflation), but to our knowledge such a study has not been performed for dark matter models. As we will see, a very rich phenomenology arises when this possibility is allowed. Of special interest is that this scenario will indeed be able to populate regions of the parameter space that can be tested in the near future, either with astrophysical observations or experimental searches. Focusing on the coupling to nucleons, the main motivation for us in this respect is threefold. First, as was already argued, we want to explore the possibility of building an ALP dark matter model with a larger such coupling. Second, we ask ourselves if these models could lie on the region of parameter space to the left of the QCD axion band in Figure 1. Finally and concerning the Big Bang Nucleosynthesis bound that seems to restrict this area of parameter space, we would like to test its robustness constraining such ALP models.
In this work we study the viability of ALPs with a non-canonical kinetic term as dark matter candidates from a purely phenomenological perspective. Let us nevertheless briefly mention some of the mechanisms that can give rise to this scenario. For instance, a nonminimal coupling of the ALP field to gravity in the so-called Jordan frame induces a noncanonical kinetic function in the usual Einstein frame. In the context of supergravity, an explicit breaking of the shift symmetry in the Kähler potential also results in non-standard kinetic terms for the ALP. Finally, in the context of compactifications, string theory a priori contains all the necessary ingredients to generate axions with non-canonical kinetic terms, caused, for example, by back-reaction effects. However, no explicit construction of the models that we consider in this work exists as of today, and this task is beyond the scope of this paper. We leave the study of the possibility of embedding this phenomenological study in a more complete framework for future work. This paper is structured as follows: in §2 we discuss the effects of non-canonical kinetic terms and set up our explicit case study. In §3 we study how this modified kinetic terms affects the cosmological evolution of the ALP field, and in §4 we analyse the isocurvature perturbations predicted in this setup. In §5 we discuss the impact of allowing for a coupling to QCD in this scenario, and conclude in §6.
Before getting started on the details we note that, although in this paper we focus mainly on the example of gluon interactions, most of our discussion is completely general and can be applied to any other coupling. Moreover, while the structure of interactions that we consider is inspired by that of pseudo-Nambu-Goldstone bosons, the essential qualitative features should also apply in the case of more general scalars and only depends on the singularities of the non-canonical kinetic terms.

Non-canonical kinetic terms
In this section we examine the effect that a non-standard kinetic term can have on the dynamics of the ALP field. Let us start with the Lagrangian where we have allowed for a general real scalar (and positive definite) function of φ, K 2 (φ), to scale the kinetic term and thus render it not canonically normalised. For definiteness, we will work with the usual periodic potential for ALP fields, We now proceed by performing a field redefinition to obtain the canonically normalised field. The formal solution is to define and thus the Lagrangian for ϕ is Being canonically normalised, ϕ is the physical (propagating) field. Let us see what kind of functions K result in ϕ being a viable dark matter candidate. The first condition is that ϕ behaves like cold dark matter in the late universe. This requires that it oscillates harmonically at late times (see, e.g. [17]). Accordingly the kinetic term should not modify the dynamics close to the origin. This is automatic if the kinetic term approaches a non-vanishing constant value close to the origin, As indicated in the equation, this constant can be chosen to be equal to 1 by a suitable choice of normalisation. So why should we now choose a non-trivial function for K and what shall we choose? As already mentioned in the introduction, we would like to find a model with larger couplings, i.e. smaller f a , that still gives a sufficient dark matter density. Roughly speaking the problem of obtaining a sufficient energy density can be understood as follows. For the potential Eq. (5) the maximal initial energy density is given by Λ 4 . This is linked to the mass m of the particle via Λ 4 = m 2 f 2 a . If f a is too small the initial and in consequence the final energy density is too small to make up all of the dark matter.
One way to avoid this problem would be to break the periodicity of the potential (5) such that the potential continues to grow for large field values, e.g. by exploiting a monodromy [39].
Here we will explore a different strategy. As long as the Hubble constant is sufficiently large the evolution of the field is frozen and the energy density is approximately constant. As discussed below the evolution and consequently the dilution of the energy starts when H 2 ∼ |V (ϕ)|. Hence, we can increase the energy density today by choosing the kinetic function K such that the potential becomes very flat for large field values 2 . A cartoon of this is shown in Figure 2. Using we see that this can be easily achieved if K has a singularity at some field value φ = a, This singular structure has an additional advantage: The non-canonically normalised field φ will never exceed φ = a during its evolution. Limits such as the one discussed in §5.4 arising from BBN that are based on a sizeable field value at some earlier epoch can thus be avoided if a is sufficiently small.
A simple function that satisfies the above requirements while keeping the periodic properties intact is, While this choice might seem rather arbitrary at first, there are some arguments that make it more general than it seems. The approach for obtaining a flattened potential for a scalar via a non-canonical kinetic term has been widely used in the context of inflationary cosmology [40,41]. Indeed over the last years, α-attractor models [44,45] have attracted special attention. In this context, [46] showed that the determining property of this class of models is the existence of a pole in the kinetic term. More precisely, it is the order and the residue of the pole that play a key role, and not so much the precise functional form of the kinetic function. We can therefore be confident that our results will not depend much on the specific choice of K. Similarly to [46], here we focus on the case of a second-order pole. As we mentioned before, this case is better motivated and may arise, for instance, as a consequence of a non-minimal coupling to gravity. Nevertheless, we check in Appendix §A that our main conclusions remain unchanged if we allow for higher-order poles.
Also, recall that the shift symmetry φ → φ + const. of the ALP field is what protects its mass from large corrections. It thus seems sensible to preserve or only slightly break this symmetry. Indeed, by our choice of potential Eq. (5), we are assuming that a small explicit breaking is present. This breaking typically occurs at the nonperturbative level [47,48] and crucially preserves the discrete shift symmetry φ → φ + 2kπf a , which allows us to retain a sufficient level of protection against quantum corrections. We would like the kinetic term to preserve, at least, this discrete shift symmetry, which requires that K(φ) is a periodic function of φ/f a . These arguments quickly lead us to Eq. (11). Once again, we stress that the fact that we are writing a specific kinetic term should not be understood as a construction of a complete model, bur rather as a benchmark for our phenomenological study.
The transformation to the canonically normalised field is given by We should note that the poles of K(φ) are located at φ/f a = π/(2N ). This means that, when doing the field redefinition (12), we are restricting the field space to φ/f a ∈ (− π 2N , π 2N ). As already mentioned above this will become important when discussing the limits arising from a gluon coupling in §5. In principle there exist a total of N different where the field could be trapped. However, the only one which has a minimum in the potential is the one closest to the origin. In other branches, the field would slow-roll towards infinity 3 , making them unappealing for the phenomenologically purposes that the we have in mind. For this reason, we focus on the phenomenologically viable region around zero.
Using the field redefinition (12) the Lagrangian for the canonically normalised field is given by By expanding about the origin, it can be checked that we indeed recover the quadratic behaviour for small field values. The potential is plotted in Figure 3 for different values of N . It indeed looks quite similar to what we imagined in Figure 2.
What about the equations of motion? Let us assume that we have a homogeneous and isotropic field, φ = φ(t) and consequently ϕ = ϕ(t). The Klein-Gordon equation for a homogeneous and isotropic field in an expanding spacetime is where H is the Hubble expansion parameter. For convenience we introduce the dimensionless field variable, in analogy to how the θ angle relates to the original axion field. Thus, we will be expressing the field value in terms of f a units. The equation of motion can then be written as where we define which corresponds to the second derivative of the physical field around the minimum at ψ = ϕ = φ = 0. m is the physical mass of the dark matter particles. Figure 3: Potential for the canonically normalised field, plotted for various values of N . Note that the potential is quadratic for small field value but flattens away from the origin.

Cosmological evolution and dark matter production
The goal of this section is to find an estimate for the dark matter density in the model defined above and compare it with the observed abundance. The energy density of the field depends on the parameters (f a , m, N ), as well as the initial conditions for the field and its cosmological evolution. For this purpose it is useful to briefly recall the misalignment mechanism [21][22][23], which gives us the basic idea of how our field evolves in a cosmological setup.

The misalignment mechanism
Here we briefly summarise how a misaligned light scalar field evolves in an expanding spacetime, closely following the description in [17]. Let us consider the simplified case of a real scalar field with Lagrangian Note that our final goal is not the harmonic case but a more complicated potential with strong anharmonicities. However, solving this simplified equation will give us helpful insights to tackle the anharmonic potential. In a homogeneous setting, the equation of motion This is the equation of a damped harmonic oscillator. There are two distinct regimes in the evolution of φ. First, at very early times when 3H m φ , the oscillator is overdamped and so the solution isφ = 0, and the field is stuck at its initial value. At a later time t 1 such that 3H(t 1 ) = m φ , the damping has decreased enough so that the field can start to oscillate. The equation of motion for the oscillating regime can then be solved using the WKB approximation: where a(t) is the scale factor. We see that the energy density, which is proportional to the amplitude of the oscillations squared, dilutes with expansion as a −3 . This means that the oscillating field behaves like pressureless matter for all processes mediated by gravitation.
In this simplified setup, the energy density in the axion field today is where is a smooth function (cf. [17]) that varies from 1 to ∼ 0.3 when T 1 ∈ (T 0 , 200 GeV). The last result assumes that the field starts oscillating during radiation domination and that the comoving entropy is conserved.

Analytical estimate of the dark matter density
After this small detour to explain the misalignment mechanism for the harmonic potential, let us go back to our case of interest: the ALP field with a non-standard kinetic term.
Recall that the equation of motion that we have obtained for the physical field ψ is We see that in the limit of small ψ, when N ψ 1, this reduces to the simplified case (19) and the evolution is exactly as we described in the simple real scalar field case. However, the situation is different in the regime N ψ 1. As we can expect by looking at Figure 3, the flatness of the potential away from the minimum at ψ = 0 will have the effect of delaying the start of the oscillations. Moreover, the oscillations, once they start, will not be harmonic until the damping has made the amplitude decrease enough to be in the small field regime. This means that the WKB approximation might not be as good in this case.
Although we suspect that the WKB approximation might break down when the amplitude of the oscillations is big due to the anharmoniticity of the potential, we will use it as a first approximation to solve the equation of motion and get an analytical estimate of the result. We will later contrast this to a more precise numerical computation. In the analytical approach, we will study the two regimes, where the damping is over-and under-critical, respectively, and build up the global evolution of the field by glueing together the solution for each regime. Our goal is to compute the current energy density of dark matter-like particles given an initial condition for the physical field.
As we saw, the first thing to do is to find the time when the oscillations start. In analogy with the simple case, where the condition was 3H = m φ , we use a generalisation of this formula for a non harmonic potential, namely In §3.3 we will see that this indeed works reasonably well to determine when the oscillations start, as it takes into account the flatness of the potential away from the origin. In the limit of large N ψ 1, the second derivative of the potential can be written as This turns out to be a very good approximation for intermediate and even small values of N ψ. One key difference with the harmonic case is that here the point in time when oscillations begin depends on the initial field value ϕ 0 . With this we already see that the oscillations are exponentially delayed for big N ψ: where we have assumed radiation domination so that H = 1/(2t). We now use this as an initial condition for the WKB approximation. In this approximation, the energy density of the physical field ϕ is where we have used the conservation of comoving entropy S = sa 3 to express it in terms of temperatures instead of scale factors. Using the expression for the Hubble constant during radiation domination we can express the current energy density of the field as a function of the initial condition ψ 0 , We can compare this density with the one corresponding to a harmonic potential. The result is ρ anh ρ harm so the energy density is exponentially enhanced 4 for large N and initial condition ψ 0 . The precise exponent that we obtain here should be taken as a very rough estimate. Indeed, a numerical computation is needed to get a precise result, which is what we will aim for in the following section.
As we can see in (30) the enhancement is exponential in N ψ. This implies that the field values required to yield the correct dark matter abundance are usually not too large. In the phenomenologically interesting region we usually do not need to have values for N ψ that are bigger than 50. The largest initial field values happen for N = 1 and are of order 50 in units of f a .
Another constraint that we have to care about is that the field is behaving like dark matter once it comes to dominate the dynamics of the universe, i.e. we do want to avoid having an additional phase of inflationary expansion driven by ψ. A sufficient condition for this is that the field has already started to oscillate at matter radiation equality. Making use of the more precise numerical estimate that we will obtain in the next section, we can estimate what region of parameter space satisfies this condition, This condition excludes the very small values of the mass and the decay constant in the upper left corner of Figure 1, which are already in tension with the nEDM experiment, BBN observations and the limits from [33].

Numerical computation
Having obtained a simple estimate of the cosmological evolution of the field, we now make use of a numerical solution of the equation of motion to have a more precise result. Our goal in this subsection is to quantify how much the solution for the nonlinear equation of motion (23) deviates from the harmonic case (19). Following the usual practice for dealing with anharmonicities in the ALP potential (see [49][50][51][52], [53] has a slightly different definition), we use an effective parametrisation in terms of an anharmonicity function f (ψ 0 ), such that where ρ is the energy density of the ALP field, computed late enough when it is already behaving as cold dark matter. This function only depends on the initial misalignment angle, and it should account for all the deviations from the harmonic solution. This approach is normally used to account for departures from the quadratic potential in the usual axion and ALP models. Our case is slightly different, mostly because we are dealing with an unbounded field range. As a consequence, the usual functional form for f (ψ 0 ) does not work here. Guided by the result obtained in the analytical approximation, we work with the following ansatz for the anharmonicity function: where b is a real parameter to be determined. This ansatz accounts for the exponential enhancement in energy density that we have found analytically. The normalisation needed is that f (ψ 0 ) → 1 when ψ 0 → 0, so as to recover the harmonic case in the small field limit.
The goal now is to fit the ansatz to a numerical computation of the energy density. To set the problem in a more straightforward way, we want to compare the numerical solution ofψ with the solution for the damped harmonic oscillator equation In this computation we use dimensionless quantities measured in units of m, denoted with a tilde:H,t,m . . . In these units, the time for the start of the oscillations in the harmonic case ist harm 1 = 3/2 (assuming radiation domination), and the period of the oscillations is 2π. We solve the equations numerically until we are well within the adiabatic regime in both cases (that is, when the amplitude of the oscillations has decreased enough so that the non-canonical potential is well approximated by the harmonic one). Then, we compute the energy density ρ = (1/2)f 2 aψ 2 + V (f a ψ) and extract the anharmonicity factor as the quotient of both energy densities. As we are within the adiabatic regime, ρ scales as ρ ∝ a −3 in both cases, so the quotient will stay constant. An example of the numerical solution can be seen in Figure 4.
This process is repeated for a large number of values of ψ 0 and N and we fit the results to the ansatz (33). We obtain a very good fit with a value of b = 0.56, as can be seen in Figure 5. One should note that we are fitting a two dimensional data sample with just one parameter, so finding a good fit confirms that we have chosen an adequate ansatz.
The anharmonicity function allows us to compute the energy density of the noncanonical ALP field in a very simple way, combining the harmonic solution (21) with the anharmonicity function (33). As long as we are within the adiabatic regime, the energy density in this approximation is given by The key difference between this equation and (27) is that here we use the well known solution of the harmonic equation of motion, instead of the full noninear one that arises in our non-canonical setup. All the information about the nonlinearity is encoded in the anharmonicity function, making it much more manageable.
In the analytical approach, we found that the quotient between non-canonical and canonical density scales as ρ NC /ρ C ∼ e (3/4)N ψ 0 . In the full numerical approach 5 we find a somewhat lower coefficient for the exponent of 0.56.
We have seen that a non-canonical kinetic term can indeed enhance the energy density of ALP dark matter. In the next few sections we will make use of the solutions for the cosmological evolution of the non-canonical ALP field to make predictions about its phenomenology, and to apply it to some particularly interesting cases.

Isocurvature perturbations
So far, we have assumed the initial misalignment angle θ 0 = f a φ 0 to be a constant value all throughout the universe, but of course we have to take into account fluctuations, e.g. those imprinted by inflation. We do this by taking the initial misalignment angle as a spatially varying quantity, and describing it in terms of its average and variance. Two very distinct scenarios arise, depending on whether the mechanism that gives rise to the ALP field turns on before or after the inflationary epoch of our Universe. If the ALP field was established, e.g. by spontaneous symmetry breaking, after inflation, the variance of the angle can be large even within our Hubble volume. The mean value will be φ 0 = 0 and the energy density is given by the fluctuations as well as other effects such as, e.g. the decay of topological defects [56,57]. In particular the latter contributions are not well understood and may also have some model dependence when going beyond the QCD axion.
To avoid this, we will focus on the scenario where the ALP field was present during inflation. Classically, if the ALP field was established before inflation, then the spatial variance of the field within a Hubble patch will be washed out as spacetime is stretched during inflation. This means that σ 2 φ → 0, and the misalignment field can take any value φ 0 in our Hubble patch.
However, this is not completely true, as any light field present during inflation will acquire quantum fluctuations (see, e.g. [58]). The power spectrum of such fluctuations for a canonically normalised scalar field is scale invariant, .
These fluctuations can be thought of as arising from a thermal spectrum at the Gibbons-Hawking temperature T GH = H I /(2π) [59]. As long as these fluctuations do not restore the spontaneously broken symmetry that gives rise to the ALPs, i.e., as long as 6 T GH < f a , this will imprint small fluctuations on top of the otherwise homogeneous ALP field. The corresponding fluctuations of the misalignment angle in Fourier space will have an amplitude of σ φ (k) = H I /(2πf a ). In real space, the fluctuations are of a size σ φ = γH I /(2πf a ), where γ ∼ O(1) is a dimensionless factor that effectively encodes the dispersive effect of the logarithmically divergent small k modes (see [50]). Its value depends on the length scales that we are interested in. Following [61] we will set γ = 2 for the CMB characteristic scale k = 0.05 Mpc −1 .
As the ALP has a negligible contribution to the total energy density of the universe during inflation, fluctuations in the field do not contribute to the usual curvature perturbations. Rather, they manifest themselves as fluctuations in the ratio of the number density of ALPs to the total entropy density, and are completely uncorrelated with the curvature perturbations. This is the reason why they are called entropy or isocurvature perturbations. As their interactions with other standard model particles are greatly suppressed, ALPs do not thermalise with the other species and their perturbations remain isocurvature [62]. At later stages of the cosmological evolution, the dark matter ALPs pick up a significant contribution to the energy density of the universe, and so they contribute to the temperature and polarisation fluctuations of the CMB as cold dark matter isocurvature modes.
Planck has set strong bounds on isocurvature perturbations [35], at 95% CL. Here, ∆ 2 φ (k ) and ∆ 2 R (k ) are the power spectrum of the axion and curvature perturbations at the pivot scale k , respectively. Once the value of ∆ 2 R (k ) is set (Planck gives ∆ 2 R (k ) = 2.1(9) × 10 −9 ), this translates into a bound on the axion isocurvature fluctuations.
To use this limits to constrain our scenario, we have to compute our prediction for that is, we need to evolve the fluctuations in the energy density until the time of emission of the CMB and compare them with the homogeneous average value. If the evolution of the field is linear, as it is in the case of canonical ALP models with a purely quadratic potential, the power spectrum is constant during the cosmological evolution. As a consequence, one can evaluate it at any point, such as right after inflation and before the onset of the oscillations in the ALP field. However, in any model that contains anharmonicities, the evolution at early times will be nonlinear, which implies that ∆ 2 φ will evolve nontrivially after inflation. Thus, to arrive at the correct prediction for the isocurvature perturbations, we have to track the evolution of the fluctuations until late times.
In addition to the limits from isocurvature fluctuations, the inflationary fluctuations 7 in the ALP field also forbid tuning the initial misalignment angle with arbitrary precision. In fact, there is an unavoidable limit to this tuning, and it is that our tuning precision cannot be better than the fluctuations, with σ θ = γH I /(2πf a ), as was argued in [42]. This has two related consequences. The first is that the initial misalignment angle cannot be infinitely close to zero. The requirement that the current ALP energy density is not bigger than the measured dark matter density Ω C h 2 ∼ 0.12 then sets a bound on the parameter space. This bound is model independent (as long as all the potentials are approximately quadratic for small θ) and roughly requires Secondly, if the field range is compact (as for the usual canonical ALP), an argument similar to the one above tells us that some regions of the parameter space will not yield enough energy density to account for all the dark matter. Indeed, it is not possible to tune the initial value of the field at the top of the potential with infinite precision, due to the presence of fluctuations. The requirement here is that π − θ 0 < γH I /(2πf a ). This particular limit will strongly depend on the anharmonicity of the potential, so it is not possible to give a more explicit expression. We discuss some particular cases in the next subsection. However, this last effect will not be relevant in our non-canonical model, as there we have an unbounded field range (our potential does not have a maximum).

Isocurvature perturbations for anharmonic potentials
We now present a general analytical expression to compute the isocurvature perturbations in general ALP models where the potential might have big anharmonicities. We do this using the anharmonicity function formalism that we presented in the previous section. An equivalent result was derived in [63] using the δN formalism. Here we provide a more straightforward derivation and extend the use of the formula to more general potentials.
To evaluate expression (39), we will use the fact that at t CMB the field should already be oscillating harmonically, as observations require it to behave as cold dark matter already by the time of matter-radiation equality. As we are already well within the adiabatic regime, the anharmonicity function approach will work well to describe the evolution of the energy density, which means that we can use equation (36). As fluctuations are small, we can work to linear order in σ φ to find 8 Note that even if this quantity is evaluated at t CMB , it directly depends only on the initial misalignment angle and the statistics of its fluctuations at inflation. All the information about the later evolution is encoded in the anharmonicity function.
We will now apply the formula (59) to both the case of the canonical ALP with a cosine potential and to our non-canonical model, and compare the results with the harmonic approximation.
For the harmonic case, where f (θ 0 ) = 1, we have the usual expression The constraints that one finds, for different values of the energy scale of inflation, are presented in Figure 6. The harmonic case in particular corresponds to the first column of plots.
Of course, the harmonic case can only be an approximation valid for small θ, as ALP models should preserve the shift symmetry θ + 2π. Among the potentials that satisfy this condition, the most commonly used is V (θ) = m 2 (1 − cos θ). The anharmonicity function that appears in this case was studied in [49][50][51][52]. After comparing with numerical simulations, we have decided to use a slightly different version of it, proposed in [53] and which provides a better fit to the numerical data, With this, it is easy to arrive to the following expression for the isocurvature perturbations, Note that this reduces to the harmonic result for small θ 0 . However, for angles close to π, the isocurvature perturbations are greatly enhanced. As expected, this function diverges at θ 0 = π, but as we have noted before, this limit is unattainable because of the fluctuations in the field. In Figure 6, we can see that the limits we can put on the parameter space are a bit stronger than in the harmonic case, in particular for low values of m and f a , which correspond to large values of the initial misalignment angle. Finally, we turn to the non-canonical case. The main difference with the canonical ALP, aside from the shape of the potential, is that here we are dealing with an unbounded field range. As the potential is asymptotically flat, it is always possible to enhance the production of ALPs by choosing a larger initial misalignment angle, as we saw in §3. This means that this model can always evade the limits related with to underproduction of dark matter. Using the anharmonicity function that we derived in the previous section, we find that the isocurvature power spectrum generated in this scenario is where b = 0.56. Again, this reduces to the harmonic case for small θ. The last column of plots in Figure 6 illustrate the limits that arise from the Planck data. Note that in the harmonic and canonical model featuring a compact field range a strong restriction on the parameter space is given by the requirement to produce enough dark matter (the limits arising from this condition are shaded in purple in Figure 6). As we have already argued, this limit is not present in our non-canonical setup, which features an unbounded field range. As a consequence, this model opens up a large region of parameter space, corresponding to low masses and decay constants, that was disfavoured until now. Finally let us remark that, as is well known, high scale inflation strongly constraints ALP models due to the generation of large isocurvature perturbations, which are not seen in the CMB. The tensor to scalar ratio r is strongly correlated with a high scale of inflation, so a detection of primordial gravitational waves would put a strong constraint on all axion and ALP dark matter models, including ours. Future experiments [64][65][66] are expected to increase the sensitivity in measuring r and thus the energy scale of inflation.

Coupling to QCD: Temperature dependent mass
So far, we have not assumed a coupling of the ALP to any other field. In what follows, we will allow for a coupling to gluons via a term θGG. We will study two distinct cases. First, we contemplate the possibility of having a non-canonical kinetic term in an otherwise QCD-axion model. Then, we add an extra term to the Lagrangian which, as we will see, allows us to construct a model of light ALPs that enjoys relatively strong gluon couplings.

The QCD axion
In this section we will focus on the QCD axion as introduced by Peccei and Quinn as a solution to the strong CP problem in quantum chromodynamics (QCD) [4][5][6].
The Lagrangian for the canonically normalised axion field is now and as usual we can define the angle θ = φ/f a , so that θ ∈ (−π, π]. For our modification with a non-canonically normalised field, we have and after we perform a field redefinition to have it canonically normalised, we find the Lagrangian There is just one difference that makes the QCD axion case particular, and it is that here the energy scale appearing in the potential is fixed by QCD to be [43] Λ It is easy to see that the mass of the axion, m a , is given by f a m a = Λ 2 QCD . It is important to note that the numerical value quoted above is only valid at zero (or very low) temperatures. Indeed, the axion potential is affected by finite temperature effects, such that the mass of the axion varies with temperature. At low temperatures below the QCD critical temperature T crit ∼ 160 − 170 MeV, the mass remains roughly constant 9 . That said, much of the dynamics that is of interest to us will happen in the early universe, at temperatures 9 The small temperature dependence can be computed using chiral perturbation theory as in [43] close or above T crit . There are different ways to compute the temperature dependence of the axion mass [42,43,49,52,67,68]. The function that controls the temperature dependence of the axion mass is the topological susceptibility χ(T ), which is usually parametrised as a power law: Here we will use 2α = −7.1 and χ 0 = 0.11, from recent lattice computations [68] that are consistent with the instanton values up to an overall normalisation factor. We see that the main effect is that the mass of the axion is approximately constant until T crit , and then it drops as a power law, so that the axion is essentially massless at high temperatures. The most important implication of the temperature dependent mass is that a smaller mass at early times can delay the start of the oscillations of the field, which in turn results in a higher energy density of axionic dark matter today. This happens both for the canonical and non-canonical axion models.

Anharmonicity function and isocurvature perturbations revisited: Temperature dependence
We have seen that coupling ALPs to QCD through φGG results in a temperature-dependent mass for the ALP, both in the canonical and non-canonical setup. This of course has an impact on its cosmological evolution, which can be of importance in computing observables such as the isocurvature perturbations that we discussed in §4. To account for this effect, we will modify the anharmonicity function formalism that we introduced in §3.3 to incorporate the temperature dependence. That is, we want to compute evaluated at a point in time late enough so that the anharmonic and temperature-dependent axion field has already entered the adiabatic regime. For definiteness, we will use the following expression for the axion mass, First of all, we note that this temperature dependence will only have an effect if the field starts oscillating before the QCD critical temperature T crit . In the harmonic limit, this means that if the mass is smaller than m * a , defined by 3H(T crit ) = m * a , the field will have acquired its late-time mass by the time it starts oscillating. Thus, the later evolution of the field will be insensitive to the temperature effects that happened earlier on. In terms of decay constants, this sets a distinct scale  (46). Both the anharmonicities of the potential and the temperature dependence of the mass are taken into account through the anharmonicity function defined in (58). Our results differ slightly from the ones obtained in [69] and [63] due to the fact that we are using the more recent data from the Planck satellite and a different anharmonicity function.
If we take into account the anharmonicities of the potential, it might happen that the start of the oscillations is delayed until after T crit , even if f a < f * a . The condition to be in this regime is that the initial misalignment angle θ 0 is larger than some value θ * 0 (f a ). This value is given for a general anharmonicity function 10 by For the case of a canonical axion with a cosine potential like in (46), we find whereas in the non-canonical case (48), we find   (47). Both the anharmonicities of the potential and the temperature dependence of the mass are taken into account through the anharmonicity function defined in (58).
For any set of decay constants and initial misalignment angles that satisfy f a < f * a and θ 0 < θ * 0 , we compute F T for a generic anharmonic potential, finding The details of the derivation of this result are given in Appendix §B. Here we see that the result depends critically on the exponent of the temperature-dependence of the axion mass at high temperatures above the QCD critical temperature. To sum up, we can write the full temperature-dependent anharmonicity function as follows, if f a < f * a and θ 0 < θ * 0 .
With this, we can use the same approach as in §4.1 to compute the isocurvature perturbations, this time using the temperature-dependent anharmonicity function, We apply this formula for both the canonical QCD axion and for our non-canonical model, and obtain the results presented in Figure 7 and Figure 8, respectively. For the canonical QCD axion, our results are an update from the ones obtained in [69] and [63], as we are using the more recent data from the Planck satellite and a better fitting anharmonicity function.

ALPs coupled to QCD with f a m << Λ 2 QCD
Let us now study the possibility of an ALP having a coupling to the GG term while satisfying f a m << Λ 2 QCD . This is an interesting region of the parameter space, as ALPs that satisfy these conditions may be found by looking for an oscillating nucleon or atomic electric dipole moment. There exist a number of proposed laboratory searches focusing on this direction [27,31,32,70] .
However, we have seen that coupling the ALP to QCD via a term proportional to GG induces an irreducible contribution to the mass, given by (49). Explicitly this contributes to the square of the axion mass as given in [43]. This contribution will also have a temperature dependence as described by (50). A priori, this irreducible contribution to the axion mass seems irreconcilable with the condition f a m << Λ 2 QCD [28]. The only known way of circumventing this caveat is to precisely cancel this contribution with an additional, fine-tuned term in the Lagrangian. Acknowledging the flaws of this ad hoc approach, we follow it and study the phenomenology of such models when allowing for a non-canonical kinetic term.
At the level of the Lagrangian, we add an extra term to the potential so that it becomes In principle there can exist a phase difference between both contributions. For our purposes, it will be necessary to require that this phase difference vanishes, so we will take α = 0. This can be viewed as equivalent to asking for a separate a solution to the strong CP problem. In principle any integer n is possible but for simplicity we will limit ourselves to the n = 1 case. In the small φ limit, this potential induces a mass for the ALP where m 0 f a = Λ 2 0 and recall that m a is completely fixed by f a as in equation (60). It is then possible to choose m 0 so that we get any zero-temperature mass for the ALP, i.e. we can set m 2 0 = m 2 a (T = 0) − m 2 . We are interested in the m 2 m 2 a (T = 0) regime. The full mass can then be expressed as Because at early times the QCD contribution is strongly suppressed, in that regime we have m 2 (T ) < 0. We will use the following simplified expression for the temperature dependent mass of the ALP Note that the negative mass does not indicate an unstable potential but only that φ = 0 is not the minimum at that time.

Canonical case
As a first step, we implement the mass subtraction and the resulting temperature dependence in an ALP model with a canonically normalised scalar field with potential given by with m(T ) defined in (64). The most relevant feature of this scenario is that before the QCD phase transition, the potential is minimised at θ = π rather than at θ = 0. Accordingly, at early times the field evolves towards its minimum at π, around which it will oscillate with damped amplitude. Then, after the QCD phase transition, the potential rapidly acquires its late-time shape, with a minimum at the origin. The field thus oscillates around its CP-conserving value θ = 0 at late times. The main role of the first set of oscillations is to set the initial condition for the second one to be close to π. We refer to Figure 9 for a cartoon explaining this evolution. It should be noted that this discussion is only valid if f a m Λ 2 QCD , that is, if we lie to the left of the QCD axion band in Figure 10. In the other limit, ie f a m Λ 2 QCD , the contribution of the QCD mass is negligible and we recover the usual constant mass ALP scenario.
Let us now be a bit more quantitative. Initially, H is large and the field is stuck at its initial value θ 0 . Then, as long as the early-time mass m a (0) overcomes the Hubble friction before the QCD phase transition, the field will oscillate around π. The condition for this to happen is roughly f a 10 17 GeV, but this value can be modified by the anharmonicities depending on the initial misalignment. These oscillations continue until the temperature decreases to T crit , at which time the amplitude is approximately given by Here, the anharmonicity function is given by (43) and T 1 is defined by 3H(T 1 ) = m a (0). The value of θ crit gives the initial condition for the oscillations that happen after the QCD phase transition, now around θ = 0 and with frequency given by the late-time mass m. Typically, θ crit is very close to π so the anharmonicites of the potential will play a key role. Taking this into account, we can compute the energy density of the oscillating scalar field as where 3H(T 2 ) = m. We can then determine in what region of parameter space the right dark matter abundance can be generated with an initial misalignment angle θ 0 of order O(1). It is possible to either enhance or suppress the energy density given in (67) by tuning the initial misalignment angle closer to zero or π. However, due to equation (59), there is an enhancement of the isocurvature perturbations each time the field gets close to a maximum of the potential, where the anharmonicity function becomes large. Because of this, the available tuning of the initial misalignment angle is very limited in this scenario due to the stringent constraints on isocurvature fluctuations. Figure 10 shows how the allowed parameter space shrinks for larger values of the Hubble scale of inflation. Despite the strong isocurvature constraints, we can see that this scenario populates some unexplored regions of parameter space to the left of the QCD axion line that could be probed by upcoming experiments looking for ALPs. All the other limits presented in Figure 1 are also applicable in this scenario, as they only depend on the dynamics of the field after the QCD phase transition.

Non-canonical case
We now want to implement the temperature dependent potential (65) in our non-canonical ALP scenario. In terms of the cosmological evolution of the field, this is effectively done by writing This is the same potential as we had before, except that for high temperatures T > T crit the mass squared will be negative and will be a function only of f a , as given in (64). This tells us that, depending on the value of the parameters m and f a , we will have two very different behaviours, which qualitatively can be understood as follows. First, if the field does not start rolling until after the QCD phase transition, then all the dynamics and the observables will not be affected at all by the features of the potential at high temperatures. This is because there is no evolution while the field is frozen by Hubble friction. Only after it has acquired its late time mass m does it start rolling, and thus the cosmological evolution is exactly as we computed in §3. However, the key difference is that now the ALP is coupled to GG, so it may be tested by observables and experiments that exploit this coupling.
The other option is, of course, that the field starts rolling before the QCD phase transition. Then the dynamics can depend strongly on the initial conditions and is rather complicated. However, we will see that this scenario leads to an overproduction of ALPs whose energy density exceeds the observed CDM one. As we are only interested in ALPS as dark matter candidates, the second scenario is not interesting for us and we just need to focus on the first one.
Let us now be more quantitative and compute what region of the parameter space allows for ALP dark matter with a non-canonical kinetic term and coupled to QCD. As we have anticipated, this ALP will only be a good dark matter candidate if its evolution is frozen until after the QCD phase transition. Then, the present ALP energy density will only depend on f a , the present mass m and the initial misalignment angle ψ 0 . The latter is given by equation (36), and satisfies where T 1 is the temperature at which the oscillations start.
We now need to find what the region of the parameter space is where the field starts oscillating only after the QCD phase transition. The QCD phase transition happens at a temperature of around T crit ∼ 160 MeV, which corresponds to a Hubble parameter of H(T crit ) ∼ 10 −11 eV. By asking that 3H(T crit ) > |V (ψ 0 )| 1/2 , we get the condition f a 10 11 GeV This region is plotted in Figure 11, together with the further cosmological and astrophysical bounds that restrict the parameter space.
Finally we still have to justify our claim that if the field starts oscillating before the QCD phase transition we always get an overproduction of ALPs. For a given (m, f a ), any initial misalignment angle bigger than the one given by (69) will lead to an energy density in ALPs greater than the observed dark matter one. But if the condition (70) is not satisfied, then the field will start rolling towards bigger ψ values, because m 2 (T ) < 0 at high temperatures. Thus, the effect of the rolling at high temperatures is to drive the field away from the required misalignment angle to give the correct dark matter abundance. This statement is independent of what misalignment angle we start with, and thus rules out ALPs in the region coloured in white in Figure 11 as dark matter candidates 11 .  Figure 11: Parameter space for ALPs with a non-canonical kinetic term of the form (11) coupled to QCD via a GG term, along with constraints coming from its cosmological evolution and searches for an oscillating EDM. Each panel represents a different value of the parameter N . This scenario can provide the right dark matter density in the yellow shaded region, while the areas excluded by overproduction of dark matter or the condition (31) to avoid a second period of inflation are coloured in white. The brown region is excluded by re-analysing data originally intended to search for a static neutron EDM in order to look for an oscillating one [30]. The dark green region in the first figure is inconsistent with the production of the observed abundance of light elements during Big Bang Nucleosynthesis [28]. This limit is effective only for N < 4 and absent in the other figures. Similarly, the limit from [33] corresponding to the ALP field being sourced at the Sun only applies for small values of N , while the Earth one stays valid in all cases. Finally, f a is (softly) bounded from above by the requirement that it does not exceed the Planck scale, and from below by the supernova limits estimated in [29].

Big Bang Nucleosythesis
those used in α-attractor models for inflation we find a significantly enlarged parameter space for dark matter. In particular, regions with larger couplings -where canonical ALPs are underproduced -now become viable, offering interesting possibilities for near future experiments.
For the production via the misalignment mechanism the key feature of the non-canonical kinetic term is that today's ALP energy density is enhanced due to a delay in the start of the oscillations. This arises because the effective potential is flattened by the growing noncanonical kinetic term, which also makes the field range of the physical field unbounded. As a consequence, any combination of mass and decay constant can generate enough ALP energy density to account for all the dark matter that we observe in the universe.
An important cosmological constraint arises from isocurvature fluctuations imprinted by inflation. To apply these constraints to our scenario we give a simple derivation of the size of isocurvature fluctuations in general models with arbitrary potential and even a temperature dependence of the potential. As a useful crosscheck we have updated the isocurvature constraints [63,69] using the newest Planck data [35] and the most recent results for the QCD topological susceptibility [68]. The result can be found in Figure 7. In our non-canonical setup the isocurvature constraints are even slightly weaker as can be seen in Figure 8.
An interesting non-trivial situation arises if the ALP is coupled to the strong interactions, i.e. via a term ∼ φGG. This is of particular interest since a number of experiments are currently searching for ALP dark matter with this coupling [30][31][32]. The coupling to gluons leads to two non-trivial features: the generation of a temperature-dependent, irreducible contribution to the ALP mass and an effective ALP field value dependent nucleons mass. The former naively makes large parts of the low mass region explored by current experiments inaccessible [28]. This can be avoided by invoking a precise cancellation with an additional term in the ALP potential (with or without non-canonical terms). The latter leads to strong constraints from Big Bang Nucleosynthesis. These are significantly weakened in our scenario with a non-canonical kinetic term. This opens up significant parameter space that can be explored in near future experiments such as Casper [31] and HeXeniA [32], as well as EDM storage rings [77].

A Effect of higher-order poles in the kinetic function
In this appendix we briefly study how our results change if we allow our non-canonical kinetic term for the ALP field to have a pole of arbitrary (even) order. We work with the Lagrangian (4), this time with the kinetic function given by Note that with this definition the order of the pole is 2p. In the main body we have focused in the p = 1 case. As opposed to the p = 1 case, for a general value of p it is not possible to find an exact analytic expression for the transformation to the canonically normalised field ψ(φ). However, we can find an approximate expression, valid close to the pole at φ/f a = π/(2N ), by expanding K(φ) in a Laurent series and keeping only the leading divergent term. With that, we can then proceed as in §3.2 and obtain an estimation for the enhancement in the relic density. The result is Due to the now unbounded field range, a significant enhancement is possible. Nevertheless, it may seem that the enhancement effect is much weaker in the p > 1 cases, which could seem counter-intuitive. However, we must note that the field redefinition ψ(φ) is different for different values of p, which means that it is not so obvious to compare the distinct cases just by looking at (76). In order to be able to compare, we can recast (76) in terms of the non-canonically normalised field φ, which avoids the problem of the p-dependent field redefinition. Doing so, we can write which is valid for all p ∈ N. Looking at (77) we can confirm that there exists an enhancement in the ALP relic density for all values of p. What's more, looking at it from the point of view of the non-canonically normalised field, the effect is stronger the higher the order of the pole, as one would naively expect.

B Temperature dependent anharmonicity function
In this section we detail how to implement the effects of the temperature dependence of the QCD axion mass into the anharmonicity function [53,69]. The approach that we follow allows us to analytically upgrade any anharmonicity function that does not include temperature effects into a full temperature dependent anharmonicity function. The derivation that we present is valid for any scalar field whose potential can be factorised as where V 0 (φ) is a temperature-independent potential that has a minimum around which the field can oscillate (possibly anharmonically), and the temperature dependence acts only as a scaling. This is the case for general axion models, including those we study in this paper. If there is no temperature dependence at all, then m(T ) ≡ m and working as in §3.3 we can express the energy density of the field with an anharmonic potential as where ρ harm is the solution to the harmonic case described in §3.1 and f (θ 0 ) is the anharmonicity function that depends on the initial value of the dimensionless field θ = φ/f a . At the effective level, we can think that the only effect of the anharmonicities is to change the time (or temperature) at which the oscillations start. This is of course not what actually happens, but with this approach we will be able to make a good estimate of the energy density of the field at late times. In this picture, we have that which follows from the dependence of the WKB solution (27) on the temperature at which oscillations start 13 . With this we can extend the equation for the condition of the start of the oscillations in the harmonic case, 3H(T harm S ) = m, to the anharmonic case, using an effective mass that encodes the effects of the anharmonicity. It reads 3H(T anh S ) = m f (θ 0 ) −2/3 .
This expression will be useful later on. We now assume that there is a temperature dependent mass that evolves as (52), as is the case for axion models. As was argued in the main body of the paper, if T anh S < T crit , the oscillations start after the QCD phase transition, when the mass has already attained its low-temperature value, and the temperature dependence has no effect on the later evolution of the field. In the harmonic case, this happens if the zero-temperature mass is smaller than m * given by 3H(T crit ) = m * , that is, m * = 3 · 1.66 g (T crit ) T 2 crit m Pl 6.6 · 10 −11 eV.
In terms of decay constants, this translates into a maximum value f * a 8.7 · 10 16 GeV above which the temperature dependence does not play a role. In the anharmonic case, it can happen that the anharmonicities delay the start of the oscillations beyond T crit , even if f a < f * a , if the initial misalignment angle is large enough. We can use equation (80) to find an expression for this critical value θ * 0 . Writing it in terms of f * a , it reads from where θ * 0 can be obtained once an explicit anharmonicity function is chosen. Finally, if both f a < f * a and θ 0 < θ * 0 , then T anh S > T crit and the temperature dependent evolution of the mass will have an impact on the oscillating behaviour of the field. First of all, we compute how the onset of the oscillations is modified. For this purpose, we can just substitute the constant mass m for the temperature dependent one m(T ) in equation (81). Using the expression for m(T ) given in (52), we have the condition To simplify the notation, we have denoted T S the temperature at which the oscillations start if we take into account both the anharmonicities of V 0 (φ) and the temperature dependence of m(T ). We can recast this equation in terms of decay constants, finding To continue, we use the expression for the energy density of an oscillating scalar field with a slowly varying mass, which can be found for instance in [17] At low temperatures below the QCD critical temperature, the quotient between this expression and the corresponding one for the harmonic case is But this is precisely what we need to define the temperature dependent anharmonicity 14 This expression is a generalisation of the WKB approximation presented before.