Gravitational Wave and CMB Probes of Axion Kination

Rotations of an axion field in field space provide a natural origin for an era of kination domination, where the energy density is dominated by the kinetic term of the axion field, preceded by an early era of matter domination. Remarkably, no entropy is produced at the end of matter domination and hence these eras of matter and kination domination may occur even after Big Bang Nucleosynthesis. We derive constraints on these eras from both the cosmic microwave background and Big Bang Nucleosynthesis. We investigate how this cosmological scenario affects the spectrum of possible primordial gravitational waves and find that the spectrum features a triangular peak. We discuss how future observations of gravitational waves can probe the viable parameter space, including regions that produce axion dark matter by the kinetic misalignment mechanism or the baryon asymmetry by axiogenesis. For QCD axion dark matter produced by the kinetic misalignment mechanism, a modification to the inflationary gravitational wave spectrum occurs above 0.01 Hz and, for high values of the energy scale of inflation, the prospects for discovery are good. We briefly comment on implications for structure formation of the universe.


INTRODUCTION
The thermal history of the very early Universe remains uncertain. It involved a sequence of eras, where each era was characterized by a certain expansion rate. The expansion rate is key to understanding the physical processes occurring during any era and is determined by ρ(a), the dependence of the energy density on the Friedmann-Robertson-Walker scale factor a. From the precise observations of the cosmic microwave background (CMB), we know that as the temperature cooled through the eV region, the universe transitioned from being dominated by radiation, with ρ(a) ∝ 1/a 4 , to being dominated by matter, with ρ(a) ∝ 1/a 3 . This matter-dominated era lasted until relatively recently when the universe entered an era apparently dominated by vacuum energy, with ρ independent of a. Furthermore, at the time of Big Bang Nucleosynthesis (BBN), when the temperature was in the MeV region, the universe was radiation-dominated, and there was likely a very early era of vacuum domination known as inflation, when ρ was independent of a. Since this is the total observational evidence we have of the very early evolution of our universe, there are clearly many possible cosmological histories, each having a different sequence of transitions between eras of differing ρ(a).
What is the underlying field theory and cosmology that leads to an era of kination? Once kination starts, it is easy to end since the kination energy density dilutes under expansion much faster than radiation, so a transition to radiation domination will occur. But how does kination begin? Going to earlier times during the kination era, the kinetic energy density of the scalar field rapidly increases. This issue is particularly important for primordial gravitational waves, whether produced from quantum fluctuations during inflation and entering the horizon during a kination era or from emission from cosmic strings during a kination era, since the spectrum for both increases linearly with frequency. This UV catastrophe must get cut off by the physics that initiates the kination era, and hence the peak of the gravitational wave distribution will have a shape determined by this physics.
Recently a field theory and cosmology for kination was proposed by two of us: axion kination [34]. An approximate U (1) global symmetry is spontaneously broken by a complex field. Early on there are oscillations in both angular and radial modes, which we call axion and saxion modes, in an approximately quadratic potential. The radial oscillation results from a large initial field displacement, and the angular mode is excited by higher-dimensional operators that break the U (1) symmetry. At some point the radial mode is damped, but by this time the "angular momentum" in the field, that is the charge density, is conserved, except for cosmic dilution, and a period of circular evolution sets in. This era is assumed to occur when the radial field is much larger than the vacuum value, f a , the symmetry-breaking scale for the axion. In fact, at first the trajectory is not quite circular, as the cosmic dilution of charge leads to a slow inward spiral of the trajectory. If the energy density of the universe is dominated by the scalar field energy, this inspiral era is a matter-dominated era with ρ(a) ∝ 1/a 3 . This era ends when the radial mode settles to f a ; the potential vanishes and the axion energy density is now entirely kinetic so that an era of kination ensues. Kination ends when the axion energy density falls below that of radiation.
The rotation of an axion field was used in [34] to generate a baryon asymmetry via Axiogenesis and in [35] to generate axion dark matter via the Kinetic Misalignment Mechanism.
In fact, such schemes did not rely on the axion field energy becoming larger than the radiation energy, so an era of kination was possible but not required. In this paper we study the implications of a kination-dominated era from this mechanism, where the era is cut off in the UV by an early matter-dominated era. We consider both the QCD axion [36][37][38][39] that solves the strong CP problem and generic axion-like particles (ALPs).
The early matter-dominated era regulates the spectrum of gravitational waves from inflation or cosmic strings. After the linear increase with frequency from the kination era, the magnitude of the gravitational wave spectrum decreases, producing a triangular peak in the spectrum. Interestingly, the shape of the peak contains information about the shape of the potential of the complex field and could reveal the origin of the spontaneous U (1) symmetry breaking resulting in an axion as a Nambu-Goldstone boson.
It is natural to expect this kination era to occur early in the cosmic history, ending well before BBN. Remarkably, for low f a and large charge density, a late era of kination domination can occur after BBN, but before matter-radiation equality. This possibility arises because axion field rotations do not generate entropy. In such scenarios, the beginning of the matter-dominated era is constrained by BBN and the end of the kination-dominated era is constrained by the CMB.
We also examine the implication of the NANOGrav signal [40] to axion kination. The signal may be explained by gravitational waves emitted from cosmic strings [41,42]. We discuss how the gravitational wave spectrum from low to high frequency is modified by axion kination and the modification, including the peak, can be detected by future observations. In Sec. 2, we discuss the above mechanism of axion kination cosmology, including two theories for the potential of the radial mode. We study the transition from the early matterdominated era to the kination-dominated era, and also the thermalization of the radial mode, as this leads to important constraints on the signals. We derive the dependence of the matter and kination transition temperatures on the axion model parameters. In Sec. 3, we analyze the constraints on axion kination from both BBN and the CMB. In Sec. 4, we review the kinetic misalignment mechanism and axiogenesis and derive predictions for the parameters of axion kination. In Sec. 5, we compute the spectrum of gravitational waves produced from inflation and cosmic strings in the axion kination cosmology. The spectra depend on and can be used to infer the kination and matter transition temperatures and hence axion parameters. We discuss whether future observations of gravitational waves can detect the imprints of axion kination. We pay particular attention to the parameter space where the axion rotation also lead to dark matter or to the baryon asymmetry and also to the case of the QCD axion. Sec. 6 is dedicated to summary and discussion.

Axion rotations
In field-theoretical realizations of an axion, the axion field φ a is the angular direction θ ≡ φ a /f a of a complex scalar field P , where S is the radial direction which we call the saxion. In the present universe, S = f a , the decay constant, spontaneously breaking an approximate U (1) symmetry. The potential for P contains small explicit U (1)-breaking terms that give a small mass to the axion. For simplicity, we take the domain wall number to be unity.
In addition, there may be explicit U (1)-breaking operators of higher dimension, V = P n M n−4 + h.c., (2.2) where M is a dimensionful parameter. In the early universe, S may take on a field value much larger than f a . The potential gradient to the angular direction given by the higherdimensional operator then gives a kick to the angular direction and the complex field begins to rotate. As the universe expands, the field value of S decreases and the higher-dimensional operator becomes ineffective. The field P continues to rotate while preserving its angular momentumθS 2 up to the cosmic expansion. Such dynamics of complex scalar fields was proposed in the context of Affleck-Dine baryogenesis [43]. The angular momentumθS 2 is nothing but the conserved charge density associated with the U (1) symmetry. It is convenient to normalize this charge density by the entropy density of the universe s, which remains constant as long as entropy is not produced. As we will see, the charge density must be large enough to obtain kination domination. This can be easily achieved in our scenario because of the large initial S.
Soon after the field rotation begins, the motion is generically a superposition of angular and radial motion and has non-zero ellipticity. Once P is thermalized, the radial motion is dissipated, while the angular motion remains because of angular momentum conservation. One may think that the angular momentum is transferred into particle-antiparticle asymmetry in the thermal bath. It is, however, free-energetically favored to keep most of the charge in the form of axion rotations as long as S T [34]. The resultant motion after thermalizaion is a circular one without ellipticity. The parameter space that leads to successful thermalization is investigated in Sec. 2.3.

Axion kination
The evolution of the energy density of the axion rotations depends on the shape of the potential of the saxion. A very interesting evolution involving kination is predicted when the saxion potential is nearly quadratic at S f a [34]. Such a potential arises naturally in supersymmetic theories, where the saxion is the scalar partner of the axion: the saxion potential may be flat in the supersymmetric limit and generated by supersymmetry breaking.
For example, the soft supersymmetry breaking coefficient of P * P = S 2 /2 may be positive at high scales but evolve under renormalization to negative values at low scales. This triggers spontaneous breaking of the U (1) symmetry, which can be described by the potential [44], which is nearly quadratic for S f a . Another example is a two-field model with soft masses, Here X is a stabilizer field that fixes the symmetry breaking field P andP on a moduli space PP = v 2 P . For P v P orP v P , the saxion potential is dominated by the soft mass m P and mP , respectively. Without loss of generality, we choose P to be initially much larger than v P and identify the saxion with the radial direction of P . We neglect possible renormalization running of m P which modifies the saxion potential only by a small amount.
For these potentials, the axion rotations evolve as follows. When S f a , the potential of S is nearly quadratic, and the equation of motion of the radial direction requiresθ 2 = V (S)/S m 2 S . The conservation of the angular momentum,θS 2 ∝ a −3 , then requires S 2 ∝ a −3 . Here a is the scale factor of the universe, not to be confused with the axion field which we denote as φ a . The potential energy ∼ m 2 S S 2 and the kinetic energy ∼θ 2 S 2 are comparable. Once S decreases and S f a , the conservation of the angular momentum requiresθ ∝ a −3 . The kinetic energy dominates over the potential energy. The scaling of the energy density in these two regimes is (2.6) The scaling of the energy density naturally leads to kination domination [34]. When S f a , the rotation behaves as matter and is red-shifted slower than radiation is, so the universe may become matter-dominated by the axion rotation. Once S approaches f a , kination domination by the axion rotation begins.
Throughout most of this paper, we adopt a piecewise approximation where ρ θ ∝ a −3 for S > f a and ρ θ ∝ a −6 for S = f a . We will comment on how the actual evolution differs from this. Within this approximation, the Hubble expansion rate H as a function of the temperature T is given by Here T RM is the temperature at which the matter domination (MD) by the axion rotation begins, T MK is the temperature at which the kination domination (KD) begins, and T KR is the temperature at which the KD ends and radiation domination (RD) begins.
In Eq. (2.7), it is assumed that P is thermalized when the rotation is still a subdominant component of the universe. It is also possible that the thermalization occurs after the rotation dominates the universe. In this case, the energy associated with the radial mode, which is comparable to or larger than the angular mode, is converted into radiation energy.
Then the universe evolves as RD → MD → RD → MD → KD → RD, and the scaling in Eq. (2.7) is applicable to the last four eras. The second RD is very short when the radial and the angular component of the initial rotation is comparable, which naturally occurs in supersymmetric theories; see [45] for details. Y θ computed after the initiation of the rotation receives entropy production from the dissipation of the radial mode, and is conserved again afterward. It is this final value of Y θ that concerns us, and we do not investigate how it is related to the UV parameters of the theory.
The cosmological progression from RD → MD → KD → RD described above is determined by three parameters, (f a , m S , Y θ ), and the three temperatures (T RM , T MK , T KR ) can be expressed in terms of these, , (2.9) . (2.10) The kination-dominated era exists if T KR < T RM , The scalings of (2.7) imply that these three temperatures are not independent, but are related by T 3 MK T RM T 2 KR . The expansion history of the universe is therefore determined by two combinations of (f a , m S , Y θ ) such as (m S f a , Y θ ), but other phenomenology depends also on the third combination. In what follows, according to convenience, we use a variety of ways of spanning the 3-dimensional parameter space. For discussion of axion physics we must include the axion mass m a as a fourth parameter. The axion mass is determined by f a for the QCD axion [36][37][38][39]. Note that T RM T KR ∝ m S f a , so one can in principle determine the product m S f a by a measurement of T RM and T KR through gravitational wave spectra discussed in Sec. 5.
Moreover, additional theoretical considerations such as the baryon asymmetry and/or axion dark matter abundance from the axion rotation, as discussed in Sec. 4, will help determine m S and f a individually. If f a and m a are also measured by axion experiments, then the theory parameters are over-constrained so that the theory can be confirmed or ruled out.
A unique feature of our kination scenario is that matter domination preceding kination domination ends without creating entropy. This is quite different from usual matter domination, where matter domination ends by dissipation of matter into radiation creating a huge amount of entropy. Because of the absence of entropy production in our scenario, matter and kination domination can occur even after BBN and before recombination.
The precise evolution of the universe differs from the piecewise approximation and depends on the saxion potential. The evolution for the one-field model of Eq. (2.4) is derived in [34] and reviewed in Appendix A, and is shown by the orange solid line in Fig. 1. Beyond the piecewise approximation, there is no sharply defined T MK , so we first define T RM and T KR by the equality of the axion energy density with the radiation energy density and then KR . The transition from matter to kination domination is not sharp, but still occurs within a temperature change of O(10). The evolution for the two-field model of Eq. (2.5) is derived in Appendix A and is shown by the blue-dashed lines. The evolution depends on the ratio m P /mP , but the transition is sharper than the one-field model and occurs within a temperature change of O(1). As we will see, this difference shows up in the spectrum of primordial gravitational waves and allows the determination of the saxion potential.

Thermalization
As discussed in the previous subsection, the motion of the field P is initiated in both angular and radial components, and the energy density associated with the radial mode is comparable to or more than the rotational energy. Since the radial mode also evolves as matter, if it is thermalized after S reaches f a , no kination-dominated era is present. Thus earlier thermalization is required. In the simplest case, we consider a Yukawa coupling between the saxion and fermions ψ andψ that couple with the thermal bath, L ⊃ y ψ Sψψ. (2.12) The simplest possibility is a Standard Model gauge charged fermion, but we may also consider a dark sector fermion. The thermalization rate is given by [46], where b is a constant and is O(0.1) when the coupling of the fermion with the thermal bath is O(1). The fermion is heavy in the early universe because of a large saxion field value, m ψ = y ψ S, while the fermions themselves need to be populated in the thermal bath in order to thermalize the saxion at the temperature T th . Such a requirement, y ψ S th ≤ T th , leads to an upper bound on the Yukawa coupling as well as on the rate (2.14) We obtain the same bound for a saxion-scalar coupling. For gauge boson couplings, which arise after integrating out charged fermions or scalars, the thermalization rate [46], so the constraints on the parameter space for this case can be obtained by putting b = 10 −5 in the following equations.
Here Y θ is determined from a fixed T KR using Eq. (2.10). In this case, the thermalization constraints can be imposed by the consistency conditions: 1) the saxion field value at thermalization must be larger than f a , i.e., S max th ≥ f a ; otherwise, thermalization would not occur or the Universe would not be kination-dominated, 2) the radiation energy density after thermalization is at least that of the saxion, i.e., π 2 g * (T th )T 4 th /30 ≥ m 2 S (S max th ) 2 , where the inequality is saturated when the saxion dominates and reheats the universe. In fact, upon assuming the existence of kination domination, T RM > T KR , condition (2) automatically guarantees that (1) is satisfied. Therefore, only condition (2) is relevant and leads to the constraint on the saxion mass m S 1.5 × 10 5 GeV b 0.1 (2.17) This thermalization constraint is shown as the green regions in the figures we will show in the following sections. For consistency with the assumption of the rotation in the vacuum potential, the thermal mass of m S must be subdominant to the vacuum one, y ψ T th < m S . However, using condition (2), one finds that this constraint becomes y ψ < (π 2 g * (T th )/30) 1/2 T th /S th , which is always weaker than the earlier constraint from requiring ψ in thermal equilibrium, When the thermalization temperature is much lower than the QCD scale, additional constraints may become important [47]. For example, the energy density deposited into the bath (or dark radiation) at late thermalization may contribute to excessive ∆N eff since this energy deposit cannot be absorbed by the Standard Model bath nor diluted by the change of g * in the Standard Model across the QCD phase transition. Effectively, the constraint from condition (2) above is strengthened by replacing g * (T th ) by 7∆N eff /4. We impose the limit ∆N eff > 0.17 from the CMB and BBN [48].
In the case when T RM O(GeV), the saxion has to be very light and is thus subject to quantum corrections, which require the coupling y ψ < 4πm S /m soft where m soft is the soft mass of ψ's scalar partner,ψ. We do not impose this constraint since one can simply assume that m soft is generated in the same way as and is of the same order of m S , in which case the constraint is trivially satisfied.
One may expect additional constraints from the relic density and warmness of ψ especially when f a is very small and ψ may be very light. However, one may consider a model where ψ has a sufficiently large vector-like mass and freezes out non-relativistically much before BBN (see Ref. [47] for details), or ψ is dark gauge-charged and effectively annihilates into massless dark gauge bosons.

COSMOLOGICAL CONSTRAINTS
The axion rotations lead to matter and kination-dominated eras. If these happen close to BBN or recombination, the modified expansion rate alters primordial light element abundances or the spectrum of the CMB. Constraints from BBN divide kination into two paradigms -"early kination" for T KR > O(MeV) and "late kination" with T RM < O(10 keV).
In this section, we discuss the constraints on the axion rotation on both early and late kination from BBN and the late kination from CMB.

BBN
When kination domination occurs before BBN, the strongest constraint comes from the helium abundance, since it is sensitive to the freeze-out of neutron-proton conversions, which occurs at an early stage of BBN. Using AlterBBN [49,50], we show the prediction on the primordial helium abundance as a function of T KR with varying the baryon abundance within values allowed by Planck 2018 (TT,TE,EE+lowE) [51], together with the constraint on the abundance [52], in the left panel of Fig. 2. The width of the prediction originates from the uncertainty in nuclear reaction rates and the neutron lifetime. From this, we obtain the constraint T KR 2.5 MeV.
Here we use Planck's allowed range for the baryon abundance with the BBN consistency condition imposed on the helium abundance for the following reason: the BBN prediction for the helium abundance with kination does not deviate from the standard prediction as much to make the helium abundance a free parameter (see Fig. 40 of [51]). Since BBN and CMB results for baryon abundance are in excellent agreement, relaxing the consistency condition and allowing the baryon abundance to range more freely gives very similar results for the allowed parameter space.
When kination domination occurs after BBN, but before recombination, the strongest constraint comes from the abundance of deuterium whose destruction freezes out at a late stage of BBN. We show the prediction on the primordial deuterium abundance as a function of T RM in the right panel of Fig. 2 from which we obtain T RM 6 keV.

CMB
The case of an early kination-dominated era with T KR > 2.5 MeV has no observable impact on the CMB. On the other hand, in the case with T RM < 6 keV the modified expansion rate of the universe can potentially lead to significant deviations in the evolution of modes on scales probed by the CMB.
The angular size of the sound horizon at the surface of last scattering, which is precisely measured, is one quantity that can be altered by a modified cosmic expansion history. We will develop some intuition for how the sound horizon is changed by assuming a kinationdominated era with T RM < 6 keV and T KR > T eq , where T eq 0.8 eV is the approximate temperature at matter-radiation equality.
The comoving sound horizon can be written as where η is the comoving horizon and c s = . Here for simplicity we will assume c s = 1 3 . We can then rewrite the integral for the comoving sound horizon at matterradiation equality as , where a eq is the scale factor at matter-radiation equality. The ΛCDM comoving sound horizon at matter-radiation equality can then be written as where a i → 0 and H i denote the scale factor and the Hubble scale deep in radiation domination and in the last equality we have used a i a eq . For the comoving sound horizon in the case of kination cosmology, we can divide the universe into two successive eras, a < a MK and a MK ≤ a ≤ a eq , where we assume that a eq a KR a RM . We consider the energy density of a scalar field that behaves like matter for a < a MK and kination for a > a MK in addition to the standard radiation density. The sound horizon in our kination cosmology up to matter-radiation equality can then be written as where in the last approximation we have used a i a eq and the identity a RM = a 3 MK a 2

KR
. Here √ 2H RM is the Hubble at the early matter radiation equality. Assuming a eq a KR and 3), the relative difference between the sound horizons in ΛCDM and kination cosmology at last scattering is approximated by where η ls is the comoving horizon at the surface of last scattering and we have assumed a ls 3a eq , T eq = 0.8 eV in the last equality. Thus for large enough T KR , the deviation in the angular scale of the sound horizon at last scattering can be minimal.  In order to quantify the bounds on the kination parameters T KR and T RM , we modify the publicly available CMB code, CLASS [53], to include kination cosmology. We also use Monte Python [54], a Markov Chain Monte Carlo (MCMC) code, along with CLASS to derive the posterior probability distribution on the cosmological parameters. We assume a sudden transition between the matter-dominated and kination-dominated era at T MK . For simplicity we also fix T RM = 5 keV, near the upper bound from BBN. We, however, find that the lower bound on T KR from the CMB is insensitive to T RM , as the CMB is mainly probing modes that enter the horizon at lower temperatures. We consider the following cos- During early matter domination and kination (i.e., between T RM and T KR ), matter perturbations grow linearly [32]. This rapid growth can result in an enhancement of the matter power spectrum on scales that were inside the horizon during the epochs of modified expansion. The excellent constraints provided by the CMB require that substantial modifications to the matter power spectrum must occur on scales k Fig. 4). Probes T KR . However, in order to accurately derive the constraint one needs to evolve the kination matter power spectrum into the non-linear regime and then use hydrodynamical simulations to derive the Lyman-α flux power spectrum to compare to experiments. This is beyond the scope of the present publication, but we will return to this in future work.

DARK MATTER AND BARYOGENESIS FROM AXION ROTATIONS
In this section, we discuss the production of axion dark matter and baryon asymmetry from axion rotations by the kinetic misalignment and axiogenesis mechanisms in the following subsections, respectively. We show the implications of these mechanisms for the

Axion dark matter from kinetic misalignment
Axion rotations can lead to a larger axion abundance today via the kinetic misalignment mechanism [35] than that from the conventional misalignment mechanism [55][56][57]. As long as the axion field velocity is much larger than its massθ m a , the axion continues to run over the potential barriers. If this motion continues past the time when the mass is equal to Hubble, then the axion kinetic energyθ 2 f 2 a is much larger than the maximum possible potential energy θ 2 i m 2 a f 2 a in the conventional case and thus the abundance is enhanced. Even when the axion field velocity is larger than the mass, the axion self-interactions can cause parametric resonance (PR) [58][59][60][61][62], which fragments the axion rotation into axion fluctuations [63][64][65]. The production of fluctuations by PR occurs at an effective rate 1 given by [65] On the other hand, the axion momentum k PR generated at the time of PR is of orderθ/2 due to the resonance condition. Therefore, the abundance of the axion is estimated as where the axion yield Y a is approximately conserved after PR. Here C is a factor that should be determined by numerical computation. In Ref. [35], C 2 was derived assuming the coherence of the axion rotation throughout the evolution. As noted in Ref. [47], the axion abundance is reduced by an O(1) factor in comparison with the estimation in Ref. [35] because of the extra energy of axions from non-zero momenta sourced by PR. Just after PR effectively occurs, the number-changing scattering rate of axion fluctuations is comparable to the Hubble expansion rate while axion fluctuations are over-occupied, so the number density may be further reduced by an O(1) factor, which should be determined by lattice computation; see also the discussion in [66][67][68]. We use the reference value C = 1 in this paper and demonstrate the impact of C < 1 on observations by showing results for C = 0.3. Requiring axion dark matter from the kinetic misalignment mechanism, we obtain 1 The effective rate is much smaller than the PR rate at the center of the first resonance band ∼ m 2 a /θ because of the narrow width of the band, the reduction of the axion velocity by the PR production [65], and the reduction of the axion momentum by cosmic expansion.
a prediction on T KR , The prediction for an ALP is shown by the purple dashed lines in the left panel of Fig. 5.
To avoid overproduction of axion dark matter by the kinetic misalignment mechanism, this prediction is also a lower bound on T KR . The bound can be avoided if the rotation is washed out at T < T KR . This is difficult for the QCD axion with the Standard Model because of the suppression of the wash-out rate by the small up Yukawa coupling [34,69] where we assume that the axion mass at T = T PR , is the same as the one in vacuum, m a (T PR ) = m a , and the saxion is at the minimum of the potential, S(T PR ) = f a . As with PR production from radial motion of the symmetry breaking field [68,70], the produced axions are initially relativistic. They may become cold enough to be dark matter by red-shifting, and will have residual warmness [47,64]. They become non-relativistic at temperature T NR 10 MeV f a 10 9 GeV  For sufficiently small m a and/or f a , these axions are too warm to be dark matter based on the current warmness constraint from the Lyman-α measurements [71], T NR > 5 keV. This constraint is shown by the red regions of Fig. 5.

Baryon asymmetry from axiogenesis
The observed cosmological excess of matter over antimatter can also originate from the axion rotation. The U (1) charge associated with the rotation defined in Eq. (2.3) can be transferred to the baryon asymmetry as shown in [34,45,72]. In the case of the QCD axion, the strong anomaly necessarily transfers the rotation into the quark chiral asymmetry, which is distributed into other particle-antiparticle asymmetry. More generically, the couplings of the QCD axion or an ALP with the thermal bath can transfer the rotation into particleantiparticle asymmetry. The particle-antiparticle asymmetry can be further transferred to baryon asymmetry via processes that violate the baryon number. We call this scheme, applicable to the QCD axion and ALPs, axiogenesis. To specifically refer to the QCD axion and ALPs, we use QCD axiogenesis and ALPgenesis, respectively.

Minimal axiogenesis
In the minimal scenario, which we call minimal axiogenesis, the particle-antiparticle asymmetry is reprocessed into the baryon asymmetry via the electroweak sphaleron processes.
If the QCD axion or an ALP has an electroweak anomaly, then the rotation can directly produce the baryon asymmetry by the electroweak sphaleron processes. The contribution to the yield of the baryon asymmetry is given by [34] where T ws is the temperature at which the electroweak sphaleron processes go out of equilibrium and is approximately 130 GeV in the Standard Model [73], and c B is a modeldependent coefficient given in Ref. [72] that parameterizes the anomaly coefficients and the axion-fermion couplings. When the transfer is dominated by axion-gauge boson couplings, c B is typically O(0.1), while if dominated by axion-fermion couplings, it can be much smaller.
For the QCD axion, to produce sufficient Y B and to avoid overproduction of dark matter by kinetic misalignment requires that f a 10 7 GeVc B /C, which is disfavored by astrophysical constraints [74][75][76][77][78][79][80] unless c B /C > 10. The baryon asymmetry can be enhanced if the electroweak phase transition occurs at a higher temperature, with both dark matter and baryon asymmetry of the universe explained by the rotation of the QCD axion when the inequality is saturated.
For an ALP, we may choose sufficiently small m a to avoid the over production without modifying the electroweak phase transition temperature. Requiring that the baryon yield of Eq. (4.6) match the observed baryon asymmetry gives a constraint onθ(T ws ); using Eq. (2.10) this can be converted to a prediction for T KR T KR = 3.5 TeV Note that this is necessarily the case when T MK > T ws . For lower T MK , S(T ws ) > f a is possible. Sinceθ m S when S > f a andθ ∝ T 3 after S = f a , we haveθ(T ws ) m S and therefore the saxion mass is predicted to be The bound is saturated when S(T ws ) > f a , which is the case if T MK < T ws .
Both dark matter and baryon asymmetry of the universe is explained by the axion rotation, which is called ALP cogenesis [72], when This prediction is shown by the green lines in Fig. 5.

B − L number violation by new physics
In the presence of an operator that violates lepton number and generates Majorana neutrino masses, the transfer of asymmetries can be more efficient. The operator creates a non-zero B − L number, which is preserved by Standard Model electroweak sphaleron processes. Since the production of B − L at high temperatures depends on whether the lepton number violating interactions are in equilibrium, the determination of the final baryon number in this scenario is sensitive to the full cosmological evolution. As an example, for the models studied in Ref. [45], the baryon asymmetry is given by  [45]. Regardless, the saxion mass is generically predicted to be O(30 − 10 4 ) TeV × (0.1/c B ) by this baryogenesis mechanism, named leptoaxiogenesis generically, QCD lepto-axiogenesis for the QCD axion, and lepto-ALPgenesis for the ALP. While m S 30 TeV appears difficult based on Eq. (4.11), the case of TeV scale supersymmetry is possible in some special cases presented in Ref. [45], using a thermal potential.
Other axiogenesis scenarios are are also considered in the literature [81,82]. Baryon asymmetry may be dominantly produced at a temperature T dec with where Γ B is the transfer rate of the axion rotation into baryon asymmetry. For models with T dec > T ws , the lower bound on m S is generically stronger than that for ALPgenesis.

Implications to axion kination parameters
In Fig

GRAVITATIONAL WAVES
In this section, we discuss how the spectrum of primordial gravitational waves is modified by eras of matter and kination domination generated from axion rotation, as discussed in Sec. 2. We consider gravitational waves created by quantum fluctuations during inflation and by local cosmic strings. In both production mechanisms, the spectrum is nearly flat in the standard cosmology with radiation domination. As we will see, the evolution of a universe with successive eras dominated by radiation, matter, kination, and back to radiation induces a triangular peak in the gravitational wave spectrum that can provide a unique signal for axion rotation and kination cosmology.

From inflation
We first discuss the primordial gravitational waves produced from quantum fluctuations during inflation [83]. In the standard cosmology with radiation domination, the spectrum is nearly flat for the following reason. After inflation, a given mode k is frozen outside the horizon, k < H. As the mode reenters the horizon, k > H, it begins to oscillate and behaves as radiation. The energy density of the mode at that point ∼ k 2 h 2 (k)M 2 Pl , where h is the metric perturbation, whose spectrum is almost flat for slow-roll inflation. Since k ∼ H at the beginning of the oscillation, the energy density of the gravitational waves normalized by the radiation energy density ∼ H 2 M 2 Pl is nearly independent of k up to a correction by the degree of freedom of the thermal bath 2 During matter or kination domination in our scenario, the energy density at the horizon crossing is still H 2 h 2 M 2 Pl , but the radiation energy density is now much smaller than H 2 M 2 Pl . As a result, the energy density of gravitational waves with a mode k is inversely proportional to the fraction of the radiation energy density to the total energy density when the mode enters the horizon. This means that the spectrum should feature a triangular peak in axion kination. The modes that enter the horizon at T > T RM are not affected and remain flat.
(See, however, a comment below). Therefore, as the horizon-crossing temperature decreases below T RM , the gravitational waves are enhanced and reach the maximal value at T MK , where the fraction of the radiation energy is minimized. The gravitational wave strength decreases again for the horizon crossing temperatures below T MK , and returns to a flat spectrum below T KR . Note that gravitational waves that enter the horizon during matter domination are also enhanced because of the absence of entropy production after matter domination. We use an analytical approximation where each mode begins oscillations suddenly at the horizon crossing. Also approximating the evolution of H by a piecewise function with kinks at the three transition, the resultant spectrum of gravitational waves is given by where V inf is the potential energy during inflation. We normalize the spectrum to match a full numerical result at f < f KR and f > f RM , which is found to be consistent with the result of the numerical computation in Ref. [85].
We comment on possible further modification of the spectrum. Axion kination relies on a nearly quadratic saxion potential, which is natural in supersymmetric theories. The degrees of freedom of the thermal bath change by about a factor of two across the superpartner mass threshold, suppressing the gravitational wave signals by a few tens of a percent at high frequency [121]. This depends on the superpartner masses, and we do not include this effect for simplicity. We also assume that the radial mode of P does not dominate the energy density of the universe. If it does, entropy is created by the thermalization of the radial mode and gravitational waves at f > f RM can be suppressed. If the initial rotation before thermalization is highly elliptical, after the thermalization the universe is radiation dominated for a long time because of the radial mode energy much larger than the angular mode energy, so the suppression occurs at f f RM . If the initial rotation is close to a circular one, the universe is radiation dominated only for a short period, so the suppression occurs right above f RM . In principle, we can learn about the very UV dynamics of axion rotations through the observations of gravitational waves.
In Fig. 8  and r = 0.001 near the sensitivity limit of future CMB observations [122], respectively. We show the spectrum for several sets of (T KR , T RM ) in different colored curves.
The spectrum shown in the solid (dotted) orange curve corresponds to the value of T KR predicted from QCD axion dark matter via the kinetic misalignment mechanism with C = 1 (C = 0.3) according to Eq. (4.3), with the maximal T RM allowed by the constraints shown in is close to the maximum, CE can also observe the signal, but the inflation scale must be almost at the present upper bound as seen in Fig. 8.
The parameter space for ALPgenesis is shown in Fig. 10   The parameter space for lepto-ALPgenesis is shown in Fig. 11. Here m S is fixed so that the observed baryon asymmetry is explained by lepto-ALPgenesis; see Eq. (4.11). f a is then fixed by Eqs. (2.9) and (2.10). The meaning of shaded regions and contours are the same as in Fig. 10.

From cosmic strings
We next discuss gravitational waves emitted from local cosmic strings [123]. Local cosmic strings are topological defects produced upon gauge symmetry breaking in the early universe, such as U (1) symmetry breaking [124]. The breaking of a local U (1) symmetry, and hence formation of a cosmic string network, arises in many theories beyond the Standard Model.
For example, one of the best motivated cases is U (1) B−L , which is the unique U (1) symmetry that does not have a mixed anomaly with the Standard Model gauge symmetry. Moreover, U (1) B−L can be embedded into SO(10) together with the Standard Model gauge group, and whose spontaneous symmetry breaking can provide the right-handed neutrino masses in the see-saw mechanism [125][126][127][128].
After production, the cosmic string network follows a scaling law with approximately O(1) long strings per Hubble volume which is maintained from the balance between conformal expansion with the universe and losses from self-intercommutation. The self-intercommutation byproducts of the long string network lead to the formation of a network of string loops with a new loop forming nearly every Hubble time and with a loop size proportional to the horizon [129]. These subhorizon loops oscillate and redshift like matter before decaying from the emission of gravitational waves. Because of the scaling law of the string network, the energy density fraction of the cosmic strings is nearly independent of temperature, and the spectrum of gravitational waves emitted from the local cosmic strings during radiation domination is nearly flat.
During kination or matter domination by axion rotation, the size of the horizon for a given temperature is smaller than it would be in a radiation dominated universe. This enhances the energy density of strings relative to the radiation density, and the spectrum of gravitational waves feature a triangular peak in our axion kination cosmology. Since the production of gravitational waves involves two steps that occur at widely separated timesthe production of string loops and their later decay-the computation is more involved than the inflation case of Sec. 5.1.
The present day gravitational wave spectrum from a stochastic background of cosmic string loops is [12,123] Here Gµ 2 P m = ΓGµ 2 m −q /ζ(q) is the power radiated by the mth mode of an oscillating string loop with Γ 50 being a constant determined from the average power over many types of string loop configurations [12,129,130]. The power index q is 4 3 , 5 3 , or 2, if the gravitational power is dominated by cusps, kinks, or kink-kink collisions, respectively [12,123,131]. We will take q = 4/3, but for now we keep it as a free parameter. The present day critical density is ρ c , and the factor C m is given by a(t) + ΓGµt (α + ΓGµ) −1 denotes the formation time of a string loop of length l that emits gravitational waves at frequency f = 2m/l at time t . The lower integration time, t scl , is the time the infinite string network reaches scaling. F ≈ 0.1 [132] characterizes the fraction of energy that is transferred by the infinite string network into loops of size l k = αt k 3 , and C eff characterizes the loop formation efficiency which depends on the equation of state of the universe at loop formation time t k . C eff can be estimated from the velocity-one-scale model of the infinite string network and is found to be [8,134] C eff (t k ) ≈ Radiation Loop Decay Era oscillation is then [9] Ω (1) where p = −3+2/(1+w 1 )+4/(3(1+w 2 )) characterizes whether the integral (5.5) is dominated at t Γ (p < 0) or the latest possible emission time t in that cosmological era (p ≥ 0) [9]. 4 The frequency dependence of Ω  Table I. In the modified cosmology under consideration, the universe transitions from being dominated by radiation to matter at T RM , to kination at T MK , and back to radiation upon merging with the standard cosmology at T KR . From Table I, we may therefore expect for sufficiently long eras of radiation, matter, and kination that Ω That is, a triangular shaped peak in spectrum.
Although the first mode dominates the total power emitted by a string loop, the sum of the contributions from all higher modes can appreciably change this power dependence of the total spectrum. The effect of the higher modes can be analytically estimated by noting that Ω GW (f /m) [42]. For example, assuming that Ω GW is a broken power law 4 In a standard radiation dominated era, Ω GW is dominantly sourced by the smallest loops in the horizon due to their greater population and the independence of gravitational wave power, ΓGµ 2 , on loop size.
The smallest loops are those about to decay and hence for the standard cosmology, p < 0. However, this is not the case in more general cosmologies.
proportional to f α for f < f 0 and f β for f ≥ f 0 , we may write the total spectrum as where the peak amplitude of the triangular spectrum, Ω GW h 2 | peak , and the peak frequency f peak , can be thought of as Ω GW h 2 | MK and f MK , since the peak is associated with loops formed at T MK . The frequency of the peak, Eq. (5.12), is set by the invariant size of the loop at the formation time t MK with the emission frequency at decay 2/l(t MK ) redshifted to the present. Similarly, from Table I, we can see that the loops that form at the transition from matter to late era radiation, t KR , are responsible for the amplitude of the lower left vertex of the axion kination triangle. Again, the frequency of these loops is the emission frequency at decay, 2/l(t KR ) redshifted to the present as given by Eq. (5.13). Last, note that the fiducial values of T KR = 100 MeV and T MK = 2 GeV, correspond to T RM ≈ 100 GeV, which corresponds to the dark purple curve of Fig. 13.
In general, for brief eras of kination and matter domination, the gravitational wave spectrum will not reach its asymptotic dependence, Ω tot nor will the kination era peak be sharply defined. Consequently, we numerically evaluate Eq. (5.4) to precisely determine Ω GW over a wide range of {T KR , T RM , T MK }. In doing so, we numerically compute the time evolution of the scale factor from the Friedmann equation + Ω k,θ a KR a(t) 6 + Ω m,θ a MK a(t) 3 1 2 (5.14) where Ω k,θ = Ω r a(t 0 ) a KR The left panel of Fig. 12 shows the imprint of the saxion potential on the stochastic string gravitational wave background. The black curve corresponds to the piecewise approximation of the ρ θ contribution to the Hubble rate as used in Eq. (5.14). The blue dotted and orange solid lines show the spectrum for the two-field model and the log potential, respectively.
Similar to the gravitational wave spectrum of Fig. 7 for inflation, the spectrum for the twofield model is close to that for the piecewise approximation, while that for the log potential deviates from them. In what follows, we use the piecewise approximation of ρ θ .
The right panel of Fig. 12 illustrates two key features that axion kination imparts to the stochastic string gravitational wave background. First, the purple curves show the m = 1 contribution to the spectrum, Ω GW , while the red curves shows Ω GW after summing over 10 4 harmonics. For the m = 1 amplitude, the triangular shaped peak approaches The sum over high modes partially flattens the right side of the kination induced peak, shifting the spectral dependence from f −1 to f −1/3 . We fix (T KR , T RM ) = (1 GeV, 10 TeV). In both panels, the second, smaller triangle at high frequencies is an additional fingerprint of axion kination and arises from loops that form in the early radiation dominated era and decay in the subsequent matter or kination dominated eras (see Table I). Both panels assume Gµ = 5 × 10 −15 , and α = 0.1.
the expected f −1 rise and f 1 fall as shown in Table I. The amplitude summed over 10 4 modes, however, demonstrates how the total amplitude deviates from the m = 1 amplitude.
Summing over higher harmonics increases the amplitude roughly by a factor of ζ(4/3), and most importantly, the contribution from the higher harmonics changes the f −1 tail on the right side of the kination induced triangle into a much shallower f −1/3 tail while leaving the f 1 decay on the left side the triangle the same. Such a long and shallow UV tail allows high frequency gravitational wave detectors to discern axion kination from the standard ΛCDM spectrum even when the triangular kination peak is at much lower frequencies. In addition, a second key feature of axion kination is shown in the second, smaller triangle at higher frequencies compared to the main triangle. The second triangular bump in the spectrum arises from loops that form in the early radiation dominated era prior to kination, and decay in the matter, kination, or late radiation era. As seen from Table I, loops that form in the early radiation era and decay in these other eras are expected to exhibit a shallower rise and fall in amplitude akin to the main triangular shaped enhancement from loops that form in the matter or kination eras. For sufficiently short eras of kination and matter domination, the smaller, second bump is visible even after summing over higher harmonics. For long eras of kination, the sum over higher harmonics can merge the main kination induced peak with this smaller second peak, as shown for instance, in the purple curves in the left panel of Fig. 13. Nevertheless, the slightly broken power law near 10 3 Hz for the solid purple and 10 5 Hz for the lighter purple contour is a remnant left over from this second triangular peak.
The observation of a broken decaying power law or the second triangular bump itself may provide a unique gravitational wave fingerprint for axion induced kination.  Similarly, the right panel of Fig. 13 shows the modified gravitational wave spectrum for late kination cosmologies. Here we show a collection of spectra that pass through the observed NANOGrav signal [40] that are consistent with CMB and BBN constraints. The dashed black contour shows Ω GW h 2 for the standard ΛCDM cosmology [41,42] while the gray contour (Gµ = 3.5 × 10 −11 ) shows Ω GW h 2 for the maximum allowed T RM (6 keV) and near the minimum allowed T KR (130 eV), producing the largest kination peak passing through NANOGrav that is consistent with CMB and BBN. As T KR and T RM converge and the kination era decreases in duration, Ω GW h 2 converges with the standard result, shown, for example, by the magenta contour (Gµ = 4.0 × 10 −11 ). Fig. 13 demonstrates that a striking difference can exist between Ω GW h 2 in the axion kination cosmologies and the standard ΛCDM cosmology.
To understand the connection between the experimental detection of axion kination via string gravitational waves and axion kination parameters, we first reduce the four dimen- For a given (T KR , T MK ), we register a detection of axion kination in a similar manner to the "turning-point" prescription of [14]: First, Ω GW h 2 must be greater than the threshold for detection in a given experiment. Second, to actually distinguish between Ω GW h 2 in the axion cosmology and the ΛCDM cosmology, we require that their percent relative difference be greater than a certain threshold within the frequency domain of the experiment. Following [14], we take this threshold at a realistic 10% and a more conservative 100% relative and 0.1, respectively, to fit the NANOGrav data [40]. For a given (T KR , T MK ), a detection is registered when the difference in amplitudes, Ω GW − Ω GW,0 is greater than 10% (solid) or 100% (dashed) of the standard cosmological amplitude, Ω GW,0 , within the sensitivity curve the detector.
In the transparent shared regions, the peak of the spectrum originated from axion kination can be detected.
to detect late axion kination cosmology. Future detectors like SKA can probe the nanohertz triangular bump. For sufficiently low T KR and high T RM , a kination signal may already be observable or excluded at NANOGrav.
Early axion kination is consistent with axiogenesis above the electroweak scale, and can be probed by laser interferometers. We show the constraints on the parameter space of minimal ALPgenesis together with the detection prospects in the upper panels of Fig. 15. Hz, thereby evading detection. Still, a good portion of the parameter space with f a 10 8 GeV can imprint signals that are detectable by future observations. Future gravitational wave detectors that can observe super-kilohertz frequencies can potentially probe earlier eras of axion kination and hence larger f a . In the transparent shaded region, the peak of the spectrum produced by axion kination can be detected. As we argued, this is a smoking-gun signature of axion kination and the detailed shape of the peak contains information about the shape of the potential of the complex field that breaks the U (1) symmetry. The lower two panels of Fig. 15 show the constraints and prospects for lepto-ALPgenesis, for values of m S used in Fig. 11. Future laser interferometers can probe much of the parameter region with low f a 10 8 GeV.

SUMMARY AND DISCUSSION
Axion fields, due to their lightness, may have rich dynamics in the early universe. In this paper, we considered rotations of an axion in field space that naturally provide kination domination preceded by matter domination, which we call axion kination. This non-standard evolution affects the spectrum of possible gravitational waves produced in the early universe.
To be concrete, we investigated gravitational waves from inflation and local cosmic strings, which have a nearly flat spectrum when they begin oscillations or are produced during radiation domination. We found that kination domination preceded by matter domination induces a triangular peak in the gravitational wave spectrum.
We studied the theory for axion kination, which involves an approximately quadratic potential for the radial mode and has three parameters: the mass of the radial mode, the axion decay constant, and the comoving charge density. We derived constraints on this parameter space from successful thermalization of the radial mode, BBN, and the CMB. We found large areas of fully realistic parameter space where the theory yields axion kination.
The allowed region splits into two pieces, one having early kination domination before BBN and the other having late kination after BBN but well before the CMB last scattering.
Introducing a mass for the axion, we found that part of the axion kination parameter space is consistent with axion dark matter by the kinetic misalignment mechanism while part is not, due to the warmness constraint on dark matter. Similarly, we showed that part of the axion kination parameter space is consistent with generating the baryon asymmetry by ALPgenesis. Furthermore, there are constrained regions with ALP cogenesis yielding both dark matter and the baryon asymmetry, and also regions with the baryon asymmetry successfully generated by lepto-ALPgenesis.
As demonstrated in Sec. 5, axion kination modifies the spectrum of possible primordial gravitational waves through the modification of the expansion history of the universe. By analyzing the spectrum, we can in principle determine the product of the radial mode mass If the inflation scale is not much below the current upper bound, future observation of gravitational waves can detect the spectrum modified by axion kination, or even the peak of the spectrum that contains information on the shape of the potential of the U (1) symmetry breaking field. In particular, if the QCD axion accounts for dark matter via the kinetic misalignment mechanism, a modification of the gravitational wave spectrum is predicted at high frequencies, f 10 −2 Hz, as shown in Fig. 8. In this case it is very interesting that this gravitational wave signal can be detected by DECIGO and BBO over most of the allowed parameter space, as shown in Fig. 9; in a significant fraction of the parameter space the gravitational wave peak will be probed. Furthermore, a signal may also be seen at CE if the inflation scale is very near the current upper bound or the sensitivity of CE is improved.
For gravitational waves from cosmic strings, for fixed axion kination cosmology parameters, the modification of the spectrum is predicted at higher frequency, so the QCD axion will not affect the spectrum observable by near future planned experiments. ALPs can affect the spectrum in an observable frequency range.
Gravitational waves from cosmic strings provide signals that can probe axion kination over a wide range of Gµ, T RM and T KR , as illustrated in Fig. 13. We examined cosmic strings with a tension suggested by NANOGrav in detail. If axion kination occurs before BBN, the NANOGrav signal can be fitted by the same cosmic strings parameters as in standard cosmology. Importantly, axion kination enhances the spectrum at higher frequencies, allowing laser interferometers to probe the kination era. The enhancement can occur in the parameter region consistent with axiogenesis scenarios, as shown in Fig. 15. If axion kination occurs after BBN, the NANOGrav signal is fitted by a smaller string tension, as shown in the left panel of Fig. 14, and a detailed examination of the spectrum will determine if axion kination is involved. The spectrum at higher frequencies is suppressed, which can be detected by laser interferometers in the parameter region shown in the right panel of Our kination era is preceded by an epoch of matter domination that ends without creating entropy. Therefore, matter and kination domination can occur even after BBN. This allows for enhancements to the matter spectrum on small scales that may be probed by observations of Lyman-α and 21 cm lines. Evolving the enhanced matter power spectrum into the non-linear regime and understanding its effects on the Lyman-α flux spectrum as well as hierarchical galaxy formation, and constraints arising from corresponding observations will be discussed in future work.
In this paper, we concentrated on gravitational waves produced by inflation or local cosmic strings and modified by axion kination. At any temperature with early matter or kination domination, the Hubble scale is larger than with radiation domination, and hence, quite generally, primordial gravitational waves are enhanced by axion kination. Furthermore, a distinctive feature appears in the spectrum, a peak or bump depending on the field potential, containing information that probes in detail the era of kination and its origin. It will be interesting to investigate other sources of primordial gravitational waves.
Note added. Soon after the current manuscript was announced on the arXiv, Ref. [135] appeared as well. While our analyses have been conducted independently, Ref. [135] also discovered the triangular peak signature of axion kination in the spectrum of the primordial gravitational waves from inflation. where W is the Lambert W function and S i is an initial field value at a scale factor of a i .
When S i f a and a is not much above a i , S 2 ∝ a −3 andθ is nearly a constant. For a a i , S f a andθ ∝ a −3 .
The dependence of the energy density so at early times, S f a , the rotation behaves as matter ρ θ ∝ a −3 , while at late times, S f a , the rotation behaves as kination ρ θ ∝ a −6 . This behavior is seen in the orange curve of Fig. 1.

Two-field model
We next consider the two-field model in Eq. (2.5). We assume that the saxion field value is much larger than the soft masses, so that we may integrate out a linear combination of P andP that is paired with X and obtain a mass ∼ S. Using the constraintP = v 2 P /P , from the kinetic and mass terms of P andP , we obtain an effective Lagrangian The potential has a minimum at |P | = √ r P v P when both m 2 P and m 2P are positive. PuttingS =Ṡ = 0, the equation of motion of S requires thaṫ The equation of motion of θ gives a conservation of the angular momentum, Without loss of generality, we assume that P v P , i.e., S v P initially. We first consider r P > 1 (m P < mP ). Eq. (A.7) meanṡ When S v P , we obtainθ 2 m 2 P . As S approaches the minimum from which we again observe that at early times, S √ 2v P , the rotation behaves as matter ρ θ ∝ a −3 , whereas at late times, S √ 2v P , the rotation behaves as kination ρ θ ∝ a −6 . This evolution is illustrated by the blue curves of Fig. 1 for various values of r P . We next consider r P = 1. Eq. (A.7) has two solutions, 1)θ 2 = m 2 P with S = √ 2v P and 2) S = √ 2v P with unrestrictedθ. For S v P , the solution is in the branch 1) and gives matter scaling. Charge conservation implies S 2 (1 + 4v 4 P /S 4 ) ∝ a −3 . As S decreases according to this scaling and reaches √ 2v P , the branch 2) should be used and charge conservation implieṡ θ ∝ a −3 , giving kination scaling. Finally, consider r P < 1. When S v P , we again obtainθ 2 m 2 P . However, as S approaches √ 2v P (before reaching the minimum at √ 2r P v P ),θ 2 derived from Eq. (A.9) diverges. This is the point at which the solution becomes unstable. Indeed, when r P < 1, for a fixed charge, it is energetically favored to have rotations inP rather than in P , so the rotation dominantly in P is at the most a meta-stable solution. When S reaches √ 2v P , the solution becomes unstable. Quantum tunneling may occur before the instability is reached.
We leave the investigation of this scenario to future work and assume r P ≥ 1.
logical input parameters for ΛCDM and the axion kination model: baryon density Ω b h 2 , DM density Ω c h 2 , spectral index n s , primordial amplitude A s , optical depth at reionization τ and the temperature at kination radiation equality T KR . We fix N eff = 3.046 and assume log flat prior for T KR between 1 eV ≤ T KR ≤ 5 keV. We fix T RM = 5 keV and use the identity T 3 MK = T 2 KR T RM , obtaining T KR > 130 eV at 95%.
For completeness, Fig. 17 shows the derived parameters: the effective redshift at reionization z reio , dark energy density Ω Λ , Hubble expansion rate today H 0 , and matter fluctuation amplitude σ 8 .