#  Discovery of faint double-peak H$\alpha$ emission in the halo of low
redshift galaxies

J. Sánchez Almeida Instituto de Astrofísica de Canarias, La Laguna, Tenerife,
E-38200, Spain  Departamento de Astrofísica, Universidad de La Laguna, Spain
J. Calhau Instituto de Astrofísica de Canarias, La Laguna, Tenerife, E-38200,
Spain  Departamento de Astrofísica, Universidad de La Laguna, Spain  C. Muñoz-
Tuñón Instituto de Astrofísica de Canarias, La Laguna, Tenerife, E-38200,
Spain  Departamento de Astrofísica, Universidad de La Laguna, Spain  A. L.
González-Morán Instituto de Astrofísica de Canarias, La Laguna, Tenerife,
E-38200, Spain  Departamento de Astrofísica, Universidad de La Laguna, Spain
J. M. Rodríguez-Espinosa Instituto de Astrofísica de Canarias, La Laguna,
Tenerife, E-38200, Spain  Departamento de Astrofísica, Universidad de La
Laguna, Spain  Instituto de Astrofísica de Andalucía, Granada, E-18008, Spain

###### Abstract

Aiming at the detection of cosmological gas being accreted onto galaxies of
the local Universe, we examined the H$\alpha$ emission in the halo of 164
galaxies in the field of view of the Multi-Unit Spectroscopic Explorer Wide
survey (MUSE-Wide) with observable H$\alpha$ (redshift $<0.42$). An exhaustive
screening of the corresponding H$\alpha$ images led us to select 118 reliable
H$\alpha$ emitting gas clouds. The signals are faint, with a surface
brightness of ${\rm 10^{-17.3\pm 0.3}\,erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$.
Through statistical tests and other arguments, we ruled out that they are
created by instrumental artifacts, telluric line residuals, or high redshift
interlopers. Around 38 % of the time, the H$\alpha$ line profile shows a
double peak with the drop in intensity at the rest-frame of the central
galaxy, and with a typical peak-to-peak separation of the order of $\pm
200\,{\rm km\,s^{-1}}$. Most line emission clumps are spatially unresolved.
The mass of emitting gas is estimated to be between one and $10^{-3}$ times
the stellar mass of the central galaxy. The signals are not isotropically
distributed; their azimuth tends to be aligned with the major axis of the
corresponding galaxy. The distances to the central galaxies are not random
either. The counts drop at a distance $>50$ galaxy radii, which roughly
corresponds to the virial radius of the central galaxy. We explore several
physical scenarios to explain this H$\alpha$ emission, among which accretion
disks around rogue intermediate mass black holes fit the observations best.

Extragalactic astronomy (506) — Circumgalactic medium (1879) — Galaxy
accretion (575) — Galaxy fountains (596) — Galaxy winds (626) — Galaxy
formation (595) — Intergalactic medium (813) — Intergalactic gas (812) —
Astrophysical black holes (98) – Intermediate-mass black holes (816)

††journal: ApJ††facilities: MUSE@VLT (Bacon et al., 2014)††software: AstroPy
(Astropy Collaboration et al., 2013, 2018), NumPy (Harris et al., 2020),
Montage (Jacob et al., 2010)

##  1 Introduction

According to theoretical considerations and current numerical simulations, the
accretion of gas from the cosmic web is the fundamental driver of star
formation in galaxies (e.g., Dekel et al., 2009; Silk & Mamon, 2012; Fox &
Davé, 2017; Péroux & Howk, 2020). The process is particularly important in
low-mass galaxies (dark matter halo mass $<10^{12}\,{\rm M}_{\odot}$; e.g.,
Dekel et al.2013), hence, in the early Universe when most galaxies had low
masses, but also in the present Universe around isolated dwarfs. The
preference for low-mass galaxies has to do with the physics of the
cosmological accretion. When the gas encounters a massive halo
($>10^{12}\,{\rm M}_{\odot}$), it becomes shock heated and requires long time
to cool down and settle. However, gas streams reach the inner halo of low-mass
galaxies directly, in a process called cold-flow accretion (Birnboim & Dekel,
2003). The cosmological gas is predicted to be tenuous, patchy, partly
ionized, multi-phase, and of very low metallicity (e.g., van de Voort &
Schaye, 2012), and it becomes intertwined with gas from the outflows produced
by stellar feedback, i.e., by the collective effect of stellar winds and SN
(super nova) explosions from short-lived stars. Inflow and outflow rates are
often comparable, which implies that the gas ending up in stars corresponds to
only a small fraction of the gas involved in large-scale gas flows (e.g., Davé
et al., 2012; Shen et al., 2012). Shocks are to be expected either when gas
accreted from the inter-galactic medium (IGM) meets the gas already existing
in the galaxy halo, or when SN driven outflows encounter the gas in the halo
(see Sánchez Almeida et al., 2014, and references therein).

The above picture is as clear in simulations as it has been difficult to
confirm observationally. The search for cosmic gas inflows has only had
partial and/or indirect success (see Sancisi et al., 2008; Sánchez Almeida et
al., 2014; Sánchez Almeida, 2017, and references therein). For example, some
of the Ly$\alpha$ and metallic line absorptions on quasar spectra seem to be
created by intervening gas in the IGM (see, e.g., Fumagalli et al., 2011;
Péroux et al., 2019). Likewise, chemical inhomogeneities in disk galaxies
appear to be the result of pristine gas accretion (Cresci et al., 2010;
Sánchez Almeida et al., 2013, 2015, 2018; Hwang & et al., 2018; Sánchez-
Menguiano et al., 2019; Scholz-Díaz et al., 2021). Cosmic web gas may explain
some of the HI filaments found in blind HI surveys (Popping & Braun, 2011), as
well as some of the high velocity clouds detected around our galaxy (e.g.,
Sancisi et al., 2008). Finally, Ly$\alpha$ emission has been detected in
cosmic web filaments (Martin et al., 2014; Cantalupo et al., 2014; Wisotzki et
al., 2016; Bacon et al., 2021), around quasars (Arrigoni Battaia et al.,
2016), and such emission becomes pervasive when going faint enough, so that
any line-of-sight intercepts emitting gas in between redshift 6 and us
(Wisotzki et al., 2018). All these findings are rather indirect, in the sense
that connecting them with the predictions of numerical simulations still
requires a significant dose of interpretation. Therefore, it is important to
develop alternative yet complementary approaches to detect and characterize
gas flowing around galaxies.

Our work is aimed at exploring a new way to detect and characterize this gas
in the local Universe through its emission in H$\alpha$. As we mentioned
above, there is much evidence for Ly$\alpha$ emission around high-redshift
galaxies tracing gas in the IGM and CGM (circum-galactic medium). The
Ly$\alpha$ emission of nearby galaxies is not accessible from ground-based
observations but, fortunately, most physical processes generating Ly$\alpha$
photons produce H$\alpha$ as well (details in App. A). The expected H$\alpha$
flux is very low, typically $\lesssim\,10^{-18}{\rm
erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$ (App. A; Table 2), but this flux limit is
within reach of the IFU Multi-Unit Spectroscopic Explorer (MUSE; Bacon et al.,
2014) in one hour integration (Bacon et al., 2015). We take advantage of the
existing MUSE-Wide ancillary data set (Urrutia et al., 2019), which has the
required depth and field-of-view (FOV; see Sect. 2) to look for faint
H$\alpha$ signals around low redshift galaxies. The FOV has to be significant
since the virial radius of one of these galaxies, where the transition between
CGM and IGM occurs, is often of the order of one arcmin around its center
(further details in Sect. 3.1). In addition, a large FOV also provides enough
central galaxies to have a proper statistics.

The paper is organized as follows: Sect. 2 describes the MUSE-Wide data set
employed in the work. We present the selected central galaxies (Sect. 2.1) and
analyze the effectiveness of the telluric line correction (Sect. 2.2). The
data analysis is detailed in Sect. 3. These days when AI-assisted data
analysis is all the rage, our signal detection algorithm was guided by human
intuition (Sect. 3.4). The reliability of the detection was assesses later on
using statistical tests. We did not have enough previous knowledge on the
expected signals to be able to devise an automated algorithm. Various sub-
sections are devoted to explain the searching area (Sect. 3.1), the extraction
of images from the original data-cubes (Sects. 3.2 and 3.3), and to discard
various potential biases (Sect. 3.5). The outcome of the search is given in
Sect. 4, where we put forward the main physical properties of the detected
H$\alpha$ emiting clumps; their number density (Sect. 4.1), line shapes (Sect.
4.2), redshifts, fluxes and sizes (Sect. 4.3), masses (Sect. 4.4), and spatial
distribution (Sects. 4.5 and 4.6). A cross-match with external catalogs is
described in Sect. 5. Section 6 is devoted to explore various physical
scenarios to explain the observation. Finally, the main results are summarized
and discussed in Sect. 7. The appendixes are devoted to estimate the magnitude
of the expected H$\alpha$ signals (App. A), the contamination by interlopers
(App. B), the shape of the residuals left by insufficient telluric line
correction (App. C), and the search for X-ray counterparts by stacking (App.
D).

Whenever it was required, we adopted the cosmological parameters
$H_{0}=70\,{\rm km\,s^{-1}\,Mpc^{-1}}$, $\Omega_{M}$ = 0.3, $\Omega_{\Lambda}$
= 0.7. None of the results presented in the paper depends significantly on
this assumption.

##  2 Description of the MUSE Wide dataset

Figure 1: Spatial distribution of all H$\alpha$ emitting candidates classified
as level 2 and 3. These are the candidates that proceed to the second stage of
scrutiny, having their spectra inspected for possible telluric contamination.
The outline illustrates the footprint of the MUSE-Wide survey in its current
state. The color code is explained in the inset.

The Multi-Unit Spectroscopic Explorer (MUSE; Bacon et al., 2014) is a second
generation instrument on the Very Large Telescope (VLT) designed for integral
field spectroscopy in the optical band (4750 Å – 9350 Å). In its wide field
mode, it has a $\rm 1\arcmin\times 1\arcmin$ FOV with a spatial sampling of
0.2″. The spectral sampling is 1.25 Å, resulting in a resolution of around
$\rm 2.5\,\AA$. It is capable of taking around 90000 spectra in a single
exposure. We direct the interested reader to Bacon et al. (2014, 2015) for an
in-depth description of the instrument.

The MUSE-Wide survey is a blind spectroscopic survey encompassing the
CANDELS/GOODS-S and CANDELS/COSMOS regions. We include the outline of the
current state of the MUSE-Wide survey in Fig. 1. In its complete form, it will
cover ${\rm 100\,arcmin^{2}}$ of the sky and will provide spectroscopic
information for thousands of galaxies. It complements the already existing
MUSE-Deep survey in the Hubble Deep Field South (Bacon et al., 2017) by
enlarging the search area at the cost of increasing noise level. Each pointing
has an integration time of 1 hours. The current release (DR1) comprises 44
${\rm 1\arcmin\times 1\arcmin}$ contiguous fields (Fig. 1) and includes
catalogues of detected galaxies with information on emission lines, stellar
masses, and redshifts. We direct the interested reader to the data release
paper by Urrutia et al. (2019).

###  2.1 Selection of target galaxies

Our target galaxies were selected for H$\alpha$ to appear within the MUSE
wavelength range. Thus, all galaxies in the MUSE-Wide DR1 catalog with
redshift $<0.42$ were chosen to search for their H$\alpha$ emission, resulting
in a total of 164 possible host galaxies with an average redshift of $\rm
0.29\pm 0.08$. Here and throughout the paper, error bars refer to the standard
deviation of the named quantity.

The properties of these galaxies are obtained by cross-matching, within 1″,
the MUSE-Wide Survey DR1 catalog with the multi-wavelength catalogue from Guo
et al. (2013). This allows us to obtain information on the radius of the
galaxies as well as on their HST fluxes in the WFC3 filters including UV,
optical, and infrared data. We also cross-match the MUSE-Wide catalog with van
der Wel et al. (2012) to obtain parameters like Sérsic index and orientation.
Neither Guo et al. (2013) nor van der Wel et al. (2012) completely overlap
with the MUSE-Wide catalog. The match includes around 80 % of the galaxies. We
attempt to recover as many galaxies as possible by cross-matching them with
other available catalogues but, when this is not possible, the missing data
are excluded from the analysis.

The galaxies selected for scrutiny have an average stellar mass ($M_{\star}$)
of $\rm 10^{8.2\pm 0.9}\,M_{\odot}$ with the most massive galaxies reaching
$\rm\sim 10^{12}\,M_{\odot}$. They present a Sérsic index of $\rm\sim 1.5\pm
1.2$, about half-way between the expected values for low-mass spirals ($\sim
1$) and massive ellipticals ($\gtrsim 4$). Their average radius and dispersion
are $\rm 5.2\pm 3.0\,kpc$. The actual values of these three parametrers are
distributed as shown in Fig. 2 (the black lines).

Figure 2:  Global physical parameters of the central galaxies considered in
our study. Top panel: distribution of stellar masses of the parent sample (the
black line), of those galaxies having H$\alpha$ emitting clumps (the red
line), and of those without emitting clumps (the blue line). All histograms
are normalized to area equals one. Middle panel: same color code as in the
upper panel but for galaxy radii. Bottom panel: same color code as in the
upper panel but for Sérsic indexes. References for the used Sérsic indexes,
masses and radii are given in the main text, Sect. 2.1.

###  2.2 Sky emission subtraction

Molecules in the terrestrial atmosphere, particularly ${\rm O_{2}}$ and ${\rm
H_{2}O}$, absorb and emit light on their own, with many lines in the red and
near-infrared. Ground based observations have therefore to be decontaminated
from these telluric lines. This correction is particularly important in our
case since we are looking for weak emission lines in the observed spectra. The
sky subtraction is handled by MUSE-Wide using principal component analysis
(PCA) through the Zurich Atmospheric Purge algorithm (zap \- Soto et al.,
2016). zap is particularly effective in our case, when the astronomical
sources are weak and small. Regions in the science data-cubes devoid of
astronomical signal are used to reconstruct the shape of the telluric line
spectrum to be removed. The shape is the same but the amplitude differs from
one spaxel to another. There is no need for a list of telluric lines to carry
out the correction, and both weak and strong emission lines are taken into
account simultaneously and self-consistently. However, the procedure has the
inherent risk of the model sky spectrum also including faint, extended, true
astronomical signals, which will be erroneously removed from the science
spectra. zap minimizes this confusion problem through spectral filtering and
by segmentation of the sky spectrum at particular wavelengths (for details,
see Soto et al., 2016). This reduces or eliminates the number of astronomical
signals included in the sky eigenspectra.

While zap does a very good job correcting for sky lines, there is still the
possibility of contamination at low-noise level, which is where we expect the
signals to be. This residual contamination is particularly deceitful since
averaging signal-free regions of the cube does not result in the residuals
being enhanced, as one would normally expect from a stacking process.
Therefore, we cannot depend on stacking to identify residuals. Nonetheless,
the existence of sky line residuals does not seem to be affecting the
conclusions of the work. The detected H$\alpha$ signals do not overlap with
significant sky emission lines (Sect. 3.5) and the shape of the detected
signals does not agree with the shape of the residuals left by zap (App. C).
Additional arguments are collected in Sect. 3.5.

##  3 Detection of Faint H$\alpha$ Emission

This section presents the steps taken to find and select faint H$\alpha$
emission around the galaxies introduced in Sect. 2.1. The searching area
around the galaxies is defined in Sect. 3.1. Section 3.2 explains how we
merged the original data-cubes as required to carry out the actual search. The
search is based on the inspection of H$\alpha$ and broad-band images produced
out of the MUSE-Wide data-cubes (Sect. 3.3). Candidates are selected as
features with excess H$\alpha$ signal and no broad-band counterpart. The
spectrum of each candidate is then inspected individually to identify bona-
fide H$\alpha$ emission clumps (Sect. 3.4). Finally, we devote Sect. 3.5 to
argue that neither artifacts from the reduction pipeline nor contamination
from telluric lines explain the detected signals, indicating that they must
have an astronomical origin.

###  3.1 Searching area

In order to identify structures with faint H$\alpha$ emission in the CGM of
the selected galaxies, we define the searching area as a disk having 100 times
the radius of the host galaxy, taken to be the radius containing 80 % of the
light measured in the HST F160W band (1.54 $\mu m$ central wavelength; Guo et
al., 2013). The transition between the CGM and the IGM is thought to occur at
distances of approximately 70% to 90% of the optical effective radius (e.g.,
Kravtsov, 2013). Therefore, searching an area equivalent to $100\,\times$ the
IR radius allows us to examine both the full CGM and the IGM closest to the
central galaxy.

While using such a large searching radius increases the chances of detection,
it complicates the exploration since the searching region is often larger than
the area covered by a single MUSE data-cube ($1\arcmin\times 1\arcmin$). Even
if the 44 MUSE-Wide data-cubes cover a contiguous FOV (Fig. 1), each one of
them has different spatial and wavelength samplings and have to be
interpolated onto a common scale to be used together. In principle, one can
either construct the H$\alpha$ and continuum images for the individual data-
cubes and merge them later on, or merge the original data-cubes and build the
images from the resulting merged super cube. We take the second approach, as
detailed in Sect. 3.2.

###  3.2 Preparation of super data-cubes

From the original MUSE-Wide data-cubes, we build new data-cubes large enough
to encompass the searching area around each central galaxies (Sect. 3.1).
These super cubes are used not only to search for H$\alpha$ emission, but also
for all subsequent analysis. The approach of building super cubes ensures that
the artifacts introduced by merging the original data-cubes is the same
throughout the analysis, making it easier to control potential biases.

The original MUSE-Wide cubes were observed at different heliocentric velocity,
so we are forced not only to match them in RA and DEC but also to interpolate
the cubes into a common wavelength grid. We define the new common grid as that
of the central cube in each super cube. In order to assemble the super cubes,
we use the standard software Montage
111http://montage.ipac.caltech.edu/index.html, which is a toolkit for
assembling FITS images and cubes into custom mosaics. We use it in the so-
called fast reprojection. Montage has the option to correct the background of
the cubes. While this makes visualization easier, it also carries the drawback
of altering the fluxes of the faintest noise-level emission. We made sure to
assemble the data-cubes with this feature disabled, thus keeping the original
background.

While merging all MUSE-Wide cubes into a single super cube would solve our
problem, it has the disadvantage of requiring prohibitively large amounts of
memory. As a compromise, we build three different super cubes, which together
encompass the totality of the MUSE-Wide FOV and guarantee the searching area
of each galaxy to be included in at least one of them. Cutouts of these cubes
centered in the galaxies are then employed to create the H$\alpha$ and broad-
band images employed to search for weak H$\alpha$ signals, as explained in the
next subsection.

Figure 3: Images derived from the MUSE-Wide data-cubes and involved in the
process of selecting line emission candidates. The upper panels show the
H$\alpha$ image and the broad band images in $g$ and $r$, as indicated in the
inset. H$\alpha$ emission candidates must have a distinguishable signal in
H$\alpha$ and no clear emission in the broad bands. The yellow circle marks
the candidate detected in this particular region, with its radius giving the
size assigned to it. The lower panel shows the candidate in context, with the
yellow circle indicating its position with respect to the corresponding
central galaxy (the green circle). 5″ scales are included in all panels.

###  3.3 Building the H$\alpha$ and broad band images

We aim at identifying faint extended H$\alpha$ emission features within the
searching radius around the target galaxies. For this, we need to find
structures showing up in H$\alpha$ images but without broad band counterpart.
Using the MUSE-Wide super data-cubes, we construct the required narrow
H$\alpha$ images tuned to the redshift of the central galaxy, as well as
pseudo broad band $r$ and $g$ images. As an additional piece of information,
we also built color images following the recipe by Lupton et al. (2004).

We compute the H$\alpha$ image by summing the data-cube within a wavelength
range of $\rm\pm\,10\,\AA$ around the H$\alpha$ wavelength of the central
galaxy. This range represents a trade off: it is narrow enough to avoid
including too much noise into the image while, at the same time, it allows the
H$\alpha$ emission to be shifted from the central galaxy with a generous
proper motion of $\pm\,400\,{\rm km\,s^{-1}}$. We make no attempt to remove
the continuum emission from the H$\alpha$ image. A similar approach was taken
with the broad bands, except that their bandpasses were defined following the
central wavelength and full width half-maximum (FWHM) of the SDSS $g$ and $r$
filters (Fukugita et al., 1996). We add up the data-cube in a range of
$\rm\pm\,FWHM/2$ around the central wavelength of the broad band filter.
Examples of H$\alpha$ and broad band images are given in Fig. 3.

The broad band images thus constructed allow us to reach a limiting surface
brightness of 27.6 and 27.2 mag arcsec-2 for the $g$ and $r$ bands,
respectively. These values correspond to three times the RMS (root mean
square) fluctuation of the total flux per unit area considering regions of the
data-cube without sources. Following Trujillo & Fliri (2016), the areas used
for the calculation were squares of $10\arcsec\times 10\arcsec$. The same
procedure applied to the H$\alpha$ image yields a limiting surface brightness
flux of $\rm 2.4\times 10^{-19}\,erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}$. These
limits are in good agreement with estimates from the literature based on MUSE
data. Fumagalli et al. (2016) obtain a $\rm 2\sigma$ limiting magnitude of
26.6 mag arcsec-2 in the $r$ band, comparing areas of 1 arcsec radius with an
integration time of 4.1 h. Using that same band and area, we find a $\rm
2\sigma$ limiting magnitude of $\rm 25.6\,mag\,arcsec^{-2}$. Assuming the
signal-to-noise ratio to scale with the square root of the integration time,
our limiting magnitude after 4.1 h should be around $\rm
26.4\,mag\,arcsec^{-2}$, in fair agreement with Fumagalli et al. (2016). We
also compare our results with Bacon et al. (2021), which analyze 120 – 140 h
of integration in the MUSE Deep-Field. MUSE-Wide has an integration time of 1
h for each data-cube so, considering the square root scaling of the signal-to-
noise ratio with integration time, this translates into our limiting surface
brightness fluxes being expected to be 10 to 12 times worse than those of
Bacon et al. (2021). As the standard to characterize noise, they use five
times the RMS fluctuations of the signal in areas of 1 arcsec2 averaging the
spectra at 7000 Å over 3.75 Å. Their detection limit turns out to be $\rm
1.3\times 10^{-19}\,erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}$. We replicate their
estimate using MUSE-Wide data to obtain a 5 $\sigma$ limiting flux of $\rm
1.4\times 10^{-18}\,erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}$, meaning that our
limiting flux is $\sim 11$ times worse than that in Bacon et al. (2021), thus
right within the expected range.

Figure 4: H$\alpha$ line profiles for several of the candidates selected by
inspection of the H$\alpha$ and broad band images. They illustrate the
differences between levels 0, 1, 2, and 3, that go from unreliable (0) to very
reliable (3). There are examples of single peak and multiple peak H$\alpha$
profiles. The upper central and upper right plots illustrate a top-hat feature
found in some of the candidates. These candidates are never classified as
higher than level 1. The red markers on the top and bottom of the plots
indicate the wavelength range covered by our narrow H$\alpha$ image (see Sect.
3.3). The vertical blue lines mark the expected position of H$\alpha$ and the
two [Nii] lines at the redshift of the central galaxy. Note that the vertical
scale has to be multiplied by $10^{-18}$.

###  3.4 Selection of H$\alpha$ emission clump candidates, assessment, and
classification

We use the images described in Sect. 3.3 to look for faint H$\alpha$ emitting
structures while excluding sources that present a broad band contribution.
Imposing lack of broad band emission allows us to discard potential
contaminants like obvious stars and background galaxies. Unfortunately, it
does not exclude line emitter interlopers (like Ly$\alpha$ emitters). We
discuss the problem of possible astronomical contaminants in Sect. 6.

The search for H$\alpha$ emitting candidates was conducted by visual
inspection of the narrow and broad band images simultaneously (Fig. 3 gives an
example). While in principle there was the possibility of designing an
automated algorithm for the search, it was discarded because the
characteristics of the H$\alpha$ emission to be identified became known only a
posteriori. Thus, we carried out a visual inspection of the H$\alpha$ and
broad band images for each host galaxy, searching for emission in the
H$\alpha$ image and selecting those candidates that did not have a counterpart
in $r$ and $g$. Sizes were estimated visually, assuming a circular geometry
for the candidates. The process typically required around 40 – 60 minutes per
galaxy amounting to around 20 8-hour days of visual inspection in order to
cover the total 164 host galaxies within the redshift range. The procedure led
to the identification of 621 potential candidates to be further inspected.

The criteria to select H$\alpha$ clumps are somewhat arbitrary, however this
fact does not downgrade our study, which is confessedly exploratory. It aims
at detecting so far undiscovered faint H$\alpha$ emission around galaxies.
This exploratory selection will allow us to characterize the properties of
these objects so that future searches can be automated using more unambiguous
criteria. In addition, all the H$\alpha$ emission candidates go through
additional screening before they are regarded as true signals and further
statistical tests on the detections support their reliability. We classify
them according to their degree of certainty through a two stage inspection, or
two passes. They are based on the visual inspection of the spectrum resulting
from averaging all the spaxels of the candidate in the spectral regions where
the emission lines H$\alpha$, H$\beta$, [Oiii]$\lambda$5007,
[Oiii]$\lambda$4959, and [Oii]$\lambda$3727 are expected. The shape and signal
around H$\alpha$ is the main driver to assign one out of 4 grades or levels,
as detailed below:

  * -

Level 0: Immediately discarded as it cannot be a valid candidate. This
generally happens when no H$\alpha$ line can be detected, which is the
absolute minimum for a candidate to be considered.

  * -

Level 1: Probably not real but retain the possibility of revisiting it. In
general, this happens when the H$\alpha$ line is not very clear, difficult to
separate from background, and there are no other lines present.

  * -

Level 2: Probably real. H$\alpha$ line appears clearly and/or there may be
suggestions of other lines as well, such as [Oiii]$\lambda$5007\. Suspected
Ly$\alpha$ emitting interlopers are also included here, as they are candidates
we definitely want to revisit.

  * -

Level 3: Real. H$\alpha$ is clearly present along with possible presence of
the other emission lines expected from emitting gas.

Figure 4 shows examples of profiles at various levels. From the original 621
candidates, we end up with 288 level 0, 167 level 1, 155 level 2, and 11 level
3.

Only candidates classified as levels 2 and 3 are selected for a second pass.
The spectra of these candidates are further inspected to see whether or not
the emission lines, particularly H$\alpha$, could be contaminated by sky line
residuals. The inspection is carried out by over-plotting the telluric lines
compiled by Hanuschik (2003) (wavelength and intensity) on the observed
spectra. If a candidate shows no clear presence of a sky line on the vicinity
of the emission line, then the possibility of a sky artifact is minimal and
the candidate is immediately accepted. Also accepted are candidates where the
spectra present telluric lines near the emission line, but there is not
serious threat of contamination since nearby much stronger telluric lines do
not show any noticeable residual in the observed spectrum. If the emission
coincides with a strong sky line and we cannot discard serious contamination,
the candidate is classified as unclear and discarded from the list of viable
candidates to be considered for further analysis. From the original 166 level
2 and 3 candidates, we ended up with 118 good candidates. We keep the level 2
and 3 classification along the manuscript but consider only the 118 candidates
validated through the second inspection. Further details on the statistical
properties are given in Sect. 4.1. In particular, the S/N at the maximum
signal is typically larger than 3 and reaches up to 20. The same figure also
holds for the S/N of the integrated line profile (i.e., for the flux).

This second inspection was blind, without knowing what the previous
classification (level 2 or 3) was, in order to minimize biases on our part. In
general, the second pass was mainly focused on H$\alpha$, as having this line
was essential for a candidate to be picked up.

We note that our selection of good candidates is quite conservative. The
overlapping of a telluric line with the putative H$\alpha$ signal does not
imply the existence of contamination. The telluric line removal procedure
(zap) does an excellent job removing strong emission lines (see Sect. 2.2 for
details).

###  3.5 Discarding contamination due to data-cube construction artifacts and
sky emission

After the identification of candidates, a concern is the possibility of
getting false positives due to artifacts arising from the merging of data-
cubes (Sect. 3.2). Limiting these artifacts was the rationale leading to
merging entire data-cubes rather than separate images (Sect. 3.2). In addition
to this concern, there is also the possibility of contamination by artifacts
from the slicer-stacker that built the original cubes themselves (Bacon et
al., 2014; Urrutia et al., 2019). We discard these two sources of artifacts by
inspecting the spatial distribution of the selected candidates (Fig. 1). They
do not follow the boundaries between merged data-cubes nor do they portray a
discernible grid-like distribution, an indication that the majority of the
sources cannot be due to artifacts resulting from the piecing together of the
individual spectra that form the data-cubes.

When it comes to sky emission, we consider the possibility of contamination
from sky line residuals that have survived the zap treatment (Sects. 2.2 and
3.4). The shape of the observed H$\alpha$ profiles (Fig. 4, to be discussed
later on in Sect. 4.2) does not present the wiggles expected if it were due to
residuals from zap, and we take this as one evidence for our H$\alpha$
detections not being artifacts due to telluric contamination (see App. C for
further details).

A second argument against the signals being caused by telluric line
contamination is the small size of the emission features. Their typical size
is around $\rm 1\arcsec$ (see Sect. 4.3 and Fig. 13). Most telluric emission
comes from the mesopause (see, e.g. van Rhijn, 1921; Noll et al., 2014), at a
height of around 100 km. Thus, for the observed emission to be produced in
there, the physical size of the emitting gas would have a characteristic
length scale of $\rm\sim 0.5\,m$, which is too small for any physically
meaningful scale at this atmospheric layer. For example, the radial width of
the mesopause is at least several km, and a typical wind of 10 m s-1 would
move a 0.5 m patch outside the MUSE FOV in just a few seconds. These order of
magnitude arguments are fully consistent with the typical extension of the
skyline emission measured to be around 1′ (e.g., Rodrigues et al., 2010a, b),
which is much too large to produce the small observed clumps.

A third sanity check was carried out comparing the observed emission line
signals and the wavelength of intense telluric lines as compiled by Hanuschik
(2003). Very often there is not overlap, as we explain in Sect. 3.4.

Finally, an extra argument against the telluric origin of the emission line
signals comes from their spatial distribution. The observed galactocentric
distances and azimuths are inconsistent with a random distribution of sources
in the FOV. Rather, they cluster toward the major axis of the central galaxy
(Sects. 4.5 and 4.6). This fact disfavors any non-astronomical origin for
them, including sky emission.

##  4 Results

Here we present the properties of the 118 good H$\alpha$ emission candidates
selected in Sect. 3.4. It starts by discussing general statistical properties
(Sects. 4.1) and the shape of the line profiles (Sect. 4.2). Fluxes, sizes,
and redshifts are presented in Sect. 4.3. In Section 4.4, we estimate the mass
of emitting gas based on the H$\alpha$ luminosity and some assumptions on the
physical properties of the gas. Sections 4.5 and 4.6 analyze the spatial
distributions of our candidates: their distance and orientation with respect
to the major axis of the central galaxy. Kolmogorov-Smirnov tests (KS tests;
see, e.g, Press et al., 1986) on various of these distributions support the
reliability of the detection.

###  4.1 Number of H$\alpha$ emitting blobs

We have examined the outskirts of 164 galaxies with H$\alpha$ within the
bandpass of MUSE. This exercise resulted in 621 candidates being selected as
possible H$\alpha$ emitting gas clouds, with an average of 3 to 4 candidates
per galaxy. After extracting 1D spectra for our candidates and applying the
first-pass classification (Sect. 3.4), the number of valid candidates
decreases considerably. Only 27 % (166 candidates) were classified as level 2
(probably real) or level 3 (secure) and were automatically selected for the
second-pass classification to discard telluric line contamination. Most of
them passed the test (118 candidates), but the rest were put aside to be
conservative because there was a potential telluric contaminant coinciding
with the emission feature.

Overall, the above numbers translate into roughly one secure candidate per
host galaxy being found on average. Before applying the second inspection,
around 87 % of the sample had between one and two candidates per host galaxy,
with only a few outliers ($\sim 13\%$) having a higher count. After passing
the two selections, a total of 56 galaxies (from the original 164, equivalent
to $\rm\sim 34\%$) are shown to host clumps consistent with our selection
criteria. Thus, a typical galaxy with clumps has 2 of them ($118/56\simeq
2.1$).

###  4.2 Line profiles

The spectra in Fig. 4 show examples of H$\alpha$ profiles, labeled according
the classification from Sect. 3.4. From level 0 and 1, considered unlikely to
be real, to levels 2 and 3, which we consider likely candidates. One
remarkable characteristic is the appearance of line profiles presenting double
peaks (Fig. 4, middle row, first panel). A few of the candidates ($\rm\sim
2\,\%$ of the original total sample) also show no clear line on the H$\alpha$
range but instead present a kind of top-hat shape, causing the candidate to
show up in the H$\alpha$ narrow band images. One such profile is shown in Fig.
4, top right panel. Some of the double peaks present a boosted blue wing,
which goes against what we expect from high redshift Ly$\alpha$ emitter
interlopers (see Sect. 6).

Figure 5: Example of a candidate with hints of H$\beta$ emission. Both
H$\alpha$ and H$\beta$ are redshifted by ${\rm\sim 90\,km\,s^{-1}}$ with
respect to their position in the central galaxy, which is marked by the
unlabeled vertical lines. The ratio H$\beta$ to H$\alpha$ is close to $1/3$,
as expected from Case B recombination.

Figure 6: Example of a candidate with [Oiii]$\lambda$5007 emission. The lines
are clearly blueshifted with respect to the central galaxy (marked by
unlabeled vertical lines). A possible H$\beta$ emission is at noise level, but
complies with the expected $\rm\sim 1/3$ flux ratio with respect to H$\alpha$.

By selection, all candidates show signal at H$\alpha$. The spectra of some of
them also suggest the presence of other characteristic emission lines. Figures
5 and 6 show two candidates with possible emission in H$\beta$ and
[Oiii]$\lambda$5007, respectively.

We approach the problem of measuring the properties and shapes of the profiles
by fitting two Gaussian functions plus a continuum to each individual
H$\alpha$ profile. Provided the fits are good, the ratio between the fluxes of
the two fitted Gaussians allows us to quantify the fraction of double peak
profiles. The top left panel in Fig. 7a displays the diagnostic diagram used
to carry out such classification. It represents a scatter plot with the flux
ratio versus the wavelength separation of the two Gaussians, color-coded with
the S/N at the maximum emission. The noise for the S/N is obtained as the
standard deviation of the spectrum in a nearby continuum window (from 6500 Å
to 6540 Å in rest-frame wavelengths). Double peak profiles are defined to be
those where the secondary component over the main component flux ratio is
$>0.5$, the peak separation is larger than the spectra resolution
($>2.5{\rm\AA}\,\equiv\,115\,{\rm km\,s^{-1}}$; see Sect. 2), and ${\rm
S/N}>3$.

Figure 7:  Properties of the two Gaussians plus continuum fit to the observed
H$\alpha$ emission line profiles. (a) Scatter plot with the flux ratio versus
the wavelength separation of the two Gaussians, color-coded with the signal to
noise ratio (S/N) at the maximum emission. Double peak profiles (the bullet
symbols) are defined as those clumps where the flux ratio is $>0.5$, the peak
separation is larger than the spectral resolution (marked by the two vertical
dashed lines), and ${\rm S/N}>3$. Those clumps that do not meet these
conditions are represented by plus symbols. The star symbols mark the location
of the two types of stacked spectra. (b) Distribution of wavelength separation
between the two Gaussians: for all clumps (the cyan line) and only those that
meet the conditions to be regarded as double peak profiles (the red line). The
black line does not show separation between peaks but difference between the
wavelength of the peak emission and H$\alpha$. Wavelength differences are
expressed in km s-1. (c) Distribution of flux ratios of the two Gaussians: for
all clumps (the cyan line) and only those that meet the conditions to be
regarded as double peak profiles (the red line). (d) Distribution of S/N for
the peak emission (the green line) and the wavelength integrated flux of the
emission line (the blue line). All histograms are normalized to area equals
one.

The double peak line profiles thus selected make up around 38 % of the
candidates (45/118). The rest of our candidates (73/118) are regarded as
single peak emission. The properties of the resulting two Gaussian fits are
summarized in Fig. 7b (distribution of wavelength separations), Fig. 7c
(distribution of flux ratios). The fits leading to Fig. 7 were carried out
fitting a single Gaussian around the maximum emission, removing this Gaussian
from the observed profiles, and then fitting a second Gaussian to the
residuals.

Figure 7d shows the distributions of S/N ratios at the peak emission and
integrated along the line profile (i.e., the flux). As we explain above, the
noise for this estimate was computed as the standard deviation of the spectrum
in a nearby continuum window to the red of H$\alpha$. In the case of the
integrated flux, the noise is found by propagating this error into the
integral using the error propagation equation (e.g., Martin, 1971). We find
that 97 % of profiles have ${\rm S/N}>3$ in either peak or flux.

As a byproduct to initialize the first Gaussian fit, we had to compute the
centroid of the maximum emission within the H$\alpha$ detection band pass. The
distribution of wavelengths corresponding to this centroid is represented in
Fig. 7b, the black line, and it provides further evidence for the detected
signals to be real. Noise is expected to fake signals at random wavelengths.
Thus, the wavelength of the peak emission produced by noise should follow a
uniform distribution within the detection band pass. However, the observed
distribution presents two conspicuous peaks (Fig. 7b, the black line). One can
use the KS test (e.g., Press et al., 1986) to assess the null hypothesis that
the wavelengths of the maximum emission are drawn from the uniform
distribution expected from noise. This exercise provides a p-value of only
0.026, thus discarding the null hypothesis with a probability of 97.4 %. Other
similar KS tests described in Sects. 4.5 and 4.6 also discard the null
hypothesis with even higher confidence.

Figure 8: Stacking all spectra of the good candidates. Spectra around
H$\alpha$, H$\beta$, and [Oiii]$\lambda$5007 obtained by stacking assuming the
redshift of the central galaxy (top row), the redshift of the maximum of the
H$\alpha$ emission feature (central row), and using the redshift of the
central galaxy but weighting each spectrum according to the noise level in the
continuum next to H$\alpha$. The red ticks on the top and bottom edges of the
H$\alpha$ panels mark the wavelength range used to search for H$\alpha$. Each
row shares a common vertical scale.

In order to increase the signal-to-noise ratio in the emission lines, we
stacked the spectra of the 118 candidates that passed the second inspection.
The spectrum of each candidate has to be shifted in wavelength and re-sampled
to a common laboratory wavelength scale, adopting the original MUSE sampling
of 1.25 Å. We create three different stacks differing in the way the
individual spectra are shifted and averaged. (1) The first stack assumes the
redshift of the candidate to be the same as the redshift of the central
galaxy. This ignores proper motions and should give us a good general idea of
the typical spectrum of a good candidate. All spectra are weighted equally. We
apply a plain averaging for this stacking and do not normalize the spectra,
which biases the final spectrum towards the strongest signals. (2) The second
stacking assumes the redshift of the candidate to be set by the wavelength of
the peak emission inside the wavelength range used for detection and
classification (observed H$\alpha$ $\pm\,10\,$Å; Sect. 3.3). This procedure
partly corrects for the relative proper motion of the H$\alpha$ emission clump
with respect to the central galaxy, thus enhancing the resulting emission line
signals. The weighting is the same as for (1). (3) This stacking uses a
redshift equals to that of the central galaxy but weights each individual
spectrum according to the noise of the continuum next to H$\alpha$ (from 6500
Å to 6540 Å). Figure 8 shows the average spectrum obtained following these
three approaches for H$\alpha$, H$\beta$, and [Oiii]$\lambda$5007\. The
stacked spectra remains noisy for all lines except H$\alpha$, which gets
significantly enhanced. Note how the H$\alpha$ stacked spectra have a small
but noticeable continuum.

The stacking using only the redshift of the central galaxy (Fig. 8, top row)
does not appear to show any feature at H$\beta$ or [Oiii]$\lambda$5007\. While
this could be readily interpreted as absence of these lines in the majority of
the stacked spectra, we note that the noise-weighted stack shows a dip with a
small emission core at the position [Oiii]$\rm\lambda$5007 (Fig. 8, bottom
right panel).

The ratio of H$\alpha$ to H$\beta$ fluxes is expected to be around 3 if it is
produced by recombination in a fully ionized medium at some $10^{4}\,$K. We
measure the ratio to be larger than approximately 6 (Fig. 8). Thus, should the
lack of H$\beta$ would be due to dust reddening, the $V$ band
extinction222Assuming a Cardelli et al. (1989) law with a MW-like extinction.
has to be $A_{V}>2.2$, which is quite significant. This high extinction would
also explain the absence of [Oiii]$\lambda$5007\. Alternatively, one could
also reproduce the large observed flux ratio if part the emission is due to
collisional excitation rather than to recombination (e.g., Raga et al., 2015).
These conditions may be met in shocks (e.g., Tseliakhovich et al., 2012),
which also produce emission where [Oiii]$\lambda$5007 is negligibly small
compared to H$\alpha$ (e.g., Alarie & Morisset, 2019). The dense environment
existing in the broad line region of AGNs also disfavors the collisionally
excited [Oiii]$\lambda\lambda$4959,5007 lines with respect to the Balmer lines
(e.g., Osterbrock, 1984; Ho, 2008).

To further characterize the emission, we also fit two Gaussians plus a
continuum to the H$\alpha$ stacks showing double peaks (Fig 8, top and bottom
left panels). In both cases, we find the two Gaussians to be separated by 8 –
9 Å, equivalent to a velocity of $\rm 350-400\,km\,s^{-1}$, and implying a
proper motion with respect to the central galaxy of the order of 150 – 200
$\rm km\,s^{-1}$. The FWHM of the Gaussians is 6.9 Å (315 km s-1) for the
plain stack and 4.6 Å (210 km s-1) for the noise-weighted stack. Despite the
uncertainties, they are much larger than the spectral resolution of MUSE (2.5
Å), implying the existence of a large velocity dispersion within the emitting
gas. The continua are $7.5\times 10^{-20}$ and $\rm 9.2\times
10^{-21}\,erg\,s^{-1}\,cm^{-2}\,\AA$ for the plain stack and the noise-
weighted stack, respectively.

The line [Nii]$\lambda$6583 is never detected. The flux ratio between
[Nii]$\lambda$6583 and H$\alpha$ is a proxy for the metallicity of the gas
(Denicoló et al., 2002), thus, the non-detection can be transformed to an
upper limit on the metallicity. Using the calibration by Pettini & Pagel
(2004) and the same flux ratio limit employed above (i.e.,
[Nii]$\lambda$6583/H$\alpha<1/6$), we get ${\rm 12+\log(O/H)<8.5}$. This limit
is not very restrictive, though. Using the solar metallicity measured by
Asplund et al. (2009), it corresponds to around 60 % of the solar metallicity.

One final remark is in order. Some of the single peak profiles may be double
peaks with one of the two peaks missing. This is suggested by a number of
issues. Most of the profiles contributing to the stacked spectra are single
peak and yet the stacks show double peaks. The distribution of wavelengths at
the maximum emission of all profiles show two peaks (Fig. 7b, the black solid
line). Finally, the stacking of only single peak spectra also results in a
double peak profile (see Fig. 9).

Figure 9:  Stacked H$\alpha$ spectra considering all 118 candidates (blue
line), only spectra classified as double peaks (orange line) and only
spectrasified as single peak (green line). Note how the mean spectra of single
peaks also shows two peaks. The H$\alpha$ wavelength and the wavelength window
used for detection are indicated as explained in the inset. The averages are
plain means as those corresponding to the top panel in Fig. 8.

###  4.3 Fluxes, luminosities, redshifts, and sizes of the H$\alpha$ emitting
clumps

Figure 10: Distribution of log flux (top panel) and log surface brightness
(bottom panel), given in cgs units as indicated in the labels. Our candidates
have a log flux of $\rm-17.3\pm 0.3$ and a log surface brightness of
$\rm-17.2\pm 0.2$.  Figure 11: Distribution of H$\alpha$ luminosity for the
candidates that passed both inspections. The detected clumps are relatively
faint with a mean luminosity around $\rm 10^{39.2}\,erg\,s^{-1}$ and an upper
limit at around $\rm 10^{40}\,erg\,s^{-1}$.

We estimate the flux of the H$\alpha$ emission features as the total flux of
the two Gaussian fit described in Sect. 4.2. Figure 10 shows the distribution
of observed fluxes (upper panel) and surface brightness (lower panel). The
final selected H$\alpha$ candidates present a $\log$ flux distribution with
mean and standard deviation of $\rm-17.3\pm 0.3$, with the flux given in cgs
units ($\rm erg\,s^{-1}\,cm^{-2}$). The distribution of $\log$ surface
brightness has mean and standard deviation $\rm-17.2\pm 0.2$, in units of
$\rm\,erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}$.

Figure 11 contains the distribution of H$\alpha$ luminosities, derived from
the observed fluxes ($F_{\rm H\alpha}$) as,

$L_{\rm H\alpha}=4\pi D_{L}^{2}F_{\rm H\alpha},$ (1)

where $D_{L}$ is the luminosity distance, estimated from the redshift of the
central galaxy and the cosmological parameters given in Sect. 1.

Figure 12: Distribution of redshift for the selected candidates. The average
redshift is $\rm 0.31\pm 0.07$, with the error bar giving the standard
deviation. The black line represents the distribution of redshifts of the
galaxies originally selected for inspection, and is included here for
reference.

Figure 12 shows the distributions of redshift of the H$\alpha$ clumps. The
candidates have an average redshift of $0.31\pm 0.07$, with the error giving
the standard deviation of the distribution. The distribution mimics the
distribution of redshifts of the central galaxies with a preference for
redshifts around 0.35 (see the black solid line in Fig. 12).

Figure 13: Distribution of the apparent diameters (in arcsec; top panel) and
physical diameters (in kpc; bottom panel). The apparent diameters are
estimated from the circle that encompasses the totality of the emitting clump,
even when the shape of the structure is not completely circular.

Figure 13 shows the distribution of projected sizes and intrinsic sizes of the
H$\alpha$ clumps. Apparent diameters are estimated from the circle that covers
the full emitting clump (Fig. 3), even when the shape of the structure is not
truly circular. The physical size follows from the apparent size assuming the
angular size distance inferred from the redshift. The average apparent
diameter is $0.9\arcsec\pm 0.4\arcsec$, with the error giving the standard
deviation of the distribution (Fig. 13, top panel). These apparent diameters
correspond to physical diameters of $4.3\pm 2.0$ kpc (Fig. 13, bottom panel).
The spatial resolution of MUSE-Wide is around 1″ (see Sect. 3.1 and Table 1 in
Urrutia et al., 2019), therefore, except for the larger H$\alpha$ clumps, the
measured sizes seem to be set by the MUSE spatial resolution. Consequently,
most of the emitting structures are spatially unresolved. Clumps significantly
smaller than the resolution (Fig. 13, top panel) can be explained by noise,
which prevents us from detecting the full extend of the underlying H$\alpha$
emitting structure.

Figure 14: Gas masses of the H$\alpha$ emitting clumps. The masses peak around
$\rm\sim 10^{5.5}\,M_{\odot}$, although the value depends on the assumed
electron density (10 ${\rm cm}^{-3}$ in this case; see Eq. [2]). If, rather
than a value for galactic HII regions, we use a value characteristic of the
CGM or IGM, all masses increase by a factor of $10^{3}$.

###  4.4 Mass of the H$\rm\alpha$ emitting gas

We estimate the mass of the H$\alpha$ emitting gas following the approach
described in Olmo-García et al. (2017; see also Carniani et al. 2015). The
emission is assumed to be produced by recombination of hydrogen in a fully
ionized plasma (i.e., in an HII region) having the solar composition. Under
these conditions, the mass of emitting gas, $M_{g}$, can be expressed as,

$M_{g}=3.15\times 10^{3}\,{\rm M}_{\odot}~{}t_{4}\,\frac{L_{{\rm
H}\alpha}/10^{38}\,{\rm erg\,s^{-1}}}{n_{e}/10^{2}\,{\rm cm^{-3}}},$ (2)

where $t_{4}$ is the temperature of the gas in units of $10^{4}$ K, $L_{{\rm
H}\alpha}$ stands for the total H$\alpha$ luminosity (Eq. [1]), and $n_{e}$
represents the electron density. The interested reader is directed to the work
by Olmo-García et al. (2017) for further details. Assuming values of $n_{e}$
and $t_{4}$ typical of Hii regions ($10\,{\rm cm}^{-3}$ and $t_{4}=1$,
respectively; see, e.g., Osterbrock 1974), and the luminosities derived in
Sect. 4.3, we obtain gas masses ranging from $\rm 10^{4}$ to $10^{6.5}\,{\rm
M}_{\odot}$, with a mean mass of $\rm 10^{5.7}\,{\rm M}_{\odot}$. The full
distribution of values is displayed in Fig. 14. If rather than $n_{e}$ for
galactic Hii regions we assume a value characteristic of the CGM or IGM (say,
$10^{-2}\,{\rm cm}^{-3}$; e.g., Sánchez Almeida et al., 2014) all masses have
to be increased by a factor of the order of $10^{3}$, becoming from $\rm
10^{7}$ to $10^{9.5}\,{\rm M}_{\odot}$, with mean mass of $10^{8.8}\,{\rm
M}_{\odot}$.

With HII region $n_{e}$, the masses are significantly smaller than the stellar
mass of the central galaxy (Sect. 2.1), with a typical mass ratio of the order
of a thousandth (on average, $\log[M_{g}/M_{\star}]\simeq-3.1$). If, on the
other hand, CGM – IGM $n_{e}$ are assumed, then $M_{g}$ and $M_{\star}$ are
comparable.

###  4.5 Azimuthal distribution of the H$\alpha$ emission with respect to the
central galaxy

As we explain in Sect. 1, the gas from cosmological accretion is expected to
prefer the plane of the disk of the central galaxy. On the contrary, outflows
from the galaxy are channeled along the direction perpendicular to the plane
of the disk. In order to investigate the existence of any preference, and to
find additional arguments for the association of the emission with the central
galaxy, we estimate the azimuthal distribution of the H$\alpha$ clumps with
respect to the central galaxy. The orientation of the major axis of the host
galaxies is taken from the catalog by van der Wel et al. (2012), which
provides the position angle (PA) of the galaxy with respect to the celestial
north. The PA of our candidates is measured taking the galaxy as pivot and
measuring the angle counter-clockwise from the north. The difference of PAs
gives the relative orientation, which we summarize in Fig. 15 as a polar-plot
histogram.

Figure 15: Distribution of the relative azimuth of the candidates with respect
to the major axis of the corresponding central galaxy. Angles are folded to
appear in the 1st quadrant. The favoured orientation seems to match the major
axis, suggesting an alignment with the galaxy plane.

The candidates prefer relative azimuths in the range of 0 – 30∘, although with
considerable dispersion (Fig. 15). The H$\alpha$ emitting clumps seems to be
roughly aligned with the plane of the galaxy, which has two important
implications. Firstly, the scenario where the emitting gas is pristine and
coming from cosmological accretion is favored. Secondly, the clumps are not
randomly distributed in azimuth. The clumps know of the existence of a central
galaxy, which adds on to the argument that the observed emission lines are
neither random artifacts nor created by interlopers.

In order to quantify the second argument, we carried out a KS test to evaluate
the null hypothesis that the clumps are randomly distributed in the sky and so
their azimuth is drawn from a uniform distribution. The resulting p-value,
i.e., the probability of getting the distribution in Fig. 15 if the signals
were randomly distributed in the sky, is only 0.002. Thus we can discard the
null hypothesis with 99.8 % confidence. Note that under the assumption that
artifacts are expected to be randomly distributed in the sky, we are
discarding with 99.8 % confidence that the signals are fake.

###  4.6 Distribution of distances from the H$\alpha$ emission to the central
galaxy

Figure 16:  Distribution of distances between the H$\alpha$ emission and their
central galaxy. Distances are represented in multiples of the galaxy radius
adopted in the work (i.e., radius containing 80 % of the light at $\sim
1.5\,\mu{\rm m}$). Most clumps seem to be at $\rm\sim 50$ galactic radii, with
a secondary peak at just under 100 radii. The black line with error bars
represents the distribution expected if the emitting clumps where randomly
distributed throughout the MUSE-Wide footprint, and it has been normalized to
the total number of candidates. The red band approximately indicates where the
transition between CGM and IGM is expected to be.

The histogram in Fig. 16 shows the distribution of distances from the line-
emitting clumps to their respective central galaxy. In order to plot together
different galaxies, distances are normalized to the radius of the galaxy
(Sect. 3.1). The distribution has a mean and standard deviation of $\rm 53\pm
25$ radii, with a drop toward 100 radii, which corresponds to the radius of
the searching area. Most of our candidates seem to be located at a distance of
40 – 50 radii from their respective host galaxies, with the number of objects
decreasing for lower and higher distances.

The drop of the histogram at $>50$ radii approximately coincides with the
expected virial radius of the central galaxy, i.e., the boundary between the
CGM and the IGM. We estimate the virial radius to be between 42 and 54 times
the radius adopted for the galaxies (Sect. 2)333Through abundance matching,
Kravtsov (2013) shows the virial radius to be approximately 70 – 90 times the
half-mass radius of a galaxy. The typical ratio between the half-light radius
and the radius containing 80 % of the light, used in this work, is of the
order of 0.6 in our galaxies. Thus, 70 – 90 times the half-light radius
corresponds to 42 – 54 times our 80 % light radius.. This range is indicated
as a red shaded area in Fig. 16.

The question arises as to whether the observed distribution of distances may
be due to a random distribution of sources in the FOV. If the emitting clumps
were randomly distributed, their number should scale with the available
searching area at each distance from the source. This area, averaged over the
164 scrutinized galaxies and normalized to the total number of detected
emission clumps, is shown in Fig. 16 as the solid line with error bars. The
error bars represent the scatter of this scaled searching area as inferred
from the scatter among the areas of individual galaxies. The available area
should increase linearly with distance if the searching area was a perfect
disk. However, the outermost parts of the searching area often extend beyond
the boundary of the MUSE-Wide footprint, due to host galaxies being close to
the edge. The effect of this can be seen in the model histogram of Fig. 16
(the black solid line), which grows linearly up to approximately 50 radii, and
then flattens up at larger radii. Since the observed distribution of distances
is so different from the expected distribution of random sources, we conclude
that the detected emission clumps tend to cluster around the central galaxy.
There is a significant excess of sources for distances $<50$ radii and a
dearth above this distance. The divide occurs at roughly the virial radius of
the central galaxy (the red stripe in Fig. 16).

As we did for the distribution of azimuths in Sect. 4.5, a KS test was carried
out to evaluate the null hypothesis that the clumps are randomly distributed
in the sky and so their true number density is described to the black line in
Fig. 16. This time the resulting p-value is of the order of $10^{-10}$, fully
discarding the null hypothesis. As a side-effect of these results, the fact
that the observed clumps know about the existence of a central galaxy provides
further evidence that they are neither artifacts nor created by randomly
occurring interlopers.

##  5 Cross-match with catalogs of X-ray and radio sources

In order to further characterize the selected H$\alpha$ emitting clumps that
have passed both our inspections, we cross-match their position with catalogs
of X-ray and radio sources. By selection, they should not overlap with bright
continuum sources (although see the follow up work mentioned at the end of
Sect. 7).

Following Urrutia et al. (2019), we cross-match our candidates with the
Chandra CDFS-7Ms observations. The large total integration time allows this
survey to reach depths in the soft (0.5–2 keV), hard (2–7 keV), and full
(0.5–7 keV) Chandra bands of $6.4\times 10^{-18}$, $2.7\times 10^{-17}$, and
$\rm 1.9\times 10^{-17}\,erg\,s^{-1}\,cm^{-2}$, respectively. We use 3 times
the X-ray positional accuracy from Luo et al. (2017) as a limit for the
uncertainty in the matching radius.

We find no good candidates separated from an X-ray source by less than three
times the positional accuracy in Luo et al. (2017). Two level 1 candidates are
at around 1 arcsec, but these are classified as probably not real H$\alpha$
signal by our screening (Sect. 3.4). Moreover, these X-ray sources seem to be
AGNs with redshifts of 0.5 and 1 (Luo et al., 2017), thus outside the redshift
range we targeted.

We also stacked X-ray images corresponding to all the good candidates and
those with double peak emission lines (App. D gives details). The procedure
does not yield any detection, but it poses and upper limit for the average
X-ray luminosity between $2.1\times 10^{39}$ and $\rm 4.6\times
10^{39}\,erg\,s^{-1}$ for the Chandra full band ($\rm 0.5-7.0\,keV$).

We also cross-match our catalog of H$\alpha$ clumps with the Very Large Array
(VLA) 1.4 GHz survey for the Extended Chandra Deep Field South (E-CDFS;
Kellermann et al., 2008; Miller et al., 2013). The second data release goes
down to an average depth of $\rm\sim 7\mu Jy$. Following Urrutia et al.
(2019), we use a fixed 1″ search radius for matching, as the resolution of the
radio images and of MUSE-Wide are comparable. We find no match with any radio
source.

Table 1: Possible physical origin of the detected signals

Mechanism | Surface | Number | Double | Central | H$\beta$/H$\alpha$ | No | Size | Spatial
---|---|---|---|---|---|---|---|---
| Brightness | Density | Peak | Drop |  | Continuum |  | Distribution
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9)
Accretion disks (IMBH) |  ✓ |  ✓ |  ✓ |  ✓ | ? |  ✓ |  ✓ |  ✓
Expanding bubbles (SN) |  ✓ |  ✓ |  ✓ |  ✓ | ? |  ✓ |  ✓ |  ✗
Cosmological gas |  ✓ |  ✓ | ? | ? |  ✗ |  ✓ |  ✓ |  ✓
Planetary nebulae |  ✗ |  ✓ |  ✓ |  ✓ | ? |  ✓ |  ✓ | ?
X-ray binaries |  ✗ |  ✓ |  ✓ |  ✓ | ? |  ✓ |  ✓ | ?
Shocks | ? |  ✓ |  ✓ |  ✗ | ? |  ✓ | ? |  ✓
Galaxy outflows |  ✓ |  ✓ |  ✓ |  ✗ |  ✓ |  ✓ | ? |  ✗
Tidal disruption events | ? |  ✗ |  ✓ |  ✓ | ? |  ✓ |  ✓ |  ✗
Interlopers |  ✓ |  ✗ | ? |  ✗ |  ✓ |  ✓ | ? |  ✗
Jets | ? | ? |  ✓ |  ✓ | ? | ? | ? |  ✗

Note. — The symbols ✓, ✗, and ? stand for yes, no, and don’t known,
respectively. (1) Potential explanation discussed in Sect. 6. (2) Does the
physical mechanism produce surface brightness signals in the observed range?
(3) Does it predict the observed number density of objects? (4) Does it
produce double peak spectra? (5) Does the drop in-between peaks appear at the
redshift of the central galaxy? (6) Does it account for a flux ratio
H$\beta$/H$\alpha\,<1/6$? (7) Does it lack continuum emission as observed? (8)
Are the expected structures spatially unresolved? (9) Does it reproduce the
radial and azimuthal distributions?

##  6 Possible astrophysical origin of the line emission signals

As we argue in previous sections, the emission line signals, even if weak and
unusual, have to be of astronomical origin. One can think of various processes
to produced them, but none of the ones we came up with is completely free from
trouble. This section puts forward several possibilities, discusses their
viability, and highlights which ones are more likely. The degree of agreement
with observations of the various physical mechanisms is summarized in Table 1,
where they are ordered according to the number of observables they may be able
to satisfy.

$\bullet$ Contamination from emission line background sources (interlopers).
Emission line sources in the background may leak into the band-pass used to
select the H$\alpha$ emission. Ly$\alpha$ is commonly the strongest UV line
and so the usual suspect to contaminate. Given the Ly$\alpha$ luminosity
function at the appropriate redshift, in App. B we work out the expected
number of Ly$\alpha$ emitters that may be detected in the FOV around the
observed local galaxies. This number is shown to be $\lesssim$0.03
contaminants per galaxy, which is much too small to explain the observed
numbers (Sect. 4.1). Other properties of the detected emission also disfavor
the Ly$\alpha$ contamination. Even if Ly$\alpha$ sometimes show double peaks,
this shape is exceptional and should not appear in a supposed Ly$\alpha$
average profile. Moreover, double peak Ly$\alpha$ profiles tend to have a blue
lobe significantly smaller than the red lobe (e.g., Matthee et al., 2021), a
feature which is not present in the observed profiles (Sect. 4.2; Figs. 4 and
8). Finally, the drop in intensity between the two peaks (Figs. 4, 6, and 8)
coincides with H$\alpha$ at the redshift of the central galaxy, which discards
the contribution to this feature of sources unrelated to the local galaxy.
Other line emitters (Ciii], [Oii], and H$\beta$) are also considered and
discarded in App. B. In particular, the double [Oii]$\lambda\lambda$3727,3729
has a separation insufficient to account for the split that we observe (Sect.
4.2).

Finally, Sects. 4.2, 4.5, and 4.6 conclude that the observed line emission
knows of the existence of the central galaxy, which would be extremely weird
if the typical signals were produced by sources at redshift around 5. A casual
coincidence can be discarded with high confidence.

$\bullet$ Cosmological gas accretion. The level of H$\alpha$ emission to be
expected from processes associated with cosmological gas accretion is in the
range from $10^{-17}$ to $10^{-21}~{}{\rm erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$
(App. A). Since this range is very wide, it does not pose a strong constraint,
but the magnitude of the detected surface densities is consistent with this
possibility (Fig. 10, top panel). We do not claim this to be the origin of the
signals since (a) it is unclear whether gas accretion can produce the double
peak profiles that we often observe and (b) there are hints of significant
dust obscuration (H$\beta$ is fainter than expected without reddening; Sect.
4.2), which disagrees with the belief that cosmological gas is metal poor
(Sect. 1). A precise determination of metallicities is needed to explore this
issue because, at present, the upper limit we can set is too loose (Sect.
4.2).

A property consistent with cosmological gas accretion is the preference for
the signals to appear aligned with major axis of the central galaxy, and
consequently, in the plane of the galaxy (Sect. 4.5). The distribution with
distance from the central galaxy (approximately uniform within the virial
radius and with a drop toward the IGM; Sect. 4.6) is also tenable if the
emission traces cosmological gas in the process of being accreted.

Conceptually, double-peak profiles can be produced by self-absorption in the
emitting gas cloud. However, H$\alpha$ self-absorption is very unexpected. It
should be tracing extremely dense gas pockets, in conditions similar to those
existing in the central regions of AGNs ($t_{4}=1$, $n_{e}\sim 10^{9}\,{\rm
cm}^{-1}$; e.g., Kwan & Krolik, 1979; Canfield & Puetter, 1981), so that
electron collisions are able to maintain a significant fraction of the neutral
H in the $n=2$ excited state (e.g., Smith, 1980).

Shocks involving cosmological gas are to be expected when gas accreted from
the IGM meets the gas pre-existing in the halo of the central galaxy (Sect.
1). These shocks can create double-peak line profiles. However, one of the
emission components should be at rest with respect to the central galaxy, a
feature that we do not observe. We discuss this issue further in the next item
on shocks.

$\bullet$ Shocks. They are expected in the CGM, where gas accreted from the
IGM encounters hot gas already existing in the CGM, or when outflows driven by
stellar winds, SNe, or AGNs collide with halo gas (see Sect. 1). Since shocks
involve gas having two distinct velocities with the transition occurring
through a sharp interface, they are expected to produce double peak emission
lines (e.g., Schultz et al., 2005; Fadeyev & Gillet, 2004). Indeed, these
double peaks are observed in shocks running through the atmospheres of Mira
variables (Gillet et al., 1983) or in solar flares (e.g., Kuridze et al.,
2015). The expected difference of speeds between the CGM and the infalling gas
can be pretty large, of the order of several hundred km s-1 (e.g., Bennett &
Sijacki, 2020), which is fully consistent with the line splitting we detect
(Sect. 4.2). The spatial distribution to be expected is uncertain but if
shocks occur where gas from the IGM meets CGM gas, then they should be
concentrated toward the plane of the disk of the central galaxy, as observed.

So far shocks seem to be consistent with all observed features, however, there
is an additional property difficulty to reconcile with observations: one of
the two peaks of the H$\alpha$ profile is created by gas at rest with respect
to the central galaxy, however, the two observed peaks are generally split
with respect to the velocity set by the central galaxy.

$\bullet$ Accretion disks around compact objects. Compact low-luminous massive
objects are expected to lurk the CGM of galaxies. From lone low-mass stars
(e.g., Helmi, 2020), to BHs arising from stellar evolution (e.g., Fender et
al., 2013), including remnants of Pop iii stars (e.g., Madau & Rees, 2001;
Filho & Sánchez Almeida, 2018), and primordial BHs (e.g., Carr & Hawking,
1974; Clesse & García-Bellido, 2015). The number density of these objects is
assumed to be high; for example, there are claims that primordial BHs account
for all the dark matter in the Universe (Clesse & García-Bellido, 2015). Even
if this is not the case (e.g., Carr & Kühnel, 2020), the expectations clearly
overwhelm the density required to account for all the emission signals that we
observe, of the order of one clump per central galaxy (Sect. 4.1). When these
compact objects are surrounded by accretion disks, they should emit in
H$\alpha$, with the rotation of the disk giving rise to two-horn H$\alpha$
profiles when observed with the appropriate viewing angle (Smak, 1969). This
kind of double peak H$\alpha$ emission is observed in X-ray binaries
(Grundstrom et al., 2007; Zamanov et al., 2013; Casares & Torres, 2018;
Monageng et al., 2020) and cataclysmic variables (e.g., Zolotukhin &
Chilingarian, 2011), with peak separations of up to hundreds of km s-1
(Zamanov et al., 2013). At a completely different mass scale, double peak
H$\alpha$ emission is sometimes observed in the broad line region of AGN,
where it is also supposed to trace a rotating disk (e.g., Eracleous & Halpern,
1994; Ho, 2008).

The double peak H$\alpha$ emission observed around stellar-mass compact
objects comes together with continuum emission (e.g., McSwain et al., 2010;
Zamanov et al., 2019). There are hints of continuum emission in our stacks of
observed spectra (Fig. 8), but the question arises at to whether this very
faint continuum is or not consistent with an accretion disk around a compact
object. A number of arguments show that the observed lack of continuum
emission does not discard the accretion disk scenario. Firstly, large
H$\alpha$ equivalent widths (EWs) are sometimes observed in X-ray binaries.
H$\alpha$ EWs larger than 100 Å are not rare during quiescence (Fender et al.,
2009; Casares & Torres, 2018), with maxima reaching 2000 Å (Muñoz-Darias et
al., 2016), which implies line to continum flux ratios enough to account for
our observations (Fig. 8). Secondly, H$\alpha$ photons are to be produced by
recombination of H atoms photo-ionized by the accretion disk (e.g., Matthews
et al., 2015). Under particular circumstances, the nebular emission produced
by a photo-ionizing source can emit an H$\alpha$ line with very little
underlying continuum. For example, young starbursts produce H$\alpha$ with EW
in excess of $10^{3}$ Å (Leitherer et al., 1999), and galaxy-integrated
spectra with EW of hundreds of Å are not uncommon (e.g., Morales-Luis et al.,
2011). The main physical ingredient for the photo-ionization source to produce
little continuum is presenting a hard spectrum with an ionization flux greatly
exceeding the flux in the optical. Accretion disks can easily match this
requirement since they can be extremely hot with spectra peaking in the far UV
and X-ray.

The small size of the observed sources, often spatially unresolved (Fig. 13),
also fit in the accretion disk scenario. As for the spatial distribution
(Figs. 15 and 16), it can be seamlessly accommodated within this explanation
as well. BHs are expected to be present everywhere in the galaxy halo,
however, they need gas to build the accretion disk. Getting gas around a lone
BH is certainly very unlikely, but the chances are higher where the gas
concentrates, i.e., near the galaxy and close to its disk.

The H$\alpha$ luminosity of accretion disks around stellar-mass compact
objects is observed to be between $10^{32}$ and $10^{35}$ erg s-1. This
estimate has been taken from the prototypical low-mass X-ray binary V404 Cygni
(Mata Sánchez et al., 2018), with the range embracing from the quiescent to
the active phases. Even during outbursts, when the brightness is largest, the
characteristic luminosities are orders of magnitude smaller than the ones we
observe (Fig. 11). Therefore, stellar-mass compact objects, such as the ones
resulting from stellar evolution described above, cannot be responsible for
the observed H$\alpha$ emission. However, the bolometric luminosity of an
accretion disk scales with the mass accretion rate, which itself scales with
the mass of the central object (e.g., Ferrarese, 2006; Spruit, 2010; Kocsis &
Loeb, 2014). If the bolometric luminosity and the H$\alpha$ luminosity scale
with each other (a reasonable ansatz), our signals could be produced by
compact objects with mass between $10^{3}$ and $10^{6}\,{\rm M}_{\odot}$,
i.e., in the realm of the IMBHs (Intermediate-mass BHs).444We note that the
mass of emitting gas estimated in Sect. 4.4, betweem $10^{3}$ and $10^{6}$ M⊙,
also discards stellar-mass BHs and seems to be consistent with the IMBH
interpretation. Rogue IMBHs in galaxies were theoretically predicted by Madau
& Rees (2001) to result from direct gravitational collapse during early phases
of the Universe. These objects keep growing through mergers as part of the
hierarchical formation of galaxies. Even though they are generated at large
mass over-densities, i.e., at the center of proto-galaxies, gravitational
recoil during BH mergers can kick them out to appear in the outskirts of
galaxies, far from their cradle (Bellovary et al., 2019, and references
therein). Numerical simulations (e.g., O’Leary & Loeb, 2009; Micic et al.,
2011; Rashkov & Madau, 2014) predict that tens to several thousands of IMBHs
should exist in a galaxy like the MW. IMBHs can also be formed by mergers and
associations of primordial BHs, a channel that outnumbers the production via
direct gravitational collapse (e.g., Carr et al., 2019; Carr & Kühnel, 2020).

Accretion disks around stellar-mass compact objects emit most of the radiation
in X-rays. The luminosity of X-ray binaries spans from $10^{31}$ to $10^{34}$
erg s-1 during quiescence and from $10^{35}$ to $10^{39}$ erg s-1 during
outbursts (e.g., Fender et al., 2009; Casares & Torres, 2018). As we explain
above, the H$\alpha$ EW is largest during quiescence, therefore, for the
accretion disk scenario to explain our H$\alpha$ luminosities the boosting
factor has to be around $10^{6}$. A naive scaling of the X-ray luminosities
during quiescence with this factor implies that the IMBH accretion disks
should have an X-ray luminosity $<10^{40}$ erg s-1, which renders a flux
density at the typical redshift of our sources ($\sim 0.3$; Fig. 12) of
$<3\times 10^{-17}$ erg s-1 cm-2. This luminosity and flux density are upper
limits, not only because they correspond to up-scaling the largest observed
X-ray luminosities but because an increase of the central object mass drops
the temperature of the accretion disk555The inner radius of the disk grows
with increasing central mass thus dropping the disk temperature. This drop is
quite large and cannot be compensated by the increase of temperate arising
from the growth of the central mass; see, e.g., Spruit (2010)., which moves
the bulk emission from X-rays to the far UV. An additional comment is in
order: some of the IMBH candidates in literature correspond to X-ray sources
with luminosities larger than $10^{40}$ erg s-1 (Mezcua, 2017), thus, larger
than the upper limit used above. This, once again, is consistent with the
scenario. These large X-ray luminosities may be the IMBH equivalent to the
stellar-mass BHs emission during extreme outbursts.

We finally note that the X-ray luminosities and fluxes predicted above are
close to but below the detection threshold of the current X-ray surveys
explored in Sect. 5.

$\bullet$ Outflows and galactic fountains. The process of star formation and
the AGN activity expel gas from galaxy disks. This ejected gas ends up in the
CGM and, depending on the depth of the gravitational well, also in the IGM.
Part of this material rains back to the galaxy in the so-called galactic
fountains (e.g., Fraternali, 2014). In principle, the outflows can be
distinguished from cosmological accretion because they are not metal-poor and
because they are channeled by the galaxy disk to break out in the direction
perpendicular to the plane of the disk (e.g., Sánchez Almeida et al., 2014).
The existence of outflows with these properties has been inferred in various
ways, for instance, as metal-rich absorption on background sources (e.g.,
Péroux et al., 2006; Fynbo et al., 2011), or as a bimodality in the azimuthal
angle and metallicity of the absorbing gas around galaxies (e.g., Kacprzak et
al., 2012; Lehner et al., 2013).

Galaxy driven outflows are able to account for the number and the surface
brightness of our H$\alpha$ emissions. However, three other properties of the
observed clumps are not well reproduced within this scenario. Firstly, the
observed clumps are concentrated toward the plane of the galaxy (Fig. 15)
rather than where the outflows should be, in the direction perpendicular to
this plane. Secondly, within the virial radius, the observed clumps are not
centrally concentrated but spread out. Thirdly, outflows are often supersonic
(e.g., Olmo-García et al., 2017) and can collide with the gas in the CGM, thus
creating a shock front that emits double-peak line profiles. However,
contrarily to observation, one of the two peaks has to be at rest with respect
to the central galaxy. Double peak line emission are indeed observed in
outflows near the disks of luminous IR galaxies, with one of the two
components at the systemic velocity of the galaxy (Heckman et al., 1990).

Outflows and the subsequent galactic fountains involve metal rich gas, where
the existence of dust grains is to be expected. Thus, the observed large
H$\alpha$ flux compared with H$\beta$ could, in principle, be explained as an
effect of dust reddening in outflows from galaxies (see Sect. 4.2).

$\bullet$ Expanding Bubbles. Expanding gas shells give rise to top-hat line
shapes (e.g., Tenorio-Tagle et al., 1996), which do not have a central dip.
Such line shape is at variance with the double peaks that we often observe,
however, if the shell is dusty then photons at the center of the line are
preferentially absorbed within the shell, resulting in two-horn profiles
(e.g., appendix in Olmo-García et al., 2017). In this case, the blue lobe
tends to be largest.

This type of shell producing double peak lines could be associated with
planetary nebulae (PNe), i.e., the ionized circumstellar gas ejected during
the late phases of evolution of solar-mass stars. It produces emission lines
with little continuum since the ionization is provided by a very hot central
star. However, the luminosity of one of these shells is not sufficient to
account to the typical $3\times 10^{38}$ erg s-1 H$\alpha$ luminosity that we
observe (PN luminosities span from $10^{34}$ to $5\times 10^{35}$ erg s-1;
Ciardullo, 2010). In addition, PNe emit more in [Oiii]$\lambda$5007 than in
H$\alpha$, a feature we do not see (Fig. 8). Thus, PNe can be discarded as
sources of the observed signals. An alternative shortcut to reach the same
conclusion is noting that the mass of emitting gas (Fig. 14), all
uncertainties notwithstanding, is much larger than the envelope that could be
ejected from a low-mass star.

Expanding bubbles driven by SN explosions do not have these drawbacks. SN
remnants tend to have a H$\alpha$ flux significantly larger than
[Oiii]$\lambda$5007 (e.g., Davis et al., 2018). The H$\alpha$ luminosity at
the explosion is much larger than the ones we observe, however, $\sim 10^{3}$
days later from the outset, the luminosity drops to the observed values (e.g.,
Benetti et al., 2016). Thus, the luminosity itself is not an issue, but the
fact that we do not observe SN remnants at their brightest point questions the
whole scenario. Why do we observe SNe only after a particular time lag from
the explosion? Double peak H$\alpha$ profiles are common in SNe and they
present a continuum consistent with observations (e.g., Benetti et al., 2016).
On the other hand, SN remnants are product of the stellar evolution and so
should mimic the spatial distribution of the stars in the central galaxy. This
fundamental characteristic of SNe is in disagreement with the spatial
distribution of observed H$\alpha$ signals which, in addition to be spread out
over large areas, they also appear in the CGM where the formation of the stars
that go off SN is unlikely (Fig. 16).

$\bullet$ Tidal disruption events (TDEs). When a star passes sufficiently
close to a super-massive BH (SMBH), it can be stretch out to disruption by
tidal forces (e.g., Carter & Luminet, 1982). A fraction of the original mass
is captured to form a transient accretion disk around the BH. The emission of
the system resembles the one from accretion disks discussed above, in
particular, double peak H$\alpha$ lines are produced (Short et al., 2020; Hung
et al., 2020) with a maximum luminosity significantly larger than the ones
that we observe ($>10^{41}\rm\,erg\,s^{-1}$; Holoien et al. 2019, Short et al.
2020). The luminosity and spectrum are variable so it is unclear whether the
double peak line shape and the right H$\alpha$ luminosity occur simultaneously
during the transient. Similarly, it is unclear whether a scaled down version
of the observed TDEs, with smaller central BH masses, can pull apart a star
and give rise to an accretion disk with the characteristics we observe. The
main handicap for the TDEs to explain our observations is the fact that the
two required ingredients, a star plus a SMBH, cluster toward the center of the
gravitational potential. Thus an encounter is far more likely to occur at the
center of the galaxies, which is at variance with the spatial distribution of
the H$\alpha$ signals (Fig. 16). In addition, it is also a short-lived event
($\lesssim 1$ yr), therefore, hard to reconcile with the relatively high
number of the observed signals.

$\bullet$ Jets. Astrophysical jets are quite common and they involve bipolar
outflows which, once spatially averaged, can give rise to double peak line
profiles like the ones that we often detect. They originate in accretion disks
around compact central objects such as BHs of all masses, neutron stars, X-ray
binaries, T Tauri stars, Herbig–Haro objects, or evolved post-AGB stars. Jets
around stellar mass objects can be discarded with the argument that the
H$\alpha$ emitting gas seems to involve between $10^{3}$ to $10^{6}\,{\rm
M_{\odot}}$ (Fig. 14). This narrows down the possibilities to IMBHs, and many
of the restrictions invoked above apply. In addition, the unresolved character
of many of the detected sources (Fig. 13) may be at variance with the large
scales expected in jets around massive BHs. Jets can carry away a significant
part of the energy released during the accretion process and so, from the
point of view of the energetics, they are similar to accretion disks (e.g.,
Romero, 2021).

##  7 Conclusions

Aiming at detecting cosmological gas being accreted onto galaxies of the local
Universe (Sect. 1), we examined the outskirts of all the 164 galaxies in the
MUSE-Wide footprint with H$\alpha$ within the bandpass of MUSE (redshift
$<0.42$). After a first visual inspection of H$\alpha$ maps around all
galaxies, an exhaustive screening of the putative faint emission line signals
led us to identify what seems to be 118 H$\alpha$ emitting gas clouds (Sects.
3.4 and 4.1). The signals are very faint, with a typical surface brightness of
${\rm 10^{-17.2\pm 0.2}\,erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$ and a typical
flux of ${\rm 10^{-17.3\pm 0.3}\,erg\,s^{-1}\,cm^{-2}}$ (Fig. 10). We
discarded that they are created by instrumental artifacts or telluric line
residuals (Sects. 3.5, 4.2, 4.5, and 4.6, and App. C). Neither are they high
redshift interlopers (Sect. 6, and App. B). KS tests on the spatial and
spectral distribution of the signals discard with high confidence ($>97$ %)
that the signals are fake.

Assuming that the emission line signals are produced by gas clumps around the
central galaxies, they have the following properties:

  * -

The H$\alpha$ line profile often shows a double peak with the drop in
intensity at the rest-frame of the central galaxy (Figs. 4, 6, and 8). The
typical peak-to-peak velocity is of the order of $\pm 200$ km s-1. The double
peak line profiles make up almost 38 % of all candidates (Sect. 4.2). Some of
the single peak profiles may be double peaks with one of the two peaks missing
(Fig. 9).

  * -

Only a few of the H$\alpha$ emitting clumps also show H$\beta$ (Fig. 5) and/or
[Oiii]$\lambda$5007 (Fig. 6). These emission lines are virtually absent in the
average stacked spectra (Fig. 8), so that their fluxes have to be smaller than
around 1/6 of the flux in H$\alpha$ (Sect. 4.2).

  * -

Using the non-detection of [Nii]$\lambda$6583, we set a loose upper limit on
the metallicity of the emitting gas, ${\rm 12+\log(O/H)<8.5}$, which roughly
corresponds to 60 % of the solar metallicity (Sect. 4.2).

  * -

Most of the detected line emission clumps are at the resolution limit of MUSE-
Wide, around 1″. We measure diameters of ${\rm 0.9\pm 0.4\,arcsec}$ (Fig. 10,
top panel), which at the redshift of the central galaxies ($0.31\pm 0.07$;
Fig. 12) correspond to ${\rm 4.3\pm 2.0\,kpc}$ (Fig. 10, bottom panel).

  * -

Using the distance of the central galaxy, the inferred H$\alpha$ luminosities
are around $\rm 10^{39}\,erg\,s^{-1}$ and they seldom exceed $\rm
10^{40}\,erg\,s^{-1}$ (Fig. 11).

  * -

The mass of emitting gas ($M_{g}$) has been estimated under the assumption
that the emission is produced by H recombination in a fully ionized medium.
With a temperature of ${\rm 10^{4}\,K}$ and $n_{e}={\rm 10\,cm^{-3}}$, $M_{g}$
goes from $10^{4}$ to ${\rm 10^{6.5}\,M_{\odot}}$ (Fig. 14). These masses are
significantly smaller than the stellar mass of the central galaxy
($M_{\star}$; Sect. 2.1), with an average mass ratio of the order of
$\log(M_{g}/M_{\star})\simeq-3.1$. $M_{g}$ scales as $n_{e}^{-1}$ so for
densities closer to those of the CGM and IGM, $M_{g}/M_{\star}\sim 1$.

  * -

The signals are not isotropically distributed with respect to the central
galaxy. The azimuth of the signals tend to be aligned with the major axis of
the corresponding central galaxy (Fig. 15).

  * -

The distances from the emitting clumps to the central galaxies are not
randomly distributed (Sect. 4.6). There is an excess of sources for distances
$<50$ galaxy radii and a dearth above this value (Fig. 16). The critical
distance roughly corresponds to the virial radius of the central galaxies,
which defines the transition between CGM and IGM.

  * -

We find no match between the position of the H$\alpha$ clumps and catalogues
of X-ray sources (Luo et al., 2017) or radio sources (Kellermann et al., 2008;
Miller et al., 2013); see Sect. 5.

We explore several physical mechanisms to explain these observations including
the cosmological gas accretion that motivated the study (Sect. 6). The degree
of agreement with observations provided by each one of them is summarized in
Table 1. According to number of observables that they satisfy, the main
possibilities are:

  * -

Accretions disks around rogue IMBHs. They may provide the observed surface
brightness, number density, double peak profiles centered at the redshift of
the galaxy, lack of continuum emission, unresolved sizes, and correct spatial
distribution. IMBHs may also be able to reproduce the observed small
H$\beta$/H$\alpha$ flux ratio.

  * -

Expanding bubbles caused by SN explosions. They account for all the
observables mentioned in the previous item except for the spatial distribution
(SN trace stellar sites, whereas the emission is spread over much larger
distances to reach out the IGM). It also remains unclear whether they can
account for H$\beta$/H$\alpha$.

  * -

The cosmological gas accretion naturally accounts for the spatial distribution
of the emitting clumps. They can also cope with their typical luminosities and
sizes. However, it has difficulty to explain the double peak in H$\alpha$. In
addition, the problem to explain H$\beta$/H$\alpha$ is particularly severe in
this scenario since the cosmological gas is metal-poor, and so, unable to
produce significant dust reddening.

  * -

Shocks driven either by cosmological gas accretion or by galaxy winds. In this
case, the main obstacle is explaining double peak H$\alpha$ profiles with the
intensity drop at the rest-frame of the central galaxy. Shocks with CGM gas
can create double peak line profiles, but with one of the peaks at the
systemic velocity. In addition, the spatial distribution is at odds with
galaxy driven winds, while shocks involving metal-poor cosmological gas has
difficulty to account for H$\beta$/H$\alpha$.

Thus, from all the possibilities, accretion disks around rogue IMBHs satisfy
the observational constraints best.

Note that the key discriminant among the posed physical scenarios is the need
to reproduce the double peak emission centered at the redshift of the galaxy.
Even if this feature is quite common (38 % or more; Sect. 4.2), most of the
observed spectra present a single peak and so are not subject to this
constraint. Thus, they can be produced by mechanisms other than the IMBHs,
including the cosmological gas accretion we were after (Sect. 1). The cause of
the detected emission could be assorted, and follow up studies will be
required to assign them individually.

The criteria we have employed to select H$\alpha$ clumps are rather fuzzy,
however, this fact does not undermine our study, which is confessedly
exploratory. It was aimed at detecting so far unknown faint H$\alpha$ emission
around galaxies. Its detection, presented in this paper, will allow us to
characterize its properties so that future searches can be automated, and so,
based on less ambiguous criteria. As an additional followup, we are in the
process of cross-matching the positions of the H$\alpha$ emissions with
existing ancillary data, in particular, with deep Hubble Space Telescope
multi-band images of the region (Whitaker et al., 2019). Thus, several faint
optical continuum counterparts of the emitting clumps have been found and
their physical properties will be reported soon (González-Morán et al. 2022,
in preparation).

Thanks are due to various colleagues that helped us along the way: Ana Monreal
Ibero, for discussions and references on may different aspects of the work.
Teo Muñoz Darias, for discussion on the interpretation of the double peak
H$\alpha$ signals based on stellar mass accretion disks. He also pointed out
the possibility of TDEs. Mar Mezcua and Lucio Mayer, for discussions and
references on the IMBH hypothesis. Peter Weilbacher, for clarifications on how
the removal of telluric lines is carried out, both using the standard approach
and ZAP. Antonio Cabrera, for discussion and references on the size and
timescales of the emitting patches if telluric. Ignacio Trujillo, for various
specific suggestions on cross matching with external catalogs and on the
parameters to be elucidated. Tanya Urrutia, for clarifications regarding the
MUSE data-cubes and their manipulation. Graham Berriman, for help with
Montage. We are particularly thankful to an anonymous referee for helping us
quantify the reliability of the detections. We acknowledge support from the
Spanish Ministry of Science and Innovation, project PID2019-107408GB-C43
(ESTALLIDOS), and from Gobierno de Canarias through EU FEDER funding, project
PID2020010050. This research made use of Montage, which is funded by the
National Science Foundation under Grant Number ACI-1440620, and was previously
funded by the National Aeronautics and Space Administration’s Earth Science
Technology Office, Computation Technologies Project, under Cooperative
Agreement Number NCC5-626 between NASA and the California Institute of
Technology.

##  Appendix A Expected H$\alpha$ signals

This appendix works out the H$\alpha$ emission to be expected in the CGM and
IGM around nearby galaxies. The estimates are very uncertain, but they all
converge to provide a level of emission that includes the surface brightness
of the signals detected in this work. Table 2 summarizes the expected signals
under the various circumstances and approximations explained in the next
paragraphs.

Table 2: H$\alpha$ signal to be expected in the local Universe.

Based on | Surface brightness | Region | Reference
---|---|---|---
(1) | (2) | (3) | (4)
Scaled Ly$\alpha$ observations | 1 – 30 | IGM | Cantalupo et al. (2014)
Scaled H$\alpha$ observations | 10 | IGM | Leibler et al. (2018)
Fluorescence | 0.01 – 10 | IGM | Kollmeier et al. (2010)
Fluorescence | 0.03 – 0.1 | IGM | Cantalupo et al. (2005)
CIB$a$$a$footnotemark:  Fluorescence | 0.001 – 0.1 | CGM | Bland-Hawthorn et al. (2017)
Mechanical Feedback | $\sim 1$ | CGM | Mori et al. (2004)
Gravity driven | 0.09 – 5 | CGM | Goerdt et al. (2010)
Gravity driven | 0.001 – 2 | CGM | Faucher-Giguère et al. (2010)
Gravity driven $z=0$ | 0.002 – 20 | IGM | Furlanetto et al. (2003)
Observed $z=0$ | $\geq 0.5$ | CGM | Cheung et al. (2016); Herenz et al. (2017)
Ly$\alpha$ blobs at $z=3$ | 1 – 30 | CGM | Matsuda et al. (2004)

aafootnotetext: Cosmic Ionizing Background.

Note. — (1) Origin of the empirical data or physical mechanism used to
estimate the H$\alpha$ signal (details in App. A). (2) In units of
$10^{-18}~{}{\rm erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$. (3) Whether the signals
correspond to the CGM or the IGM. (4) Main reference used for the estimate.

As we put forward in Sect. 1, diffuse Ly$\alpha$ emission around galaxies is
routinely detected at high $z$. Most Ly$\alpha$ emission mechanisms ultimately
rely on the recombination of electrons and protons, which grants the
production of H$\alpha$ photons together with the Ly$\alpha$ photons.
Detecting these H$\alpha$ photons rather than the Ly$\alpha$ photons has pros
and cons. Among the cons, the flux ratio between Ly$\alpha$ and H$\alpha$
($\epsilon$) is expected to be around 8.7 if the lines are produced by
recombination (e.g., Miller, 1974). It is observed to be between 0.5 and 15
(Baldwin, 1977; Hu, 1992; Atek et al., 2009; Hayes et al., 2010; Leibler et
al., 2018). The emissivity of the lines is expected to increase with
increasing density (e.g., Osterbrock, 1974), therefore, another con has to do
with the evolution of the cosmic web gas that decreases in density
($\rho_{g}$) with decreasing redshift ($z$) following the global expansion of
the universe (e.g., van de Voort & Schaye, 2012), i.e,

$\rho_{g}\propto(1+z)^{3}.$ (A1)

Among the pros, the faint Ly$\alpha$ signals have been detected around high
redshift galaxies whereas H$\alpha$ can be observed in local galaxies that are
closer and thus less affected by cosmological dimming. From redshift 0 to $z$,
there is a gain in surface brightness that scales as $(1+z)^{4}$ (e.g., Goerdt
et al., 2010; Hogg, 1999). A final advantage of H$\alpha$ is the fact that
Ly$\alpha$ photons are more easily destroyed than H$\alpha$ photons (e.g., by
dust absorption or collisional de-excitation) and so $\epsilon$ tends to be
reduced with respect to the theoretical recombination value. For the sake of
having a rough estimate of the expected signal, we will assume that the H
emission roughly drops with the mean gas density of the Universe. Then the
ratio between the Ly$\alpha$ surface brightness at redshift $z$ ($S\\!B_{\rm
Ly\alpha}[z]$) and the H$\alpha$ surface brightness if the same structure is
observed in the local Universe ($S\\!B_{\rm H\alpha}[0]$) would be,

$S\\!B_{\rm H\alpha}(0)=\frac{(1+z)^{4}}{\epsilon(1+z)^{3}}S\\!B_{\rm
Ly\alpha}(z)=\frac{1+z}{\epsilon}S\\!B_{\rm Ly\alpha}(z),$ (A2)

which yields

$S\\!B_{\rm H\alpha}(0)\sim S\\!B_{\rm Ly\alpha}(z),$ (A3)

provided that $z$ is in the range between 2 and 4 and $\epsilon$ remains
around 5, which is reasonable considering the range of redshifts where
extended Ly$\alpha$ emission has been detected (Sect. 1 and the references
below) and the observed values for $\epsilon$. To sum up, the two fluxes are
expected to be similar (Eq. [A3]), a result that depends on two key
assumptions, namely, (1) $\epsilon$ is not far from the value provided by
recombination and (2) the line emission scales with the mean density of the
universe.

$\bullet$ Emission expected scaling the signal found by Cantalupo et al.
(2014). The work by Cantalupo et al. is used as a reference to quantify the
Ly$\alpha$ signals observed at $z$ between 2 and 3, but similar results are
obtained using any of the other works mentioned in Sect. 1 (e.g., Martin et
al., 2014; Wisotzki et al., 2016, 2018). The peak surface brightness of the
Ly$\alpha$ structure detected by Cantalupo et al. (2014) is $S\\!B_{\rm
Ly\alpha}(3)\simeq{3\times 10^{-17}}~{}{\rm
erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$. The Ly$\alpha$ signal presents a range of
values of around a factor of 30, therefore, considering Eq. (A3) $F_{\rm
H\alpha}(0)$ is expected to be in between $1$ and $30$ in units of
$10^{-18}~{}{\rm erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$.

Leibler et al. (2018) measured the H$\alpha$ emission corresponding to the
Ly$\alpha$ nebula detected by Cantalupo et al. (2014). H$\alpha$ is in the
near infrared at the redshift of the source, so, they use the spectrograph
MOSFIRE666https://irlab.astro.ucla.edu/instruments/mosfire/ on the Keck
telescope. The H$\alpha$ surface brightness is around $3\times 10^{-18}\,{\rm
erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$, with a Ly$\alpha$ to H$\alpha$ flux ratio
of 5.5. This observational result is consistent with the above estimates. If
we correct for the variation with redshift, the H$\alpha$ signal would be
around $10^{-17}\,{\rm erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$, as listed in Table
2.

$\bullet$ Fluorescence of UV photons. UV photons from the cosmological
background or from nearby sources are expected to photo-ionize the gas clouds
existing in the CGM and IGM. The recombination of the photo-ionized H produces
Ly$\alpha$ photons that, after internal scattering, can escape producing a
diffuse Ly$\alpha$ signal that traces the cloud. In optically thick clouds,
the mechanism very effectively transforms the UV photons into Ly$\alpha$
photons, with an efficiency approaching one (e.g., Gould & Weinberg, 1996).
This whole process is often refereed to as fluorescence (e.g., Cantalupo et
al., 2005; Gould & Weinberg, 1996). Since photons are finally produced through
recombination, fluorescence produces all the other recombination lines of H,
including H$\alpha$.

There are works estimating the Ly$\alpha$ signals expected due to the
illumination by the cosmic UV background and by nearby quasars. Using
cosmological simulations of gas clouds at $z\simeq 3$, Cantalupo et al. (2005)
work out the Ly$\alpha$ signals, that turn out to be in the range between
$3\times 10^{-20}~{}{\rm erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$ and a few times
$10^{-19}~{}{\rm erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$, the latter corresponding
to the presence of a nearby quasar. Kollmeier et al. (2010) carry out a
similar exercise rendering larger signals that can reach a value of
$10^{-17}{\rm erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$ if the quasar illumination
is strong enough. The signal excited by the UV background is in the range of
$10^{-20}~{}{\rm erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$. For more references on
the estimates see, e.g., Cantalupo et al. (2012). The scaling between the
fluxes in H$\alpha$ and Ly$\alpha$ in Eq. (A3) can also be used in this case
to render the range of expected H$\alpha$ fluxes included in Table 2.

Bland-Hawthorn et al. (2017) model the H$\alpha$ emission in the outer gas
disk of MW-like galaxies produced by fluorescence of photons coming from the
cosmic ionizing background (CIB). The estimated signals are in the range
between $10^{-19}$–$10^{-21}~{}{\rm erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$.

$\bullet$ Emission driven by mechanical feedback from galaxies. The mechanical
energy injected by galaxy winds may also lead to shocks that are powerful
sources of ionizing photons produced in the hot post shock plasma as it cools
down (e.g., Dopita & Sutherland, 1995). They may be responsible for some of
the extended Ly$\alpha$ emission found at $z\sim 2$ – 3 (e.g., Taniguchi &
Shioya, 2000; Ohyama et al., 2003). Mori et al. (e.g., 2004) model the winds
produced by SN explosions in primordial galaxies at $z=3$, and work out the
emission properties assuming an optically thin gas in collisional ionization
equilibrium. The total flux in Ly$\alpha$ is around $10^{-16}~{}{\rm
erg\,sec^{-1}\,cm^{-2}}$, and it is emitted in region of around 10 arcsec,
which renders a surface flux density of some $10^{-18}~{}{\rm
erg\,sec^{-1}\,cm^{-2}\,arcsec^{-2}}$. According to Eq. (A3), this is also the
level of signal to be expected if observing the shocks in H$\alpha$ in the
local Universe.

$\bullet$ Gravity driven emission. The potential energy of the gas falling
into a gravitational well is transformed into kinetic and thermal energy, and
has to be released for the gas to become gravitationally bound. This energy is
partly radiated away in the form of H recombination lines (Dijkstra & Loeb,
2009).

There are several studies in the literature that model this process and the
resulting line emission. We take the work by Goerdt et al. (2010) as
reference. They use cosmological numerical simulations of galaxy formation to
work out the Ly$\alpha$ flux to be expected from the gas streams that fed
galaxies at high redshift. UV background excitation, collisions with free
electrons, and dust attenuation are included when computing the Ly$\alpha$
emission. They emission is mainly driven by the excitation produced by
collisions with electrons since the filaments are dense enough to be shielded
from the UV background. Most of the Ly$\alpha$ emission comes from extended
(50 – 100 kpc) narrow, partly clumpy, inflowing, cold streams of
$\sim\,10^{4}$ K that feed the growing galaxies. The predicted morphology is
irregular, with dense clumps and elongated extensions. The typical surface
brightness increases with decreasing distance from the halo centre. The peak
surface brightness is $2\times 10^{-17}$ ${\rm
erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$ for MW-like halos at $z=2.5$, with typical
signals in the range from $3\times 10^{-17}$ to $5\times 10^{-19}\,{\rm
erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$. No estimate for the H$\alpha$ flux is
given, however, they mention that other H emission lines are expected to be
two orders of magnitude less luminous in our model. H$\alpha$ is expected to
be the strongest one so we assume $\epsilon$ (i.e. the ratio Ly$\alpha$ and
H$\alpha$) to be 20. Using this value for $\epsilon$ in Eq. (A2), we end up
with the expected H$\alpha$ surface brightness range reported in Table 2.

Faucher-Giguère et al. (2010) carry out a similar simulation but with more
sophisticated radiative transfer. As in the work by Goerdt et al. (2010), the
signals are concentrated towards the center of the dark matter halos, but some
also emerge from the IGM. Depending on the treatment of the radiative
transfer, the predicted Ly$\alpha$ signals can vary a lot, going from
$10^{-17}$ to $5\times 10^{-21}~{}{\rm erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$.
Re-scaling to H$\alpha$ and to the local Universe, as we did for Goerdt et al.
signals, the expected H$\alpha$ surface brightness is in the range from
$2\times 10^{-18}$ to $10^{-21}\,{\rm erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$
(Table 2).

Furlanetto et al. (2003) worked out predictions directly for the local
Universe ($z\sim 0.15$), therefore, there is no need to apply high-$z$
corrections in this case. Their Fig. 2 shows structures producing between
$10^{2}$ to $10^{6}$ photons cm-2 s-1 sr-1, which properly transformed into
the usual units render $S\\!B_{\rm Ly\alpha}(0)$ in the range between $4\times
10^{-20}$ and $4\times 10^{-16}~{}{\rm erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$.
Correcting only for the ratio of emissivities between H$\alpha$ and
Ly$\alpha$, $S\\!B_{\rm H\alpha}(0)$ turns out to be between $2\times
10^{-21}$ and $2\times 10^{-17}~{}{\rm erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$.

Likewise, Rosdahl & Blaizot (2012) model the Ly$\alpha$ emission in numerical
simulations. The emitting gas is at the center of the dark matter halos
radiating away the excess of gravitational energy explicitly modeled in the
previous works. The expected signals are similar to the ones resulting from
the previous simulations. Updated numerical simulations by Witstok et al.
(2021) at $z>2$ give results consistent with the previous ones.

$\bullet$ Observed diffuse H$\alpha$ emission in the local Universe. The
existing H$\alpha$ surveys in the local Universe have a surface brightness
limit insufficient to detect the faint signals to be expected (Table 2): equal
or fainter than $10^{-17}~{}{\rm erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$. This
level is not within reach of the existing surveys (e.g., Moss et al., 1988;
James et al., 2004; Epinat et al., 2008; Dale et al., 2008; Jaiswal & Omar,
2016; Greimel et al., 2021). Moreover, some of them are focused on the
emission of the MW, either diffuse (e.g., Gaustad et al., 2001) or its point-
like sources (e.g. Barentsen et al., 2014).

However, diffuse H$\alpha$ emission have been detected around individual
galaxies of the local Universe when the observation is deep enough. As part of
the LARS777Ly$\alpha$ Reference Sample (Östlin et al., 2014). project to study
the visibility, strength, and escape fraction of the Ly$\alpha$ photons in 14
galaxies the local universe, Östlin et al. (2014) and Hayes et al. (2014)
obtained deep H$\alpha$ maps hinting at the existence of extended emission.
Although the sensitivity is good (around $10^{-18}~{}{\rm
erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$; Hayes et al., 2013), the FOV includes
only a few effective radii around the galaxies, therefore, missing most of the
CGM and the full IGM. Flux-wise, the integrated emission observed in H$\alpha$
and Ly$\alpha$ are comparable.

Herenz et al. (2017) detect filaments of ionised gas in the CGM of the
extremely metal-poor galaxy SBS 0335-052E. The features are detected in
H$\alpha$ and [Oiii]$\lambda$5007 down to a limiting surface brightness of
$5\times 10^{-19}~{}{\rm erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$. The galaxy
diameter is around 3 kpc, with the filaments extending more than 9 kpc and
connecting seamlessly in velocity space to the galaxy. Extended H$\alpha$
emission at a level of $5\times 10^{-18}~{}{\rm
erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$ have been found around 3 dwarf galaxies by
Lee et al. (2016). Cheung et al. (2016) find a lone H$\alpha$ emitting blob
next to an evolved galaxy with a signal around $3\times 10^{-18}~{}{\rm
erg\,s^{-1}\,cm^{-2}\,arcsec^{-2}}$. The flux representative of these
observations is included in Table 2.

##  Appendix B Contamination by background sources

Assume that emission line sources with luminosity function $\Phi(L,z)$ are
leaking into the band-pass of the filter $f(\lambda)$ used to detect
H$\alpha$,

$f(\lambda)=\Pi\Big{(}\frac{\lambda-\lambda_{obs}}{\Delta\lambda}\Big{)},$
(B1)

with $\lambda_{obs}$ and $\Delta\lambda$ representing the central wavelength
and the width of the filter, respectively. In our case, the filter has a top-
hat shape,

$\Pi(x)=\cases{1&$|x|\leq 1/2$,\cr 0&elsewhere,}$ (B2)

with $\Delta\lambda=20$ Å and $\lambda_{obs}$ the wavelength of H$\alpha$
($\lambda_{\rm H\alpha}$) at the redshift of the central galaxy $z_{g}$, i.e.,
$(1+z_{g})\,\lambda_{\rm H\alpha}$. If the contaminants are Ly$\alpha$
emitters, the window function defining the range of redshifts $z$ where they
contribute is

$W(z)=f\Big{[}(1+z)\lambda_{Ly\alpha}\Big{]}.$ (B3)

Given any window function, the total number of targets leaking in is (e.g.,
Peebles 1993, Eqs. [13.60] and [13.61]; Fried et al. 2001, Eq. [3]),

$N=\Omega\int_{0}^{\infty}\,W(z)\Big{(}\int_{Lmin(z)}^{Lmax(z)}\,\Phi(L,z)\,dL\Big{)}\,{{dV(z)}\over{dz}}\,dz,$
(B4)

where $dV(z)/dz$ represents the variation with redshift of the comoving volume
of the universe per unit solid angle, and $\Omega$ is the solid angle of the
field covered by the measurement. If $F_{min}$ is the minimum flux density
that we can observe (with $[F_{min}]=$ erg s-1 cm-2), then

$L_{min}(z)=4\,\pi\,D_{L}^{2}(z)\,F_{min},$ (B5)

with $D_{L}$ the luminosity distance (e.g., Peebles, 1993; Hogg, 1999). On the
other hand, we do not include in our sampling sources brighter than $F_{max}$,
otherwise they would be detected as individual galaxies and excluded from the
search. The value of $F_{max}$ sets $L_{max}$, explicitly,

$L_{max}(z)=4\,\pi\,D_{L}^{2}(z)\,F_{max}.$ (B6)

Since the redshift window $W(z)$ is very narrow, all other functions of $z$ in
Eq. (B4) can be regarded as constant and so can be pulled out of the integral,
rendering

$N\simeq C\,\Omega\,\frac{\Delta\lambda}{\lambda_{\rm Ly\alpha}}{{dV(z_{\rm
Ly\alpha})}\over{dz}}\,\int_{Lmin(z_{\rm Ly\alpha})}^{Lmax(z_{\rm
Ly\alpha})}\,\Phi(L,z_{\rm Ly\alpha})\,dL,$ (B7)

with

$z_{\rm Ly\alpha}=\frac{\lambda_{\rm H\alpha}}{\lambda_{\rm
Ly\alpha}}(1+z_{g})-1.$ (B8)

A global fudge factor $C$ has to be included in Eq. (B7) to account for the
completeness of detections. This factor can be far from unity when the
expected signals are close to the noise level, which is our case (e.g.,
Blanton et al., 2005; Herenz et al., 2019).

All in all, the number of Ly$\alpha$ emitters leaking in ($N$ in Eq. [B7])
depends on the luminosity function $\Phi$ at redshift $z_{\rm Ly\alpha}$, the
flux limits ($F_{min}$ and $F_{max}$), the angular extend of the field-of-view
around each central galaxy ($\Omega$), and the cosmological parameters that
define the Universe (through $dV/dz$ and $D_{L}$).

Figure 17a uses Eq. (B7) to work out the expected number of contaminants per
galaxy assuming typical values for the parameters describing our observation:
$F_{min}=10^{-17.9}$ erg sec-1 cm-2, $F_{max}=4\times F_{min}$ (Fig. 10),
$C=0.2$ (Herenz et al., 2019), and $\Omega=\pi r_{obs}^{2}$ with $r_{obs}=100$
arcsec (Sect. 3). The derivative of the comoving volume ($dV/dz$) and $D_{L}$
have been computed using the astropy cosmology package (Astropy Collaboration
et al., 2013, 2018) and the cosmological parameters adopted in the text (Sect.
1). The Ly$\alpha$ luminosity functions used in the estimate are represented
in Fig. 17b (Henry et al., 2012; Zheng et al., 2013; Sobral et al., 2018).
Figure 17b also shows the range of luminosities used to estimate the number of
contaminants (the blue band). Only the luminosity function worked out by
Sobral et al. (2018) reaches faint enough luminosities for the estimate to be
based on actual Ly$\alpha$ emitter counts (i.e., the blue solid line and the
blue band in Fig. 17b overlap). This luminosity function predicts less than
0.03 contaminants per central galaxy at $\lambda_{obs}\simeq 8500$ Å(Fig.
17a), which corresponds to $z_{g}\simeq 0.3$ where most signals appear (Fig.
12). This level of contamination is insufficient to account for the number of
observed signals, of the order of one per galaxy (Sect. 4.1).

Figure 17: (a) Expected number of background sources contaminating our
signals. We give the number of contaminants per observed galaxy as function of
the observed wavelength ($\lambda_{obs}=[1+z_{g}]\lambda_{\rm H\alpha}$). The
axis on top shows the redshift if the contaminants were Ly$\alpha$ emitters.
The dashed lines correspond to estimates based on luminosity functions
extrapolated outside the observed range of luminosities, whereas the solid
lines use the range of observed luminosities. (b) Luminosity functions used to
compute the number of contaminants in (a). They are either Schechter functions
or Saunders et al. functions based on luminosities observed within the range
shown as solid lines, with the luminosity function outside this range shown as
dashed lines. The color code is the same as in panel (a). The colored regions
indicate the range of luminosities used to estimate the number of
contaminants. As indicated in the inset, we include several high redshift
Ly$\alpha$ luminosity functions (Henry et al., 2012; Zheng et al., 2013;
Sobral et al., 2018), a Ciii]$\lambda$1909 luminosity function (Stroe et al.,
2017, at redshift 1), a [OII]$\lambda\lambda 3726,3729$ luminosity function
(Comparat et al., 2016, at redshift 0.8), and a H$\beta$ luminosity function
(Comparat et al., 2016, at redshift 0.3). As usual, rather than the actual
luminosity function $\Phi(L,z)$, the plot shows the luminosity function per
$\log L$ interval ($=\ln[10]\times L\times\Phi[L,z]$).

Note that above formalism holds for any other emission line replacing $\Phi$
and $\lambda_{\rm Ly\alpha}$ with the corresponding luminosity function and
rest-frame wavelength of emission. We expected that these other contaminants
are less important than the Ly$\alpha$ emitters because they are less numerous
for the same $L$ and because they sample a smaller redshift range (i.e., they
sample a smaller volume of Universe). To check it out, we also include in Fig.
17 the luminosity function and the expected number of counts for
Ciii]$\lambda$1909, [Oii]$\lambda\lambda 3726,3729$, and H$\beta$ emitters
(Stroe et al., 2017; Comparat et al., 2016). The largest number of expected
counts correspond to [Oii] emitters, and it is similar to the (low number of)
Ly$\alpha$ contaminants and so insufficient to account for observation
presented in the paper.

##  Appendix C Line shapes expected from partly corrected telluric lines

Figure 18: Monte Carlo simulation showing the residuals to be expected under
insufficient telluric line contamination correction. Top left panel: mock
emission lines used to construct the PCA database that reproduce them all
(black solid lines). Right panels: first 4 eigenvectors of the resulting PCA
database. Bottom left panel: residual left when the profile shown as the red
solid line in the top left panel is represented in the PCA database using
different sets of eigenvectors (as indicated in the insets). Intensities and
wavelengths are given in arbitrary units.

Telluric line contamination is one of the main concerns to interpret the
detected signals (Sect. 3.5). This appendix describes a schematic Monte Carlo
simulation carried out to show the shapes to be expected if telluric line
residuals were present in the observed spectra. The telluric line removal is
based on a principal component analysis (PCA) decomposition of the spectra in
the part of the FOV assumed to have no astronomical signal and thus assumed to
represent sky emission (Soto et al., 2016). Simply put, such sky emission
signal is used to compute the eigenvectors that characterize the telluric
contamination. Then the spectrum of each spaxel of the full FOV is decomposed
in this PCA database, and the sky emission spectrum resulting for every spaxel
is subtracted out from the observed spectrum. In a very schematic
representation to illustrate the procedure, we carried out the PCA
decomposition of $10^{4}$ mock Gaussian emission line spectra with various
strengths, widths, and shifts: the grey lines in the top left panel of Fig.
18. They mock any of the telluric emission lines present in the MUSE FOV. The
first four eigensvectors of this decomposition are shown in Fig. 18, right
panels. Changes in line width are described by odd eigensvectors whereas
changes in line shift are represented by even eigenvectors. This PCA database
is used to represent the problem profile shown as the red solid line in Fig.
18, top left panel. The residuals left when 1, 2, 3, and 4 eigenvectors are
included in the representation of the problem profile are shown in the bottom
left panel of Fig. 18. As the number of eigenvectors increases, the residuals
have more wiggles. The lobes of the residuals are positive and negative.
Importantly, the observed H$\alpha$ shapes (Sect. 4) do not have any of the
properties of these mock residuals. Although this Monte Carlo simulation is
very schematic, and in reality the PCA decomposition of the sky is based on
many emission lines, the shape of the true residuals do not vary much from our
simple modeling (see Figs. 2 and 3 in Soto et al., 2016). Thus, we use the
difference in shape between the observed H$\alpha$ profile and the expected
residuals as an additional argument that disfavor the observed signals to be
caused by telluric line contamination. Similar results and identical
conclusions are reached when using other Monte Carlo realizations and problem
profiles.

##  Appendix D X-ray stacking

The cross-matching of our sources with the available Chandra X-ray catalogue
for the CDFS field did not yield positive results (Sect. 5). In order to
improve the X-ray sensitivity, we resort to stacking techniques. We make use
of the images publicly available from Luo et al. (2017) and create stacks of
our candidates using mean and median statistics.

We start by creating cut-outs of $\rm 100\times 100\,px^{2}$ centered around
each gas candidate, then we stack the associated exposure maps to create
images with the mean and the median counts s-1 in every pixel (Fig. 19). The
signal is measured as the central count rate in the stack using an aperture of
$\sim 1.5\arcsec$ ($\sim 3$ px) radius, which is large enough to account for
the sizes of the majority of our clumps (Sect. 4.3). We estimate the noise by
randomly sampling the signal in the stacked image with the same aperture and
avoiding sources. The distribution of count rates of 2000 such random
apertures is used to compute the noise. Specifically, we use as noise three
times the difference between the mean and the 16th percentile of the
distribution, a value equivalent to the 3-$\sigma$ error bar for a Gaussian
distribution.

We stack the X-ray cutouts of our candidates using mean and median statistics
for two different groups. The first one uses all level 2 and 3 candidates
whereas the second one includes only the candidates who present evidence of
multiple peaks. We find no evidence of X-ray emission in any of the stacks
(Fig. 19 shows the result of stacking all sources with mean and median
statistics). However, we are still in a position to provide upper limits for
the X-ray flux and luminosity of our stacked candidate.

We determine the conversion factor between count rate and X-ray flux (in the
Chandra 0.5-7 keV band) using all the sources in Luo et al. (2017). The mean
of the ratio between the flux of each source and its count rate is used as
conversion factor888 $\rm(1.5\pm 0.5)\times
10^{-11}\,\,erg\,s^{-1}\,cm^{-2}\textordmasculine\,(counts\,s^{-1})^{-1}$..
Using this calibration, the flux corresponding to the 3-$\sigma$ noise level
of the mean stack considering all candidates and those with multiple peaks are
$\rm 1.18\times 10^{-17}\,erg\,s^{-1}\,cm^{-2}$ and $\rm 6.5\times
10^{-18}\,erg\,s^{-1}\,cm^{-2}$, respectively The values using median stacking
are similar. The corresponding X-ray luminosity is computed using Eq. (1),
with the luminosity distance $D_{L}$ set by the median redshift of the clumps
in the stack. The luminosity equivalent to the noise level ($L_{X}$) is of the
order of $\rm 4.6\times 10^{39}\,erg\,s^{-1}$, for the average stacking of the
full sample, and $\rm 2.1\times 10^{39}\,erg\,s^{-1}$, for the stacking using
only multiple peaks.

Figure 19:  The result of stacking our candidates in full band (0.5-7 keV) 7Ms
X-ray Chandra images. The left panel shows mean stacking whereas the right
panel shows median stacking. The mean stacking reveals the presence of some
clear sources, likely in the vicinity of several of our clumps, but no actual
detection in the position of our candidates (marked with a circle), a lack of
signal also supported by the median stack. Our mean stacking reveals an X-ray
non-detection with an upper limit luminosity of $\rm L_{X}<4.6\times
10^{39}\,erg\,s^{-1}$. The images have been scaled to their minimum (black)
and $\rm 0.85\times maximum~{}(white)$ values.

## References

  * Alarie & Morisset (2019) Alarie, A., & Morisset, C. 2019, Rev. Mexicana Astron. Astrofis., 55, 377, doi: 10.22201/ia.01851101p.2019.55.02.21
  * Arrigoni Battaia et al. (2016) Arrigoni Battaia, F., Hennawi, J. F., Cantalupo, S., & Prochaska, J. X. 2016, ApJ, 829, 3, doi: 10.3847/0004-637X/829/1/3
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481, doi: 10.1146/annurev.astro.46.060407.145222
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Atek et al. (2009) Atek, H., Kunth, D., Schaerer, D., et al. 2009, A&A, 506, L1, doi: 10.1051/0004-6361/200912787
  * Bacon et al. (2014) Bacon, R., Vernet, J., Borisova, E., et al. 2014, The Messenger, 157, 13
  * Bacon et al. (2015) Bacon, R., Brinchmann, J., Richard, J., et al. 2015, A&A, 575, A75, doi: 10.1051/0004-6361/201425419
  * Bacon et al. (2017) Bacon, R., Conseil, S., Mary, D., et al. 2017, A&A, 608, A1, doi: 10.1051/0004-6361/201730833
  * Bacon et al. (2021) Bacon, R., Mary, D., Garel, T., et al. 2021, A&A, 647, A107, doi: 10.1051/0004-6361/202039887
  * Baldwin (1977) Baldwin, J. A. 1977, MNRAS, 178, 67P, doi: 10.1093/mnras/178.1.67P
  * Barentsen et al. (2014) Barentsen, G., Farnhill, H. J., Drew, J. E., et al. 2014, MNRAS, 444, 3230, doi: 10.1093/mnras/stu1651
  * Bellovary et al. (2019) Bellovary, J., Brooks, A., Colpi, M., et al. 2019, BAAS, 51, 175.  https://arxiv.org/abs/1903.08144
  * Benetti et al. (2016) Benetti, S., Chugai, N. N., Utrobin, V. P., et al. 2016, MNRAS, 456, 3296, doi: 10.1093/mnras/stv2811
  * Bennett & Sijacki (2020) Bennett, J. S., & Sijacki, D. 2020, MNRAS, 499, 597, doi: 10.1093/mnras/staa2835
  * Birnboim & Dekel (2003) Birnboim, Y., & Dekel, A. 2003, MNRAS, 345, 349, doi: 10.1046/j.1365-8711.2003.06955.x
  * Bland-Hawthorn et al. (2017) Bland-Hawthorn, J., Maloney, P. R., Stephens, A., Zovaro, A., & Popping, A. 2017, ApJ, 849, 51, doi: 10.3847/1538-4357/aa8f45
  * Blanton et al. (2005) Blanton, M. R., Lupton, R. H., Schlegel, D. J., et al. 2005, ApJ, 631, 208, doi: 10.1086/431416
  * Canfield & Puetter (1981) Canfield, R. C., & Puetter, R. C. 1981, ApJ, 243, 390, doi: 10.1086/158606
  * Cantalupo et al. (2014) Cantalupo, S., Arrigoni-Battaia, F., Prochaska, J. X., Hennawi, J. F., & Madau, P. 2014, Nature, 506, 63, doi: 10.1038/nature12898
  * Cantalupo et al. (2012) Cantalupo, S., Lilly, S. J., & Haehnelt, M. G. 2012, MNRAS, 425, 1992, doi: 10.1111/j.1365-2966.2012.21529.x
  * Cantalupo et al. (2005) Cantalupo, S., Porciani, C., Lilly, S. J., & Miniati, F. 2005, ApJ, 628, 61, doi: 10.1086/430758
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245, doi: 10.1086/167900
  * Carniani et al. (2015) Carniani, S., Marconi, A., Maiolino, R., et al. 2015, A&A, 580, A102, doi: 10.1051/0004-6361/201526557
  * Carr et al. (2019) Carr, B., Clesse, S., Garcia-Bellido, J., & Kuhnel, F. 2019, arXiv e-prints, arXiv:1906.08217.  https://arxiv.org/abs/1906.08217
  * Carr & Kühnel (2020) Carr, B., & Kühnel, F. 2020, Annual Review of Nuclear and Particle Science, 70, 355, doi: 10.1146/annurev-nucl-050520-125911
  * Carr & Hawking (1974) Carr, B. J., & Hawking, S. W. 1974, MNRAS, 168, 399, doi: 10.1093/mnras/168.2.399
  * Carter & Luminet (1982) Carter, B., & Luminet, J. P. 1982, Nature, 296, 211, doi: 10.1038/296211a0
  * Casares & Torres (2018) Casares, J., & Torres, M. A. P. 2018, MNRAS, 481, 4372, doi: 10.1093/mnras/sty2570
  * Cheung et al. (2016) Cheung, E., Stark, D. V., Huang, S., et al. 2016, ApJ, 832, 182, doi: 10.3847/0004-637X/832/2/182
  * Ciardullo (2010) Ciardullo, R. 2010, PASA, 27, 149, doi: 10.1071/AS09022
  * Clesse & García-Bellido (2015) Clesse, S., & García-Bellido, J. 2015, Phys. Rev. D, 92, 023524, doi: 10.1103/PhysRevD.92.023524
  * Comparat et al. (2016) Comparat, J., Zhu, G., Gonzalez-Perez, V., et al. 2016, MNRAS, 461, 1076, doi: 10.1093/mnras/stw1393
  * Cresci et al. (2010) Cresci, G., Mannucci, F., Maiolino, R., et al. 2010, Nature, 467, 811, doi: 10.1038/nature09451
  * Dale et al. (2008) Dale, D. A., Barlow, R. J., Cohen, S. A., et al. 2008, AJ, 135, 1412, doi: 10.1088/0004-6256/135/4/1412
  * Davé et al. (2012) Davé, R., Finlator, K., & Oppenheimer, B. D. 2012, MNRAS, 421, 98, doi: 10.1111/j.1365-2966.2011.20148.x
  * Davis et al. (2018) Davis, B. D., Ciardullo, R., Feldmeier, J. J., & Jacoby, G. H. 2018, Research Notes of the American Astronomical Society, 2, 32, doi: 10.3847/2515-5172/aab045
  * Dekel et al. (2013) Dekel, A., Zolotov, A., Tweed, D., et al. 2013, MNRAS, 435, 999, doi: 10.1093/mnras/stt1338
  * Dekel et al. (2009) Dekel, A., Birnboim, Y., Engel, G., et al. 2009, Nature, 457, 451, doi: 10.1038/nature07648
  * Denicoló et al. (2002) Denicoló, G., Terlevich, R., & Terlevich, E. 2002, MNRAS, 330, 69, doi: 10.1046/j.1365-8711.2002.05041.x
  * Dijkstra & Loeb (2009) Dijkstra, M., & Loeb, A. 2009, MNRAS, 400, 1109, doi: 10.1111/j.1365-2966.2009.15533.x