Cosmic birefringence from the Axiverse

We revisit the evidence for CMB birefringence in the context of a rich Axiverse. Using probability density functions (PDFs) for various axion parameters, such as the mass and axion decay constant, we construct the PDF for the cosmic birefringence angle and investigate its properties. By relating the observed value of the birefringence angle to the mean or standard deviation of the constructed PDF, we constrain the shape of the input PDFs, providing insights into the statistical distribution of the Axiverse. We focus on three different types of axion potentials: cosine, quadratic, and asymptotically linear axion monodromy. Our analysis showcases the potential of cosmic birefringence in constraining the distribution of axion parameters and uncovering possible correlations among them. We additionally offer predictions for “birefringence tomography”, anticipating future measurements of birefringence from lower multipoles, and show how it can be used to rule out simpler versions of the Axiverse. Our findings contribute to the ongoing exploration of the Axiverse and its implications for cosmic birefringence.

Several analyses of the CMB data have found an intriguing signal of Cosmic Birefringence [1][2][3][4], which measures the angle of rotation of the CMB radiation, with the most recent value being β = 0.342 +0.094  −0.091 , excluding the null result by more than 3σ.This angle refers to the static and isotropic rotation of the polarization of photons emitted at the "last scattering" (LS) surface which shows up in CMB experiments, specifically in the EB cross-correlation signal C EB l ∝ sin(4β)(C EE,CMB l − C BB,CMB l ) (see Ref. [5] for a recent review).This measurement has been notoriously complicated because the birefringence angle β is degenerate with the miscalibration of the polarimeters, but the new technique developed in Refs.[6,7] made it possible to distinguish the two and achieve better precision.The measurement appears robust under instrumental systematics, but potential EB contribution from galactic dust currently represents the biggest limitation of the analysis [8].
The birefringence signal finds an elegant explanation in the context of axion physics where a light scalar field ϕ, which can be also tied to dark matter or dark energy, couples to the electromagnetic field strength F µν via the Chern-Simons interaction ϕF F .This additional term in the Lagrangian breaks the parity of the electromagnetic sector and changes the dispersion relation between the left and right-handed photons.This effectively makes the Universe, filled with such scalar field ϕ, behave as a birefringent material with a different refraction index for the right and left-handed polarization.Thus the linearly polarized light travelling in such a "medium" experiences a rotation of the direction of polarization, which is called Cosmic Birefringece.A new way to understand such a phenomenon is as a "photon chiral memory effect" [9].This memory effect measures the permanent change of the spin angular momentum of electromagnetic fields by chiral symmetry violating processes (ϕF F in our case) occurring in the bulk.This reveals an interesting connection between the birefringence effect, electromagnetic memory effect, and the asymptotic symmetries/charges at null infinity1 [9].
Within the axion framework, the cosmic birefringence angle depends only on the field displacement between the times of emission and observation and it is frequency independent [10][11][12][13], which is consistent with current results [3].The angle of rotation of CMB light is thus significant when the scalar field has a non-trivial evolution between LS and today.In particular, the dominant contribution comes from fields that start oscillating in this time window, corresponding to axion masses of H 0 ∼ 10 −33 eV ≤ m ≤ H LS ∼ 10 −29 eV.We focus on axion fields produced via the misalignment mechanism with decay constant f a greater than the energy scale of inflation H I , such that each field is mostly uniform in our observable patch of the Universe [14][15][16][17].This is the most common scenario discussed in the literature [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33] to explain the static and isotropic nature of the signal, which is generated through the evolution of the axion background.In this case, the unavoidable fluctuations that would be generated during inflation are suppressed by a factor of ∼ H I /f a with respect to the isotropic component, in line with the current absence of evidence of an anisotropic counterpart [34][35][36].Note that in the whole mass range that could give rise to the signal (10 −40 eV ≲ m ≲ 10 −26 eV [27]), the axion field can be responsible for dark energy or comprise a subdominant fraction of dark matter.With future experiments, it will be possible to access higher masses by examining the different effects on β of the CMB polarization at early and late times [26].The signal could be also generated by axions forming a long-lived network of domain walls or cosmic strings; this scenario would be distinguishable through a characteristic imprint on the anisotropic spectrum of cosmic birefringence [37][38][39][40][41].
However, one can argue for the existence of many axion fields or axion-like particles (ALPs), that were dynamical in various stages of cosmic evolution.In the context of the String Axiverse [17,42,43][ [44][45][46], many axion fields exist (hundreds or thousands) and are spread over several orders of magnitude in mass 2 .We thus explore the implications of having β ∼ 0.3 deg in the presence of many dynamical axions or ALPs.Previous work has explored the case of a single field [27,28] where constraints on the axion-photon coupling g ϕγ are inferred by inverting the relation β = − 1 2 g ϕγ [ϕ(t LS ) − ϕ(t 0 )], where t 0 and t LS correspond to the present and "last scattering" time and the initial condition is fixed by the maximum axion abundance in the relevant mass range [47,48].In this work, we follow a different approach: we consider the field's initial value and the axion decay constant as random variables drawn from theoretically motivated distributions.We can then calculate the statistical distribution of the resulting birefringence angle as a function of the number of axions and check how many are needed to match the observed signal within one standard deviation, for different distributions of the input parameters (axion mass, decay constant etc).
Interestingly, we find that combining the results for the cosmic birefringence and the axion abundance gives a generic upper bound for the central value of the decay constant, independently of the other parameters for the quadratic potential.Regarding the mass distribution, so far the cosmic birefringence cannot constrain it directly, but we show that the projected abundance at higher masses gives important information on the maximum allowed mass.Moreover, the future 2 Interestingly, cosmic birefringence provides a complementary view of the Axiverse, compared to other axion-driven phenomena, such as supperadiance, due to the different axion masses involved.However, CMB birefringence requires the axions to couple to the SM photon, which the supperadiance mechanism does not.measurement of the birefringence from lower multiples, corresponding to photons emitted around reionization [49][50][51], can be used as a probe of the mass distribution which is typically uniformly distributed in logarithmic space.We also discuss the impact of correlations, especially between (ϕ in , f a ) and (m a , f a ), that can emerge from the underlying UV theory [43,52].In the second part of the paper, we repeat the study for axions with monodromy potential that have a different monomial dependence for field values larger than some new transition scale and we show how the constraining region changes in these cases.
The paper is structured as follows: In Section II we introduce the model, describe the single field result for axion-driven CMB birefringence and introduce the notion of a statistical treatment in the presence of many dynamical axion fields in the universe.In Section III we present the manyfield results in the case of the usual cosine potential.In Section IV we analyze the case of many axions with independent quadratic potentials.This allows for an in-depth analysis and for several results to be derived analytically.In Section V we use axion monodromy potentials, focusing on potentials that are asymptotically linear at large field values.We discuss our findings and propose future directions in Section VI.

II. MODEL
The Lagrangian of an axion-like particle (ALP) ϕ coupled to a U (1) gauge field, in this case the Standard Model photon, can be written as where c is the anomaly coefficient and α EM ≃ 1/137 is the fine-structure constant.This leads to the well-known equation of motion for the two circular polarizations of the gauge field in an expanding universe with a scale factor a: where dots denote derivatives with respect to cosmic time and H = ȧ/a is the Hubble parameter.It is evident from the dispersion relation ω 2 ± ≡ k 2 a 2 ∓ k a g ϕγ φ that the two circular polarizations (±) are affected differently by the rolling pseudo-scalar field ϕ.This can have interesting phenomenology, including efficient transfer of energy to gauge fields, the fragmentation of ϕ and the generation of gravitational waves [53][54][55][56].In the current work, we are instead interested in the different phases of the two polarizations of CMB photons, as they propagate through the Universe under the effect of a slowly rolling ALP [10,12,18].The total birefringent angle β affecting the CMB photons generated by the intervening axion field is given by the integral over the line of sight of the difference ω − − ω + which is shown to be [12] g ϕγ φ/2 and can be separated into two parts [27]: The first contribution comes from the evolution of the background field whereas the second contribution is sourced by the primordial, quantum in origin, perturbations.We consider axion fields with misalignment production and the U (1) symmetry broken before or during inflation, such that during LS the axion fields are homogeneous over the whole sky, except for perturbations of the order δϕ ∼ H I /f a ≪ ϕ LS [57].Thus the evolution of the background field and its local perturbations ∆ϕ and δϕ 0 produce a uniform rotation of the polarization plane over the whole sky called isotropic cosmic birefringence.On the other hand, the fluctuations at LS lead to a spatially dependent rotation called anisotropic cosmic birefringence, because the corresponding values of fluctuations in various directions of the sky are uncorrelated at distant positions of the LS surface [58][59][60][61].In our analysis, we focus on axions with masses H 0 ≤ m ≤ H LS because, for similar initial conditions, the field displacement of these fields dominates the signal by orders of magnitude.
Indeed, the contribution from lighter fields m ≪ 10 −33 eV is suppressed by (m/H 0 ) 2 ≪ 1, this is because they are just starting to roll and their evolution deviates from the initial value by a factor (mt 0 ) 2 ≪ 1 where t 0 ∼ H −1 0 is the age of the Universe.On the other side, heavier fields m ≫ 10 −29 eV experience a wash-out effect due to the several oscillations over the finite duration of the LS [26][42].In the next section, we briefly review the evolution of a rolling axion field during the matter-dominated era (the complete analysis can be found in Ref. [17]) that corresponds to axions with masses in the window of interest, allowing us to use the corresponding equations in our analysis.Subsequently, we discuss some theoretical motivations for choosing the form of the probability density functions (PDFs) for the axion properties (including the mass and axion decay constant), drawing inspiration from the String Axiverse scenario.

A. ALP background motion as a source of CMB birefringence
The axion field evolves under the Klein-Gordon equation in an expanding universe, which in the quadratic approximation of the potential is written as If we consider the ansatz a(t) ∝ t p , where p = 1/2 in radiation-domination and p = 2/3 in matter domination (we assume that the scalar field does not dominate the expansion rate), the scalar field evolution is expressed in terms of Bessel functions as: with n = (1/2) 9p 2 − 6p + 1.We neglect the solution Y n since it is divergent at early times, whereas we know that the initial value is determined by the misalignment angle.The parameter A is fixed by the initial field amplitude.The main features of the evolution are captured by J n , which is approximately constant when mt ≪ 1, meaning that the field is "frozen", and it oscillates as cos(mt) when mt ≫ 1, corresponding to dark matter behaviour.The transition between the two behaviours is well approximated by mt = 1 or, using H = p/t, H osc = pm.After the transition, one can easily show that for a osc ≳ a LS , i.e. in matter domination, the field amplitude scales like ϕ(t) ≃ ϕ i [a osc /a(t)] 3/2 and ϕ osc ∼ ϕ i .Thus the axion displacement, in the mass range of interest, is ∆ϕ ≃ −ϕ LS because ϕ 0 ≪ ϕ LS and the latter is related to the present field value as: Because the field is frozen before LS, we can approximate ϕ LS ≃ ϕ in and relate it to the current abundance as: Therefore in this mass regime, the field displacement takes the value Note that in this mass regime, the current density is more constrained compared to higher or smaller axion masses from linear cosmology as studied in Ref. [47] and more recently in Ref. [48].
In the single-field axion case, the birefringence angle is given by: Combining Eqs. ( 9) and (10) leads to Using the measurement β obs ∼ 0.3 deg, one finds f a /c ≤ 10 −2 M Pl , which fits well with the theoretical predictions from string theory [42,[44][45][46].We later show how these results extend in the case of a large number of fields.

B. Statistics of many ALP's
In this work, we consider the existence of many axions whose parameters are distributed following some probability density function (PDF).These PDFs can be either derived through some underlying UV theory or taken as phenomenological "parameters", which can be constrained through experiments and observations.Indeed, in String Theory, the axions are naturally linked to the geometry of the compactification that determines the distribution of the mass m and decay constant f a [17,42,45,46,62,63].Generically, instantons generate a potential for the axions that is exponentially suppressed by the instanton action, thus m ∝ e −S ins .This last quantity is linked to the parameter of compactification which is likely to be uniformly distributed rather than being concentrated around one particular scale [42], giving a uniform mass distribution in logarithmic space (see App E of Ref. [17]).This distribution has been found in different compactification scenarios, for example [46,62,63], and numerically in Ref. [43], at least far from the Planck scale.The extremes of the mass distribution are model-dependent, thus we take an agnostic point of view and consider the mass bounded from above by the Planck mass m max ≲ M Pl .Determining minimum mass is more difficult, since it is also connected to dark energy [64].We follow Refs [42,43,62] and choose m min ∼ H 0 .
String theory, in general, predicts a non-flat distribution for the decay constant which has the following parametric dependence f a ∼ M Pl /S ins for a single field [45].For the Axiverse, it has been shown that in several cases it follows a Gaussian distribution in logarithmic space with a decreasing mean value for an increasing number of axions [43,52,65,66].In particular, in Ref. [66] it is explained how the log-normal distribution of the decay constant, computed as f 2 a = m 2 /λ, as the product distribution of the mass and the axion self-interaction coupling λ that follows a similar distribution along with a strong correlation between the two which make the distribution peak at a certain scale.For comparison, in section IV C we discuss the results for the log-normal and log-uniform distribution for the decay constant.
The initial value of an axion in a cosine potential can be taken to be uniformly distributed in the fundamental domain of the cosine.In a non-periodic potential, like a quadratic or axion monodromy potential, de-Sitter fluctuations during inflation can generate a large field value.For simplicity, and in order not to include specific inflationary dynamics in our analysis, we assume that the distribution of initial field values follows a Gaussian distribution with a spread σ ϕ < M Pl .
We explore these distributions in the context of CMB birefringence by varying their mean and standard deviation.Furthermore, we also consider the existence of correlations between the different parameters , as correlations have been shown to arise in the Axiverse [43,52].For example, a correlation between the initial field value and the decay constant can be understood as the initial field value given by the product of the initial misalignment angle and the decay constant ϕ = f a θ.This last quantity differs from the scale entering in the coupling g ϕγ ∝ c/f a by the anomaly coefficient c, thus we can imagine that in the presence of axion mixing the relation ϕ and f a is not exact, but they can be strongly correlated.Ref. [30] provides a concrete example of axion's mixing and the implication for cosmic birefringence.Similarly, a correlation between the decay constant and the axion mass, which is exact for the QCD axion m ∝ f −1 a , can arise for string axions when axion's mixing is present [30,43,52].

III. COSINE POTENTIAL
We start with the cosine potential, as this is the usual potential for axions, arising from instanton effects (see e.g.Ref. [67]).Furthermore, it allows us to make contact with the estimate given in Ref. [42], where β ∝ √ N .If there is no mixing between the different axion fields, such that the decay constants appearing in the potential ϕ i /f a,i and in the couplings g ϕγ,i ∼ f −1 a,i are the same for each axion, the two cancel in Eq. ( 10), since we expect ∆ϕ i ∝ f a,i .The birefringence angle thus depends only on the initial field position within the fundamental domain of the cosine.Considering ϕ i /f a,i uniformly distributed in the interval [−π, π], the average value of β is zero and its variance is controlled by the variance of the initial field value (ϕ in,i /f a,i ) 2 = π 2 /3.Neglecting the effects of the anomaly coefficient Requiring that the standard deviation of the birefringence angle matches the observed value β ≃ 0.3 deg, leads to N ≃ 25 fields, which is a large enough number to justify a statistical treatment.

A. Anomaly coefficient and correlations
We go beyond this simple result by introducing the anomaly coefficient c i that can be different for each axion, drawn itself from a PDF, while we allow it to be correlated to the initial field amplitude θ i ≡ ϕ in,i /f a,i .The correlation coefficient ρ is defined as where σ c , σ θ are the standard deviations of the two distributions.The average value of β in the general case is We typically consider distributions that are symmetric with respect to zero, meaning that they have zero mean.For example, in order to "bias" the distribution of initial field amplitudes θ i , some parity-violating process must occur in the early universe.This is not impossible, but requires some concrete model-building.If at least one of the average values of the field amplitude and anomaly coefficient vanishes, then the average value of β scales linearly both with the number of fields, as well as with the correlation coefficient, becoming zero for uncorrelated c i and θ i .
The calculation of the variance of a product of random variables is more complicated; we present all relevant formulas in Appendix A. In the case of Gaussian random variables with mean µ c and µ θ (which we can use as a qualitative guide for our case), the standard deviation becomes When both random variables have zero average values, the above equation becomes We see that this is similar to Eq. ( 13) for σ θ , σ c , ρ ∼ 1.Most importantly however, the different scaling of ⟨β⟩ with N and σ β with √ N , means that ⟨β⟩ ≫ σ β in the many-field limit, unless |ρ| ≪ 1.If ⟨c i ⟩⟨θ i ⟩ = 0 and ρ ̸ = 0, the mean of the birefringence angle will dominate over its standard deviation.This corresponds to the case of aligned initial conditions discussed in Ref. [43], because there is an overall tendency for the axions to give either positive or negative contribution to β.We do not pursue this case further here and instead focus on independent distributions for the various parameters that determine β.We will return to the subject of correlated random variables in Section IV C, for the case of a quadratic axion potential.

B. Testing the mass distribution with birefringence tomography
Subsequently, we can ask about the total number of axions needed to explain the observed signal β obs ∼ 0.3 deg for a given PDF of the axion mass.We consider a log-flat distribution As discussed in section II B, we take leaving some freedom for the lower mass cutoff.With this (simple) prior for the mass distribution, the probability that an axion mass falls within the interesting parameter space for birefringence is given by P(−33 ≤ log(m/eV) ≤ −29) = 4/(27 + 33 log 10 (m min /H 0 )).Taking m min ∼ H 0 and N ≃ 25 axions in the interesting mass range, as follows from Eq. ( 13), we arrive at to N dec ≃ 6 axions per decade of mass and N tot ∼ 360 total axions with any mass.This number is compatible with the typical expectation of a few hundred axions (see e.g.Ref. [65]).However, postulating a large number of axion fields can lead to inconsistencies with the dark matter abundance.Using Eq. ( 8), which is strictly only true for the quadratic case (and well approximates a cosine, when the field does not start very close to the maximum of the potential), and using ϕ in = O(f a ), we estimate that the contribution of each axion to the dark matter abundance will be Ω ϕ = O(f 2 a /M 2 Pl ).This is very small for axions with a small decay constant and can be catastrophically large for f a = O(M Pl ).
We will return to this point in Section IV B, where we compute the ALP contribution to the DM abundance for the case of a quadratic axion potential.Note that one can increase or decrease the inferred number of axions by changing the mean value of the anomaly coefficient, which can be computed in a definite axion model.For our purposes, we treat it as a phenomenological parameter, whose properties (value, correlations, PDF) can be constrained from the birefringence measurement.
In principle, we can constrain the mass distribution, since it gives a clear prediction for the ratio of the birefringence angle associated with recombination β rec and reionization β rei , that can be extracted by studying the angular dependence of the polarization data [50,51].The angle β rec measures how many axions have rolled to the minimum of their respective potential between recombination and the present time whereas β rei is only sensitive to those axions that started rolling after reionization (H rei ∼ 10 −31 eV).Therefore the ratio between the two gives information about the relative amplitude of the mass parameter space explored in the two time intervals: whereas in the "aligned" case β rei ≃ 0.5β rec .Note that although the result depends on there being a large enough number of axions in the relevant mass window to justify the statistical treatment, it does not depend on the lower-mass cutoff of the mass prior and therefore is formally independent of the total number of axions across all masses.
We can probe the mass distribution of axions using Eq. ( 19) by assuming that the mass distribution over the interesting range does not follow Eq. ( 18), but instead has a preference for larger or smaller masses.A phenomenological distribution with this property is where n m is a dimensionless number.For n m = 0 we recover the log-uniform distribution.Note that in Eq. ( 20) we neglect the overall normalization for simplicity and because it drops out of Eq. ( 19).
Using this "tilted log-uniform" distribution, we can compute β rei /β rec ≃ 2 −(nm+1)/2 .For example a ±20% difference from the result of the log-uniform distribution requires a tilt of n m ≃ −0.5 or n m ≃ 0.7.A larger difference from log-uniform will arise if the mass distribution has a stronger scale dependence, e.g. if it follows a log-normal distribution centered around some preferred mass range.Of course, this result can be degenerate with changes in the potential or introduction of correlations between the various parameters of the model, such as the mass, decay constant and anomaly coefficient (see e.g.Section.IV D).However Eq. ( 19) provides a clean expectation for the log-flat prior of the axion mass.A future detection that is inconsistent with this result would rule out this simple Axiverse model.

IV. QUADRATIC POTENTIAL
We now move to the case of axions with a quadratic potential.This simplification, allows us to derive several results analytically and gain physical intuition by testing a variety of PDFs.
A. Fixed decay constant f a In the case of a quadratic potential, the total birefringence is given by the sum of the individual contributions of axions in the mass range −33 ≤ log 10 (m/eV) ≤ −29, as given by Eq. ( 10): where we reabsorbed the anomaly coefficient into f a .By fixing the axion decay constant and taking the field amplitude as a random variable with a PDF that is symmetric around zero, the standard deviation of the birefringence angle is Therefore, by using the observed value β obs ≃ 0.3 deg and equating that to the standard deviation σ β ≃ β obs , we find a simple dependence of the number of axions on the standard deviation of the initial field value σ ϕ and the fixed decay constant f a : From this expression, we see that the number of axions depends quadratically on the decay constant and the standard deviation of the initial field value, thus by varying these parameters the required number of axions changes significantly.Finally, in order for our statistical analysis to be viable, we must require that N ≫ 1 which gives f a ≫ σ ϕ /10.
It is worth remembering that the (standard deviation of the) birefringence angle scales as √ N , where N is the number of axions3 (with N ≫ 1), whereas for the axion dark matter abundance the corresponding scaling is Ω ϕ ∝ N .This leads to a constraint on the number of axions predicted from the measurement of β, so as not to exceed the allowed DM abundance.Even in the simplest case with fixed f a , we derive non-trivial constraints.Extending Eq. ( 8) for the case of N axions, the total axion abundance is: The sum of the squares of N independent Gaussian variables with zero mean is a χ 2 N −distribution with N degrees of freedom, where both the mean ⟨χ 2 N ⟩ = N and the variance (which is equal to 2N ) scale with the number of degrees of freedom. 4The mean value of the total abundance becomes: The combination σ 2 ϕ N can be substituted using Eq. ( 23), leading to: Requiring Ω ϕ ≤ 0.003 (1% of the DM abundance in this mass range [48]), we find a constraint for the decay constant f a ≤ 2.2 × 10 16 GeV.It is worth comparing this to the constraint on the coupling g ϕγ ≳ 10 −20 GeV derived in Ref. [28] using a single rolling axion in a quadratic potential, which gives f a ≲ 10 17 GeV.We see that the many-axion analysis leads to a somewhat stricter constraint on the axion decay constant.
Let us now briefly discuss the "aligned case", where the distribution of initial field values has a non-zero mean.Assuming that the mean dominates over the standard deviation, the birefringence becomes leading to N = 10f a /⟨ϕ in ⟩, by requiring that β obs ∼ 0.3 deg.In order for the standard deviation of the birefringence angle to be negligible compared to its mean, we require Similarly, the average value of the axion abundance is By using the relation for the birefringence, we get We can derive a loose bound, by only considering the second term, leading to The last relation is a condition on the geometric mean of the axion decay constant and the alignment, as given by the average field value.Another interesting limit is one where all dynamics is controlled by a single scale, the axion decay constant.By taking ⟨ϕ in ⟩ ∼ f a ∼ σ β , the mean birefringence is larger than its standard deviation, by a factor of 3, so our assumption for neglecting the standard deviation is borderline valid.The constraint in this case, arising from avoiding overproduction of axion dark mater, is f ≲ 0.01M Pl .
The main result of this section can be summarized as follows: When ϕ in,i has no preferred sign the DM abundance gives a maximum effective displacement of the axion field ∆ϕ 2 ≃ i ϕ 2 in,i that, combined with Eq. ( 22), leads to an upper bound on the decay constant given by Eq. ( 26) 5 .This relation will change if we change the relation between the field value and the abundance, as we see in Section V for the family of axion monodromy potentials.10 -28 10 -26 10 -24 10 -22 10 -20 10 -18 FIG. 1: The projected abundance for each decade in mass given by Eq. (31).

B. Projected abundance to higher masses
In the axiverse picture, it's interesting to ask whether the expected number of axions we infer from the birefringence measurement, given in Eq. ( 22), is consistent with the abundance constraints at higher masses.For instance, we can test the assumption that the same dynamics, which sets the distribution of the initial field displacement and decay constant, also determines the axion abundance at higher masses.For axions with masses m > H eq ∼ 10 −28 eV, the present abundance for a single axion field is given by [47] where Ω r is the present radiation density.Eq. ( 30) leads to the following abundance for each decade in mass which, using N dec = 25(f a /σ ϕ ) 2 , depends solely on the axion decay constant.We can check that the average abundance, given in Eq. ( 31), does not exceed the current cosmological upper bound on the axion abundance reported for example in Refs.[48] and [68].In Figure . 1 we show the projected abundance and the corresponding upper limit computed as the linear interpolation, in decade of mass, between the following values: The first (stronger) constraint, for axions with m ≤ 10 −28 eV, comes from CMB and Large Scale Structure analyses, which relaxes for m ≳ 10 −25 eV [48].At higher masses, the constraint from the Lyman-α forest becomes stronger, which bounds the axion to be at most 30% of the total dark matter abundance for m < 10 −21 eV [68].
It is interesting to note that the assumption of a unique decay constant for each decade in mass is in conflict with the expectations shown in Figure 1, where the corresponding abundance exceeds the cosmological bounds.For instance, as shown in Figure 1, this occurs at m ∼ 10 −24 eV for f a ∼ 10 16 GeV, whereas the abundance is consistent up to m ∼ 10 −20 eV for f a ≲ 10 15 GeV.This would indicate, for example, that axions in different mass ranges have different decay constants or that they have different production mechanisms.These upper bounds assume that we can do statistics within each decade of mass, which is a rather fine-tuned version of the axiverse.A better estimator is the sum of all axions, whose masses are drawn from a distribution We take the average value of this, leading to where we took the mass and initial field amplitude to be uncorrelated.Taking ⟨ϕ in,i ⟩ = 0 implies that ⟨ϕ 2 in,i ⟩ = σ 2 ϕ .Considering again a log-uniform distribution for the masses, it is straightforward to compute where we took m max ≫ m min .If we assume that the mass distribution is log-uniform up to some large cut-off, we can compute the upper value of this cutoff m max from the constraint on Ω ϕ .The total number of axions is also related to this cutoff, since N tot = N dec log 10 m max /m min .
Interestingly, this makes the final expression for the dark matter abundance depend on m max only through the square root, as Therefore, the final axion abundance is determined by the product of (f a /M Pl ) 2 and m max /H 0 .
Note that Eq. ( 36) is almost identical to Eq. ( 31) with the substitution of m max , indicating that the total abundance is dominated by the more massive axions, as shown in Figure 1.Indeed, because of the Ω ϕ ∼ √ m dependence, the contribution from lighter axions becomes quickly negligible for a fixed decay constant.
C. Joint PDF for the decay constant and the initial displacement We now add one more source of complexity in our calculation, by introducing a distribution for the decay constant f (f a ) and a possible correlation with the initial field amplitude ϕ in .Keeping the axion mass in the range −33 ≤ log 10 (m/eV) ≤ −29 and expressing ϕ in in units of M Pl , the total birefringence angle is given by where z is a random variable uniformly distributed between [−1, 1].
We want to focus on the size of the initial field value rather than the sign, especially when we consider correlations, thus we take as marginal distribution for ϕ in a half-normal distribution which is defined only for positive values, whereas for f a we choose a log-normal distribution: with corresponding average values ⟨ϕ in ⟩ = µ ϕ and ⟨f a ⟩ = exp (µ a + σ 2 a /2) 6 .Despite focusing on these marginalised distributions, we checked that the resulting distribution of the total birefringence angle does not change significantly if we take the distribution of the decay constant uniform in logarithmic space 7 , as shown in Figure 2. Of course, this result depends on the dispersion of f a and, for increasing values of σ a , the result starts diverging for the normal and uniform distribution as can be seen in Figure 16  Right: The birefringence PDF for the log-normal distribution with different values of σ a .We present the curves as smooth instead of histograms to make them visually easier to compare.
The first important result is shown in Figure 3, which displays the number of axions needed to match β obs ∼ 0.3 deg as a function of the mean values of the distributions of (f a , ϕ in ) and the corresponding total abundance.We explain the details of the numerical simulations in Appendix C. The main message is that the number of axions is very sensitive to the mean values of (ϕ in , f a ), exhibiting an approximately quadratic dependence as follows from Eq. ( 23), and equal-N lines positively correlate the two parameters.The bottom-left corner of the parameter space, corresponding to the low decay constant and large initial field value, leads to N < 1, meaning that not even one axion is allowed in this regime.Note that for the statistical treatment to be valid one should demand N ≳ 10, which corresponds to the second contour line in Figure 3 and in the following ones.Interestingly, we find that the constraints on the abundance, whose dependence on the number of axions is shown on the right panel of Figure 3, translate into a constraint on the decay constant which is independent of the initial field value, as was previously found for the fixed f a case.Figures 4 and 5 show the same results when positive and negative correlations are introduced between the two input parameters.
Even in these cases, the results do not qualitatively change and the equal abundance contours remain straight, but their position shifts to the left or the right of the zero-correlation value by an O(1) factor, quantitatively changing the constraint for f a .We can understand this behaviour from analytic estimations of the propagation of the variance for a multivariable function with FIG.4: Same as Figure 3, but for positive correlation ρ = 0.9 between ϕ in and f a .
correlations that we review in Appendix A. In this case, the variance of β is8 : where we neglected terms proportional to ρ 2 as well as higher-order correlators 9 .It is clear that a positive correlation requires a greater number of axions in order to give σ β ∼ 0.3 deg compared to the uncorrelated case, which in turn gives a stronger upper bound on f a .The opposite is true in the case of negative correlation.These results are summarized in Table I.Note that the previous expansion is not accurate when the distribution has long tails, since one should include more terms in Eq. ( 40), thus the actual shape of the distribution matters and the result is not solely determined by the means and variances of the field amplitude and axion decay constant.In particular, we found that a broader spread of β leads to an overall shift of the expected N to lower values.This suppression relaxes the upper bound on the decay constant.
We can compare the previous results to the case of "aligned" axions where β scales linearly with the number of axions, thus a smaller value of N is needed for given initial field values and decay constant, as shown in Figure 6.In this case, the lines of constant Ω ϕ depend non-trivially on the initial value ϕ in and so does the maximum value of f a .The numerical details of the "aligned" case are also described in Appendix C. The corresponding total axion abundance.

D. Birefringence tomography
In this section we discuss the expectations for the birefringence tomography [50,51], i.e. the ratio of birefringence angles coming from reionization and recombination β rei /β rec , in the Axiverse scenario and compare it to the single-field case.In the single-field case, we can easily estimate this quantity using the analytical solution for the axion field evolution in the matter domination (MD) epoch ϕ(t) ∝ sin(mt)/mt [17], which gives This result is shown in Figure 7 as a function of the axion mass 10 .We also show the field evolution for four different axion masses.Notice that the dependence on the mass is rather strong when 10 Here we take t rei t 0 = a rei a 0

3/2
= (1 + zrei) −3/2 ≃ 10 −2 for zrei ∼ 8 and trec t 0 ≃ 10 −4 for zrec ∼ 1000.m ∼ H rei ∼ 20H 0 , which is why future detectors can use the difference in the birefringence between recombination and reionization to probe axions in the window 10 ≲ m/H 0 ≲ 100 [50].In the multiple-axion case, we can write where in the regime of interest ϕ rec ≃ ϕ in and we absorbed again the anomaly coefficient into f a,i .
To make contact with the result given in Section III B, notice that the contribution of axions with mt rei ≫ 1 is suppressed, simply meaning that the axions that have started oscillating before reionization do not contribute to β rei .From this and the fact that the average value of the birefringence is zero, we conclude that β rei /β rec ∼ σ rei β /σ rec β = 1/ √ 2 ≃ 0.7 when the axions have a uniform mass distribution.This result can change significantly in the presence of correlations.For example, a correlation between the mass and the decay constant, that we discussed in SectionII B11 , will weigh differently the contribution of axions rolling before and after reionization, enhancing the signal in one of the two regimes.As an example, for a strong negative correlation between the mass and the decay constant, the contribution of lighter axions is suppressed, which changes the variance of the distributions of β rec and β rei .This is shown in Figure 8 for ρ(m, f a ) = −0.9, which leads to smaller than what we expect for zero correlation.On the other hand, for a positive correlation, the contribution from more massive axions is suppressed leading to β rei /β rec ≃ 1 because most of the birefringence is generated at later times.Notice that even with this type of correlation, the birefringence ratio is typically between zero and one, similar to the single-field case, as can be seen from the probability density functions shown in Figure 9 for three different choices of ρ(m, f a ) and from the monotonic dependence of its mean with the increase of the correlation shown in Figure 10.
Interestingly, in the multi-field case, it is also possible that |β rei /β rec | > 1, but this requires a cancellation between fields that roll after and before reionization.This can happen in the "aligned" scenario with a correlation between the initial field value, randomly distributed around zero, and the mass.In this case, a positive correlation implies that more massive axions have preferably positive initial values whereas lighter ones have negative initial values and the opposite for a negative correlation.Because of this cancellation between positive and negative displacements, the mean values of the distributions of β rei and β rec change as shown in the left panel of Figure 8.The shift of the mean values of the distribution of β rec and β rei are related to each other as can be intuitively expected.Indeed, approximating the distribution of ϕ in and m with a Binormal with ρ ̸ = 0, the average value of the birefringence gets shifted proportionally to the correlation Inserting this into Eq.( 42), it is straightforward to see that the ratio of the birefringence angles does not depend on the correlation where ⟨m⟩ ≃ 10 2 H 0 for uniform distributions of masses within H 0 < m < 10 4 H 0 .This simple result is surprisingly accurate when compared to the numerical calculation of ⟨β rec ⟩/⟨β rei ⟩ as shown in Figure 10 for sufficiently large correlation, when the ratio is dominated by the mean and not the variance; this condition breaks down for ρ ∼ 0. Note that in Figure 10 we take ⟨β rec ⟩/⟨β rei ⟩ and not the inverse because for ⟨β rec ⟩ < ⟨β rei ⟩ the numerical error is smaller.It is worth noting that a similar situation was recently discussed in Ref. [69] for a two-axion model to motivate the search of birefringence originated at late-time through the recently proposed polarized Sunyaev Zel'dovich effect [70,71].This search is further motivated in our axiverse scenario where the signal is generated at all cosmic times.

V. AXION MONODROMY
We now move to the more complicated case of axion monodromy potentials [72][73][74][75][76][77][78][79][80].Our goal is to explore the parameter space of the Axiverse for this type of potential since altering the largefield dynamics of the rolling axion field can affect both the birefringence angle as well as the axion abundance.We model axion monodromy as  44).This is however only valid when the ratio is dominated by the mean of β and not by its variance, as is the case close to ρ = 0.
where p is a free parameter that controls the large-amplitude slope of the potential.The value p = 1/2 leads to the usual asymptotically linear potential, but values of p < 1/2 have appeared in the literature in the context of inflation [81].Furthermore, the mass-scale M defines the transition from a quadratic potential for |ϕ| < M to a potential that grows like |ϕ| 2p for |ϕ| > M .Choosing M close to ϕ in will largely reproduce the results shown in the previous section, as the dynamics will be dominated by the quadratic part of the potential.The pre-factor of the potential is chosen, such that m is the free-particle mass of the axion.Axion monodromy potentials of the form of Eq. ( 45) introduce two new parameters, compared to the quadratic case, the tranisiton scale M and the power-law exponent p.
We start by examining the axion abundance Ω ϕ as a function of M , ϕ in and p.The potential of Eq. ( 45) can be approximated as: Because at large field values, the potential is flatter than in the quadratic case, the onset of oscillations is delayed.This can be estimated as in Ref. [82], by comparing the two relevant timescales, that of cosmic expansion t H ≃ H −1 and that of the motion driven by the potential t ϕ = ω −1 ϕ .Defining ω ϕ ≃ |V ,ϕ /ϕ| leads to H osc ≃ 3ω ϕ and, by using the scaling a(t) ∝ t 2/3 in the MD era, we find that the field starts rolling, and subsequently oscillating, around: where φ = M (H/3m) 1 p−1 is the field value evaluated at that time.For the remainder of this section, we work with non-dimensional units t → H 0 t, m → m/H 0 and ϕ → ϕ/M , unless otherwise stated.
It is important to note that the parameter space of interest for having cosmic birefringence, i.e.
the field rolls significantly between recombination and today 10 4 < t osc < t 0 , does not depend only on the axion mass, as is the case for a quadratic potential, but is very sensitive to the initial field value and the power-law parameter p.In particular, the maximum axion mass that can contribute to the signal is lifted to higher values as m max ≃ 10 4 (ϕ in ) 1−p ≫ 10 4 (in units of H 0 ) thus opening the available parameter space to masses m > 10 4 .The cosmological evolution of the background field in these potentials is more complicated since at large field amplitude |ϕ| ≫ 1, the abundance redshifts as [83] ρ and later when the field only explores the quadratic region around the minimum, it dilutes as standard dark matter with w = 0.In what follows we focus on the usual p = 1/2 case with linear behaviour at large field values.We study in detail how the background field redshifts and discuss how the phenomenological constraints from birefringence change in this case.

A. Asymptotically linear potential
We now focus on the case of p = 1/2, where the potential can be approximated as: and we briefly reintroduced the mass scale M for clarity.Taking the initial condition as ϕ in ≫ 1, the axion field starts evolving in the linear potential where the analytic solution is [84] After the field reaches unity ϕ = 1 in a time-scale of a first "transition" of dynamical behaviour occurs.The subsequent evolution, in particular the damping of the oscillations, depends on the value of ϕ osc , which defines when the potential term becomes dominant compared to the cosmological term 3H(t), which acts as a time-dependent viscosity.
When ϕ in ≫ 1 the field abundance redshifts as in Eq. ( 48) with w = −1/3, until the field amplitude is sufficiently damped (ϕ < 1) and the abundance redshifts as w = 0. Therefore, depending on the initial condition, we can distinguish three types of behaviour: rolling in a linear potential, oscillating and probing the linear part of the potential and finally oscillating in a quadratic potential.We identify the transition between oscillation in linear and quadratic potential with the time t 1 corresponding to field amplitude of ϕ = 1 as discussed later.In what follows we identify the relative parameter space for each behaviour to occur: (i) Oscillations in the quadratic potential t tr > t osc : When the initial value is smaller than unity from the beginning, the field never probes the linear part of the monodromy potential.
Actually, this regime extends until ϕ in < 3M because at the time the field starts oscillating ϕ < M and the field evolution continues as in the quadratic case.This can be seen from the linear evolution given in Eq. ( 50) evaluated at the time when the oscillations start (47): leading to ϕ osc = ϕ in /3.The abundance is given by Eq. ( 8), as in the quadratic case.
(ii) Oscillations in the linear potential and no transition t osc < t tr < t 0 < t 1 : This happens when ϕ in > 1, thus the evolution of the axion abundance includes the rolling in the linear part of the axion potential and the subsequent oscillations, which also probe the linear part of the potential.The evolution of the abundance during the first part of the motion can be accurately computed from the analytic solution of Eq. ( 50) for t < t tr and from Eq. ( 48) for t > t tr . 12The final value of the abundance becomes (we neglect numbers smaller than O(1) compared to ϕ in ): During the evolution, the amplitude of the oscillations gets reduced and the field will eventually probe only the quadratic part of the potential.Numerically we find that ϕ(t 1 ) = 1 well approximates this transition.For p = 1/2 we find that the amplitude of the scalar field scales as the density from Eq. ( 48) ϕ amp (t) ∼ ϕ osc (t osc /t) 4/3 , which gives Given the values of m and ϕ in , this transition hasn't happened yet if t 1 > t 0 ∼ 1 (in units of (iii) Oscillations in the linear part and transition to quadratic potential t osc < t tr < t 1 < t 0 : Following the previous discussion, in this regime, the evolution of the abundance has three different scalings: rolling and oscillations in the linear potential and subsequent oscillations in the quadratic region.In Figure 11 we show that the analytical scaling in these three phases closely follows the numerical evolution and gives the current abundance Remarkably, as for the quadratic potential case, the final abundance does not depend on the mass of the field.
The regions of parameter space corresponding to the three different types of evolution are shown in the left panel of Figure 11.

Probability of rolling after recombination
A natural question one can ask is: how likely is it that the fields start rolling after recombination, i.e. 10 −4 < t osc < 1, given an initial probability distribution of ϕ in and m and how much does it change for the asymptotically linear potential?
In the case of a cosine or quadratic potential, the onset of oscillations depends solely on the mass (assuming that the initial field amplitude is not fine-tuned to be close to the maximum of the cosine potential), thus the number of active axions after recombination depends on how many populate the mass range 1 ≲ m/H 0 ≲ 10 4 .If we take the mass uniformly distributed in log-space over several orders of magnitude, e.g.log 10 (m/eV) ∈ [−33 log 10 (m min /H 0 ), 27], this is given by N tot × 4/60 for m min = H 0 .For the monodromy case, Eq. (47) shows that the onset of oscillations depends also on the initial distribution of ϕ in .In particular, the time of the onset of oscillations in logarithmic space is: thus its distribution is given by the convolution of the individual distributions of log(ϕ in ) and log(m).When the initial field has a distribution peaked at a certain scale, we found that the probability of having 10 −4 < t osc < 1 does not change much compared to the previous case, since the mass is distributed over many orders of magnitude O (4/(log 10 (m max /m min )).As can be seen from Figure 11, the masses of axions active between recombination and today are still around four orders of magnitude, but, compared to the quadratic case, the mean value gets shifted to higher masses as ϕ in /M increases.This means that the probability is not affected if the mass distribution is flat, but it is certainly altered if the distribution is peaked at some scale or if it is tilted towards higher or lower masses.
The main goal of this work is to use the birefringence signal and the DM fraction to determine the interesting parameter space of (ϕ in , f a , M, m) where the signal can be explained in the Axiverse scenario without conflicting with constraints on the current axion abundance.In Section V A we derived the relation between (ϕ in , M, m) and the present axion abundance, which gives an upper bound on the transition scale in the different regimes.On the other hand, the birefringence angle connects the field amplitude to the decay constant f a and from the astrophysical bound on f a we can infer a lower bound on M .To our knowledge, one of the tightest constraints on f a comes from the absence of evidence of spectral distortions of X-rays flux coming from quasars (e.g.Ref. [85]) which leads to f a ≳ 10 9 GeV.This bound is derived assuming one axion and not multiple light axions, as in our case.In the latter case, the photons can convert to all axions lighter than m a ≪ 10 −12 eV [85] which can be treated as massless.Following Ref. [65], we can define a single axion that couples to F µν and has an effective coupling g eff ϕγ = i g ϕγ,i .If the coupling of different axions is uncorrelated and is distributed around ⟨g ϕγ ⟩, one naively expects that the coupling is enhanced by a factor √ N 12 where N 12 is the fraction of axions with mass below the threshold mass m a ≪ 10 −12 eV 13 .We confirm numerically that taking the distribution of f a given in Eq. ( 38) the effective coupling scales as √ N14 , therefore we take g eff ϕγ = N dec × log 10 (m 12 /m min )α EM /2π⟨f a ⟩ where m 12 = 10 −12 eV and m min = 10 −33 eV, thus log 10 (m 12 /m min ) = 21.With this, we proceed with the scan of allowed values for the scale of the transition in the various regimes: (i) Quadratic regime (⟨ϕ in ⟩ ≪ M ): In this case, using Eq. ( 8) we get (ii) Linear regime: In this regime, there is a minimum initial value shown as the bottom corner of the green region of Figure 11 that is ⟨ϕ in ⟩ ≃ 10M for m ∼ H 0 , which gives the upper bound M ≲ 0.01M Pl from the abundance constraint.However, there is no minimum value for the transition scale that would come from a maximum value of ⟨ϕ in ⟩, as we can keep t osc constant when ϕ in increases by increasing m.Moreover, one can always adjust the transition scale to not exceed the abundance constraints, M 2 /M 2 Pl ≪ 1 as follows from Eq. ( 53).
(iii) Linear plus transition: The maximum transition scale for which the field starts in the linear regime is given by the abundance constraint when ⟨ϕ in ⟩ ∼ 3M , which gives √ 4N dec M/M Pl ≲ 1.On the other hand, the minimum transition scale can be related to the maximum initial value ⟨ϕ in ⟩ ≲ 10 6 M (corresponding to the upper corner of the orange region in Figure 11) and the birefringence through leading to 10 4 GeV ≲ M ≲ M Pl / √ 4N dec , where we used the lower bound for the decay constant.Eq. ( 61).As the initial field value increases, the mean value of the axion mass contributing also increases, as shown in Figure 14.

C. Implications from Cosmic Birefringence
In direct analogy to Figure 3 for the quadratic potential, we present the results for the number of axions and the corresponding abundance for the "linear plus transition" and "linear" cases.We continue using the half-normal distribution for the initial field value and the log-normal distribution for the axion decay constant, as in Eq. ( 38), and introduce a random sign for each axion, to capture the effects of positive or negative initial field amplitude.As we discussed in the previous section, the bounds on the transition scale M are different in the different regimes and so is the constraining region of (⟨ϕ in ⟩, ⟨f a ⟩).In Figures 12 and 13 we fix M = 10 10 GeV and show the results for the "linear plus transition" and "linear" case respectively.Note that the underlying contours for the numbers of axions in the (⟨ϕ in ⟩/M, ⟨f a ⟩/M ) plane are the same for both cases since we approximate ∆ϕ ≃ −ϕ in , but the constraints from the abundance change its final value in the two cases.This is also clear from the different abundance contours of the right panels of Figures 12 and 13.We can understand the slope of the constant Ω ϕ -lines from Eqs. ( 53) and ( 55) through substituting the scaling N ∼ (⟨f a ⟩/⟨ϕ in ⟩) 2 .The "linear plus transition" case is similar to the quadratic one, since the final abundance does not depend on the mass, thus we arrive at the following scaling Since the last factor of Eq.( 60) grows with ϕ in , equal ⟨Ω ϕ ⟩-lines show anti-correlation between ⟨ϕ in ⟩ and ⟨f a ⟩, as can be seen in Figure 12.On the other hand, the "linear" case is more complicated, because the final abundance, which scales as depends on the axion mass, whose mean value varies with the initial condition according to Eq. (47).
In Figure 14 we show the average value (m/H 0 ) 2/3 as a function of the initial condition 15 for three different values of the transition scale M .Note that the mean mass increases for decreasing M , because higher values of ϕ in are allowed (without exceeding Ω ϕ,max ) therefore more massive fields can contribute.Because of this dependence of the average mass on the initial condition, the final abundance shows a change of behaviour that shows up as a "knee" in Figure 13.After this "knee" the average value of the mass grows slower with ⟨ϕ in ⟩/M and the behaviour of the abundance is determined by the last factor in Eq. ( 61), which decreases with ⟨ϕ in ⟩/M .Then Ω ϕlines show a positive correlation between ⟨ϕ in ⟩ and ⟨f a ⟩, as opposed to the case with the transition.
Therefore, it is important to understand how different field evolution affects the degeneracy lines in the (⟨ϕ in ⟩, ⟨f a ⟩) parameter space and modifies the relative constraints.
We now examine the dependence of our results on the transition scale M .Note that the dynamics of the scalar field, discussed in Section V A, does not depend on its absolute value, but 15 We compute the mean value of the mass for corresponding initial value by constructing many realizations of the pairs (m/H0, ϕin/M ) taken from the green region of the allowed parameter space shown in Figure 11.For each ⟨ϕin/M ⟩ we select those values of m/H0 within the 3σ ϕ band, and construct a histogram of the corresponding masses that we then use for computing (m/H0) 2/3 .II.We see that changing the transition scale M by two orders of magnitude, mildly affects the other parameters.

VI. SUMMARY AND DISCUSSION
In this work, we explore the expectations and the implications for the multidimensional probability density functions of the axions' parameters as the initial field value, the decay constant and the mass in the Axiverse scenario.We consider three different types of axion potential: cosine, quadratic and asymptotically linear monodromy potential.In most cases, the average value of the birefringence angle vanishes, since axion contributions with different sign cancel each other out on average.Therefore, except in special "aligned" cases, we consider the standard deviation of the distribution of values of the birefringence angle as the characteristic value.
In the case of the usual cosine potential, we note that more than O(20) axions are needed to explain the observed value of the birefringence angle, thus a statistical treatment is justified.In the multidimensional parameter space resulting from string compactifications, correlations might appear and we studied the effect of those on the number of inferred axions.In the Gaussian approximation, we give analytical formulas for the variance of the birefringence angle with a correlation between the initial misalignment and the anomaly coefficient that can change the dependence on the number of axions from √ N to linear in N , where N is the number of axions that contribute to the birefringence signal.We also showed how to use future observational data, providing the ratio between birefringence from reionization and recombination, to probe the distribution of axion masses, given a model for the axion potential and couplings.In the case of a log-uniform distribution, there is a clear expectation of β rei /β rec ≃ 0.7.Deviations from this result can be parameterized using a tilt in the mass distribution, which can be compared to expectations in the context of the string Axiverse.
The case of a quadratic potential allows for significant analytical treatment since the axion evolution can be derived analytically and thus we can connect the final abundance to the initial field amplitude.Limits on the present abundance provide an effective maximum displacement since each field contribution adds up linearly as opposed to the variance of the birefringence.Because of the different dependence on N , the abundance provides an upper bound on the value of the decay constant.If we take the decay constant to be the same for all axions, (defined in Eq. ( 1)) This work represents an initial attempt of using the cosmic birefringence data to constrain the properties of the Axiverse.In order to keep the analysis tractable, we limited ourselves to specific potentials and considered a limited number of input random variables in each case.Despite the expectation of a universe filled with (pseudo)scalar fields, most analyses of axion phenomenology have focused either on single or few-field models (with important exceptions, e.g.Refs.[43,66], focusing on black hole superradiance).On the contrary, multi-field dynamics has been extensively applied to inflationary model-building [87][88][89][90][91][92], thus providing interesting parallels between the two areas of physics.A similar research program is necessary, in order to fully appreciate the implications of an axiverse (whether of string theory origin or not).Recently, machine learning methods have found applications in physics [93], specifically in explorations of the String landscape [94,95], as well as in reconstructing the inflationary potential from cosmological data [96].Building on these advances presents an interesting avenue for a more extensive evaluation of the multi-dimensional model space that Axiverse models possess, including axions with different potentials and correlations between several of the parameters.for the initial field value σ ϕ i = 10 −2 .In the second and third panels from the left, we can clearly see the difference between choosing a positive (ρ = 0.9) or a negative correlation (ρ = −0.9) between the two variables.
• The birefringence is computed using Eq. ( 37) for N axions over several realizations.For each point of the grid, we then compute the value of σ 2 β ∝ N from interpolating the variance of the resulting distribution of β for an increasing number of axions.
• We find the number of axions needed for each parameter combination (grid point) by requiring that σ β (µ a , µ ϕ , N ) = 0.3 deg.We then construct the contour plot coming from the 2D interpolation of N (µ a , µ ϕ , β obs ) over the grid.
• Concerning the abundance, we interpolate the mean value of Ω ϕ , for each µ ϕ , as a function of an increasing number of axions.Subsequently, we compute the corresponding abundance for each grid point by evaluating Ω ϕ (µ ϕ , µ a , N ).
In the "aligned" case we take β ∼ ⟨β⟩ ± σ β and instead of using the fit β ∼ N , we find the parameters of the more complicated fit We confirmed that this model recovers well the scaling of β in this case.The following steps are the background motion as a source of CMB birefringence B. Statistics of many ALP's III.Cosine potential A. Anomaly coefficient and correlations B. Testing the mass distribution with birefringence tomography IV.Quadratic potential A. Fixed decay constant f a B. Projected abundance to higher masses C. Joint PDF for the decay constant and the initial displacement D. Birefringence tomography V. Axion Monodromy A. Asymptotically linear potential 1. Probability of rolling after recombination B. Range of transition scale values for Monodromy potentials C. Implications from Cosmic Birefringence VI.Summary and discussion I. INTRODUCTION FIG.1:The projected abundance for each decade in mass given by Eq. (31).The maximum masses corresponding to the values f a /M Pl = {10 −2 , 10 −3 , 10 −4 } shown in the figure are m max = {10 −24 , 7 × 10 −20 , 10 −15 } eV.

FIG. 2 :
FIG. 2: Left: The Probability Density Function of the birefringence angle for the log-normal and log-uniform marginal distribution of the decay constant.In both cases we take N = 50, σ ϕ ≃ 10 −2 , µ a = 16.5 and σ a = 0.25 and the resulting standard deviation of the birefringence angle is σ β ∼ 0.3 deg.In red we show the Gaussian distribution with σ = 0.3 that captures the emergent birefringence distribution.

1 FIG. 3 :
FIG. 3: Left: Number of axions needed to saturate σ β ∼ 0.3 deg as a function of the µ ϕ and µ a parameters of the PDFs (38) in the case of zero correlation ρ = 0. Blank regions are excluded because N < 1 whereas the high f a region is excluded from abundance constraints.Right: Abundance of the total axions whose number comes from the left panel.The red dashed line corresponds to N = 1 coming from the left panel whereas the black line corresponds to Ω ϕ,max = 3 × 10 −3 .Note that the "waviness" of the lines in both panels can be attributed to sampling noise.

1 FIG. 6 :
FIG. 6:Left: Similar to Figure3, but for the "aligned" case with ρ = −0.8.Compared to the previous scenarios, fewer axions are needed and their dependence is very different with respect to µ ϕ and µ a .Right:

FIG. 7 :
FIG. 7: Left:The ratio β rei /β rec for a single rolling axion as a function of its mass, normalized by H 0 .Right: Four examples of rolling axions with m/H 0 = 10, 140, 190, 430 (green, blue, red and black respectively).The rolling before and after reionization is shown in solid and dashed respectively.The units for the field value are arbitrary, since the evolution does not depend on the initial field amplitude for a quadratic potential.

FIG. 8 :
FIG. 8:The PDF of the birefringence angle from recombination and reionization in the presence of different types of correlations.Left: An example of a large positive correlation between the mass and the initial field value.This leads to |⟨β rei ⟩| > |⟨β rec ⟩| because of the partial cancellation of the birefringence from recombination from axions with different masses.Right: An example of a large negative correlation between the mass and the decay constant (with a log-normal distribution).This mainly changes the variance of the distribution and not its mean value as in the previous case.

FIG. 10 :
FIG. 10: Left: The monotonic dependence of σ rei β /σ rec β on the correlation between the mass and the decay constant.Right: We show the mean value ⟨β rei ⟩ and ⟨β rec ⟩ as a function of the correlation between the mass and the initial value.Contrary to the left panel the angle from reionization is bigger in amplitude than that from recombination, therefore, we take ⟨β rec ⟩/⟨β rei ⟩.At large correlation, this oscillates around −0.4 as estimated by Eq. (44).This is however only valid when the ratio is dominated by the mean of β

FIG. 11 :
FIG. 11: Left: Regions of the parameter space of the initial values and masses corresponding to the different evolution for the asymptotic linear potential explained in Section V A with 10 −4 < t osc < 1.The green region corresponds to fields that oscillate exploring the large field values and the transition from oscillation in linear to quadratic potential hasn't occurred yet.Fields in the orange region experience such a transition within the present time and are oscillating in the quadratic region around the minimum of the potential.Evolution in the blue region is quadratic corresponding to ϕ in ≤ 3M .Right: An example of the background evolution of an axion field that experiences a transition from the oscillations in the linear to the quadratic regime of the potential.The chosen values are m = 10 4 and the time-scales t osc = 10 −2.6 and t 1 = 0.1 are denoted by the dashed red and black vertical lines respectively.
leading to the upper bound M < 0.1M Pl which gets reduced by 1/ √ 4N dec in the case of multiple axions.Inverting Eq. (40) for the variance of β leads to the scaling √ 4N dec ⟨ϕ in ⟩ ≃ 10⟨f a ⟩, thus reintroducing the transition scale ⟨ϕ in ⟩ ≃ 10⟨f a ⟩/ √ 4N dec < M and the upper bound on the decay constant ⟨f a ⟩/ √ 21N dec ≳ 10 9 GeV leads to 10 10 GeV < M < 0.1M Pl / 4N dec .

2 FIG. 12 :
FIG. 12: Left: Number of axions needed to achieve σ β = 0.3 deg as a function of µ ϕ and µ a and corresponding excluded region for the "linear plus transition" case, i.e. orange region in Figure 11.The transition scale is fixed to M = 10 10 GeV.Compared to the quadratic case the maximum initial field value is lower by an order of magnitude ⟨ϕ in ⟩ ≲ 10 16 GeV.Right: Corresponding total abundance for the number of axions coming from the left panel.Note the negative correlation between the initial field value and decay constant of regions with constant abundance.

1 FIG. 13 :
FIG.13: Analogous to Figure12for initial conditions leading to just linear evolution, i.e. the green region in Figure11.Note that in this regime, the abundance contours present a change of behaviour, similar to a "knee", because of the dependence of the mass on the formula of the final axion abundance in

3 FIG. 14 :
FIG. 14: Average value of (m/H 0 ) 2/3 as a function of the mean initial field value ⟨ϕ in /M ⟩ for three selected values of the transition scale M = 10 12 , 10 10 , 10 8 GeV.It is evident that lowering the transition scale increases the average value of the mass of the axions that become active after recombination.

8 FIG. 15 :
FIG. 15: The analogous plots to the left panels of Figures 12 (left) and 13 (right) of the "linear plus transition" and "linear" regime for a transition scale of M = 10 12 GeV.

9 FIG. 17 :
FIG. 17: Joint PDF of the decay constant, whose marginal is a Gaussian distribution in log-space, and initial condition, whose marginal is a half Gaussian distribution.The marginal distribution of the decay constant is taken with log 10 (⟨f a ⟩/GeV) = 16 and dispersion σ a = 0.25,