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

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.


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 L = − 1 4 g aγ a F µν Fµν , where a is the ALP field, F µν is the electromagnetic field strength, Fµν its dual, and g aγ 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 [15,16,[19][20][21][22].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.

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 II1 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×1 arcmin 2 field of view with a spatial sampling of 0.2 arcsec pixel −1 .The data cover a wavelength range of 4700 − 9350 Å, sampled at a resolution of 1.25 Å.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 Å (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 Å 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 σ 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 whitelight 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σ, 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.

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 λ produced by decays of ALPs from a given di- rection identified by θ can be computed as: The decay rate Γ a depends on the ALP mass m a and the effective ALP-two-photon coupling g aγ .In natural units, it reads Γ a = g 2 aγ m 3 a /(64π).The wavelength of emission can be computed as λ em = c/ν em with ν em = m a /(4π).In Eq. 3.1, we neglect the velocity dispersion of ALPs in the dwarf halo, since it is typically ≲ 10 −4 c, namely smaller than the spectral resolution of MUSE, σ λ /λ which varies between 1.2 × 10 −4 and 2.7 × 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 λ obs is then given by λ obs = λ 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(Ω).
Under the assumption of spherical symmetry, the DM density ρ 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 ℓ 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]: where ρ 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, ρ 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γ and m a .In total, there are four free parameters describing the expected flux from ALPs.

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 L j Jeans depending only on ρ 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): where S i th is the theoretical estimate for the flux density in the pixel i, S i obs is the observed flux density and σ i rms 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 λ,f lat 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 ′′ of radius.The number of pixels within the MUSE angular beam is N F W HM pix , which has a size given by the aforementioned FWHM.We define a likelihood L dif f which depends only on ALP parameters by profiling out S j λ,f lat 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: and a combined likelihood considering all five targets simultaneously: To compute L j dif f , we scan a three-dimensional logarithmically spaced grid of values of the following parameters: k = g aγ √ ρ s , r s and S λ,f lat .k is sampled in the range [10 −15 , 10 −10 ] GeV −1 (10 8 M ⊙ kpc −3 ) 1/2 .The values of g aγ and ρ s are then separated by sampling g aγ in the range [10 −13.6 , 10 −10.5 ] GeV −1 , while the ranges of variation of ρ s and r s are different for every galaxy and are determined from the 95% C.L. contour of L j Jeans .Lastly, for each channel, we first compute the best-fit value of S λ,f lat for a vanishing g aγ , and then we scan S λ,f lat in a range from 1/10 to 10 times this best-fit value.
We assume that λ c (g aγ ) = −2 ln[L(g aγ , ⃗ ρ lbf s , ⃗ r lbf s )/L(g b.f.aγ , ⃗ ρ gbf s , ⃗ r gbf s )] follows a χ 2distribution with one d.o.f. and with one-sided probability given by P = ∞ √ λc dχ e −χ 2 /2 / √ 2 π, where g b.f.aγ 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γ = g b.f.aγ , whilst lbf denotes the best-fit of ⃗ ρ s and ⃗ r s for that given g aγ .For the analysis on a single dwarf j, one has just to replace (⃗ ρ s , ⃗ r s ) with (ρ j s , r j s ) in the expression of the estimator λ c .The 95% C.L. upper limit on g aγ at mass m a is obtained from λ c = 2.71.
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γ is constrained at a level around 10 −12 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 3 a , mitigated by an opposite energy dependence of the observational capabilities (angular and energy resolutions, foreground).
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 ′′ ).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.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γ , 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.For what concerns possible evidence of an ALP signal, we define, again at any given mass, Due to the imperfect subtraction of emission lines from the night sky, many channels present large values of √ λ 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". 2As 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 This work Leo T Figure 5: Left: bounds derived in this work assuming an NFW profile exclude g aγ 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 γ-ray attenuation [39] (cyan dotted line).Right: comparison with the bounds derived by the MUSE observation of Leo T in [14].
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 (λ) = f true (λ) + N (0, σ i (λ) 2 ), where f i (λ) is the flux measured at a given wavelength λ, f true (λ) is the true flux at a given wavelength λ, and N (0, σ i (λ) 2 ) represents the normally distributed random noise at each wavelength, with variance σ i (λ) 2 .We use the spectra of 220 stars in the Leo T field of view, and, for each star, we compute f true (λ) 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 (λ) should be nearly constant.If the uncertainty estimates σ i (λ) are reliable, then the difference between the observed flux and the true flux, normalized by σ i (λ), should be distributed approximately as a Gaussian with zero mean and standard deviation Σ(λ) = 1.Next, we superimpose the sky spectrum of Leo T with Σ(λ).We observe that around the frequencies corresponding to atmospheric emission lines, Σ deviates significantly from 1.To determine the characteristic range of wavelengths affected, we look at three well-known isolated oxygen emission lines (λ P = 5577.3,6300.3,6363.8Å).We find that Σ(λ) − 1 > 0.05 in a range [λ P − δλ, λ P + δλ], with δλ = 5.7σ λ and λ P being the wavelength of the peak and σ λ being the spectral resolution as in Eq. (3.1).
We thus conduct the evidence search excluding channels that fall into the ±δλ 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 λ d defined above, combining all targets.We find no evidence for ALP DM in our data, i.e., no case with λd > 5, where λd has the same meaning of λ d but corrected for the LEE.

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.   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.

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

Figure 2 :
Figure 2: 95% C.L. upper limits on the effective ALP-two-photon coupling g aγ 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.

Figure 3 :
Figure3: 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.

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

Figure 8 :
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 Figure4(blue dots).