11institutetext: IRAP, Université de Toulouse, CNRS, CNES, 9 avenue du Colonel
Roche, 31028 Toulouse, France 22institutetext: Institut d’Estudis Espacials de
Catalunya, Universitat de Barcelona (ICC-UB), Martí i Franquès 1, 08028
Barcelona, Spain
e-mail<EMAIL_ADDRESS>
# Statistical study of a large and cleaned sample of ultraluminous and
hyperluminous X-ray sources

Hugo Tranin Statistical study of a large and cleaned sample of ultraluminous
and hyperluminous X-ray sourcesStatistical study of a large and cleaned sample
of ultraluminous and hyperluminous X-ray sourcesStatistical study of a large
and cleaned sample of ultraluminous and hyperluminous X-ray sourcesStatistical
study of a large and cleaned sample of ultraluminous and hyperluminous X-ray
sources Natalie Webb and Olivier Godet Statistical study of a large and
cleaned sample of ultraluminous and hyperluminous X-ray sourcesStatistical
study of a large and cleaned sample of ultraluminous and hyperluminous X-ray
sourcesStatistical study of a large and cleaned sample of ultraluminous and
hyperluminous X-ray sourcesStatistical study of a large and cleaned sample of
ultraluminous and hyperluminous X-ray sources

###### Abstract

Context. Ultra-/hyperluminous X-ray sources (ULX and HLX) could constitute
interesting laboratories to further improve our understanding of the
supermassive black hole growth through super-Eddington accretion episodes and
successive mergers of lighter holes. ULXs are thought to be powered by super-
Eddington accretion onto stellar-mass compact objects, while HLXs are of
unknown nature, but could be good candidates for accreting intermediate mass
black holes (IMBH). However, a significant portion of the sample of ULX/HLX
candidates derived from catalogue searches are in fact background active
galactic nuclei (AGN).

Aims. Here we build samples of ULXs and HLXs from the three largest X-ray
catalogues available, compiled from XMM-Newton, Swift-XRT and Chandra
detections, and the GLADE catalogue containing 1.7 million galaxies at D¡1000
Mpc. We aim to characterise the frequency, environment, hardness and
variability of ULXs and HLXs in order to better assess their differences and
understand their populations.

Methods. After a thorough classification of these X-ray sources, we were able
to remove $\sim$42% of $S/N>3$ sources which were shown to be contaminants, to
obtain the cleanest sample of ULX/HLXs to date. From a sample of 1342 ULXs and
191 HLXs detected with a signal to noise $>3\sigma$, we study the occupation
fraction, hardness, variability, radial distribution and preferred environment
of the sources. We build their Malmquist-corrected X-ray luminosity functions
(XLF) and compare them with previous studies. Thanks to the unprecedented size
of the sample, we are able to statistically compare ULXs and HLXs and assess
the differences in their nature. The interpretation of HLXs as IMBHs is
investigated.

Results. A significant break is seen in the XLF at$\leavevmode\nobreak\ \sim
10^{40}$ erg s-1. With our ULX sample, containing $\lesssim 2$% of
contaminants, we are able to confirm that ULXs are located preferentially in
spiral galaxies and galaxies with higher star-formation rates, with an excess
seen in dwarf galaxies. While X-ray binaries (XRBs), ULXs and most HLXs share
common hardness and variability distributions, a fraction of HLXs appear
significantly softer. Unlike ULXs, HLXs seem to reside equally in spiral and
lenticular/elliptical galaxies. 35% of the HLX candidates have an optical
counterpart, and we estimate the mass of 120 of them to be in the range of
$2\times 10^{3}-10^{5}\mathrm{M}_{\odot}$. Most of the HLX population is found
to be consistent with an accreting massive black hole in a dwarf galaxy
satellite of the primary host. This diverse nature needs to be confirmed with
deeper optical and infrared observations, as well as upcoming X-ray
facilities.

###### Key Words.:

Catalogs – X-rays: general – X-rays: binaries – X-rays: galaxies – Methods:
statistical

##  1 Introduction

Supermassive black holes (SMBH, in the mass range
$10^{6}-10^{10}\mathrm{M}_{\odot}$) are found in the centres of the most
massive galaxies observed so far, and are thought to play a major role in
galaxy evolution notably during their violent episodes of accretion and
ejection (Fabian 2012; Padovani et al. 2017). However, the formation
mechanisms of these objects are yet to be fully understood. Some SMBH observed
at $z>6$, less than 1 Gyr after the Big Bang, are more massive than
$10^{9}\,\mathrm{M}_{\odot}$ (e.g. Mortlock et al. 2011; Wang et al. 2021),
implying rapid and efficient growth mechanisms that may be a mixture of
intermediate-mass black hole (IMBH, $10^{2}-10^{5}\mathrm{M}_{\odot}$) mergers
and sustained super-Eddington accretion onto these seeds (Volonteri et al.,
2008; Haiman, 2013; Pacucci & Loeb, 2022). However, compelling evidence for
the existence of IMBH is still missing, with the exception of some recent
discoveries through gravitational waves and X-ray observations (e.g. Farrell
et al. 2009; Mezcua 2017; Abbott et al. 2020).

The search and study of X-ray sources in nearby galaxies, especially in the
luminosity range $10^{39}-10^{42}$ erg s-1, could be a promising way to gain
insights into the SMBH growth mechanisms. Indeed, for an object spherically
accreting hydrogen gas, the luminosity must be lower than the Eddington limit
– $L_{Edd}=1.26\cdot 10^{38}(\frac{M}{\mathrm{M}_{\odot}})$ erg s-1 – for the
accretion flow to be sustained. It is interesting to note, however, that in
rare cases there appears to be physical mechanisms that help stabilise the
accretion flow in the super-Eddington regime, as suggested by both recent
models and observations (e.g. Inayoshi et al. 2016; Massonneau et al. 2023;
Lin et al. 2017; Belfiore et al. 2020). However, using $L<L_{Edd}$ is one way
to estimate a lower limit on the accretor mass. Extra-galactic and off-nuclear
ultraluminous X-ray sources (ULXs) with X-ray luminosities $L_{X}>10^{39}$ erg
s-1 were thus first thought to be good IMBH candidates (e.g. Colbert &
Mushotzky 1999; Liu & Bregman 2005). However, the spectral curvature seen in
ULXs, hinting at an accretion regime different than for X-ray binaries (e.g.
Bachetti et al. 2013), and the discovery of pulsating ULXs (e.g. Bachetti et
al. 2014; Fürst et al. 2016; Israel et al. 2017; Carpano et al. 2018; Quintin
et al. 2021) later revealed that some of these sources are instead neutron
stars accreting above the Eddington limit. In ULX spectra, the inner
temperature of the disc is cooler than usually found in black hole X-ray
binaries ($kT=0.1-0.3$ keV), while the powerlaw is steeper ($\Gamma=2-4.5$)
and cut off at much lower energies ($2-7$ keV compared to $\gtrsim 60$ keV in
black hole X-ray binaries). The absorption is intermediate, with column
densities generally in the range $1-3\times 10^{21}$ cm-2. This led the
community to imagine a different accretion state than in X-ray binaries
(XRBs), namely the ultraluminous state (Gladstone et al., 2009; Kaaret et al.,
2017). ULXs are also likely to cause important feedback on their environment,
with optical and radio signatures of a “bubble” surrounding them, inflated by
powerful winds and intense photo-ionisation (e.g. Abolmasov et al. 2008;
Berghea et al. 2020; Gúrpide et al. 2022), which may seem inconsistent with a
narrow beaming of X-ray emission that could explain a super-Eddington
luminosity (Pakull & Mirioni, 2002; Berghea et al., 2010).

On the other hand, the more luminous and much rarer hyperluminous X-ray
sources (HLXs, $L_{X}>10^{41}$ erg s-1) may still remain excellent candidates
to look for IMBH, as evidenced with ESO 243–49 HLX–1 (Farrell et al., 2009),
one of the best IMBH candidates known, reaching an X-ray luminosity of
$10^{42}$ erg s-1 in the 0.3–10 keV band. It shows spectral evolution similar
to that observed in Galactic black holes (Godet et al., 2009; Servillat et
al., 2011; Godet et al., 2012). The black hole in ESO 243–49 HLX–1 is thought
to be fed by episodes of mass-transfer, induced by repetitive partial
stripping of a white dwarf-like star when passing at periapsis (Godet et al.,
2014). This black hole may be embedded in a stellar cluster or be the central
black hole in a dwarf galaxy stripped by or in interaction with ESO 243–49
(Webb et al., 2010, 2017; Farrell et al., 2012). Very few other HLX candidates
are known so far, most often located at high distance (Swartz et al., 2011;
Zolotukhin et al., 2016), where confusion with a nuclear or background source
is more likely. In a recent study, Barrows et al. (2019) identified a large
sample of 169 HLX candidates thanks to the resolving power and high
sensitivity of two surveys, namely Chandra and the SDSS galaxy catalogue.
However, most of their candidates are at $D>1000$ Mpc where a background
source is difficult to detect in the deepest large optical surveys. Their
median unabsorbed luminosity of $4.6\times 10^{42}$ erg s-1 in the band $2-10$
keV, as well as the significant fraction of HLXs having an optical counterpart
(28%), could also be explained by a significant number of objects being
background sources, although they interpret them as active galactic nuclei
(AGN) in satellites of the host galaxy and report a background contamination
rate of $7-8$% estimated from the cosmic X-ray background curves of Moretti et
al. (2003). In contrast, NGC 5907 ULX-1 occasionally exceeded $10^{41}$ erg
s-1 in luminosity while hosting a neutron star, appearing as an exception in
this class of objects (Israel et al., 2017) . In the following, HLX refers to
objects whose average broad-band luminosity over all X-ray detections by a
given instrument exceeds $10^{41}$ erg s-1, which is not the case for this
one.

Studying ULX samples has led to a statistical picture of the ULX population,
using many catalogues of ULX candidates from old (Swartz et al. 2004; Liu &
Bregman 2005; Liu & Mirabel 2005; Swartz et al. 2011; Walton et al. 2011; Wang
et al. 2016) and recent (Earnshaw et al. 2019; Kovlakas et al. 2020; Inoue et
al. 2021; Bernadich et al. 2022; Walton et al. 2022) releases of large X-ray
catalogues. Focusing either on the hardness (e.g. Earnshaw et al. 2019) and
variability of ULXs (e.g. Sutton et al. 2013; Bernadich et al. 2022) or on
their environment and counterparts (e.g. Kovlakas et al. 2020), or on their
X-ray luminosity function (XLF, Mineo et al. 2012; Wang et al. 2016), these
studies showed several features for the ULX population. ULXs are more
prominent in spiral, star-forming galaxies. Their frequency increases with
galaxy mass and their luminosity function may be explained as a mere powerlaw
extension of the high-mass X-ray binary luminosity function, for which a break
is noted at a few $10^{38}$ erg s-1 in elliptical galaxies (e.g. Kim &
Fabbiano 2004, Wang et al. 2016). A break or cutoff at $1-2\times 10^{40}$ erg
s-1 may also be present (Swartz et al., 2011; Mineo et al., 2012) in the XLF
of spiral galaxies, however this is still under debate (see e.g. Walton et al.
2011; Wang et al. 2016). Moreover, no break is systematically observed in
these galaxies at the Eddington luminosity of a $2\mathrm{M}_{\odot}$ neutron
star ($L_{Edd}\sim 2\times 10^{38}$ erg s-1) or a $100\mathrm{M}_{\odot}$
stellar black hole ($L_{Edd}\sim 10^{41}$ erg s-1), as would be expected if
the Eddington luminosity was a hard limit (Fabbiano, 2006; Kaaret et al.,
2017). ULXs in spiral (resp. elliptical) galaxies are found to be consistent
with a dominant population of high-mass (resp. low-mass) X-ray binaries (HMXB
and LMXB), in agreement with the expected environment for these two classes.
HMXBs and LMXBs consist of a black hole or a neutron star accreting matter
from a close and more massive (respectively less massive) star, either because
of stellar winds (Bondi-Hoyle accretion, Bondi & Hoyle 1944; Bondi 1952) or by
the star overflowing its Roche lobe and forming an accretion disc (e.g. van
den Heuvel & De Loore 1973). In HMXBs, because mass transfer by Roche lobe
overflow is unstable when the binary has a large mass ratio, accretion is
generally driven by stellar winds (e.g. Iben & Livio 1993; Nelemans et al.
2000; however see Pavlovskii et al. 2017 where this point is discussed). LMXBs
reside near the galactic bulge or in globular clusters, while HMXBs, for which
the compact object is most often identified as a neutron star (e.g. Sidoli &
Paizis 2018), are located close to the galactic plane and at larger
separations from the bulge (Grimm et al., 2002; Repetto et al., 2017), where
massive stars are more frequent due to more recent star formation episodes.

Some differences were still reported between XRB and ULX populations. However,
only a handful of the brightest (and often closest) ULXs have been studied in
detail. Many of the conclusions on the population of ULXs are often drawn from
simple model fitting, which does not give a physical description of the data
(e.g. Swartz et al. 2004), or on easily computed quantities such as the
hardness ratio between two bands. A large sample approach was carried out by
Walton et al. (2011) who reported a small offset of ULXs towards harder
accretion states when residing in spiral galaxies. Earnshaw et al. (2019)
found no significant difference between the hardness ratio distributions of
ULXs and XRBs. Likewise, the larger dataset of Bernadich et al. (2022) showed
similar locations for these two populations in the hardness-hardness diagrams.
In contrast, ULXs seem less prone to exhibit high variability on the timescale
of years compared to XRBs, and the most variable ones are hosted in spiral
galaxies (Bernadich et al., 2022). Confirmed neutron star ULXs were reported
to be the hardest ULXs, with highly variable hard radiation (Gúrpide et al.,
2021). Some ULXs are also variable on the timescale of a few kiloseconds
(within a single observation), as found by Earnshaw et al. (2019) and studies
of individual ULXs suggest the presence of winds, outflows or pulsations (e.g
Sutton et al. 2013; Koliopanos et al. 2019).

Nevertheless, ULX population studies suffer from a moderately high
contamination fraction due to foreground stars and (mainly) background AGN,
often estimated to be $\sim$20% from the $\log N-\log S$ diagram of cosmic
X-ray background sources (e.g. Walton et al. 2011), increasing with X-ray
luminosity (e.g. Bernadich et al. 2022). Despite evidence for ULXs being
super-Eddington stellar-mass accretors, the physical mechanism behind the
apparent super-Eddington luminosities is still poorly constrained. Modelling
over the last twenty years suggests different possible scenarios such as
anisotropic (beamed) emission (e.g. King 2009; Wiktorowicz et al. 2019),
supercritical accretion discs around stellar-mass black holes with radiation-
driven winds (e.g. Poutanen et al. 2007; Middleton et al. 2015) or around
neutron stars (e.g. Erkut et al. 2019; Kuranov et al. 2020), thermal-timescale
mass transfer - in particular from a Helium-burning secondary (e.g.
Wiktorowicz et al. 2015; Pavlovskii et al. 2017), accretion flows around
pulsating highly-magnetised neutron stars (e.g. Bachetti et al. 2014;
Mushtukov et al. 2015; Israel et al. 2017; Mushtukov et al. 2017; Koliopanos
et al. 2017) and accretion of clumpy stellar wind enhanced by X-ray ionization
(Krtička et al., 2022).

In this work, we develop a novel approach to build cleaner ULX and HLX
samples, using a general-purpose classification of X-ray sources developed in
Tranin et al. (2022). We assess the added value of such a clean sample, and we
conduct a statistical study of our sample – the largest to date – to offer a
more complete view of ULX and HLX properties.

In Section 2 we explain our ULX selection method, involving the use of a
classification of X-ray sources, and the way we validate the sample and build
a complete subsample. The sample of contaminants and the selection of HLX
candidates are also described. In Section 3 we present a statistical study of
our clean sample, in terms of X-ray luminosity function, ULX rate evolution
with the environment, hardness and variability properties. Samples of XRBs,
ULXs and HLXs are compared. We discuss these results in Section 4 and compare
them to catalogues and results obtained in previous works. The nature of HLXs
is also investigated. We summarise our study in Section 5. Unless stated
otherwise, errors are quoted at the 1-$\sigma$ level.

##  2 The sample

###  2.1 X-ray catalogues

In order to retrieve the largest ULX sample possible, we consider the three
largest X-ray catalogues to date, generated from observations with XMM-Newton,
the Chandra X-ray observatory and the Swift X-ray telescope (XRT).

From the largest of these catalogues, 4XMM-DR11 (Webb et al., 2020), we select
the 496645 point-like sources that have a detection with a reasonable
detection flag $SUM\\_FLAG\leq 1$ (source parameters may be affected, but
spurious detections are unlikely). A few extended sources are considered as
point-like by the XMM-Newton pipeline ($SC\\_EXTENT=0$) and remain in this
selection. They have typically $SC\\_EXT\\_ML>100$ and $SC\\_SUM\\_FLAG>1$ or
$SC\\_EXT\\_ML>10^{4}$, so we remove the 561 corresponding sources. We also
remove sources below the detection likelihood $SC\\_DET\\_ML=10$ threshold, to
limit the false source rate while still keeping a high fraction of the
catalogue sources (83%).

The Chandra CSC2 catalogue (Evans et al., 2010, 2019) is treated in a similar
way: removing sources flagged as extended or ambiguous (as given in the source
name), saturated or overlapping a bright streak. Unlike in the XMM-Newton
pipeline, extended sources seem to be all flagged as extended in the Chandra
pipeline, but some genuinely point-like sources may also be flagged likewise:
from manual inspection, we choose to keep 1217 sources flagged as extended but
not being in a confused area and whose median major axis across energy bands
(computed in the Chandra pipeline) is at most 0.5 arcsec (i.e. the pixel scale
of ACIS cameras). We keep only sources in the likelihood class TRUE and with a
non-zero flux, having detections with a $conf\\_code$ lower than 256 (to
remove sources with unreliable parameters due to an overlapping extended
source); 214757 sources remain after this filtering process. As detailed in
Section 3.1, 199 sources flagged by the Chandra pipeline and not included at
this stage were added a posteriori to the ULX catalogue after manual
inspection.

Likewise, from the Swift-XRT catalogue 2SXPS (Evans et al., 2020) we select
sources having $detflag=0$ and $fieldflag\leq 1$ to remove possibly spurious
sources and sources in polluted fields, as well as sources with a zero flux.
Since Swift-XRT has the largest PSF among the three X-ray telescopes, some
supposedly point-like sources may actually be extended at this stage, and no
extent estimate is made available. This was notably confirmed by a visual
inspection of some fields containing galaxy clusters. We thus flagged 2SXPS
sources matching the extent of an XMM-Newton or Chandra extended source,
unless they also match a reliable source at less than 10 arcsec. In galaxy
clusters, because the cluster is extended, its Simbad and XRT position are
sometimes separated by $\sim$20 arcsec: 188 sources at less than 20 arcsec
from a Simbad galaxy cluster were conservatively flagged as well, as possibly
due to the hot gas extended emission. This leaves a sample of 130162 clean
sources.

###  2.2 X-ray matches

In order to assess the number of unique ULXs left in our final sample, and to
be able to probe long-term variability of these sources, we perform a
crossmatch of the 3 catalogues with each other using the TOPCAT software (Tool
for Operations on Catalogues and Tables, Taylor 2005). The Sky with Errors
algorithm is used with the 3–$\sigma$ position error of each catalogue.
Grouped matches are found, meaning the ambiguous association of a source to at
least two sources from the other catalogue, either in crowded fields (like
galaxies) or because of spatial resolution issues (it is not unusual that 2
close Chandra sources are confused into a single XMM-Newton or Swift one,
since these instruments cannot resolve them). We choose to remove these
associations and flag the corresponding source, unless they become ungrouped
when a crossmatch using the 1–$\sigma$ position error is performed. In this
way, ambiguous associations in which one of the potential associations is
clearly favoured are well retrieved, while more ambiguous associations are
flagged and removed.

From the 412242 XMM-Newton, 214757 Chandra and 130162 Swift unique sources,
the cross-correlation retrieves 687291 unique sources including 65366 with
detections in at least two facilities and 4229 with a grouped match flag. In
particular, 1.7% of 2SXPS sources matching a CSC2 source are in ambiguous
association. The median separation between the two counterparts of an
ambiguous 2SXPS-CSC2 (4XMM-CSC2) association is 8 arcsec. Consequently, for
sources associated with a galaxy, this fraction of ambiguous associations
increases with galaxy distance and becomes $\sim$10% in galaxies at 40 Mpc
(Figure 1). Similarly, 0.7% of 4XMM-CSC2 associations are ambiguous, their
median separation is 5.5 arcsec, and $\sim$5% of 4XMM-CSC2 associations are
ambiguous in galaxies at 40 Mpc. At greater distances, some Chandra sources
must be confused as well. This source confusion issue is important in ULX
studies because tight groups of XRB (resp. ULX) can be mistaken for ULXs
(resp. HLXs) (e.g. Wolter et al. 2015).

###  2.3 The galaxy sample

Figure 1: Fraction of 2SXPS-CSC2 and 4XMM-CSC2 ambiguous matches as a function
of the host distance given by GLADE.

The first published ULX catalogues were constructed by cross-correlating X-ray
catalogues with rather small catalogues of bright galaxies (in particular RC3,
de Vaucouleurs et al. 1991). More recently, in order to achieve more complete
ULX catalogues, larger galaxy catalogues were considered, such as HyperLEDA
(Paturel et al., 2003), the Catalogue of Neighbouring Galaxies (Karachentsev
et al., 2004) or HECATE (Heraklion Extragalactic Catalogue, Kovlakas et al.
2021). Here we use a recent compilation of galaxy catalogues, GLADE (Galaxy
List for the Advanced Detector Era, Dálya et al. 2018), intended to help
locate the origin of gravitational wave events. It is 100% complete in bright
galaxies (defined as those accounting for half of the integrated Schechter
luminosity function, Dálya et al. 2018) up to 91 Mpc and more than 90%
complete at 200 Mpc. Unlike more recent versions, its 2016 version, which is
based on HyperLEDA, GWGC (the Gravitational Wave Galaxy Catalogue, White et
al. 2011) and 2MASX (the 2MASS extended source catalogue, Skrutskie et al.
2006), contains information on the extent of each galaxy. It contains more
than 270000 entries at D¡200 Mpc, $\sim$30% more than HECATE for the same
distance range. We found that some extended infrared sources from 2MASX,
considered as galaxies in GLADE, are actually young stellar objects or diffuse
emission from nebular regions in our Galaxy: from a manual inspection of the
ones matching an X-ray source, we found that their B magnitude, inferred in
GLADE, is generally brighter than 14, or that their Gaia colours are distinct
from actual galaxies (whose $G$ and $BP$ Gaia magnitudes follow $G-BP>0$). We
therefore removed the corresponding 2MASX entries. We also remove entries
matching a stellar object in Simbad within 10 arcsec (typical angular size of
the closest young stellar objects), unless they are also categorised as
galaxy. Since the major axis of 2MASX entries is missing in the catalogue, we
retrieve the $r\\_fe$ column native of 2MASX. Last but not least, because
GLADE is a compilation of catalogues, some duplicate entries are found. We
retrieve 6725 such galaxies from a Sky internal crossmatch of 10 arcsec.
Although their distances in each catalogue are generally close to each other,
large differences can occur when the distance is photometrically estimated. In
these cases, we favour the distance of HyperLEDA or GWGC over the one of
2MASX, as they are more consistent with each other (Dálya et al., 2018). To
remove Galactic globular clusters present in GLADE, where we do not expect the
presence of ULXs, we limit the sample to distances above 1 Mpc. The resulting
galaxy sample is composed of $\sim$1.7 million galaxies essentially at
$B_{mag}$ brighter than 19 (median $B_{mag}=17.3$).

Unlike HECATE, GLADE does not contain any information on the galaxy
morphology, star-formation rate (SFR) or stellar mass. The latter is however
estimated from infrared integrated luminosity in the latest release of GLADE,
GLADE+ (Dálya et al., 2022), so we retrieve this column from this catalogue.
The Hubble type $t$ of galaxies is given in HyperLEDA, so this column is
retrieved as well and provides morphology information for $\sim$460000
galaxies. Such information is necessary to distinguish the study of ULX
populations in spiral and elliptical galaxies, which have significant
differences as mentioned in Section 1. Similar to the definition used by
Earnshaw et al. (2019) and Bernadich et al. (2022) , we define spiral and
elliptical galaxies as those with $t\geq 0$ and $t<0$, respectively. Irregular
galaxies are thus included in the spiral sample in the same way as Walton et
al. (2011) and Walton et al. (2022) . To get a more complete census of spiral
and elliptical galaxies, we cross-correlate GLADE with the catalogue of galaxy
morphologies inferred by machine learning on PanSTARRS images (Goddard &
Shamir, 2020). As a result, we obtain stellar masses and morphologies for
$\sim$1.5 million and $\sim$1 million galaxies, respectively. Last but not
least, as done in HECATE and introduced by Cluver et al. (2017) (see also
Kennicutt & Evans 2012), we cross-correlate GLADE with the WISE (Wide-field
Infrared Survey Explorer) catalogue of infrared sources (Cutri et al., 2012)
to estimate the star-formation rate from the W3-band (12 $\mu$m) magnitude and
the galaxy distance. More than 99% of our GLADE subset has a WISE counterpart.
This method provides consistent values of SFR for late-type galaxies, while
the heated dust content generally dominates in early-type galaxies, often
leading to an SFR overestimation (e.g. Galliano et al. 2018; Kovlakas et al.
2021).

###  2.4 ULX candidates

To form the ULX sample, the first step is to identify X-ray sources matching a
galaxy, but outside its nuclear region to avoid AGN contaminants. We used the
matching tool of TOPCAT to this end, with the Sky Ellipses algorithm to
retrieve all sources whose X-ray error circle overlap the ellipse representing
the galaxy area. The major axis of the galaxy is chosen to be the Holmberg
diameter $D_{Holm}=1.26\times D_{25}$ as in Walton et al. 2022, where $D_{25}$
is the isophotal diameter at surface brightness 25 mag arcsec-2, in order to
retrieve more candidates and study the spatial distribution of ULXs out to
large radii. The resulting samples of X-ray matches are used later as input
for source classification. They contain 18506, 13055 and 7243 sources for
Chandra, XMM-Newton and Swift, respectively.

To exclude the central region, we select sources satisfying the criterion

${\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
d>3(POSERR+0.5)}$ (1)

where $d$ is the angular separation to the galaxy centre and $POSERR$ is the
X-ray position error at $1\sigma$. A minimal offset of 3 arcsec is required.
The 0.5 arcsec term is applied to correct for astrometric errors on the galaxy
centre reported in GLADE – it is the typical offset between GLADE and SDSS
positions. The mean observed X-ray luminosity over the instrument broad band
(0.2-12 keV for XMM-Newton, 0.5-7 keV for Chandra and 0.3-10 keV for Swift) is
computed from the mean X-ray flux given in the X-ray catalogues and the galaxy
distance, and must strictly exceed $L_{X}>10^{39}$ erg s-1 for the source to
be considered as a ULX candidate. Sources brighter than $10^{41}$ erg s-1 are
kept in the same sample at this stage.

Figure 2: Illustration of ULX selection on the galaxy NGC 1365. (Left) XMM-
Newton X-ray image (Right) DSS Optical image. The $D_{25}$ ellipse of the
galaxy is shown in grey. Blue, orange, green and red circles represent sources
classified as AGN, soft sources, XRB and initial ULX candidates, respectively.
The dashed red circle corresponds to a candidate classified as a contaminant,
and is indeed a background AGN. An interactive view of this figure is
available at https://xmm-
ssc.irap.omp.eu/claxson/xray_analyzer.php?srcquery=53.4019%20-36.1406

###  2.5 Source classification to filter contaminants

One main issue in ULX studies is the high rate of contaminants in ULX samples,
due to the presence of foreground (essentially stars in our Galaxy) and
background (AGN) sources. Some central AGN may also pass the criterion cited
above because of bad astrometry. Previous attempts to remove these
contaminants mainly focused on removing known AGN and stars (e.g. Kovlakas et
al. 2020; Bernadich et al. 2022) using the large existing datasets for these
types. However, this is not sufficient to remove all contaminants, since many
of them remain uncatalogued. Another option, applied in Bernadich et al.
(2022), is to remove all sources that are too bright in optical to be a ULX.
To this end, they use Gaia and PanSTARRS magnitudes obtained by a positional
crossmatch with X-ray sources to compute the X-ray to optical flux ratio: any
source with $\log(F_{X}/F_{Opt})<-2.2$ is considered as a star, and other
sources with $\log(F_{X}/F_{Opt})<0$ were manually inspected to flag likely
background sources while keeping the bright HII star-forming regions where
ULXs are preferentially found.

While efficient, this process is tedious and also misses all possible
background AGN having $F_{X}>F_{Opt}$. To overcome this issue, we make use of
the automated probabilistic classification of X-ray sources we recently
developed (Tranin et al., 2022). In a nutshell, the dataset is first prepared
by enriching the X-ray catalogue by completing X-ray sources with the
following information: their optical and infrared counterparts from large
ground-based surveys, identified with the Bayesian cross-matching tool Nway
(Salvato et al., 2018); galaxies hosting the source, as explained in Section
3.3; X-ray variability ratio between multi-instrument observations; and source
identification given by external catalogues, to form the training sample. The
classification scheme starts with a Naive Bayes Classifier based on the
densities of the training sample for all source properties obtained after this
enrichment. The scheme is then fine-tuned to maximise the classification
performance of a chosen class (here X-ray binaries), by increasing the weight
of the most discriminating properties. The probabilities of each class as well
as the final class (giving maximum probability) are computed. Here, the X-ray
sample of sources matching a GLADE galaxy is divided into a training sample of
known sources and a test sample to classify. The classification is based on
three classes: AGN (corresponding to background contaminants and most sources
located in galaxy centres), X-ray binaries (the sample of in-situ genuine
sources) and soft sources (a miscellaneous class containing stars, supernovae
and supernova remnants, which are in-situ and foreground contaminants). The
properties used in the classification process and the list of catalogues used
to retrieve AGN, stars and XRB are detailed in Tables 1 and 2 of Tranin et al.
(2022). Supernovae and supernova remnants are retrieved through a match with
Simbad. The parameters of the classifier for each catalogue are summarised in
Table 1, while the sample size of the training samples of XMM-Newton, Chandra
and Swift are detailed in Table 2, as well as the classification results in
the training samples. The classifier is very efficient at retrieving AGN, and
retrieves more than 80% of X-ray binaries. False positive rates are also low
($\lesssim$15%) and can be further decreased by applying a probability
threshold on the class under consideration.

Figure 3: Densities of some source properties for each class of the 4XMM-DR11
training sample. $SC\\_HR2$: second hardness ratio of 4XMM-DR11 (between bands
0.5-1 and 1-2 keV), from EPIC-pn and MOS cameras. $\log(F_{X}/F_{W1})$:
logarithm of the X-ray to infrared (W1-band, 3.4 $\mu$m) flux ratio. $L_{X}$:
mean 0.2–12 keV observed luminosity. Separation/Gal. radius: galactocentric
distance, i.e. the source separation to the host centre, in units of the host
radius at the source position angle.

As a sanity check, we obtain that 99% of sources at $L_{X}>10^{42}$ erg s-1
and 94% of sources matching the galaxy centre (see Equation (1)) are
classified as AGN. The vast majority of other central sources are classified
as soft sources, mostly being absorbed (Seyfert II) AGN or hot gas.

Table 1: Classification parameters for 4XMM, CSC2 and 2SXPS.

4XMM
---
Priors $\mathcal{P}_{AGN},\leavevmode\nobreak\ \mathcal{P}_{Soft},\leavevmode\nobreak\ \mathcal{P}_{XRB}$ | 67%, 8%, 25%
Weights $\alpha_{loc},\leavevmode\nobreak\ \alpha_{spe},\leavevmode\nobreak\ \alpha_{ctp},\leavevmode\nobreak\ \alpha_{var}$ | 4.8, 9.1, 1.0, 5.5
CSC2
Priors $\mathcal{P}_{AGN},\leavevmode\nobreak\ \mathcal{P}_{soft},\leavevmode\nobreak\ \mathcal{P}_{XRB}$ | 50%, 10%, 40%
Weights $\alpha_{loc},\leavevmode\nobreak\ \alpha_{spe},\leavevmode\nobreak\ \alpha_{ctp},\leavevmode\nobreak\ \alpha_{var}$ | 6.9, 5.3, 3.6, 4.2
2SXPS
Priors $\mathcal{P}_{AGN},\leavevmode\nobreak\ \mathcal{P}_{soft},\leavevmode\nobreak\ \mathcal{P}_{XRB}$ | 67%, 8%, 25%
Weights $\alpha_{loc},\leavevmode\nobreak\ \alpha_{spe},\leavevmode\nobreak\ \alpha_{ctp},\leavevmode\nobreak\ \alpha_{var}$ | 9.6, 9.1, 1.3, 4.5
111For each catalogue are shown (first line) the prior proportions applied to
each class and (second line) the weights for properties related to location,
spectrum, multiwavelength counterparts and variability, once optimised to
maximise the XRB classification performance. See Table 1 and Equation (3) of
Tranin et al. (2022) for further details.

Table 2: Confusion matrixes resulting from the classification of 4XMM-DR11,
CSC2 and 2SXPS training samples. 4XMM
---
| AGN | Soft | XRB | Precision
$\rightarrow$AGN | 1482 | 26 | 10 | 98%
$\rightarrow$Soft | 15 | 33 | 18 | 50%
$\rightarrow$XRB | 29 | 20 | 320 | 87%
Recall | 97% | 42% | 92% |
CSC2
| AGN | Soft | XRB | Precision
$\rightarrow$AGN | 996 | 27 | 58 | 92%
$\rightarrow$Soft | 25 | 140 | 27 | 73%
$\rightarrow$XRB | 81 | 80 | 830 | 84%
Recall | 90% | 57% | 91% |
2SXPS
| AGN | Soft | XRB | Precision
$\rightarrow$AGN | 1568 | 31 | 14 | 97%
$\rightarrow$Soft | 8 | 43 | 13 | 67%
$\rightarrow$XRB | 65 | 36 | 363 | 78%
Recall | 96% | 39% | 93% |

To remove contaminants, the probability that the source is an XRB is simply
the probability that it is not a contaminant (indeed
$P_{XRB}=1-P_{AGN}-P_{Soft}$). It is then compared to the fraction of
contaminants in the galaxy area within the source separation. This fraction is
obtained with the formula:

$\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
f_{cont}=2\frac{\pi abn_{cont}sep^{2}}{N_{ULX}}$ (2)

where $a$ and $b$ are the galaxy semi-major and semi-minor axes, $n_{cont}$ is
the density of contaminants computed with the analytical formula of the
log(N)-log(S) relation in Moretti et al. (2003), for the hard cosmic X-ray
background, $sep$ is the galactocentric distance (i.e. the source separation
in units of the galaxy radius at its position angle), and $N_{ULX}$ is the
number of ULX candidates in the considered galaxy. Regardless of this
quantity, sources with a probability to be a contaminant higher than $\sim$95%
are also classified as contaminant: it is notably the case of many foreground
sources, and many of the few spurious sources remaining in the sample.
Selected candidates are thus those following $P_{XRB}>max(0.05,f_{cont})$.
Here, $P_{XRB}$ should be understood as the probability that the ULX candidate
is not a background or foreground contaminant, rather than the probability
that the ULX is actually an XRB, which requires further study to be affirmed.
Four known ULXs are contained in this training sample of XRBs. From manual
inspection of 150 sources, we estimate that more than 90% of reliable ULXs are
retrieved (as detailed in Section 2.9) and that at most 15% of selected
candidates are compatible with background contaminants. Removed candidates
represent $\sim$42% of the initial candidates. This high fraction is the
result of three factors: first, the selection out to $1.26D_{25}$ increases
the galaxy area by 60% and the number of background contaminants by the same
amount. Second, foreground contaminants and spurious sources (resulting from
the X-ray detection pipelines) are also removed in this process. Third, a few
valid ULX candidates are removed as well, because their properties used in
this work are also compatible with an AGN (e.g. NGC 3921 ULX X-2, Jonker et
al. 2012).

###  2.6 Filtering remaining contaminants

We manually verified a large number of our selected ULX candidates, which
greatly helped to develop and assess the filtering pipeline described above.
To this end, we used the virtual observatory tools Aladin
Lite222https://aladin.u-strasbg.fr/AladinLite/ (Bonnarel et al., 2000; Boch &
Fernique, 2014), the Simbad database333https://simbad.cds.unistra.fr/simbad/
and the VizieR catalogue access tool444https://vizier.cds.unistra.fr/viz-
bin/VizieR (Ochsenbein et al., 2000), developed at CDS, Strasbourg
Observatory, France.

Some contaminants remain after the filtering process described above: to
remove them, we match our sample with Simbad (3 arcsec) and exclude objects of
types AGN or stellar objects (75% of these objects were already successfully
identified by the classification). We also inspect manually all sources with
an optical counterpart and selected as ULX, having optical colour $b-r>0.5$,
due to the redder nature of background AGN (about 700 sources). There are 135
sources from XMM-Newton, 139 from Chandra and 53 from Swift that are discarded
in this process, bringing the expected contamination rate to about 2% (Section
4.2). In the following, we refer to the resulting cleaned sample as ”selected
ULXs”.

In the same way as the contamination rate, the false negative rate is expected
to increase with luminosity, because the likelihood of the AGN class is
enhanced (Figure 3, lower left panel).

###  2.7 HLX sample

Hyperluminous candidates were selected as described in Section 2.4, as off-
nuclear sources with a mean observed X-ray luminosity $>10^{41}$ erg s-1 in
the broad energy band. To limit the fraction of spurious sources, only
candidates having a $S/N>3$ were kept. Respectively 195, 360 and 110 sources
from CSC2, 4XMM and 2SXPS satisfy these criteria. Each of these HLX candidates
was manually inspected to remove spurious sources (notably due to source
confusion issues) and contaminants such as stars, hot gas overdensities and
jet hotspots. This step led to the removal of 157 sources. Likewise, 33
sources were removed as being mistakenly associated to GLADE galaxies, whose
extent was manifestly overestimated.

Unlike ULXs, which have well-constrained properties that are essentially
similar to XRBs for the classifier, selecting HLXs from the classification
results can induce important biases. In particular, bright, reliable HLX
candidates having a Gaia (Gaia Collaboration et al., 2021), PanSTARRS
(Chambers et al., 2016) or DES (Dark Energy Survey, Abbott et al. 2018)
optical counterpart are misclassified as background AGN. A more robust
approach to eliminate background AGN is to use the redshift (spectroscopic or
photometric) of the optical counterpart (as done by Barrows et al. 2019).
Recently, in addition to spectroscopic measurements of the redshift for a few
million sources (e.g. in the SDSS-BOSS survey, Bolton et al. 2012), large and
deep surveys led to the release of several $10^{9}$ photometric redshifts. In
SDSS, this work was notably performed by Beck et al. (2016), inferring the
redshift of 208 million sources up to $z\sim 0.6$ from their (u,g,r,i,z)
magnitudes. Tarrío & Zarattini (2020) computed the photometric redshift of 1.1
billion sources using the (g,r,i,z,y) PanSTARRS bands, extending up to
redshift $z\sim 1$. Zou et al. (2022) performed a similar work using the
(g,r,i,z,Y) and (g,r,z,W1,W2) bands from the Dark Energy Survey and the DESI
Legacy Survey, respectively. These three studies provide highly reliable
photometric redshifts, with a typical accuracy better than $\sigma(\Delta
z_{norm})=0.03$. Besides, most recently, identification as quasars or galaxies
and redshifts were provided for 7.8 million Gaia DR3 sources (Gaia
Collaboration et al., 2022), using the low-resolution optical Gaia spectra.
This set of redshift catalogues was used as input to identify background and
foreground contaminants: X-ray sources at $<3$ arcsec of an optical source
with a redshift measurement inconsistent with the distance of the assumed host
(i.e. $|z_{host}-z|>z_{err}$, with $z$ and $z_{err}$ the redshift value and
error in the external catalogue) were discarded. This led to the removal of
283 background and 5 foreground sources, respectively. In contrast, sources
having a distance match ($|z_{host}-z|<z_{err}$) have more chance to be indeed
located in their associated galaxy: this is the case of 37 sources. At the end
of this selection process, 115, 75 and 13 HLX candidates are found in 4XMM-
DR11, CSC2 and 2SXPS respectively, totalling 191 unique sources including 11
observed by multiple X-ray instruments at the HLX level according to their
mean luminosity (for X-ray sources observed with several instruments, the
entry of smallest position error is kept). From them, 63% have no counterpart
in Gaia, PanSTARRS and DES catalogues, 24% have a counterpart but no redshift
measurement and 13% have a redshift consistent with their host.

###  2.8 Complete subsample

As introduced for example in Walton et al. (2011), ULX properties and
luminosity function have to be assessed on a complete sample, so that the
contribution of brighter ULXs is not overestimated. Indeed, the limiting
sensitivity of the instrument defines the minimum flux that a source should
have to be detected; for extragalactic sources, this translates into a
limiting luminosity for each distance bin. Usually, ULX surveys keep only
galaxies with limiting luminosity below $10^{39}$ erg s-1 to ensure that all
ULXs are well-detected. For instance, Kovlakas et al. (2020) (resp. Bernadich
et al. 2022) keep galaxies closer than 40 Mpc (resp. 29 Mpc). This leads to
the removal of most bright ULXs located in more distant galaxies. In this
study, we keep all sources with luminosities above the galaxy limiting
luminosity $L_{lim}=4\pi F_{lim}D^{2}$. As further detailed in Section 3.1,
the contribution of ULXs in each luminosity bin is later weighted by the
number of galaxies having $L_{lim}$ below this luminosity.

The only prerequisite is thus to compute the sensitivity $F_{lim}$ for each
galaxy in the sky coverage of each instrument. We use here a data-driven
approach, for simplicity and because sensitivity maps are not accessible for
all 3 instruments. For each X-ray catalogue, we infer the flux – effective
exposure time relation from subsets of $\sim 10^{5}$ detections with signal-
to-noise close to 3. Chandra sources show the highest deviation from a single
powerlaw model, because a third parameter, the off-axis angle $\theta$, is
determinant in the sensitivity value. We empirically find that if we consider
a Chandra effective exposure time $t^{\prime}=2t/max(\theta,2\text{ arcmin})$,
the three X-ray subsets show a $\pm$0.5 dex deviation from a single powerlaw
model. The resulting $F_{lim}$ relations are as follow:

$\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
F_{lim}=\begin{cases}\big{(}\frac{1.2\times
10^{-15}}{t^{\prime}}\big{)}^{0.75}&\text{for CSC2}\\\ \big{(}\frac{1\times
10^{-13}}{t}\big{)}^{0.8}&\text{for 4XMM}\\\ \big{(}\frac{5\times
10^{-13}}{t}\big{)}^{0.8}&\text{for 2SXPS}\end{cases}$ (3)

The limiting X-ray sensitivities of Chandra, XMM-Newton and Swift are
illustrated in Figure 4, showing the cumulative distribution of the limiting
X-ray flux for GLADE galaxies in their sky coverage. The difference between
on-axis and off-axis Chandra galaxies is to be noted. Once applied to X-ray
sources matching a GLADE galaxy, the $F>F_{lim}$ cut keeps 70 to 95% of
sources with $S/N>3$ and removes 80 to 90% of sources with $S/N<3$. This
empirical sensitivity can thus be considered as a proxy for the $3\sigma$
sensitivity.

Figure 4: Cumulative distribution of GLADE galaxies X-ray 3$\sigma$
sensitivities, computed from Equation (3) as a function of flux. $\theta$
refers to the Chandra off-axis angle.

###  2.9 Contaminants sample

In order to validate the automatic filtering method of the test sample, we
compiled a sample of objects classified as contaminants from the three
catalogues to compare it with the compilation of selected ULX. They are ULX
candidates following $d>3(POSERR+0.5)$ and $P_{XRB}<max(0.05,f_{cont})$, and
having no identification as one of the three classes AGN, soft source or XRB
in Simbad or other catalogues, in order not to bias the classification
assessment. This results in a sample of 1431 sources (1331 unique) with 494,
715 and 222 sources from CSC2, 4XMM and 2SXPS, respectively.

We manually inspected 100 sources having S/N¿3 from this sample: 83% were
indeed contaminants (mostly background AGN), 7% may be reliable ULXs missed by
our classification, and 10% were ambiguous cases, either lying far from the
galaxy optical extent but having no optical counterpart, or having source
confusion issues leading to unreliable properties.

The distributions of galactocentric distance, mean X-ray luminosities, GLADE
distance and Hubble type are shown in Figure 5 for both selected ULXs and
contaminant samples. These results agree with what is reported in previous
works: sources further from the galaxy centre, of higher luminosities or in
elliptical galaxies are more prone to be contaminants. In particular, before
manual inspection of the HLXs, all but 6 sources above $10^{42}$ erg s-1 were
classified as contaminants. One of the 6 sources is 2CXO J115324.3+493104, an
AGN jet hotspot (Sambruna et al., 2006), proving the presence of unexpected
types of contaminants in the ULX sample. We thus removed it from the HLX
sample, as well as 2CXO J003703.9-010904, another jet hotspot in the HLX
luminosity range (Martel et al., 1998; Kataoka et al., 2003).

Figure 6 shows the fraction of sources classified as contaminants among
initial ULX candidates as a function of the mean X-ray luminosity: in
particular, for a search within the Holmberg diameter 1.26$\times D_{25}$,
more than 30% of candidates above $10^{40}$ erg s-1 are contaminants, even
after removing known contaminants. This evolution is quantitatively consistent
with the 70% contamination rate observed for HLX candidates (Zolotukhin et
al., 2016; Kaaret et al., 2017), although this is dependent on the apparent
size of selected galaxies. Consequently, previous analyses of the bright end
of the X-ray luminosity function of ULXs (e.g. Swartz et al. 2011; Mineo et
al. 2012; Wang et al. 2016) may be severely affected by contamination.

Figure 5: Normalised distributions of some properties of selected ULX
candidates and candidates classified as contaminants, from the compilation of
CSC2, 4XMM-DR11 and 2SXPS. Figure 6: Fraction of contaminants per luminosity
bin, for candidates in the $D_{25}$ and 1.26$\times D_{25}$ ellipses, among
the initial ULX candidates (before filtering).

##  3 Results

The final ULX samples are made up of 1234, 667 and 304 sources from Chandra,
XMM-Newton, and Swift, totalling a compiled sample of 1901 unique sources
(i.e. after removing duplicate entries for sources observed by several
instruments; this sample is hereafter called the compiled sample). The rates
of HLX candidates in each sample are 6%, 16% and 1.9%, respectively. This
difference is expected because of the different capabilities of these X-ray
facilities, as detailed in Section 3.1. Table 5 summarises the number counts
of ULXs, complete ULXs (i.e. with $L_{X}>L_{lim}$), HLXs, and galaxies of
different morphologies in each ULX sample. To our knowledge, our work provides
both the cleanest (about 2% of contaminants, Section 4.2) and the largest
census of ULXs, just above the recently published multi-mission catalogue of
Walton et al. (2022) (1843 sources). Only a handful of well-known ULXs are
absent from our sample. For instance, out of the 17 well-studied ULXs compared
in Gúrpide et al. (2021), only three ULXs are missed: Holmberg IX X-1 and
Circinus X-5, because the $1.26\times D_{25}$ ellipse given in GLADE are
somewhat smaller than the actual galaxy area ; and NGC 55 X-1, because its
mean X-ray luminosity is just below $10^{39}$ erg s-1 for each instrument. All
the eight pulsating ULXs known so far, including the recent pulsating ULX
(PULX) candidate of Quintin et al. (2021), are retrieved in our sample. ESO
243-49 HLX-1 (along with other HLXs, Section 3.4) is retrieved as well.
Conversely, supersoft ULXs studied by Urquhart & Soria (2016) are not in our
sample because of having $\langle L_{X}\rangle<10^{39}$ erg s-1.

Table 3: Samples of ULXs, HLXs and galaxies under study

| CSC2 | 4XMM | 2SXPS | Total (unique)
---|---|---|---|---
ULX | 1234 | 667 | 304 | 2205 (1901)
cULX | 925 | 470 | 266 | 1661 (1448)
HLX | 75 | 115 | 13 | 203 (191)
cHLX | 63 | 92 | 13 | 168 (158)
Galaxies | 741 | 599 | 207 | 1547 (1303)
Spiral | 377 | 320 | 171 | 868 (713)
Elliptical | 206 | 112 | 19 | 337 (325)
555Number counts of ULXs, complete ULXs (i.e. with $L_{X}>L_{lim}$, Section
2.8), HLXs, complete HLXs and number of galaxies that host them (all types,
spiral and elliptical), for CSC2, 4XMM-DR11, 2SXPS and their sum (in
parenthesis: minus duplicates). In the compilation of unique sources, for
X-ray sources observed with several instruments, the entry of smallest
position error is kept. By construction of the HLX sample, the sources all
have $S/N>3$ and thus a higher fraction of them are complete.

###  3.1 Malmquist-corrected XLF

Previous works on the X-ray luminosity function of ULXs only kept galaxies
with limiting X-ray luminosities $L_{lim}<10^{39}$ erg s-1 . While this method
gives the intrinsic shape of the XLF, by probing the same volume in all bins,
it also cuts out most of the detected bright sources, located in further
galaxies. Here, we keep all sources brighter than the limiting luminosity of
their galaxy $L_{X}>L_{lim,host}$. Thus, the distribution of $L_{X}$ is the
XLF convolved with the survey volume, each luminosity bin
$[L_{X},L_{X}+dL_{X}]$ comes from a different survey volume (Malmquist bias,
also known as the Eddington bias), and a 1/$V_{max}$ correction has to be
applied to each bin. This volume $V_{max}$ is here given by the number of
galaxies complete to the luminosity of the bin, that is to say having
$L_{lim}<L_{X,bin}$. The deconvolved differential XLF is computed by dividing
the luminosity distribution of ULXs by the cumulative limiting luminosity
distribution of galaxies, $N(L_{lim}<L_{X,bin})$. Although this approach has
been extensively used to constrain the luminosity functions of stars (e.g.
Leggett & Hawkins 1988; Stobie et al. 1989; Kroupa 1995), galaxies (e.g.
Binggeli et al. 1988; Loveday et al. 1992; Bouwens et al. 2011) and quasars
(e.g. Lynden-Bell 1971; Ueda et al. 2003; Aird et al. 2010), it seems absent
from the literature of XRB and ULX statistical studies at the time of writing.
The underlying assumption is that a universal XLF exists and is the same in
the different survey volumes, which seems reasonable since we study the very
local Universe $z<0.16$ (however see Section 4.4 for further discussion).

Figure 7 shows the deconvolved cumulative XLF (complete in each bin) for
spiral (top panel) and elliptical (bottom panel) galaxies. It includes
selected ULXs below $10^{41}$ erg s-1 from the complete sample, and complete
robust HLXs for the bright end. Poisson errors are assumed in each bin on the
number of X-ray sources and the number of complete galaxies. The three
datasets are in good agreement (consistent at the 90% level except at faint
luminosities where the different resolution capabilities lead to substantially
different numbers of detections, Section 4.3) after making a correction on the
CSC2 sample. Indeed, we noted a deficit of Chandra sources at medium
luminosities ($\gtrsim 10^{40}$ erg s-1 ), which after a thorough inspection
of 4XMM-CSC2 intersecting fields was found to be caused by genuine ULXs that
were flagged in CSC2 as extended or confused, or not detected by the Chandra
pipeline which is less sensitive close to the edges of the field of view.
After manual inspection of 480 flagged sources, we retrieved 199 such sources
including 51 at $L_{X}>10^{40}$ erg s-1 .

We verified that both deconvolved and ULX-complete XLF had the same shape.
Notably, we considered sources at low luminosities $10^{39}-10^{40}$ erg s-1
in spiral galaxies, so that the sample becomes comparable to that of Wang et
al. (2016) (who have only 20 candidates above this luminosity). Using a single
powerlaw model, as they did in this luminosity range, we retrieve an XLF slope
consistent with their result, $\alpha=0.93\pm 0.03$ instead of $\alpha=0.96\pm
0.05$ (1$\sigma$ errors). In elliptical galaxies, we generally find flatter
slopes than in the literature; this point is discussed in Section 4.4.2.

Each XLF is then fitted with two models, a single powerlaw

$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
n(>L)=n_{39}(L/10^{39})^{\alpha}$

and a broken powerlaw

$\displaystyle n(>L)=n_{39}(L/10^{39})^{\alpha_{1}}\leavevmode\nobreak\
\mathrm{if\leavevmode\nobreak\ }L<L_{b}$ $\displaystyle
n(>L)=n_{39}(L/10^{39})^{\alpha_{2}}\leavevmode\nobreak\
\mathrm{\leavevmode\nobreak\ otherwise}$

where $n_{39}$ is the total ULX rate and $L_{b}$ is the break luminosity.
Parameters of the XLF fits and their uncertainties are probabilistically
estimated in a bayesian framework, by sampling the XLF with 40000 Monte-Carlo
trials using observed values and errors. Flat priors are applied on each
parameter, in the ranges $n_{39}\in[0,1]$, $\alpha\in[0,2]$ (single powerlaw)
and $n_{39}\in[0,1]$, $\alpha_{1}\in[0.3,2]$, $\alpha_{2}\in[0,2]$,
$\log(L_{b})\in[39,41]$ (broken powerlaw). Results from the single powerlaw
and broken powerlaw fits of the deconvolved cumulative XLF are detailed in
Tables 6 and 7. Unlike the XLF of elliptical galaxies, in spiral galaxies, the
broken powerlaw fit is always preferred over a single powerlaw fit. However,
the different catalogues disagree on the exact location of the break, even at
the three sigma level. This discrepancy is further discussed in Section 4.4.
However, it is to be noted that all three datasets are well-fitted by a
powerlaw breaking at $L_{X,break}=5\times 10^{39}$ erg s-1 (fixed parameter),
with $\chi^{2}_{r}<1$.

Unlike the lower luminosity break seen in the XLF of elliptical galaxies at a
few $10^{38}$ erg s-1 which is consistent with the Eddington limit of neutron
star binaries (Kim & Fabbiano, 2010), the physical origin of this higher
luminosity break has been poorly understood. Previous works consider it as the
suggestion for a different class of objects above the break luminosity, in
particular intermediate mass black holes, as it corresponds to the Eddington
luminosity of a $\sim 80\mathrm{M}_{\odot}$ black hole (Swartz et al., 2011;
Wang et al., 2016; Kaaret et al., 2017). The reasons for this break are
further examined in Section 4.4.1.

Table 4: Parameters of the fits of the deconvolved XLF with a single powerlaw
model.

Sample | $n_{39}$ | $\alpha$ |  $\chi^{2}_{r}$ (dof)
---|---|---|---
CSC2, ETG | $0.31^{+0.01}_{-0.02}$ | $0.99^{+0.02}_{-0.02}$ | 0.54 (28)
4XMM, ETG | $0.16^{+0.02}_{-0.03}$ | $0.99^{+0.05}_{-0.04}$ | 0.09 (27)
2SXPS, ETG | $0.03^{+0.02}_{-0.01}$ | $1.02^{+0.30}_{-0.24}$ | 0.68 (9)
compiled, ETG | $0.29^{+0.01}_{-0.01}$ | $1.06^{+0.02}_{-0.02}$ | 0.67 (28)

666 “ETG” refer to the sample of complete ULXs in elliptical galaxies.
$n_{39}$ and $\alpha$ are the amplitude and slope of the powerlaw model,
respectively. See the text for more details.

Table 5: Parameters of the fits of the deconvolved XLF of spiral galaxies with
a broken powerlaw model.

Sample | $n_{39}$ | $\alpha_{1}$ | $\alpha_{2}$ | $\log(L_{break})$ |  $\chi^{2}_{r}$ (dof)
---|---|---|---|---|---
CSC2, LTG | $0.45^{+0.02}_{-0.02}$ | $0.67^{+0.09}_{-0.10}$ | $1.36^{+0.06}_{-0.05}$ | $39.45^{+0.07}_{-0.06}$ | 0.37 (27)
4XMM, LTG | $0.35^{+0.02}_{-0.02}$ | $0.81^{+0.06}_{-0.07}$ | $1.24^{+0.05}_{-0.04}$ | $39.73^{+0.11}_{-0.11}$ | 0.19 (33)
2SXPS, LTG | $0.29^{+0.02}_{-0.02}$ | $0.80^{+0.04}_{-0.05}$ | $1.73^{+0.12}_{-0.10}$ | $39.97^{+0.07}_{-0.08}$ | 0.33 (24)
compiled, LTG | $0.41^{+0.01}_{-0.01}$ | $0.86^{+0.04}_{-0.06}$ | $1.31^{+0.03}_{-0.03}$ | $39.66^{+0.08}_{-0.09}$ | 0.77 (33)
CSC2, LTG39 | $240^{+12}_{-12}$ | $0.80^{+0.08}_{-0.10}$ | $1.90^{+0.07}_{-0.14}$ | $39.67^{+0.08}_{-0.11}$ | 1.13 (12)
4XMM, LTG39 | $168^{+10}_{-10}$ | $0.71^{+0.08}_{-0.09}$ | $1.80^{+0.14}_{-0.17}$ | $39.74^{+0.10}_{-0.10}$ | 0.17 (14)
2SXPS, LTG39 | $150^{+10}_{-10}$ | $0.65^{+0.10}_{-0.12}$ | $1.39^{+0.25}_{-0.17}$ | $39.70^{+0.20}_{-0.19}$ | 0.13 (18)
compiled, LTG39 | $424^{+15}_{-15}$ | $0.89^{+0.05}_{-0.06}$ | $1.97^{+0.03}_{-0.06}$ | $39.71^{+0.05}_{-0.06}$ | 0.60 (12)
777“LTG” refer to the samples of complete ULXs in spiral galaxies. “LTG39”
stands for the samples of spiral galaxies having $L_{lim}<10^{39}$ erg s-1;
they correspond to the classical, non-renormalised XLF. $n_{39}$,
$\alpha_{1}$, $\alpha_{2}$ and $\log(L_{break})$ are the amplitude, faint-end
slope, bright end slope and break luminosity of the broken powerlaw model,
respectively. See the text for more details.

Figure 7: Cumulative X-ray luminosity functions of ULX candidates from 4XMM,
CSC2 and 2SXPS. (Top panel) Deconvolved XLF of spiral galaxies, and residuals
from the single (middle) and broken (bottom) powerlaw fits detailed in Table
7. (Bottom panel) Deconvolved XLF of elliptical galaxies, and residuals from
the single powerlaw fits detailed in Table 6. For comparison, we overlay the
cumulative distributions of sources matching their host nucleus, and of ULX
candidates classified as background contaminants.

###  3.2 ULX rates in different environments

Generally speaking, precursory studies found that spiral galaxies were more
prone to host ULXs, in particular in their star-forming regions (e.g. King et
al. 2001; Swartz et al. 2004; Liu & Bregman 2005). Elliptical galaxies also
host a significant ULX population, especially those that have undergone a
recent star-formation event (Kim & Fabbiano, 2004). The luminosity function of
ULXs and their location in their host (distance to the centre) seem to be
consistent with the extrapolation of HMXB in spiral galaxies and LMXB in
elliptical galaxies, at higher luminosities (Gilfanov et al., 2004; Swartz et
al., 2011; Mineo et al., 2012; Kovlakas et al., 2020). Other demographic
studies showed an excess of ULXs in dwarf galaxies and in low-metallicity
galaxies (e.g. Swartz et al. 2008, Kovlakas et al. 2020). The scaling
relations linking ULX rates with galaxy mass, star-formation rate and
metallicity were even calibrated in recent works (Anastasopoulou et al., 2019;
Kovlakas et al., 2020; Lehmer et al., 2021).

Figure 8 shows the evolution of ULX rates $n_{ULX}$ with various galaxy
parameters. The qualitative trends reported in the literature are also present
in our cleaned ULX sample:

  * •

ULX rates globally increase with the galaxy stellar mass, regardless of the
star-formation rate; however an excess is visible for dwarf galaxies
($10^{8}<M_{gal}/\mathrm{M}_{\odot}<10^{9}$)

  * •

ULX rates increase with the star-formation rate, at least in spiral galaxies

  * •

ULX rates are higher in spiral galaxies, regardless of the star-formation
rate; a significant ULX population does exist in elliptical galaxies

Modelling this ULX rate – SFR relation with a powerlaw,
$\log(n_{ULX})=\alpha\log(\mathrm{SFR})+\beta$ where $\alpha$ is the slope and
$\beta$ the normalization, we obtain $\alpha=0.43\pm 0.04$ for spiral galaxies
shown in the middle panel of Figure 8. This is in excellent agreement with the
value obtained by Kovlakas et al. (2020) on average over all types of spiral
galaxies: $\alpha=0.45^{+0.06}_{-0.09}$. Elliptical galaxies seem to present a
sharper scaling relation, with a slope $\alpha=0.84\pm 0.07$ at
$\mathrm{SFR}<1\leavevmode\nobreak\ \mathrm{M}_{\odot}$ yr-1 and a drop above
this value. However, the SFR estimator we used is unreliable for elliptical
galaxies, being degenerate with the dust mass, so this result is inconclusive.
In the top panel of Figure 8, we can see an excess of ULXs in dwarf galaxies
and a tenfold increase of the ULX rate over five orders of magnitude in
stellar masses covered by the sample: the specific ULX rate $n_{ULX}/M_{*}$
thus decreases with stellar mass, in agreement with Walton et al. (2011) and
Kovlakas et al. (2020).

Figure 8: Rates of ULX as a function of different properties characterising
galaxy environment.

We find tentative evidence that ULXs and HLXs do not share exactly the same
environment: the latter tend to be hosted in equal rates in spiral and
elliptical galaxies (Figure 9). Note that this trend was reported in Bernadich
et al. (2022), but not found in Kovlakas et al. (2020). This trend can be
interpreted in terms of the mean slope of the XLF of ULXs, which (in our
sample) is shallower for elliptical galaxies than for spiral galaxies. This is
in tension with the trend reported in the literature (e.g. Wang et al. 2016).
This point is discussed in Section 4.4.2.

Figure 9: Evolution with X-ray luminosity of the ratio between ULX rates
$N(L_{X})$ in spiral and elliptical galaxies.

The radial distribution of ULXs in their host also gives constraints on their
favourable environments. However, it cannot be studied in galaxies of small
angular size or large distance, due to source confusion issues and the
exclusion of ULXs close to the nucleus. We thus examine the number of ULXs in
spiral and elliptical galaxies having a major axis $D_{25}>30$ arcsec, as a
function of their galactocentric distance. The result is shown in Figure 10,
where the first bin is not shown since the ULX census in this inner part of
the galaxy ($\lesssim 0.15\times D_{25}$) is incomplete. Regardless of the
sample in use, the surface density of ULXs in spiral galaxies smoothly
decreases towards larger separations, while it becomes flatter in elliptical
galaxies before being cut off beyond $D_{25}$. In the subsample of spiral
galaxies, containing enough ULXs to probe several luminosity bins, the form of
the distribution seems independent of the source luminosity. This is in
agreement with the results of Kovlakas et al. (2020), performing the same
analysis in galaxies closer than 40 Mpc. However, our ULX selection extending
beyond the $D_{25}$ ellipse as well as our filtering of contaminants allow us
to observe for the first time the drop in ULX density in the outer part of
elliptical galaxies.

Figure 10: Radial distribution of ULXs in elliptical (top) and spiral (middle,
bottom) galaxies, showing their surface density as a function of
galactocentric distance. The last panel shows the surface density of different
luminosity bins in the compilation of Chandra, XMM-Newton and Swift selected
ULXs.

###  3.3 Comparison of XRB and ULX

It is now widely accepted that most ULXs below $\sim 10^{41}$ erg s-1 contain
either neutron stars or stellar-mass black holes accreting above the Eddington
limit. However, the geometry of their accretion flow, as well as the
mechanisms explaining their spectral states are still poorly understood. From
the few nearby individual ULXs that have been observed multiple times with
large exposures, two decades of detailed spectral study led to the results
cited in the introduction.

For the three X-ray catalogues studied in this work, XRB are well-detected
only in nearby galaxies ($D\lesssim 20$ Mpc) which have larger apparent sizes
and thus overlap numerous background sources. For this reason, and because
they correspond to a source type actually present in the training sample, we
select them using a different criterion to avoid background and foreground
contaminants: sources following $P_{XRB}>0.7$, $sep<0.8$ and $L_{X}<10^{39}$
erg s-1 make up our sample of XRB candidates. The contamination rate among
these selected candidates is 10–12% from manual inspection.

####  3.3.1 Hardness

Figure 11 shows three hardness – hardness diagrams comparing ULX and XRB from
each survey. In all three samples, their distributions are mostly overlapping,
consistent with previous studies. Except for the minority of nearby ULX
sources observed with long exposure times, this hardness information is not
sufficient to probe differences in the spectral shape. We also plot ULX
candidates classified as background ($P_{AGN}>P_{Soft}$) and soft
($P_{AGN}<P_{Soft}$) contaminants. The locus of sources classified as
background contaminants is significantly offset from these populations, and as
expected, it matches precisely the locus of AGN. This gives further credit to
their classification as contaminants, and is in agreement with previous
studies on AGN colours. Soft contaminants are fewer and dispersed at different
loci of this parameter space.

More conclusive results may be found by exploring the hardness evolution when
the flux varies. Indeed, hardness-luminosity studies of black hole XRBs have
shown a hysteresis cycle between two canonical states, high/soft and low/hard
(Shakura & Sunyaev, 1973; Remillard & McClintock, 2006) while neutron star and
some black hole ULXs globally harden when they brighten (Kaaret et al., 2017;
Gúrpide et al., 2021). In ULXs, a softening occurring when the ULXs brighten
may indicate the presence of a black hole (e.g. Narayan et al. 2017), as found
in black hole XRBs; however this does not apply to at least some ULXs (e.g.
NGC 5907 ULX1, Gúrpide et al. 2021). We thus looked at the hardness evolution
of a sample of XRB and ULXs having varied significantly between detections.
While such a study would require a detailed spectral modelling of each source,
we implement a simple approach to identify sources showing large hardness
variations. In order to combine the detections coming from all three
instruments, we built a flux hardness ratio $HR67$ between the bands 0.2–2 and
2–12 keV:

$HR67=\frac{F_{2-12}-F_{0.2-2}}{F_{2-12}+F_{0.2-2}}$

This corresponds to existing XMM-Newton energy bands, and an extrapolation of
Swift and Chandra energy bands, for which we used fixed conversion factors. To
this end, although this is a severe approximation, we assume an absorbed
powerlaw model with ($\Gamma=1.7$, $n_{H}=3\times 10^{20}$ cm-2) as already
done in CSC2, 4XMM and part of Swift catalogues to convert count rates into
fluxes. From the ULX and XRB samples in all three surveys, we select 59 unique
sources that varied by a factor of at least 10 during their follow-up, and we
consider their 799 detections that have $HR67\\_err<0.25$. This selection
targets high S/N sources having a good follow-up, thus most of them are
already well-studied (32 have an entry in Simbad). Sources having their peak
flux below (resp. above) $10^{39}$ erg s-1 are considered as XRB (resp. ULX).
Figure 12 shows the detected luminosities of these 59 sources, with hardness
being colour-coded.

While a subsequent fraction ($\sim 50\%$) of sources do not present a
significant hardness evolution ($\Delta HR<0.2$), most of the other sources –
which have variable hardness – follow the expected trend: from the lower
luminosity state to the higher luminosity state, a significant fraction of
XRBs become softer and most ULXs harder. A few outliers are to be noted: some
ULXs become softer instead, such as M101 ULX-1 (ObjID=38) which is a well-
known supersoft ULX thought to host a black hole (Liu et al., 2013; Shen et
al., 2015), or ESO 243-49 HLX-1 (ObjID=59) having a well-studied high/soft
state (e.g. Servillat et al. 2011; Godet et al. 2009, 2012). NGC 4490 ULX-3
(ObjID=57) has a very low HR at peak (Swift) detection, which is unreliable
($L_{X,det}>10^{41}$ erg s-1 ), probably due to its confusion with the ambient
hot gas and the vicinity of ULX-1 which is softer – indeed the Swift
coordinates are offset in this direction, by 5 arcsec from the Chandra source.
While 4XMM J095524.8+690113 (ObjID=24) is always very soft, it matches SN
1993J. NGC 5907 ULX-1 (ObjID=55) and 4XMM J022239.1+422328 (ObjID=25) are
always very hard (this is confirmed by an inspection of their XMM-Newton
spectra), as expected from their locations right inside a gas-rich, edge-on
host galaxy suggesting an important absorption. This sample of 59 sources
represents a large sample of uniformly selected, highly variable extragalactic
XRBs and ULXs.

Figure 11: Hardness–hardness diagrams of XRB and ULX candidates from 4XMM (top
panel), CSC2 (middle panel) and 2SXPS (bottom panel). Contours are shown to
ease visualisation. The greyscale density in the background refers to known
AGN. Candidates classified as foreground and background contaminants are shown
as orange circles and pale blue contours, respectively. The hardness ratios
plotted here correspond to the following energy ranges (in keV): $0.5-1/1-2$
(SC_HR2), $1-2/2-4.5$ (SC_HR3), $0.5-1.2$/$1.2-2$ (hard_ms), $1.2-2/2-7$
(hard_hm), $0.3-1/1-2$ (HR1) and $1-2/2-10$ (HR2). Figure 12: Luminosity-
hardness evolution of the subset of 59 variable XRB and ULX. $HR67$ is the
hardness ratio between $0.2-2$ and $2-12$ keV bands. All detections having
$HR67\\_err<0.25$ are shown. To ease visualisation, for each source, the
correlation, anti-correlation, or absence of correlation between
$L_{X,detection}$ and $HR67$ are highlighted at the bottom of the graph by the
”+”, ”$-$” and ”=” signs, respectively.

####  3.3.2 Variability

The variability within an observation is hard to characterise on a systematic
basis, as it takes various forms depending on the source type and is strongly
affected by the binning. Additionally, automated processing of the detailed
light curve is time consuming and subject to biases. However, some variability
indicators are made available as columns of the X-ray catalogues: for
instance, the fractional variability $F_{var}$ represents the dispersion of
the flux between snapshots, weighted by the flux error, and is given in 4XMM-
DR11. It can be used to isolate the most significantly variable objects. The
probability that the source is constant within the observation is available as
well in all 3 catalogues, using various tests. 2SXPS records the maximum rate
among snapshots in each band, which is valuable to spot short-term outbursts.
For each 2SXPS source, we compute the ratio between the maximum and the median
rate among all snapshots in the most variable band. Since Swift snapshots are
shorter (20-30 minutes) than the typical duration of our multi-mission
observations ($\sim 10$ks), this indicator probes variability on shorter
timescales.

We find no significant difference between the two populations for XMM-Newton
and Chandra variability indicators. However, 2SXPS ULXs seem less variable
than XRBs between observations and (especially) snapshots (Figure 13, y-axis).
On average, the former varied by 0.71 dex between snapshots, while the latter
varied by 1.15 dex (the standard deviation are 0.4 and 0.5 dex, respectively).
This difference may be intrinsic, due to the physics behind these populations;
but it may also be an observer bias, because unlike other telescopes, the
Swift monitoring of an XRB is generally performed to track its particular
variability. To remove this doubt, we also computed the variability ratio
between snapshots for XMM-Newton sources, by rebinning the 4XMM light curves
to 20 minute long bins, and recording the peak rate of each band. The result
is shown in the bottom panel of Figure 13: the same offset is visible, with
ULXs varying of 0.56 dex between snapshots on average and XRB varying of 0.93
dex.

However, this indicator is subject to at least two biases: first, the exposure
time, because episodes of high flux variations are more likely to be detected
in longer exposures. XRB have to be at lower distances to be detected, in
particular in better-studied galaxies, targeted by longer or more numerous
observations. This is why they are detected at lower fluxes in our samples. In
the Swift sample, 112 out of 293 XRB are located in just 4 galaxies observed
for more than 200 ks: M51, M81, M101 and NGC 300. After removal of these
galaxies, the ratio between the mean variability of XRBs and ULXs is reduced
by 30%. The second important parameter is the flux, because the noise
contribution to the source flux is more significant at low flux and equal
exposure. At equal fluxes and exposure time, both populations have similar
variability distributions in both surveys. This is tentative evidence that
both populations undergo a short-term flux variability of equal amplitude,
however this result would need a detailed treatment of the light curves to be
confirmed, which is beyond the scope of our study.

Figure 13: Variability between snapshots of XRB and ULX candidates from 2SXPS
(top panel) and 4XMM (bottom panel), as a function of their flux. Contours are
shown to ease visualisation.

###  3.4 The hyperluminous X-ray sample

In our work, thanks to the high completeness of GLADE even at several hundred
Mpc, we are able to retrieve 13, 115 and 75 HLX candidates from Swift, XMM-
Newton and Chandra samples, respectively. Of these 191 unique candidates, 76
sources have a bright optical counterpart and are split in two categories:
“galaxy pair candidates” (22 objects), if the source matches the nucleus of a
known galaxy in interaction with the host; and “weak HLX candidates” (54
objects), for other point-like optical sources or sources considerably offset
from the galaxy area . The latter could indeed be background AGN, as suggested
by the identification of 4 such sources as mid-infrared-detected AGN (Secrest
et al., 2015). However, they could also be dwarf satellite galaxies of the
host, which are a favourable environment to look for IMBH (Webb et al., 2010;
Greene et al., 2020; Barrows et al., 2019); and actually 7 of them have
photometric redshifts consistent with the host distance. The 115 remaining
HLXs with no optical counterpart, or a faint counterpart overlapping the area
of the galaxy, are qualified as “robust HLX candidates”. Figure 14 shows a
random selection of 15 candidates belonging to one of these three categories,
taken from the three X-ray surveys. For clarity, we summarise the number
counts of each subsample of initial and selected HLX candidates in Table 8 .

The distributions of the distances and mean observed 0.5–10 keV X-ray
luminosities of the 191 HLX candidates are shown in Figure 15. As usually
found, HLXs are detected much further away than ULXs. Most of them have
luminosities in the range $10^{41}-10^{42}$ erg s -1, with galaxy pair
candidates being the most luminous. By construction, all candidates have
$S/N>3$. However, we note that only 17% of the sample of 169 (robust+weak)
candidates have $S/N>10$. The median $S/N$ for this sample is 4.9.

Eight candidates that were removed, but are present in the literature are
discussed here. Some sources are simply not present in our sample: NGC 2276
ULX-1 is located outside the GLADE $1.26\times D_{25}$ ellipse, which does not
fit the actual galaxy area (this is the case for a small proportion of GLADE
galaxies). Reported in Zolotukhin et al. (2016) as a reliable candidate,
XMM0838+24 is here located outside the extent of its GLADE association, its
separation being 1.5 times the galaxy radius at its position angle. The other
reliable candidate cited by this study, XMM1226+12, was discarded during
manual inspection as lying visibly outside the extent of the host. M82 X-1 is
flagged as confused, or extended, in all three surveys. Thus these two sources
are not present in our sample. The candidate in IC 4320 (an elliptical galaxy
at $\sim$93 Mpc), proven to be a background AGN (Sutton et al., 2015), was
discarded due to a PanSTARRS photometric redshift higher than the host
($z_{ph}=0.18\pm 0.06$, Tarrío & Zarattini 2020). Because we use the mean
luminosity in each survey for selection criterion, the candidates located in
NGC 5907, NGC 4077, UGC 6697 and Cartwheel are not in our HLX sample, due to
their $\langle L_{X}\rangle<10^{41}$ erg s-1. We still retrieve some well-
known candidates in our work: they include ESO 243-49 HLX-1, NGC 470 HLX-1
(Walton et al., 2011; Sutton et al., 2012) and 3XMM J161604.0-223726 in IC
4596 (Earnshaw et al., 2019).

Figure 14: Images of a random sample of new HLX candidates. The images are
structured as follows: column (1,2), (3,4) and (5) display CSC2, 4XMM-DR11 and
2SXPS candidates, respectively. (odd numbers) and (even numbers) rows show the
X-ray and optical image at the same scale – the white line corresponds to 10
arcsec. Rows (1,2), (3,4) and (5,6) correspond to “galaxy pair”, “robust” and
“weak” candidates, respectively. These candidates have the following names:
1–2CXO J124208.4+331854, 2–2CXO J025921.5+132913, 3–4XMM J123441.0+020846,
4–4XMM J104444.7-012018, 5–2SXPS J164652.0+234011, 6–2CXO J082215.9+210535,
7–2CXO J105210.4+552243, 8–4XMM J062447.9-372122, 9–4XMM J125708.6-044144,
10–2SXPS J111416.1+481833, 11–2CXO J005151.7+474019, 12–2CXO J131133.3-011656,
13–4XMM J083235.2-225804, 14–4XMM J145753.6-113959 and 15–2SXPS
J115109.2+570340. X-ray images correspond to the Chandra and XMM-EPIC images
displayed in ESASky, or the science exposures of Swift-XRT where the green
circle has a 10 arcsec radius. Optical images come from the PanSTARRS survey
(Chambers et al., 2016), except candidate 8 (Digitized Sky Survey, McLean et
al. 2000) as it lies in a region not covered by PanSTARRS. The 3$\sigma$ X-ray
position error circle is shown in cyan. Table 6: Subsamples of HLXs under
study

| CSC2 | 4XMM | 2SXPS | Total (unique)
---|---|---|---|---
Initial candidates | 195 | 360 | 110 | 665 (619)
With redshift | 96 | 187 | 41 | 324 (298)
Background | 82 | 165 | 35 | 282 (259)
Foreground | 2 | 3 | 0 | 5 (5)
Distance match | 12 | 19 | 6 | 37 (34)
Selected candidates | 75 | 115 | 13 | 203 (191)
Robust | 51 | 66 | 5 | 122 (115)
Weak | 15 | 35 | 7 | 57 (54)
Galaxy pair | 9 | 14 | 1 | 24 (22)
With opt. counterpart | 29 | 63 | 8 | 100 (94)
Soft ($HR67<-0.6$) | 2 | 5 | 1 | 8 (6)
Hard ($HR67>-0.6$) | 67 | 110 | 12 | 189 (179)

888Number counts of initial and selected HLXs in different subsamples (see
Section 2.7 and Section 3.4, respectively, for details on the subsamples). The
sample of “soft HLXs” is defined in Section 3.5.1, and is the complementary of
“hard HLXs” (note that 6 Chandra candidates are not included in these samples,
because their $HR67$ could not be calculated).

Figure 15: (Left) Distribution of distances for the selected samples of ULXs
and HLXs. ULXs are detected up to 500 Mpc and HLXs up to 1100 Mpc. (Right)
Distribution of the 0.5–10 keV observed mean X-ray luminosity of the 191 HLXs
candidates, including the 22 galaxy pair candidates.

###  3.5 Comparison of ULXs and HLXs

Because of the low number of known HLXs, few studies have attempted to
characterise their typical environment. The location of ESO 243-49 HLX-1
outside the disk of its host may be due to it being embedded in a stellar
cluster or the stripped core of a dwarf galaxy (e.g. Webb et al. 2010) but ESO
243-49 HLX-1 in itself seems to be an outlier among HLXs (e.g. Sutton et al.
2015). Other HLX discoveries led to the suggestion that HLXs and ULXs share
the same type of environment, namely star-forming, spiral galaxies (Sutton et
al., 2015; Barrows et al., 2019). This is at odds with the ratio of our HLX
rates between spiral and elliptical (Figure 9), suggesting an equal rate of
robust HLXs in these two morphologies. We compare the behaviour followed by
ULX and HLX rates as a function of three environmental parameters (Figure 16):
the galaxy stellar mass, its star-formation rate and its Hubble type. Unlike
ULXs, our HLXs tend indeed to be hosted equally in (regular) spiral and
elliptical galaxies. However we note that only 38 HLXs have a Hubble type
value provided in HyperLEDA, and that 40% of them have errors on this type
$e_{t}\geq 3$ (equivalent to 1.5 major graduations in Figure 8). On the other
hand, similarly to the trends obtained on ULXs, we find a positive correlation
between the HLX rate and the stellar-mass and star-formation rate.

Figure 16: Rates of ULX and HLX as a function of different properties
characterising galaxy environment.

####  3.5.1 Hardness

We investigate the locus of the 191 HLX candidates in a hardness–hardness
diagram. Figure 17 shows the diagrams obtained from XMM-Newton, Chandra and
Swift data, after discarding sources with $HR=-1$ or 1 to limit unreliable
measurements. The HLX and ULX distributions are mainly overlapping, with only
a fraction of XMM-Newton HLX candidates looking harder above 4.5 keV
($SC\\_HR4$). They consist mainly of 15 robust HLX candidates having
$SC\\_HR4>0$ on average among detections; their low signal-to-noise ratio
($S/N<10$) does not allow a more in depth spectral study. However, the energy
range covered by $SC\\_HR4$ essentially probes the photon index of the
spectrum if it were an absorbed powerlaw spectrum, and in this case such
values would only be produced by $\Gamma\lesssim 2$.

Besides, in the Swift sample, ESO 243-49 HLX-1 is an outlier located far in
the lower left of the distribution. Three other outliers are found in 4XMM:
4XMM J215022.4-055109, a well-known IMBH candidate which underwent a tidal
disruption event (Lin et al., 2018), 4XMM J161534.3+192707, a variable HLX
possibly consistent with another offnuclear tidal disruption event (Tranin et
al., in prep), and 4XMM J085253.8+180110, a weak candidate located in a
cluster of galaxies. Another soft outlier (not visible in Figure 17 because of
having $SC\\_HR4=-1$) is 4XMM J231818.7-422237, a variable soft source already
cited as a candidate IMBH in Lin et al. (2014). It could actually be
associated with two GLADE galaxies, and would be of ULX luminosity if
belonging to the one closer to us. It is however softer than all other ULXs
($SC\\_HR2=-0.73$, $SC\\_HR3=-0.85$). The outlier nature of such soft HLXs is
even more visible in the distribution of $HR67$, the hardness ratio between
bands $0.2-2$ and $2-12$ keV, where they form the peak close to -1
($HR67<-0.6$) in Figure 18. All but two of them are also seen as outliers by
the source classification, with an outlier measure $>12$ (Tranin et al.,
2022), higher than 99% of ULX candidates.

Figure 17: Hardness–hardness diagrams of ULX and HLX candidates from Chandra
(top), XMM-Newton (middle) and Swift (bottom). Grey and red dots correspond to
ULXs and HLXs, respectively. Contours at 25th, 50th and 75th percentiles are
shown to ease visualisation. Figure 18: Distribution of the hardness ratio
$HR67$ between the 0.2-2 and 2-12 keV energy bands for selected HLX
candidates.

####  3.5.2 Variability

Some extreme ULXs were seen to approach or overcome the luminosity threshold
of HLXs for transient periods. Super-Eddington radiation well above the
Eddington limit seems possible for the timescale of an observation (Israel et
al., 2017), and in this case the Eddington limit cannot give accurate mass
estimations of the central object. HLX candidates detected only once could
thus be powered by stellar-mass accretors. However, we find that out of the
191 HLX candidates, 96 were detected several times (up to 20 times, and even
more for ESO 243-49 HLX-1, NGC 470 HLX-1 and the AGN in a galaxy pair XMMU
J134736.6+173403). Only 26 of them have a detection below $10^{41}$ erg s-1 –
including the known candidates in ESO 243-49, IC 4596 and NGC 470 – suggesting
that a majority of our HLX sample is actually persistent. The time between the
first and the last observation is $>1$ yr for 75% of persistent candidates.
Figure 19 shows the distribution of their variability ratio between
observations, along the one of ULXs. No significant variability excess is seen
in either of these two categories.

Figure 19: Distribution of the variability ratio between observations
($logFratio=\log\big{(}\frac{F_{max}}{F_{min}}\big{)}$), for ULXs and HLXs.

##  4 Discussion

###  4.1 Comparison to other studies

Recent releases of large X-ray surveys and complete galaxy catalogues enabled
the discovery of large samples of ULXs and HLXs. We discuss our findings in
the context of previous studies on this subject.

Kovlakas et al. (2020) found 629 ULX candidates in 309 galaxies in the local
40 Mpc volume, from Chandra CSC2 catalogue matched with HECATE (Kovlakas et
al., 2021). They estimated a contamination rate of 20%, using the $\log N-\log
S$ method. Because of the high Chandra sensitivity, their sample is likely to
be complete. In our work, despite the use of a different galaxy catalogue and
slightly different selection criteria, we find a similar number of Chandra ULX
candidates (567) in galaxies closer than 40 Mpc. However, we remove 58
candidates of this sample from their classification as contaminants, 16 others
from their identification as another type in the literature, and 41 are also
flagged during manual inspection. This leaves a sample of 442 sources located
within 40 Mpc, representing 78% of the initial sample. Adding Swift and XMM-
Newton candidates, the size of the sample is almost doubled (839). This
outlines the need to cross-correlate several X-ray catalogues to find these
sources, since focusing on one catalogue decreases the survey coverage and
misses some candidates as they may have been fainter than $10^{39}$ erg s-1 at
the time of the observation. Besides, we were able to confirm key results
present in their work, namely the predominance of contaminants in uncleaned
samples at luminosities $L_{X}>10^{40.5}$ erg s-1, the evolution of the ULX
rate with SFR and the decrease of specific ULX rate with galaxy stellar mass.

Bernadich et al. (2022) used a recent release of the 4XMM catalogue, 4XMM-
DR10, matched with HECATE galaxies, to build a clean ULX sample. Similarly to
our work, their cleaning pipeline relies on removing the known AGN and
foreground contaminants, as well as manual inspection of the sources having
optical counterparts (from PanSTARRS and Gaia). However, without an automatic
classification, they were not able to remove contaminants having
$F_{X}>F_{Opt}$ or no optical counterparts. Besides, we were able to improve
this cleaning by manual inspection of faint counterparts in the southern
hemisphere, using the DES survey. Overall, their sample of 779 candidates is
somewhat larger than our cleaned sample (667 candidates). They however
consider a source as ULX as soon as $L_{X}+L_{X,err}>10^{39}$ erg s-1 in at
least one detection, and only 632 of their candidates have $\langle
L_{X}\rangle>10^{39}$ erg s-1.

Their work also focuses on the statistical study of bright ULXs, having
$L_{X}+L_{X,err}>5\times 10^{40}$ erg cm-2 s-1 in at least one detection.
While this bright sample contains 94 such sources, it is significantly
different from our sample of HLXs for several reasons: first, because we use
the mean luminosity to select them, and because we use the more conventional
threshold of $10^{41}$ erg $s^{-1}$ (only 25 of their bright candidates have
$\langle L_{X}\rangle>10^{41}$ erg $s^{-1}$). A large fraction of their sample
is thus likely to be similar to luminous ULXs powered by a stellar-mass
accretor. Second, because we use GLADE, which is more complete than HECATE as
it is a compilation of catalogues (including HyperLEDA), and because HECATE is
limited to galaxies at $D<200$ Mpc, our HLX sample contains significantly more
candidates. As a result, we find significant differences in the preferred
environments and counterparts of ULXs and HLXs (Section 4.5). We do not find
evidence in our sample for the softer average spectra for HLXs reported in
Bernadich et al. (2022).

Walton et al. (2022) is similar to our work in the sense that they included
the same three X-ray surveys. However, their galaxy sample is different,
compiling HyperLEDA and the Catalogue of Neighbouring Galaxies. While they
used the same matching limit as we do, $1.26\times D_{25}$, the contamination
rate of their sample must be high (Figure 6) given that they only cleaned out
known contaminants. Their final sample contains 1843 unique candidates,
compiled from the 4XMM-DR10, 2SXPS and CSC2 samples containing 641, 501 and
1031 sources, respectively. This total is similar to the 1901 ULXs contained
in our compilation, with however slightly different proportions from each
survey (667, 304 and 1185). Our slightly larger sample, despite our strict
cleaning, was possible thanks to the use of GLADE. Our sample thus constitutes
both the largest and the cleanest ULX catalogue to date.

We examine the classification of X-ray sources that we have in common with
these three studies. We find 601 sources of Kovlakas et al. (2020) in our
sample, of which only 423 pass the offnuclear and $\langle
L_{X}\rangle>10^{39}$ erg s-1 criteria. The rest are mostly discarded because
of insufficient luminosities, due to their use of CSC2 fluxes from the PSF
90%-energy width ($flux\\_aper90\\_b$) instead of the source region aperture
($flux\\_aper\\_b$). 359 sources are selected as ULXs after classification and
manual screening, hence $\sim$15% of candidates are discarded as contaminants.
Similarly, in the clean ULX catalogue of Earnshaw et al. (2019) claiming a
contamination rate of 24%, we find 17% of candidates removed in our pipeline.
From the ULX sample of Walton et al. (2022), 1734 sources are included in this
study, of which 1291 pass the first criteria (their luminosity criterion being
less stringent). 27% of these sources are discarded as contaminants, higher
than their estimated background contamination rate (18–23%). Most of the ULXs
of Bernadich et al. (2022) are present in our sample (746), and 574 have ULX-
consistent positions and luminosities in our sample. However, 27% of these
sources are discarded by the classification and manual screening, considerably
higher than expected from their contamination estimation ($\sim$2%). Their
estimation is however calculated for their entire X-ray sample (including non-
ULXs), perhaps explaining this discrepancy.

Last but not least, Barrows et al. (2019) produced a catalogue of 169 HLX
candidates from the CSC2 catalogue and SDSS images. This unprecedented sample
size is allowed by their detailed treatment of optical images to detect
galaxies up to greater distances and fainter surface brightness; as a result,
most of their candidates are located at distances $>1000$ Mpc, and present
X-ray luminosities $>10^{42}$ erg s-1. Few of their galaxies are present in
the GLADE catalogue: as a result, only 10 of their candidates are present in
our CSC2 sample. Two of them have $\langle L_{X}\rangle<10^{41}$ erg s-1, and
one is discarded when excluding the galaxy central region. Interestingly, 2CXO
J122148.2+491131 is flagged as background contaminant due to its photometric
redshift $z_{ph}=0.46\pm 0.1$ (Beck et al. 2016, also consistent with
PanSTARRS and DESI estimates), leaving only 6 candidates out of the possible
10 in our sample. Whilst their study has the advantage of a low contamination
rate (7–8%) and probes the HLX population to greater distances in the
Universe, our sample is more likely to contain candidates well-suited for a
follow-up as well as HLXs of lower luminosity, probing lower IMBH masses and
possibly earlier stages of black hole growth. These points are further
discussed in Section 4.5.

###  4.2 Number of contaminants

The cosmic X-ray background (CXB) has been extensively studied for several
decades by several X-ray missions, revealing that it is essentially due to
radiation from unresolved point-like sources, in particular AGN (e.g., Moretti
et al. 2003 and references therein). As such, its flux distribution has often
been used in ULX and HLX studies to assess the number of background/foreground
contaminants behind a given galaxy. For instance, the $\log N-\log S$ relation
of Moretti et al. (2003), gathering wide-field and pencil-beam surveys probing
six orders of magnitude in flux, was used by Walton et al. (2011), Zolotukhin
et al. (2016) and Barrows et al. (2019) to compute their number of
contaminants. To this aim, they converted the broad band flux of their sample
to the hard energy band of Moretti et al. (2003) (2–10 keV) which yields the
most contaminants, using an absorbed powerlaw spectrum with $\Gamma=1.69$ &
$n_{H}\sim 3\times 10^{20}$ cm-2. However, a possible caveat of this method is
that CXB sources present a variety of spectra (Mushotzky et al., 2000) that
cannot be accurately represented by this model. To alleviate this problem, we
also compute the number of contaminants using $\log N-\log S$ distributions
calibrated in an energy band similar to the broad bands of our sample, namely
0.2–12 keV (XMM-Newton), 0.3–10 keV (Swift) and 0.5–7 keV (Chandra). This is
the case of the relations of Georgakakis et al. (2008) and Kim et al. (2007),
respectively used in the ULX studies of Zolotukhin et al. (2016) and Kovlakas
et al. (2020). After converting our fluxes in the proper energy band, we find
excellent agreement between the results obtained from each relation.

The number of contaminants in each galaxy is computed by multiplying the
cumulative flux density, $\log(N(>S))$ at a given flux $S$, by the extent of
the galaxy used for finding ULXs (i.e. after the centre exclusion criterion
has been applied). The flux to consider is directly given by the galaxy
distance (such that the contaminant is found in the ULX luminosity range) and
the X-ray sensitivity (so that the contaminant is detectable), as
$\min(f_{39},f_{lim})$ where
$f_{39}=\frac{10^{39}\text{erg\leavevmode\nobreak\ s}^{-1}}{4\pi D^{2}}$.

Because our centre exclusion radius $3(POSERR+0.5)$ depends on the source
position error, we estimate the number of contaminants using the 10%, 50% and
90% percentiles of the position error distribution in each sample of
candidates. Results are given in Table 9 for the ULX and HLX contaminants,
where error bars correspond to the variations in position error. For ULXs, it
is comparable to the number of candidates known or classified as AGN, or
removed by manual inspection of the optical source. For HLXs, it is similar to
the number of background sources (using photometric redshifts) added to the
number of HLX candidates flagged as “weak”. This is a further indication that
the selection of ULX and HLX candidates, and thus the automatic classification
and the photometric redshifts estimates, are reliable and robust.

Except in the 4XMM sample, our filtering process identifies more ULX
candidates as background contaminants than expected from the CXB. Using the
CXB curves of Moretti et al. (2003), yielding the larger contamination rates,
the 41 remaining background contaminants in our sample thus represent 2.2% of
the sample of selected ULXs.

Table 7: Number of contaminants estimated using different published $\log
N-\log S$ relations, in the ULX and HLX samples.

Survey | Moretti+03 | Georgakakis+08 | Kim+07 | Identified
---|---|---|---|---
2SXPS ULXs | $140^{+7}_{-15}$ | $141^{+7}_{-14}$ | $139^{+7}_{-15}$ | 171
2SXPS HLXs | $13^{+2}_{-6}$ | $13^{+3}_{-5}$ | $13^{+2}_{-6}$ | 25
4XMM ULXs | $582^{+55}_{-73}$ | $558^{+53}_{-70}$ | $554^{+53}_{-69}$ | 541
4XMM HLXs | $136^{+34}_{-44}$ | $137^{+33}_{-44}$ | $135^{+33}_{-43}$ | 121
CSC2 ULXs | $483^{+38}_{-67}$ | $451^{+35}_{-62}$ | $452^{+35}_{-62}$ | 482
CSC2 HLXs | $75^{+17}_{-31}$ | $73^{+17}_{-30}$ | $73^{+17}_{-30}$ | 64

999Positive and negative errors correspond to the estimates at 10% and 90% of
the exclusion radius ($3(POSERR+0.5)$) distribution. The last column gives the
number of candidates filtered out by our pipeline as background contaminants,
or identified as HLX “weak candidates”.

###  4.3 Heterogeneity of the sample and instrumental biases

Because of the different capabilities of Swift-XRT, XMM-Newton and Chandra, we
do not expect to probe exactly the same ULX populations in each survey.
Chandra on-axis observations are more sensitive and have a higher angular
resolving power, so it can probe larger distances where the galaxies look
smaller. XMM-Newton and Swift have covered larger sky areas, and the superior
follow-up capabilities of Swift can also reveal populations of variable ULXs.
The design of Swift as a gamma-ray burst monitor has led to a more uniform sky
coverage, less biased towards galaxies of special interest, so we can expect a
lower observer bias of its ULX population. However, the centre-excluding
criterion (Equation (1)) leads to exclude a significant fraction of actual
ULXs at low galaxy angular separation. This is particularly the case of Swift
HLX candidates, which are all offset from their host centre by $>10$ arcsec,
while half of XMM-Newton and Chandra HLX candidates have lower separations.
This probably causes the steeper slope of the 2SXPS XLF above $L_{X}>10^{40}$
erg s-1.

Another bias in our selection of ULXs and HLXs is the significant difference
between Swift, XMM-Newton and Chandra energy bands: in particular, the
narrower energy band of the latter probably removes valid candidates just
below the ULX luminosity threshold. Assuming a phenomenological absorbed
powerlaw model for ULX spectra, with typical parameters $\Gamma=2.4$,
$n_{H}=2\times 10^{21}$ cm-2 (Gladstone et al., 2009), the conversion factor
from the XMM-Newton 0.2–12 keV band to the Chandra 0.5–7 keV band is 0.85, and
becomes as low as 0.7 for a photon index $\Gamma\lesssim 2$. If this
conversion is valid on average, the Chandra XLF break in spiral galaxies would
be shifted to up to $\sim 4\times 10^{39}$ erg s-1, better matching the other
XLFs.

Nevertheless, results obtained from the XLF and the comparison with XRBs, ULXs
and HLXs are comparable between instruments, showing the robustness of both
our cleaning pipeline and our statistical study with respect to instrumental
biases.

The XLF inferred from 2SXPS, 4XMM and CSC2 are compatible at the $2-\sigma$
level except at the faint end $L_{X}<5\times 10^{39}$ erg s-1 . At these
luminosities, most of the sample is made up of sources in galaxies in the
distance range of 5–60 Mpc: those sources are most often resolved by Chandra,
whereas a substantial fraction of XMM-Newton and Swift-XRT sources are
confused (Figure 1). Moreover, because of the higher noise level and the
poorer spatial resolution, sources hidden in diffuse X-ray emission are harder
to detect with those facilities. Consequently, the intrinsic slope of the XMM-
Newton and Swift XLF are substantially underestimated at the faint end.

A similar problem occurs at fainter luminosities for Chandra sources: as
explained in Wang et al. (2016), sources fainter than $10^{38}$ erg s-1 are
more likely to be below the detection threshold, in particular when located in
diffuse emission. This leads to the flattening of the XLF at these
luminosities, so that diagnostics on the XLF faint end of XRB must rely on
deep observations of close individual galaxies, as done in Mineo et al. (2012)
or include a completeness correction as in Wang et al. (2016).

Another moderately significant difference is the steeper slope in the bright
end of 2SXPS XLFs, above $\sim 2\times 10^{40}$ erg s-1, while 4XMM and CSC2
are consistent with each other. A closer look at each sample above this
luminosity reveals that half of 4XMM and CSC2 candidates are detected in
galaxies of apparent diameter $D_{25}<20$ Mpc, whereas 2SXPS candidates in
such galaxies are systematically discarded by the centre exclusion criterion
due to their larger position error. There is thus a lack of HLXs in the 2SXPS
sample. As a result, the XLF based on Chandra sources is probably the most
representative of the intrinsic, “universal” XLF of spiral galaxies.

###  4.4 X-ray luminosity function

####  4.4.1 XLF break

Our work shows significant evidence of the existence of a break between
$3\times 10^{39}$ and $10^{40}$ erg s-1 . This somewhat agrees with the
findings of Swartz et al. (2011) and Mineo et al. (2012), suggesting a cutoff
luminosity of $1-2\times 10^{40}$ erg s-1 . In contrast, Wang et al. (2016)
reported that this break is spuriously caused by the lack of statistics in
these studies, and stated that no break is needed at this point of their XLF.
We see two reasons for this inconsistency: first, their sample is made up by
searching X-ray sources within twice the $D_{25}$ of the galaxy, and they do
not implement any procedure to remove the abundant contaminants this implies.
These contaminants have an important influence at high luminosities (Figure 6)
and thus can flatten the XLF around $10^{40}$ erg s-1 . Second, they only
considered the differential luminosity function in their XLF fits, showing
higher dispersion, so that a break can be buried in this ”noise”. In their
Figure 9 showing the cumulative XLF, there might be an hint for a break around
$10^{40}$ erg s-1 .

The precise location of the break varies between samples, with 2SXPS sources
yielding the largest value. We argue that the lack of low-luminosity ULXs due
to source confusion issues, given the lower spatial resolution of Swift-XRT,
and of HLXs due to the large Swift exclusion radius, are likely to distort the
XLF and shift the break towards higher luminosities. This effect was already
pointed out by Wang et al. (2016) when considering more distant elliptical
galaxies, increasing the population of unresolved sources. Similarly, if ULXs
and HLXs have different nature (see Section 4.5), adding HLXs in the ULX
cumulative XLFs will distort it so that the break will shift towards fainter
luminosities and there will be a flattening at the high-luminosity end.

We test how the location of the break is changed when sources $>10^{41}$ erg
s-1 are removed before fitting the XLF: while it is unchanged for 2SXPS ULXs,
we obtain a best-fit break of $(5.7\pm 1.2)\times 10^{39}$ erg s-1 and
$(9.1\pm 2.6)\times 10^{39}$ erg s-1 for the CSC2 and 4XMM samples,
respectively. The single powerlaw fit is still excluded. Once converted to the
XMM-Newton broad energy band of 0.2–12 keV (assuming an absorbed-powerlaw
spectrum with $\Gamma=1.7$, $n_{H}=3\times 10^{20}$ cm-2), the CSC2 break is
further moved to $\sim 8\times 10^{39}$ erg s-1 , consistent with the break
fitted to the 4XMM sample.

More generally, although the use of cumulative XLFs provides higher statistics
and may reveal more features, it is also subject to more biases. In
particular, the faint end of the cumulative XLF is impacted by all the
selection biases, leading to its flattening, while having the smallest
poissonian error bars. This is likely to force the faint end of the XLF fit to
flatten, potentially shifting the fitted break luminosity towards fainter
luminosities to alleviate the tension between the model and the data. When the
differential XLF is fitted (in the form $dN/d\ln L_{X}$ so that the slopes of
the cumulative version are conserved), the break is shifted to $(6.4\pm
2.4)\times 10^{39}$, $(6.5\pm 4.2)\times 10^{39}$, $(1.21\pm 0.27)\times
10^{40}$ and $(7.9\pm 2.5)\times 10^{39}$ erg s-1 for the samples of CSC2,
4XMM, 2SXPS and all ULXs, respectively. As a result, the break is probably the
cutoff feature reported by Mineo et al. (2012) and Swartz et al. (2011) for
spiral galaxies. Most recently, this feature was also found in the XLF from
deep Chandra observations of the Virgo cluster (Soria et al., 2022), with a
normalisation jump occurring in the range $4-5\times 10^{39}$ erg s-1 .

We investigate the physical reasons for this break. The first possibility is
that a different galaxy population is probed in each luminosity bin of the
XLF. Indeed, the XLF shape was found to significantly depend on the galaxy
properties, with e.g. a flatter slope for more star-forming galaxies (Wang et
al., 2016). We find a moderate increase of the mean star formation rate of ULX
hosts in bins of increasing luminosity, with typical values $<0.5$ M⊙ yr-1
below $5\times 10^{39}$ erg s-1 and $>3$ M⊙ yr-1 above $10^{41}$ erg s-1.
Extrapolating the result of Wang et al. (2016) to higher luminosities, this
should produce a flattening of our XLF with luminosity, whereas we find the
opposite behaviour. The stellar mass follows the same trend as the SFR, so
that each bin probed similar specific SFRs. The XLF break is thus unlikely to
be a selection effect of different environments.

Another possibility is that the luminosity break corresponds to the Eddington
luminosity of a given class of objects, such as neutron stars or a certain
population of black holes. For example, the Eddington luminosity of $\sim
2\mathrm{M}_{\odot}$ neutron stars was argued to be consistent with the
$5\times 10^{38}$ erg s-1 break seen in the luminosity function of LMXB, when
accounting for beaming effects and the accretion of helium-rich material (Kim
& Fabbiano, 2004; Wang et al., 2016). With the same arguments, we could expect
the break found in this study to be due to the Eddington luminosity of $\sim
25$ M⊙ black holes. Recent modelling works show that radiative processes allow
a luminosity as high as $10^{40}$ erg s-1 from supercritically accreting
$20\mathrm{M}_{\odot}$ stellar-mass black holes (Krtička et al., 2022).
However, the $>50$ M⊙ black holes detected by gravitational wave facilities
rule out this break to match the maximum mass of stellar-mass black holes,
unless the mass distribution of stellar-mass black holes is bimodal (which
could be the case if black holes in XRBs and black holes in mergers do not
share the same formation channels, Fishbach & Kalogera 2022, but see
Belczynski et al. 2021). Moreover, the same luminosity could be reached by the
most magnetised neutron stars (Mushtukov et al., 2015), but the discovery of
pulsars reaching even greater luminosities (e.g. Israel et al. 2017) may have
challenged this interpretation.

On the other hand, the source spectrum, and hardness ratios, may allow better
diagnostics of the source nature. Using a sample of 17 ULXs including 6 PULXs,
Gúrpide et al. (2021) find that the secure neutron stars of their sample are
among the hardest sources. In order to probe which type of accretor dominate
at different luminosities, and if they all present the same XLF features, we
divide our samples in three bins of $HR67$: soft ($HR67<0.2$), medium
($0.2<HR67<0.5$) and hard ($HR67>0.5$) ULXs. The result is plotted in Figure
20, showing no significant difference between each XLF. The possible reason is
that $HR67$ is inefficient to distinguish neutron star and black hole ULXs, or
perhaps more interestingly, that a similar break should exist in both
populations even if for different physical reasons. For instance, since the
contribution of LMXB in spiral galaxies is known to be non-negligible even at
ULX luminosities (e.g. Mineo et al. 2012; Lehmer et al. 2019), the break could
be the manifestation of a normalisation jump in the XLF due to a luminosity
cutoff at a few $10^{39}$ erg s-1 for LMXBs, as mentioned in Soria et al.
(2022).

Alternatively, the HMXB XLF was found to present a large dispersion in the
range $10^{38}-10^{40}$ erg s-1 related to galaxy metallicity (Lehmer et al.,
2021), with a high luminosity cutoff located between $3\times 10^{39}$ erg s-1
(low metallicity galaxies, 12+log(O/H)$\simeq$7.6) and $10^{40}$ erg s-1
(galaxies of solar metallicity, 12+log(O/H)$\simeq$8.7 Allende Prieto et al.
2001). The metallicity information is absent from GLADE, but we find 106
galaxies in our ULX sample that have a match in HECATE with a metallicity
value: 85% of them are in the range 8.4$<$12+log(O/H)$<$8.9, close to solar
metallicity, supporting this explanation for the break. In this scenario, the
XLF break we obtained may be due to the combination of an excess of luminous
HMXB in low metallicity environments followed by a cutoff, and an increasingly
large contribution of another type of object at the bright end, possibly
intermediate-mass black holes, corresponding to the bulk of HLXs.

Figure 20: Deconvolved XLFs of Chandra, XMM-Newton and all selected ULXs in
spiral galaxies, after removal of HLX candidates. Three bins of hardness ratio
$HR67$ (see the text for details) are also shown, representing soft
($HR67<0.2$), medium ($0.2<HR67<0.5$) and hard ($HR67>0.5$) sources.

####  4.4.2 XLF slopes

Throughout this Section, we compare the best-fit slopes obtained on cumulative
XLFs to slope values reported in the literature. Since the latter are often
obtained on differential XLF, e.g. $dN/dL_{X}\propto L_{X}^{-\alpha}$ for a
simple powerlaw model, we convert these values to the cumulative case whenever
needed ($N(>L_{X})\propto L_{X}^{-\alpha^{\prime}}$ where
$\alpha^{\prime}=\alpha-1$).

As discussed in the previous Section, the fits to the cumulative XLFs are
subject to a number of biases which affect the fitted parameters, including
the XLF slopes. For example, the inclusion of HLXs flattens the bright end of
the distribution. Analysing the cumulative XLF with HLXs removed, we obtain
best-fit high-luminosity slopes ($\alpha_{2}$) of $1.82\pm 0.1$, $1.66\pm
0.13$, $1.79\pm 0.14$ and $1.80\pm 0.10$ for CSC2, 4XMM, 2SXPS and all ULXs,
respectively, in better agreement with the values inferred on the ULX-complete
galaxies and close to the value $1.73^{+1.58}_{-0.54}$ obtained by Mineo et
al. (2012) for HMXBs in the same luminosity range.

Besides, resolution effects are also likely to distort the XLF, flattening it
at the faint end when two low-luminosity ULXs are mistaken for a unique,
brighter source. This also explains why the number of ULXs per galaxy
decreases between CSC2, 4XMM and 2SXPS samples, with decreasing spatial
resolution. Similarly, given the centre exclusion criterion, the larger
position errors of Swift lead to the removal of some bona fide HLXs,
steepening the high-luminosity slope of the XLF.

Some of these biases can be alleviated by fitting differential XLFs, rather
than their cumulative version. The slopes we obtain are indeed different from
the cumulative XLF fits, with $(\alpha_{1},\alpha_{2})$ being $(0.64\pm
0.16,1.69\pm 0.17)$, $(0.88\pm 0.27,1.16\pm 0.12)$, $(0.58\pm 0.11,1.70\pm
0.18)$ and $(0.82\pm 0.09,1.38\pm 0.09)$ for CSC2, 4XMM, 2SXPS and all ULXs,
respectively. These low luminosity slope values better agree with the slope of
the HMXB XLF in Mineo et al. (2012) and Lehmer et al. (2019), however the high
luminosity slopes are not well-constrained.

Unlike spiral galaxies, the XLF we obtained for elliptical galaxies are in
disagreement with the steeper ones reported in the literature. For instance,
Walton et al. (2011) find a best-fit slope $\alpha=1.5\pm 0.4$ for their
elliptical galaxies, while Wang et al. (2016) obtain steeper slopes $1.80\pm
0.03$, $2.74\pm 0.10$ and $2.29\pm 0.18$ for elliptical, lenticular and all
early-type galaxies, respectively. Note that the definition of elliptical
galaxies in the present paper includes lenticular galaxies. In nearby
galaxies, the XLF of LMXBs is found to present a break at a few $10^{38}$ erg
s-1 followed by a powerlaw decrease of slope $\sim 1.8$ (Humphrey & Buote
2008; Kim & Fabbiano 2010).

The flatter slope of the XLF for elliptical galaxies in our study is examined.
One possibility is that the hot gas content, prominent in these galaxies, is
drowning the signal of ULXs well above the limiting luminosity we compute,
which does not depend on the background level. In that case the census of ULXs
is largely incomplete, up to a higher luminosity. For instance, M86 has a CSC2
limiting luminosity of $2\times 10^{38}$ erg s-1 from Equation (3) but the
only sources detected within the $D_{25}/2$ ellipse have $L_{X}>2\times
10^{39}$ erg s-1 and $S/N<3$. At fixed hot gas density, the hot gas flux in a
typical aperture increases with galaxy distance, affecting the detection of
increasingly bright sources. To test this effect, we fit the deconvolved XLFs
for samples with a range of distance limits. Figure 21 shows a clear
flattening of the XLF when extending the sample to further distances. Below 70
Mpc, the XLF of 4XMM ULXs becomes flatter than that of CSC2, probably because
of the lower spatial resolution of XMM-Newton. Even the closest elliptical
galaxies have few low-luminosity ULXs detected by Swift-XRT because of its
lower spatial resolution, hence the flat faint end of the XLF and the
inconsistent slope. For a distance limit of $\sim 40$ Mpc, as adopted in
previous studies of this XLF (e.g. Wang et al. 2016), we find a slope in the
range $1.6-1.8$, in good agreement with values reported in the literature.

Actually, Mineo et al. (2012) noted a large diversity of XLFs from different
spiral galaxies. This diversity of XLF is arguably due to different star
formation histories, metallicities and initial mass functions. Ongoing efforts
to constrain these quantities for a large fraction of galaxies is thus
valuable in this regard. The assumption of a universal XLF thus cannot hold
regardless of the galaxy environment, and ideally the XLF should be considered
in specific bins of (morphology, specific SFR, $M_{*}$, metallicity), which is
beyond the scope of this study.

Figure 21: Best-fit slope of the deconvolved cumulative XLF of elliptical
galaxies, as a function of the distance limit of the sample.

###  4.5 Nature of HLXs

Taking advantage of the large sky coverage and sensitivities of the catalogues
we used, we were able to constitute a sample of 191 HLX candidates. Unlike
ULXs, their classification by the classifier (mostly as AGN) cannot be
interpreted in terms of contaminants. Moreover, Barrows et al. (2019)
(hereafter B19) found that their HLX population is at least partly consistent
with accreting black holes at the centre of dwarf galaxies, satellites of the
matched host. Similarly, their high luminosities and galactocentric distances
tend to bias the classifier towards the AGN class, particularly given that no
HLX is present in the training sample of the classifier. On the other hand,
because some reliable HLXs are seen to have soft spectra, with most events
detected below 1 keV, some of which being explained by the total or partial
tidal disruption of a star by an intermediate-mass black hole (Godet et al.,
2014; Lin et al., 2018), their classification as soft sources is expected. The
filtering of HLXs thus relies on manual inspection and the removal of known
background and foreground sources.

A few contaminants may be left in the sample of selected HLXs. From the $\log
N-\log S$ estimations reported in Table 9, the fraction of contaminants in the
2SXPS, 4XMM and CSC2 samples are $\sim 20$%, $\sim 70$% and $\sim 60$%, after
removal of spurious sources and before removal of known background and
foreground sources. The shallower sensitivity of 2SXPS results in HLX
candidates in nearby galaxies only, hence the lower contamination rate. The
high contamination rates obtained for other surveys agree with the 70%
estimate of Zolotukhin et al. (2016), but are in contrast with the 7% of
contaminants reported in B19. As they explain, this difference is mainly
driven by the lower angular sizes of their galaxies, most of them being
located beyond 1000 Mpc, a distance range which is not probed by our sample.
By using large and deep photometric redshift surveys, we were able to
drastically reduce our high initial contamination fraction (background sources
represent 85-95% of sources reported in the last column of Table 9). Following
the reasoning of Section 4.2 on ULX candidates, we conservatively expect the
contamination rate among robust candidates (resp. robust and weak candidates)
to be $\sim 20\%$ (resp. $\sim 30$%). Thus, our slightly larger HLX sample is
the result of the larger sky coverage of GLADE, Swift-XRT and XMM-Newton, to
provide local analogues to their HLX population, better suited for future
individual follow-ups.

We investigate the possible nature of our HLX candidates. First of all, the
bimodality of $HR67$ as shown in Figure 18 suggests a difference of nature
between soft and hard HLXs, the former being rarer and sometimes transient,
with at least some of them consistent with partial or tidal disruption events
(e.g. Godet et al. 2014; Lin et al. 2018). No such candidates were present in
the sample of B19, rather showing intrinsically hard spectra or important
absorption. By construction, some AGN in galaxy pairs are also included in our
HLX sample, and flagged as “galaxy pair” candidates. Since such objects are
discarded when their true host is in the GLADE sample, the remaining sources
represent a minority of our HLX candidates, but their luminosity distribution
is consistent with the one of “nuclear” sources (Figure 7). In contrast,
robust HLXs candidates (having no clear optical counterpart and visibly
overlapping the area of the galaxy) and weak HLX candidates (having an optical
counterpart in Gaia, PanSTARRS, or DES, or visibly offset from the area of the
galaxy), taken together or separately, present a steeper luminosity function,
consistent with the bright end slope fitted to ULXs. This maybe supports a
common mechanism behind the X-ray emission of these objects, such as accreting
black holes of stellar to intermediate mass. However, following the method
described in Section 3.2 for ULXs, we find a moderately significant difference
between the radial distribution of ULXs and HLXs. Figure 22 illustrates this
discrepancy: unlike ULXs, HLXs may not follow typical light profiles and may
be equally distributed in the extent of their host. This behaviour further
supports the scenario in which they are the central accretor of a galaxy
satellite or a globular cluster.

Figure 22: Radial distribution of all ULXs (filled circles) and robust HLXs
(empty squares), showing their surface density as a function of galactocentric
distance. Examples of typical galaxy light profiles are shown in solid lines.

Unlike ULXs, a significant fraction ($\sim 35$%) of HLXs is found to have an
optical counterpart. The X-ray to optical flux ratio can often give clues on
the source nature. For instance, in AGN, the optical-UV emission of the disc
scales with the soft X-ray emission, so that their ratio is generally in the
range 0.1–10 (e.g. Maccacaro et al. 1988; Barrows et al. 2019; Tranin et al.
2022). Following the methodology of B19, we first convert the B- and R-band
fluxes given different matched optical surveys (PanSTARRS, DES, Gaia, SDSS and
the Beijing Arizona Sky Survey) to the V-band, using the conversion of Jester
et al. (2005). X-ray broad band fluxes are converted to the 0.3–3.5 keV band
using the same fixed spectrum as above. While this is surely a limitation of
our study, the heterogeneity of our sample and the low signal-to-noise levels
do not allow for a more detailed spectral study. For HLXs present in both B19
and this work, we find a mean bias of our CSC2 fluxes of -0.2 dex with respect
to their rest-frame unabsorbed fitted fluxes ($\pm$0.3 dex deviation), so this
estimation is still reasonable for a population analysis. Finally, for X-ray
sources having no counterpart, we attempt to compute an upper limit on their
optical flux. B19 estimated the local background level directly from optical
images, using an annulus of width being twice the radius of the host. However,
we argue that such a large area is likely to contain other optical sources and
thus overestimate the optical sensitivity at the HLX location, as suggested by
the bright locus of flux upper limits in their $F_{X}-F_{V}$ plane. We adopt a
simpler approach in which this sensitivity is assimilated to the flux of the
faintest optical source having $S/N>5$ within a radius of 15 arcsec around the
HLX candidate. Flux limits obtained with this method are found to scale well
with the limits computed in B19 ($\pm 0.5$ dex deviation), but are on average
ten times fainter. Figure 23 shows the $F_{X}-F_{V}$ plane resulting from this
analysis. In agreement with the finding of B19, most of $F_{X}/F_{V}$ values
in our sample are consistent with expectations for AGN. The use of more
optical surveys allows us to obtain more stringent constraints on the optical
flux upper limit.

Figure 23: $F_{X}-F_{V}$ plane of robust and weak HLX candidates, showing
their observed 0.3–3.5 keV X-ray flux as a function of the V-band flux of
their optical counterpart (filled circles), or the upper limit on $F_{V}$ at
this location (arrows). The grey solid lines indicate ratios of
$F_{X}/F_{V}=0.1$, 1 and 10.

Considering the scenario that HLX candidates are indeed accreting black holes
associated with their host, we estimate their mass and the mass of their
stellar counterpart. Black hole masses are estimated using bolometric
corrections of Marconi et al. (2004) and a fixed Eddington ratio. Unlike B19
who compute bolometric luminosities from their unabsorbed hard (2–10 keV)
X-ray flux obtained after spectral fitting, we rather use the softer (0.5–2
keV) part of the spectrum which better matches the union of several energy
bands in CSC2, 4XMM and 2SXPS, and is less affected by background emission.
Since we do not correct for absorption, this approach is likely to yield some
underestimated masses, but we estimate this bias to be $\lesssim 0.6$ dex from
a comparison of CSC2 fluxes with B19 fluxes. The derived bolometric
corrections are in the range 10–30 depending on the source luminosity (for
comparison, the band 2–10 keV would yield corrections in the range 7–20, and
B19 use the fixed value 10). Assuming spectral properties similar to their
sample, we adopt a value $f_{Edd}=0.24$ to derive black hole masses, and the
standard formula $L_{Edd}\mathrm{[erg/s]}\simeq 1.26\times
10^{38}M_{BH}[\mathrm{M}_{\odot}]$. The resulting masses are in the range
$2\times 10^{3}-2\times 10^{6}\mathrm{M}_{\odot}$ for robust and weak HLXs,
with a median value of $4\times 10^{4}\mathrm{M}_{\odot}$. In particular, 122
of these 169 (robust+weak) candidates fall in the IMBH mass range, below
$10^{5}\mathrm{M}_{\odot}$. Being more luminous, HLXs in galaxy pairs have an
estimated black hole mass $\gtrsim 10^{4}\mathrm{M}_{\odot}$ with a median at
$2\times 10^{5}\mathrm{M}_{\odot}$.

Concerning galaxy stellar masses, as done in B19, we use the mass to light
ratio calibrated on optical colours (Bell et al., 2003). Point-like optical
counterparts of HLXs are assimilated to their stellar associations (we note as
a caveat that this neglects the optical contribution from the putative AGN,
while on the contrary we previously interpreted $F_{V}$ as the AGN flux).
However, since some optical surveys in our sample do not provide $i$-band
magnitudes, the relation giving $M_{*}/L_{r}$ as a function of the $(g-r)$
colour is used instead, before the computation of $M_{*}$. Another maybe
important bias is that the point source flux $F_{r}$ is used in this step,
instead of the flux integrated over the Sérsic component. Except for HLXs in
pairs, for which $M_{*}$ may well be largely underestimated, we do not expect
this effect to have a large impact on the computed stellar masses. Most of
stellar counterparts of robust and weak HLXs are in the dwarf regime,
$M_{*}<3\times 10^{9}\mathrm{M}_{\odot}$.

The resulting $M_{BH}-M_{*}$ plane is shown in Figure 24, also showing the
scaling relation calibrated by Reines & Volonteri (2015) on local AGN in low-
mass galaxies. Upper limits corresponding to the stellar mass of the host are
shown for HLXs having no optical counterpart detected in both $g$ and $r$
bands. While some robust and weak HLXs appear under-massive, most of them are
hard sources, with $F_{0.5-2\leavevmode\nobreak\ \mathrm{keV}}$ accounting for
less than 30% of their broad-band flux: their black hole mass is likely to be
underestimated. A second estimation using the extrapolated band 2–10 keV gives
black hole masses in the range $2\times 10^{3}-2\times
10^{7}\mathrm{M}_{\odot}$ for robust and weak HLXs (median $9\times
10^{4}\mathrm{M}_{\odot}$), with these sources all above $3\times
10^{4}\mathrm{M}_{\odot}$. Some HLXs in pairs are largely under-massive,
possibly because of an underestimated $M_{BH}$ (occurring for instance if
$f_{Edd}<0.24$), an overestimated $M_{*}$ (occurring for instance if the
optical source is dominated by the AGN), or a different nature for these
candidates (such as a tight group of XRBs or ULXs, mimicking a point source).
Interestingly, we find a significant proportion ($\sim 25$%) of overmassive
black holes among robust and weak HLXs, supporting a scenario in which these
black holes have undergone a merger-free growth and/or that their stellar
counterpart have been tidally stripped by the primary host, or is instead a
globular cluster.

This contrasts with B19 results who find a broad alignment of their HLXs to
the relation of Reines & Volonteri (2015), maybe because our work is focusing
on the more local Universe and thus probes lower black hole masses and an