Cosmological window onto the string axiverse and the supersymmetry breaking scale

In the simplest picture, the masses of string axions populating the axiverse depend on two parameters: the supersymmetry breaking scale $M_{\rm susy}$ and the action $S$ of the string instantons responsible for breaking the axion shift symmetry. In this work, we explore whether cosmological data can be used to probe these two parameters. Adopting string-inspired flat priors on $\log_{10}M_{\rm susy}$ and $S$, and imposing that $M_{\rm susy}$ be sub-Planckian, we find $S=198\pm28$. These bounds suggest that cosmological data complemented with string-inspired priors select a quite narrow axion mass range within the axiverse, $\log_{10}(m_a/{\rm eV}) = -21.5^{+1.3}_{-2.3}$. We find that $M_{\rm susy}$ remains unconstrained due to a fundamental parameter degeneracy with $S$. We explore the significant impact of other choices of priors on the results, and we comment on similar findings in recent previous literature.


I. INTRODUCTION
The theory of QCD allows for a currently unobserved CP-violating interaction [1][2][3][4]. A plausible solution to the so-called strong-CP problem [5,6] predicts the existence of the QCD axion a ¼ aðxÞ [7,8], a pseudoscalar Goldstone boson that couples to the number density of QCD instantons via Here, the trace is taken in the three-dimensional representation of SU (3), and f a is the spontaneous symmetrybreaking scale (or axion decay scale). In viable "invisible" axion models, stellar object cooling considerations provide a bound f a ≳ 10 9 GeV [9][10][11][12], with the exact value of the bound depending on the axion model considered. The shift symmetry a → a þ const, which holds at the classical level, is explicitly broken by the same QCD instanton effects which are also responsible for generating the axion potential. Owing to these nonperturbative effects, the QCD axion acquires a mass m a ¼ Λ 2 QCD =f a , where Λ QCD ≃ 75.5 MeV [7,8,13]. It is well known that the axion can constitute the cold dark matter (CDM) [14][15][16], with the axion decay constant being as large as the grand unification theory (GUT) energy scale in the so-called anthropic axion window [17][18][19][20][21][22][23][24][25][26][27][28][29][30][31].
Along with the QCD axion, other "axionlike" particles (ALPs) arise either from the breaking of "accidental" symmetries [32][33][34][35][36][37][38] or from manifold compactification within string theory [27,[39][40][41][42][43][44][45][46][47][48]. Both these scenarios feature a symmetry-breaking scale Λ a and an ALP decay constant f a , with the axion field acquiring the mass m a ¼ Λ 2 a =f a . However, unlike the case for the QCD axion, the ALP energy scale Λ a is not tied to the QCD energy scale, so the mass m a and the decay scale f a can effectively be treated as two independent parameters. Although somewhat fundamentally less motivated than the QCD axion, ALPs are potentially suitable dark matter candidates [49][50][51][52]. Of particular interest are ultralight axions (ULAs) [53][54][55][56][57][58][59][60][61][62][63][64], the mass of which resides in the range m a ∈ ½10 −27 ; 10 −18 eV. These ULAs manifest their wavelike behavior by suppressing the matter power spectrum at astrophysical scales. It has been argued that this suppression of power could be the key to addressing a number of controversies (e.g., the "missing satellites" and the "cuspcore" problems) arising in the standard ΛCDM cosmology on galactic and subgalactic scales (see Ref. [65] for a review). Hereafter, we shall refer to axionlike particles interchangeably as "ALPs" or "axions" when we are not interested in their origin (e.g., from the string axiverse), as opposed to "string axions" which originate specifically from the string axiverse (to be discussed later).
It has long been noticed that axions arise naturally within string theory compactifications as Kaluza-Klein zero modes of antisymmetric tensor fields [39,41]. Zero modes originate from noncontractable cycles on the compactified manifold. The number of zero modes is fixed by the topology of the compactification manifold itself and is generally in the order of hundreds (for instance, for compactifications on Calabi-Yau manifolds, the number of zero modes is given by the Hodge number of the manifold itself). Notice that this assumes that the size of the extra dimensions is finite, unlike, e.g., the case of Randall-Sundrum models [66][67][68][69]. A fraction of these zero modes is expected to acquire a mass through nonperturbative string instanton effects [70], which can be characterized by the (dimensionless) action of the instantons S (which scales with the volume of the corresponding cycles) [41], as well as a nonperturbative UV cutoff scale μ. This nonperturbative scale is related in turn to the supersymmetry (SUSY-) breaking scale M susy [41,42] (with SUSY almost inevitably appearing in any realistic string theory) as μ ∝ ffiffiffiffiffiffiffiffiffiffiffi M susy p , if the axion potential arises from the superpotential generated by string instantons (which is often the case). This suggests that a reasonable way to characterize string axions is by exploring the μ-S (or equivalently M susy -S) parameter space.
Away from the swampland (Ref. [71]; see Refs.  for recent developments), the landscape of string vacua [40,41] give rise to a plethora of light axions, known in the literature as the string axiverse [42], along with various other massless modes [100]. Since for the QCD axion the θ parameter is constrained to be smaller than 10 −10 , string corrections are negligible compared to those given by the QCD potential. On the contrary, the mass of lighter axions is primarily fixed by the nonperturbative string contributions. The exploration of the string axiverse scenario, wherein a multitude of axions populates various orders of magnitude in masses, is being pursued by a variety of searches, including a rotation in the polarization of the cosmic microwave background (CMB) radiation spectrum (for masses ranging from 10 −33 to 4 × 10 −28 eV), a suppression in the power spectrum of density perturbations at small scales (for masses ranging from 10 −28 to 10 −18 eV), the altered dynamics of rotating black holes due to the effect of superradiance (for masses ranging from 10 −22 to 10 −10 eV) [42], and various laboratory searches (for masses larger than about 10 −15 eV) like ABRACADABRA [101], ADMX [102], KLASH [103], QUAX [104], X3 [105], CULTASK [106], MADMAX [107], ARIADNE [108], IAXO [109], and CASPEr [110]; see, e.g., Ref. [111] for a review. Axions lighter than the present Hubble rate H 0 ∼ 10 −33 eV are still frozen today and do not contribute to the present matter content of the Universe.
In the present work, it is our aim to address the following question: what can cosmology tell us about the string axiverse and its parameters? As discussed previously, we address the question by focusing on the nonperturbative scale μ (or equivalently the SUSY-breaking scale M susy ) and the dimensionless symmetry-breaking instanton action S as parameters characterizing the string axiverse. Focusing for definiteness on the case in which the axions are present during inflation, we will also consider the initial misalignment angle θ i and the primordial isocurvature fraction β as additional parameters. For an incomplete list of other works examining the axiverse, and especially its cosmology, see, e.g., Refs. [112][113][114][115][116][117][118][119][120][121][122][123][124][125][126]. We do not necessarily assume that the axion is the totality of the dark matter, so we assume a mixed dark matter scenario [127][128][129][130][131].
We therefore characterize the string axiverse by the fourdimensional parameter space spanned by the parameters ðμ; S; θ i ; βÞ [alternatively ðM susy ; S; θ i ; βÞ] and explore how these parameters can be constrained by cosmological data. On this matter, a caveat/warning is in order at this point. As it is, to the best of our knowledge, the first time an attempt to constrain the parameters M susy and S is made, our goal is not to provide a full-fledged analysis utilizing all available cosmological data (e.g., the full CMB temperature and polarization anisotropy power spectra or measurements of galaxy power spectra) but rather to get a feel for whether cosmology can actually provide information on M susy and S and, if so, which region of parameter space is selected and what the physical motivation is for such a region being chosen. To this end, we include the following information/ requirements from cosmology: (i) constraints on the tensor-to-scalar ratio r from the Planck satellite; (ii) constraints on the primordial isocurvature fraction from the Planck satellite; (iii) the requirement that the energy density in axions should not exceed the dark matter energy density measured by Planck. We obtain constraints on the model parameters ðM susy ; S; θ i ; βÞ by performing Bayesian parameter inference in light of the aforementioned cosmological data, with string-motivated priors for the model parameters (flat in log 10 M susy and S), which we discuss in detail in Sec. IV. In this first analysis, we find S ¼ 198 AE 28 at 68% C.L.. On the other hand, we find μ and hence M susy to be poorly constrained, due to a fundamental parameter degeneracy with S. To break this degeneracy, we perform for purely illustrative purposes two additional analyses in which we fix the instanton action to S ¼ 198 and S ¼ 153, corresponding to the central value and the 2σ lower bound on S, respectively. When fixing S ¼ 198, we find that the SUSYbreaking scale is of the order of M susy ∼ 10 8 TeV, whereas for S ¼ 153, we find that M susy ≲ 10 4 TeV at 95% C.L., which could lead to potentially interesting signatures at the proposed 100 TeV collider [132].
Within the string axiverse, theoretical results make it possible to link the quantities M susy and S to the more familiar axion mass and decay constant, m a and f a [39][40][41]. We exploit these results to convert our previous constraints on M susy and S to constraints on m a and f a for the string axiverse case (hence adopting string-inspired priors on the parameters). Perhaps surprisingly, we find that the string axiverse results select a rather tight range of axion masses, around m a ≈ 10 −22 eV, thus favoring the ULA interpretation. Interestingly, this appears to be in agreement with theoretical work which suggests that most of the axions originating from the string axiverse would be of the ULA type [42].
The previous result, selecting a very specific range for m a in the string axiverse case, raises the following question: what if we were to repeat the analysis without focusing on string axions, hence without adopting string-inspired priors for μ and S? We address this question by performing a fourth and final analysis in which we consider the parameter space spanned by the parameters ðm a ; f a ; θ i ; βÞ, as per previous discussions, with phenomenology-inspired priors on m a and f a (flat in log 10 m a and log 10 f a ). In this case, we find that our analysis selects a rather broad region in the m a parameter space, with log 10 ðm a =eVÞ ≈ −11 þ11 −6 at 68% C.L. and ≈ − 11 þ14 −15 at 95% C.L., slightly preferring the heavier mass region, albeit at a very mild significance. This suggests that the choice of prior distributions (string inspired vs not string inspired) plays a relevant role in our results, a fact worth keeping in mind when reading our paper. Our main results are shown in Fig. 1, in which we plot the posterior probability distributions we obtain for the string axion (first analysis, solid black line) and the ALP (fourth analysis, dashed black line) cases, with detection techniques for axions of various mass ranges overlain on the plot. The results of the four analyses we perform are briefly summarized in Table I. In principle, additional readily available datasets can be added to the analysis performed, to further constrain different parts of the axion parameter space range, including: (i) limits on the axion-photon coupling from direct searches for the axion in the lab through lightshining-through-wall experiments, like OSQAR [133,134] and ALPS [135,136]; (ii) cavity searches [137,138], the sensitivity of which depends on the local density of axion dark matter, like ADMX [102,139] and YWL [105,140]; (iii) axion helioscopes [141], like CAST [142,143]; (iv) astrophysical constraints on the axion-photon coupling from the supernova 1987A [144]; (v) constraints on the axion-electron coupling from the cooling of white dwarves [145,146]; (vi) constraints on the couplings of the axion to electrons and photons from consideration on the branch of red giant stars [147]; (vii) searches for gravitational waves in relation to black hole superradiance [42]; (viii) constraints on the mass of ultralight bosons in relation to the Lyman-α forest [148][149][150]; (ix) the matter power spectrum [151,152]. A list of experiments and astrophysical observations has been recently used to constrain the parameter space of the QCD axion in Ref. [153], while precision cosmological data have been used to explore the axiverse in Refs. [151,152]. We discuss these findings in comparison with our results in Sec. VI. This paper is organized as follows. In Sec. II, we first review some basics of axion cosmology. We review the FIG. 1. Marginalized posterior probability distributions for log 10 ðm a =eVÞ (with m a the axion mass) normalized to their maximum values, for two of the analyses considered in our work. In the first analysis (solid black curve), we take string-inspired flat priors on log 10 M susy and S, with M susy the supersymmetrybreaking scale and S the dimensionless action of the instantons responsible for breaking the axion shift symmetry and generating the axion mass; see Sec. III. In the latter analysis (dashed black curve), we assume phenomenology-inspired flat priors on log 10 m a and log 10 f a , with m a the ALP mass and f a the axion decay scale. In both cases, we apply the bounds on the parameter space described in Sec. II. Hashed areas represent the various detection techniques for different mass ranges of the string axiverse.
TABLE I. Summary of the parameters used and inferred values for log 10 ðm a =eVÞ and log 10 ðM susy =TeVÞ obtained for each analysis. Intervals of the form μ AE σ are 68% C.L. intervals, whereas quoted upper limits are 95% C.L. upper bounds. We consider a parameter to be unconstrained when its 95% C.L. interval is almost as wide as the parameter prior.

Analysis
Parameters log 10 ðm a =eVÞ log 10 ðM susy =TeVÞ Unconstrained string theory axion and the string axiverse in Sec. III. In Sec. IV, we discuss in more detail the cosmological data used and the statistical analysis performed. In Sec. V, we discuss the constraints we obtain on the relevant string parameters, M susy and S. We then compare the constraints obtained on m a (axion mass) and f a (axion decay constant) in the string axion case (by converting the bounds we obtained on M susy and S in the first analysis), as well as in the ALP case (fourth analysis). We also discuss the results obtained by fixing S for purely illustrative purposes, in order to resolve the M susy -S degeneracy. Finally, in Sec. VI, we provide concluding remarks and discuss the implications of our results for future axion and supersymmetry searches.

II. AXION COSMOLOGY
In this section, we review axion cosmology without focusing on the possible string nature of the axion itself. See, e.g., Ref. [154] for a recent comprehensive review. As stated previously, we consider for definiteness the case in which the axion decay constant is larger than the Hubble rate during inflation H I , implying that the axion is present during inflation. The initial value of the axion misalignment angle θ i is drawn randomly from the unit circle, with the randomness arising from the Peccei-Quinn mechanism. The axion field is frozen in its initial configuration until the Universe's expansion rate has slowed down to a value comparable to the axion mass. At this point, the axion field starts to oscillate about the minimum of its potential and the axion number density redshifts as expected for a nonrelativistic matter component. This transition occurs when the scale factor takes a value a osc obtained by requiring m a ≈ 3Hða osc Þ, where the value of the Hubble rate as a function of the scale factor a is given in terms of the Hubble parameter and scale factor at the present time, H 0 and a 0 , as follows: Here, we denote by Ω i ≡ ρ i =ρ crit the ratio of the current energy density of species i ∈ fΛ; m; rg (with Λ, m, and r corresponding to dark energy, matter, and radiation, respectively) to the current critical energy density ρ crit ¼ 3H 2 0 M 2 Pl , with M Pl the reduced Planck mass. The present energy density of cold axions, ρ a ¼ ð1=2ÞΛ 4 a θ 2 i ða osc =a 0 Þ 3 , must not exceed the present CDM energy density measured by Planck ρ CDM ∼ 10 −47 GeV 4 [155].
As for any nearly massless scalar field, axions inherit quantum fluctuations from the inflationary period, with a standard deviation σ θ . Primordial quantum fluctuations later develop into isocurvature perturbations [156], which modify the number density of axions, since the gaugeinvariant entropy perturbation S a is nonzero [157][158][159].
Cold axions that spectate inflation differ from thermally produced weakly interacting massive particles because of these imprints from isocurvature fluctuations, the amplitude of which is related to the energy scale of inflation. In this work, we focus on single-field, slow-roll inflation for which the amplitude of the bispectrum is related to the spectral tilt of the power spectrum [160], and a relation between the Hubble rate at the end of inflation H I and the tensor-to-scalar ratio r emerges [161][162][163]. In this framework, the standard deviation of the axion field in units of the decay constant is σ θ ¼ H I =2πf a . In the following, we assume that there are no couplings between the axion and the inflaton field other than gravity. Other scenarios have been discussed in Ref. [164] in relation to axion isocurvature fluctuations. Isocurvature fluctuations can also be suppressed by coupling the axion to a hidden sector [165], which we do not take into consideration here.
We have parametrized the power spectrum of the isocurvature fluctuations at the scale k 0 according to [166][167][168] The Planck Collaboration [169][170][171][172][173] constrains both the primordial isocurvature fraction β and the tensor-to-scalar ratio r, while the curvature power spectrum is measured to be Δ 2 R ðk 0 Þ ≈ 2.2 × 10 −9 [155]. In terms of the axion physics quantities, the ratio of the axion energy density to the CDM energy density today, ω, and the tensor-toscalar ratio, r, read An important caveat is in order at this point. In the following, we will conservatively require ω ≤ 1 so that the current axion energy density does not exceed the current CDM energy density. In reality, this is an approximate requirement for two reasons. The first is that, for m a ≲ 10 −27 eV, the resulting ULA actually has a dark energy-like rather than dark matter-like behavior, and so a more appropriate requirement would be ρ a =ρ Λ < 1. Moreover, an earlier analysis using precision cosmological data (including measurements of the CMB temperature and polarization anisotropy power spectra, galaxy power spectrum, and baryon acoustic oscillations) showed that, in the region where 10 −32 eV ≲ m a ≲ 10 −26 eV, ω ≲ 0.05 is required [151] (see also Refs. [150,174]).
In our work, we conservatively choose to only require ω ≤ 1 instead, for a variety of reasons. First, a posteriori, we do not expect the choice of setting ω ≲ 0.05 for 10 −32 eV ≲ m a ≲ 10 −26 eV to affect our bounds substantially, since the posterior distributions in the region of m a parameter space in question are already quite suppressed for both the string axion and the ALP analyses (see Fig. 1). Second, our focus in this paper is on the string axion case, and it is worth noting that the mathematical relations connecting M susy and S to m a and f a [Eqs. (7) and (8) to be discussed later] are uncertain to a factor of a few (especially the relation between f a and S), so actually worrying about modeling the exact constraints on ω when there are these other uncertainties at play appears to some extent incongruous. We also remind the reader that our goal is not to provide a full-fledged analysis with precision cosmological data but to get a feel for whether cosmology can actually provide information on M susy and S; we believe that, to this end, modeling the exact constraints on ω is not essential. Of course, for future work aiming to obtain more robust constraints on the string axiverse from precision cosmology data, dealing with the aforementioned problem will be of utmost importance, and we plan to return to this in a follow-up work using the modified Boltzmann solver axionCAMB [151].

III. AXIONS IN STRING THEORY
In this section, we consider axions originating from string theory. Starting from the ten-dimensional low-energy Lagrangian of the heterotic string [175] and reducing to four dimensions by compactifying a six-dimensional manifold Z of which the volume is V Z , an effective Lagrangian describing the field ϕ ¼ a=f a is found to be [39,41,176] where, in various string compactification models [41,177], the axion decay constant satisfies f a ≲ xM Pl =S, with x a factor of Oð1Þ 1 and S the (dimensionless) action of the instantons that break the axion shift symmetry, generating the axion potential. As stated previously, we expect S ≈ Oð200Þ. The dependence on the type of string theory (the string scale l s , the asymptotic expansion parameter g s , and the volume of Z) enters only through the string definition of the reduced Planck scale M Pl ¼ ffiffiffiffiffiffiffiffiffiffiffi 4πV Z p =g s l 4 s and S ¼ 2πV Z =g 2 s l 6 s . For instance, nonperturbative world sheet or membrane instantons violate the shift symmetry [41,177], leading to Λ 4 a ¼ μ 4 e −S , where μ is an UV nonperturbative scale which can be as high as the Planck energy scale.
To connect to more physical quantities, the resulting axion mass is then given by m a ¼ Λ 2 a =f a . Supersymmetry, which almost inevitably appears in any realistic string theory, suppresses the UV nonperturbative energy scale by a factor ðM susy =M Pl Þ 2 , where M susy is the SUSY-breaking scale. We then expect μ 2 ≈ M susy M Pl [178]. In the following, we assume f a ¼ xM Pl =S with x ¼ 1. Since the axion decay constant enters Eqs. (4) and (5) only through the combination f a θ i , choosing lower values x < 1 should not alter our results significantly since the factor x would be reabsorbed into a different value of the misalignment angle θ i (which we are not particularly interested in). Admittedly, this does nonetheless introduce an uncertainty of order unity when converting from S to f a , which should be kept in mind. When performing the analysis, we convert from the parameters fμ; Sg (or alternatively fM susy ; Sg) to the parameters fm a ; f a g through In terms of the string parameters, Eqs. (4) and (5) read

IV. ANALYSIS
In this section, we describe the method used to obtain observational constraints on the string axiverse parameters. Our aim is to obtain observational constraints on the four parameters μ, S, θ i , and β, jointly denoted by Θ 1 , in light of observational data d. To do so, we perform a Bayesian analysis, for which we need to specify priors for the parameters Θ 1 .
We begin by discussing our choice of priors. For the initial misalignment angle θ i , we choose a uniform prior over the region ½−π; π, consistent with the parameter being drawn randomly from the unit circle. Our prior on β is given by the posterior distribution for the same parameter obtained by the Planck Collaboration analyzing the Planck TT,TE,EE+lowP dataset (including temperature as well as large-scale and small-scale polarization data) at the wave number k 0 ¼ 0.002 Mpc −1 and provided in Ref. [173]. The fractional primordial contribution of isocurvature modes at the comoving wave number considered is constrained as β < 0.02 at 95% C.L. The distribution of β peaks roughly at β ≈ 0; i.e., it is consistent with an upper limit and not a detection. We have also tested the results =2 in the context of the modelindependent heterotic string. against a different prior based on the TT þ lowP dataset, obtaining that results are not sensitive to the different choice on the prior on β. We choose a flat prior for S ∈ ½50; 450. The range is chosen to match the theoretical expectation S ≈ 200, and we verify a posteriori that it is broad enough to not cut the posterior where the latter is significantly nonzero. Finally, we impose a flat prior for log 10 ðμ=TeVÞ ∈ ½7; 16, where the choice for the range considered corresponds roughly to M susy ≈ μ 2 = M Pl ∈ ½100 GeV; M Pl , i.e., to the SUSY-breaking scale lying between the electroweak scale and the Planck scale. The upper limit in the choice on the range for M susy conforms to theoretical expectations from realistic string theories, whereas the lower limit is consistent with the nonobservation of supersymmetric partners at colliders.
On top of the priors we discussed, we include one further prior on the set of parameters Θ 1 , conditioned on the value of the tensor-to-scalar ratio r computed from Θ 1 ≡ fμ; S; θ i ; βg through Eq. (10). We choose the prior in such a way that it reflects constraints on the tensor-to-scalar ratio obtained by the Planck Collaboration by analyzing the Planck TT; TE; EE þ lowP dataset. Operationally, for each point in parameter spaceΘ 1 selected by our Markov chain Monte Carlo algorithm (to be discussed briefly later), we first computer using Eq. (10); then, to this point in parameter space, we assign a prior probability the value of which is numerically equivalent to the posterior probability distribution for the tensor-to-scalar ratio r determined by the Planck Collaboration and evaluated atr (or, equivalently, we reweigh each point in our Markov chain by this value). 2 This distribution peaks at r ≈ 0 and constrains r < 0.12 at 95% C.L.; i.e., the distribution is consistent with a nondetection, similarly to what discussed for the prior on β.
On top of these priors, we further restrict the available parameter space by requiring the following: (i) The axion decay rate into two photons is smaller that the present expansion rate of the Universe; this imposes further cuts in the m a -f a subspace, so in turn on the μ-S subspace through Eqs. (7) and (8). (ii) The axion is present during inflation, so we demand f a > H I =2π, with H I the Hubble rate at the end of inflation. Since H I can be expressed as a function of r and hence of all four parameters in Θ 1 , this condition imposes further cuts over the whole parameter space. (iii) The current axion energy density does not exceed the CDM energy density, i.e., that ω ≤ 1, with ω given by Eq. (9). Since ω can be expressed as a function of μ, S, and θ i , this condition restricts the available μ-S-θ i subspace, without affecting β. Our discussion so far has been concerned with string axions and on the parameters Θ 1 ≡ fμ; S; θ i ; βg. As discussed in the Introduction, after converting the resulting bounds on μ and S to bounds on m a and f a using Eqs. (7) and (8), we are a posteriori brought to also consider a more generic analysis in which we sample directly on the latter two parameters, with phenomenology-inspired priors on the two. We therefore perform a separate analysis wherein we consider the parameter space spanned by Θ 2 ≡ fX; Y; θ i ; βg, with X ≡ log 10 ðm a =eVÞ and Y ≡ log 10 ðf a =M Pl Þ. This choice is driven by our expectation that m a and f a evenly span various orders of magnitude [42]; hence, it is more appropriate to work with their logarithms (however, see Refs. [179]). We impose uniform priors on X and Y within the ranges ½−40; 8 and ½−10; 0, respectively. We also further impose the same bounds as discussed for the string axion analysis (i.e., the bounds concerning the axion decay rate into two photons, the axion spectating inflation, and the axion energy density not exceeding the CDM energy density today). The range chosen for X is very broad and conservative and allows for values of m a < H 0 ∼ 10 −33 eV, for which the axion has yet to begin oscillating. The range for Y is chosen such that the decay constant f a is sub-Planckian and is bounded from below by the negative results of the CAST searches [143].
To sample the posterior distribution, we use Markov chain Monte Carlo (MCMC) methods. We use the Metropolis-Hastings sampler implemented in the cosmological MCMC package MONTEPYTHON [180], which we configure to act as a generic sampler. From the generated MCMC chains, we compute the joint and marginalized posterior probability distributions of the four parameters and, in particular, of the axion mass and the axion decay constant. From now on, we quote credible regions for the parameters at 68% C.L. unless otherwise stated, whereas all upper limits are quoted at 95% C.L. to conform to standard practices.

V. RESULTS
In this section, we describe the results obtained using the methodology outlined previously. We begin by considering the string axion case, with parameter space described by Θ 1 ≡ fμ; S; θ i ; βg (or equivalently fM susy ; S; θ i ; βg). Our first analysis yields a poor constraint on μ (and hence M susy ), which remains basically unconstrained due to the strong degeneracy between μ and S which data are unable to lift, while we find S ¼ 198 AE 28. One might at first glance be surprised by the fact that S is quite well constrained whereas μ is not. However, an explanation for this puzzling observation is readily found by examining Eqs. (9) and (10), in which we see that S enters exponentially in the observables, and hence it is not possible to vary S too much without then spoiling agreement with observations. In particular, high values of S are excluded by the limits on r, since r ∝ e 2S , whereas low values are excluded by the requirement that The degeneracy between μ and S is easy to understand if we glance at Eqs. (9) and (10). There, we see that μ and S enter the observables through the combination μ 4 e −S , and the data are unable to break this degeneracy (this is somewhat similar to the case of the CMB temperature power spectrum, in which the overall amplitude depends on the combination A s e −2τ , with A s and τ thus strongly correlated and the degeneracy only being partially broken by including polarization data; see, e.g., the discussion in Ref. [181]). We thus expect a strong correlation between the two parameters. This is confirmed by our analysis, wherein we find a correlation coefficient of 0.93 between the two parameters, which are thus close to being perfectly correlated. In Fig. 2, we show a triangular plot in the log 10 μ-S space, which clearly shows this degeneracy.
Using Eqs. (7) and (8), we translate the obtained bounds on μ and S into bounds on m a and f a . We do so by using the aforementioned relations and converting each fμ; Sg sample in our MCMC chains into an fm a ; f a g sample so that the latter two are effectively derived parameters. Notice that doing so retains the information on the priors we used on μ and S. The resulting marginalized posterior probability distribution for m a , normalized to its maximum value, is given by the solid black line in Fig. 1. We find log 10 ðm a =eVÞ ¼ −21.5 þ1.3 −2.3 and f a ¼ ð5.1 þ0.5 −0.9 Þ×10 −3 M Pl . The reason for m a and f a being quite well constrained despite μ being unconstrained is due to the fact that S is quite well constrained. Note that the values obtained for m a and f a are centered around what was originally proposed for ULAs [64], corresponding to an axion decay constant of the order of the GUT scale and with a Compton wavelength of order m −1 a ≈ 0.1 pc, thus with a de Broglie wavelength of galactic size. We therefore conclude that cosmological data complemented with our choice of string-inspired priors favor the region of the string axiverse that comprises an ULA of mass m a ≈ 10 −22 eV. This result is perhaps somewhat artificial, as it depends on the choice of priors and in particular on our choice of having a sub-Planckian M susy . We will discuss in much more detail the impact of this choice on our results in the final paragraphs of this section.
The strong degeneracy between μ (or M susy ) and S suggests that an improved determination of the former parameter(s) would be greatly enhanced by an improved determination of the latter. For purely pedagogical purposes, we explore how the constraints on M susy (currently unconstrained) would change if we were to fix S to a selected value. We choose two representative values of S, namely, S ¼ 198 (corresponding to the inferred central value of S) and S ¼ 153 (corresponding to the 2σ lower limit on S). One could, for instance, imagine considering a specific string realization wherein the value of S is so well determined from theoretical considerations that it is for all intents and purposes fixed. In the former case, we find log 10 ðM susy =TeVÞ ¼ 7.8 þ1.4 −2.3 , whereas in the latter case, we find M susy ≲ 8 × 10 3 TeV at 95% C.L., indicating that within this scenario SUSY could possibly be within reach of the proposed 100 TeV collider [132]. We show the marginalized posterior distribution, normalized to its maximum value, of log 10 ðM susy =TeVÞ when S is fixed to 198 and 153 in Fig. 3, with the dashed vertical line corresponding to M susy ¼ 100 TeV.
We now perform our fourth analysis, in which we consider ALPs (thus disregarding their fundamental origin) with the parameter space described by Θ 2 ≡ fm a ; f a ; θ i ; βg and flat priors on log 10 m a and log 10 f a . The marginalized posterior probability distribution for m a , normalized to its maximum value, is given by the dashed black line in Fig. 1. In particular, we find X ¼ log 10 ðm a =eVÞ ¼ −11.1 þ11.5 −5.9 , mildly favoring higher values for m a and mildly disfavoring ULAs, although not at a statistically significant level. Two notable features are discernible from the posterior distribution. The first is the complete loss of sensitivity (flat posterior) at very low values of m a ≲ H 0 , due to the fact that axions which are extremely light have yet to begin oscillating. In the very high mass range [m a ≳ OðeVÞ], instead, the posterior distribution is sharply cut by the requirement that the axion is present during inflation, hence f a ≳ H I . The requirements that ω ≤ 1 and that r ≪ 1 also cut the posterior distribution at low and high masses, respectively, as evident from Eqs. (4) and (5).
For the axion decay constant, we find Y ¼ log 10 ðf a =M Pl Þ ¼ −5.5 þ1.8 −3.0 . The central value corresponds to f a ≈ 10 13 GeV, which is slightly higher than what is expected for the QCD axion to be the CDM particle. Moreover, we find that X and Y are strongly anticorrelated, with a correlation coefficient of −0.92. We show this Fig. 4, in which we plot the joint two-dimensional posterior distribution in the X-Y plane along with the 68% C.L. (dark blue) and the 95% C.L. (light blue) contours. The degeneracy is approximately along the axis defined by the relation m a f 2 a ≈ const. This can be understood by combining Eqs. (4) and (5) for a radiation-dominated universe (m a ≳ 10 −27 eV), for which we obtain r ∝ ðm a f 2 a θ 2 i Þ −1 β. Overlaid on the same figure is the m a -f a relation for the QCD axion (black dashed line), which instead lies along the axis defined by m a f a ¼ Λ 2 QCD , as well as the central value of the analysis by Klaer and Moore (which for the first time included large short-distance contributions to the axionic string tension in their numerical simulations), which yields m a ¼ ð26.2 AE 3.4Þ μeV [182].
If we interpret the results obtained in Fig. 4 in terms of the string axiverse parameters μ and S, the preferred value for the decay constant forces S ∼ M Pl =f a ≈ 10 6 . This is far from the theoretically preferred value S ≈ Oð200Þ. Moreover, it leads to a trans-Planckian value of M susy ¼ m a e S=2 =S for any reasonable value of the axion mass, due to the exponential sensitivity of the SUSY scale to S. In fact, demanding that M susy be at most equal to M Pl , we find that the size of the cycle is at most S ¼ 145 (S ¼ 240) for m a ¼ 10 −2 eV (m a ¼ 10 −22 eV), a value which is consistent with Oð200Þ. 3 In fact, the consistency of the relation between m a , M susy , and S in Eq. (7) with the theoretical prior M susy < M Pl leads to S ∼ 200 to within an uncertainty of about ∼50. The role of the data is that of further shrinking the error bar, leading to the result S ¼ 198 AE 28 we have reported. The role of the data in further improving the determination of S relies on the fact that the upper limits on ω and r squeeze the allowed region of parameter space from opposite directions; see Eqs. (9) and (10), with the limit on ω excluding the tail at low values of S, whereas the bounds on r exclude higher values of S.

VI. DISCUSSIONS AND CONCLUSIONS
In this work, we have for the first time attempted to constrain fundamental parameters describing the string axiverse, using cosmology. We have in particular focused on the SUSY-breaking scale M susy and the dimensionless  4. Two-dimensional joint posterior distribution in log 10 ðm a =eVÞ-log 10 ðf a =M Pl Þ parameter space, obtained in our axion analysis with flat priors on log 10 m a and log 10 f a , with 68% C.L. and 95% C.L. credible regions corresponding to the dark and light blue regions, respectively. Also shown is the relation defining the QCD axion (dashed black line) and the value of the axion mass predicted by Klaer and Moore [182] assuming that the axion makes up the totality of the dark matter (red diamond), m a ¼ ð26.2 AE 3.4Þ μeV.
action of the string instantons responsible for breaking the axion shift symmetry, S. Imposing string-inspired uniform priors on log 10 M susy and S and using current observational bounds on the tensor-to-scalar ratio, the primordial isocurvature fraction, and the energy density of dark matter, we have performed a Bayesian inference analysis to constrain the parameters characterizing the string axiverse.
We have found that M susy is essentially unconstrained (due to a strong parameter degeneracy with S), while we find S ¼ 198 AE 28 at 68% C.L. which is partly due to consistency with the theoretical priors as discussed in the previous section. When interpreting these results in terms of the more familiar axion mass and decay constant, m a and f a , through Eqs. (7) and (8), we find that the lower range of m a is somewhat artificially favored, with log 10 ðm a =eVÞ ¼ −21.5 þ1.3 −2.3 and f a ¼ ð5.1 þ0.5 −0.9 Þ × 10 −3 M Pl at 68% C.L.. These values lie within the ultralight axion range; see the solid black line in Fig. 1.
To break the M susy -S degeneracy, for purely pedagogical purposes, we explore the impact of fixing S ¼ 198 or S ¼ 153, respectively, corresponding to the central value and 2σ lower limit from the first analysis. This respectively gives log 10 ðM susy =TeVÞ ¼7.8 þ1.4 −2.3 or M susy ≲ 8 × 10 3 TeV. The latter bound is particularly interesting, since it suggests that such a scenario could hypothetically be probed in a future 100 TeV collider [132].
We then explored the impact of our choice of stringinspired priors on the inferred values of m a and f a . We have done so by performing a different analysis in which we directly sample the m a -f a parameter space, with flat priors in log 10 m a and log 10 f a . We have found log 10 ðm a =eVÞ ¼ −11.1 þ11.5 −5.9 and log 10 ðf a =M Pl Þ ¼ −5.5 þ1.8 −3.0 at 68% C.L., which mildly disfavors the lighter end of the axion mass spectrum and is consistent with the value of the mass inferred by Klaer and Moore assuming that the axion makes up all the dark matter [182]; see the dashed black line in Fig. 1. Clearly, the difference between the two analyses is at least partly due to the different choices of priors, which, however, we argued arise quite naturally when considering the two different theoretical or phenomenological approaches. In the string-inspired analysis, our results are consistent with theoretical works suggesting that axions originating from string theory preferentially populate the low-mass end of the string axiverse [42], although we find that only a relatively narrow band around the value m a ≈ 10 −22 eV is consistent with a combination of theoretical priors and data; see the solid black line in Fig. 1. Ultralight axions with mass of order 10 −22 eV have been ruled out as the dark matter component by various methods: matching the profile of the Fornax dwarf [183], Jeans analysis on numerically simulated dark matter halos [184][185][186], and the Lyman-α forest [148,149].
Recently, the work in Ref. [153] appeared, and a global fit that uses a Bayesian analysis technique to explore the parameter space of the QCD axion based on the code GAMBIT [187] and its module DARKBIT [188] is presented. In particular, Ref. [153] considers the scenario in which the Peccei-Quinn symmetry breaks during a period of inflation, although, at this stage, a fit that accounts for the inflation module has not been implemented yet, so the analysis reported does not take into account bounds from isocurvature fluctuations. A further difference with our analysis lies in the relation m a f a ¼ Λ 2 QCD ¼ ð75.5 MeVÞ 2 , which holds for the QCD axion. The likelihood analysis with GAMBIT takes into account various results, including laboratory experiments from light-shining-through-wall, helioscopes, and cavity searches as well as astrophysical observations for the distortions of gamma-ray spectra, supernovae, horizontal branch stars, and the hint from the cooling of white dwarves. Given these differences, the marginalized posterior distribution obtained in Ref. [153] when demanding that the totality of dark matter is in axions gives the range 0.12 μeV ≤ m a ≤ 0.15 meV at the 95% equal-tailed confidence interval. Their results are also dependent on the choice for the prior on f a as well as those on the axion-toelectron coupling and the anomaly ratio.
In Refs. [151,152], a Bayesian technique using a MCMC sampling has been performed using the cosmological parameter space fA s ; n s ; τ; Ω b h 2 ; h; Ω c h 2 ; Ω a h 2 ; m a ; H I g, where the amplitude of scalar fluctuations A s and the scalar spectral index n s are defined at the pivotal scale k 0 ¼ 0.05 Mpc −1 , τ is the optical depth to reionization, Ω b h 2 is the baryon density, and Ω c h 2 is the density in dark matter other than axions, so that Ω CDM ¼ Ω a þ Ω c . The analysis in the paper is performed using the axionCAMB code [151] and reveals no evidence for an axion component in the mass range 10 −33 eV ≤ m a ≤ 10 −24 eV, extending similar previous findings [189].
In conclusion, in this work, we provide for the first time a new window into the string axiverse and the SUSYbreaking scale from cosmology. There are plenty of avenues for follow-up work. First, it would be worth performing a full-fledged analysis carefully taking into account precision cosmology data, such as data from the CMB temperature and polarization anisotropy spectra, as well as from galaxy surveys (galaxy power spectrum and/or baryon acoustic oscillations) and weak lensing surveys. Such an analysis requires a suitable modification of the axionCAMB code [151] and would be especially relevant given that our analysis of the string-inspired axion constrains the axion mass to reside within a narrow region where the axion affects the matter power spectrum at small, but still cosmologically relevant, scales. In addition, it would be useful to revisit the theoretical uncertainties entering Eqs. (7) and (8) relating M susy and S to m a and f a , in order to strengthen our analysis. We plan to return to these and other issues in future work.