ALP Anarchy

String theory models generically predict the existence of multiple axion-like particle (ALP) fields, yet the majority of both theoretical and experimental works have assumed only one ALP. In this paper, we discuss the phenomenology of systems with multiple ALPs that can undergo oscillations akin to neutrino oscillations. Motivated by this effect, we extend the 'anarchy' framework, which has been used to predict neutrino oscillation parameters, to generate the parameters of many ALP systems. We explore the phenomenology of these ALP anarchy models in some of the leading ALP search strategies, including the CERN Axion Solar Telescope, magnetic white dwarfs and the gamma-ray spectra of distant blazars. We include both the ALP-photon and the ALP-electron coupling. We find that ALP anarchy models predict drastically different results than single ALP models.


Introduction
Axion-like particles (ALPs) have emerged as a leading candidate for physics Beyond the Standard Model.There are three key physics motivations for the existence of ALPs.First, they behave as cold dark matter [1][2][3][4][5].Second, the QCD axion offers a promising solution to the strong CP problem [6,7].Finally, string theory suggests the existence of a large number of ALPs [8][9][10][11][12][13].Unlike the QCD axion, ALPs need not couple to gluons; therefore, their mass and Standard Model couplings are not necessarily related.ALPs are pseudo-scalar fields and are singlets under the Standard Model gauge group.This determines the form of their interactions with the Standard Model fields.For example, the Lagrangian of a single ALP ϕ interacting with photons and electrons is given by where g γ is the coupling of ϕ to the electromagnetic field strength tensor, F µν , g e is the dimensionless ALP-electron coupling, m e is the electron mass and ψ denotes the electron field.Models containing many ALPs may differ qualitatively in their phenomenology from those containing a single axion or ALP.Recent works have explored many ALP scenarios in a range of contexts [14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31].The goal of this paper is to explore the interaction of many ALP models with the Standard Model.As pointed out in Ref. [26], this can be understood by considering the misalignment between the ALPs' mass basis and the basis in which only one linear combination of ALPs couples to electromagnetism.As shown below, this results in oscillations between the electromagnetically coupled ALP and a number of orthogonal hidden ALP states where the physics at work is analogous to the physics of neutrino oscillations.This effect has also been studied in the context of Kaluza-Klein axions [32].This paper will expand on the results of Ref. [26] by considering several new aspects of these ALP oscillations.In particular, we will introduce ALP anarchy models, similar to the anarchy approach that has been used in neutrino physics to explain the structure of the leptonic mixing matrix.We will show that these many ALP models display dramatically different phenomenology to single ALP models with the same effective couplings.The framework developed here also applies to any system with multiple axion or ALP fields.
This paper is structured as follows.In Section 2 we will introduce many ALP models and the oscillation effect.In Section 3 we will introduce ALP anarchy models.In Section 4, Section 5 and Section 6 we will calculate the phenomenology of ALP anarchy models in the Cern Axion Solar Telescope, in magnetic white dwarfs and in the gamma-ray spectra of distant blazars respectively.Finally, in Section 7, we will discuss our results further and conclude.

ALP oscillations
We will consider a model containing N ALPs, each coupling to both photons and electrons: where ϕ i is an ALP field with mass m i , g γ i is the coupling of ϕ i to the electromagnetic field strength tensor, F µν , g e i is the dimensionless ALP-electron coupling, m e is the electron mass and ψ denotes the electron field.In this paper, we will consider only ALP interactions with photons and electrons, but similar considerations would apply to other couplings.We note that Eq. (2.1) is in the mass basis.In this work, we will consider string-motivated ultralight ALPs, and hence, we will assume that the QCD axion (which may also emerge from string theory) is too heavy to contribute to the scenarios considered here.Furthermore, string ALPs with masses heavier than the energy scales considered below will also not contribute.
It will be convenient to rotate to the electromagnetic basis, in which only a single ALP (the 'electromagnetic ALP') couples to photons: where g γ = i g γ i 2 is the total effective ALP-photon coupling, we have chosen ϕ 1 = i g γ i ϕ i /g γ as the electromagnetic ALP and ϕ i and g e i have been appropriately redefined.This basis has the advantage that only ϕ 1 interacts directly with the photon.The other ALPs, ϕ {2...N } , are 'hidden' with respect to the electromagnetic interaction.
For example, any ALP production via the electromagnetic interaction in cases where the mass is irrelevant will produce the electromagnetic ALP state ϕ 1 .The production rate can, in this case, be calculated simply by considering a single ALP with coupling g γ to photons.However, this does not necessarily mean we can treat the system as though there is only one ALP with coupling to photons g γ .As we assume that the ALP mass states differ and that the mass and electromagnetic basis are misaligned, ϕ 1 may mix with the hidden states, ϕ {2...N } .This is an analogous effect to neutrino oscillations resulting from misalignment between the neutrinos' interaction and mass eigenbases.
Similarly, we can define the electronic basis in which only one ALP (the 'electronic ALP') couples to electrons: where g e = i g e i 2 and we have now chosen ϕ 1 = i g e i ϕ i /g e as the electronic ALP and ϕ i and g γ i have been appropriately redefined.As with the electromagnetic ALP, this basis has the advantage that only ϕ 1 interacts directly with the electron while the other ALPs ϕ {2...N } are 'hidden' with respect to the electron interaction.Note that the electronic and electromagnetic ALP states are generally neither orthogonal nor colinear.Hence, in scenarios where both ALP-photon and ALP-electron interactions are relevant, we must potentially consider three different bases -the mass, the electromagnetic, and the electronic.If other interactions between the ALPs and the Standard Model are relevant, these may introduce further relevant bases.
As shown below, ALP oscillations are phenomenologically significant for many observations, particularly when we hope to detect an ALP that has propagated a large distance.However, ALP search strategies that rely only to the disappearance of Standard Model particles or energy into ALP degrees of freedom, such as stellar cooling bounds on ALPs [33][34][35], are not significantly effected by ALP oscillations.Therefore, comparison between ALP searches is substantially more complicated in many ALP systems than if we assume only a single axion or ALP.

Anarchy models
String theory provides a framework to understand ALP properties, including their mixing matrices that parameterise the misalignment between the interaction and mass bases.Calculating these mixing matrices is a highly non-trivial task [36].Nonetheless, the ALP photon coupling has been modelled in a range of string axiverse scenarios [13,37].Such modelling of the electronic ALP properties from string theory has not been undertaken.
This current lack of knowledge of the mixing matrices, from a first principle string theory calculation, presents a challenge when motivating the choice of ALP mixing parameters and couplings.In this work, we circumvent this issue and remain agnostic to the ALPs' particular ultraviolet physics by considering a large set of randomly sampled mixing matrices.This framework, also known as anarchy, has been applied in neutrino physics and refers to the postulate that the neutrino mass matrix has no particular structure but that its elements are randomly chosen O(1) parameters [38][39][40][41][42][43][44][45][46][47][48][49].Randomness in O(1) coupling constants is expected in sufficiently complicated models or with many fields mixing with each other.While it remains unclear if string theory predicts an anarchy-like mixing pattern, in this work, we use anarchy to explore the general properties of multi-ALP phenomenology and their relation to the number of ALP mass eigenstates.
To implement this anarchical approach and determine the mixing matrices between the hidden and visible ALPs, we assume that the non-diagonal ALP mass matrices, M γ and M e , are real and can, therefore, be related to the mass basis states as follows: where we assume U α ∈ SO(N ) and D = diag (m 1 , m 2 , . . ., m N ) is a real diagonal matrix with m i denoting the mass of ALP field ϕ i .Following our assumption that the electromagnetic and electronic bases are misaligned, a given model will be described by a set of two rotations U α -one for each interaction basis.We sample SO(N ) such that elements are uniformly distributed over the group manifold to generate these mixing matrices.This can be achieved using the Haar measure, which describes the density of elements in a Lie group.In spherical coordinates, the Haar measure for SO(N ) is: where there are N (N − 1)/2 mixing angles, {θ ij } 1≤i≤j≤N −1 with θ 1,j ∈ [0, 2π] and θ i,j ∈ [0, π] for i, j > 1. Sampling uniformly in dV yields the desired distribution of the angles that parameterise the mixing matrices U γ and U e and we outline our numerical procedure for this task in Appendix A. Mixing matrices, although providing a simple means to understand a given parametrisation, contain a large degree of redundancy.In practice, and in what follows, we are more interested in the relationship between the couplings in the mass basis -{g γ i } and {g e i }.Given, for example, the EM coupling in the EM-ALP basis, the mass basis couplings are given by: An analogous statement is also true for the electron-ALP couplings.From Eq. (3.3), it is clear that the top row of U γ parameterises the mixing of the electromagnetic ALP with the hidden states, which will have observable phenomenological consequences.The remaining N − 1 rows determine the mixing between the hidden states, which is not observationally relevant.We note that the mixing matrices U parameterise only the misalignment between the electromagnetic, electronic and mass bases and not the magnitude of the total effective couplings g γ and g e .This allows us to distinguish between the effects of basis misalignment and the effects of simply varying the total coupling, which is also present in the single ALP case.In the following sections, we will explore the phenomenology of ALP anarchy models.In particular, we will compare anarchy models with different numbers of ALPs to single ALP models with the same total effective couplings.We note here that, in a string axiverse scenario, the fundamental parameters are the couplings g γ i and g e i of the mass eigenstates and the effective couplings of the electromagnetic and electronic ALP states emerge from these.In this case, the possible values of the effective couplings g γ and g e are determined by the string model.

The Cern Axion Solar Telescope
In this section, we will outline the production of electromagnetic and electronic solar ALPs in the Sun and their detection by CAST.We will discuss how oscillations, within the framework of an anarchical mixing pattern, influence the number of electromagnetic solar ALP states that reach the CAST detector.Finally, we will explain how we reinterpret the CAST experimental constraint on the single ALP parameter space for our multi-ALP scenario.Of the potential astrophysical sources, the Sun is one the most accessible to the search for ALPs as the internal dynamics of the Sun are well understood and can be accurately modelled by a weakly coupled plasma.Within this setting, it is possible to find analytical forms for the emitted ALP flux, and these can be used to set competitive bounds on ALP couplings for mass ranges relevant to solar processes.ALPs are produced predominantly through three mechanisms: axio-bremsstrahlung (the Primakoff process) [50,51], axio-recombination and axio-deexcitation [52][53][54], and Compton scattering [55][56][57], see Ref. [58] for a comprehensive overview.Of these processes, shown in Fig. 1, only the Primakoff process depends on the ALP-photon coupling with the others arising from the ALP-electron interaction.The CERN Axion Solar Telescope (CAST) is a helioscope designed to use the Primakoff effect to scatter axions or ALPs produced by the Sun into photons.CAST consisted of a 9.2 m long evacuated cylinder with a sustained magnetic field oriented transverse to the ALP propagation direction.A bound on g γ and g e can be placed based on the lack of an ALP signal.In this work, we consider how the multi-ALP scenario will alter this bound.This modification occurs as the electromagnetic and electronic ALP states produced in the Sun may oscillate into hidden ALP states which CAST cannot detect.

Solar ALP emission and detection
If ALPs couple to the electromagnetic field strength tensor or an electronic current, as in Eq. (2.1), they can be produced in the Sun.The relevant ALP production processes are shown in Fig. 1, and the dominant ALP production mechanisms are Bremmstrahlung (B), Compton Scattering (C) and the Primakoff process (P ).In the well-studied single ALP case, with no oscillation into hidden ALPs, the fluxes at Earth, in units of m generated by these mechanisms are given by [59]: where the ALP energy, ω, is in keV.A more detailed study of the solar ALP flux can be found in [60].Integrating over the energy range (0.8 − 6.8 keV) of the CAST analysis [61], the total fluxes can be found for Φ B , Φ C and Φ P : The CAST experimental constraints assume a single ALP that couples to electromagnetism and electrons.We will now consider the CAST signal from models with multiple ALP mass eigenstates.

Oscillation: The two ALP case
We will expand upon the results of [26] by discussing the basics of ALP oscillations relevant to CAST in a simplified two ALPs scenario where ϕ e and ϕ γ are linear combinations of only two massive ALP states, ϕ 1 and ϕ 2 , which have masses m 1 and m 2 respectively.The Lagrangian in the mass basis is We can rotate to a basis where a single ALP field, the electromagnetic ALP, couples to electromagnetism while the other hidden electromagnetic ALP field does not.These are given by respectively.Analogously, the electronic ALP and the hidden electronic ALP are the following orthogonal combinations: respectively.Note that the electromagnetic and electron ALPs ϕ γ and ϕ e are generally neither colinear nor orthogonal.Therefore, ϕ γ h will in general have some non-zero coupling to electrons and ϕ e h will in general have some non-zero coupling to photons.The following unitary rotation matrices relate the mass basis and the electromagnetic and electron bases: with where α = γ, e.
We will assume that all couplings are real, and their values are determined using the anarchical approach outlined in Section 3. We will also assume that m 1 and m 2 are much less than other relevant energy scales.In particular, we will assume m 1 , m 2 < 10 −2 eV, corresponding to CAST bounds for evacuated magnet bores.In this mass range, the ALP masses are also much lower than their production energy in the Sun and can be treated in the relativistic limit.This means their effects on ALP production may be neglected, so axio-recombination, axio-de-excitation and Compton scattering produce the state ϕ e while Primakoff production produces the state ϕ γ .Furthermore, in this mass range, the CAST sensitivity is independent of mass, and therefore, CAST will detect the state ϕ γ as a single signal.
As CAST aims to detect electromagnetic ALPs, ϕ γ , we are interested in the probability that a solar electronic or electromagnetic ALP oscillates to an electromagnetic ALP after travelling a distance L. We can calculate these probabilities using the fact that the mass eigenstates propagate as |ϕ i (L)⟩ = e −i m 2 i L 2ω |ϕ i (0)⟩, where L is the distance travelled.The former is given by ) where ∆m 2 = m 2 2 −m 2 1 is the mass squared splitting between the ALP mass states.Rewriting Eq. (4.12) in terms of couplings, g α , yields In addition to the electronic ALPs produced in the Sun oscillating to electromagnetic ALPs when they reach CAST, we must consider the survival probability of the solar electromagnetic ALP which is the usual two-state survival probability familiar from neutrino physics: Several simplifications can be made for solar ALP oscillations: first, the matter potential induced by the solar electron background is negligibly small and, therefore, does not affect the electron ALP propagation through the Sun.This can be estimated from the fact that the potential experienced by electron neutrinos from electrons in the Sun is V ≈ G F N e ∼ 10 −12 eV where G F is Fermi's constant and N e is the number density of electrons in the Sun's core.In contrast, the potential experienced by the electronic ALP (with coupling g e 2me = 10 N e ∼ 10 −29 eV.Likewise, the potential induced by the Sun's magnetic field is negligible, so vacuum oscillation between the ALP states will be applied.Second, for Sun-Earth distances, keV solar ALP energies with ∆m 2 > 10 −12 eV 2 , which we assume, the oscillation probability of Eq. (4.12) and Eq.(4.13) averages when integrated over the CAST energy range yielding:

Oscillation: The many ALP case
We now turn to the case where many ALP mass eigenstates couple to electrons and photons: such that the electromagnetic and electronic ALPs produced in the Sun are linear combinations of the mass states: where we again assume that all ALP masses considered are m i < 10 −2 eV.Any mass eigenstates with m i > 10 −2 eV would not contribute to the signal considered here, as they would not produce a signal in CAST with an evacuated bore.Again, we will assume that ∆m 2 > 10 −12 eV 2 so that the oscillation probabilities average when integrated over the CAST energy range.Under these conditions, it can be shown that the electromagnetic ALP survival probability is where N is the number of ALP mass eigenstates, is the coupling of the electromagnetic ALP to photons, and VAR({g γ i 2 }) is the variance of the ALP-photon couplings squared in the mass basis which is given by Likewise, the probability that an electronic ALP oscillates to an electromagnetic ALP is where g e = g e i 2 is the coupling of the electronic ALP to electrons and is the covariance of the ALP-photon and ALP-electron couplings squared and is

Reintepretation of CAST results
The non-observation of an excess number of photons allows CAST to place constraints on the (g γ , g e ) parameter space, which bounds the solar ALP fluxes on Earth.In the multi-ALP case, the fluxes of the electromagnetic and electronic ALP on Earth, namely the oscillated fluxes, are, respectively: where Φ γ = Φ P and Φ e = Φ B + Φ C and P a→b is the probability of species a oscillating into species b during propagation as derived in Section 4.3.Since CAST makes use of a strong magnetic field to convert ALPs, the relevant flux to consider for the recasting from the single to the multi-ALP scenario is Φ osc γ .CAST has placed bounds on (g γ ) 2 as a function of ALP mass and on g γ as a function of g e , for small ALP mass m a , for a single ALP.Here we determine an upper bound on g γ as a function of g e for multi-ALP scenarios with m a ≤ 10 −2 eV and ∆m 2 > 10 −12 eV 2 for all relevant ALP mass eigenstates; for which the bound becomes mass independent [61].We do this by considering the maximum flux compatible with the non-detection of ALPs by CAST.Given the original single ALP bound on the ALP-electromagnetic coupling, which we denote as g γ N =1 (g e ), (shown by the black line in Fig. 2), the maximum flux of electromagnetic ALP on Earth is bounded to be: where Φ N =1 = Φ P + Φ B + Φ C is the flux at the detector in the single ALP case.With Φ osc γ (g e , g γ ) in hand, we now perform a grid scan over the coupling parameter space, (g e , g γ ), to determine the new bound.For a given point in the ALP anarchy parameter space to be allowed by existing CAST data, the total EM-ALP flux at the detector given in Eq. (4.22) must be less than Φ Max .
To determine the bound on the total effective couplings in the ALP anarchy scenario, it is therefore necessary to compute P γ→γ and P e→γ .From Section 4.3, these are given by Eq. (4.18) and Eq.(4.20), respectively, and can be written in terms of the relationship between the individual ALP couplings (g e i and g γ i ).It is this relationship that encodes the mixing between the various ALP states.We consider 10 4 realisations of {g γ i } and {g e i } for each overall coupling pair (g γ , g e ).For each point in (g γ , g e ) space, we determine the proportion of viable realisations such that Φ osc γ < Φ Max (g e , g γ ).These results are shown in Fig. 2; for a detailed outline of this numerical procedure, see Appendix B. The black line indicates the original N = 1 bound (ϕ e ≡ ϕ γ ) where in the mass independent region log 10 g γ [GeV −1 ] ≲ −10.The boundary between regions where 0 and 100% of realisations are viable can be interpreted as the bound on g γ as a function of g e in the ALP anarchy scenario.From Fig. 2, it can be seen that increasing the number of ALPs decreases the competitiveness of the bound.The left plot of Fig. 2 shows the bound with N = 2 ALP fields, and we observe that the effect of an additional state is marginal; however, for N = 30, we observe that the bound on g γ is relaxed by almost half an order of magnitude with log 10 (g γ [GeV −1 ]) ≲ −9.6.As the number of hidden states increases, the oscillated flux of electromagnetic ALPs on Earth decreases since they can oscillate into the hidden ALP states that are not detectable by CAST.Hence, the effective coupling g γ can increase to compensate for this decrease in the detectable flux.To quantify this relationship, we consider ALP multiplicities N ∈ [2,30].For each N in this set, we determine the value of g γ , g γ 50 (N ) for which 50% of mixing realisations satisfying Φ osc γ < Φ Max (g e , g γ ) at a point in the horizontal region of (Fig. 2).We fix g e = 10 −15 .For this low ALP-electron coupling, the production processes for ϕ e , namely bremsstrahlung and Compton scattering, are ineffective, and the production of ϕ γ in the Sun dominates.g γ 50 (N ) can be interpreted as an approximate bound on g γ in the ALP anarchy scenario.This result is shown in Fig. 3, which we fit the following function to: where we find m = 0.25 and c = −23.This fitting function is also shown in Fig. 3.We can understand the dependence of the bound on Ng γ 50 (N ) ∝ N 1/4 .We first note that when electromagnetic ALP production dominates, CAST is sensitive to g γ 4 , as the ALP must be produced and detected via this coupling to photons.In the many ALP scenario, the bound placed on g γ 4 is weakened in proportion to the electromagnetic ALP survival probability P γ→γ given in Eq. (4.18).For large N , the variance term becomes negligible and we have P γ→γ ∼ 1 N .We therefore obtain g γ 50 (N ) ∝ N 1/4 , as found numerically.We have seen that if the ALP-photon and ALP-electron interactions are an effect from multiple ALP mass eigenstates, the CAST bounds on the total effective couplings may be somewhat reduced.We have calculated this reduction in the ALP anarchy scenario for mass differences ∆m 2 > 10 −12 eV 2 .As shown in [26], for ∆m 2 < 10 −14 eV 2 , there is no significant oscillation into hidden states.For intermediate mass differences, there is a non-trivial oscillation structure depending on ω.

Magnetic white dwarfs
Having considered the constraints from the CAST experiment, we now examine how the limits on the single ALP electromagnetic and electronic coupling from observations of magnetic white dwarfs (MWDs) can be applied to our multi-ALP scenario.MWDs produce ω ∼ keV energy ALPs via axio-bremsstrahlung, which can convert to X-rays in the magnetosphere surrounding the MWD.Subsequently, searches for observable X-ray signals provide one of the most stringent constraints on the ALP parameter space, see e.g.Ref. [63].More specifically, in the single ALP scenario, where N = 1 and ϕ e ≡ ϕ γ , the (g e , g γ ) parameter space is constrained from the non-observation of astrophysical X-ray emission from the MWD RE J0317-853 by the Suzaku telescope [64].The flux of ALP-induced X-ray photons on Earth, in the low ALP mass regime, is approximately proportional to Φ X-ray ∝ (g e g γ ) 2 as the ALP luminosity is proportional to g e2 while the probability the ALP transitions to X-ray photons is proportional to g γ 2 .The non-observation of excess X-rays provides an upper bound on Φ X-ray and hence a corresponding bound on (g e , g γ ).
In the multi-ALP case, since the radius of the MWD is relatively small (less than a percent of the Sun's radius [65]), the oscillations between the electronic ALP and hidden states do not have time to develop, assuming ∆m 2 ≲ 4R = 10 −10 eV 2 , where R is the MWD's radius.Moreover, we assume that all ALP masses considered are m i < 10 −2 eV.To perform a simple recast of the single to the multi-ALP case, we compute the conversion probability (P e→γ ) and scale the N = 1 bound on g e g γ , denoted as b N =1 , as follows: where b N >1 is the new bound on g e g γ assuming the existence of N ALP states.A decrease in the conversion probability will allow for a greater possible coupling strength.In the MWD setting, in the case where oscillations do not have time to develop, P e→γ is given by the scalar product of ϕ e and ϕ γ : Note that this is a different limit to that considered in Section 4 as the distance propagated is much lower.The resulting bound as a function of number of axions is shown in Fig. 4 for values: L = 4 × 10 −3 Solar Radii [65], E = 5 keV and masses distributed logarithmically in [10 −9 , 10 −6 ] eV.The conversion probability, accounting for the presence of hidden ALPs (in the simple two-ALP case, see Eq. (4.12)), is given by which describes the projection of the electronic ALP state onto the EM ALP state.In the limit θ e = θ γ , naturally, the probability P e→γ = 1 because the electron and EM ALP states are completely aligned, and the N = 1 case is recovered.If these angles are very different, then the probability P e→γ can be significantly smaller and make the bound proportionally weaker.In Fig. 4, we show how the limit on |g γ g e |, in the low ALP mass regime, from Ref. [63] changes as the number of hidden states is increased.We observe that as the number of hidden states is increased, the total allowed coupling strength of the ALP to EM and electrons increases.This occurs because increasing the number of hidden states leads to an increase in the typical orthogonality between the electromagnetic and electronic ALP states, decreasing the chance that an ALP produced in an electron interaction will convert into a photon.We find that the bound decreases almost three orders of magnitude as the number of ALP mass states increases from 2 to 18.

Very high energy blazars
In this final section, we examine how our multi-ALP scenario can be constrained by the observation of very high-energy gamma rays.We begin by outlining the simulation of the propagation of these high-energy photons emitted by blazars toward Earth, considering the possibility of photon mixing with multiple ALP states.We detail the density matrix equations and the model of the magnetic fields used in the simulation.Additionally, we discuss the selection of blazars and how these sources can be utilised to constrain the electromagnetic coupling as a function of the number of ALPs in our multi-ALP scenario.Blazars produce a large flux of Very High Energy (VHE) photons.These TeV scale photons can scatter off the isotropic extragalactic background light (EBL) as they propagate to Earth, producing positron-electron pairs.The probability with which this scattering occurs increases with energy.As such, we expect significant attenuation of high-energy photons travelling through intergalactic space.Telescope observations suggest that the Universe may be more transparent than expected to VHE photons [66,67], although the evidence for this effect is not conclusive [68,69].Several phenomenological studies have introduced ALPs to explain this discrepancy [70][71][72][73].VHE photons can oscillate into ALPs in the magnetic field of the blazar or the intergalactic medium and thus travel unimpeded through the Universe.These oscillations may conspire for appropriate masses and couplings, such that the ALPs reconvert into photons in the Milky Way, allowing for their detection on Earth.In this case, the measured flux of VHE photons on Earth will be amplified, accounting for the observations.In this section, we will explore this effect in the context of our multi-ALP model.Following and extending the analysis of Ref. [73], we consider an arbitrarily mass-mixed set of ALP states.We note that the ALP-electron coupling does not play a significant role in this process, and therefore will not be considered in this section.

Simulation
In the following subsections, we describe the approach by which we simulate the propagation of VHE photons from their blazar source to Earth.We consider both photon-EBL scattering and ALP-photon mixing.As the ALP mass is relevant for this effect, we work in the mass basis where each ALP couples separately to the photon.Photon-EBL scattering is dissipative (VHE photons can scatter off EBL photons, creating positron-electron pairs) and introduces non-unitary into the evolution.To account for this effect, we use a density matrix formalism.
Our simulation can be broken down into three spatial regions: propagation through the galaxy cluster hosting the blazar; propagation through the intergalactic medium (IGM); and propagation through the Milky Way (MW).The principal difference between each region is the form of the magnetic field present.Photon to ALP conversion occurs most readily in the galaxy cluster where the VHE-photon density and magnetic fields are large.As the VHE photon/ALP beam propagates out of the cluster, the relatively strong cluster fields give way to a much weaker intergalactic field, largely suppressing any ALP-photon processes.The vast IGM is instead responsible for the attenuation of photons by EBL scattering.Finally, once the blazar jet reaches the Milky Way, we expect to see ALP-photon back-conversion induced by the strong galactic magnetic fields.
The effect of each propagation region is depicted in Fig. 5 which shows the photon survival probability of a photon with E = 400 GeV as a function of the distance of propagation from the blazar.The solid black line indicates the scenario without ALPs (N = 0), while the red and blue shows the effect of N = 1 and N = 20 ALP states.In the no-ALP scenario, the survival probability remains unity until propagation through the IGM where, due to EBL scattering, it decreases to ∼ 0.2.For the cases that include ALPs, a number of model unknowns associated with oscillation lead us to consider a large set possible P γ→γ -each element generated with a unique B-field structure and, for N > 1, an anarchical choice of {g γ i }.The central third of these {P γ→γ } are indicated by the filled regions in the figure.In both the one and 20 ALP cases, we see the effect of photon to ALP conversion in the galaxy cluster.In the one ALP case, we see a significant increase in the photon survival probability as the beam travels through the Milky Way, corresponding to the reconversion of ALPs to photons.However, this increase is not present in the 20 ALP case.This is because the electromagnetic ALP produced in the galaxy cluster is very likely to oscillate into a hidden ALP before it reaches the Milky Way.Due to oscillation into hidden ALPs, if the ALP-photon coupling is spread over 20 ALPs, we always expect a decrease in the blazar luminosity rather than the increase seen in the single ALP case.

The galaxy cluster
The simulation begins with a jet of pure photons at the site of the blazar.These photons are unpolarised, hence we take as our initial state to be the following the density matrix: where the first two diagonal elements correspond to the photons' two polarisations, and the remaining diagonal elements correspond to the ALP mass eigenstates; off-diagonal elements are associated with a superposition of states.To facilitate comparison between the multi-ALP and single ALP effects, we follow [73] in our description of the magnetic fields.We take the cluster field to have a domain-like structure with a radially dependent magnitude given by: where, the core radius is r core = 200 kpc and the central field strength is B C 0 = 10 µG and η = 0.5.Within a given magnetic domain, the field is assumed to be constant.The direction of the magnetic field is randomised in each domain.The coherence of the magnetic field, therefore, depends on the domain length, which is taken to be ∆L c = 10 kpc.The total radius of the cluster is 2 Mpc.
Galaxy clusters are home to a large population of charged particles that will give an effective mass to the photon.To accurately describe these, we consider the electron density given by: where n 0 el = 10 −2 cm −3 .With this description of the cluster, we evolve the state matrix of Eq. (6.1) by constructing an evolution operator (G) for each domain: where H denotes the propagation Hamiltonian with the following parameters: here, E is the photon/ALP energy (in GeV), n e is the electron density and B x and B y are the x and y components of the magnetic field strength (in µG) [74].The mass of the i th ALP state, denoted m i , is taken to be distributed logarithmically within 10 −8 , 10 −5 eV.The various components of H can be interpreted as follows: ∆ ϕ i are the ALP mass terms; ∆ pl is a photon effective mass terms, induced by thermal effects in the electron plasma; (∆ ϕγ x ) i and (∆ ϕγ y ) i are the ALP-photon couplings for each photon polarisation; and, ∆ QED implements vacuum polarisation effects.As we are working in the mass basis, each ALP is independent of all other ALP states so H is diagonal in the ALP sector.We construct a new evolution operator, G i , for each domain using the central values of B and n el .The state of the system after propagation through the cluster is then given by: where i runs over each domain in the cluster.

Intergalactic space
The intergalactic magnetic field is modelled with a similar domain-like structure to that of the AGN cluster, albeit with a significantly lower overall strength and a much larger domain length ∆L IG = 50 Mpc.The intergalactic medium has a luminosity redshift (z) dependent field strength given by: where B IG 0 = 1 nG.The intergalactic electron density is approximated with a constant value of n IG el = 10 −7 cm −3 .To account for VHE photon scattering off the EBL we introduce a non-unitary decay matrix (D) for each domain: where τ is the optical depth associated with propagation through a given domain.In general, τ depends on the photon energy and the redshift of the domain.We use the EBL opacity model given in Ref. [75].After computing D(τ i ) for each domain, the ALP-photon state is propagated as follows: where G IG i is calculated as before using Eq.(6.4) with i running over each domain in the IGM.

The Milky Way
Compared to the previous two field scenarios, the Milky Way field has a significantly more complex structure that comprises a halo field that surrounds a spatially compact disk field.Following Ref. [76], we describe the disk by a series of logarithmic spirals and the halo field by a superposition of piece-wise functions extending relatively far above and below the galactic plane, see Fig. 6.It is this extensive halo field that is responsible for the majority of ALP-Photon conversion.A detailed description of the field model used is provided in [76].We note that in our work, a small alteration has been made to correct the function describing boundaries between the consecutive log-spiral regions in the disk: where α open = 11.5 o is the opening angle of the log-spiral and r j are the radii at which each spiral boundary crosses the negative x-axis.Finally, we use a constant electron density of 0.1 cm −3 to construct the MW evolution operator.The final state of the system is then given by the product over domains i: where again we discretise our computation of the evolution operator (G i ) -in this case, using an approximate-coherence length of 100 pc.We begin the Milky Way propagation stage when the ALP/photon reaches a radial distance of 20 kpc from the galactic centre.The probability that an emitted photon is detected as a photon on Earth is found by projecting the final state (ρ M W out ) onto the two possible photon states: .12)

Results
Due to observational limitations, the field orientations in each domain of the galaxy cluster and intergalactic medium are unknown and small changes thereto could have a large effect on the final value of P γ→γ .It is thus necessary to simulate with a large set of different field directions.P γ→γ now represents a stochastic quantity that we use to marginalise over the unknown field states.As with the rest of this work, we also encounter an ambiguity in the choice of mixing parameters ({g γ i }), which we again account for using the ALP anarchy model described above.Each simulation run has a unique set of coupling parameters and magnetic field orientations.The final results described below were obtained using a sample of N samp = 1000 simulation runs and a set of 6 Blazar sources.These sources (listed in Table 1) were selected based on the availability of their data and on their inclusion in other works on this subject [73].We plan to perform a more detailed follow-up study using a larger sample of sources and more recent data.
For a given choice of blazar source and grid point (N, g γ ) in model parameter space, the result of our simulation is a set of N samp survival probability spectra P γ→γ (E), where E is the photon energy.To marginalise over these, we use the students p-statistic (p t ), following [73]; the calculation of p t is described in Appendix C. For each grid point, we obtain a distribution of N samp p-statistics.We collate these data by determining the value of p t corresponding to the field configuration, resulting in a better agreement between the model and corrected spectra than 95% of other configurations.We denote this as p 95 .In this way, we arrive at a two-dimensional parameter scan for each source that preserves more of the underlying p t distribution than could be achieved, for example, by simply averaging the p t over N samp .This method facilitates comparison with Ref. [73].Note that this method results in relatively conservative bounds on ALP scenarios as we are choosing a B field configuration that agrees rather well with the data.
Finally, to construct an overall parameter scan we determine p 95 on the union of the sample data for each source.The resulting grid scan is shown as a heat map in Fig. 7. Lower values of − log 10 (p 95 ) correspond to the model better reproducing the observations.We see a significant improvement over the Standard Model with the addition of a single ALP state, this corroborates the findings of Ref. [73]; the introduction of ALP states can alleviate tensions in standard intergalactic opacity models.Interestingly, and novel to this work, we also observe that increasing N provides a worse fit to the VHE data.As the number of ALP states increases, so too does the number of hidden states.Working in the interaction basis, states in models with many ALPs will exist less frequently as the interaction state.We might expect, therefore, a suppression in the back conversion probability.This can be seen explicitly in Fig. 5, where in the 20 ALP case, we do not observe any increase in P γ→γ in the MW.This effect can be very significant, in this case, completely negating the opacity decrease caused by reduced EBL scattering; P γ→γ is significantly lower in the 20 ALP case than it is for no ALPs (black line in the Fig. 5).It should be noted that the results depicted in Fig. 5, having been computed for a single source and energy, may not be representative.In contrast, Fig. 7 reflects our statistical analysis using multiple blazars and energies, where the single ALP state shows the greatest agreement with observed spectra.In addition to favouring fewer ALP states, we also see a lower degree of tension at greater ALP couplings.The lower regions of Fig. 7 behave as no ALP models, and there is no significant reduction in EBL scattering.

Discussion and conclusions
In this paper, we have explored the phenomenology of ALP anarchy models comprising many axion-like particles whose masses and Standard Model couplings are related by random matrices.String compactifications typically generate many ALPs; therefore, the phenomenology of many ALP systems is an important direction in studying physics beyond the Standard Model.The ALP anarchy scenario provides a benchmark for this phenomenology.
We have shown that a key feature of many ALP phenomenology is oscillations between the ALP states that couple to the Standard Model and hidden ALP states that do not.A given ALP model can be characterised by the total ALP photon and ALP electron couplings g γ 2 = g γ i 2 and g e2 = g e i 2 .However, oscillation into hidden ALP states can significantly reduce the signal in some ALP searches, such as CAST and IAXO.Other ALP search strategies sensitive only to photon disappearance into ALP degrees of freedom, such as arguments from stellar cooling, will be largely insensitive to the number of ALP fields for a given total ALPphoton coupling.Still, other search strategies such as ADMX [77] and other Dark Matter haloscopes rely on a mass resonance and would, therefore, be sensitive to each ALP mass eigenstate individually rather than to the total effective ALP couplings.If the total ALPphoton coupling and dark matter density is shared over many ALP states, the expected signal in haloscope experiments for a given mass would be correspondingly reduced.
As discussed in section Section 6, ALPs have been proposed as a solution to increase the observed transparency of the Universe to very high energy photons.However, we found that for many ALP systems, the photon survival probability instead decreases due to the presence of ALPs -the opposite effect to that observed for a single ALP.We, therefore, conclude that many ALP phenomenology leads to a number of new effects not captured by consideration of a single ALP field.

B Producing CAST bounds
This section provides a detailed overview of the approach used to recast the CAST bounds in Section 4.4.The original CAST bound can be generalised by considering the detector's sensitivity.From any point (g γ N =1 , g e N =1 ) on the original bound, the black line shown in Fig. 2, we can extract a proxy for the maximum signal (σ max ) that is consistent with non-detection: where Φ tot is the total single ALP flux at the detector; computed using Eqs.For any point in the coupling parameter space, the total flux at the detector must, therefore, be less than: Φ max (g γ ) = σ max /g γ 2 (B.3) We note that to obtain the total flux at the detector we integrated over the differential fluxes of Eqs.(4.1), (4.2) and (4.3) using a lower integration boundary of ω = 1 keV (corresponding to the detector threshold) and an upper bound of ω = 10 keV.
The recast bounds in Fig. 2 are shown as heat maps.The weight assigned to each set of couplings (g e , g γ ), corresponds to the proportion of our sample of survival and conversion probabilities, {P γ→γ } and {P e→γ }, that are consistent with Eq. (B.3).That is to say that, for a given (g e , g γ ) and probability P l e→γ ∈ {P e→γ }, we find the proportion (η) of {P γ→γ } that satisfy: Φ max (g γ ) ≥ Φ P (g γ )P γ→γ + (Φ B (g e ) + Φ C (g e ))P e→γ (B.4) At this point, η could be found by Monte Carlo sampling {P γ→γ } and determining the number of elements that satisfy Eq. (B.4).This approach is, however, very computationally expensive.A significant saving can be achieved by instead considering the critical value (P Crit ) of P γγ that achieves equality in Eq. (B.4): P Crit (g e , g γ , P l eγ ) = 1 Φ P (g γ ) Φ max (g γ ) − (Φ B (g e ) + Φ C (g e ))P l e→γ , (B.5) η is now given by the cumulative sum -evaluated at P Crit -of the histogram of {P γ→γ }; η = C(P Crit ).Repeating this process for every P l e→γ ∈ {P e→γ }, the total weight associated with a coupling pair is given by: W (g e , g γ ) = 1 N ae l C(P Crit (g e , g γ , P l e→γ )) , (B.6) where N ae is the size of {P e→γ } -taken here to be N ae = N aγ = 10 4 .W (g e , g γ ) is equivalent to the fraction of {P γ→γ } and {P e→γ } that satisfy Eq. (B.4) and will saturate to 1 for sufficiently low couplings.
C Fit quality for VHE spectra and statistical approach to multi-ALP parameter space To assess the fit of our multi-ALP model to the VHE Blazar spectra, we employ the same statistical methods as Refs.[73,78].In particular, we use the observed spectra, Φ obs i , shown in Table 1.We note that Φ obs i contains data on the blazar source plus all ambient backgrounds.We compute the photon survival probability for a given point in the theory parameter space: where E is the photon energy, and this probability can be computed on a bin-by-bin basis.
The absorption-corrected flux is given by Φ i (g γ , g e ) = ⟨P γ→γ ⟩ −1 i (g γ , g e )Φ obs i , (C.2) where again, we have kept the theory parameter dependence explicit.For a given point in the model parameter space, (g γ , g e ), we are tasked with understanding its compatibility with the observed spectra.For this purpose, we must quantify the observed flux in the absence of new physics and do so by fitting Φ i (g γ , g e ) to a power law: f (E) = N 0 (E/E 0 ) −Γ , p fit ⩾ 0.05 N 0 (E/E 0 ) −(Γ+βc ln(E/E 0 )) , otherwise (C.3) where p fit denotes the fit probability.The first power law contains three parameters: N 0 , E 0 , and Γ while in the second (which is a parabola in log-log space), there are four fit parameters, N 0 , E 0 , β c and Γ.The power law function is used unless the baseline (no ALP) fit probability (as calculated by Ref. [73]) is less than 0.05.We note that, when computing the fit parameters, we perform a Chi-Squared analysis, taking the spectral uncertainties to be statistical.For each energy bin (E i ), the fit residual is calculated: where σ i denotes the statistical measurement uncertainty at 68% confidence level.The value of χ 2 can now easily be obtained by summing the residuals over all energy bins and squaring.
Assuming that P γ→γ accurately predicts the Universe's opacity to VHE γ rays, we would expect the residuals in the optically thick regime to follow a Gaussian distribution with a mean of zero.To test this hypothesis, we employ the t-test, where we calculate the test statistic, t, as follows: where χ denotes the mean and σ χ the variance of the residual distribution, which comprises a total of N χ data points.The p t test statistic we use to quantify the fit of our model to the observed spectra can be obtained using a standard t-test procedure.

Figure 1 :
Figure 1: Top left image shows the Primakoff process that produces the electromagnetic ALP, and the middle and top right images show the Compton scattering and axio-recombination processes that produce the electronic ALP.The bottom row shows axio-deexcitation and Bremsstrahlung processes that produce the electronic ALP.

Figure 2 :
Figure 2: The fraction of realisations consistent with non-detection as a function of coupling (g e and g γ ) shown for two different values of N -Left: N = 2; Right: N = 30.The grey region indicates the excluded region from solar neutrinos [62].

Figure 3 :
Figure 3: Bounds on photon ALP coupling (g γ ) as a function of number of ALPs (N ) for g e = 10 −15 .The scatter plot depicts the minimum upper bound consistent with 50% of mixing matrix realisations.The fit to that is shown as a line plot.

Figure 4 :
Figure4: The bounds on g γ g e as a function of N for the low mass (m a < 10 −6 eV), low propagation limit of the MWD.The bounds were computed for a set 1000 coupling pairs (g e and g γ ).The central 90% of bounds lie within the red band, with the central third being encompassed by the blue band.The one ALP bound was recast from Ref.[63].

Figure 5 :
Figure 5: Photon survival probability against propagation distance for a photon energy of 400 GeV produced by 1ES0414+009.The zero ALP case is shown in black, with the central third of samples shown in red and blue for the 1 ALP and 20 ALP cases respectively.

Figure 6 :
Figure 6: The magnitude of the Milky Way field in two perpendicular planes through the galactic centre.Left: cross-section at z = 0 looking upon the galactic disk; Right: crosssection perpendicular to disk at y = 0.The sign of the field indicates the orientation of the azimuthal component: positive -anticlockwise (looking down negative z); negative -clockwise.The location of the earth, (−8.5, 0, 0)kpc, is marked by the black circle.

Figure 7 :
Figure 7: A heat map showing p 95 computed using 6 sources and N samp = 1000 samples per source.Lower values (higher p 95 ) indicate better agreement between the model and the data.

Table 1 :
List of VHE gamma-ray sources used in our analysis and their corresponding fit functions; Log Parabola (LP) or Power Law (PL), see Appendix C.