# Updated bounds on Axion-Like Particle Dark Matter with the optical MUSE-
Faint survey

###### Abstract

Bounds are derived on the axion-like particle (ALP) to two-photon coupling in
the mass range $2.65-5.27$ eV. The bounds are obtained by searching for the
signal from ALP decay in the Multi Unit Spectroscopic Explorer (MUSE)
observations of five dwarf spheroidal galaxies, under the assumption that ALPs
constitute the dark matter component of the haloes. These bounds are of the
same order and improve on the robustness of those of Reference [1], and
currently represent the strongest bounds within the considered mass range.

###### keywords:

ALP dark matter, dwarf spheroidal galaxy, ALP decay into photons

Elisa Todarello

Dipartimento di Fisica, Università di Torino
via Giuria 1, I–10125 Torino, Italy
Istituto Nazionale di Fisica Nucleare, Sezione di Torino
via P. Giuria 1, I–10125 Torino, Italy
Email Address<EMAIL_ADDRESS>

##  1 Introduction

Axion-like particles (ALPs), just like QCD axions, couple to photons through
the interaction Lagrangian

$\mathcal{L}=-\frac{1}{4}gaF_{\mu\nu}\tilde{F}^{\mu\nu}\enspace,$ (1)

where $\tilde{F}^{\mu\nu}$ is the dual of the electromagnetic field strength
tensor $F^{\mu\nu}$, $a$ is the ALP (or axion) field, and $g$ is a coupling
constant with dimensions of inverse energy. The interaction described by
Equation (1) is used in a wide range of laboratory, cosmological, and
astrophysical searches for ALPs, see References [2, 3, 4] for reviews. One
process that can occur as a result of this interaction is the decay of an ALP
into two photons, which we use in this work to establish upper bounds on $g$.

We search for the diffuse photon emission from ALP decay in five dwarf
galaxies targeted by the MUSE-Faint survey [5]: Eridanus 2, Grus 1, Hydra II,
Leo T, and Sculptor. Dwarf satellite galaxies are ideal targets for dark
matter searches, thanks to their high light-to-mass ratio and proximity to
Earth. We selected these particular five dwarfs because information on the
dark matter profile, needed to calculate the expected flux, is available
either from the MUSE collaboration itself [6] or can be found in the
literature.

The photons from ALP dark matter decay have an energy of about half the ALP
mass. Since the MUSE instrument measures in the optical band [7], we can
constrain $g$ for masses between 2.65 and 5.27 eV. Other bounds in a similar
mass range have been derived from stellar evolution in globular clusters [8,
9], other telescope searches for ALP decay focused on galaxy clusters [10],
and anisotropies in the cosmic optical background [11, 12].

In this work, we update the bounds presented in Reference [1], which were
obtained from the MUSE observation of only one dwarf, namely Leo T. Using
multiple targets enables us to establish more reliable limits. All five dwarfs
provide comparable bounds of the order of $g<10^{-12}$ GeV-1, increasing our
confidence in the validity of such a constraint. We also provide combined
bounds, which are slightly more stringent than those from any of the single
dwarfs.

The analysis presented here does not take into account the uncertainty on the
dark matter profile, whose parameters we fix to the best-fit values form [6]
or, in the case of Sculptor, from a new analysis to be presented in an
upcoming publication [13]. While Ref. [13] will include the likelihood of the
dark matter profile in the statistical analysis, the present manuscript offers
a more pedagogical discussion of the bounds and their derivation.

##  2 Observations and data from the MUSE instrument

The Multi Unit Spectroscopic Explorer (MUSE) located at the Very Large
Telescope (VLT) is a large-field medium-resolution Integral Field
Spectrograph. It provides data in the optical band in the wavelength range
$470-935$ nm, with a sampling of 1.25 Å for a total of 3721 channels. The
spectral resolution has a full width at half maximum $\Delta\lambda$ varying
with wavelength between $2.5$ and $3.0$ Å. The data is in the form of a
“cube”, the first two dimensions being the angular coordinates on the sky and
the third dimension being the wavelength. The field of view is $1\times 1$
arcmin2 with a spatial resolution of $\sim 0.6$ arcsec.

Known emission lines from the night sky are subtracted from the data using the
MUSE Data Reduction Software (DRS) [5]. The subtraction cannot be executed
perfectly, something we take into account in our analysis by adding an
arbitrary spatially flat term to the expected flux, as emission from the night
sky is not expected to vary over the angular scales of the image. The DRS also
provides the r.m.s. error estimate that we’ll use for our data analysis.

For each galaxy, we mask in every channel the brightest sources, such as
stars. To create the mask we run SExtractor [14] on the white-light image,
i.e. the image obtained by integrating the cube over the wavelength axis.

##  3 Expected flux

The rate of ALP decay into two photons is

$\Gamma=\frac{g^{2}m^{3}}{64\pi}\enspace,$ (2)

where $m$ is the ALP mass. Assuming spherical symmetry, the expected flux
density at the detector from a direction $\theta$ on the sky is

$S_{\lambda}(\theta)=\frac{\Gamma}{4\pi}\,\frac{1}{\sqrt{2\pi}\sigma_{\lambda}}e^{-\frac{(\lambda-\lambda_{obs})^{2}}{2\sigma_{\lambda}}}\int
d\Omega\,d\ell\,\rho[r(\theta,\Omega,\ell)]\,B(\Omega)\enspace,$ (3)

where $\rho(r)$ is the dark matter energy density as a function of the
distance from the center of the dwarf, $\lambda_{obs}$ is twice the Compton
wavelength for the ALP mass under consideration, Doppler shifted as described
below, $\sigma_{\lambda}=\Delta\lambda/(2\sqrt{2\ln(2)})$ is the spectral
resolution, and $B(\Omega)$ is the detector’s angular response. The integral
is taken over the line of sight $\ell$ and the position inside the angular
beam $\Omega$.

In the ALP rest frame, the photons are emitted back to back each carrying
energy $m/2$. In the lab frame, the emission line is smeared out due to the
dark matter velocity dispersion in the dwarf, and Doppler-shifted because of
the relative velocity between the dwarf and Earth. For the galaxies we
consider, the velocity dispersion is ${\delta v\lesssim 10^{-4}}$, while the
resolution in wavelength of the MUSE instrument is
$\sigma_{\lambda}/\lambda\geq 1.2\times 10^{-4}$. We thus neglect the velocity
dispersion and assume that all decay photons have a frequency of $m/2$ when
emitted.

We cannot instead neglect the radial velocities of the dwarfs. These
velocities are reported in Table 1 of the Appendix. When converted to natural
units, we see that for all dwarfs except Leo T, the radial velocity is well
above the best spectral resolution. Therefore, photons from ALP decay would be
detected in a different frequency channel than the one with frequency $m/2$.
To account for this, in our calculation of the bounds presented below, we
apply a correction for each dwarf

$\nu_{em}=\nu_{obs}(1+v_{radial})\enspace,$ (4)

where $\nu_{obs}$ is the frequency of the MUSE channel being considered, and
$\nu_{em}$ is the corresponding frequency at emission. We then report the
bounds as a function of the ALP mass $m=4\pi\nu_{em}$.

##  4 Data analysis and results

We assume a cored dark matter density profile for the five dwarfs. The
parameters are taken to be the best-fit values from the GravSphere Jeans
analysis of [6] for Eridanus 2, Grus 1, Hydra II, and Leo T, while for
Sculptor we use the best-fit values from an analysis to be presented in [13].
This new analysis makes use of data from [15, 16]. The expression for the
energy density as a function of radius, as well as the complete set of
parameter values, are presented in the Appendix.

For our statistical analysis, we use a methodology completely analogous to
that of [1]. For each galaxy, and for each frequency channel, we compute the
likelihood

$\mathcal{L}=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}\enspace.$
(5)

The sum runs over pixels in the image. It is divided by the number of pixels
contained within the detector’s angular resolution $N_{pix}^{FWHM}$, as the
fluxes in these pixels are correlated. We assume that the theoretically
predicted flux density $S_{th}^{i}$ at a wavelength $\lambda$ consists of a
flat component $S_{\lambda,flat}$ constant across the image plus
$S_{\lambda}^{i}$ from Equation (3). $S_{\lambda,flat}$ accounts for imperfect
subtraction of emission from the night sky. The noise estimate
$\sigma_{rms}^{i}$ is derived through the MUSE DRS. For each channel, we scan
a grid of values of $\Gamma$ and $S_{\lambda,flat}$.

Figure 1: Illustration of the profile likelihood methodology used to obtain
upper bounds on $g$, shown for the channel of wavelength 5325 Å of Leo T as an
example. On the bottom panel, the color gradient represents the value of
$\tilde{q}(\Gamma,S_{\lambda,flat})$. On the vertical axis,
$S^{(0)}_{\lambda,flat}$ is the best-fit value of $S_{\lambda,flat}$ in the
absence of dark matter. We use $S^{(0)}_{\lambda,flat}$ as an initial guess to
determine the range of values to scan. The red line marks, for each $\Gamma$,
the lowest value of $\tilde{q}$, while the red dot marks the best fit
$\hat{\Gamma},~{}\hat{S}_{\lambda,flat}$. Notice that, for this channel,
$S^{(0)}_{\lambda,flat}$ is negative. As we increase $\Gamma$, larger negative
values of $S_{\lambda,flat}$ are preferred to compensate for the fact that the
diffuse emission from dark matter is not a good fit. The top panel shows the
resulting profile likelihood $q(\Gamma)$ (solid blue). Again, the dot marks
the best-fit value $q(\Gamma)=0$, while the intersection between the thin
dotted line and the blue line corresponds to the 95% C.L. upper bound on
$\Gamma$ and thus on $g$. Figure 2: Upper bounds on the ALP-photon coupling
$g$ at the 95% C.L. derived from MUSE-Faint observations of five dwarf
spheroidal galaxies assuming a ‘coreNFWtides’ profile, and combined bounds
(solid black, bottom-right panel). On the combined panel we show the QCD axion
band [17] in yellow. We also show bounds from star evolution in globular
clusters from [9] (red dash-dotted line) and from observation of galaxy
clusters [10] (blue solid line). The grey line on the combined panel
represents the bounds of [1] derived assuming an NFW profile for Leo T and
using the best-fit D-factor from [18] (see main text for details).

To establish our bounds, we choose a frequentist approach that does not
require incorporating prior knowledge on the allowed parameter ranges. We are
concerned with constraining the ALP-photon coupling $g$, while we are not
interested in $S_{\lambda,flat}$. A commonly employed technique to deal with
such nuisance parameters within the frequentist approach is the profile
likelihood. The profile likelihood function is obtained by maximizing the
likelihood with respect to the nuisance parameter for each fixed value of the
parameter of interest. This yields a one-dimensional slice through the
likelihood surface, as illustrated on the bottom panel of Figure 1 for the
channel of wavelength 5325 Åof Leo T. On the axes of the bottom panel, we see
the grid of scanned values. The color gradient represents the value of
${\tilde{q}(\Gamma,S_{\lambda,flat})=-2\ln\left(\mathcal{L}(g,S_{\lambda,flat})/\mathcal{L}(\hat{g},\hat{S}_{\lambda,flat})\right)}$,
where a hat indicates the best-fit value. For each value of $\Gamma$, we
select the lowest value of $\tilde{q}(\Gamma,S_{\lambda,flat})$ corresponding
to any of the values of $S_{\lambda,flat}$. Such values are marked by the red
line. The result of the profiling procedure is shown on the top panel of
Figure 1. The blue line represents the likelihood ratio

$q(g)=-2\ln\left(\frac{\mathcal{L}(g,\hat{\hat{S}}_{\lambda,flat})}{\mathcal{L}(\hat{g},\hat{S}_{\lambda,flat})}\right)\enspace,$
(6)

where double hat indicates that the parameter has been profiled out.

We use the following test statistics

$t(g)=\begin{cases}q(g)&\qquad g\geq\hat{g}\\\ 0&\qquad
g<\hat{g}\enspace.\end{cases}$ (7)

To obtain upper bounds on $g$, we assume that the probability density function
for $t(g)$ is a half $\chi^{2}$ with one degree of freedom, meaning that the
probability of $t(g)$ being larger than a given value $t_{c}$ is
$P(t_{c})=\int_{t_{c}}^{\infty}dt^{\prime}e^{-t^{\prime}/2}/(2\sqrt{2\pi
t^{\prime}})$. To set the 95% C.L. upper bound, we search for the lowest value
of $g$ for which $t(g)>2.71$. This is also depicted on the top panel of Figure
1. For the combined bound, we multiply the likelihoods for the single
galaxies, here labeled by $j$, yielding the combined likelihood

$\mathcal{L}^{tot}(g)=\prod_{j=1}^{5}\mathcal{L}^{j}(g,\hat{\hat{S}}_{\lambda,flat})\enspace,$
(8)

where $S_{\lambda,flat}$ has been profiled out for each dwarf independently.

Upper bounds on $g$ at the 95% C.L. are shown in Figure 2. All five dwarfs
give similar bounds in the $10^{-12}$ GeV-1 ballpark. The large oscillations
for low masses are due to bright atmospheric emission lines that saturate the
detector. Since these lines are the same in the observations of all galaxies,
they do not compensate and appear in the combined bound as well. Bounds cannot
be found for $4.12~{}\mathrm{eV}<m<4.27~{}\mathrm{eV}$ because a filter is
applied to block light from the sodium laser of the adaptive optic system.

We notice that the bounds from Leo T alone, as well as the combined bounds,
are less stringent than those of [1]. This is due to the more conservative
dark matter profile assumed here. In Ref. [1], the profile considered was a
Navarro-Frenk-White (NFW), and the integral over the line of sight of the dark
matter density of Leo T was derived from the D-factor reported in [18], while
here we’re using a ‘coreNFWtides’ with parameters derived from the MUSE-Faint
observations of [6]. These observations reach distances from the center of the
dwarf well below the half-light radius (see for example Figure 2 of [6]),
while the data from [18] are on angular scales of order or larger than the
half-light radius. This is one of the factors that makes the present bounds
more robust than those of [1].

Figure 3:  Illustration of how our bounds change when the dark matter profile
is modified in the central region of the galaxy. On the left panel, we see two
possible profiles for the Leo T galaxy, a cuspy NFW and a cored
‘coreNFWtides’. The maximum radius shown corresponds to our integration limit
of 1′ from the center of the galaxy. The two profiles are similar for the
largest radii. This results in the bounds, shown on the right panel, being
almost identical (the blue line is covered by the pink one). See the main text
for more details.

In order to understand how the bounds change when a cuspy profile is
considered, we plot in Figure 3 a comparison between the ‘coreNFWtides’ and
NFW profiles for Leo T. The parameters for the NFW profile are taken from
(Table A.1) of [6] and reported in the Appendix for completeness. On the left
panel, we plot both profiles up to our maximum integration radius $r_{max}$,
corresponding to an angular size of $1^{\prime}$. From Equations (3) and (5),
we see that for fixed $S^{i}_{obs}$, if we increase $\rho$ in a given pixel,
the fit will tend to prefer lower values $\Gamma$, up to variations in
$S_{\lambda,flat}$. Next, we notice that the largest the radius, the more
pixels associated with it, their number increasing proportionally to $r$. Thus
we expect the larger radii to determine the typical variation in the bound on
$g$ as $\rho$ is modified. In the example of Figure 3, the average value of
$\rho_{NFW}/\rho_{coreNFWtides}$, weighted by $r$, is approximately 1.7. So we
should expect the bounds on $g$ using the cored profile to be worse by a
factor of about 1.3. As a matter of fact, the bounds using the two profiles
are basically indistinguishable, and we can attribute the bulk of the
discrepancy with our estimate to the effect of the nuisance parameter
$S_{\lambda,flat}$, which can also compensate for the change in $\rho$. We
conclude that if the core is much smaller than $r_{max}$, and if the profiles
are similar outside of the core region, the bounds will be roughly unaffected.
For larger cores, the bounds will be less stringent for the cored profile by a
factor of order $\lesssim\sqrt{\rho_{NFW}(r_{max})/\rho_{core}(r_{max})}$.

Let us finally address the question of whether the data shows any evidence for
ALP dark matter. For example, the Leo T channel of Figure 1 contains no
evidence, as $q$ does not increase significantly for values of $\Gamma$
approaching zero. However, this is not the case for all channels. In fact, for
all dwarfs, we find numerous channels in which the null hypothesis is strongly
disfavored. This is due, just like the oscillations in the value of the upper
bounds, to bright atmospheric emission lines that cannot be perfectly
subtracted from the data. This renders the search for evidence of ALPs
complicated, as the ‘fake’ evidence due to the atmosphere has to be removed.
Such an analysis is in preparation and will appear in [13].

##  5 Conclusions

We have provided upper bounds on the axion-like particle to photon coupling by
searching for a signal from radiative decay. Our analysis is based on
observations of five dwarf spheroidal galaxies, namely Eridanus 2, Grus 1,
Hydra II, Leo T, and Sculptor, obtained from the MUSE-Faint survey. Unlike
Reference [1], which relied on observations of a single galaxy, we used data
from multiple targets, making our constraints more reliable and conclusive.
Additionally, information on the dark matter profile is taken, except for
Sculptor, from the very same MUSE-Faint observations [6], extending to
distances significantly shorter than the half-light radius of the dwarf. This
contributes to the enhanced reliability of the current constraints compared to
those presented in the previous work of Ref. [1].

We have also discussed in detail the method used for our analysis, namely the
profile likelihood, and the effect using a cuspy versus a cored profile has on
the bounds.

The bounds reported in this manuscript do not consider the uncertainty on the
dark matter profile parameters. In a paper currently in preparation, we will
take into account the likelihood of the dark matter profile of the dwarfs in
the analysis.

Acknowledgements

The author would like to thank the members of the team working on this
project: Marco Regis, Javier Reynoso, and Marco Taoso from University of Turin
and INFN Turin, Jarle Brinchmann, Daniel Vaz, and Sebastiaan Zoutendijk from
the MUSE collaboration. The author also thanks the Galileo Galilei Institute
for Theoretical Physics for its hospitality during the completion of this
work. The author acknowledges support from the ‘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.

##  Appendix A Dark matter profile parameters

| Eri 2 | Gru 1 | Hya II | Leo T | Scu
---|---|---|---|---|---
$\rho_{s}~{}[10^{8}~{}M_{\odot}~{}\mathrm{kpc}^{-3}]$ | 0.171 | 0.345 | 0.0388 | 0.221 | 0.154
$r_{s}~{}[\mathrm{kpc}]$ | 1.43 | 0.581 | 23.1 | 1.51 | 1.81
$r_{c}~{}[\mathrm{kpc}]$ | 0.0827 | 0.0112 | 0.018 | 0.0217 | 0.0493
$n~{}$ | 0.854 | 0.364 | 0.631 | 0.707 | 0.775
$r_{t}~{}[\mathrm{kpc}]$ | 1.09 | 3.73 | 16.3 | 9.78 | 1.67
$\delta$ | 3.35 | 3.79 | 3.94 | 4.20 | 4.22
$d~{}[\mathrm{kpc}]$ | 366 | 125 | 151 | 417 | 86
$v_{radial}~{}[\mathrm{km/s}]$ | 76 | -141 | 303 | 38 | 111
Table 1: Values of the ‘coreNFWtides’ parameters used for the analysis. They
are the best-fit values from the likelihood of the Jeans analysis of [6] using
GravSphere or from an analysis using data from [15, 16] in the case of
Sculptor. The table also shows the value of the distance of the dwarfs $d$,
taken from [6, 19, 20] and the radial velocities from [21].

The bounds shown in Figure 2 are derived assuming the ‘coreNFWtides’ dark
matter density profile used for the GravSphere analysis of [6]. This profile
is constructed starting from an NFW progenitor

$\rho_{\rm{NFW}}(r)=\frac{\rho_{s}}{\left(\frac{r}{r_{s}}\right)\left(1+\frac{r}{r_{s}}\right)^{2}}\enspace.$
(9)

and redistributing some of the mass from the central regions to larger radii
as follows

$\rho_{\mathrm{cNFW}}(r)=f^{n}\rho_{\mathrm{NFW}}(r)+\frac{nf^{n-1}\Bigl{(}1-f^{2}\Bigr{)}}{4\pi
r^{2}r_{\mathrm{c}}}M_{\mathrm{NFW}}(<r)\enspace,$ (10)

where $M_{\mathrm{NFW}}(<r)$ is the mass contained within a radius $r$ if the
density profile is NFW and $f=\tanh(r/r_{\mathrm{c}})$. Moreover, beyond a
tidal radius $r_{t}$, the energy density drops abruptly. The expression for
the full profile is

$\rho_{\mathrm{cNFWt}}(r)=\begin{cases}\rho_{\mathrm{cNFW}}(r)&r<r_{\mathrm{t}}\enspace,\\\
\rho_{\mathrm{cNFW}}(r_{t})\,(r/r_{\mathrm{t}})^{-\delta},&r>r_{\mathrm{t}}\enspace.\end{cases}$
(11)

The values of the NFW profile parameters of Leo T used to produce Figure 3 are

$\rho_{s}=0.227\times 10^{8}~{}M_{\odot}~{}\mathrm{kpc}^{-3}\qquad
r_{s}=1.45~{}\mathrm{kpc}\enspace,$ (12)

taken from Table A.1 of [6].

## References

  * [1] M. Regis, M. Taoso, D. Vaz, J. Brinchmann, S. L. Zoutendijk, N. F. Bouché, M. Steinmetz,  _Phys. Lett. B_ 2021, _814_ 136075\.
  * [2] I. G. Irastorza, J. Redondo,  _Prog. Part. Nucl. Phys._ 2018, _102_ 89\.
  * [3] L. Di Luzio, M. Giannotti, E. Nardi, L. Visinelli,  _Phys. Rept._ 2020, _870_ 1\.
  * [4] P. Sikivie,  _Rev. Mod. Phys._ 2021, _93_ , 1 015004.
  * [5] S. L. Zoutendijk, J. Brinchmann, L. A. Boogaard, M. L. P. Gunawardhana, T.-O. Husser, S. Kamann, A. F. Ramos Padilla, M. M. Roth, R. Bacon, M. den Brok, S. Dreizler, D. Krajnovi´c,  _A &A_ 2020, _635_ A107.
  * [6] S. L. Zoutendijk, M. P. Júlio, J. Brinchmann, J. I. Read, D. Vaz, L. A. Boogaard, N. F. Bouché, D. Krajnović, K. Kuijken, J. Schaye, M. Steinmetz,  _arXiv e-prints_ 2021, arXiv:2112.09374.
  * [7] R. Bacon, S. Conseil, D. Mary, J. Brinchmann, M. Shepherd, M. Akhlaghi, P. M. Weilbacher, L. Piqueras, L. Wisotzki, D. Lagattuta, B. Epinat, A. Guerou, H. Inami, S. Cantalupo, J. B. Courbot, T. Contini, J. Richard, M. Maseda, R. Bouwens, N. Bouché, W. Kollatschny, J. Schaye, R. A. Marino, R. Pello, C. Herenz, B. Guiderdoni, M. Carollo,  _A &A_ 2017, _608_ A1.
  * [8] A. Ayala, I. Domínguez, M. Giannotti, A. Mirizzi, O. Straniero,  _Phys. Rev. Lett._ 2014, _113_ , 19 191302.
  * [9] M. J. Dolan, F. J. Hiskens, R. R. Volkas,  _JCAP_ 2022, _10_ 096\.
  * [10] D. Grin, G. Covone, J.-P. Kneib, M. Kamionkowski, A. Blain, E. Jullo,  _Phys. Rev. D_ 2007, _75_ 105018\.
  * [11] K. Nakayama, W. Yin,  _Phys. Rev. D_ 2022, _106_ , 10 103505.
  * [12] P. Carenza, G. Lucente, E. Vitagliano 2023.
  * [13] E. Todarello, et al.,  _in preparation_ .
  * [14] E. Bertin, S. Arnouts,  _Astron. Astrophys. Suppl._ 1996, _117_ 393\.
  * [15] M. G. Walker, M. Mateo, E. W. Olszewski,  _Astron.J._ 2009, _137_ , 2 3100.
  * [16] M. G. Coleman, G. S. D. Costa, J. Bland-Hawthorn,  _The Astronomical Journal_ 2005, _130_ , 3 1065.
  * [17] L. Di Luzio, F. Mescia, E. Nardi,  _Phys. Rev. Lett._ 2017, _118_ , 3 031801.
  * [18] V. Bonnivard, C. Combet, M. Daniel, S. Funk, A. Geringer-Sameth, J. A. Hinton, D. Maurin, J. I. Read, S. Sarkar, M. G. Walker, M. I. Wilkinson,  _Monthly Notices of the Royal Astronomical Society_ 2015, _453_ , 1 849.
  * [19] S. L. Zoutendijk, J. Brinchmann, N. F. Bouché, M. d. Brok, D. Krajnović, K. Kuijken, M. V. Maseda, J. Schaye,  _Astron. Astrophys._ 2021, _651_ A80.
  * [20] A. W. McConnachie,  _The Astronomical Journal_ 2012, _144_ , 1 4.
  * [21] M. Wenger, et al.,  _Astron. Astrophys. Suppl. Ser._ 2000, _143_ 9\.


A_{t}(\bm{s}_{t})\\}}\cdot\mathbbm{1}_{\\{\bm{c}_{t}\geq\bm{a}_{j}\\}}\cdot\left(\sum_{\bm{s}^{\prime}\in\mathcal{S}}p_{t}(\bm{s}_{t},\bm{s}^{\prime})\cdot\sum_{i\in\mathcal{A}_{j}}\frac{1}{C_{i}}\cdot\sum_{j^{\prime}\in\mathcal{B}_{i}}\nu^{t+1}_{j^{\prime}}(\bm{s}^{\prime})\right)\right]$
$\displaystyle\geq$
$\displaystyle\mathbb{E}\left[\sum_{t\in[T]}\sum_{j\in\hat{A}_{t}(\bm{s}_{t})}\mathbbm{1}_{\\{j\sim\hat{A}_{t}(\bm{s}_{t})\\}}\cdot\mathbbm{1}_{\\{\bm{c}_{t}\geq\bm{a}_{j}\\}}\cdot\left(\sum_{\bm{s}^{\prime}\in\mathcal{S}}p_{t}(\bm{s}_{t},\bm{s}^{\prime})\cdot\sum_{i\in\mathcal{A}_{j}}\frac{1}{C_{i}}\cdot\sum_{j^{\prime}\in\mathcal{B}_{i}}\nu^{t+1}_{j^{\prime}}(\bm{s}^{\prime})\right)\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\sum_{t\in[T]}\sum_{j\in\hat{A}_{t}(\bm{s}_{t})}\phi_{j}(\hat{A}_{t}(\bm{s}_{t}),\bm{s}_{t})\left(\left[r_{j}-\sum_{\bm{s}^{\prime}\in\mathcal{S}}p_{t}(\bm{s}_{t},\bm{s}^{\prime})\sum_{i\in\mathcal{A}_{j}}\frac{1}{C_{i}}\sum_{j^{\prime}\in\mathcal{B}_{i}}\nu^{t+1}_{j^{\prime}}(\bm{s}^{\prime})\right]^{+}\sum_{i\in\mathcal{A}_{j}}\frac{1}{C_{i}}\sum_{t^{\prime}=1}^{t-1}\mathbbm{1}_{\\{j^{\prime}\text{\leavevmode\nobreak\
purchased at time\leavevmode\nobreak\
}t^{\prime},i\in\mathcal{A}_{j^{\prime}}\\}}\right)\right]$ $\displaystyle=$
$\displaystyle\sum_{t\in[T]}\sum_{\bm{s}\in\mathcal{S}}P(\bm{s}_{t}=\bm{s})\cdot\sum_{j\in\hat{A}_{t}(\bm{s})}\phi_{j}(\hat{A}_{t}(\bm{s}),\bm{s})\cdot\left(\left[r_{j}-\sum_{\bm{s}^{\prime}\in\mathcal{S}}p_{t}(\bm{s},\bm{s}^{\prime})\cdot\sum_{i\in\mathcal{A}_{j}}\frac{1}{C_{i}}\cdot\sum_{j^{\prime}\in\mathcal{B}_{i}}\nu^{t+1}_{j^{\prime}}(\bm{s}^{\prime})\right]^{+}\cdot\sum_{i\in\mathcal{A}_{j}}\frac{\mathbb{E}[b_{i,t}|\bm{s}_{t}=\bm{s}]}{C_{i}}\right).$

We apply Markov’s inequality to show that for each $i\in[m]$, it holds

$\frac{\mathbb{E}[b_{i,t}|\bm{s}_{t}=\bm{s}]}{C_{i}}\geq P(b_{i,t}\geq
C_{i}|\bm{s}_{t}=\bm{s})=P(c_{i,t}\leq 0|\bm{s}_{t}=\bm{s}).$

Then, by the union bound, we have

$\sum_{i\in\mathcal{A}_{j}}\frac{\mathbb{E}[b_{i,t}|\bm{s}_{t}=\bm{s}]}{C_{i}}\geq\sum_{i\in\mathcal{A}_{j}}P(c_{i,t}\leq
0|\bm{s}_{t}=\bm{s})\geq P\left(\exists
i\in\mathcal{A}_{j}:c_{i,t}<a_{i,j}|\bm{s}_{t}=\bm{s}\right)=1-P(\bm{c}_{t}\geq\bm{a}_{j}|\bm{s}_{t}=\bm{s}).$

Therefore, we know that

$\displaystyle\mathbb{E}[\text{term II}]$ (56) $\displaystyle=$
$\displaystyle\sum_{t\in[T]}\sum_{\bm{s}\in\mathcal{S}}P(\bm{s}_{t}=\bm{s})\cdot\sum_{j\in\hat{A}_{t}(\bm{s})}\phi_{j}(\hat{A}_{t}(\bm{s}),\bm{s})\cdot\left(\left[r_{j}-\sum_{\bm{s}^{\prime}\in\mathcal{S}}p_{t}(\bm{s},\bm{s}^{\prime})\cdot\sum_{i\in\mathcal{A}_{j}}\frac{1}{C_{i}}\cdot\sum_{j^{\prime}\in\mathcal{B}_{i}}\nu^{t+1}_{j^{\prime}}(\bm{s}^{\prime})\right]^{+}\cdot\sum_{i\in\mathcal{A}_{j}}\frac{\mathbb{E}[b_{i,t}|\bm{s}_{t}=\bm{s}]}{C_{i}}\right)$
$\displaystyle\geq$
$\displaystyle\sum_{t\in[T]}\sum_{\bm{s}\in\mathcal{S}}P(\bm{s}_{t}=\bm{s})\cdot\sum_{j\in\hat{A}_{t}(\bm{s})}\phi_{j}(\hat{A}_{t}(\bm{s}),\bm{s})\cdot\left(\left[r_{j}-\sum_{\bm{s}^{\prime}\in\mathcal{S}}p_{t}(\bm{s},\bm{s}^{\prime})\cdot\sum_{i\in\mathcal{A}_{j}}\frac{1}{C_{i}}\cdot\sum_{j^{\prime}\in\mathcal{B}_{i}}\nu^{t+1}_{j^{\prime}}(\bm{s}^{\prime})\right]^{+}\cdot\left(1-P(\bm{c}_{t}\geq\bm{a}_{j}|\bm{s}_{t}=\bm{s})\right)\right).$

Combining (53), (55), and (56), we have

$\displaystyle\mathbb{E}\left[\sum_{t\in[T]}\sum_{j\in
A_{t}(\bm{s}_{t})}\mathbbm{1}_{\\{j\sim A_{t}(\bm{s}_{t})\\}}\cdot
r_{j}\right]$ (57) $\displaystyle\geq$
$\displaystyle\sum_{t\in[T]}\sum_{\bm{s}\in\mathcal{S}}P(\bm{s}_{t}=\bm{s})\cdot\sum_{j\in\hat{A}_{t}(\bm{s})}\phi_{j}(\hat{A}_{t}(\bm{s}),\bm{s})\cdot\left(\left[r_{j}-\sum_{\bm{s}^{\prime}\in\mathcal{S}}p_{t}(\bm{s},\bm{s}^{\prime})\cdot\sum_{i\in\mathcal{A}_{j}}\frac{1}{C_{i}}\cdot\sum_{j^{\prime}\in\mathcal{B}_{i}}\nu^{t+1}_{j^{\prime}}(\bm{s}^{\prime})\right]^{+}\right).$

Finally, from the bid price computing rule (20) and by noting
$\nu^{T+1}_{j}(\bm{s})=0$ for any $j\in[n]$ and $\bm{s}\in\mathcal{S}$, we
also have

$\sum_{\bm{s}\in\mathcal{S}}p_{1}(\bm{s})\cdot\sum_{j\in[n]}\nu^{1}_{j}(\bm{s})=\sum_{t\in[T]}\sum_{\bm{s}\in\mathcal{S}}P(\bm{s}_{t}=\bm{s})\cdot\sum_{j\in\hat{A}_{t}(\bm{s})}\phi_{j}(\hat{A}_{t}(\bm{s}),\bm{s})\cdot\left(\left[r_{j}-\sum_{\bm{s}^{\prime}\in\mathcal{S}}p_{t}(\bm{s},\bm{s}^{\prime})\cdot\sum_{i\in\mathcal{A}_{j}}\frac{1}{C_{i}}\cdot\sum_{j^{\prime}\in\mathcal{B}_{i}}\nu^{t+1}_{j^{\prime}}(\bm{s}^{\prime})\right]^{+}\right).$
(58)

Comparing (57) and (58), we have

$\mathbb{E}\left[\sum_{t\in[T]}\sum_{j\in
A_{t}(\bm{s}_{t})}\mathbbm{1}_{\\{j\sim A_{t}(\bm{s}_{t})\\}}\cdot
r_{j}\right]\geq\sum_{\bm{s}\in\mathcal{S}}p_{1}(\bm{s})\cdot\sum_{j\in[n]}\nu^{1}_{j}(\bm{s})$

which completes our proof.

$\Halmos$

##### Proof of Lemma 12.2.

We first show that
$\\{\hat{\beta}^{t}_{i}(\bm{s}),\hat{\theta}^{t}(\bm{s}),\forall
i\in[m],\forall t\in[T],\forall\bm{s}\in\mathcal{S}\\}$ defined in (16) and
(17) is feasible to LP (19) with the bid price $\nu^{t}_{j}(\bm{s})$ defined
in (20). Note that from the definition in (20), we must have
$\nu^{t}_{j}(\bm{s})\geq\sum_{\bm{s}^{\prime}\in\mathcal{S}}p_{t}(\bm{s},\bm{s}^{\prime})\cdot\nu^{t+1}_{j}(\bm{s}^{\prime})$
for any $j\in[n]$, $t\in[T]$ and $\bm{s}\in\mathcal{S}$, which implies that

$\hat{\beta}^{t}_{i}(\bm{s})\geq\sum_{\bm{s}^{\prime}\in\mathcal{S}}p_{t}(\bm{s},\bm{s}^{\prime})\cdot\hat{\beta}^{t+1}_{i}(\bm{s}^{\prime}).$

Therefore, in order for the constraint (19b) to be satisfied, it is sufficient
to show that

$\hat{\theta}^{t}(\bm{s})-\sum_{\bm{s}^{\prime}\in\mathcal{S}}p_{t}(\bm{s},\bm{s}^{\prime})\cdot\hat{\theta}^{t+1}(\bm{s}^{\prime})\geq\sum_{j\in
A}\phi_{j}(A,\bm{s})\cdot\left[r_{j}-\sum_{\bm{s}^{\prime}\in\mathcal{S}}p_{t}(\bm{s},\bm{s}^{\prime})\cdot\sum_{i\in[m]}a_{i,j}\cdot\hat{\beta}^{t+1}_{i}(\bm{s}^{\prime})\right]^{+},$

holds for any $\forall A\in\mathcal{F}$. Note that

$\displaystyle\hat{\theta}^{t}(\bm{s})-\sum_{\bm{s}^{\prime}\in\mathcal{S}}p_{t}(\bm{s},\bm{s}^{\prime})\cdot\hat{\theta}^{t+1}(\bm{s}^{\prime})$
$\displaystyle=\sum_{j\in[n]}\nu^{t}_{j}(\bm{s})-\sum_{\bm{s}^{\prime}\in\mathcal{S}}p_{t}(\bm{s},\bm{s}^{\prime})\cdot\sum_{j\in[n]}\nu^{t}_{j}(\bm{s}^{\prime})$
$\displaystyle=\sum_{j\in\hat{A}_{t}(\bm{s})}\phi_{j}(\hat{A}_{t}(\bm{s}),\bm{s})\cdot\left[r_{j}-\sum_{\bm{s}^{\prime}\in\mathcal{S}}p_{t}(\bm{s},\bm{s}^{\prime})\cdot\sum_{i\in\mathcal{A}_{j}}\frac{1}{C_{i}}\cdot\sum_{j^{\prime}\in\mathcal{B}_{i}}\nu^{t+1}_{j^{\prime}}(\bm{s}^{\prime})\right]^{+}$
$\displaystyle\geq\sum_{j\in
A}\phi_{j}(A,\bm{s})\cdot\left[r_{j}-\sum_{\bm{s}^{\prime}\in\mathcal{S}}p_{t}(\bm{s},\bm{s}^{\prime})\cdot\sum_{i\in\mathcal{A}_{j}}\frac{1}{C_{i}}\cdot\sum_{j^{\prime}\in\mathcal{B}_{i}}\nu^{t+1}_{j^{\prime}}(\bm{s}^{\prime})\right]^{+}$
$\displaystyle=\sum_{j\in
A}\phi_{j}(A,\bm{s})\cdot\left[r_{j}-\sum_{\bm{s}^{\prime}\in\mathcal{S}}p_{t}(\bm{s},\bm{s}^{\prime})\cdot\sum_{i\in\mathcal{A}_{j}}\beta^{t+1}_{i}(\bm{s}^{\prime})\right]^{+}$
$\displaystyle=\sum_{j\in
A}\phi_{j}(A,\bm{s})\cdot\left[r_{j}-\sum_{\bm{s}^{\prime}\in\mathcal{S}}p_{t}(\bm{s},\bm{s}^{\prime})\cdot\sum_{i\in[m]}a_{i,j}\cdot\hat{\beta}^{t+1}_{i}(\bm{s}^{\prime})\right]^{+}$

for all $A\in\mathcal{F}$, where the first inequality follows from the
definition of $\hat{A}_{t}(\bm{s})$ in (21). Therefore, we show that
$\\{\hat{\beta}^{t}_{i}(\bm{s}),\hat{\theta}^{t}(\bm{s}),\forall
i\in[m],\forall t\in[T],\forall\bm{s}\in\mathcal{S}\\}$ defined in (16) and
(17) with the bid price $\nu^{t}_{j}(\bm{s})$ defined in (20) is feasible to
LP (19).

The objective value of LP (19) under the feasible solution
$\\{\hat{\beta}^{t}_{i}(\bm{s}),\hat{\theta}^{t}(\bm{s}),\forall
i\in[m],\forall t\in[T],\forall\bm{s}\in\mathcal{S}\\}$ can be upper bounded
as follows.

$\displaystyle\sum_{\bm{s}\in\mathcal{S}}p_{1}(\bm{s})\cdot\left(\hat{\theta}^{1}(\bm{s})+\sum_{i\in[m]}C_{i}\cdot\hat{\beta}^{1}_{i}(\bm{s})\right)$
$\displaystyle=\sum_{\bm{s}\in\mathcal{S}}p_{1}(\bm{s})\cdot\sum_{j\in[n]}\nu^{1}_{j}(\bm{s})+\sum_{\bm{s}\in\mathcal{S}}p_{1}(\bm{s})\cdot\sum_{i\in[m]}\sum_{j\in\mathcal{B}_{i}}\nu^{t}_{j}(\bm{s})$
$\displaystyle=\sum_{\bm{s}\in\mathcal{S}}p_{1}(\bm{s})\cdot\sum_{j\in[n]}\nu^{1}_{j}(\bm{s})+\sum_{\bm{s}\in\mathcal{S}}p_{1}(\bm{s})\cdot\sum_{j\in[n]}\sum_{i\in\mathcal{A}_{j}}\nu^{t}_{j}(\bm{s})$
$\displaystyle\leq(1+L)\cdot\sum_{\bm{s}\in\mathcal{S}}p_{1}(\bm{s})\cdot\sum_{j\in[n]}\nu^{1}_{j}(\bm{s})$

where it follows from definition (13) that the number of elements in the set
$\mathcal{A}_{j}$ is upper bounded by $L$ for any $j\in[n]$. Our proof is thus
completed.

$\Halmos$