# Robust bounds on ALP dark matter from dwarf spheroidal galaxies
in the optical MUSE-Faint survey

Elisa Todarello  Marco Regis  Javier Reynoso-Cordova  Marco Taoso  Daniel Vaz
Jarle Brinchmann  Matthias Steinmetz  Sebastiaan L. Zoutendijk

###### Abstract

Nearby dwarf spheroidal galaxies are ideal targets in the search for indirect
dark matter (DM) signals. In this work, we analyze MUSE spectroscopic
observations of a sample of five galaxies, composed of both classical and
ultra-faint dwarf spheroidals. The goal is to search for radiative decays of
axion-like particles (ALPs) in the mass range of 2.7-5.3 eV. After taking into
account the uncertainties associated with the DM spatial distribution in the
galaxies, we derive robust bounds on the effective ALP-two-photon coupling.
They lie well below the QCD axion band and are significantly more constraining
than limits from other probes, in the relevant mass range. We also test the
possible presence of a positive signal, concluding that none of the wavelength
channels selected for this analysis, i.e., not affected by large background
contamination, is exhibiting such evidence.

##  1  Introduction

Axion-like particles (ALPs) are pseudo Nambu-Goldstone bosons that arise in
extensions of the Standard Model and that can act as cold dark matter (DM)
candidates [1, 2, 3]. A particularly well-motivated example is the QCD axion,
which is associated with the Peccei-Quinn symmetry solution to the strong CP
problem [4, 5, 6, 7]. Generically, ALPs couple to photons through the operator
$\mathcal{L}=-\frac{1}{4}g_{a\gamma}\,a\,F_{\mu\nu}\tilde{F}_{\mu\nu},$ where
$a$ is the ALP field, $F_{\mu\nu}$ is the electromagnetic field strength,
$\tilde{F}_{\mu\nu}$ its dual, and $g_{a\gamma}$ the coupling constant. This
interaction term leads to a variety of possibilities to detect ALPs in
laboratory experiments or with astrophysical and cosmological probes, see e.g.
[8, 9, 10] for reviews. In astrophysical environments, an almost monochromatic
photon emission is produced by the radiative decay of ALP DM, and, for ALP
masses in the eV range, this photon line falls in the optical and near-
infrared bands. In this frequency range, several upper bounds on this signal
have been derived from observations [11, 12, 13, 14, 15, 16, 17, 18]. In
particular, ref. [14] derived the currently most stringent constraints on ALP
radiative decays for masses between 2.7 and 5.3 eV, improving previous bounds
by more than an order of magnitude. Interestingly, in recent years, ALPs
masses in the eV mass range have been invoked to explain excesses in the
measured cosmic near-infrared background and its angular anisotropies [19, 20,
21, 22, 15, 16]. The upper limits of [14] severely challenge some of these
scenarios [21].

The analysis in [14] is based on spectroscopic observations of the Leo T dwarf
spheroidal galaxy obtained with the Multi Unit Spectroscopic Explorer (MUSE)
at the Very Large Telescope (VLT) [23]. Dwarf spheroidal galaxies are ideal
targets for searching for DM decay signals because they contain large DM
densities and are relatively close to us. Moreover, measurements of the line-
of-sight velocity of the stars in these objects allow us to infer the
underlying DM distribution along with its uncertainty. This information is
instrumental in order to reliably predicting the ALP signal, and deriving
robust bounds on the ALP decay lifetime.

With the present work, we extend and improve the analysis in [14] in two ways:
we enlarge the dataset exploiting recent MUSE observations of other dwarf
galaxies, and we implement a more detailed treatment of the DM distribution
and its uncertainty, taking advantage of recent analyses. More specifically,
in addition to Leo T, we consider MUSE observations of the dwarf spheroidal
galaxies Sculptor, Eridanus 2, Grus 1, and Hydra II. Then, for Eridanus 2,
Grus 1, Hydra II, and Leo T, we make use of the recent determination of the DM
content in these objects performed by the MUSE collaboration [24]. Concretely,
we consider two parametrizations for the DM density, namely the
Navarro–Frenk–White (NFW) model [25] and a cored profile. We account for the
uncertainty on these DM distributions including the corresponding likelihoods
derived in [24] in our statistical analysis; see Sec. 4 for details. For
Sculptor, we follow the same procedure but we derive the DM distribution and
the relevant likelihood ourselves. This is accomplished by means of a Jeans
analysis, using data from [26, 27] and employing the same method as [24].

For each target, we perform a search of ALP decay signals in the MUSE data,
and then we combine the individual bounds in a global analysis. We find upper
limits on the ALP lifetime similar to, but slightly weaker than, those in
[14]. Finally, excluding channels severely contaminated by background, we do
not find significant evidence for an ALP signal.

The structure of this paper is as follows. The data from MUSE observations are
presented in Sec. 2. The calculation of the ALP decay signal is discussed in
Sec. 3. In Sec. 4, we discuss the statistical analysis and results. We
conclude in Sec. 5. The Jeans analysis for Sculptor is discussed in Appendix
A. The caveats of our analysis assuming a cored dark matter profile are
addressed in Appendix B.

##  2  Observations and data reduction

As part of MUSE-Faint, a GTO survey of faint dwarf galaxies (PI Brinchmann),
Leo T, Sculptor, Eridanus 2, Grus 1, and Hydra II111Sculptor, which contains a
significantly larger number of stars than the other targets, is commonly
considered to be a “classical” dwarf, while Leo T, Eridanus 2, Grus 1, and
Hydra II are classified as “ultra-faint” dwarf galaxies. were observed with
MUSE, a large-field medium-resolution Integral Field Spectrograph installed on
the VLT. We use multiple exposures of 900 s of each galaxy, with a total
exposure time of 3.75 hours on Leo T (one field), 3 hours on Sculptor (one
field), 21.5 hours on Eridanus 2 (five fields), 4 hours on Grus 1 (one field),
and 14.75 hours on Hydra II (four fields).

The data were taken in the Wide Field Mode with adaptive optics (WFM-AO),
which provides a $1\times 1\leavevmode\nobreak\ \mathrm{arcmin^{2}}$ field of
view with a spatial sampling of $0.2\leavevmode\nobreak\
\mathrm{arcsec\leavevmode\nobreak\ pixel^{-1}}$. The data cover a wavelength
range of $4700-9350\leavevmode\nobreak\ \rm{\AA}$, sampled at a resolution of
$1.25\leavevmode\nobreak\ \rm{\AA}$. A blocking filter was used to remove the
light from the sodium laser of the adaptive optics system to avoid
contamination. This filter blocked light in the $5820-5970\leavevmode\nobreak\
\rm{\AA}$ ($2.13-2.08$ eV) range, which appears as a gap in the constraints
presented in the following.

For data reduction, we refer the reader to Ref. [24] for details, while here
we provide a brief summary. We performed the standard data reduction procedure
using the MUSE Data Reduction Software (DRS; version 2.8 [28]). Flux
calibration was carried out using flux standards observed during the night,
while atmospheric emission lines were removed by accounting for Raman
scattering caused by the laser light of the adaptive optics system. We
subtracted emission lines from the night sky that have well-known wavelengths
and result in increased noise at those wavelengths. We measured a spatial
resolution (full-width half maximum) of 0.61, 0.50, 0.53, 0.67, and 0.40
arcsec for Leo T, Sculptor, Eridanus 2, Grus 1, and Hydra II, respectively, at
a wavelength of $7000\leavevmode\nobreak\ \rm{\AA}$ in the reduced datacubes.

To ensure an accurate analysis of the data cubes, it is crucial to have a
reliable estimate of the noise. Previous studies (e.g., Ref [29]) have shown
that the MUSE Data Reduction Software (DRS) underestimates uncertainties in
the final data cube. Therefore, we proceeded as in Ref. [14] and re-estimated
the pixel-to-pixel variance directly from each individual exposure data cube
using the method described in Ref. [29], creating mask images using SExtractor
[30]. We then combined all single-exposure data cubes using MPDAF [31], to
create two final data cubes that were used in the subsequent analysis, one for
the r.m.s. noise $\sigma_{rms}$ and one for the observed flux density
$S_{obs}$.

The data contain numerous stellar sources within the field of view, both from
the dwarf galaxy and also from some likely foreground stars from the Milky
Way, as well as some galaxies. To minimize the impact of these sources on the
final results, we have identified and masked the brightest ones. This was
achieved by following the same approach as in Ref. [14] and involved two
steps. First, we generated a white-light image by summing over the wavelength
axis in each datacube. Next, we ran SExtractor on this white-light image, with
a detection threshold of $3\sigma$, resulting in a segmentation map that was
used to mask sources. Therefore, we consider only pixels where no sources are
detected in the white-light image. In Fig. 1 we show, for all the dwarf
galaxies under consideration, the flux density per beam solid angle averaged
over all the unmasked pixels.

Figure 1: Flux density per beam solid angle averaged over the unmasked pixels
of the map as a function of the wavelength of observation.

##  3  ALP signal

We model the DM halo in a dwarf galaxy as a spherical system. The possibility
of a non-spherical halo in Leo T and Sculptor has been investigated, e.g., in
[32]. They found negligible impact in the derived D-factor uncertainties. We
are not aware of non-spherical analysis for Eridanus 2, Grus 1, and Hydra II.
For simplicity, and given that we generically expect lower triaxiality for
lower-mass DM halos [33], we assume spherical symmetry for all the targets.

The flux density at wavelength $\lambda$ produced by decays of ALPs from a
given direction identified by $\theta$ can be computed as:

$S_{\lambda}(\theta)=\frac{\Gamma_{a}}{4\pi}\,\frac{1}{\sqrt{2\pi}\sigma_{\lambda}}\exp{\left[-\frac{(\lambda-\lambda_{obs})^{2}}{2\sigma_{\lambda}^{2}}\right]}\int
d\Omega\,d\ell\rho_{a}[r(\theta,\Omega,\ell)]\,B(\Omega)\;.$ (3.1)

The decay rate $\Gamma_{a}$ depends on the ALP mass $m_{a}$ and the effective
ALP-two-photon coupling $g_{a\gamma}$. In natural units, it reads
$\Gamma_{a}=g_{a\gamma}^{2}\,m_{a}^{3}/(64\pi)$.

The wavelength of emission can be computed as $\lambda_{em}=c/\nu_{em}$ with
$\nu_{em}=m_{a}/(4\pi)$. In Eq. 3.1, we neglect the velocity dispersion of
ALPs in the dwarf halo, since it is typically $\lesssim
10^{-4}\leavevmode\nobreak\ c$, namely smaller than the spectral resolution of
MUSE, $\sigma_{\lambda}/\lambda$ which varies between $1.2\times 10^{-4}$ and
$2.7\times 10^{-4}$. However, we cannot neglect the heliocentric radial
velocity of dwarfs, except for Leo T, and we correct for the Doppler shift
when deriving the results below. The observed wavelength $\lambda_{obs}$ is
then given by $\lambda_{obs}=\lambda_{em}(1+v_{radial})$.

We assume a Gaussian behavior for both the energy and the angular responses of
the detector, with FWHM as a function of the wavelength taken from Ref. [29]
and normalized to the value at 7000 Å mentioned in Sec. 2. The angular beam is
denoted by $B(\Omega)$.

Under the assumption of spherical symmetry, the DM density $\rho_{a}(r)$ is a
function only of the radial distance $r$ from the center of the dwarf, which
can be expressed in terms of the coordinate along the line of sight $\ell$ and
the angle of observation. Our reference scenario for the description of the DM
density profile is given by a “cuspy” distribution given by the NFW functional
form [25]:

$\rho_{\rm{NFW}}(r)=\frac{\rho_{s}}{\left(\frac{r}{r_{s}}\right)\left(1+\frac{r}{r_{s}}\right)^{2}}\;,$
(3.2)

where $\rho_{s}$ and $r_{s}$ are respectively the scale density and radius. We
also consider a second parameterization, dubbed “coreNFWtides”, which modifies
the cusp by allowing for a central core (e.g., due to star-formation feedback)
and includes a decrease in density beyond a tidal radius [34]. This
parametrization is not completely independent from the NFW one, since it
builds on it, but it has the advantage to allow an easy assessment of the
effect a core and a tidal radius have on our bounds. The coreNFWtides profile
is described in more detail in the Appendix A.

We constrain the DM profiles of the dwarf galaxies through a Jeans analysis of
the velocity dispersion of their stars. For Eridanus 2, Grus 1, Hydra II, and
Leo T, we used the likelihood derived in [24] (GravSphere method). In the case
of Sculptor, we derived the likelihood ourselves, using the same approach as
in [24], but with data from [26, 27]. The results for the Sculptor case are
reported in Appendix A. In the rest of our analysis, $\rho_{s}$ and $r_{s}$
are treated as free parameters of the model. For the coreNFWtides profile, the
additional parameters describing the DM distribution, see Appendix A, are
fixed to their global best-fit values from the just mentioned analyses of
dispersion velocities, in order to reduce the number of free parameters. We
expect that a full scan of all the parameters of the coreNFWtides distribution
would lead to a mild weakening of our bounds, as discussed in Appendix B.

The other parameters that enter Eq. 3.1 and that will be sampled in our scans
are $g_{a\gamma}$ and $m_{a}$. In total, there are four free parameters
describing the expected flux from ALPs.

##  4  Methods and Results

In our statistical analysis, we consider two types of data. On one side, we
have the dispersion velocities of the stellar component in the dwarf galaxy,
which allow us to infer the DM spatial distribution via Jeans analysis. From
the likelihood defined in [24], we derive a profile likelihood
$\mathcal{L}_{Jeans}^{j}$ depending only on $\rho_{s}$ and $r_{s}$, where all
the other “nuisance” parameters are profiled out (and, in the case of the
coreNFWtides profile, the additional parameters are set to their global best-
fit values). The index $j$ stands for the dwarf considered.

The second type of dataset we consider is the diffuse emission probed by MUSE
observations in the direction of the dwarf galaxies. As done in [14], we
compare the expected ALP signal with the observed data in each dwarf by means
of a Gaussian likelihood (omitting the index $j$ for simplicity):

$\mathcal{L}_{diff}=e^{-\chi^{2}/2}\;\;\;{\rm
with}\;\;\;\chi^{2}=\frac{1}{N_{pix}^{FWHM}}\sum_{i=1}^{N_{pix}}\left(\frac{S_{th}^{i}-S_{obs}^{i}}{\sigma_{rms}^{i}}\right)^{2}\;,$
(4.1)

where $S_{th}^{i}$ is the theoretical estimate for the flux density in the
pixel $i$, $S_{obs}^{i}$ is the observed flux density and $\sigma_{rms}^{i}$
is the r.m.s. error, both described in Sec. 2. The theoretical estimate is
given by Eq. 3.1 along with an additional spatially flat term
$S_{\lambda,flat}$ that we incorporate in the fit to each individual map at
every wavelength to account for incomplete sky subtraction. We consider this
flat term a nuisance parameter. $N_{pix}$ is the total number of pixels in the
area under investigation, which we chose to be a circle of $60^{\prime\prime}$
of radius. The number of pixels within the MUSE angular beam is
$N_{pix}^{FWHM}$, which has a size given by the aforementioned FWHM.

We define a likelihood $\mathcal{L}_{diff}$ which depends only on ALP
parameters by profiling out $S_{\lambda,flat}^{j}$ from the likelihood in Eq.
4.1. Then, assuming the two types of datasets to be independent, we can
define, at any given mass $m_{a}$, a global likelihood for each dwarf $j$:

$\mathcal{L}^{j}(g_{a\gamma},\rho_{s}^{j},r_{s}^{j})=\mathcal{L}_{diff}^{j}(g_{a\gamma},\rho_{s}^{j},r_{s}^{j})\times\mathcal{L}^{j}_{Jeans}(\rho_{s}^{j},r_{s}^{j})\;,$
(4.2)

and a combined likelihood considering all five targets simultaneously:

$\mathcal{L}^{all}(g_{a\gamma},\vec{\rho}_{s},\vec{r}_{s})=\prod_{j=1}^{5}\mathcal{L}^{j}(g_{a\gamma},\rho_{s}^{j},r_{s}^{j})$
(4.3)

To compute $\mathcal{L}_{diff}^{j}$, we scan a three-dimensional
logarithmically spaced grid of values of the following parameters:
$k=g_{a\gamma}\sqrt{\rho_{s}}$, $r_{s}$ and $S_{\lambda,flat}$. $k$ is sampled
in the range [$10^{-15}$, $10^{-10}$] GeV${}^{-1}(10^{8}\leavevmode\nobreak\
M_{\odot}\leavevmode\nobreak\ \mathrm{kpc}^{-3})^{1/2}$. The values of
$g_{a\gamma}$ and $\rho_{s}$ are then separated by sampling $g_{a\gamma}$ in
the range [$10^{-13.6}$, $10^{-10.5}$] GeV${}^{-1},$ while the ranges of
variation of $\rho_{s}$ and $r_{s}$ are different for every galaxy and are
determined from the 95% C.L. contour of $\mathcal{L}^{j}_{Jeans}$. Lastly, for
each channel, we first compute the best-fit value of $S_{\lambda,flat}$ for a
vanishing $g_{a\gamma}$, and then we scan $S_{\lambda,flat}$ in a range from
1/10 to 10 times this best-fit value.

We assume that
$\lambda_{c}(g_{a\gamma})=-2\ln[\mathcal{L}(g_{a\gamma},\vec{\rho}_{s}^{\,lbf},\vec{r}_{s}^{\,lbf})/\mathcal{L}(g_{a\gamma}^{b.f.},\vec{\rho}_{s}^{\,gbf},\vec{r}_{s}^{\,gbf})]$
follows a $\chi^{2}$-distribution with one d.o.f. and with one-sided
probability given by
$P=\int^{\infty}_{\sqrt{\lambda_{c}}}d\chi\,e^{-\chi^{2}/2}/\sqrt{2\,\pi}$,
where $g_{a\gamma}^{b.f.}$ denotes the best-fit value for the coupling at a
specific ALP mass. The superscript $gbf$ indicates the global best-fit, i.e.
for $g_{a\gamma}=g_{a\gamma}^{b.f.}$, whilst $lbf$ denotes the best-fit of
$\vec{\rho}_{s}$ and $\vec{r}_{s}$ for that given $g_{a\gamma}$. For the
analysis on a single dwarf $j$, one has just to replace
$(\vec{\rho}_{s},\vec{r}_{s})$ with $(\rho_{s}^{j},r_{s}^{j})$ in the
expression of the estimator $\lambda_{c}$. The 95% C.L. upper limit on
$g_{a\gamma}$ at mass $m_{a}$ is obtained from $\lambda_{c}=2.71$.

Figure 2: 95% C.L. upper limits on the effective ALP-two-photon coupling
$g_{a\gamma}$ as a function of the ALP mass, derived in this work for the five
dwarf spheroidal galaxies under investigation, and in the combined case (last
panel). The limits are computed assuming an NFW spatial profile for the ALP
DM.

Results are shown in Figs. 2 (NFW profile) and 4 (coreNFWtides). We see that
the different targets provide similar bounds, which also further motivates us
to perform the combined analysis. The coupling $g_{a\gamma}$ is constrained at
a level around $10^{-12}\leavevmode\nobreak\ \rm{GeV}^{-1}$ with significant
fluctuations between adjacent masses, due to the noise from the process of
subtracting the foreground emission lines. Such rapid variation is more
pronounced at lower masses/longer wavelengths reflecting the presence of
strong OH emission lines from the night sky in this wavelength range. The
bounds improve slightly from low to high masses, which is due to the scaling
of the decay rate with $m_{a}^{3}$, mitigated by an opposite energy dependence
of the observational capabilities (angular and energy resolutions,
foreground).

Figure 3: Bounds for the Leo T galaxy as a function of the integration radius
for 20 channels with wavelengths equally spaced in the range [5324.57,
5562.07] Å, and using a coreNWFtides profile. The vertical line marks the
best-fit value of the core radius.

By comparing Figs. 2 and 4, we see that the reduction of constraining power
for the coreNFWtides profile is very limited. This is because we are probing a
relatively large portion of the targets, implying that the majority of the
pixels entering the statistical analysis are not from the central region,
where the two profiles differ, but from distances where the two profiles
basically coincide. To further illustrate this point, we show in Fig. 3 how
the bounds change as a function of the integration radius $r_{int}$,
corresponding to the radius of the circular area used to compute the bound (we
remind that the default value in our analysis corresponds to an angular radius
of $60^{\prime\prime}$). Each of the 20 lines corresponds to one channel,
while the vertical line marks the best-fit value of the core radius $r_{c}$.
As a general trend, the bound becomes more stringent once regions larger than
the core radius are considered.

Figure 4: Same as Fig. 2 but for the coreNFWtides DM profile.

The robustness of our results against different masking and error estimates is
tested in the same way as discussed in [14]. We find negligible differences in
the derived bounds from the alternative analyses, with results very similar to
what shown in Fig. 3 of [14].

In the left panel of Fig. 5, we summarize our findings for the NFW profile and
include the bound derived in Ref. [13] from the observation of clusters, in
Ref. [35] from the ratio of horizontal branch (HB) to Asymptotic Giant Branch
stars in globular clusters and, for reference, the preferred region for the
QCD axion [36]. In the wavelength/mass range covered by our analysis, we can
confidently exclude the QCD axion, which is also in tension with other
astrophysical and laboratory probes associated with couplings different from
$g_{a\gamma}$, see e.g. [37], and the possible interpretation of near-infrared
background anisotropies in terms of ALP dark matter [19].

In the right panel of Fig. 5, we compare the results of our combined analysis
to Ref. [14], obtained from the MUSE data of Leo T. The current analysis is
typically more conservative, even though bounds are at a comparable level, and
this is mainly due to the treatment of the DM profile. Indeed, in Ref. [14],
the profile was derived by extrapolating results from a Jeans analysis at
larger radii [38], while here it is derived directly from data. We found that,
in the case of Leo T, the extrapolation slightly overshoots the real DM
profile. On top of that, the uncertainty associated with the profile
determination is now taken into account in a more rigorous statistical way, as
described above.

Figure 5: Left: bounds derived in this work assuming an NFW profile exclude
$g_{a\gamma}$ in the blue shaded region. The yellow region marks the QCD axion
band [36]. We show also bounds from star evolution in globular clusters from
[35] (red dash-dotted line), from observation of galaxy clusters [13] (green
solid line), and the 99% C.L. constraint from $\gamma$-ray attenuation [39]
(cyan dotted line). Right: comparison with the bounds derived by the MUSE
observation of Leo T in [14].

For what concerns possible evidence of an ALP signal, we define, again at any
given mass,
$\lambda_{d}=2\ln[\mathcal{L}(g_{a\gamma}=0,\vec{\rho}_{s}^{\,lbf},\vec{r}_{s}^{\,lbf})/\mathcal{L}(g_{a\gamma}^{b.f.},\vec{\rho}_{s}^{\,gbf},\vec{r}_{s}^{\,gbf})]$.
The ALP discovery would occur if $\sqrt{\lambda_{d}}>5$.

Due to the imperfect subtraction of emission lines from the night sky, many
channels present large values of $\sqrt{\lambda_{d}}$. In order to identify
possible emission peaks due to the presence of an ALP, we need to remove this
spurious evidence from our data, i.e. determine which channels are
“unreliable”.222Note that, in those channels, the bounds described above are
weakened and not strengthened by the presence of such residual emission. As a
first step, we look at the sky spectrum of Leo T, obtained from a data cube
without sky subtraction, by summing the flux over an area with no bright
sources. The sky spectrum presents large peaks above a non-zero frequency-
dependent baseline level. We search for peaks with heights above a threshold
that we choose to be five times the standard deviation of fluctuations about
the baseline in a region without large emission lines. With this criterion, we
identify 225 peaks, all of which correspond to known atmospheric emission
lines, except two, which we will not consider as sources of fake evidence in
the following. Next, in order to determine how many channels are affected by a
bright emission line, we construct an estimate of the reliability of the
errors used in our analysis.

For this purpose, we use the individual Leo T exposures. We reduce the data
for 11 individual exposures, following the data reduction process described in
Sec. 2. By assuming that the data are perfectly aligned across exposures, we
have:
${f_{i}(\lambda)=f_{true}(\lambda)+N(0,\sigma_{i}(\lambda)^{2})}$, where
$f_{i}(\lambda)$ is the flux measured at a given wavelength $\lambda$,
$f_{true}(\lambda)$ is the true flux at a given wavelength $\lambda$, and
$N(0,\sigma_{i}(\lambda)^{2})$ represents the normally distributed random
noise at each wavelength, with variance $\sigma_{i}(\lambda)^{2}$. We use the
spectra of 220 stars in the Leo T field of view, and, for each star, we
compute $f_{true}(\lambda)$ as the average of the 11 observations (therefore
$i$ runs over 220 stars and 11 exposures). It is important to note that this
modeling approach is only suitable for sources, as the sky background may
vary. On the other hand, the fluxes from stars are relatively stable, and
$f_{true}(\lambda)$ should be nearly constant. If the uncertainty estimates
$\sigma_{i}(\lambda)$ are reliable, then the difference between the observed
flux and the true flux, normalized by $\sigma_{i}(\lambda)$, should be
distributed approximately as a Gaussian with zero mean and standard deviation
$\Sigma(\lambda)=1$. Next, we superimpose the sky spectrum of Leo T with
$\Sigma(\lambda)$. We observe that around the frequencies corresponding to
atmospheric emission lines, $\Sigma$ deviates significantly from 1. To
determine the characteristic range of wavelengths affected, we look at three
well-known isolated oxygen emission lines ($\lambda_{P}=5577.3,$ $6300.3,$
$6363.8\leavevmode\nobreak\ \rm{\AA}$). We find that $\Sigma(\lambda)-1>0.05$
in a range $[\lambda_{P}-\delta\lambda,\lambda_{P}+\delta\lambda]$, with
$\delta\lambda=5.7\sigma_{\lambda}$ and $\lambda_{P}$ being the wavelength of
the peak and $\sigma_{\lambda}$ being the spectral resolution as in Eq. (3.1).

We thus conduct the evidence search excluding channels that fall into the
$\pm\delta\lambda$ region around all background emission peaks. With this
procedure, approximately 2000 channels out of 3719 are excluded. Finally, we
correct for the Look Elsewhere Effect (LEE) by dividing the p-value by a
factor $N_{trials}$ equal to the total number of used channels divided by the
number of channels falling within the spectral resolution. We compute the
p-value from the test statistics $\lambda_{d}$ defined above, combining all
targets.

We find no evidence for ALP DM in our data, i.e., no case with
$\sqrt{\tilde{\lambda}_{d}}>5$, where $\tilde{\lambda}_{d}$ has the same
meaning of $\lambda_{d}$ but corrected for the LEE.

##  5  Conclusions

Most ALP models predict a coupling between photons and ALPs. This implies that
we expect a monochromatic photon flux generated by ALP decays inside
astrophysical structures. Nearby dwarf spheroidal galaxies are ideal targets
for this search since they are DM-dominated and are relatively close to us.
Assuming ALPs to constitute all the DM in galaxy halos, we analyzed MUSE
spectroscopic observations of five dwarf spheroidal galaxies to search for ALP
radiative decays in the mass range 2.7-5.3 eV. The excellent spectral
resolution and sensitivity of the spectroscopic observations obtained with the
MUSE instrument at the VLT allowed us to probe quite faint and diffuse
monochromatic line emissions.

We tested the possible presence of an ALP DM signal, concluding that none of
the channels selected for this analysis, i.e., not affected by large
background contamination, is exhibiting a detection. We derived robust bounds
on the effective ALP-two-photon coupling consistent among the five galaxies
considered. We took into account the uncertainties associated with the DM
spatial distribution in each dwarf galaxy by including a profile likelihood
depending on the DM profile parameters derived from the Jeans analysis of Ref.
[24]. We considered two possible profiles: an NFW and a cored profile,
yielding comparable bounds. A different description of the DM profile from
what considered here might lead to slightly different bounds. However, we do
not expect large modifications since the constraining power comes from dSph
regions where the determination of the DM density from kinematic data appears
to be robust and so weakly dependent on the parameterization.

The resulting bounds lie well below the QCD axion band, and are significantly
more constraining than limits from other probes, in the relevant mass range.

## Acknowledgements

MT acknowledges support from the research grant ‘The Dark Universe: A Synergic
Multimessenger Approach’ No. 2017X7X85K funded by MIUR. MR, JR, MT and ET
acknowledge support from the project “Theoretical Astroparticle Physics
(TAsP)” funded by the INFN.

MR, JR and ET acknowledge support from ‘Departments of Excellence 2018-2022’
grant awarded by the Italian Ministry of Education, University and Research
(miur) L. 232/2016 and Research grant ‘From Darklight to Dark Matter:
understanding the galaxy/matter connection to measure the Universe’ No.
20179P3PKJ funded by miur. JB and DV acknowledge support by Fundação para a
Ciência e a Tecnologia (FCT) through the research grants UIDB/04434/2020 and
UIDP/04434/2020 and through grant PTDC/FIS-AST/4862/2020. JB acknowledges work
contract 2020.03379.CEECIND and DV acknowledges support from the Fundação para
a Ciência e a Tecnologia (FCT) through the Fellowship 2022.13277.BD.

Based on observations made with ESO Telescopes at the La Silla Paranal
Observatory under programme IDs 0100.D-0807, 0101.D-0300, 0102.D-0372 and
0103.D-0705.

## References

  * [1] J. Preskill, M. B. Wise, and F. Wilczek, Cosmology of the Invisible Axion,  Phys. Lett. B120 (1983) 127–132. [,URL(1982)].
  * [2] L. F. Abbott and P. Sikivie, A Cosmological Bound on the Invisible Axion,  Phys. Lett. B120 (1983) 133–136. [,URL(1982)].
  * [3] M. Dine and W. Fischler, The Not So Harmless Axion,  Phys. Lett. B120 (1983) 137–141. [,URL(1982)].
  * [4] R. D. Peccei and H. R. Quinn, CP Conservation in the Presence of Instantons,  Phys. Rev. Lett. 38 (1977) 1440–1443. [,328(1977)].
  * [5] R. D. Peccei and H. R. Quinn, Constraints Imposed by CP Conservation in the Presence of Instantons,  Phys. Rev. D16 (1977) 1791–1797.
  * [6] S. Weinberg, A New Light Boson?,  Phys. Rev. Lett. 40 (1978) 223–226.
  * [7] F. Wilczek, Problem of Strong p and t Invariance in the Presence of Instantons,  Phys. Rev. Lett. 40 (1978) 279–282.
  * [8] P. W. Graham, I. G. Irastorza, S. K. Lamoreaux, A. Lindner, and K. A. van Bibber, Experimental Searches for the Axion and Axion-Like Particles,  Ann. Rev. Nucl. Part. Sci. 65 (2015) 485–514, [1602.00039].
  * [9] I. G. Irastorza and J. Redondo, New experimental approaches in the search for axion-like particles,  Progress in Particle and Nuclear Physics 102 (Sept., 2018) 89–159, [1801.08127].
  * [10] L. Di Luzio, M. Giannotti, E. Nardi, and L. Visinelli, The landscape of QCD axion models,  Phys. Rept. 870 (2020) 1–117, [2003.01100].
  * [11] M. T. Ressell, Limits to the radiative decay of the axion,  Phys. Rev. D44 (1991) 3001–3020.
  * [12] M. A. Bershady, M. T. Ressell, and M. S. Turner, Telescope search for multi-eV axions,  Phys. Rev. Lett. 66 (1991) 1398–1401.
  * [13] D. Grin, G. Covone, J.-P. Kneib, M. Kamionkowski, A. Blain, et al., A Telescope Search for Decaying Relic Axions,  Phys. Rev. D75 (2007) 105018, [astro-ph/0611502].
  * [14] M. Regis, M. Taoso, D. Vaz, J. Brinchmann, S. L. Zoutendijk, et al., Searching for light in the darkness: Bounds on ALP dark matter with the optical MUSE-faint survey,  Physics Letters B 814 (Mar., 2021) 136075, [2009.01310].
  * [15] K. Nakayama and W. Yin, Anisotropic cosmic optical background bound for decaying dark matter in light of the LORRI anomaly,  Phys. Rev. D 106 (2022), no. 10 103505, [2205.01079].
  * [16] P. Carenza, G. Lucente, and E. Vitagliano, Probing the Blue Axion with Cosmic Optical Background Anisotropies, 2301.06560.
  * [17] T. Bessho, Y. Ikeda, and W. Yin, Indirect detection of eV dark matter via infrared spectroscopy,  Phys. Rev. D 106 (2022), no. 9 095025, [2208.05975].
  * [18] W. Yin and K. Hayashi, Indirect Detection of Decaying Dark Matter with High Angular Resolution: Case for axion search by IRCS at Subaru Telescope, 2305.13415.
  * [19] Y. Gong, A. Cooray, K. Mitchell-Wynne, X. Chen, M. Zemcov, et al., Axion decay and anisotropy of near-IR extragalactic background light,  Astrophys. J. 825 (2016), no. 2 104, [1511.01577].
  * [20] O. E. Kalashev, A. Kusenko, and E. Vitagliano, Cosmic infrared background excess from axionlike particles and implications for multimessenger observations of blazars,  Phys. Rev. D 99 (Jan., 2019) 023002, [1808.05613].
  * [21] A. Caputo, A. Vittino, N. Fornengo, M. Regis, and M. Taoso, Searching for axion-like particle decay in the near-infrared background: an updated analysis,  JCAP 05 (2021) 046, [2012.09179].
  * [22] J. L. Bernal, G. Sato-Polito, and M. Kamionkowski, Cosmic Optical Background Excess, Dark Matter, and Line-Intensity Mapping,  Phys. Rev. Lett. 129 (2022), no. 23 231301, [2203.11236].
  * [23] R. Bacon, S. Conseil, D. Mary, J. Brinchmann, M. Shepherd, et al., The MUSE Hubble Ultra Deep Field Survey. I. Survey description, data reduction, and source detection,  Astron. Astrophys. 608 (Nov., 2017) A1, [1710.03002].
  * [24] S. L. Zoutendijk, M. P. Júlio, J. Brinchmann, J. I. Read, D. Vaz, et al., The MUSE-Faint survey. III. No large dark-matter cores and no significant tidal stripping in ultra-faint dwarf galaxies, arXiv e-prints (Dec., 2021) arXiv:2112.09374, [2112.09374].
  * [25] J. F. Navarro, C. S. Frenk, and S. D. M. White, The Structure of cold dark matter halos,  Astrophys. J. 462 (1996) 563–575, [astro-ph/9508025].
  * [26] M. G. Walker, M. Mateo, and E. W. Olszewski, Stellar Velocities in the Carina, Fornax, Sculptor, and Sextans dSph Galaxies: Data From the Magellan/MMFS Survey,  Astron.J. 137 (Feb., 2009) 3100–3108, [0811.0118].
  * [27] M. G. Coleman, G. S. D. Costa, and J. Bland-Hawthorn, The absence of extratidal structure in the sculptor dwarf spheroidal galaxy,  The Astronomical Journal 130 (sep, 2005) 1065.
  * [28] P. M. Weilbacher, R. Palsa, O. Streicher, R. Bacon, T. Urrutia, et al., The Data Processing Pipeline for the MUSE Instrument, arXiv e-prints (June, 2020) arXiv:2006.08638, [2006.08638].
  * [29] R. Bacon and et al., The MUSE Hubble Ultra Deep Field Survey. I. Survey description, data reduction, and source detection,  Astron. Astrophys. 608 (Nov., 2017) A1, [1710.03002].
  * [30] E. Bertin and S. Arnouts, SExtractor: Software for source extraction.,  Astron. Astrophys. Suppl. 117 (Jun, 1996) 393–404.
  * [31] L. Piqueras, S. Conseil, M. Shepherd, R. Bacon, F. Leclercq, et al., MPDAF - A Python package for the analysis of VLT/MUSE data, arXiv e-prints (Oct., 2017) arXiv:1710.03554, [1710.03554].
  * [32] K. Hayashi, K. Ichikawa, S. Matsumoto, M. Ibe, M. N. Ishigaki, et al., Dark matter annihilation and decay from non-spherical dark halos in galactic dwarf satellites,  Mon. Not. R. Astron. Soc. 461 (2016), no. 3 2914–2928, [1603.08046].
  * [33] M. D. Schneider, C. S. Frenk, and S. Cole, The shapes and alignments of dark matter halos,  JCAP 2012 (May, 2012) 030, [1111.5616].
  * [34] J. I. Read, M. G. Walker, and P. Steger, The case for a cold dark matter cusp in Draco,  Mon. Not. Roy. Astron. Soc. 481 (2018), no. 1 860–877, [1805.06934].
  * [35] M. J. Dolan, F. J. Hiskens, and R. R. Volkas, Advancing globular cluster constraints on the axion-photon coupling,  JCAP 10 (2022) 096, [2207.03102].
  * [36] L. Di Luzio, F. Mescia, and E. Nardi, Redefining the Axion Window,  Phys. Rev. Lett. 118 (2017), no. 3 031801, [1610.07593].
  * [37] Particle Data Group, M. Tanabashi et al., Review of Particle Physics,  Phys. Rev. D 98 (2018), no. 3 030001.
  * [38] V. Bonnivard et al., Dark matter annihilation and decay in dwarf spheroidal galaxies: The classical and ultrafaint dSphs,  Mon. Not. R. Astron. Soc. 453 (2015), no. 1 849–867, [1504.02048].
  * [39] J. L. Bernal, A. Caputo, G. Sato-Polito, J. Mirocha, and M. Kamionkowski, Seeking dark matter with $\gamma$-ray attenuation,  Phys. Rev. D 107 (2023), no. 10 103046, [2208.13794].
  * [40] J. Binney, The radius-dependence of velocity dispersion in elliptical galaxies,  Monthly Notices of the Royal Astronomical Society 190 (04, 1980) 873–880, [https://academic.oup.com/mnras/article-pdf/190/4/873/18194217/mnras190-0873.pdf].
  * [41] J. Binney and G. A. Mamon, M/L and velocity anisotropy from observations of spherical galaxies, or must M87 have a massive black hole?,  Monthly Notices of the Royal Astronomical Society 200 (09, 1982) 361–375, [https://academic.oup.com/mnras/article-pdf/200/2/361/9403099/mnras200-0361.pdf].
  * [42] J. I. Read and P. Steger, How to break the density-anisotropy degeneracy in spherical stellar systems,  MNRAS 471 (Nov., 2017) 4541–4558, [1701.04833].
  * [43] J. Reynoso-Cordova, M. Regis, and M. Taoso, Upper limits on the dark matter content in globular clusters,  JCAP 10 (2022) 038, [2203.13735].
  * [44] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, emcee: The MCMC Hammer,  Publ.Astron.Soc.Pac. 125 (Mar., 2013) 306, [1202.3665].
  * [45] A. Genina, J. I. Read, C. S. Frenk, S. Cole, A. Benítez-Llambay, et al., To $\beta$ or not to $\beta$: can higher order Jeans analysis break the mass-anisotropy degeneracy in simulated dwarfs?,  MNRAS 498 (Oct., 2020) 144–163, [1911.09124].
  * [46] M. L. M. Collins, J. I. Read, R. A. Ibata, R. M. Rich, N. F. Martin, et al., Andromeda XXI - a dwarf galaxy in a low-density dark matter halo,  Mon.Not.Roy.Astron.Soc. 505 (Aug., 2021) 5686–5701, [2102.11890].

##  Appendix A Jeans analysis of Sculptor

Figure 6: Posterior distribution of the dark matter parameters entering the
NFW profile in Eq. 3.2. Units are
$\log_{10}(\rho_{s}/(M_{\odot}/\rm{kcp}^{3}))$ and
$\log_{10}(r_{s}/\rm{kpc})$. Figure 7: The same as in Fig. 6 but using a
coreNFWtides profile. Units are $\log_{10}(M_{200}/(M_{\odot}))$, $c_{200}$,
$\log_{10}(r_{c}/[\mathrm{kpc}])$, $n$, $\log_{10}(r_{t}/[\mathrm{kpc}])$,
$\delta$.

The data we have used in our Jeans analysis for Sculptor are the line-of-sight
velocities reported in [26] alongside the photometric measurements presented
in [27] for surface density data. We model the stellar dynamics assuming a
spherically symmetric and non-collisional Jeans equation [40, 41]:

$\frac{1}{\nu(r)}\frac{\partial}{\partial
r}\left(\nu(r)\sigma_{r}^{2}\right)+\frac{2\beta(r)\sigma_{r}^{2}}{r}=\frac{GM(<r)}{r^{2}},$
(A.1)

where $\nu(r)$ is the stellar density, $\sigma_{r}$ the radial velocity
dispersion, $\beta(r)$ the velocity anisotropy and $M(<r)$ is the total
enclosed mass within a radius $r$ from the center of the target. Furthermore,
we model the radial stellar density profile as a sum of three Plummer spheres
[41]:

$\nu(r)=\sum_{j=1}^{3}\frac{3M_{j}}{4\pi
a_{j}^{3}}\left(1+\frac{r^{2}}{a_{j}^{2}}\right)^{-5/2},$ (A.2)

$M_{j}$ and $a_{j}$ are parameters which can be constrained through the data.
The expression A.2 can be seen as a density expansion, analogous to a Gaussian
decomposition. The velocity anisotropy $\beta(r)$ is parametrized as:

$\beta(r)=\beta_{0}+(\beta_{\infty}-\beta_{0})\,\frac{1}{1+\left(r_{a}/r\right)^{\eta}},$
(A.3)

where $\beta_{0}$, $\beta_{\infty}$, $r_{a}$ and $\eta$ are also free
parameters describing respectively the inner and outer orbital anisotropy, the
radius and the sharpness of the transition.

The mass distribution is the sum of the DM component and the stellar
contribution, which is modeled as in Eq. A.2 but with a free parameter fixing
the overall normalization (in practice this corresponds to a free mass-to-
light ratio). As explained in Sec. 3, we consider two options for the DM
distribution, namely the NFW profile in Eq. 3.2, and the coreNFWtides model of
[34], which modifies the NFW distribution allowing for the presence of a
central core, and a reduced density beyond a tidal radius. More specifically,
the central core is implemented by modifying the NFW mass distribution as
follows:

$M_{\rm{cNFW}}(<r)=M_{\rm{NFW}}(<r)f^{n},$ (A.4)

with $0\leq n\leq 1$ and

$f=\tanh{\left(\frac{r}{r_{c}}\right)},$ (A.5)

being $r_{c}$ the size of the core. The associated density profile is:

$\rho_{\rm{cNFW}}(r)=f^{n}\,\rho_{\rm{NFW}}(r)+\frac{nf^{n-1}\left(1-f^{2}\right)}{4\pi
r^{2}r_{c}}M_{\rm{NFW}}(<r).$ (A.6)

Furthermore, since the galaxies that we are describing experience strong
gravitational interactions with the host galaxy, tidal stripping is expected
in their external regions. Such an effect can be modeled by further modifying
the enclosed mass beyond a tidal radius $r_{t}:$

$M_{\rm{cNFWt}}(<r)=\begin{cases}M_{\rm{cNFW}}(<r)&\text{if $r<r_{t}$},\\\
M_{\rm{cNFW}}(r_{t})+4\pi\rho_{\rm{cNFW}}(r_{t})\frac{r_{t}^{3}}{3-\delta}\left[\left(\frac{r}{r_{t}}\right)^{3-\delta}-1\right],&\text{if
$r>r_{t}$},\end{cases}$ (A.7)

which in terms of the density profile reads

$\rho_{\rm{cNFWt}}(r)=\begin{cases}\rho_{\rm{cNFW}}(r)&\text{if $r<r_{t}$},\\\
\rho_{\rm{cNFW}}(r_{t})\left(\frac{r}{r_{t}}\right)^{-\delta},&\text{if
$r>r_{t}$},\end{cases}$ (A.8)

and the external slope $\delta$ is taken to be $\delta\geq 3.$

Given all these ingredients, we can use the radial velocity dispersion
obtained by solving equation A.1 to compute the line-of-sight velocity
dispersion:

$\sigma_{\rm{L.O.S}}^{2}(R)=\frac{2}{\Sigma_{*}(R)}\int_{R}^{\infty}\left(1-\beta\frac{R^{2}}{r^{2}}\right)\frac{\nu(r)\sigma_{r}^{2}r}{\sqrt{r^{2}-R^{2}}}dr,$
(A.9)

being $\Sigma_{*}(R)$ the projected stellar surface density, which can be
expressed as:

$\Sigma_{*}(R)=\sum_{j=1}^{3}\frac{M_{j}}{\pi
a_{j}^{2}}\left(1+\frac{R^{2}}{a_{j}^{2}}\right)^{-2}.$ (A.10)

Finally, we have compared the model to the data by considering the Gaussian
likelihood
$-2\ln{\mathcal{L}}=\chi^{2}_{\rm{L.O.S}}+\chi_{\Sigma_{*}}^{2}+\chi^{2}_{\rm{VSP}_{1}}+\chi^{2}_{\rm{VSP}_{2}}$.
The first two terms correspond to the chi-squared for the line-of-sight
velocity dispersion and surface density data, respectively. The last two
contributions are constructed from the fourth moments of the velocities
distribution, known as the virial shape parameters $\rm{VSP1}_{1}$ and
$\rm{VSP}_{2},$ which can be computed using equations [20-23] from Ref. [42].
See Refs. [42, 43] for more details. The model is based on a total of 13 (17)
free parameters: 2 (6) for NFW (coreNFWtides), 4 for the velocity anisotropy,
and 7 for the stellar component. We have explored this parameter space through
a Markov Chain Monte Carlo (MCMC) simulation. The numerical analysis is
performed through the emcee [44] sampler which is implemented by using public
python code pyGravSphere [45]. Prior to the MCMC, the stellar model is fitted
to the surface density data. Then, in the MCMC the corresponding parameters
$a_{j}$ and $M_{j}$ are allowed to vary in a 50% range around their previously
determined surface density best-fit values. The dark matter profile A.8 was
added to the already existent profiles in pyGravSphere using the same priors
as in [24].

We present our results for the case of an NFW profile in Fig. 6, where we show
the posterior probability distributions of
$\log_{10}(\rho_{s}/(M_{\odot}/\rm{kpc}^{3}))$ and
$\log_{10}(r_{s}/\rm{kpc}),$ alongside their 0.16, 0.5 and 0.84 percentile
values.

Analogous results are shown in Fig. 7 for the coreNFWtides distribution. In
this case, in order to follow the numerical implementation of the public code
GravSphere [42, 34, 46], we parametrize the NFW profile through the
concentration $c_{\Delta}$ and the mass $M_{\Delta},$ which are related to
$\rho_{s}$ and $r_{s}$ by the equations:

$c_{\Delta}=r_{\Delta}/r_{s}\;\;\;,\;\;\;\rho_{s}=\frac{\Delta\rho_{c}c_{\Delta}^{3}}{3\left(\ln{(1+c_{\Delta})-\frac{c_{\Delta}}{1+c_{\Delta}}}\right)},$
(A.11)

where $\rho_{c}$ is the critical density of the Universe, and $\Delta=200$,
which in turn defines $r_{200}$ as the radius at which the DM density in the
halo is 200 times the critical density of the Universe

$r_{200}=\left(\frac{3}{800}\frac{M_{200}}{\pi\rho_{c}}\right)^{1/3}.$ (A.12)

Taking into account the additional parameters described in Eqs. A.4-A.8, the
coreNFWtides profile is defined by six parameters:
$\log_{10}(M_{200}/[M_{\odot}])$, $\log_{10}(r_{200}/[\rm{kpc}])$,
$\log_{10}(r_{c}/[\rm{kpc}])$, $n$, $\log_{10}(r_{t}/[\rm{kpc}])$ and
$\delta$. We show their posterior distribution from our analysis in Fig.7.
Given the large number of parameters, this case shows significant degeneracy
among the model parameters. On the other hand, the two most relevant ones,
$M_{200}$ and $c_{200}$ are suitably well constrained.

##  Appendix B Caveats of the coreNFWtides analysis

Figure 8: Comparison between bounds obtained assuming a coreNFWtides profile
for the Leo T galaxy performing a full scan of the profile parameter space
(red diamonds) and a partial scan with four of the profile parameters fixed to
their best-fit values as in Figure 4 (blue dots).

For the derivation of the bounds of Figure 4 under the assumption of the
coreNFWtides profile Eq. (A.6), we have fixed $r_{c}$, $n$, $r_{t}$, and
$\delta$ to their best-fit values from the Jeans analysis of Ref. [24].
Although, in principle, all the parameters defining the profile distribution
should be varied, the computational cost of scanning the eight-dimensional
(six profile parameters, $g_{a\gamma}$ and $S_{flat}$) parameter space would
render the analysis very expensive. In particular, the bottleneck is the
calculation of the likelihood from the MUSE data Eq. (4.1), which involves a
sum over a large number of pixels.

As discussed in Section 4, most of the constraining power comes from
intermediate regions, located outside the core radius of the coreNFWtides
profile as shown in Figure 3. In these regions the difference between the
density profiles accounted for with our simplified procedure and with a full
scan of all the parameters is limited.

To further check the validity of this argument, we perform a full parameter
space scan for 6 channels of Leo T. The results are shown in Figure 8. As
expected, the change in the bounds compared to Figure 4 is minor.