# ALP Anarchy

Francesca Chadha-Day  James Maxwell and  Jessica Turner

###### Abstract

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.

Preprint: IPPP/23/64

##  1 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]. Second, the
QCD axion offers a promising solution to the strong CP problem [3, 4].
Finally, string theory suggests the existence of a large number of ALPs [5, 6,
7, 8, 9, 10]. 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 $\phi$ interacting with photons
and electrons is given by

$\mathcal{L}\supset-\frac{1}{2}\partial^{\mu}\phi\partial_{\mu}\phi-\frac{1}{2}m^{2}\phi^{2}-g^{\gamma}\phi\tilde{F}^{\mu\nu}F_{\mu\nu}+\frac{g^{e}}{2m_{e}}\bar{\psi}\gamma^{\mu}\gamma_{5}\psi\partial_{\mu}\phi\,,$
(1.1)

where $g^{\gamma}$ is the coupling of $\phi$ to the electromagnetic field
strength tensor, $F^{\mu\nu}$, $g^{e}$ is the dimensionless ALP-electron
coupling, $m_{e}$ is the electron mass and $\psi$ 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 [11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25]. The goal of this paper is to explore the interaction of
many ALP models with the Standard Model. As pointed out in Ref. [23], 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 [26]. This paper will expand on the results of Ref. [23] 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.

##  2 ALP oscillations

We will consider a model containing $N$ ALPs, each coupling to both photons
and electrons:

$\mathcal{L}\supset\sum_{i=1}^{N}\left(-\frac{1}{2}\partial^{\mu}\phi_{i}\partial_{\mu}\phi_{i}-\frac{1}{2}m_{i}^{2}\phi_{i}^{2}-g^{\gamma}_{i}\phi_{i}\tilde{F}^{\mu\nu}F_{\mu\nu}+\frac{g^{e}_{i}}{2m_{e}}\bar{\psi}\gamma^{\mu}\gamma_{5}\psi\partial_{\mu}\phi_{i}\right)\,,$
(2.1)

where $\phi_{i}$ is an ALP field with mass $m_{i}$, $g^{\gamma}_{i}$ is the
coupling of $\phi_{i}$ to the electromagnetic field strength tensor,
$F^{\mu\nu}$, $g^{e}_{i}$ is the dimensionless ALP-electron coupling, $m_{e}$
is the electron mass and $\psi$ 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 ultra-light
ALPs, and hence, we will assume that the QCD axion is too heavy to contribute
to the scenarios considered here. It will be convenient to rotate to the
electromagnetic basis, in which only a single ALP (the ‘electromagnetic ALP’)
couples to photons:

$\mathcal{L}\supset-\sum_{i}\frac{1}{2}\partial^{\mu}\phi_{i}\partial_{\mu}\phi_{i}-\sum_{i,j}\frac{1}{2}M^{\gamma}_{ij}\phi_{i}\phi_{j}+\sum_{i}\frac{g^{e}_{i}}{2m_{e}}\bar{\psi}\gamma^{\mu}\gamma_{5}\psi\partial_{\mu}\phi_{i}-g^{\gamma}\phi_{1}\tilde{F}^{\mu\nu}F_{\mu\nu}\,,$
(2.2)

where $g^{\gamma}=\sqrt{\sum_{i}{g^{\gamma}_{i}}^{2}}$ is the total effective
ALP-photon coupling, we have chosen
$\phi_{1}={\sum_{i}g^{\gamma}_{i}\phi_{i}}/{g^{\gamma}}$ as the
electromagnetic ALP and $\phi_{i}$ and $g^{e}_{i}$ have been appropriately
redefined. This basis has the advantage that only $\phi_{1}$ interacts
directly with the photon. The other ALPs, $\phi_{\\{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
$\phi_{1}$. The production rate can, in this case, be calculated simply by
considering a single ALP with coupling $g^{\gamma}$ 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^{\gamma}$. As we assume that the ALP mass
states differ and that the mass and electromagnetic basis are misaligned,
$\phi_{1}$ may mix with the hidden states, $\phi_{\\{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:

$\mathcal{L}\supset-\sum_{i}\frac{1}{2}\partial^{\mu}\phi_{i}\partial_{\mu}\phi_{i}-\sum_{i,j}\frac{1}{2}M^{e}_{ij}\phi_{i}\phi_{j}-\sum_{i}g^{\gamma}_{i}\phi_{i}\tilde{F}^{\mu\nu}F_{\mu\nu}+\frac{g^{e}}{2m_{e}}\bar{\psi}\gamma^{\mu}\gamma_{5}\psi\partial_{\mu}\phi_{1}\,,$
(2.3)

where $g^{e}=\sqrt{\sum_{i}{g^{e}_{i}}^{2}}$ and we have now chosen
$\phi_{1}={\sum_{i}g^{e}_{i}\phi_{i}}/g^{e}$ as the electronic ALP and
$\phi_{i}$ and $g^{\gamma}_{i}$ have been appropriately redefined. As with the
electromagnetic ALP, this basis has the advantage that only $\phi_{1}$
interacts directly with the electron while the other ALPs $\phi_{\\{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 [27, 28, 29], 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.

##  3 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 [30]. Nonetheless, the ALP photon coupling has been modelled in a
range of string axiverse scenarios [31, 10]. 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 $\mathcal{O}(1)$ parameters [32, 33, 34, 35].
Randomness in $\mathcal{O}(1)$ coupling constants is expected in sufficiently
complicated models or with many fields mixing with each other. 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^{\gamma}$ and $M^{e}$, are real and can, therefore, be related to
the mass basis states as follows:

$M^{\alpha}=U^{\alpha}D{U^{\alpha}}^{T}\,\quad\alpha=e,\gamma\,,$ (3.1)

where we assume $U^{\alpha}\in SO(N)$ and
$D=\operatorname{diag}\left(m_{1},m_{2},\ldots,m_{N}\right)$ is a real
diagonal matrix with $m_{i}$ denoting the mass of ALP field $\phi_{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^{\alpha}$ \- 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:

$dV=\left(\prod_{i\in[1,N-1],j\in[i,N-1]}\sin^{i-1}\theta_{i,j}\right)d\theta_{1,1}\cdots
d\theta_{N-1,N-1}\quad\,,$ (3.2)

where there are $N(N-1)/2$ mixing angles, $\left\\{\theta_{ij}\right\\}_{1\leq
i\leq j\leq N-1}$ with $\theta_{1,j}\in[0,2\pi]$ and $\theta_{i,j}\in[0,\pi]$
for $i,j>1$. Sampling uniformly in $dV$ yields the desired distribution of the
angles that parameterise the mixing matrices $U^{\gamma}$ 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}^{\gamma}\\}$ and $\\{g_{i}^{e}\\}$. Given, for example, the EM
coupling in the EM-ALP basis, the mass basis couplings are given by:

$\begin{pmatrix}g_{1}^{\gamma}\\\ g_{2}^{\gamma}\\\
\vdots\end{pmatrix}=U^{\gamma}\begin{pmatrix}g^{\gamma}\\\ 0\\\
\vdots\end{pmatrix}.$ (3.3)

An analogous statement is also true for the electron-ALP couplings. From Eq.
(3.3), it is clear that the top row of $U^{\gamma}$ 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^{\gamma}$ 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.

##  4 The Cern Axion Solar Telescope

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) [36, 37], axio-
recombination and axio-deexcitation [38, 39, 40], and Compton scattering [41,
42, 43], see Ref. [44] 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\,\rm{m}$ long evacuated cylinder with a sustained magnetic field oriented
transverse to the ALP propagation direction. A bound on $g^{\gamma}$ 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.

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.

###  4.1 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
$\mathrm{m}^{-2}\,\text{year}^{-1}\,\mathrm{keV}^{-1}$, generated by these
mechanisms are given by [45]:

$\displaystyle\left.\frac{\mathrm{d}\Phi_{a}}{\mathrm{\leavevmode\nobreak\
d}\omega}\right|_{B}$ $\displaystyle=8.3\times
10^{20}\left(\frac{g_{e}}{10^{-13}}\right)^{2}\frac{\omega}{1+0.667\omega^{1.278}}e^{-0.77\omega}\,,$
(4.1)
$\displaystyle\left.\frac{\mathrm{d}\Phi_{a}}{\mathrm{\leavevmode\nobreak\
d}\omega}\right|_{C}$ $\displaystyle=4.2\times
10^{18}\left(\frac{g_{e}}{10^{-13}}\right)^{2}\omega^{2.987}e^{-0.776\omega}\,,$
(4.2)
$\displaystyle\left.\frac{\mathrm{d}\Phi_{a}}{\mathrm{\leavevmode\nobreak\
d}\omega}\right|_{P}$ $\displaystyle=2.0\times
10^{18}\left(\frac{g_{\gamma}}{10^{-12}\mathrm{GeV}^{-1}}\right)^{2}\omega^{2.450}e^{-0.829\omega}\,,$
(4.3)

where the ALP energy, $\omega$, is in $\mathrm{keV}$. Integrating over the
energy range ($0.8-6.8$ keV) of the CAST analysis [46], the total fluxes can
be found for $\Phi_{B}$, $\Phi_{C}$ and $\Phi_{P}$:

$\displaystyle\frac{\Phi_{B}}{\mathrm{m}^{-2}\,\text{year}^{-1}}=4.1\times
10^{46}{g^{e}}^{2}\,,$ (4.4)
$\displaystyle\frac{\Phi_{C}}{\mathrm{m}^{-2}\,\text{year}^{-1}}=5.2\times
10^{45}{g^{e}}^{2}\,,$ (4.5)
$\displaystyle\frac{\Phi_{P}}{\mathrm{m}^{-2}\,\text{year}^{-1}}=1.0\times
10^{43}\left({\frac{g^{\gamma}}{\mathrm{GeV}^{-1}}}\right)^{2}\,.$ (4.6)

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.

###  4.2 Oscillation: The two ALP case

We will expand upon the results of [23] by discussing the basics of ALP
oscillations relevant to CAST in a simplified two ALPs scenario where
$\phi_{e}$ and $\phi_{\gamma}$ are linear combinations of only two massive ALP
states, $\phi_{1}$ and $\phi_{2}$, which have masses $m_{1}$ and $m_{2}$
respectively. The Lagrangian in the mass basis is

$\displaystyle\mathcal{L}\supset$
$\displaystyle-\frac{1}{2}\partial^{\mu}\phi_{1}\partial_{\mu}\phi_{1}-\frac{1}{2}\partial^{\mu}\phi_{2}\partial_{\mu}\phi_{2}-\frac{1}{2}m_{1}^{2}\phi_{1}^{2}-\frac{1}{2}m_{2}^{2}\phi_{2}^{2}$
(4.7) $\displaystyle-
g_{1}^{\gamma}\phi_{1}\tilde{F}^{\mu\nu}F_{\mu\nu}-g_{2}^{\gamma}\phi_{2}\tilde{F}^{\mu\nu}F_{\mu\nu}+\frac{g^{e}_{1}}{2m_{e}}\bar{\psi}\gamma^{\mu}\gamma_{5}\psi\partial_{\mu}\phi_{1}+\frac{g^{e}_{2}}{2m_{e}}\bar{\psi}\gamma^{\mu}\gamma_{5}\psi\partial_{\mu}\phi_{2}\,.$

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

$\phi_{\gamma}=\frac{g_{1}^{\gamma}\phi_{1}+g_{2}^{\gamma}\phi_{2}}{\sqrt{g_{1}^{\gamma
2}+g_{2}^{\gamma
2}}}\,,\quad\phi_{\gamma_{h}}=\frac{g_{2}^{\gamma}\phi_{1}-g_{1}^{\gamma}\phi_{2}}{\sqrt{g_{1}^{\gamma
2}+g_{2}^{\gamma 2}}}\,,$ (4.8)

respectively. Analogously, the electronic ALP and the hidden electronic ALP
are the following orthogonal combinations:

$\phi_{e}=\frac{g_{1}^{e}\phi_{1}+g_{2}^{e}\phi_{2}}{\sqrt{g_{1}^{e2}+g_{2}^{e2}}}\,,\quad\phi_{e_{h}}=\frac{g_{2}^{e}\phi_{1}-g_{1}^{e}\phi_{2}}{\sqrt{g_{1}^{e2}+g_{2}^{e2}}}\,,$
(4.9)

respectively. Note that the electromagnetic and electron ALPs $\phi_{\gamma}$
and $\phi_{e}$ are generally neither colinear nor orthogonal. Therefore,
$\phi_{\gamma_{h}}$ will in general have some non-zero coupling to electrons
and $\phi_{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:

$\displaystyle\left(\begin{array}[]{l}\phi_{\gamma}\\\
\phi_{h_{\gamma}}\end{array}\right)=\left(\begin{array}[]{cc}\cos(\theta^{\gamma})&\sin(\theta^{\gamma})\\\
-\sin(\theta^{\gamma})&\cos(\theta^{\gamma})\end{array}\right)\left(\begin{array}[]{l}\phi_{1}\\\
\phi_{2}\end{array}\right)\,,$ (4.10)
$\displaystyle\left(\begin{array}[]{l}\phi_{e}\\\
\phi_{h_{e}}\end{array}\right)=\left(\begin{array}[]{cc}\cos(\theta^{e})&\sin(\theta^{e})\\\
-\sin(\theta^{e})&\cos(\theta^{e})\end{array}\right)\left(\begin{array}[]{l}\phi_{1}\\\
\phi_{2}\end{array}\right)\,,$

with

$\cos\left(\theta^{\alpha}\right)=\left(\frac{g_{1}^{\alpha}}{\sqrt{g_{1}^{\alpha
2}+g_{2}^{\alpha
2}}}\right)\,,\quad\sin\left(\theta^{\alpha}\right)=\left(\frac{g_{2}^{\alpha}}{\sqrt{g_{1}^{\alpha
2}+g_{2}^{\alpha 2}}}\right)\,,$ (4.11)

where $\alpha=\gamma,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}\,\rm{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 $\phi_{e}$ while Primakoff production produces the state
$\phi_{\gamma}$. Furthermore, in this mass range, the CAST sensitivity is
independent of mass, and therefore, CAST will detect the state $\phi_{\gamma}$
as a single signal.

As CAST aims to detect electromagnetic ALPs, $\phi_{\gamma}$, 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 $\ket{\phi_{i}(L)}={\rm
e}^{-i\frac{m_{i}^{2}L}{2\omega}}\ket{\phi_{i}(0)}$, where $L$ is the distance
travelled. The former is given by

$\displaystyle P(\phi_{e}\to\phi_{\gamma})\equiv P_{e\to\gamma}$
$\displaystyle=\frac{1}{2}\left(1+\cos(2\theta^{e})\cos(2\theta^{\gamma})+\sin(2\theta^{e})\sin(2\theta^{\gamma})\cos\left(\frac{\Delta
m^{2}L}{2\omega}\right)\right)\,,$ (4.12)

where $\Delta 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^{\alpha}$,
yields

$P_{e\to\gamma}=\frac{\left({g_{1}^{e}}^{2}g_{1}^{\gamma
2}+g_{2}^{e2}g_{2}^{\gamma
2}\right)}{\left({g_{1}^{e}}^{2}+g_{2}^{e2}\right)\left(g_{1}^{\gamma
2}+g_{2}^{\gamma
2}\right)}\left(1+\frac{2g_{1}^{e}g_{1}^{\gamma}g_{2}^{e}g_{2}^{\gamma}}{\left({g_{1}^{e}}^{2}g_{1}^{\gamma
2}+g_{2}^{e2}g_{2}^{\gamma 2}\right)}\cos\left(\frac{\Delta
m^{2}L}{2\omega}\right)\right)\,.$ (4.13)

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:

$\displaystyle P(\phi_{\gamma}\to\phi_{\gamma})\equiv P_{\gamma\to\gamma}$
$\displaystyle=1-\sin^{2}(2\theta^{\gamma})\sin^{2}\left(\frac{\Delta
m^{2}L}{4\omega}\right)$ (4.14) $\displaystyle=1-4\frac{g_{1}^{\gamma
2}g_{2}^{\gamma 2}}{\left(g_{1}^{\gamma 2}+g_{2}^{\gamma
2}\right)^{2}}\sin^{2}\left(\frac{\Delta m^{2}L}{4\omega}\right)\,.$

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\approx G_{F}N_{e}\sim
10^{-12}\,\mathrm{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
$\frac{g^{e}}{2m_{e}}=10^{-11}\,\mathrm{GeV}^{-1}$) is
$V\approx\left(\frac{g^{e}}{2m_{e}}\right)^{2}N_{e}\sim 10^{-29}\rm{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 $\Delta
m^{2}>10^{-12}\,\mathrm{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:

$P_{e\to\gamma}=\frac{\left({g_{1}^{e}}^{2}g_{1}^{\gamma
2}+g_{2}^{e2}g_{2}^{\gamma
2}\right)}{\left({g_{1}^{e}}^{2}+g_{2}^{e2}\right)\left(g_{1}^{\gamma
2}+g_{2}^{\gamma 2}\right)}\,,\quad
P_{\gamma\to\gamma}=\frac{{g^{\gamma}_{1}}^{4}+{g^{\gamma}_{2}}^{4}}{\left(g_{1}^{\gamma
2}+g_{2}^{\gamma 2}\right)^{2}}\,.$ (4.15)

###  4.3 Oscillation: The many ALP case

We now turn to the case where many ALP mass eigenstates couple to electrons
and photons:

$\mathcal{L}\supset\sum_{i}^{N}\left(\frac{1}{2}\partial^{\mu}\phi_{i}\partial_{\mu}\phi_{i}-\frac{1}{2}m_{i}^{2}\phi_{i^{2}}-g_{i}^{\gamma}\phi_{i}\tilde{F}^{\mu\nu}F_{\mu\nu}+g_{i}^{e}\bar{\psi}\gamma^{\mu}\gamma_{5}\psi\partial_{\mu}\phi_{i}\right)\,,$
(4.16)

such that the electromagnetic and electronic ALPs produced in the Sun are
linear combinations of the mass states:

$\phi_{\gamma}=\frac{\sum_{i}g_{i}^{\gamma}\phi_{i}}{\sqrt{\sum_{i}{g_{i}^{\gamma}}^{2}}}\,,\quad\phi_{e}=\frac{\sum_{i}g_{i}^{e}\phi_{i}}{\sqrt{\sum_{i}{g_{i}^{e}}^{2}}}\,,$
(4.17)

where we again assume that all ALP masses considered are
$m_{i}<10^{-2}\,\rm{eV}$. Any mass eigenstates with $m_{i}>10^{-2}\,\rm{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 $\Delta
m^{2}>10^{-12}\,\mathrm{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

$P_{\gamma\to\gamma}=\frac{1}{N}\left(1+\frac{N^{2}\mathrm{VAR}\left(\\{{g_{i}^{\gamma}}^{2}\\}\right)}{{g^{\gamma}}^{4}}\right)=\frac{\sum_{i}^{N}{g^{\gamma}_{i}}^{4}}{\left(\sum_{i}^{N}{g^{\gamma}_{i}}^{2}\right)^{2}}\,,$
(4.18)

where $N$ is the number of ALP mass eigenstates,
$g^{\gamma}=\sqrt{\sum{g_{i}^{\gamma}}^{2}}$ is the coupling of the
electromagnetic ALP to photons, and $\mathrm{VAR}(\\{{g^{\gamma}_{i}}^{2}\\})$
is the variance of the ALP-photon couplings squared in the mass basis which is
given by

$\mathrm{VAR}\left(\\{{g_{i}^{\gamma}}^{2}\\}\right)=\frac{\sum_{i}^{N}\left[{{g^{\gamma}_{i}}}^{2}-\frac{{{g^{\gamma}}}^{2}}{N}\right]^{2}}{N}\,.$
(4.19)

Likewise, the probability that an electronic ALP oscillates to an
electromagnetic ALP is

$P_{e\to\gamma}=\frac{1}{N}\left(1+\frac{N^{2}\mathrm{COVAR}\left(\\{{g_{i}^{\gamma}}^{2}\\},\\{{g_{i}^{e}}^{2}\\}\\}\right)}{{g^{\gamma}}^{2}{g^{e}}^{2}}\right)=\frac{\sum_{i}^{N}{g^{e}_{i}}^{2}{g^{\gamma}_{i}}^{2}}{\sum_{i}^{N}{g^{e}_{i}}^{2}\sum_{i}^{N}{g^{\gamma}_{i}}^{2}}\,,$
(4.20)

where $g^{e}=\sqrt{\sum{g_{i}^{e}}^{2}}$ is the coupling of the electronic ALP
to electrons and
$\mathrm{COVAR}\left(\\{{g_{i}^{\gamma}}^{2}\\},\\{{g_{i}^{e}}^{2}\\}\\}\right)$
is the covariance of the ALP-photon and ALP-electron couplings squared and is

$\mathrm{COVAR}\left(\\{{g_{i}^{\gamma}}^{2}\\},\\{{g_{i}^{e}}^{2}\\}\\}\right)=\frac{\sum_{i}^{N}\left[{g^{\gamma}_{i}}^{2}-\frac{{g^{\gamma}}^{2}}{N}\right]\left[{g^{e}_{i}}^{2}-\frac{{g^{e}}^{2}}{N}\right]}{N}\,.$
(4.21)

###  4.4 Reintepretation of CAST results

The non-observation of an excess number of photons allows CAST to place
constraints on the $(g^{\gamma},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:

$\displaystyle{\Phi^{\rm{osc}}_{\gamma}}$
$\displaystyle=P_{\gamma\to\gamma}\Phi_{\gamma}+P_{e\to\gamma}\Phi_{e}\,,$
(4.22) $\displaystyle\Phi_{e}^{\rm{osc}}$ $\displaystyle=P_{\gamma\rightarrow
e}\Phi_{\gamma}+P_{e\rightarrow e}\Phi_{e}\,,$

where $\Phi_{\gamma}=\Phi_{P}$ and $\Phi_{e}=\Phi_{B}+\Phi_{C}$ and
$P_{a\rightarrow 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
$\Phi_{\gamma}^{\rm{osc}}$. CAST has placed bounds on $(g^{\gamma})^{2}$ as a
function of ALP mass and on $g^{\gamma}$ as a function of $g^{e}$, for small
ALP mass $m_{a}$, for a single ALP. Here we determine an upper bound on
$g^{\gamma}$ as a function of $g^{e}$ for multi-ALP scenarios with
$m_{a}\leq\,10^{-2}\,{\rm eV}$ and $\Delta m^{2}>10^{-12}\,{\rm eV}^{2}$ for
all relevant ALP mass eigenstates; for which the bound becomes mass
independent [46]. 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^{\gamma}_{N=1}(g^{e})$,
(shown by the black line in Fig. 2), the maximum flux of electromagnetic ALP
on Earth is bounded to be:

$\Phi^{\mathrm{Max}}(g^{e},g^{\gamma})=\frac{(g^{\gamma}_{N=1}(g^{e}))^{2}}{{g^{\gamma}}^{2}}\Phi^{N=1}(g^{e},g^{\gamma}_{N=1}(g^{e}))\,,$
(4.23)

where $\Phi^{N=1}=\Phi_{P}+\Phi_{B}+\Phi_{C}$ is the flux at the detector in
the single ALP case. With ${\Phi_{\gamma}^{\rm{osc}}}(g^{e},g^{\gamma})$ in
hand, we now perform a grid scan over the coupling parameter space,
($g^{e},g^{\gamma}$), 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
$\Phi^{\mathrm{Max}}$.

To determine the bound on the total effective couplings in the ALP anarchy
scenario, it is therefore necessary to compute $P_{\gamma\to\gamma}$ and
$P_{e\to\gamma}$. 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_{i}^{e}$ and $g_{i}^{\gamma}$). It is this
relationship that encodes the mixing between the various ALP states. We
consider $10^{4}$ realisations of $\\{g^{\gamma}_{i}\\}$ and $\\{g^{e}_{i}\\}$
for each overall coupling pair ($g^{\gamma},g^{e}$). For each point in
($g^{\gamma},g^{e}$) space, we determine the proportion of viable realisations
such that $\Phi_{\gamma}^{\rm{osc}}<\Phi^{\mathrm{Max}}(g^{e},g^{\gamma})$.
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
($\phi_{e}\equiv\phi_{\gamma}$) where in the mass independent region
$\log_{10}\left(g^{\gamma}\,[\rm{GeV}^{-1}]\right)\lesssim-10$. The boundary
between regions where $0$ and $100\%$ of realisations are viable can be
interpreted as the bound on $g^{\gamma}$ as a function of $g^{e}$ in the ALP
anarchy scenario.

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

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^{\gamma}$ is
relaxed by almost half an order of magnitude with
$\log_{10}(g^{\gamma}\,[\rm{GeV^{-1}}])\lesssim-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^{\gamma}$ can increase to
compensate for this decrease in the detectable flux.

To quantify this relationship, we consider ALP multiplicities
$N\in\left[2,30\right]$. For each $N$ in this set, we determine the value of
$g^{\gamma}$, $g^{\gamma}_{50}(N)$ for which $50\%$ of mixing realisations
satisfying $\Phi_{\gamma}^{\rm{osc}}<\Phi^{\mathrm{Max}}(g^{e},g^{\gamma})$ 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 $\phi_{e}$,
namely bremsstrahlung and Compton scattering, are ineffective, and the
production of $\phi_{\gamma}$ in the Sun dominates. $g^{\gamma}_{50}(N)$ can
be interpreted as an approximate bound on $g^{\gamma}$ in the ALP anarchy
scenario. This result is shown in Fig. 3, which we fit the following function
to:

$g^{\gamma}=e^{m\log N-c},$ (4.24)

where we find $m=0.25$ and $c=-23$. This fitting function is also shown in
Fig. 3.

Figure 3: Bounds on photon ALP coupling ($g^{\gamma}$) 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.

We can understand the dependence of the bound on $N$ \-
$g^{\gamma}_{50}(N)\propto N^{1/4}$. We first note that when electromagnetic
ALP production dominates, CAST is sensitive to ${g^{\gamma}}^{4}$, as the ALP
must be produced and detected via this coupling to photons. In the many ALP
scenario, the bound placed on ${g^{\gamma}}^{4}$ is weakened in proportion to
the electromagnetic ALP survival probability $P_{\gamma\rightarrow\gamma}$
given in Eq. (4.18). For large $N$, the variance term becomes negligible and
we have $P_{\gamma\rightarrow\gamma}\sim\frac{1}{N}$. We therefore obtain
$g^{\gamma}_{50}(N)\propto 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 $\Delta
m^{2}>10^{-12}\,\rm{eV}^{2}$. As shown in [23], for $\Delta
m^{2}<10^{-14}\,\rm{eV}^{2}$, there is no significant oscillation into hidden
states. For intermediate mass differences, there is a non-trivial oscillation
structure depending on $\omega$.

##  5 Magnetic white dwarfs

We now consider limits on ALP couplings from observations of magnetic white
dwarfs (MWDs). MWDs produce $\omega\sim$ 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. [48].
More specifically, in the single ALP scenario, where $N=1$ and
$\phi_{e}\equiv\phi_{\gamma}$, the $(g^{e},g^{\gamma})$ parameter space is
constrained from the non-observation of astrophysical X-ray emission from the
MWD RE J0317-853 by the Suzaku telescope [49]. The flux of ALP-induced X-ray
photons on Earth, in the low ALP mass regime, is approximately proportional to
$\Phi_{\text{X-ray}}\propto(g^{e}g^{\gamma})^{2}$ as the ALP luminosity is
proportional to ${g^{e}}^{2}$ while the probability the ALP transitions to
X-ray photons is proportional to ${g^{\gamma}}^{2}$. The non-observation of
excess X-rays provides an upper bound on $\Phi_{\text{X-ray}}$ and hence a
corresponding bound on $(g^{e},g^{\gamma})$.

In the multi-ALP case, since the radius of the MWD is relatively small (less
than a percent of the Sun’s radius [50]), the oscillations between the
electronic ALP and hidden states do not have time to develop, assuming $\Delta
m^{2}\lesssim 4R=10^{-10}\,\rm{eV}^{2}$, where $R$ is the MWD’s radius.
Moreover, we assume that all ALP masses considered are
$m_{i}<10^{-2}\,\rm{eV}$. To perform a simple recast of the single to the
multi-ALP case, we compute the conversion probability
($P_{e\rightarrow\gamma}$) and scale the $N=1$ bound on $g^{e}g^{\gamma}$,
denoted as $b_{N=1}$, as follows:

$b_{N>1}=b_{N=1}/P_{e\rightarrow\gamma}$ (5.1)

where $b_{N>1}$ is the new bound on $g^{e}g^{\gamma}$ 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\rightarrow\gamma}$ is given by
the scalar product of $\phi_{e}$ and $\phi_{\gamma}$:

$P_{e\rightarrow\gamma}=\left|\frac{\sum_{i}g_{i}^{e}g_{i}^{\gamma}}{g^{e}g^{\gamma}}\right|^{2}$
(5.2)

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\times 10^{-3}$ Solar Radii [50],
$E=5\,\mathrm{keV}$ and masses distributed logarithmically in
$[10^{-9},10^{-6}]\,\mathrm{eV}$. The conversion probability, accounting for
the presence of hidden ALPs (in the simple two-ALP case, see Eq. (4.12)), is
given by

$P_{e\to\gamma}=\mathrm{cos}^{2}(\theta^{e}-\theta^{\gamma})\,,$ (5.3)

which describes the projection of the electronic ALP state onto the EM ALP
state. In the limit $\theta^{e}=\theta^{\gamma}$, naturally, the probability
$P_{e\to\gamma}=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\to\gamma}$ can be significantly smaller and make
the bound proportionally weaker.

Figure 4: The bounds on $g_{\gamma}g_{e}$ as a function of $N$ for the low
mass ($m_{a}<10^{-6}\mathrm{eV}$), low propagation limit of the MWD. The
bounds were computed for a set $1000$ coupling pairs ($g^{e}$ and
$g^{\gamma}$). 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. [48].

In Fig. 4, we show how the limit on $\lvert g^{\gamma}g^{e}\rvert$, in the low
ALP mass regime, from Ref. [48] 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 decreases the typical
overlap between the electromagnetic and electronic ALP states, therefore
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$.

##  6 Very high energy blazars

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 [51, 52], although the evidence
for this effect is not conclusive [53, 54]. Several phenomenological studies
have introduced ALPs to explain this discrepancy [55, 56, 57, 58]. 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. [58], we consider an arbitrarily mass-mixed set of ALP
states.

###  6.1 Simulation

Figure 5: Photon survival probability against propagation distance for a
photon energy of $400\,\rm{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.

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\,\rm{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 $\sim 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_{\gamma\rightarrow\gamma}$ – each element generated with a unique B-field
structure and, for $N>1$, an anarchical choice of $\\{g^{\gamma}_{i}\\}$. The
central third of these $\\{P_{\gamma\rightarrow\gamma}\\}$ 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.

####  6.1.1 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:

$\rho_{\gamma}=1/2\begin{bmatrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&0&0\\\
0&0&0&\ddots\,\end{bmatrix}\,,$ (6.1)

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 [58] 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:

$B^{C}(r)=B_{0}^{C}(1+(r/r_{\mathrm{core}})^{2})^{-\eta}\,,$ (6.2)

where, the core radius is $r_{\mathrm{core}}=200\,\mathrm{kpc}$ and the
central field strength is $B_{0}^{C}=10\,\mu\mathrm{G}$ and $\eta=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 $\Delta
L_{\mathrm{c}}=10\,\mathrm{kpc}$. The total radius of the cluster is
$2\,\mathrm{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:

$n_{\mathrm{el}}(r)=n_{\mathrm{el}}^{0}(1+r/r_{\mathrm{core}})^{-1}\,,$ (6.3)

where $n_{el}^{0}=10^{-2}\mathrm{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:

$\displaystyle G$ $\displaystyle=e^{iH\Delta L}\,,\quad
H=\begin{bmatrix}\Delta^{\mathrm{pl}}+2\Delta^{\mathrm{QED}}&0&(\Delta^{\phi\gamma}_{x})_{1}&(\Delta^{\phi\gamma}_{x})_{2}&(\Delta^{\phi\gamma}_{x})_{3}\\\
0&\Delta^{\mathrm{pl}}+\frac{7}{2}\Delta^{\mathrm{QED}}&(\Delta^{\phi\gamma}_{y})_{1}&(\Delta^{\phi\gamma}_{y})_{2}&(\Delta^{\phi\gamma}_{y})_{3}\\\
(\Delta^{\phi\gamma}_{x})_{1}&(\Delta^{\phi\gamma}_{y})_{1}&\Delta_{1}^{\phi}&0&0\\\
(\Delta^{\phi\gamma}_{x})_{2}&(\Delta^{\phi\gamma}_{y})_{2}&0&\Delta_{2}^{\phi}&0\\\
(\Delta^{\phi\gamma}_{x})_{3}&(\Delta^{\phi\gamma}_{y})_{3}&0&0&\ddots\\\
\end{bmatrix}\,,$ (6.4)

where $H$ denotes the propagation Hamiltonian with the following parameters:

$\displaystyle\Delta^{\mathrm{pl}}$ $\displaystyle=\frac{-1.1\times
10^{-10}\rm{GeV}}{E\times 10^{-3}}\frac{n_{e}}{10^{-3}c\rm{m}^{3}}$ (6.5)
$\displaystyle\Delta^{\mathrm{QED}}$
$\displaystyle=4.1\times\frac{10^{-6}\rm{GeV}}{E\times
10^{-3}}\frac{B_{x}^{2}+B_{y}^{2}}{\mu\rm{G}^{2}}$
$\displaystyle\Delta^{\phi}_{i}$ $\displaystyle=\frac{-7.8\times
10^{-3}\rm{GeV}}{E\times
10^{-3}}\left(\frac{m_{i}}{10^{-8}\rm{GeV}}\right)^{2}$
$\displaystyle(\Delta^{\phi\gamma}_{x})_{i}$ $\displaystyle=7.6\times
10^{-2}\frac{g^{\gamma}_{i}\rm{GeV}}{5\times 10^{-11}}\frac{B_{y}}{\mu\rm{G}}$
$\displaystyle(\Delta^{a\gamma}_{y})_{i}$ $\displaystyle=7.6\times
10^{-2}\frac{g^{\gamma}_{i}\rm{GeV}}{5\times
10^{-11}}\frac{B_{x}}{\mu\rm{G}}\,,$

here, $E$ is the photon/ALP energy (in $\rm{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 $\mu\rm{G}$) [59]. The mass of the $i^{\rm{th}}$ ALP state,
denoted $m_{i}$, is taken to be distributed logarithmically within
$\left[10^{-8},10^{-5}\right]\rm{eV}$. The various components of $H$ can be
interpreted as follows: $\Delta^{\phi}_{i}$ are the ALP mass terms;
$\Delta^{pl}$ is a photon effective mass terms, induced by thermal effects in
the electron plasma; $(\Delta^{\phi\gamma}_{x})_{i}$ and
$(\Delta^{\phi\gamma}_{y})_{i}$ are the ALP-photon couplings for each photon
polarisation; and, $\Delta^{\rm{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_{\mathrm{el}}$. The state of the system after propagation through the
cluster is then given by:

$\displaystyle\rho_{\mathrm{out}}^{C}=$
$\displaystyle\tilde{G}^{C}\rho_{\gamma}\left(\tilde{G}^{C}\right)^{\dagger}\,,$
(6.6) $\displaystyle\tilde{G}^{C}=$ $\displaystyle\prod_{i}G_{i}^{C}\,,$

where $i$ runs over each domain in the cluster.

####  6.1.2 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 $\Delta
L^{IG}=50\,\mathrm{Mpc}$. The intergalactic medium has a luminosity redshift
($z$) dependent field strength given by:

$B^{IG}(z)=B_{0}^{IG}(1+z)^{2}\,,$ (6.7)

where $B_{0}^{IG}=1\,\mathrm{nG}$. The intergalactic electron density is
approximated with a constant value of
$n_{\mathrm{el}}^{IG}=10^{-7}\,\mathrm{cm}^{-3}$. To account for VHE photon
scattering off the EBL we introduce a non-unitary decay matrix ($D$) for each
domain:

$D(\tau)=\begin{bmatrix}\exp(-\tau/2)&0&0\\\ 0&\exp(-\tau/2)&0\\\
0&0&\ddots\end{bmatrix}\,,$ (6.8)

where $\tau$ is the optical depth associated with propagation through a given
domain. In general, $\tau$ depends on the photon energy and the redshift of
the domain. We use the EBL opacity model given in Ref. [60]. After computing
$D(\tau_{i})$ for each domain, the ALP-photon state is propagated as follows:

$\displaystyle\rho_{\mathrm{out}}^{IG}$
$\displaystyle=\tilde{G}^{IG}\rho_{\mathrm{out}}^{C}\left(\tilde{G}^{IG}\right)^{\dagger},\,$
(6.9) $\displaystyle\tilde{G}^{IG}$
$\displaystyle=\prod_{i}D(\tau_{i})G_{i}^{IG}\,,$

where $G_{i}^{IG}$ is calculated as before using Eq. (6.4) with $i$ running
over each domain in the IGM.

####  6.1.3 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. [61], 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 [61]. 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:

$r_{\mathrm{bound}}=r_{j}e^{(\theta-\pi)\cot\theta}\,,\quad\quad\theta=\frac{\pi}{180}(90-\alpha_{\mathrm{open}})\,,$
(6.10)

where $\alpha_{\mathrm{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\,\mathrm{cm}^{-3}$
to construct the MW evolution operator. The final state of the system is then
given by the product over domains $i$:

$\displaystyle\rho_{\mathrm{out}}^{MW}$
$\displaystyle=\tilde{G}^{MW}\rho_{\mathrm{out}}^{IG}\left(\tilde{G}^{MW}\right)^{\dagger}\,,$
(6.11) $\displaystyle\tilde{G}^{MW}$ $\displaystyle=\prod_{i}G_{i}^{MW}\,,$

where again we discretise our computation of the evolution operator ($G_{i}$)
- in this case, using an approximate-coherence length of $100\,\mathrm{pc}$.
We begin the Milky Way propagation stage when the ALP/photon reaches a radial
distance of $20\,\mathrm{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 ($\rho_{\mathrm{out}}^{MW}$) onto the two possible photon states:

$\displaystyle P_{\gamma\rightarrow\gamma}$
$\displaystyle=\operatorname{Tr}\left(2\rho_{\gamma}\mathcal{T}_{\mathrm{out}}^{MW}\right)$
(6.12)
$\displaystyle=\operatorname{Tr}\left(2\rho_{\gamma}\left[\tilde{G}^{MW}\tilde{G}^{IG}\tilde{G}^{C}\right]\rho_{\gamma}\left[\tilde{G}^{MW}\tilde{G}^{IG}\tilde{G}^{C}\right]^{\dagger}\right)\,.$

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: cross-section 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)\mathrm{kpc}$, is marked by the black circle.

###  6.2 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_{\gamma\rightarrow\gamma}$.
It is thus necessary to simulate with a large set of different field
directions. $P_{\gamma\rightarrow\gamma}$ 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^{\gamma}_{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_{\mathrm{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 [58]. We plan to perform a more
detailed follow-up study using a larger sample of sources and more recent
data.

$j$ | Source | Experiment | Fit Function
---|---|---|---
1 | Mkn 501 | HEGRA | LP
2 | 1ES0414+009 | H.E.S.S | PL
3 | 1ES0229-200 | H.E.S.S. | PL
4 | Mkn421 | H.E.S.S. | LP
5 | 1ES1101-232 | H.E.S.S | PL
6 | 1ES0347-121 | H.E.S.S | PL
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.

For a given choice of blazar source and grid point ($N,g^{\gamma}$) in model
parameter space, the result of our simulation is a set of $N_{\mathrm{samp}}$
survival probability spectra $P_{\gamma\rightarrow\gamma}(E)$, where $E$ is
the photon energy. To marginalise over these, we use the students p-statistic
($p_{t}$), following [58]; the calculation of $p_{t}$ is described in Appendix
C. For each grid point, we obtain a distribution of $N_{\mathrm{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_{\mathrm{samp}}$. This method facilitates comparison with Ref.
[58]. 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}\left(p_{95}\right)$
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. [58]; 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_{\gamma\rightarrow\gamma}$ in the
MW. This effect can be very significant, in this case, completely negating the
opacity decrease caused by reduced EBL scattering;
$P_{\gamma\rightarrow\gamma}$ is significantly lower in the $20$ ALP case than
it is for no ALPs (black line in the Fig. 5). 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.

Figure 7: A heat map showing $p_{95}$ computed using $6$ sources and
$N_{\mathrm{samp}}=1000$ samples per source. Lower values (higher $p_{95}$)
indicate better agreement between the model and the data.

##  7 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^{\gamma}}^{2}=\sum{g_{i}^{\gamma}}^{2}$ and
${g^{e}}^{2}=\sum{g_{i}^{e}}^{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 ALP-photon
coupling. Still, other search strategies such as ADMX [62] 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 ALP-photon 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.

## Acknowledgements

We thank Anthony Brown, Christopher Dessert, James Halverson, Sebastian Hoof
and David Marsh for valuable conversations. This work used the DiRAC@Durham
facility managed by the Institute for Computational Cosmology on behalf of the
STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS
capital funding via STFC capital grants ST/P002293/1, ST/R002371/1 and
ST/S002502/1, Durham University and STFC operations grant ST/R000832/1. DiRAC
is part of the National e-Infrastructure. This work was partly performed in
part at the Aspen Center for Physics, which the National Science Foundation
supports grant PHY-2210452. The authors are supported by STFC grant
ST/T001011/1. FCD is supported by EPSRC Stephen Hawking Fellowship
EP/T01668X/1.

## Data Access Statement

The code used to simulate VHE photon propagation from blazars is hosted in the
following repository:
https://gitlab.dur.scotgrid.ac.uk/James_Maxwell/blazarphotonaxion

##  Appendix A Appendix A: Statistical approach to producing mixing matrices
using the Haar measure

Here we can outline how the mixing matrices are inductively produced.
Following Section 3, one can generate a sample of $N-1$ mixing angles,
$\\{\theta_{ij}\\}$, in spherical polar coordinates. As a simplifying
assumption, we take the mixing matrices, $U\in SO(N)$. Since $SO(N)$ is
parametrised by $N(N-1)/2$ mixing angles, $\left\\{\theta_{ij}\right\\}_{1\leq
i\leq j\leq N-1}$, the mixing matrix for $SO(2)$ can be parametrised by a
single angle, $\theta_{11}$:

$U_{2}=\begin{bmatrix}\cos\theta_{11}&\sin\theta_{11}\\\
-\sin\theta_{11}&\cos\theta_{11}\end{bmatrix}\,.$ (A.1)

The Haar measure of this matrix is given by Eq. (3.2) and explicitly is

$dV=d\theta_{11}\,,$ (A.2)

which informs us we sample over $\theta_{11}$ uniformly in $[0,2\pi]$.
Generalising to higher dimensions, the structure of the matrices becomes much
less trivial. The most straightforward approach to writing down a higher
$SO(N)$ matrix is as a product of matrices:

$U_{N}=U_{N}^{\prime}S_{N}(U_{N-1})$ (A.3)

where $S_{N}(U_{N-1})$, given by:

$\displaystyle S_{N}(U_{N-1})=\begin{bmatrix}U_{N-1}&\bm{0}\\\
\bm{0}^{T}&1\end{bmatrix}\,,$ (A.4)

where $\bm{0}$ is a $(N-1)$ dimensional zero vector and $S_{N}(U_{N-1})$
represents rotations in the hyperplane orthogonal to the $N$th dimension.
$U_{N}^{\prime}$ is the $N$-dimensional analog of $U_{2}$, describing
rotations that involve the $N$th coordinate:

$\displaystyle{U_{N}^{\prime}}_{1,1}$
$\displaystyle=\cos\left(\theta_{1,N-1}\right)$ (A.5) $\displaystyle
U_{N1,i+1}^{\prime}$
$\displaystyle=\cos\left(\theta_{i+1,N-1}\right)\prod_{m=1}^{i}\sin\left(\theta_{m,N-1}\right)\quad\forall\quad
i\in[1,N-2]$ $\displaystyle{U_{N}^{\prime}}_{1,N}$
$\displaystyle=\prod_{m=1}^{N-1}\sin\left(\theta_{m,N-1}\right)$
$\displaystyle{U_{N}^{\prime}}_{i,j}$ $\displaystyle=0\quad\forall\quad
j<i-1\And i>1$ $\displaystyle{U_{N}^{\prime}}_{i,j}$
$\displaystyle=-\sin\left(\theta_{j,N-1}\right)\quad\forall\quad j=i-1\And
i>1$ $\displaystyle{U_{N}^{\prime}}_{i+1,i+1}$
$\displaystyle=\cos\left(\theta_{i,N-1}\right)\cos\left(\theta_{i+1,N-1}\right)\quad\forall\quad
i\in[1,N-2]$ $\displaystyle{U_{N}^{\prime}}_{N,N}$
$\displaystyle=\cos\left(\theta_{N-1,N-1}\right)$

Explicitly, the $SO(3)$ mixing matrix has the following form:

$U_{3}=\left[\begin{array}[]{ccc}\cos\theta_{12}&\sin\theta_{12}\cos\theta_{22}&\sin\theta_{12}\sin\theta_{22}\\\
-\sin\theta_{12}&\cos\theta_{12}\cos\theta_{22}&\cos\theta_{12}\sin\theta_{22}\\\
0&-\sin\theta_{22}&\cos\theta_{22}\end{array}\right]\left[\begin{array}[]{ccc}\cos\theta_{11}&\sin\theta_{11}&0\\\
-\sin\theta_{11}&\cos\theta_{11}&0\\\ 0&0&1\end{array}\right]\,.$ (A.6)

Hence we produce the symbolic form of the N-dimensional mixing matrices
inductively. We then sample the angles according to the Haar measure as shown
in Eq. (3.2), and we can visualise these samples as populating an
$N$-dimensional sphere uniformly.

##  Appendix 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^{\gamma}_{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
($\sigma^{\mathrm{max}}$) that is consistent with non-detection:

$\sigma_{\mathrm{max}}=\Phi_{\rm
tot}(g^{e}_{N=1},g^{\gamma}_{N=1})(g^{\gamma}_{N=1})^{2}\,,$ (B.1)

where $\Phi_{\rm tot}$ is the total single ALP flux at the detector; computed
using Eqs. (4.4), (4.5) and (4.6):

$\Phi_{\rm
tot}(g^{e}_{N=1},g^{\gamma}_{N=1})=\Phi_{P}(g^{\gamma}_{N=1})+\Phi_{B}(g^{e}_{N=1})+\Phi_{C}(g^{e}_{N=1})\,.$
(B.2)

For any point in the coupling parameter space, the total flux at the detector
must, therefore, be less than:

$\Phi_{\mathrm{max}}(g^{\gamma})=\sigma_{\mathrm{max}}/{g^{\gamma}}^{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 $\omega=1\,\rm{keV}$ (corresponding to the detector threshold) and
an upper bound of $\omega=10\,{\rm keV}$.

The recast bounds in Fig. 2 are shown as heat maps. The weight assigned to
each set of couplings $(g^{e},g^{\gamma})$, corresponds to the proportion of
our sample of survival and conversion probabilities,
$\\{P_{\gamma\rightarrow\gamma}\\}$ and $\\{P_{e\rightarrow\gamma}\\}$, that
are consistent with Eq. (B.3). That is to say that, for a given
$(g^{e},g^{\gamma})$ and probability
$P_{e\rightarrow\gamma}^{l}\in\\{P_{e\rightarrow\gamma}\\}$, we find the
proportion ($\eta$) of $\\{P_{\gamma\rightarrow\gamma}\\}$ that satisfy:

$\Phi_{\mathrm{max}}(g^{\gamma})\geq\Phi_{P}(g^{\gamma})P_{\gamma\rightarrow\gamma}+(\Phi_{B}(g^{e})+\Phi_{C}(g^{e}))P_{e\rightarrow\gamma}$
(B.4)

At this point, $\eta$ could be found by Monte Carlo sampling
$\\{P_{\gamma\rightarrow\gamma}\\}$ 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_{\mathrm{Crit}}$) of $P_{\gamma\gamma}$ that achieves
equality in Eq. (B.4):

$P_{\mathrm{Crit}}(g^{e},g^{\gamma},P_{e\gamma}^{l})=\frac{1}{\Phi_{P}(g^{\gamma})}\left[\Phi_{\mathrm{max}}(g^{\gamma})-(\Phi_{B}(g^{e})+\Phi_{C}(g^{e}))P_{e\rightarrow\gamma}^{l}\right]\,,$
(B.5)

$\eta$ is now given by the cumulative sum – evaluated at $P_{\mathrm{Crit}}$ –
of the histogram of $\\{P_{\gamma\rightarrow\gamma}\\}$;
$\eta=C(P_{\mathrm{Crit}})$. Repeating this process for every
$P_{e\rightarrow\gamma}^{l}\in\\{P_{e\rightarrow\gamma}\\}$, the total weight
associated with a coupling pair is given by:

$W(g^{e},g^{\gamma})=\frac{1}{N_{ae}}\sum_{l}C(P_{\mathrm{Crit}}(g^{e},g^{\gamma},P_{e\rightarrow\gamma}^{l}))\,,$
(B.6)

where $N_{ae}$ is the size of $\\{P_{e\rightarrow\gamma}\\}$ – taken here to
be $N_{ae}=N_{a\gamma}=10^{4}$. $W(g^{e},g^{\gamma})$ is equivalent to the
fraction of $\\{P_{\gamma\rightarrow\gamma}\\}$ and
$\\{P_{e\rightarrow\gamma}\\}$ that satisfy Eq. (B.4) and will saturate to $1$
for sufficiently low couplings.

##  Appendix 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. [58, 63]. In particular, we use the
observed spectra, $\Phi_{i}^{\text{obs }}$, shown in Table 1.

We note that $\Phi_{i}^{\text{obs }}$ 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:

$\left\langle
P_{\gamma\rightarrow\gamma}\right\rangle_{i}(g^{\gamma},g^{e})=\frac{1}{\Delta
E_{i}}\int_{\Delta E_{i}}\mathrm{\leavevmode\nobreak\
d}EP_{\gamma\rightarrow\gamma}(E)\,,$ (C.1)

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

$\Phi_{i}(g^{\gamma},g^{e})=\left\langle
P_{\gamma\rightarrow\gamma}\right\rangle_{i}^{-1}(g^{\gamma},g^{e})\Phi_{i}^{\mathrm{obs}}\,,$
(C.2)

where again, we have kept the theory parameter dependence explicit. For a
given point in the model parameter space, $(g^{\gamma},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 $\Phi_{i}(g^{\gamma},g^{e})$ to a power law:

$f(E)=\begin{cases}N_{0}\left(E/E_{0}\right)^{-\Gamma},&p_{\text{fit
}}\geqslant 0.05\\\
N_{0}\left(E/E_{0}\right)^{-\left(\Gamma+\beta_{c}\ln\left(E/E_{0}\right)\right)},&\text{
otherwise }\end{cases}$ (C.3)

where $p_{\text{fit }}$ denotes the fit probability. The first power law
contains three parameters: $N_{0}$, $E_{0}$, and $\Gamma$ while in the second
(which is a parabola in log-log space), there are four fit parameters,
$N_{0}$, $E_{0}$, $\beta_{c}$ and $\Gamma$. The power law function is used
unless the baseline (no ALP) fit probability (as calculated by Ref. [58]) 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:

$\chi_{i}=\frac{\Phi_{i}-f\left(E_{i}\right)}{\sigma_{i}}\,,$ (C.4)

where ${\sigma_{i}}$ denotes the statistical measurement uncertainty at 68$\%$
confidence level. The value of $\chi^{2}$ can now easily be obtained by
summing the residuals over all energy bins and squaring.

Assuming that $P_{\gamma\rightarrow\gamma}$ accurately predicts the Universe’s
opacity to VHE $\gamma$ 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:

$t=\frac{\bar{\chi}}{\sqrt{\sigma_{\chi}/N_{\chi}}}\,.$ (C.5)

where $\bar{\chi}$ denotes the mean and $\sigma_{\chi}$ the variance of the
residual distribution, which comprises a total of $N_{\chi}$ 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.

## References

  * [1] D. J. E. Marsh, Axion Cosmology, Phys. Rept. 643 (2016) 1–79, [arXiv:1510.07633].
  * [2] F. Chadha-Day, J. Ellis, and D. J. E. Marsh, Axion dark matter: What is it and why now?, Sci. Adv. 8 (2022), no. 8 abj3618, [arXiv:2105.01406].
  * [3] R. D. Peccei and H. R. Quinn, CP Conservation in the Presence of Instantons, Phys. Rev. Lett. 38 (1977) 1440–1443.
  * [4] F. Wilczek, Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons, Phys. Rev. Lett. 40 (1978) 279–282.
  * [5] M. B. Voloshin, Channel coupling in e+ e- annihilation into heavy meson pairs at the D* anti-D* threshold, hep-ph/0602233.
  * [6] P. Svrcek and E. Witten, Axions In String Theory, JHEP 06 (2006) 051, [hep-th/0605206].
  * [7] M. Cicoli, M. Goodsell, and A. Ringwald, The type IIB string axiverse and its low-energy phenomenology, JHEP 10 (2012) 146, [arXiv:1206.0819].
  * [8] I. Broeckel, M. Cicoli, A. Maharana, K. Singh, and K. Sinha, Moduli stabilisation and the statistics of axion physics in the landscape, JHEP 08 (2021) 059, [arXiv:2105.02889]. [Addendum: JHEP 01, 191 (2022)].
  * [9] M. Demirtas, N. Gendler, C. Long, L. McAllister, and J. Moritz, PQ Axiverse, arXiv:2112.04503.
  * [10] N. Gendler, D. J. E. Marsh, L. McAllister, and J. Moritz, Glimmers from the Axiverse, arXiv:2309.13145.
  * [11] M. J. Stott and D. J. E. Marsh, Black hole spin constraints on the mass spectrum and number of axionlike fields, Phys. Rev. D 98 (2018), no. 8 083006, [arXiv:1805.02016].
  * [12] V. M. Mehta, M. Demirtas, C. Long, D. J. E. Marsh, L. McAllister, and M. J. Stott, Superradiance in string theory, JCAP 07 (2021) 033, [arXiv:2103.06812].
  * [13] M. J. Stott, D. J. E. Marsh, C. Pongkitivanichkul, L. C. Price, and B. S. Acharya, Spectrum of the axion dark sector, Phys. Rev. D 96 (2017), no. 8 083510, [arXiv:1706.03236].
  * [14] M. Reig, The stochastic axiverse, JHEP 09 (2021) 207, [arXiv:2104.09923].
  * [15] J. E. Kim, H. P. Nilles, and M. Peloso, Completing natural inflation, JCAP 01 (2005) 005, [hep-ph/0409138].
  * [16] S. Dimopoulos, S. Kachru, J. McGreevy, and J. G. Wacker, N-flation, JCAP 08 (2008) 003, [hep-th/0507205].
  * [17] P. Agrawal, J. Fan, M. Reece, and L.-T. Wang, Experimental Targets for Photon Couplings of the QCD Axion, JHEP 02 (2018) 006, [arXiv:1709.06085].
  * [18] N. Kitajima and F. Takahashi, Resonant conversions of QCD axions into hidden axions and suppressed isocurvature perturbations, JCAP 01 (2015) 032, [arXiv:1411.2011].
  * [19] I. Obata, Implications of the cosmic birefringence measurement for the axion dark matter search, JCAP 09 (2022) 062, [arXiv:2108.02150].
  * [20] T. Higaki, N. Kitajima, and F. Takahashi, Hidden axion dark matter decaying through mixing with QCD axion and the 3.5 keV X-ray line, JCAP 12 (2014) 004, [arXiv:1408.3936].
  * [21] N. Glennon, N. Musoke, and C. Prescod-Weinstein, Simulations of multifield ultralight axionlike dark matter, Phys. Rev. D 107 (2023), no. 6 063520, [arXiv:2302.04302].
  * [22] B. Gavela, P. Quílez, and M. Ramos, Multiple QCD axion, arXiv:2305.15465.
  * [23] F. Chadha-Day, Axion-like particle oscillations, JCAP 01 (2022), no. 01 013, [arXiv:2107.12813].
  * [24] Z. Chen, A. Kobakhidze, C. A. J. O’Hare, Z. S. C. Picker, and G. Pierobon, Phenomenology of the companion-axion model: photon couplings, Eur. Phys. J. C 82 (2022), no. 10 940, [arXiv:2109.12920].
  * [25] S. Gasparotto and E. I. Sfakianakis, Cosmic Birefringence from the Axiverse, arXiv:2306.16355.
  * [26] K. R. Dienes, E. Dudas, and T. Gherghetta, Invisible axions and large radius compactifications, Phys. Rev. D 62 (2000) 105023, [hep-ph/9912455].
  * [27] A. Ayala, I. Domínguez, M. Giannotti, A. Mirizzi, and O. Straniero, Revisiting the bound on axion-photon coupling from Globular Clusters, Phys. Rev. Lett. 113 (2014), no. 19 191302, [arXiv:1406.6053].
  * [28] G. G. Raffelt, Axion Constraints From White Dwarf Cooling Times, Phys. Lett. B 166 (1986) 402–406.
  * [29] S. I. Blinnikov and N. V. Dunina-Barkovskaya, The cooling of hot white dwarfs: a theory with non-standard weak interactions, and a comparison with observations, Mon. Not. Roy. Astron. Soc. 266 (1994) 289–304.
  * [30] J. Halverson and F. Ruehle, Computational Complexity of Vacua and Near-Vacua in Field and String Theory, Phys. Rev. D 99 (2019), no. 4 046015, [arXiv:1809.08279].
  * [31] J. Halverson, C. Long, B. Nelson, and G. Salinas, Towards string theory expectations for photon couplings to axionlike particles, Phys. Rev. D 100 (2019), no. 10 106010, [arXiv:1909.05257].
  * [32] L. J. Hall, H. Murayama, and N. Weiner, Neutrino mass anarchy, Phys. Rev. Lett. 84 (2000) 2572–2575, [hep-ph/9911341].
  * [33] N. Haba and H. Murayama, Anarchy and hierarchy, Phys. Rev. D 63 (2001) 053010, [hep-ph/0009174].
  * [34] A. de Gouvea and H. Murayama, Statistical test of anarchy, Phys. Lett. B 573 (2003) 94–100, [hep-ph/0301050].
  * [35] A. de Gouvea and H. Murayama, Neutrino Mixing Anarchy: Alive and Kicking, Phys. Lett. B 747 (2015) 479–483, [arXiv:1204.1249].
  * [36] G. G. Raffelt, ASTROPHYSICAL AXION BOUNDS DIMINISHED BY SCREENING EFFECTS, Phys. Rev. D 33 (1986) 897.
  * [37] L. M. Krauss, J. E. Moody, and F. Wilczek, A STELLAR ENERGY LOSS MECHANISM INVOLVING AXIONS, Phys. Lett. B 144 (1984) 391–394.
  * [38] S. Dimopoulos, G. D. Starkman, and B. W. Lynn, Atomic Enhancements in the Detection of Weakly Interacting Particles, Phys. Lett. B 168 (1986) 145–150.
  * [39] S. Dimopoulos, J. A. Frieman, B. W. Lynn, and G. D. Starkman, Axiorecombination: A New Mechanism for Stellar Axion Production, Phys. Lett. B 179 (1986) 223–227.
  * [40] M. Pospelov, A. Ritz, and M. B. Voloshin, Bosonic super-WIMPs as keV-scale dark matter, Phys. Rev. D 78 (2008) 115012, [arXiv:0807.3279].
  * [41] K. O. Mikaelian, Astrophysical Implications of New Light Higgs Bosons, Phys. Rev. D 18 (1978) 3605.
  * [42] M. Fukugita, S. Watamura, and M. Yoshimura, Light Pseudoscalar Particle and Stellar Energy Loss, Phys. Rev. Lett. 48 (1982) 1522.
  * [43] M. Fukugita, S. Watamura, and M. Yoshimura, Astrophysical Constraints on a New Light Axion and Other Weakly Interacting Particles, Phys. Rev. D 26 (1982) 1840.
  * [44] J. Redondo, Solar axion flux from the axion-electron coupling, JCAP 12 (2013) 008, [arXiv:1310.0823].
  * [45] I. G. Irastorza et al., Towards a new generation axion helioscope, JCAP 06 (2011) 013, [arXiv:1103.5334].
  * [46] K. Barth et al., CAST constraints on the axion-electron coupling, JCAP 05 (2013) 010, [arXiv:1302.6283].
  * [47] P. Gondolo and G. G. Raffelt, Solar neutrino limit on axions and keV-mass bosons, Phys. Rev. D 79 (2009) 107301, [arXiv:0807.2926].
  * [48] C. Dessert, A. J. Long, and B. R. Safdi, X-ray Signatures of Axion Conversion in Magnetic White Dwarf Stars, Phys. Rev. Lett. 123 (2019), no. 6 061104, [arXiv:1903.05088].
  * [49] A. Harayama, Y. Terada, M. Ishida, T. Hayashi, A. Bamba, and M. S. Tashiro, Search for non-thermal emissions from an isolated magnetic white dwarf, euve j0317 855, with suzaku, Publications of the Astronomical Society of Japan 65 (2013), no. 4 73.
  * [50] B. Kulebi, S. Jordan, E. Nelan, U. Bastian, and M. Altmann, Constraints on the origin of the massive, hot, and rapidly rotating magnetic white dwarf RE J 0317-853 from an HST parallax measurement, Astron. Astrophys. 524 (2010) A36, [arXiv:1007.4978].
  * [51] G. I. Rubtsov and S. V. Troitsky, Breaks in gamma-ray spectra of distant blazars and transparency of the Universe, Soviet Journal of Experimental and Theoretical Physics Letters 100 (Nov., 2014) 355–359, [arXiv:1406.0239].
  * [52] K. Kohri and H. Kodama, Axion-like particles and recent observations of the cosmic infrared background radiation, Physical Review D 96 (sep, 2017).
  * [53] J. Biteau and D. A. Williams, THE EXTRAGALACTIC BACKGROUND LIGHT, THE HUBBLE CONSTANT, AND ANOMALIES: CONCLUSIONS FROM 20 YEARS OF TeV GAMMA-RAY OBSERVATIONS, The Astrophysical Journal 812 (oct, 2015) 60.
  * [54] A. Domínguez and M. Ajello, Spectral analysis of Fermi-LAT blazars above 50 GeV, Astrophys. J. Lett. 813 (2015), no. 2 L34, [arXiv:1510.07913].
  * [55] A. De Angelis, G. Galanti, and M. Roncadelli, Relevance of axion-like particles for very-high-energy astrophysics, Phys. Rev. D 84 (2011) 105030, [arXiv:1106.1132]. [Erratum: Phys.Rev.D 87, 109903 (2013)].
  * [56] M. Simet, D. Hooper, and P. D. Serpico, Milky way as a kiloparsec-scale axionscope, Physical Review D 77 (mar, 2008).
  * [57] M. A. Sanchez-Conde, D. Paneque, E. Bloom, F. Prada, and A. Dominguez, Hints of the existence of Axion-Like-Particles from the gamma-ray spectra of cosmological sources, Phys. Rev. D 79 (2009) 123511, [arXiv:0905.3270].
  * [58] M. Meyer, D. Horns, and M. Raue, First lower limits on the photon-axion-like particle coupling from very high energy gamma-ray observations, Physical Review D 87 (feb, 2013).
  * [59] E. Masaki, A. Aoki, and J. Soda, Photon-axion conversion, magnetic field configuration, and polarization of photons, Physical Review D 96 (aug, 2017).
  * [60] A. Franceschini, G. Rodighiero, and M. Vaccari, Extragalactic optical-infrared background radiation, its time evolution and the cosmic photon-photon opacity, Astronomy & Astrophysics 487 (jun, 2008) 837–852.
  * [61] R. Jansson and G. R. Farrar, A NEW MODEL OF THE GALACTIC MAGNETIC FIELD, The Astrophysical Journal 757 (aug, 2012) 14\.
  * [62] ADMX Collaboration, C. Boutan et al., Piezoelectrically Tuned Multimode Cavity Search for Axion Dark Matter, Phys. Rev. Lett. 121 (2018), no. 26 261302, [arXiv:1901.00920].
  * [63] M. Meyer, The Opacity of the Universe for High and Very High Energy $\gamma$-Rays.  PhD thesis, Hamburg U., 2013.