Moduli Stabilisation and the Statistics of Axion Physics in the Landscape

String theory realisations of the QCD axion are often said to belong to the anthropic window where the decay constant is around the GUT scale and the initial misalignment angle has to be tuned close to zero. In this paper we revisit this statement by studying the statistics of axion physics in the string landscape. We take moduli stabilisation properly into account since the stabilisation of the saxions is crucial to determine the physical properties of the corresponding axionic partners. We focus on the model-independent case of closed string axions in type IIB flux compactifications and find that their decay constants and mass spectrum feature a logarithmic, instead of a power-law, distribution. In the regime where the effective field theory is under control, most of these closed string axions are ultra-light axion-like particles, while axions associated to blow-up modes can naturally play the role of the QCD axion. Hence, the number of type IIB flux vacua with a closed string QCD axion with an intermediate scale decay constant and a natural value of the misalignment angle is only logarithmically suppressed. In a recent paper we found that this correlates also with a logarithmic distribution of the supersymmetry breaking scale, providing the intriguing indication that most, if not all, of the phenomenologically interesting quantities in the string landscape might feature a logarithmic distribution.


Introduction
The Peccei-Quinn mechanism is without any doubt the most elegant solution to the strong CP problem. It postulates the existence of an anomalous global U (1) PQ symmetry which is spontaneously broken at f a . The corresponding Goldstone boson is the so-called QCD axion a which enjoys a continuous shift symmetry. QCD instantons lift the axionic direction and provide a minimum where CP is conserved. The QCD axion develops a mass of order m a ∼ Λ 2 QCD /f a and naturally contributes to the dark matter (DM) abundance. The phenomenologically allowed window for the axion decay constant f a is given by 10 9 GeV f a 10 12 GeV, where the lower bound is due to astrophysical and direct observations while the upper bound comes from the requirement to avoid DM overproduction if the initial misalignment angle takes natural O (1) values.
This scenario for the solution of the strong CP problem relies on some assumptions which have to be checked in a UV complete embedding. Some crucial questions which need to be answered are: (i) What is the origin of the axion shift symmetry?; (ii) What dynamics breaks U (1) PQ spontaneously and sets the value of f a ?; (iii) Is f a related to other important physical quantities like the Planck scale M p , the string scale M s , the GUT scale M GU T , the Kaluza-Klein scale M KK or the scale of soft supersymmetry breaking terms M soft ?; (iv) what dynamics breaks U (1) PQ explicitly and sets the value of m a ?; (v) Is m a generated by QCD instantons or by other effects?: (vi) How many axion-like particles (ALPs) can arise from UV physics?; (vii) What is the allowed range of f a and m a for these ALPs?
Several studies performed during the last 15 years revealed that string theory can provide a successful answer to many, if not all, of the previous questions [1][2][3][4]. However string theory admits a plethora of 4D solutions which goes under the name of string landscape. Even if all 4D string vacua share some generic features about axion physics, the number of axions and the corresponding values of f a and m a take different values in different string vacua. In order to make contact with observations, it is therefore crucial to perform a statistical analysis of the distribution in the string landscape of phenomenologically relevant quantities like f a and m a which determine the axion DM abundance.
As stressed in our recent paper [5] where we derived the distribution of the supersymmetry breaking scale in the string landscape, these statistical studies need to be based on a solid understanding of moduli stabilisation. In the case of axion physics, the motivation is the following. We shall focus on the type IIB flux compactifications which provide a well-defined subset of the string landscape. A model-independent origin of 4D axions is provided by the higher-dimensional gauge form C 4 which gives rise to pseudoscalars with a continuous shift symmetry when reduced on internal 4-cycles Σ i 4 : θ i = Σ i 4 C 4 . These axions are the imaginary parts of the Kähler moduli T i = τ i +i θ i whose real parts τ i control the volume of Σ i 4 in string units. Due to a combination of supersymmetry, scale invariance and the axionic shift symmetry, at tree-level each T i is a flat direction [6]. The axionic directions θ i are lifted by instantons which preserve only a discrete shift symmetry. On the other hand, the saxions τ i can be stabilised either at perturbative or at non-perturbative level. Let us comment on the implications of these two situations for axion physics: • If a given saxion τ is fixed by non-perturbative physics as in KKLT models [7], the stabilisation is at leading order supersymmetric, implying m θ ∼ m τ ∼ m 3/2 . Given that the absence of any cosmological moduli problem requires m τ O(50) TeV [8] and m 3/2 sets the mass of the superpartners which cannot be lower than the TeVscale, in this case the axion θ is generically very heavy, and so cannot play the role of the QCD axion [4].
• If τ is stabilised by perturbative physics (such as α and/or string loop corrections to the Kähler potential), at this level of approximation m τ ∼ m 3/2 while m θ = 0.
In the regime where the effective field theory (EFT) is under control, i.e. where non-perturbative contributions are exponentially suppressed with respect to pertur-bative terms, instanton effects will lift θ while inducing negligible corrections to the stabilisation of τ . This would produce the mass hierarchy m θ m τ ∼ m 3/2 which identifies θ as a promising QCD axion candidate with θ a/f a .
Besides focusing on models where τ is fixed at perturbative level, the other conditions to be checked to get a viable QCD axion are that θ couples to the QCD sector coming from stacks of D7-branes, 1 and that stringy instantons generate a mass m θ,str for θ which is smaller than the one developed by QCD instantons, i.e. m θ,str Λ 2 QCD /f a . If all these conditions are satisfied, one has still to derive the value of f a which determines all the main phenomenological properties of the QCD axion: its mass, its couplings and its contribution to the DM abundance. Depending on the topology of the 4-cycle Σ 4 , f a can be either of order M KK for bulk cycles, or of order M s for blow-up modes [2,4]. In a given moduli stabilisation framework, these two fundamental scales can be explicitly written down in terms of the underlying parameters (like the string coupling g s and the vacuum expectation value of the tree-level superpotential W 0 ) which depend on flux quanta. By exploiting the known distributions of g s and W 0 in the flux landscape [10], one can therefore derive the distribution of f a . We shall perform this analysis by focusing on the Large Volume Scenario (LVS) [11,12] which fixes some Kähler moduli at perturbative level, and find that f a features a logarithmic distribution.
We shall consider two possible realisations of the QCD axion: (i) axions associated to blow-up modes, and (ii) axions associated to bulk cycles. In case (i) f a is independent on the Standard Model (SM) gauge coupling, while in case (ii) the decay constant is fixed around the GUT scale by the requirement of reproducing the observed visible sector gauge coupling α SM . Hence, once we focus on phenomenologically relevant vacua with α −1 SM ∼ O(10-100), only axions associated to blow-up modes feature a logarithmic distribution of f a . However we consider this to be the generic situation for vacua where the EFT is under control since axions from bulk cycles require an anisotropic shape of the extra dimensions which corresponds to a tuned situation for moduli stabilisation. The reason is the interplay between two conflicting conditions: the low-energy 4D EFT can be trusted only for large values of the internal Calabi-Yau (CY) manifold, while α −1 SM ∼ O(10-100) implies that the 4-cycle supporting the SM brane system cannot be too large.
This result confirms the naive expectation that a generic 4D string model is characterised by a QCD axion with a GUT scale decay constant which would overproduce DM if the initial misalignment angle θ in is not tuned close to zero. However it also shows that string vacua with a QCD axion with an intermediate scale f a and an O(1) value of θ in are not so rare since the number of flux vacua grows with f a only as a logarithm, instead of a power-law.
Interestingly, in [5] we found that also the distribution of the gravitino mass in the flux landscape is logarithmic, 2 providing the intriguing indication that most, if not all, of the 1 Notice that when the QCD sector lives on D3-branes at singularities, C4-axions are eaten up by anomalous U (1) and the QCD axion arises from open string modes [9]. 2 To be more precise, in [5] we concluded that the distribution of the gravitino mass is logarithmic in LVS models and power-law in KKLT scenarios. However recent explicit constructions of KKLT models [13][14][15] might indicate a logarithmic distribution also for the KKLT case.
phenomenologically interesting quantities in the string landscape might feature a logarithmic distribution. Once the distributions of more than one phenomenological quantity are known, it is important to look at potential correlations among them. In our case, we find that vacua with an intermediate scale f a are also characterised by TeV-scale soft-terms, as typical of LVS models [16]. Let us finally mention that a generic CY gives rise to many Kähler moduli in the 4D EFT. If several of them are stabilised by perturbative effects, only one of them will play the role of the QCD axion while all the others would behave as ALPs which tend to be ultra-light in the regime where the computational control over the EFT is solid. These ALPs have interesting applications to DM [17], dark radiation [18][19][20][21] and astrophysics [22,23]. We shall therefore derive also the distribution in the flux landscape of the decay constants, the mass spectrum and the DM contribution of stringy ALPs, finding again a logarithmic dependence.
One may wonder whether our findings provide a trustable representation of the generic situation for axion physics in the flux landscape since they are based on the LVS framework while other moduli stabilisation mechanisms at perturbative level have been proposed [24]. However recent studies of the Kähler cone of CY manifolds with a large number of Kähler moduli h 1,1 revealed that, in the regime where the volume of each holomorphic curve is larger than the string scale so that the α expansion is under control, the overall volume in string units grows as V (h 1,1 ) 7 [25]. This clearly implies that for a generic CY with h 1,1 ∼ O(100), the EFT can be under control only if the moduli are fixed at V O (10 14 ). Given that only LVS models yield an exponentially large CY volume which can naturally account for such a large value of V, we believe that the genericity of our results is rather robust. Ref. [26] presented an explicit LVS moduli stabilisation procedure which can lead to an exponentially large CY volume for arbitrarily large h 1,1 exploiting instantons on del Pezzo divisors and O(α 3 ) corrections at O(F 2 ) and O(F 4 ) (where F denotes an F-term). This moduli stabilisation scenario leads to several ultra-light axions in agreement with the expectation of [25].
Moreover, ref. [27,28] derived the distributions of the axion decay constants and masses for different values of h 1,1 but at a given point in the moduli space, focusing in particular on the tip of the so-called stretched Kähler cone, i.e. the point closest to the origin which allows to keep the EFT under control. Interestingly, they found that the mean value of f a decreases as h 1,1 increases. Our results are complementary to the ones of [27,28] since we included moduli stabilisation and worked out the distribution of f a and m a as a function of flux quanta, i.e. moving in the moduli space at fixed h 1,1 . The results of [27,28] can be integrated with ours since they provide the boundaries of the region in moduli space where the EFT is under control and our logarithmic distributions can be trusted, i.e. our logarithmic distributions are valid for f a f a,max (h 1,1 ) or m a m a,max (h 1,1 ) with f a,max (h 1,1 ) and m a,max (h 1,1 ) as given in [27,28] as a function of the number of Kähler moduli h 1,1 .
Similar considerations apply to the comparison of our findings with the ones of [29] which noticed that, in the presence of N 1 ALPs which are effectively massless, there is just a linear combination of them which couples to photons. Ref. [29] derived the distribution of the corresponding ALP-photon coupling g aγγ as a function of N (with N ∼ h 1,1 ) at a fixed point in moduli space, choosing again the tip of the stretched Kähler cone. For type IIB flux vacua, they found g aγγ (N ) ∼ 10 −21 N 4 GeV −1 which, according to our previous considerations, can be considered as a lower bound for a logarithmic distribution of g aγγ as a function of different flux vacua at fixed N , when moduli stabilisation is taken into account along the lines of our paper. This paper is organised as follows. In Sec. 2 we discuss in depth the interplay between axion physics and moduli stabilisation. We first describe in detail an example with h 1,1 = 4 where the QCD axion can arise from either a bulk or a blow-up cycle, and then we discuss a more general example with arbitrarily large h 1,1 . In Sec. 3 we derive the distribution in the type IIB flux landscape of several quantities of axion physics relevant for phenomenology: decay constants, masses, DM abundance, axion couplings to gauge bosons and axion dark radiation in Fibre Inflation models [26,[30][31][32][33][34][35][36]. We discuss our results and present our conclusions in Sec. 4. Three appendices are devoted to provide technical details: App. A gives the details of the axion canonical normalisation; App. B provides a few benchmark points which reproduce the observed fuzzy DM abundance for ultra-light stringy ALPs; App. C shows the distribution of additional quantities relevant for phenomenology, like moduli masses and the reheating temperature from moduli decay, which also feature a logarithmic distribution in the flux landscape.

Axions and moduli stabilisation
As explained in Sec. 1, axions can be light (i.e. much lighter than the gravitino and the soft terms) only if supersymmetry is broken and the corresponding saxions are fixed at perturbative level. Moreover models with a large number of Kähler moduli require a huge CY volume to keep control over the α expansion. These two considerations single out type IIB LVS models as the best framework to study the interplay between axion physics and moduli stabilisation.
We shall now describe moduli stabilisation for a toy-model which can feature up to 3 light axions. This model is, at the same time, simple enough to perform moduli stabilisation in full detail, and rich enough to be a good representative of a more generic situation. In fact, it has 1 axion which becomes as heavy as the gravitino because of non-perturbative stabilisation, 1 ultra-light bulk axion which plays the role of an ALP, and 2 QCD axion candidates arising from the reduction of C 4 over a bulk or a local 4-cycle.

The geometry
The total number of Kähler moduli is h 1,1 (X) = 4 and the CY X features a K3 or T 4 divisor D 1 fibred over a P 1 base contained in a second divisor D 2 , and two additional rigid divisors D 3 and D 4 with only self-intersections. The Kähler form can be expanded in a basis of (1, 1)-forms as J = t 1D1 + t 2D2 − t 3D3 − t 4D4 where the t i are 2-cycle volumes and the negative signs have been chosen to ensure that all 2-cycle volumes are positive (in particular those dual to rigid divisors). The only non-vanishing intersection numbers are k 122 , k 333 and k 444 . Explicit examples with these properties can be found in [31,32,37]. Thus the CY volume form looks like: The 4-cycle moduli τ i = 1 2 XD i ∧ J ∧ J become: These relations can be inverted and V can be written in terms of 4-cycle moduli as: k 444 . Before dwelling on the details of moduli stabilisation, let us outline the main features of this representative model. We assume that the SM is built on stacks of D7-branes wrapping a 4-cycle in the geometric regime. As typical of LVS models, the internal volume is stabilised at exponentially large values. On the other hand, the 2 blow-up modes τ 3 and τ 4 are fixed at small values, and so the volume can be approximated as V α √ τ 1 τ 2 . Given that V is controlled by 2 moduli, we can consider 2 different regimes in moduli space: 1. Isotropic limit with SM on a local cycle: In this case τ 1 ∼ τ 2 τ 3 ∼ τ 4 . Both τ 1 and τ 2 are exponentially large, and so none of them can be wrapped by the SM D7-stack since the corresponding gauge coupling α −1 SM = τ i , i = 1, 2, would be hyper-weak. Hence the SM lives on a D7-stack wrapping the local divisor D 3 . τ 4 and θ 4 are fixed by instantons which make both of them as heavy as the gravitino. τ 1 and τ 2 are fixed by a combination of α and g s effects, and so the corresponding axions θ 1 and θ 2 are ultra-light ALPs. τ 3 is stabilised by a combination of D-terms, F-terms of matter fields and string loops. The associated axion θ 3 plays the role of the QCD axion with a decay constant of order M s which is around the intermediate scale for TeV-scale soft terms.
2. Anisotropic limit with SM on a bulk cycle: In this case τ 2 τ 1 ∼ τ 3 ∼ τ 4 . τ 1 and τ 2 are again frozen by perturbative corrections to the Kähler potential, and τ 3 by non-perturbative contributions to the superpotential. Contrary to the previous case, τ 3 is instead stabilised by non-perturbative physics. Given that τ 1 is hierarchically smaller than τ 2 , the underlying CY has an anisotropic shape with 2 extra dimensions much larger than the other 4. Thus the SM can live on the bulk divisor D 1 . θ 1 becomes the QCD axion with a decay constant set by the Kaluza-Klein scale associated to the fibre divisor D 1 which turns out to be of order the GUT scale. The mass of θ 3 and θ 4 is around m 3/2 , whereas θ 1 plays again the role of an ultra-light ALP.

Moduli stabilisation: leading results
The model-independent closed string moduli involve the axion dilaton S = e −φ + iC 0 , h 1,2 (X) complex structure moduli U a , and h 1,1 (X) = 4 Kähler moduli T i = τ i + i θ i . Dimensional reduction yields the following tree-level Kähler potential: 3 (2.4) S and the U -moduli are fixed at tree-level by turning on the 3-form flux G 3 = F 3 + iSH 3 which generates the superpotential: where Ω(U ) is the U -dependent holomorphic (3, 0)-form of X. The complex structure moduli and the axio-dilaton develop a mass of order m 3/2 , and so the associated axions are too heavy to be relevant for low-energy phenomenology.
Because of the no-scale cancellation, the T -moduli remain flat at semi-classical order. These directions are lifted by including non-perturbative corrections to (2.5) and perturbative corrections to (2.4). Focusing just on the Kähler sector and including only the leading order α effects and instanton contributions, K and W become: 3 with χ(X) the Euler number of X and ζ(3) 1.2, W 0 is the vacuum expectation value of (2.5), A 4 ∼ O(1) and a 4 = 2π for an ED3 wrapping D 4 , while a 4 = 2π/n 4 for gaugino condensation on a stack of n 4 D7-branes on D 4 .
Plugging (2.6) into the standard form of the 4D N = 1 supergravity F-term scalar potential, we end up with (up to an overall S and U -dependent factor): Minimising (2.7) with respect to the 3 moduli V, τ 4 and θ 4 results in: (2.8) This vacuum is AdS but there exist several mechanisms to uplift it to a dS solution [38][39][40][41]. The order of magnitude of the induced moduli masses is: showing that the axion θ 4 becomes too heavy to be relevant for low-energy phenomenology since m 3/2 sets also the order of magnitude of the soft terms, M soft m 3/2 , which cannot be below the TeV-scale. At this level of approximation all the other Kähler moduli are still flat.

Moduli stabilisation: subleading results and axion physics
Let us now describe the stabilisation of the remaining Kähler moduli by including additional contributions to K and W which are subdominant with respect to those considered in Sec. 2.2. The small parameters controlling the g s and α expansions are respectively e φ 1 and V −1/3 1. We shall consider the isotropic and anisotropic limits separately.

Isotropic limit with SM on a local cycle
In this case the SM lives on D7-branes wrapped around the 'small' rigid divisor D 3 . Because of the well-known tension between chirality and non-perturbative effects [42], τ 3 cannot be stabilised by instantons, and so θ 3 , contrary to θ 4 , remains light and can play the role of the QCD axion. Let us see this important issue in detail.

Moduli stabilisation
The total world-volume fluxes on the SM D7-stack and an ED3 instanton (similar considerations apply to gaugino condensation) on D 3 look like: where f SM ∈ Z and the half-integer contributions are due to Freed-Witten anomaly cancellation on non-spin divisors [43,44]. In order to obtain an O(1) instanton which contributes to W , the B-field has to chosen as B = 1 2D 3 so that F ED3 = 0. This, in turn, gives: The number of chiral intersections between the ED3 and the SM D7-stack is then given by: These zero modes can kill the ED3 contribution to W if they develop vanishing vacuum expectation values, as expected for visible sector fields in order not to break any of the SM gauge symmetries at high energies. In fact, the gauge flux F SM in (2.11) induces the following U (1)-charge for the modulus T 3 : and W has to be gauge invariant. This implies that the prefactor of the non-perturbative W has also to depend on charged matter fields. Considering for simplicity just a single open string field φ, the relevant U (1) transformations are: (2.14) Thus the non-perturbative superpotential (including the possibility of gaugino condensation): transforms under the anomalous U (1) as: implying that W ED3 can be gauge invariant if n = q T 3 /(n 3 q φ ). As can be clearly seen from (2.15), A 3 = 0 if φ = 0 (for n > 0). This problem comes along with the following correlated issue. A non-zero gauge flux on the D7-stack generates also a moduli-dependent Fayet-Iliopoulos term of the form [45,46]: If φ = 0, a vanishing D-term potential requires ξ SM = 0 which, in turn, implies τ 3 → 0, causing the collapse of the divisor D 3 to a singularity. This shrinking can be avoided in 2 ways: (i) by considering a slightly different geometry where the 2 rigid divisors D 3 and D 4 intersect each other so that ξ SM = 0 would just fix τ 3 in terms of τ 4 ; (ii) by considering the case where φ is a SM gauge singlet (like a right handed sneutrino) which can develop a non-zero vacuum expectation value by D-term cancellation.
In what follows we shall focus on the option (ii) since in the case (i) the anomalous U (1) would become massive by eating up a combination of the θ 3 and θ 4 , leaving no light closed string axions to behave as the QCD axion.
The D-term potential reads (taking, without loss of generality, φ as a canonically normalised field): A vanishing D-term potential then fixes the open string field at: The anomalous U (1) becomes massive by eating up a combination of θ 3 and the phase θ φ of φ = |φ| e iθ φ . Its mass is given by: where: is the decay constant of the open string axion θ φ , while f cl is the decay constant of the closed string axion θ 3 . This last quantity can be derived from the kinetic terms: where a θ 3 f cl is the canonically normalised axion, implying: (2.23) -9 -Comparing (2.21) with (2.23), it is easy to see that f op f cl for τ 3 1, signaling that the combination of θ 3 and θ φ eaten up by the anomalous U (1) is mostly given by the open string axion θ φ since the largest contribution to M U (1) in (2.20) comes from f op [9]. Thus θ 3 survives in the low-energy theory and can play the role of the QCD axion a. The corresponding saxion τ 3 develops a potential via 2 effects: 1. The F-term potential of the matter field φ generated by supersymmetry breaking effects, after writing |φ| in terms of τ 3 using (2.19): 2. The potential generated by string loop corrections to the Kähler potential due to the exchange of Kaluza-Klein modes between the D7-stack wrapped around D 3 and O7-planes or D3-branes [47,48]: where c loop is expected to be an O(1−10) coefficient which depends on the U -moduli.
The potential V matter +V loop admits a minimum at τ 3 = c loop /c ∼ O(10) which reproduces the correct order of magnitude of the SM gauge coupling g −2 SM τ 3 . It can be easily checked that the saxion τ 3 develops a mass of order m 3/2 similarly to τ 4 and θ 4 .
Notice that T 3 -dependent instanton corrections to W as in (2.15) would generate a potential of the form (using (2.19) and setting θ φ = 0 and n 3 = 1): (2.26) For n = q T 3 /q φ ≥ 2 and 2πτ 3 1 (i.e. the limit where higher-order instanton corrections can be neglected), the potential (2.26) is exponentially suppressed with respect to V matter + V loop , and so it produces just a tiny shift of the minimum for τ 3 . On the other hand, it would generate a mass for θ 3 of order: which for n ≥ 2, m 3/2 O(10 10 ) GeV and τ 3 = α −1 SM 25 is always subdominant with respect to the contribution from QCD instantons m a Λ 2 QCD /f a for any possible f a M p . This guarantees that θ 3 is a good QCD axion candidate.
The only saxion which remains to be fixed is the fibre modulus τ 1 . This field develops a potential via string loop corrections which experience an 'extended no-scale cancellation' that suppresses them with respect to the leading α correction [49]. The resulting scalar potential for τ 1 reads [30]: 4 where A and B are flux-dependent parameters. The minimum of (2.28) is located at: This result has 3 important implications: 1. For α λ O(1) and g s O(0.1), τ 1 is roughly of the same order as τ 2 , implying that the CY volume is isotropic. Without loss of generality, in what follows we shall consider αλ 3/2 g 2 s = 1, i.e. τ 1 = τ 2 .
2. The scalar potential (2.28) scales as V −10/3 , and so for V 1 it is indeed suppressed with respect to the leading order LVS potential (2.7) which scales as V −3 .
3. For V 1 the SM cannot live either on τ 1 or on τ 2 since the resulting gauge coupling would be too small. Hence the SM has to be supported by τ 3 .
Notice that the axions θ 1 and θ 2 are lifted only by tiny non-perturbative corrections to the superpotential of the form: , and a i = 2π/n i for i = 1, 2. Given that τ 1 τ 2 1, these effects would make θ 1 and θ 2 2 ultra-light, i.e. almost massless, ALPs.

Mass spectrum and decay constants
The mass spectrum of the 3 moduli fixed at leading order, V, τ 4 and θ 4 has been given in (2.9). The mass of the remaining moduli turns out to be: Notice that α effects are under control when V τ 3/2 1 τ 3/2 2 10 3 since the corresponding expansion parameter is V −1/3 0.1. In this regime the 2 ALPs θ 1 and θ 2 are almost massless, m θ 1 ∼ 0 and m θ 2 ∼ 0, since their mass would turn out to be smaller than the present value of the Hubble constant, H 0 10 −33 eV, for n 1 = n 2 = 1. Larger values of n 1 and n 2 can however raise m θ 1 and m θ 2 above H 0 , with interesting application to fuzzy DM [17]. In what follows we shall therefore consider θ 1 and θ 2 as ultra-light.
The decay constants of the QCD axion θ 3 and the 2 ALPs θ 1 and θ 2 can be derived from canonical normalisation and take the generic form (see [4] and App. A.1): where λ i is the i-th eigenvalue of the Kähler metric and n i determines the periodicity of the cosine potential which enjoys a discrete shift symmetry (with n 3 = 1 for the QCD axion).
33) where the c i 's are moduli-independent coefficients: Notice that the QCD axion decay constant scales as the string scale, f a M s M p / √ V, while the decay constants of the 2 ultra-light ALPs behave as the Kaluza-Klein scale, The order of magnitude of all these mass scales is set by the overall volume V. An interesting regime in moduli space is the one where

Axion couplings to gauge bosons
Other quantities which are relevant for phenomenology are the couplings of the axions to the gauge fields of the visible and hidden sectors. We shall focus just on the couplings of the QCD axion and the 2 ultra-light ALPs which we will express in terms of the corresponding canonically normalised fields a 3 ≡ a, a 1 and a 2 (see App. A for the details of canonical normalisation). The visible sector lives on D 3 while the hidden sector involves 2 intersecting stacks of D7-branes wrapped around D 1 and D 2 . Knowing that the gauge kinetic function of each sector is given by the corresponding unnormalised Kähler modulus, and denoting the field strengths of the canonically normalised gauge bosons respectively as F vis , F 1 and F 2 , we obtain: Notice that the QCD axion a has a stronger than Planckian coupling to the visible gauge bosons while its coupling to hidden sector degrees of freedom is very suppressed. On the other hand, the 2 ALPs have a standard O(1/M p ) coupling to hidden gauge bosons but they are decoupled from the visible sector. These results are due to the combination of two effects: (i) the visible sector lives on a shrinkable del Pezzo D 3 which has no intersection with the bulk divisors D 1 and D 2 ; (ii) the axions θ 1 and θ 2 are in practice massless.

Anisotropic limit with SM on a bulk cycle
In LVS scenarios the SM can be realised on a stack of D7-branes wrapped around a bulk cycle only if the underlying geometry has an anisotropic shape. In this case the overall volume can be exponentially large in agreement with a SM gauge coupling which is not too small.

Moduli stabilisation
We focus on the case where the SM lives on the K3 or T 4 fibre D 1 . Hence the visible sector gauge coupling is given by large, the internal geometry needs to be anisotropic with 2 extra dimensions much larger than the other 4. This can be achieved via the following moduli stabilisation procedure: • V, τ 4 and θ 4 are stabilised as in (2.8) and the CY volume becomes exponentially large in string units.
• Given that D 3 is not wrapped by the SM D7-stack, the non-perturbative superpotential (2.15) is not suppressed anymore due to chiral intersections with visible sector states. Hence the freezing of τ 3 and θ 3 is completely similar to the stabilisation of τ 4 and θ 4 . Contrary to the isotropic scenario, in this case θ 3 acquires a mass of order m 3/2 and plays no role for low-energy physics.
• The fibre divisor τ 1 is stabilised by string loop corrections as in (2.29) but with αλ 3/2 = 1 and g s 1. 5 This results in the following hierarchy: (2.36) • T 1 -dependent non-perturbative corrections to W would be very suppressed due to: (i) the presence of chiral intersections as for T 3 -dependent instantons in the isotropic case; (ii) the fact that D 1 is a non-rigid cycle with extra fermionic zero modes which tend to kill instanton contributions. Therefore the closed string axion θ 1 is a perfect QCD axion candidate which becomes massive via standard QCD instantons.
• The remaining closed string axion θ 2 is an almost massless ALP which develops a tiny mass via non-perturbative corrections to W which are exponentially suppressed in terms of the large 4-cycle τ 2 . 5 To be more precise, we envisage a situation similar to the explicit CY cases discussed in [32] where the volume (neglecting blow-up modes) is V √ τ1τ2τ2. Due to the intersection between τ2 andτ2, the Fayet-Iliopoulos term induced by gauge fluxes fixes τ2 ∝τ2 for vanishing VEVs of open string fields. An appropriate combination of closed string axions is eaten up by the corresponding anomalous U (1). Substituting τ2 ∝τ2 in V, one obtains effectively the same expression that we are considering: V √ τ1τ2.

Mass spectrum and decay constants
The mass of V, τ 4 and θ 4 is again given by (2.9). The mass spectrum of the other moduli instead reads: The decay constants of the QCD axion θ 1 and the ALP θ 2 now become: Moreover, for g s 0.1 and τ 1 10, (2.36) yields τ 2 10 3 which, in turn, implies that the ALP θ 2 is ultra-light, i.e. m θ 2 ∼ 0. The decay constant of this ALP is set by the Kaluza-Klein scale of the effective 6D theory since: Interestingly, this scale is one order of magnitude above the gravitino mass since: The decay constant f θ 2 can take different values depending on the order of magnitude of τ 2 . Varying τ 2 corresponds to varying V which is mainly controlled by g s , as can be seen from (2.8). The string coupling affects also the relation (2.36) where however τ 1 has to remain fixed to get the right SM gauge coupling. This can be achieved by varying W 0 as well. Hence (2.8), (2.36) and τ 1 = α −1 SM ∼ O(10 − 100) imply an interesting relation between W 0 and g s (ignoring O(1) numerical factors): . We conclude that this scenario naturally predicts a QCD axion decay constant around the GUT scale, an almost massless ALP with decay constant around 10 14 GeV and supersymmetry at intermediate scales.

Axion couplings to gauge bosons
Let us now focus on the coupling of the canonically normalised QCD axion a 1 ≡ a and ALP a 2 to gauge bosons belonging to the visible sector on D 1 and hidden sectors on D 2 and D 3 . In fact, SM particles are not charged under the gauge symmetries of the D7-stack wrapping D 3 since D 3 does not intersect with D 1 . Similar considerations apply to D 4 , and so we shall ignore the possibility of a hidden sector on D 4 since it would have the same features of the hidden sector on D 3 . On the other hand, the SM degrees of freedom can be charged under the gauge group on D 2 since there is an intersection among D 1 and D 2 . However this would still be a hidden sector since τ 2 is a big cycle, and so the corresponding gauge coupling would be hyper-weak. Thus the relevant couplings are (see App. A for the details of canonical normalisation): where the µ i 's are O(1) numerical coefficients. Notice that the QCD axion a has a standard Planckian coupling to visible gauge bosons since it arises from a bulk cycle. On the other hand its coupling to the hidden gauge bosons on D 2 is V-suppressed, while a is essentially decoupled from the hidden sector on D 3 since (m a /m θ 3 ) 2 ∝ (Λ QCD /M p ) 4 10 −76 . The ALP a 2 features instead an O(1/M p ) coupling to hidden gauge bosons on D 2 but it is decoupled from the other sectors. These results are again due to the fact that a 2 is essentially massless and D 3 has no intersection with D 1 and D 2 .

An example with arbitrary h 1,1
A generic CY threefold is characterised by hundreds of Kähler moduli, i.e. h 1,1 ∼ O(100), and so one may wonder whether the axion physics of this more complicated case would display features similar to the ones of the relatively simple case with h 1,1 = 4 analysed above. As we have stresses, this depends on the details of moduli stabilisation. In this section we shall describe how to freeze all Kähler moduli for arbitrary h 1,1 following [26]. We will obtain an LVS vacuum where V can be taken large enough to trust the α expansion.
The only requirement on the geometry is the presence of 2 blow-up modes, the first, D SM , to host the SM and the second, D np , to support non-perturbative effects. This condition is not too restrictive since del Pezzo divisors arise very frequently in CY constructions. We shall therefore consider an internal volume of the form: As explained in Sec. 2.2, the leading contributions to the scalar potential in a large-V expansion arise from O(α 3 ) corrections to K and T np -dependent non-perturbative corrections to W as in (2.6), which stabilise τ np ∼ g −1 s , θ np ∼ π/a np and V ∼ e 1/gs . Similarly to the isotropic case studied in Sec. 2.3.1, the SM cycle τ SM is instead fixed by the interplay of D-terms, F-terms of matter fields and τ SM -dependent loop corrections. This stabilisation procedure ensures that the internal volume can be exponentially large while τ SM = α −1 SM ∼ O(10 − 100) can reproduce the observed value of the SM gauge coupling. Moreover the axion θ SM behaves as a perfect QCD axion candidate with a decay constant of order the string scale.
At this level of approximation, there are still (N − 1) saxionic and N axionic flat directions (without considering the QCD axion θ SM which we assume to be lifted by QCD instantons). All the (N − 1) flat saxions can be lifted by including subdominant α effects. In 10D the first higher derivative corrections which modify the 4D scalar potential upon dimensional reduction, arise at O(α 3 ) and scale as G 2 3 R 3 . In 4D they generate the term proportional to ξ in (2.7). Additional 10D O(α 3 ) terms scale as G 4 3 R 2 , G 6 3 R and G 8 3 , and they give rise in 4D to higher F-term contributions to the scalar potential which scale respectively as F 4 , F 6 and F 8 [50]. When the superspace derivative expansion is under control [51], these terms represent just negligible corrections to the LVS potential (2.7). However they can be the leading effects to lift any remaining flat direction. In particular, the form of O(α 3 ) F 4 corrections for an arbitrary CY X has been determined to be [50] (ignoring the dependence on del Pezzo moduli): where λ ∝ g −1/2 s is a positive coefficient [52] and the Π i 's are O(1) topological quantities which can be expressed in terms of the second Chern class c 2 as Π i = X c 2 ∧D i [50]. Notice that Π i ≥ 0 ∀i = 1, ..., N in a basis of the Kähler cone where t i ≥ 0. The total potential can thus be written schematically as: where we have highlighted the moduli-dependence of each contribution. Extremising with respect to the 2-cycle moduli, we obtain: Using t i τ i = 3V, it is easy to realise that t i ∂ t i V tot = 0 implies: Plugging this result in (2.46) we find: This relation fixes (N − 1) moduli in terms of one of them, say τ N , as: The positivity of λ and the Π j 's ensures that this is a well-behaved minimum [26]. Substituting (2.49) in (2.43), we obtain: Given that this stabilisation is purely perturbative, at this level of approximation N axions are still flat. They can be lifted by including non-perturbative corrections to W which however tend naturally to give rise to axion masses below the present Hubble constant for V (h 1,1 ) 7 10 14 . Let us stress that for such a large value of V the SM is naturally expected to be supported on a blow-up mode since matching τ * = α −1 SM ∼ O(10-100) for a bulk cycle τ * would need from (2.49) a very unnatural hierarchy between Π * and Π N of order 10 −8 for τ N ∼ V 2/3 ∼ 10 10 .

Mass spectrum and decay constants
The mass of the 3 moduli fixed at leading order, V, τ np and θ np is given in (2.9). The mass spectrum of the remaining moduli becomes: where all the ALPs θ i 's are essentially massless. The decay constants scale as in the isotropic case with h 1,1 = 4 analysed in Sec. 2.3.1:

Axion couplings
We assume that the SM can be realised with a stack of magnetised D7-branes wrapped around D SM . On the other hand, the 'big' divisors D i , i = 1, ..., N , can in principle host several hidden sectors. The coupling of the QCD axion and the N ultra-light ALPs to visible and hidden gauge bosons can be derived from the moduli-dependence of the corresponding gauge kinetic functions. Denoting the canonically normalised QCD axion as a, the ALPs as a i (the results of App. A can be easily generalised to the isotropic case with many bulk Kähler moduli), and the field strengths as F vis and F i , we end up with: where again λ SM , theλ i 's and theλ i 's are numerical O(1) coefficients. The coupling of the QCD axion a to visible gauge bosons is enhanced with respect to 1/M p while the coupling to hidden degrees of freedom is V-suppressed. This is due again to the fact that D SM is a shrinkable del Pezzo divisor and the vanishing of the mass of the N ALPs.

Statistics of axion physics in the flux landscape
Building on our previous results [5], we investigate the statistical distribution in the type IIB flux landscape of various quantities of axion physics which are phenomenologically interesting. To this end, we first express these quantities in terms of the microscopic parameters, as we did in Sec. 2 via moduli stabilisation, and we then exploit their distributions. In particular the shall consider the following distributions for the underlying flux-dependent parameters g s and W 0 , and the rank n of the condensing gauge group which generates non-perturbative corrections to the superpotential: • The distribution of the string coupling g s is taken to be uniform. This result was explicitly checked in [5] for rigid CY manifolds and is believed to hold for more general cases as well [53]. Hence in the following we shall take dN dg s where N is the number of flux vacua.
• Based on the seminal work [10], the tree-level superpotential W 0 is assumed to be uniformly distributed as a complex variable, resulting in dN |W 0 |d|W 0 |. Note that this distribution might be different in regions where |W 0 | is exponentially small since recent constructions of KKLT vacua obtained |W 0 | ∼ e −1/gs [13][14][15]. However, as explained in Sec. 1, KKLT vacua feature only heavy axions with a mass of order m 3/2 , and so we shall focus just on regions where |W 0 | ∼ O(1-10) where its distribution can be taken as uniform.
• The distribution of the rank of the condensing gauge group n in the type IIB flux landscape is still poorly understood. All globally consistent type IIB CY models which have been constructed so far feature contributions to the superpotential which arise from just gaugino condensation in a pure SO(8) sector (corresponding to n = 6) and ED3 instantons (with n = 1). It is therefore still unclear if an actual distribution of n exists. If so, we argue that it should scale as dN −n −r dn with r > 0, since the number of flux vacua N is expected to decrease when n increases as D7-tadpole cancellation is easier to satisfy for smaller values of n.

Axion decay constants
Let us start with the axion decay constants. After evaluating the decay constants at the minimum of the scalar potential, we compute their distributions in the flux landscape using the scaling of the number of vacua N with the underlying parameters g s , W 0 and n.

Isotropic limit
The axion decay constants in the isotropic limit are given in (2.33). Being exponentially large, the main quantity which controls their distribution is the overall volume V. Using (2.8), we can therefore approximate the axion decay constants as: Notice that, at leading order, the decay constants do not depend on W 0 . Hence we can vary them with respect to just g s and n 4 , obtaining: where f can be any of the 3 decay constants, f a , f θ 1 and f θ 2 , and in the last step we have introduced Planck units. Using dg s dN and dN −n −r 4 dn 4 with r > 0, (3.2) takes the form: wherec = c for f a andc = 4c/3 for f θ 1 and f θ 2 . Ignoring subdominant logarithmic effects, we therefore obtain that the distributions of the decay constants of both the QCD axion and the 2 ultra-light ALPs scale as: Let us stress 3 important points: 1. Isotropic models with the SM localised on a blow-up cycle feature only a mild logarithmic preference for higher values of the axion decay constants.
2. As can be seen from (3.3), for f M p and 0 < r ≤ 2, the distribution of f is driven mainly by the distribution of g s . Moreover the final result (3.4) is unchanged if the distribution of g s is taken to be power-law.
3. As can be seen again from (3.3), the unknown distribution of n 4 would start being important only for r > 2. However it is reassuring to notice that it would affect only the form of subleading logarithmic corrections to (3.4).

Anisotropic limit
For the anisotropic case the axion decay constants are given in (2.38). As already observed in Sec. 2.3.2, the QCD axion decay constant is fixed around the GUT scale, f a ∼ M GU T , by the need to match the correct SM gauge coupling. If this phenomenological condition is dropped, however f a would feature a logarithmic distribution as in (3.4). The same is true for the distribution of the decay constant of the ALP θ 2 . However in this case the phenomenological requirement α −1 SM = τ 1 O(10-100) still leaves some freedom to vary f θ 2 in the flux landscape since from (2.38) we have f θ 2 n 2 g 2 s M GU T . Notice that the ALP decay constant does not depend on W 0 which has however to respect the relation (2.41) to keep τ 1 = α −1 SM constant when g s is varied. Differentiating f θ 2 with respect to g s and n 2 we thus obtain: Using again dg s dN and dN −n −r 2 dn 2 with r > 0, (3.5) reduces to: For 0 < r ≤ 1, n 2 ≥ 1 and g s 1, the second term in brackets in (3.6) is always smaller than unity. This term could instead become larger than 1 for r > 1. However, for g s 0.1, this would require large values of n 2 which are hard to realise in explicit examples (the largest value obtained so far is n 2 = 6 for gaugino condensation in a pure SO(8) gauge theory which would however still yield a second term of O(1) for r ≤ 3). We shall therefore consider the term in brackets in (3.6) of order unity, and obtain: Let us make 2 important observations: 1. Interestingly, we obtained now a power-law distribution for the ALP decay constant whose scaling is however very similar to the logarithmic case due to the mild square root dependence.
2. The distribution (3.7) holds as long as W 0 can be tuned to satisfy the relation (2.41) which implies W 0 ∼ e −1/gs . In the absence of a dynamical mechanism which fixes the flux superpotential in terms of dilaton-dependendent non-perturbative effects, this relation would however not hold anymore when g s is taken very small. A good estimate for the lowest value of the ALP decay constant for which (3.7) still applies, can be obtained for g s 0.01 which would give f θ 2 10 12 GeV.
Model with arbitrary h 1,1 As shown in Sec. 2.4, the results of the isotropic case with the SM on a blow-up cycle can be generalised to models with an arbitrarily large number of Kähler moduli where all saxions can be explicitly stabilised by α corrections to the scalar potential. In this case the axion decay constants are given in (2.52), and they scale with the CY volume as in the isotropic case discussed above. Hence we expect again a logarithmic distribution in the type IIB flux landscape as in (3.4): Let us comment on the regime of validity of these distributions. They hold at fixed h 1,1 when moving in the Kähler moduli space by varying microscopic parameters like g s after the decay constants are written in terms of them thanks to moduli stabilisation. These results are complementary to the ones of [27,28] (see also [29] for qualitatively similar findings) which found an approximate log-normal distribution for the axion decay constants of a given CY model focusing at the tip of the stretched Kähler cone. Moreover they found that the mean value of the f θ i 's decreases when h 1,1 increases. Given that the tip of the stretched Kähler cone corresponds to the smallest values of the Kähler moduli which are compatible with a controlled α expansion, and the axion decay constants are inversely proportional to 4-cycle volumes, the values of the f θ i 's obtained by [27,28] represent the largest values of the axion decay constants compatible with a trustable EFT. These values would therefore provide an upper bound for the regime of validity of our logarithmic distributions (3.8) which can be integrated with the results of [27,28] to describe how the number of flux vacua changes has a function of both f θ i and h 1,1 .
As an illustrative example, we consider the distribution N (f, h 1,1 ) where f is the mean value of the axion decay constants. As derived in [25], the requirement to trust the α expansion implies that the volume of each 4-cycle grows with h 1,1 as τ i (h 1,1 ) 3 ∀i = 1, ..., h 1,1 (at least for basis elements obtained from generators of the cone of effective divisors). On the other side, as we have seen in Sec. 2, the axion decay constants scale as f θ i M p /τ i . Combining the two results gives a qualitative understanding of the fact the mean value f of the log-normal distribution found in [27,28] decreases as h 1,1 increases. Moreover, we can obtain an explicit estimate of the upper bound for our logarithmic distributions: where, for a given h 1,1 , f can take different values by moving in the stretched Kähler cone in a way compatible with moduli stabilisation, and f = f max at the tip. Hence we expect the following distribution for the number of type IIB flux vacua as a function of f and h 1,1 (see Fig. 1): (3.10)

Axion masses
Let us now compute the distribution of axion masses in the flux landscape. As in the previous section we first compute the differential of the masses and then use the known scaling of the parameters g s , W 0 and n in order to determine the distribution.

Isotropic limit
The mass spectrum of the isotropic case with SM on a blow-up cycle is summarised in (2.31). Using (3.3) which can be approximated as df a f a dN , the distribution of the On the other hand, the distributions of the masses of the two ultra-light ALPs can be easily derived by noticing that their masses can be expressed in terms of the corresponding decay constants as: This result implies: which yields the following distribution (neglecting subdominant logarithmic corrections): (3.14) Notice that we find again a logarithmic distribution for the mass of both the QCD axion and the 2 ultra-light ALPs. However (3.11) and (3.14) have a different sign, implying that, in the QCD axion case, the type IIB flux landscape has a mild logarithmic preference for low scale masses, while the ALP case features more vacua at large mass values.

Anisotropic limit
The masses for the anisotropic geometry with SM on the bulk divisor D 1 are summarised in (2.37). As already pointed out in the previous section, this model is strongly constrained by the requirement to match the observed SM coupling. This sets the QCD axion decay constant around the GUT scale, f a M GU T , without a distribution. Thus the QCD axion mass would also be fixed at m a Λ 2 QCD /f a 1 neV. The mass of the ALP θ 2 can instead take different values in the flux landscape with a distribution which is again logarithmic. This result can be easily inferred by first writing m θ 2 as in (3.12) and then differentiating as below: where we used (3.5) approximated as: Barring subleading logarithmic effects, (3.15) therefore implies: This distribution, similarly to the one of f θ 2 derived in (3.7), holds as long as W 0 can be tuned to satisfy the condition (2.41) which keeps τ 1 = α −1 SM fixed at the right SM gauge coupling. In the case of the ALP decay constant, we estimated that its distribution is valid for f θ 2 10 12 GeV. Using (3.12), this gives however a lower bound for the ALP mass which can be safely ignored since it would be much smaller than today's Hubble constant: Model with arbitrary h 1,1 The mass spectrum for the model with a generic number of Kähler moduli is given in (2.51). Following the discussion of the distribution of the axion decay constants, the results for the distribution of the axion masses for the model with arbitrary h 1,1 would again be qualitatively similar to the isotropic case. Hence we expect logarithmic distributions of the form: As for the case of the axion decay constants discussed above, the results of [27,28] can be combined with ours to give an upper bound for the regime of validity of the distributions of the ALP masses as a function of h 1,1 , i.e. m θ i m max θ i (h 1,1 ).

Dark matter abundance
Let us now study the distribution of the axion DM abundance produced via the standard misalignment mechanism. We distinguish between the case where the DM particle is the QCD axion and the case where it is an ultra-light ALP. In the QCD axion case, the DM abundance is given by: while for the case of an ALP θ i , it reads: (3.20) For natural O(π) values of the initial misalignment angles θ in and θ i,in , the QCD axion can reproduce the observed DM adundance for f a 10 11 GeV, while an ALP would require m θ i 5 · 10 −21 eV for f θ i 10 16 GeV (see App. B for a detailed scan through the underlying parameter space and some benchmark points).

Isotropic limit and model with arbitrary h 1,1
In the isotropic case with SM on a blow-up cycle, the distribution of QCD axion DM abundance can be computed deriving (3.19) with respect to f a and then using the result (3.3) which, at first approximation, implies df a f a dN . Hence we end up with:  (3.20) is mainly controlled by m θ i . We can therefore derive (3.20) just with respect to m θ i and obtain: (3.23) where we used (3.13) approximated as dm θ i m θ i dN . This implies for both θ 1 and θ 2 : (3.24) Thus we realise that also the distribution of the ALP DM abundance features a logarithmic behaviour. As already explained, this result should apply also to the distribution of the DM abundance of each ultra-light ALP of the model with arbitrarily large h 1,1 .

Anisotropic limit
In the anisotropic limit with the SM on the fibre divisor, the value of the QCD axion decay constant is fixed around the GUT scale once we focus just on vacua which match the SM coupling. Hence this would represent a typical stringy case which has been considered as 'anthropic' since for f a ∼ M GU T ∼ 10 16 GeV (3.19) would reproduce the correct DM abundance only for θ in ∼ 0.001π. The DM abundance associated to the ultra-light ALP θ 2 would instead be distributed as: where we used dm θ 2 m θ 2 dN from (3.15). Similarly to the isotropic case, we find again a logarithmic distribution with however the difference, as already pointed out, that in the anisotropic case all expressions, (3.20) included, are valid only for f θ 2 10 12 GeV (while we have seen that any value of m θ 2 is allowed).

Axion couplings to gauge bosons
Let us now study the distribution of the couplings between axions and gauge fields. Following the analysis of the previous sections, we first evaluate the couplings at the minimum of the scalar potential and we then determine their distribution in the flux landscape.

Isotropic limit and model with arbitrary h 1,1
We start with the couplings in the isotropic case which are summarised in (2.35). The axion couplings to visible and hidden gauge bosons feature a similar behaviour also in the model with arbitrarily many Kähler moduli, as can be seen from (2.53). Hence the results which we will obtain for the isotropic case can be directly extended to the more general case with arbitrarily large h 1,1 and SM built with a stack of D7-branes wrapped around a local del Pezzo divisor.
Interestingly, each ultra-light ALP couples in practice just to the corresponding hidden gauge fields with a coupling that is fixed at Planckian strength without a distribution in the landscape. This is a typical stringy behaviour, as expected for the imaginary part of a standard bulk modulus. On the other hand, the coupling between the QCD axion a and SM gauge fields γ is controlled by the string scale M s ∼ M p / √ V since it is inversely proportional to f a : Thus the distribution of g aγγ takes the form: where we used (3.3) approximated as df a f a dN . Notice the mild logarithmic preference for smaller couplings. The coupling of the QCD axion to hidden gauge bosons γ h living on stacks of D7-branes wrapped around the bulk divisors D 1 and D 2 , is instead weaker than Planckian since these two divisors do not intersect with the del Pezzo 4-cycle D 3 : Hence, using df a f a dN g aγ h γ h dN , the distribution of the coupling between the QCD axion and hidden gauge bosons scales as: Contrary to the coupling to visible gauge fields, in this case the flux landscape features a logarithmic preference for larger couplings.
where the CY volume takes the same form as in (2.3) and the fibre modulus τ 1 plays the role of the inflaton. The inflationary potential is generated by perturbative corrections to the Kähler potential and the CY volume is fixed around V 10 3 -10 4 by the need to reproduce the observed amplitude of the density perturbations generated by inflaton fluctuations during inflation. In order to have an efficient production of SM degrees of freedom at reheating, the SM D7-stack has to wrap the fibre divisor D 1 . Hence a viable realisation of Fibre Inflation models requires to focus on the anisotropic case.
The inflaton τ 1 is the lightest Kähler modulus and its perturbative decay after the end of inflation produces SM particles together with the QCD axion θ 1 and the ultra-light ALP θ 2 which are both relativistic and yield a g s -dependent contribution to the effective number of relativistic species N eff [54]. One can thus exploit the known distribution of g s to derive the distribution of extra dark radiation in the flux landscape of Fibre Inflation models. The amount of extra dark radiation is parameterised by ∆N eff which is determined by the ratio of the inflaton branching ratio into hidden and visible degrees of freedom [54]: where the parameter γ controls the coupling of the inflaton to visible sector gauge bosons and depends on the string coupling: where we have used (2.36) with α = 1 in the the volume form (2.3). Let us stress that in Fibre Inflation models the CY volume is fixed around V 10 3 by the need to reproduce the observed amplitude of the density perturbations generated by inflaton fluctuations during inflation. Hence in (3.34) V should be considered as constant. When varying g s , this can be achieved by an appropriate choice of W 0 (see (2.8)). Moreover g s should be varied by keeping the SM coupling fixed at its phenomenological value. Given that α SM reads: where h(F) ≥ 0 is a non-negative function of the intersection numbers and the gauge flux F on the SM D7-brane stack, this implies that any variation of γ (by varying g s ) should be compensated by a suitable change of h(F) by considering a different choice of F (if this is allowed by the discreteness of the gauge flux quanta and by tadpole cancellation). Notice that h(F) vanishes for F = 0, implying from (3.35) γ = 1 and ∆N eff fixed at ∆N eff 0.6 [54]. However h(F) > 0 for F = 0, and so in this case ∆N eff features a distribution in the flux landscape due to its dependence on g s . We can estimate the regime of validity of this distribution by setting in (3.34) α −1 SM = 25, V = 5 · 10 3 and g s 0.25 to trust perturbation theory, which gives ∆N eff 0.17. We can also obtain an upper bound on ∆N eff by requiring τ 1 ≥ α −1 SM from (3.35) since h(F) > 0. This gives γ ≥ 1 from (3.34), and so ∆N eff 0.6.
Varying now (3.33) with respect to γ using (3.34) and dg s dN , we obtain: Interestingly we find that the flux landscape of Fibre Inflation models features more vacua around ∆N eff 0.17 which helps to satisfy current bounds on extra relativistic species. We stress again that this distribution is valid only for values of ∆N eff corresponding to values of g s which are compatible with a choice of h(F) that keeps α SM constant.

Discussion and conclusions
In this paper we studied the statistics of axion physics in the type IIB flux landscape focusing on the model-independent case of closed string axions coming from the dimensional reduction of C 4 . We argued that a proper understanding of moduli stabilisation is crucial in order to derive the main features of the low-energy phenomenology of stringy axions.
In KKLT-like scenarios all axions are as heavy as the corresponding saxions due to non-perturbative stabilisation. If the saxion masses are larger than O(50) TeV in order to avoid cosmological problems, each axion is thus too heavy to behave as the QCD axion or as a very light ALP for fuzzy DM. On the contrary, moduli stabilisation schemes which rely on perturbative corrections are characterised by axion masses which are exponentially suppressed with respect to saxion masses. This singles out LVS models as the best case scenarios for analysing axion physics since they also yield an exponentially large CY volume which allows to keep the EFT under control even for a large number of Kähler moduli.
Hence we focused on an LVS model with h 1,1 = 4 which is simple enough to perform moduli stabilisation in full detail but, at the same time, rich enough to show all the main features of axion physics which we consider to be valid in general for models with more Kähler moduli. We considered two regimes: (i) the isotropic limit with the SM on D7branes wrapping a local blow-up cycle, and (ii) the anisotropic limit where the SM lives on a D7-stack wrapped around a bulk divisor. In both cases all phenomenologically interesting quantities, like axion decay constants, axion masses, contributions to the DM abundance and axion-gauge bosons couplings, feature a logarithmic distribution in the flux landscape. In the isotropic case however, the request to reproduce the correct SM gauge coupling selects a subset of the underlying parameter space where some distributions turn into a mild power-law behaviour.
Regarding the QCD axion, in the isotropic case it comes from the reduction of C 4 on a blow-up mode, whereas in the anisotropic case it is associated to a bulk cycle. In the last case its decay constant is fixed around the GUT scale by the need to match α SM . On the other hand, in the first case f a is distributed logarithmically with just a mild preference for GUT-scale values in comparison with cases where f a is around intermediate scales. We consider this case to be more generic in the string landscape since realisations of the QCD axion from bulk cycles require an anisotropic moduli fixing which would require a good amount of tuning for relatively large values of V, while the case of a blow-up QCD axion can work with either isotropic or anisotropic models. We therefore conclude that what has been so far claimed to be the typical stringy situation with a GUT-scale QCD axion decay constant and a tuned initial misalignment angle to avoid DM overproduction, could be not so predominant in the flux landscape with respect to more natural cases where f a ∼ O(10 11 ) GeV and θ in ∼ O(π).
On top of the QCD axion, the isotropic and anisotropic scenarios feature either 1 or 2 ultra-light ALPs. In agreement with previous studies [25,[27][28][29], we argued that the presence of several ultra-light ALPs is a general characteristic of 4D string models where the EFT is under control, as we have shown explicitly in a model with arbitrary h 1,1 where full moduli stabilisation can be achieved by exploiting higher derivative α corrections following [26]. Interestingly, we found that the decay constants, the mass spectrum and the contribution to the DM abundance of all these ultra-light ALPs are also logarithmically distributed in the type IIB flux landscape.
In our recent paper [5] we found that the the number of flux vacua is also a logarithmic function of the gravitino mass and the supersymmetry breaking scale. Moreover in App. C we showed that other quantities relevant for phenomenology, as the moduli masses and the reheating temperature from moduli decay, share the same statistical properties. We are therefore tempted to argue that most, if not all, of the low-energy properties of the string theory landscape seem to obey a logarithmic distribution once moduli stabilisation is properly taken into account. Apart from the particular case of extra dark radiation in Fibre Inflation models, the only exception which we have encountered so far seems to be the supersymmetry breaking scale in KKLT scenarios which might be characterised by a power-law distribution. However its statistical significance is still unclear since this result relies on the assumption that W 0 is uniformly distributed [10] also in the exponentially small regime where however it is very hard to built explicit examples. The only ones which have been constructed so far feature W 0 = 0 and a flat direction at perturbative level [13][14][15]. The flat direction is lifted by non-perturbative physics which generates dynamically an exponentially small W 0 ∼ e −1/gs . Exploiting the uniform distribution of the string coupling, this relation would again produce a logarithmic distribution of the gravitino mass.
It is worth stressing that these distributions follow from moduli stabilisation which applies only to corners of the string landscape where the EFT is under control thanks to supersymmetry and weak couplings. In order to judge their genericity one would have therefore to be able to control the EFT beyond the regime of validity of these approximations. Despite the difficulty to achieve this goal, scaling arguments and approximate symmetries inherited from the 10D theory [6] could be used as a powerful guideline to shed light on larger portions of the string landscape. This top-down analysis of the statistical properties of quantities relevant for phenomenology is crucial to provide more theoretical guidance to recent bottom-up approaches to understand naturalness and string theory predictions for several observables [55][56][57][58][59][60].
We finally comment on the fact that the relative flatness of logarithmic distributions in the string landscape might be seen at first sight as an indication of a difficulty to make sharp predictions from string theory. However a key-feature of string theory is the correlation between different low-energy phenomenological quantities due to the underlying UV framework. It is this interplay which should be used to sharpen the predictions of the string landscape. As an example, we mention the fact that in LVS models an intermediate scale QCD axion decay constant would correlate with TeV-scale soft terms and a volume mode mass around 1 MeV. Thus, in the absence of a mechanism to avoid cosmological problems associated to the presence of such a light modulus [61], a natural QCD axion DM situation with f a ∼ O(10 11 ) GeV and θ in ∼ O(π) would not be viable even if the number of vacua with these features is only logarithmically suppressed with respect to the number of vacua with a GUT-scale decay constant. We leave the important study of the UV correlation between different particle physics and cosmological observables with logarithmic distributions for future work.

A Canonical normalisation
In this appendix we shall perform the canonical normalisation of the axion fields.

A.1 A single axion
Let us start with the simple case with a single closed string modulus T = τ + iθ where it is easy to identify the correct definition of the axion decay constant and periodicity. We start with the following Lagrangian: where b is a non-Abelian index and the gauge kinetic function is given by f = T /(2π).
Expressing L in terms of the canonically normalised axion a = 2K TT θ and Yang-Mills field strength G b µν = Re(f ) F µν b , we end up with: where we used the fact that τ = α −1 b . This expression suggests the definition of the axion decay constant f a as (inserting the appropriate power of M p ): since L would simplify to the standard expression:   Table 2: Benchmark points which match the observed ALP DM abundance for the anisotropic case setting A 2 = A 4 = 1, θ in,2 = π, α = 1/6 and ξ = 0.46. All benchmark points satisfy the phenomenological constraint τ 1 = α −1 SM = 25.

C Other distributions relevant for phenomenology
In this appendix we shall show that other phenomenologically interesting quantities feature also a logarithmic distribution in the type IIB flux landscape.

C.1 Moduli masses
Let us investigate the distribution of moduli masses in the flux landscape. For all cases, the isotropic and anisotropic cases with h 1,1 = 4 and the model with arbitrarily large h 1,1 , the mass of each Kähler modulus scales with the CY volume as: with p i > 0 ∀i = 1, ..., h 1,1 . (C.1) Following the same logic as in Sec. 3, we find again a logarithmic distribution for each m τ i since these masses are controlled by the exponentially large volume V: N (m τ i ) ∼ ln m τ i M p , ∀i = 1, ..., h 1,1 .

(C.2)
For the anisotropic case one has just to make sure that the bound g s 0.01 (coming from the ability to tune W 0 to satisfy (2.41)) does not set a lower bound on m τ i for the regime of validity of the distribution (C.2). However this bound is negligible since combining (2.8) with (2.41) one would find m 3/2 10 −45 M p 10 −16 eV for g s 0.01.

C.2 Reheating temperature
Using the moduli masses we can study the distribution of the reheating temperature coming from moduli decay [54,62]. The reheating temperature due to the perturbative decay of the i-th Kähler modulus is given by [54]: T rh,i = 40c vis c tot π 2 g * (T rh ) where c vis and c vis control the strength the interaction of the modulus τ i with the visible and the hidden sector respectively, c tot = c vis + c hid , and the decay rate Γ τ i looks like: (C.4) Thus using (C.2) we obtain again a logarithmic distribution for all cases: