# The impact of primordial black holes on the stellar mass function of ultra-
faint dwarf galaxies

Nicolas Esser1, Sven De Rijcke2, and Peter Tinyakov1
1Service de Physique Théorique, Université Libre de Bruxelles (ULB), CP225
Boulevard du Triomphe, B-1050 Bruxelles, Belgium
2Department of Physics and Astronomy, Ghent University, Krijgslaan 281, 9000
Ghent, Belgium E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

If primordial black holes constitute the dark matter, stars forming in dark-
matter dominated environments with low velocity dispersions, such as ultra-
faint dwarf galaxies, may capture a black hole at birth. The capture
probability is non-negligible for primordial black holes of masses around
$10^{20}$ g, and increases with stellar mass. Moreover, infected stars are
turned into virtually invisible black holes on cosmologically short
timescales. Hence, the number of observed massive main-sequence stars in
ultra-faint dwarfs should be suppressed if the dark matter was made of
asteroid-mass primordial black holes. This would impact the measured mass
distribution of stars, making it top-light (i.e. depleted in the high-mass
range). Using simulated data that mimic the present-day observational power of
telescopes, we show that already existing measurements of the mass function of
stars in local ultra-faint dwarfs could be used to constrain the fraction of
dark matter composed of primordial black holes in the – currently
unconstrained – mass range of $10^{19}-10^{21}$ g.

###### keywords:

cosmology : dark matter – stars : luminosity function, mass function –
galaxies : dwarf – galaxies : stellar content – methods : statistical

††pubyear: 2023††pagerange: The impact of primordial black holes on the
stellar mass function of ultra-faint dwarf galaxies–8

##  1 Introduction

Since the discovery of gravitational waves from black hole mergers (Abbott et
al., 2019, 2021), the hypothesis of primordial black holes (PBHs) has
attracted a lot of attention in the scientific community. These hypothetical
black holes, first proposed in the late 1960s (Zel’dovich & Novikov, 1967;
Hawking, 1971), could have been produced in the early Universe and potentially
explain part, or the totality, of the dark matter (DM). However, many bounds
have been placed on the fraction $f=\Omega_{\text{PBH}}/\Omega_{\text{DM}}$ of
DM that can be made of PBHs, and currently only PBHs of masses between $\sim
10^{17}$ g and $\sim 10^{23}$ g can still make up the entirety of the DM (Carr
et al., 2021). The black holes in this range of masses are particularly
difficult to constrain as they are massive enough to avoid Hawking
evaporation, but still of microscopic size, and of masses below the numerous
lensing constraints that arise for heavier PBHs.

One way to constrain PBHs in this mass range is to use their capture by stars,
which leads to the accretion of matter and conversion of the infected stars
into sub-solar mass black holes on cosmologically short time scales.
Observationally, this effect would manifest itself as the disappearance of
stars, since the resulting sub-solar mass black holes are basically invisible.
This process has been considered in application to neutron stars and white
dwarfs by Capela et al. (2013a, b, 2014) and to main-sequence stars by Esser &
Tinyakov (2023). The advantage of the latter as PBH detectors is that they
themselves are much easier to observe.

The probability of PBH capture by a star is generally very small. It may only
reach values of order one in exceptional DM-rich environments such as ultra-
faint dwarf galaxies (UFDs) which are known to have high DM density and low
velocity dispersion – a combination that provides the most favourable
conditions for capture. Requiring that no more than half of the stars in the
most promising UFD, Triangulum II, has been converted into black holes would
limit the fraction of PBHs $f$ to at most 40% (Esser & Tinyakov, 2023).

Based on the present-day observations, it is not easy to reliably constrain
the fraction of stars that have been converted into black holes in a given
dwarf galaxy over its history, as one would have to rely on modelling a UFD’s
evolution, which incurs large uncertainties. The purpose of the present paper
is to demonstrate that constraining PBHs through star destruction is
nevertheless feasible. The key observation is that the PBH capture, and hence
the star destruction is significantly more efficient for more massive stars,
and therefore overall this process suppresses the stellar mass function at
larger masses. The stellar mass functions of many UFDs have been measured. We
show that the already existing and near-future data may be sufficient, in
quantity and quality, to constrain the impact of PBHs and, consequently, their
abundance.

In this paper, we follow a proof-of-concept strategy and do not use real data.
Instead, we generate a mock stellar sample which is similar in size, mass
distribution, and other properties to the samples collected for real dwarf
galaxies. We will refer to this sample as the “mock data”. By construction,
this sample is unaffected by PBHs. We then use these mock data to estimate the
constraining power of the stellar mass function with respect to possible
effects of PBHs. To this end, we take as a model a general parametrization of
the stellar mass function, either a broken power law or a log-normal
distribution, each depending on a number of free parameters, and apply the
effect of star destruction by PBHs. The fraction $f$ of PBHs is treated as a
free parameter of the model. We use a Bayesian approach to fit the model to
the mock data, and to determine the likely values of the model parameters, and
in particular of the PBH-DM fraction $f$. With this approach, we find that it
is possible to exclude the value $f=1$ for the PBH-free mock data with
reasonable confidence.

The rest of this paper is organized as follows. In Sec. 2, we summarize the
observational facts about dwarf galaxies relevant for this work. In Sec. 3, we
calculate the probability of the PBH capture by a star, as a function of the
star mass. In Sec. 4, we discuss the initial mass function and construct a
corresponding model of the present-day mass function that includes the effect
of PBHs. In Sec. 5 we give the details of the mock data sample, and perform a
Bayesian fit of the model to the mock data. In Sec. 6 we cross-check our
results with the frequentist approach. We will finally conclude in Sec. 7.

##  2 Dark matter and stars in dwarf galaxies

The ultra-faint dwarf galaxies are the most DM-dense systems known, which
makes them ideal candidates for constraining the abundance of PBHs. It was
shown that the star formation in such systems was shut down at the time of
reionization, around 1 billion years after the big bang (Brown et al., 2014).
To a good approximation, the star population of such galaxies can thus be
considered as isochrone, i.e. all the stars in a given UFD were born at the
same time $\sim 12.8$ billion years ago. Over this time, stars more massive
than $\sim 0.8M_{\odot}$ would have evolved off the main sequence. Due to the
faintness of their remnants, the impact of PBHs on the population of stars
above this mass is complicated to measure, and we will concentrate on stars
below this threshold in what follows.

Multiple photometric studies of low-mass stars in ultra-faint dwarf galaxies
have been carried out in the last two decades. Using the stellar luminosities
measured in different colour bands, the colour-magnitude diagram of these
galaxies was obtained in order to study their chemical and dynamical
properties (see e.g. Conn et al., 2018; Filion et al., 2020; Baumgardt et al.,
2022). While all of these works involved star counting in the observed
galaxies (typically containing between a few hundred and a few thousand
individually resolved stars), a fit of the mass distribution of stars in these
galaxies was only performed in Feltzing et al. (2002); Geha et al. (2013);
Gennaro et al. (2018a, b); Filion et al. (2022). Due to the limited
sensitivity of telescopes, it was impossible for a long time to detect
individual stars with masses below $\sim 0.5M_{\odot}$ in UFDs. However,
recent observations were able to resolve stars with masses as low as $\sim
0.2M_{\odot}$ in some of these galaxies. Preliminary results from the James
Webb Space Telescope even showed a resolution that reached $0.09M_{\odot}$ in
the dwarf galaxy Draco II, located at a distance of $\sim 20$kpc (Weisz et
al., 2023).

Other studies, as reviewed by Simon (2019), put bounds on different properties
of these galaxies, such as their metallicity and velocity dispersion, both of
which turn out to be extremely low in these systems. The masses of many UFDs
were also estimated. From these measurements, one finds that the UFDs are
completely dominated by dark matter, and obtains an estimate of the DM
density. In the context of this paper, two parameters are of major interest:
the dark matter density $\rho_{\text{DM}}$, and its velocity dispersion,
$\sigma_{v}$, which is assumed to be the same as the stellar velocity
dispersion. Their combination $\rho_{\text{DM}}/\sigma_{v}^{3}$ controls the
probability of the PBH capture by a star (see next section). When normalized
to the reference values $\rho_{\text{DM}}=100$ GeV/cm3 and $\sigma_{v}=7$ km/s
it gives the “merit factor” $\eta$ introduced in Esser & Tinyakov (2023),

$\eta=\dfrac{\rho_{\rm
DM}}{100~{}\text{GeV}/\text{cm}^{3}}\left(\dfrac{7~{}\text{km/s}}{\sqrt{2}\sigma_{v}}\right)^{3}.$
(1)

Among currently observed UFDs, the galaxy Triangulum II was found to have the
highest merit factor of $\eta=0.95$, while another 4 galaxies (Tucana III,
Draco II, Segue I, Grus I) have merit factors exceeding 0.35.

##  3 Star destruction by PBHs

Main sequence stars capture PBHs most efficiently at the time of their
formation. The capture is an uncorrelated random process, and therefore the
number of PBHs captured by a given star follows a Poisson distribution
characterized by the mean captured PBH number,
$\bar{N}=\bar{M}_{\text{cap}}/m$, $m$ being the mass of the PBH and
$\bar{M}_{\text{cap}}$ the mean PBH mass captured by the star. If a star
captures at least one PBH, it will inevitably be destroyed; the probability of
survival is thus $P_{S}=\exp(-\bar{N})=\exp(-\bar{M}_{\text{cap}}/m)$.

The mean captured mass $\bar{M}_{\text{cap}}$ depends on both the PBH mass $m$
and the star mass $M$, as well as on the local DM density $\rho_{\rm DM}$, the
velocity dispersion $\sigma_{v}$, and on the fraction of PBH $f$ in the total
amount of DM. In a wide mass range $m=10^{19}-10^{21}$ g, which we will focus
on in what follows, the dependence on most of these parameters is trivial. The
captured mass is obviously proportional to the PBH fraction $f$. It is also
proportional to the combination $\rho_{\rm DM}/\sigma_{v}^{3}$ which is, up to
normalization, the merit factor $\eta$ (1). Finally, as has been shown in
Esser & Tinyakov (2023) (cf. their figure 1), in this mass range it is also
proportional to the PBH mass $m$. Therefore, one can write the mean captured
number $\bar{N}$ as

$\bar{N}=\bar{M}_{\text{cap}}/m=f\eta\,\nu(M)$ (2)

where $\nu(M)$ is the mean captured number $\bar{M}_{\text{cap}}/m$ calculated
for the reference values of the DM parameters $\eta=1$ and $f=1$. The function
$\nu(M)$ only depends on the star mass $M$. The star survival probability can
thus be written as

$P_{S}(f,M)=\exp[-f\eta\,\nu(M)].$ (3)

The function $\nu(M)$ plays a key role in our analysis. We show this function
in Fig. 1, computed numerically for $12.8$ Gyr old stars with stellar masses
in the range $M\in[0.2,0.8]M_{\odot}$ using the method described in Esser &
Tinyakov (2023). There, the mean mass $\bar{M}_{\text{cap}}$ of PBHs captured
by a star at birth was computed. The time spent by the PBH to sink towards the
star and acquire a trajectory which is enclosed in the star was taken into
account in this calculation. The possibility that the PBH could get pushed out
of its trajectory by other astrophysical objects in the vicinity of the star
was also considered. Finally, it was checked that, once the PBH is completely
trapped inside, it destroys the star in a finite time.

Figure 1:  Mean captured PBH number $\nu$ for the reference values $\eta=f=1$
as a function of the stellar mass $M$. The statistical uncertainties of the
simulations are too small to be visible on the plot. The blue line shows the
fit with the function (4).

Note that $\nu(M)$ was computed for uniform dark matter and stellar
distributions within the dwarf galaxy. We however checked that the results
remain similar in case of non-uniform profiles. Since $\nu(M)$ is independent
of the PBH mass in the interval $[10^{19},10^{21}]$ g, these results apply not
only to the case where all PBHs have the same mass, but also to any non-
monochromatic population of black holes whose masses lie between these two
values.

For the sake of numerical efficiency, we fitted the results of the
computations with an analytical function

$\nu(M)=a\left({\frac{M}{M_{\odot}}}\right)^{b}$ (4)

with the parameters $a=2.95$ and $b=2.13$. This fit is shown in blue in Fig.
1.

##  4 Impact of PBHs on the stellar mass function

###  4.1 Present-day mass function model

According to the isochrone hypothesis, the stellar population of UFD galaxies
is born during one initial star-formation burst, and is unchanged after that.
Hence, in the absence of PBHs, and for stars below $\sim 0.8M_{\odot}$ that
did not turn into compact objects yet, the present-day observed mass
distribution is the same as the initial mass function (IMF).

The presence of PBHs changes this picture, as a fraction of the stars in the
galaxy are destroyed at or after the initial star-formation burst. The effect
of PBHs on the present-day stellar mass function can therefore be easily
calculated from the IMF and the probability of PBH capture.

The IMF, i.e. the mass distribution of stars at birth, has been one of the
leading topics of research in stellar astrophysics for several decades. While
there exist some standard and universal fits to the IMF (Kroupa, 2001;
Chabrier, 2003), there likewise exists ample evidence for IMF variations in
dwarf galaxies (Geha et al., 2013; Gennaro et al., 2018a, b).

More generally, the IMF can be modelled with a broken power-law (BPL)
distribution

$\dfrac{dN(\alpha_{1},\alpha_{2},M)}{dM}\propto\begin{cases}&M^{\alpha_{1}}\text{
for }0.2M_{\odot}\leq M<0.5M_{\odot}\\\ &kM^{\alpha_{2}}\text{ for }M\geq
0.5M_{\odot},\end{cases}$ (5)

with the continuity constant $k=(0.5M_{\odot})^{(\alpha_{1}-\alpha_{2})}$ and
$\alpha_{1,2}$ the two power-law exponents. Here we have adopted
$M=0.5M_{\odot}$ for the break mass between the two power laws, as it is the
commonly used value, and current observations (cf. e.g. Gennaro et al.
(2018b)) show little deviation from this number. Note that for
$\alpha_{1}=-1.3$ and $\alpha_{2}=-2.3$, the Kroupa (2001) IMF is recovered.
In case one solely considers stars of masses above $0.5M_{\odot}$, the only
relevant parameter left is $\alpha_{2}$. We will call this particular case the
single power-law (SPL) distribution.

Another model which is often used to describe the IMF is the log-normal (LN)
distribution,

$\dfrac{dN(M_{c},\sigma_{\text{LN}},M)}{dM}\propto\dfrac{1}{M}\exp\left(-\dfrac{\left(\log_{10}(M/M_{c})\right)^{2}}{2\sigma_{\text{LN}}^{2}}\right),$
(6)

where the two free parameters are the characteristic mass $M_{c}$ and the
width $\sigma_{\text{LN}}$ of the distribution. The (single-star) Chabrier
(2003) IMF is recovered for $M_{c}=0.08M_{\odot}$ and
$\sigma_{\text{LN}}=0.69$.

With the isochrone hypothesis in mind, one may write the normalized present-
day mass function (PDMF) of the stars in a given UFD, including the
destruction by PBHs, as follows

$f_{\text{PDMF}}(f,\theta,M)=C\dfrac{dN}{dM}P_{S}$ (7)

with $C$ the normalization factor

$C^{-1}(f,\theta)=\int_{M_{\text{min}}}^{M_{\text{max}}}\dfrac{dN}{dM}P_{S}dM$
(8)

such that the PDMF is normalized to $1$. Here $\theta$ stands for the
parameters of the chosen IMF which are $(\alpha_{1}$, $\alpha_{2})$ for the
BPL; $(M_{c}$, $\sigma_{\text{LN}})$ for the LN; and $\alpha_{2}$ for the SPL.
The factor $P_{S}(f,M)$ is the survival probability of a star against its
destruction by PBHs over the last $12.8$ billion years. Its presence in eq.
(7) quantifies the impact of PBHs on the stellar mass function. Note that
without the isochrone approximation, $P_{S}$ needs to be computed individually
for each star, taking into account its age. To better visualize the impact of
PBHs on the mass function, in Fig. 2 we show, as an example, the Kroupa IMF
(BPL with $\alpha_{1}=-1.3$ and $\alpha_{2}=-2.3$) and the resulting PDMF
after destruction by PBHs. More massive stars are destroyed with higher
probability, which translates into a PDMF with a steeper slope than the IMF.

Figure 2:  Kroupa IMF (blue) and associated PDMF (orange) after the
destruction of a fraction of the stars by PBHs. A logarithmic scale is used on
the vertical axis. The normalization of the IMF, and hence the scale on the
vertical axis, is arbitrary.

###  4.2 Remark on binary systems

The comparison of model mass functions with observations is in general
complicated by the presence of binary systems, which in most cases cannot be
resolved. These unresolved systems appear as a single bright object, and thus
tend to make the observed mass function more top-heavy than what the real
stellar mass distribution is. In most observational mass function studies,
these unresolved systems are taken into account through one extra parameter,
the binary fraction.

The PBH capture by stars in binary stellar systems has not yet been considered
in the literature. It involves three-body dynamics, which may modify the
capture rate. On the one hand, the presence of a companion may deviate the
trajectories that are passing through the star and prevent the capture. On the
other hand, the trajectories not passing through the star may start doing so
as a result of the perturbation by the companion. The net effect of a
companion star, if any, is therefore unclear.

One can take a conservative approach, assuming that all the trajectories that
are sufficiently perturbed are not captured, and calculate the fraction of the
PBH trajectories that are not impacted by the presence of companion stars.
Using the database of trajectories from the previous simulations of Esser &
Tinyakov (2023), we can determine the probability distribution of the PBH
apastrons in their orbits around the star. Assuming that the typical distance
$d$ between the partners in a binary is distributed as $1/d$ for
$d\in[15,10^{5}]$AU (De Rijcke & Dejonghe, 2002) — which again is a
conservative assumption for newly formed stars around which the PBHs orbit —
one can use the 3D analogue of eq. (6) of Esser & Tinyakov (2023) to find the
fraction of PBH orbits that are captured even in the presence of a companion.
We found that this fraction varies between 15% and 65% depending on the
combination of stellar masses and on the PBH mass, being generally larger for
less massive PBHs.

To give an example, in dwarf galaxies Coma Berenices and Boötes I, for the
star mass range similar to what we are using, the binary fractions were found
to be $\sim 0.35$ and $\sim 0.7$, respectively (Gennaro et al., 2018b; Filion
et al., 2022), with $\sim 50$% uncertainties. So, assume for simplicity that
half of the stars are in binaries that are not resolved and appear as single
objects. Then of all observed “stars”, 2/3 are actual stars and capture the
PBH with the rate set by Fig. 1, while 1/3 are unresolved binary systems that
capture the PBH at a rate reduced by, say, 50% on the basis of the previous
discussion. Clearly, the resulting overall loss in the capture rate of $\sim
15\%$ is a correction that for the purposes of this paper can be neglected. We
will thus disregard the issue of binary stars in what follows.

##  5 Bayesian analysis of the PDMF model

###  5.1 Mock data

To mimic the real data, we generate a mock sample of stars (referred to as
“mock data” in what follows) with parameters that resemble the observed ones.
We will later run it through the analysis as if it were real data. Since we
aim at constraining PBHs, we do not include any PBH effects in this sample.

Specifically, the mock isochrone stellar population is generated using the
MIST (Dotter, 2016; Choi et al., 2016) isochrone models and the ArtPop python
package (Greco & Danieli, 2022). All the stars in the sample have a
metallicity $[\text{Fe/H}]=-2.2$ and an age of $12.8$ billion years, which are
values typical of ultra-faint dwarf galaxies (Simon, 2019). They follow the
Kroupa (2001) IMF, i.e. the distribution (5) with $\alpha_{1}=-1.3$ and
$\alpha_{2}=-2.3$.

Based on the review of the existing measurements of individual stars in local
dwarf galaxies, cf. Sec. 2, we choose the mock sample to contain $N=1000$
stars with masses between $M_{\text{min}}=0.2M_{\odot}$ and
$M_{\text{max}}=0.79M_{\odot}$. The number of stars in the mock sample is of
the same order as what has been measured in different UFDs. The upper bound
$M_{\text{max}}=0.79M_{\odot}$ is the mass of the most massive stars, within
the MIST model, that are still on the main sequence in a 12.8 Gyr old galaxy.
The lower bound of $0.2M_{\odot}$ corresponds to the smallest stellar mass
that can currently be resolved in dwarf galaxies. One may expect that recent
and future telescopes, such as the James Webb Space Telescope, will resolve
stars of even lower masses in UFDs in the future (Weisz et al., 2023).

The stellar population of the mock sample is represented on a colour-magnitude
diagram, Fig. 3. The absolute magnitude is displayed on the left vertical
axis, while the associated masses are displayed on the right. The F606W and
F814W filters (from the Advanced Camera for Surveys (ACS) on board of the
Hubble Space Telescope (HST)) are used to describe the luminosity and colour
of the stars. Note that for simplicity, the population is considered at a
distance of $10$pc, so that the apparent and absolute magnitudes are the same.
Reddening and extinction are neglected. These mock data are obviously
idealized, and real observations would show a more scattered distribution
around the isochrone curve. To assess the impact of noise on the data, we have
also generated a – more realistic – scattered mass sample and performed the
analysis described in Appendix 8. The results were very similar to those
obtained with the idealized data depicted in Fig. 3.

Figure 3:  Colour-magnitude diagram of the mock data. The absolute magnitudes
(in the F606W and F814W filters described in the text) are displayed on the
left vertical axis, while the stellar masses are shown on the right vertical
axis. They are a function of the colours of the stars (i.e. the magnitude
difference between the two bands), displayed on the horizontal axis.

###  5.2 Likelihood

Given the current observational capabilities to resolve stars in UFDs as
encoded in the mock data, we want to evaluate if the modifications to the
stellar mass function of these dwarf galaxies due to PBHs are detectable. For
this, we use a Bayesian approach. Given a sample of $N$ stars with individual
masses $M_{i}$, either observed or simulated, one may define a binless
likelihood $L$ that quantifies how well this set follows the distribution (7)
by taking the product of $f_{\text{PDMF}}(f,\theta,M)$ over all stars in the
sample. Its logarithm

$\displaystyle\log
L(f,\theta|\\{M_{i}\\})=\sum_{i=1}^{N}\log(f_{\text{PDMF}}(f,\theta,M_{i}))$
(9)
$\displaystyle=N\log(C(f,\theta))+\sum_{i=1}^{N}\log\left({\dfrac{dN(\theta,M_{i})}{dM}}\right)+\sum_{i=1}^{N}\log\left(P_{S}(f,M_{i})\right),$

is a standard quantity used for parameter estimation in statistical analysis.

As defined above, the log-likelihood function depends on a number of
parameters: the fraction of PBHs in the total amount of the dark matter $f$
and the parameters of the stellar IMF ($(\alpha_{1}$, $\alpha_{2})$ for BPL;
$(M_{c}$, $\sigma_{\text{LN}})$ for LN; $\alpha_{2}$ for SPL). Our primary
interest is the estimation of the fraction $f$. The other, IMF-dependent
parameters are considered as nuisance parameters. As for the merit factor, it
was fixed for simplicity to the highest measured value of $\eta=0.95$. Note
however that from eq. (3) $\eta$ and $f$ are completely degenerate parameters.
Hence, one can easily rescale $\eta$ to any value by adapting accordingly the
result obtained for $f$.

###  5.3 Bayesian analysis

We compare our mock data sample to three IMF models: BPL, LN and SPL. In each
of these models, we include the star destruction by PBHs, considering their
fractional abundance $f$ as a free parameter to be determined from the fit
together with other parameters of the model. The goal is to see how strongly
the value of $f$ is constrained by the mock data which resembles the existing
observations and, by construction, has no PBH effects in it. .

To estimate the model parameters, including $f$, we use the likelihood $L$
defined in (9). For the BPL and LN cases, the log-likelihood depends on 3
parameters: the fraction $f$ and the two IMF parameters: $\alpha_{1}$ and
$\alpha_{2}$ for the BPL distribution, and $M_{c}$ and $\sigma_{\text{LN}}$
for the LN one. For the SPL case, the likelihood depends only on $f$ and
$\alpha_{2}$. In the latter case, the range of stellar masses is restricted to
$M\geq 0.5M_{\odot}$ in order to mimic the data in which stars below
$0.5M_{\odot}$ are not resolved. Keeping only the stars above this threshold
reduces the size of our mock data sample for the SPL case from $N=1000$ to
$N=233$ stars.

Based on the existing measurements of the IMFs of ultra-faint dwarfs galaxies
with very low merit factors $\eta$, in which the destruction of stars by PBHs
is expected to be negligible (Feltzing et al., 2002; Geha et al., 2013;
Gennaro et al., 2018a, b; Filion et al., 2022), we constrained the power-law
exponents $\alpha_{1,2}$ to values $\in[-2.6,-0.8]$ and the LN parameters
$M_{c}\in[0.08,0.6]$ and $\sigma_{\text{LN}}\in[0.5,0.7]$, as there are very
few measured parameters outside these bounds. Note that these bounds are
conservative regarding the presence of unresolved binary systems, as they
encompass some observations which fit the single-star mass function and others
that fit the system mass function. On the other hand, the fraction $f$ is by
definition in the range $[0,1]$. Using uniform priors $\Pi(f,\theta)$ for all
the parameters within their bounds, the posterior probability density function
(pdf) can be computed for every IMF model as

$\text{pdf}(\\{M_{i}\\}|f,\theta)=\dfrac{\Pi(f,\theta)L(f,\theta|\\{M_{i}\\})}{\int_{f^{\prime},\theta^{\prime}}\Pi(f^{\prime},\theta^{\prime})L(f^{\prime},\theta^{\prime}|\\{M_{i}\\})},$
(10)

where $\Pi$ is in this case a flat window function for each parameter and
$\theta$ stands for all the parameters of the chosen IMF except $f$.

One can see from eq. (10) that the pdf is equal to the likelihood, to within a
normalization factor. In order to avoid calculating this factor, which is
computationally expensive for a sample size of $N=1000$ stars, we mapped the
pdf numerically by rejection sampling. The method is straightforward: a random
point $(f_{j},\theta_{j})$ is generated within the prior boundaries with
uniform probability. A random uniform number $x_{j}\in[0,\max(L)]$ is also
drawn. If $L(f_{j},\theta_{j}|\\{M_{i}\\})>x_{j}$, the point is kept. It is
discarded otherwise. The process is repeated until $n=10^{6}$ points are
obtained. The distribution of these points in the parameter space then allows
one to map the posterior probability density function.

In the case where the values of the (nuisance) parameters $\theta$ are
unknown, one can marginalize the pdf distribution with respect to these
parameters in order to get the probability density function of $f$ alone,
$\text{pdf}_{\text{marg}}(\\{M_{i}\\}|f)=\int_{\theta^{\prime}}\text{pdf}(\\{M_{i}\\}|f,\theta^{\prime})d\theta^{\prime}$.
This corresponds to drawing a histogram of the distribution of the points
$f_{j}$ previously obtained, regardless of the values of the other coordinates
$\theta_{j}$ of these points. For the BPL, LN and SPL cases respectively,
these histograms are shown on the top left, top right, and bottom panels of
Fig. 4.

Figure 4:  Marginalized probability density function of the PBH fraction $f$
for the BPL (top left panel), LN (top right panel) and SPL (bottom panel)
cases. The green region corresponds to the values that are within $1\sigma$
(i.e. $68.3\%$ of realizations), the blue region corresponds to points between
$1\sigma$ and $2\sigma$ (i.e. $95.5\%$ of realizations) and the orange region
corresponds to points between $2\sigma$ and $3\sigma$ (i.e. $99.7\%$ of
realizations).

One can read from these histograms that the most likely values of $f$ (maximum
of the pdf) are located at (or close to) $f=0$. Our mock data are thus best
described by the model $f_{\text{PDMF}}(f,M)$ for the value $f=0$. This is to
be expected, since the mock data were generated as if the environment is
devoid of PBHs. Moreover, and most importantly, we see that in the BPL and LN
cases, a fraction of DM made of PBHs $f\gtrsim 0.8$ is excluded with a
confidence of $3\sigma$, while $f\gtrsim 0.55$ is excluded at $2\sigma$. In
the SPL case, $f\gtrsim 0.95$ is excluded at $3\sigma$ and $f\gtrsim 0.7$ is
excluded at $2\sigma$.

The first thing we would like to point out regarding these results is that the
exclusion regions resulting from the BPL and LN cases are very similar, even
though the pdf maximum of these two distributions is slightly different. This
is consistent with the fact that these mass functions are very close to one
another. The exclusion resulting from the SPL case is less stringent,
expressing the fact that the latter analysis only considered stars above
$0.5M_{\odot}$, which represents $\sim 1/4$ of the full sample of $1000$
stars. Note however that even this analysis of “poorer data” results in the
exclusion of the case $f=1$, i.e. all DM consisting of PBHs, at the $3\sigma$
confidence level.

Secondly, we recall that the sensitivity depends on the chosen priors for the
nuisance parameters. More stringent constraints would arise from narrower
window functions on $\theta$. The chosen windows used to produce these results
are benchmark values based on observational data, as we lack consistent
theoretical models for the mass function of stars in dwarf galaxies. Future
improvements on the modelling of these parameters would impact the sensitivity
obtained here.

Finally, we would like to emphasize once more that these results are specific
to the generated mock dataset. Another data sample, real or computer-
generated, could give results that fluctuate from what has been obtained here.
We however analysed several other randomly generated mock datasets and found
that $f=1$ was always excluded at $3\sigma$ for $N=1000$.

##  6 Consistency check : frequentist analysis

The goal of this section is to check the robustness of the results obtained
with the Bayesian analysis by comparing them with results from a frequentist
approach. For simplicity, we will focus on the BPL model (eq. (5)) for the IMF
and fix $\alpha_{1}=-1.3$ (as in the Kroupa IMF) in what follows. Since PBHs
destroy more high-mass than low-mass main-sequence stars, the range of stellar
masses that is most impacted by PBHs is the one above $0.5M_{\odot}$, which is
the reason we are keeping $\alpha_{2}$ free. The full pdf distribution of the
present-day mass function then only depends on $f$ and $\alpha_{2}$ and can be
represented easily on a plot.

First, using the Bayesian approach described in the previous section, we
mapped the corresponding pdf, cf. Fig. 5. The maximum of the pdf (MP) located
at $(f=0.19,\alpha_{2}=-1.8)$ is indicated with a purple dot. We then chose a
specific point (call it CP, for “comparison point”), represented with a red
dot, which corresponds to the parameters $(f=0.22,\alpha_{2}=-2.3)$. The
choice of this point is arbitrary. However, note that it is located on the
boundary of the $2\sigma$-exclusion region of this pdf.

Figure 5:  Probability density function for the BPL case with
$\alpha_{1}=-1.3$ fixed. The purple dot corresponds to the location of the pdf
maximum (MP) of the real sample, while the red dot is an arbitrary point used
for comparison (CP). Green, blue and orange regions are as described on figure
4.

Now, in the frequentist approach, we generate many (ideally, an infinite
number) trials of a specific model, with fixed values of the parameters
$(f,\alpha_{2})$ corresponding to the comparison point, and check how often
the pdf maximum of the real sample (in our case, the mock data) is recovered
given this model. In particular, choosing the BPL model with values of the
parameters fixed to the ones of the comparison point, i.e. $f=0.22$ and
$\alpha_{2}=-2.3$, we expect the maximum of the pdf of the real sample to fall
around the $2\sigma$ boundary of the new maxima distribution, in order for the
Bayesian and frequentist approaches to be consistent with each other.

To verify that this is the case, we generated $10^{6}$ fake populations, using
the ArtPop package with the hypotheses described in Sec. 5.1 and with specific
values of the input parameters $(f=0.22,\alpha_{1}=-1.3,\alpha_{2}=-2.3)$. We
computed the maximum likelihood for each of these populations. The density
distribution of these maxima in the parameter space is shown on figure 6.

Figure 6:  Distribution of the maxima of the likelihoods for the BPL case with
$\alpha_{1}=-1.3$ fixed. The fake populations were generated with $f=0.22$ and
$\alpha_{2}=-2.3$. The purple dot corresponds to the location of the pdf
maximum (MP) of the real sample, while the red dot is an arbitrary point used
for comparison (CP). Green, blue and orange regions are as described on figure
4.

We see that, obviously, the point CP is now the maximum likelihood of this new
distribution. More importantly, the point MP falls around the $2\sigma$
border, as expected. The slight offset from this boundary is due to the finite
number of stars in the populations, and the finite number of fake populations
that were generated in the frequentist case.

We have repeated this whole exercise with other IMF models and specific
parameter values, always finding the consistency between Bayesian and
frequentist analyses. Hence, our results are robust and independent of the
statistical approach employed.

##  7 Conclusion

The main conclusion from this work is that it is possible to constrain
primordial black holes in the mass range $10^{19}-10^{21}$ g using already
existing photometric observations of ultra-faint dwarf galaxies. The
constraints, perhaps marginal with the existing data, could be improved in the
future with observations from higher resolution telescopes and a better
understanding of the dynamical properties of UFDs.

In order to demonstrate this, we first computed, based on the results of Esser
& Tinyakov (2023), the probability of stars of different masses being
destroyed by PBHs in the last $12.8$ Gyr. This probability depends on the
fraction of DM made of PBHs $f$, as well as the dynamical properties of the
dark matter in a given UFD. The key observation is that this probability
increases with the mass of the star.

Using the isochrone hypothesis and the fact that the initial and present-day
mass function of stars below $\sim 0.8M_{\odot}$ should be equal in the
absence of PBHs, we computed the impact of PBH destruction on the present-day
mass function. Because of preferential destruction of heavier stars, the
resulting PDMF becomes more top-light, i.e. is suppressed on the high-mass
end. In this way, we computed the model present-day mass functions
corresponding to three different parameterizations of the initial mass
functions.

We then generated a mock dataset, using the ArtPop stellar population
generation package, which mimics current observations of ultra-faint dwarf
galaxies. This mock dataset was not affected by PBHs. We performed a Bayesian
analysis of this dataset, given the PDMF model computed previously, in order
to determine the best-fit values of the PDMF model parameters, including $f$.
We found that the value $f=1$ (all the DM is in the form of PBH) could be
excluded at more than $3\sigma$ for all three IMF considered, even for the
single-power law with poorer data quality. Finally, we checked, by generating
fake populations using the ArtPop package, that the Bayesian results are
consistent with a frequentist approach.

In this paper we have improved on the results of Esser & Tinyakov (2023) in
two respects. First, we have identified a more robust observable sensitive to
the effect of PBH — the UFD stellar mass function — that is less prone to
astrophysical uncertainties than a mere fraction of destroyed stars. Second,
we have developed the Bayesian approach to quantify the statistical
significance of the constraints on the PBH fraction.

As a final remark, we note that the same approach can be used to discover PBHs
rather than constrain their fraction. When applied to the real observations,
one may in principle find that the value $f=0$ is not compatible with the
data, which would be an indication of the PBH existence and a motivation for
dedicated searches for these hypothetical objects.

## Acknowledgements

S.D.R. acknowledges support from grant Segal ANR-19-CE31-0017 of the French
Agence Nationale de la Recherche. N. E. is a FRIA grantee of the Fonds de la
Recherche Scientifique – FNRS. P.T. is supported in part by the Institut
Interuniversitaire des Sciences Nucléaires (IISN) Grant No. 4.4503.15.
Computational resources have been provided by the Consortium des Équipements
de Calcul Intensif (CÉCI), funded by the Fonds de la Recherche Scientifique de
Belgique (F.R.S.-FNRS) under Grant No. 2.5020.11 and by the Walloon Region.
Tag : ULB-TH/23-17.

## References

  * Abbott et al. (2019) Abbott B. P., et al., 2019, Physical Review X, 9, 031040
  * Abbott et al. (2021) Abbott R., et al., 2021, Physical Review X, 11, 021053
  * Baumgardt et al. (2022) Baumgardt H., Faller J., Meinhold N., McGovern-Greco C., Hilker M., 2022, MNRAS, 510, 3531
  * Brown et al. (2014) Brown T. M., et al., 2014, ApJ, 796, 91
  * Capela et al. (2013a) Capela F., Pshirkov M., Tinyakov P., 2013a, Phys. Rev. D, 87, 023507
  * Capela et al. (2013b) Capela F., Pshirkov M., Tinyakov P., 2013b, Phys. Rev. D, 87, 123524
  * Capela et al. (2014) Capela F., Pshirkov M., Tinyakov P., 2014, Phys. Rev. D, 90, 083507
  * Carr et al. (2021) Carr B., Kohri K., Sendouda Y., Yokoyama J., 2021, Reports on Progress in Physics, 84, 116902
  * Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
  * Choi et al. (2016) Choi J., Dotter A., Conroy C., Cantiello M., Paxton B., Johnson B. D., 2016, ApJ, 823, 102
  * Conn et al. (2018) Conn B. C., Jerjen H., Kim D., Schirmer M., 2018, ApJ, 852, 68
  * De Rijcke & Dejonghe (2002) De Rijcke S., Dejonghe H., 2002, MNRAS, 329, 829
  * Dotter (2016) Dotter A., 2016, ApJS, 222, 8
  * Esser & Tinyakov (2023) Esser N., Tinyakov P., 2023, Phys. Rev. D, 107, 103052
  * Feltzing et al. (2002) Feltzing S., Wyse R., Houdashelt M., Gilmore G., 2002, in Grebel E. K., Brandner W., eds, Astronomical Society of the Pacific Conference Series Vol. 285, Modes of Star Formation and the Origin of Field Populations. p. 223 (arXiv:astro-ph/0101370), doi:10.48550/arXiv.astro-ph/0101370
  * Filion et al. (2020) Filion C., Kozhurina-Platais V., Avila R. J., Platais I., Wyse R. F. G., 2020, ApJ, 901, 82
  * Filion et al. (2022) Filion C., Platais I., Wyse R. F. G., Kozhurina-Platais V., 2022, ApJ, 939, 38
  * Geha et al. (2013) Geha M., et al., 2013, ApJ, 771, 29
  * Gennaro et al. (2018a) Gennaro M., et al., 2018a, ApJ, 855, 20
  * Gennaro et al. (2018b) Gennaro M., et al., 2018b, ApJ, 863, 38
  * Greco & Danieli (2022) Greco J. P., Danieli S., 2022, ApJ, 941, 26
  * Hawking (1971) Hawking S., 1971, MNRAS, 152, 75
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Simon (2019) Simon J. D., 2019, ARA&A, 57, 375
  * Weisz et al. (2023) Weisz D. R., et al., 2023, The JWST Resolved Stellar Populations Early Release Science Program. II. Survey Overview (arXiv:2301.04659), doi:10.3847/1538-4365/acdcfd
  * Zel’dovich & Novikov (1967) Zel’dovich Y. B., Novikov I. D., 1967, Soviet Ast., 10, 602

##  Appendix A Scattered sample analysis

In reality, when trying to infer the mass function of an isochrone stellar
population, multiple sources of uncertainties have to be taken into account.
Apart from unresolved binaries, which we already mentioned in the main text,
the stars do not have exactly the same age or metallicity. Hence, they do not
follow exactly the same isochrone curve, which results in some extra scatter
on the color-magnitude diagram. Furthermore, there is for each star an
associated photometric uncertainty. For example, in Filion et al. (2022), the
relative photometric error on the magnitude in the F606W and F814W bands is of
the order of $\sim 0.2\%$. To take these uncertainties into account — the
photometric one, but also the extra scatter due to the non-exactitude of the
isochrone hypothesis — we generate a sample similar to the one used for the
main study of this paper, but with a Gaussian noise of $1\%$ on the F606W and
F814W magnitudes of each star. The amplitude of the Gaussian noise is chosen
in such a way that the resulting color-magnitude diagram is similar to real
diagrams obtained from observational studies (cf. Geha et al. (2013); Gennaro
et al. (2018a, b); Filion et al. (2022)). This color-magnitude diagram is
shown on fig 7.

Figure 7:  Colour-magnitude diagram of the noisy mock data. The absolute
magnitudes (in the F606W and F814W filters described in the text) are
displayed on the left vertical axis, while the stellar masses are shown on the
right vertical axis. They are a function of the colours of the stars,
displayed on the horizontal axis.

We performed the analysis described in section 5.3 on this new sample, and
obtained the results depicted in fig 8. One can see that the exclusion limits
of the fraction $f=1$ are very close to the original ones. Hence, we conclude
that adding scatter to the mock data sample does not change the constraining
power of the method significantly.

Figure 8:  Marginalized probability density function of the PBH fraction $f$
for the BPL (top left panel), LN (top right panel) and SPL (bottom panel)
cases with the scattered sample. The green region corresponds to the values
that are within $1\sigma$ (i.e. $68.3\%$ of realizations), the blue region
corresponds to points between $1\sigma$ and $2\sigma$ (i.e. $95.5\%$ of
realizations) and the orange region corresponds to points between $2\sigma$
and $3\sigma$ (i.e. $99.7\%$ of realizations).


# Retrieval Augmented Generation and Representative Vector Summarization for
Large Unstructured Textual Data in Medical Education

S. S. Manathunga
Department of Pharmacology
University of Peradeniya
Sri Lanka
<EMAIL_ADDRESS>
&Y. A. Illangasekera
Department of Pharmacology
University of Peradeniya
Sri Lanka
<EMAIL_ADDRESS>

###### Abstract

Large Language Models are increasingly being used for various tasks including
content generation and as chatbots. Despite their impressive performances in
general tasks, LLMs need to be aligned when applying for domain specific tasks
to mitigate the problems of hallucination and producing harmful answers.
Retrieval Augmented Generation (RAG) allows to easily attach and manipulate a
non-parametric knowledgebases to LLMs. Applications of RAG in the field of
medical education are discussed in this paper. A combined extractive and
abstractive summarization method for large unstructured textual data using
representative vectors is proposed.

_K_ eywords AI for Medicine $\cdot$ Representative Vectors

##  1 Introduction

Large language models (LLMs) are demonstrating amazing zero-shot learning
capabilities. They have many applications including content generation and
general question answering. These models learn the data they were trained on
very well and the knowledge is retained encoded in the model weights as a
parameterized knowledgebase. This also makes it difficult to revise or update
the knowledge of an LLM later, once it is trained on a large corpus of
training data [1]. The models can hallucinate, producing factually inaccurate,
plausible-sounding answers. Retrieval Augmented Generation (RAG) tries to
mitigate these issues by attaching a non-parametric knowledgebase to the LLM,
in the form of a vector database [2]. Knowledge in the vector database can
easily be modified and different kinds of unstructured textual data can be
vectorized and stored in a single vector database.

RAG-assisted summarization is an important aspect of document intelligence as
well. It poses more challenges than retrieval. A well-written summary should
highlight the key points and provide an overview of the document’s topics.
Generating summaries of large documents directly using LLMs is difficult
because the input content is often too large to fit into the model’s context
window. Performance tends to decrease with increasing context lengths even in
large context window models. Furthermore, the ’Lost in the Middle’ problem,
the phenomenon that the LLMs tend to ignore facts in the middle of the context
and assign more importance to the beginning and the end of the context becomes
worse with large input contexts [3].

To overcome these challenges, RAG-assisted Representative Vector Summarization
(RVS) is introduced in this paper. RVS selects a pre-defined number ($k$) of
representative text chunks from the non-parametric knowledgebase and applies a
combined abstractive and extractive summarizing workflow to generate the final
summary. The value of parameter $k$ is determined based on the maximum token
limit that can be afforded. The model can identify keywords of the text chunks
and their relative importance and create visual representations of the overall
content of the document using word clouds and scatter plots.

RVS is implemented in `docGPT`, a document intelligence program written in
Python using `langchain` framework and the source is available at
https://github.com/ssm123ssm/docGPT-pharm [4].

##  2 Methods

Described here is the workflow implemented in docGPT.

###  2.1 Retrieval

The model extracts text from unstructured sources including PDFs, text
documents, spreadsheets and slide presentations. It can extract text from
images and scanned PDF files using optical character recognition (OCR) [5].
The extracted text is initially stored in memory as a Document and then the
Document is split into text chunks using recursive character text splitting.
The chunks are then embedded in a 1536-dimensional vector space using OpenAI
`text-embedding-ada-002` embedding engine and stored in Facebook AI Similarity
Search (`FAISS`) vector database [3]. The `FAISS` vector database is used to
find the $k$-most similar chunks to a given query at the query time. The
original query, combined with the retrieved chunks is compiled into a prompt
and passed to the LLM for generating the answer. The model workflow and the
steps excecuted at build-time and at query-time are depicted in Figure 1.

Figure 1: Retrieval workflow of docGPT

###  2.2 Summarization

####  2.2.1 Representative Vector

A maximum affordable token limit ($T$) is initially defined. The target is to
select $k$ number of chunks which is the maximum number of chunks that can be
selected without their total token size exceeding $T$, from the vector
database of $n$ chunks. Assuming each chunk has an average token size of $s$,
$k$ is simply obtained by optimizing for $\max_{k\leq n,k\times s\leq T}k$

Once $k$ is calculated, the chunks in the high-dimension vector space are
quantized using the k-means clustering algorithm. $k$-number of clusters which
minimize the within-cluster sum of squared Euclidian distances from the
corresponding centroids are created within the vector space. Since the
distribution of chunks in the vector space is based on their contextual
similarity, we assume each cluster captures different aspects/semantics of the
original document.

After this quantization, one representative chunk, which is the closest to the
corresponding centroid from each cluster is extracted.

Let vectors be an $n\times d$ matrix, where $n$ is the number of chunks and
$d$ is the dimensionality of the embeddings. Let $centroids[i]$ denote the
centroid of $i$-th cluster. The Euclidian distance between the $m$-th data
point and the centroid of the $i$-th cluster is calculated by

$distance_{m}=\sqrt{\sum_{j=1}^{d}(vectors[m][j]-centroids[i][j])^{2}}$ (1)

and stored in the $distances[i]$ array.

Next, the index of the data point that has the minimum distance to the
centroid of the $i$-th cluster, denoted as $closest\\_index_{i}$ , is
calculated by:

$closest\\_index_{i}=\operatorname*{argmin}{distances_{i}}$ (2)

Finally, the $closest\\_index_{i}$ values for all clusters are stored in a
list.

$closest\\_indices=[closest\\_index_{0},closest\\_index_{1},\dots,closest\\_index_{k-1}]$

Once the indices of the representative chunks are identified, they are stored
in a separate representative Document list.

####  2.2.2 Keyword generation and mapping

Even though the token size of the representative text chunks obtained with
this method is relatively smaller compared to the original document, it can
still be too large to fit into the LLM’s context window. For example, `docGPT`
uses a default maximum affordable token limit of 10,000 tokens. The default
LLM that `docGPT` calls is OpenAI `gpt-3.5-turbo`, which has a maximum context
window of 4k. Therefore, the representative Document list undergoes an
intermediate step of extractive summarization (keyword generation and mapping)
for each chunk in it.

Three keywords are generated for each representative chunk and the keywords
are distributed among all the other members of the same cluster. A word cloud
is generated from all the keywords from the n chunks. The relative frequency
that each keyword appears is reflected by the size of the word in the word
cloud.

The original high-dimensional vectors are then reduced to two dimensions using
t-distributed Stochastic Neighbor Embedding (t-SNE) and plotted on a scatter
plot with colors corresponding to their clusters and keywords. This enables
the identification of the distribution of the contents of the entire document
at a glance.

Finally, the mapped summaries are used to create a final abstractive summary
and to create a list of key points from the documents.

##  3 Evaluation

###  3.1 Retrieval

Responses for queries on clinical medicine and general pharmacology from LLM
without a non-parametric knowledgebase and a RAG model with vector databases
built from Kumar and Clark’s Clinical Medicine 10th Edition and British
National Formulary 82 were compared and checked for accuracy [6, 7]. The
responses generated by `chatGPT` and RAG implementation in `docGPT` are
tabulated with the excerpts from the reference books in the supplementary
material named Queries. For each query, `docGPT` generated more targetted and
accurate answers, while `chtGPT` answers were more generic.

###  3.2 Summarization

Performances of RVS and retrieve-stuff-summarize methods were tested using
Kumar and Clark Clinical Medicine 10th edition and BNF 82 ebooks. Results are
summarized in supplementary material named Summaries.

Kumar and Clark had 1508 pages with 13024 text chunks, each having an average
of 789 tokens. The maximum affordable token limit was set at 15,000 and 19
representative chunks were selected. The word cloud and the t-SNE
visualization are depicted in Figures 2 and 3.

Figure 2: Word cloud for Kumar and Clark Clinical Medicine

Figure 3: t-SNE visualization of clustering of chunks in embedding space for
Kumar and Clark Clinical Medicine

BNF had 1805 pages with 7278 text chunks with an average token size of 486.
The model chose 10 representative chunks under the constraints of 5000 maximum
affordable tokens. The word cloud and t-SNE are depicted in Figures 4 and 5.

Figure 4: Word cloud BNF 82

Figure 5: t-SNE visualization of clustering of chunks in embedding space for
BNF

##  4 Implementation

RVS is implemented in `docGPT`, a document intelligence program written in
Python using `langchain` framework and the source is available at
https://github.com/ssm123ssm/docGPT-pharm.

##  5 Conclusion

Clinical medicine is a knowledge-intensive domain. Both clinicians and medical
students would benefit from efficient methods for retrieving information
quickly from large knowledgebases. We believe the proposed retrieval augmented
generation workflow and representative vector summarization for large
documents would be of help in this context. Even though the workflow was
tested on medical reference books and use cases related to medical education,
the concept of RAG and RVS can be adopted by other domains as well.

## References

  * [1] Karan Singhal, Shekoofeh Azizi, Tao Tu, S. Sara Mahdavi, Jason Wei, Hyung Won Chung, Nathan Scales, Ajay Tanwani, Heather Cole-Lewis, Stephen Pfohl, Perry Payne, Martin Seneviratne, Paul Gamble, Chris Kelly, Nathaneal Scharli, Aakanksha Chowdhery, Philip Mansfield, Blaise Aguera y Arcas, Dale Webster, Greg S. Corrado, Yossi Matias, Katherine Chou, Juraj Gottweis, Nenad Tomasev, Yun Liu, Alvin Rajkomar, Joelle Barral, Christopher Semturs, Alan Karthikesalingam, and Vivek Natarajan.  Large language models encode clinical knowledge.
  * [2] Patrick Lewis, Ethan Perez, Aleksandra Piktus, Fabio Petroni, Vladimir Karpukhin, Naman Goyal, Heinrich Küttler, Mike Lewis, Wen-tau Yih, Tim Rocktäschel, Sebastian Riedel, and Douwe Kiela.  Retrieval-augmented generation for knowledge-intensive NLP tasks.
  * [3] Jeff Johnson, Matthijs Douze, and Hervé Jégou.  Billion-scale similarity search with GPUs.
  * [4] Harrison Chase.  LangChain.  original-date: 2022-10-17T02:58:36Z.
  * [5] R. Smith.  An overview of the tesseract OCR engine.  In Ninth International Conference on Document Analysis and Recognition (ICDAR 2007) Vol 2, pages 629–633. IEEE.  ISSN: 1520-5363.
  * [6] Kumar and clark’s clinical medicine, 10th edition - 9780702078682.
  * [7] British national formulary | BNF publications.


# The MURAVES experiment: study of the Vesuvius Great Cone with Muon
Radiography

M. D’Errico  M. D’Errico1,2 F. Ambrosino 1,2 A. Anastasio2 S. Basnet 3 L.
Bonechi 4 M. Bongi4,5 A. Bross 6 R. Ciaranfi4 L. Cimmino1,2 C. Ciulli4,5 R.
D’Alessandro4,5 A. Giammanco3 F. Giudicepietro7 S. Gonzi4,5 R. Karnam3 G.
Macedonio7 V. Masone2 N. Mori4,5 M. Moussawi3 M. Orazi7 G. Passeggio2 R.
Peluso7 A. Pla- Dalmau6 C. Rendon8 A. Samalan8 G. Saracino1,2 G. Scarpato7 P.
Strolin1,2 M. Tytgat8 E. Vertechi7 L.Viliani4,5

( 1 University of Naples Federico II, Naples, Italy
2 INFN sez. di Napoli, Naples, Italy
3 Centre for Cosmology, Particle Physics and Phenomenology, Université
catholique de Louvain, Louvain-la-Neuve, Belgium
4 INFN sez. di Firenze, Florence, Italy
5 University of Florence, Florence, Italy
6 Fermilab, Batavia, IL, USA
7INGV, Osservatorio Vesuviano, Naples, Italy
8 Department of Physics and Astronomy, Ghent University, Ghent, Belgium
)

###### Abstract

The MURAVES experiment aims at the muographic imaging of the internal
structure of the summit of Mt. Vesuvius, exploiting muons produced by cosmic
rays. Though presently quiescent, the volcano carries a dramatic hazard in its
highly populated surroundings. The challenging measurement of the rock density
distribution in its summit by muography, in conjunction with data from other
geophysical techniques, can help the modeling of possible eruptive dynamics.
The MURAVES apparatus consists of an array of three independent and identical
muon trackers, with a total sensitive area of 3 square meters. In each
tracker, a sequence of 4 XY tracking planes made of plastic scintillators is
complemented by a 60 cm thick lead wall inserted between the two downstream
planes to improve rejection of background from low energy muons. The apparatus
is currently acquiring data. Preliminary results from the analysis of a first
data sample are presented.

##  1 Introduction

Mt. Vesuvius is an explosive volcano, still active though quiescent since
about 80 years and located near the city of Naples. Possible eruptions would
be particularly dangerous due to the dense urbanization even in its close
surroundings. The history of Vesuvius is characterized by several violent
eruptions. The best known is the Plinian eruption of 79 b. C. that destroyed
the city of Pompeii and other settlements around it.
In the course of millennia the morphology of the volcano strongly changed due
to its activity. In particular, several explosive eruptions in the last period
of activity caused the collapse of the caldera and the formation of the
present crater. The interior of the cone is thought to be composed of
materials having different densities, probably distributed in a layered
structure.
So far muon radiography (or ”muography”) was applied to volcanoes requiring
muon penetration through a rock thickness below 1 km. However, exploring the
mass distribution inside the summit of Vesuvius’ cone roughly requires
doubling the extent of the muon range, which implies coping with very low muon
rates and unexperienced backgrounds. This is the challenge of the MURAVES
(MUon RAdiography of VESuvius) project that is described in the following,
concluding with the presentation of some preliminary results.

##  2 Technology and Infrastructure

Figure 1: One of the three independent muon trackers.

The MURAVES muon telescope has a total sensitive cross-sectional area of 3 m2,
provided by an array of three identical muon trackers. Each muon tracker is
formed by a sequence of four tracking stations distributed over a length of 2
meters, with a 60 cm thick lead wall inserted between the two downstream
tracking stations for the purpose of background reduction, in particular
tagging and rejecting low energy muons. The tracking stations, realized on the
basis of the  Mu-Ray technology [1] [2] [3], consist of two adjacent planar
arrays of plastic scintillator bars, orthogonally oriented to provide the
horizontal and vertical coordinates of the muon impact point. Figure 1 shows
one of the three muon trackers. The bars have a cross-section shaped as an
isosceles triangle with 1.7 cm height and 3.3 cm basis (Figure 2) and are
assembled with alternate up-down orientation of the triangle. The space
resolution of the corresponding coordinate is improved through the weighted
average of the signal in adjacent bars. The plastic scintillator bars were
extruded at the FERMILAB-NICADD facility from a bulk of polystyrene with the
addition of PPO and POPOP scintillating dopants, emitting in a blue wavelength
spectrum centered on a 420 nm wavelength. They are of the same type as those
used in the D0 [4] and Minerva [5] experiments. Multi-clad Kuraray Y11 S-35
wavelength shifting (WLS) fibers with a 1.2 mm diameter run in a co-extruded
central hole of 1.5 $\pm$ 0.1 mm diameter. Their absorption and emission
spectra are in the blue wavelength range 400-470 nm and in the green range
470-550 nm, respectively. The WLS fibers are read-out by ASD-RGB1C-P Silicon
Photo-Multipliers (SiPM), produced by the Advansid company [6].

Figure 2: Top: The 32 scintillator bars and WLS fibers assembled in a half of
a planar array. Bottom: triangular section of the scintillator bars.

Figure 3: The temperature control system: the components (left); the system
installed on the detector (right).

Each muon tracker has 8 planar arrays of plastic scintillator bars. The 32
SiPMs of each half array (16 half-arrays for each muon tracker) are hosted on
a hybrid circuit and are connected to a 32-channel front-end electronics (FEE)
based on the ASIC chip EASIROC [7]. The setting of the trigger and the data
acquisition of each muon tracker is provided by 16 slave boards managed by a
master board. The readout and data acquisition system is described in [8]. The
experimental apparatus and the data transfer are remotely controlled.
As they are solid state devices, the SiPMs operation is affected by
temperature variations. A custom designed temperature control system based on
Peltier cells cools or heats the photosensors according to the requirement.
For each group of 32 SiPMs it consists of two Peltier cells, a copper strip
providing an uniform thermal conduction and a couple of fans for heat
dissipation. A custom circuit controls the cells through Pulse Width
Modulation drivers. The temperature can be set within 5-7 ∘C maximum
difference with respect to the ambient temperature. The bias voltage and the
temperature determine the working point of the SiPMs. The stability of the
SiPMs performance is thus obtained by changing the bias voltage according to
the temperature, so as to maintain a constant overvoltage. The temperature
control system is shown in Figure 3. Anyhow, the SiPMs must be operated at a
temperatures resulting in a sufficiently low power consumption and far enough
from the dew point, to prevent damage from condensation.
The MURAVES telescope is housed inside a container, installed on the south-
west slope of Mt. Vesuvius at 600 m a.s.l. and 1500 m distance from the
summit. The electric power is supplied by a solar panel system located on the
roof of the container and connected to an array of batteries, which ensures
continuity at night or cloudy days. Four concrete platforms inside the
container rest directly on the ground and can support lead walls of up to 90
cm thickness (at present 60 cm). One of the platforms is intended to take
calibration data with trackers in turn pointing to free sky opposite to the
mountain. Figure 4 shows the container and a scheme of the muon trackers’
arrangement inside the container.

Figure 4: Left: The container installed at Mt. Vesuvius. Right: Scheme of the
arrangement of the muon trackers inside the container: the not-colored tracker
indicates the position of a muon tracker while taking calibration data.

##  3 Measurement and sensitivity

The experiment aims at mapping the mean density of the matter crossed by muons
in the traversal of the volcano, through the measurement of the muon flux that
reaches the detector. Its ratio with the muon flux measured in calibration
runs with the muon tracker pointing to open sky in the opposite direction
(figure 5) gives the muon transmission:

$T(\theta,\phi)=\frac{N_{\mu}^{v}(\theta,\phi)/\Delta
t^{v}}{N_{\mu}^{fs}(\theta,\phi)/\Delta t^{fs}}=\frac{\epsilon^{v}\cdot
S_{eff}(\theta,\phi)\int_{E_{min(\rho)}}^{\infty}\Phi(\theta,\phi;E)dE}{\epsilon^{fs}\cdot
S_{eff}(\theta,\phi)\int_{E_{0}}^{\infty}\Phi(\theta,\phi;E)dE}$ (1)

where $\Phi(\theta,\phi;E)$ is the differential muon flux,
$N_{\mu}^{v}(\theta,\phi)$ and $N_{\mu}^{fs}(\theta,\phi$ are the muons
observed in a time interval $\Delta t^{v}$ and $\Delta t^{fs}$ respectively,
with the apices $v$ and $fs$ indicating volcano and free-sky data.
$E_{min(\rho)}$ is the energy needed to survive the rock and to be seen by the
detector; $E_{0}$ is the minimum energy needed to not be absorbed in the
detector itself. The measured fluxes can be factorized as indicated on the
r.h.s. so that with a sufficiently good approximation the efficiency
$\epsilon$ and the effective sensitive area $S_{eff}$ of the muon tracker
cancel out, isolating the ratio of the muon fluxes $\Phi(\theta,\phi;E)$ as a
function of the elevation and azimuthal angles $\theta$ and $\phi$.

Figure 5: Scheme of the configuration for the measurement of the muon flux
transmitted through the volcano (top figure) and of the free-sky calibration
flux (bottom figure). Detector not in scale.

Figure 6 shows the rock thickness to be traversed by muons, as evaluated
exploiting a Digital Elevation Model (DEM) of Mt. Vesuvius [9]. The rock
thickness ranges from some hundreds meters at the uppermost part of the cone
to almost 5000 meters at its basis. The sensitivity of the experiment depends
on the statistics which can be accumulated and on the background which is left
in spite of the tools for its reduction.

Figure 6: Map of the rock thickness to be traversed by muons through Mt. Vesuvius as a function of the elevation and horizontal angles. The horizontal range reflects the detector’s geometrical acceptance.  $\Delta$y |  $\Delta$x |  $\bar{L}$ = 500 m |  $\bar{L}$ = 1000 m |  $\bar{L}$ = 3000 m
---|---|---|---|---
9 m | 9 m | 8 months | 3 years | 100 years
9 m | 26 m | 3 months | 1 year | 33 years
9 m | 130 m | 15 days | 2.5 months | 6 years
26 m | 130 m | 5 days | 1 month | 2 years
52 m | 260 m | 2 days | 6 days | 16 months
Table 1: Exposure times expected to be required to measure the mean density
with a 10 % statistical uncertainty for a set of values of the muon path
length in the rock and of horizontal and vertical space resolutions $\Delta$x
and $\Delta$y.

The expected muon flux through the volcano was evaluated using simulations.
Table 1 gives estimates of the exposure times required to measure the mean
rock density $\rho$ with a statistical uncertainty $\Delta\rho/\rho=10\%$, for
a set of volumes defined by $\Delta x$ and $\Delta y$ and the rock thicknesses
$L$, averaged to $\bar{L}$. The required exposure times have been obtained as
described in [10]. The nominal density has been taken to be $\rho$ = 2.65 g
cm-3. In a first approximation, data are assumed to be background free and,
presumably of a lesser importance, the detector efficiency were not taken into
account. Looking at the thickness map of Mt. Vesuvius in Figure 6 in the light
of Table 1 one can deduce that a detector of a few m2 area like MURAVES can
face the challenge of exploring the summit of the cone of Mt. Vesuvius. A
larger area is required for a path length in the 3 km range.

##  4 Datasets and data analysis

Data adequate for a preliminary analysis was acquired by two of the three muon
trackers of the MURAVES telescope, named NERO and ROSSO. Because of the
variations in the environmental conditions, the SiPMs were operated at two
different working points corresponding to temperatures 15∘C and 20∘C.
Correspondingly, data is subdivided in the four groups shown in Table 2. The
table gives the respective effective exposures. The exposures pointing to the
volcano quoted in the table have an order of magnitude of 1-2 months, to be
compared to a few years of data taking foreseen for MURAVES, with the full
muon telescope. The table gives also the effective exposure already available
for each free-sky calibration of the datasets. The free-sky exposures are
shorter, but thanks to the high muon flux they give a comparable or even lower
statistical uncertainty.

Dataset | Vesuvius | Free-sky
---|---|---
ROSSO wp 15∘C | 51 days | 9.5 days
ROSSO wp 20∘C | 40 days | 14.3 days
NERO wp 15∘C | 43 days | 10 days
NERO wp 20∘C | 26 days | 17 days
Table 2: Exposure times of the four datasets with the two muon trackers
pointing to Mt. Vesuvius or to free-sky for calibration.

A preliminary analysis was performed separately for the four datasets. After a
selection of the good muon events, the data analysis proceeded by dividing the
angular range corresponding to the upper part of the cone into regions large
enough to contain a reasonable statistics, as shown in Figure 7. These regions
were chosen to be larger in azimuthal angle than in elevation, so as to obtain
a better visibility of a possible stratification of materials of different
densities. For a better balance of the rates, the layers cover an increasing
range of elevation angles going down from the top of the volcano. The three
layers were further subdivided in the left and right parts as shown in the
figure (west and a east side, respectively) in order to show possible
asymmetries in principle visible from the detector location.

Figure 7: Left: Map of the muon path length through the volcano, with
indication of the angular regions where the average density is at present in
the process to be evaluated. Right: map of the Vesuvius area with the
indication of the detector orientation.

The number of events registered in each of the six chosen angular regions are
listed in Table 3 separately for each of the four datasets. The statistical
uncertainty is maintained below 10% in most cases, in a few cases is larger
but not exceeding 20%.

| N events | Stat. unc. (%)
---|---|---
Dataset | left | right | left | right
Layer 1
ROSSO wp 15 ∘C | 428 | 439 | 0.05 | 0.05
ROSSO wp 20 ∘C | 346 | 323 | 0.05 | 0.05
NERO wp 15 ∘C | 231 | 258 | 0.05 | 0.05
NERO wp 20 ∘C | 128 | 258 | 0.07 | 0.06
Layer 2
ROSSO wp 15 ∘C | 164 | 140 | 0.08 | 0.08
ROSSO wp 20 ∘C | 106 | 109 | 0.10 | 0.10
NERO wp 15 ∘C | 78 | 79 | 0.11 | 0.11
NERO wp 20 ∘C | 61 | 63 | 0.13 | 0.13
Layer 3
ROSSO wp 15 ∘C | 61 | 76 | 0.12 | 0.11
ROSSO wp 20 ∘C | 58 | 63 | 0.13 | 0.13
NERO wp 15 ∘C | 47 | 47 | 0.14 | 0.15
NERO wp 20 ∘C | 27 | 30 | 0.19 | 0.18
Table 3: Events and statistical uncertainties in the six angular regions,
separately for the four datasets.

The expected muon flux through the volcano was evaluated using the PUMAS
software [11], a backward Montecarlo simulation program especially thought for
muography applications. The muon absorption in the volcano and the lead wall
energy threshold were taken into account. Figure 8 compares the expected muon
flux to the measured flux given by one of the datasets, to be taken as an
example. Within uncertainties there is agreement. The simulation chain is
still under development and a higher level of details is being included. In
particular, the detector responses and the effects of analysis cuts for
background reduction, such as those on the scattering in the detector and on
the muon time-of-flight for forward/backward discrimination. The expected
transmission will be evaluated with a dedicated software chain, that is
currently under construction. Finally the density will be measured in each
angular region by comparing the measured transmission to the expected
transmission.

Figure 8: Expected Flux, evaluated with PUMAS (left); Measured Flux obtained
from one of the datasets in Table 2.

##  5 Conclusions

The MURAVES experiment exploits muography to investigate, in combination with
other techniques, the rock density distribution inside the summit of Mt.
Vesuvius, a very dangerous active volcano in the Southern Italy. A muon
telescope of 3 m2 total sensitive area, subdivided into three identical and
independent muon trackers, were installed on the south-west slope of the
volcano inside a container. Four initial datasets are being analyzed, each
corresponding to exposures ranging from 26 to 51 days by single muon trackers.
Within uncertainties the measured muon agrees with the muon flux expected from
a simulation that accounts for the muon absorption in the path through the
volcano and for the energy threshold set by the detector. Calibration
datasets, taken in the same conditions, exist for each of the four datasets. A
first measurement of the mean density of the rock in six angular regions
projected on the volcano is under way. Data taking is foreseen to continue for
a few years.

## References

  * [1] G. Ambrosi et al. Nuclear Inst. and Methods in Physics Research, A 628, 120-123 (2011).
  * [2] A. Anastasio et al. Nuclear Inst. and Methods in Physics Research, A 732, 423-426 (2013).
  * [3] F. Ambrosino et al. JINST 9, C02029 (2014)
  * [4] P. Baringer et al. Nuclear Instruments and Methods in Physics Research Section A 469, p. 295 (2001)
  * [5] http://minerva.fnal.gov/
  * [6] http://www.advansid.com/
  * [7] S. Callier et al Physics Procedia 37, (2012)
  * [8] L. Cimmino et al. Annals of Geophysics, 60, 0104 (2017)
  * [9] INGV-OV 16 June 2006 Digital Terrain Model of the Campania region at 1 : 5000 scale, educational licence by Campania region n∘ 2
  * [10] N. Lesparre et al. Geophys. J. Int. 183, 1248-1361 (2010)
  * [11] V. Niess et al. Computer Physics Communications 229, 54-57 (2018)