# Statistical study of a large and cleaned sample of ultraluminous and hyperluminous X-ray sources

Hugo Tranin, Natalie Webb, Olivier Godet

► To cite this version:

Hugo Tranin, Natalie Webb, Olivier Godet. Statistical study of a large and cleaned sample of ultra-luminous and hyperluminous X-ray sources. Astronomy and Astrophysics - A&A, 2024, 681, pp.A16. 10.1051/0004-6361/202244952. [hal-04093516](https://hal.science/hal-04093516)

HAL Id: hal-04093516

https://hal.science/hal-04093516

Submitted on 5 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Statistical study of a large and cleaned sample of ultraluminous and hyperluminous X-ray sources*

Hugo Tranin<fnref n="1" />,<fnref n="2" />,✉, Natalie Webb<fnref n="1" />, Olivier Godet<fnref n="1" />, and Erwan Quintin<fnref n="1" />

Received 12 September 2022 / Accepted 12 April 2023

## ABSTRACT

**Context.** Ultraluminous and hyperluminous X-ray (ULX and HLX) sources 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 an unknown nature, but they could be good candidates for accreting intermediate mass black holes (IMBHs). However, a significant portion of the sample of ULX and 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 to better assess their differences and understand their populations.

**Methods.** After a thorough classification of these X-ray sources, we were able to remove ~42% of sources with a signal-to-noise ratio ($S/N$) > 3 which were shown to be contaminants, to obtain the cleanest sample of ULXs and HLXs to date. From a sample of 1342 ULXs and 191 HLXs detected with a $S/N > 3\sigma$, we study the occupation fraction, hardness, variability, radial distribution, and preferred environment of the sources. We built their Malmquist-corrected X-ray luminosity functions (XLFs) and compared them with previous studies. Thanks to the unprecedented size of the sample, we were 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 ~$10^{40}$ erg s$^{-1}$. With our ULX sample, containing ≤2% of contaminants, we are able to confirm that ULXs are located preferentially in spiral galaxies and galaxies with higher star formation rates. 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 as well as lenticular and elliptical galaxies. We note that 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 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.** black hole physics – X-rays: galaxies – X-rays: binaries – methods: observational – accretion, accretion disks

## 1. Introduction

Supermassive black holes (SMBHs, in the mass range $10^6-10^{10} M_{\odot}$) have been found in the centres of the most massive galaxies observed so far and they 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 have yet to be fully understood. Some SMBHs observed at $z > 6$, which is less than 1 Gyr after the Big Bang, are more massive than $10^9 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 M_{\odot}$) mergers and sustained super-Eddington accretion onto these seeds (Volonteri et al. 2008; Haiman et al. 2013; Pacucci & Loeb 2022). However, compelling evidence for the existence of IMBHs 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 insight into the SMBH growth mechanisms. Indeed, for an object spherically accreting hydrogen gas, the luminosity must be lower than the Eddington limit – $L_{\text{Edd}} = 1.26 \times 10^{38} (M_{\odot}/10^{10}) \text{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_{\text{Edd}}$ is one way to estimate a lower limit on the accretor mass. Extra-galactic and off-nuclear ultraluminous X-ray (ULX) sources with X-ray luminosities $L_{\text{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

<footnote n="1">¹ IRAP, Université de Toulouse, CNRS, CNES, 9 Avenue du Colonel Roche, 31028 Toulouse, France</footnote>

<footnote n="2">² Institut d'Estudis Espacials de Catalunya, Universitat de Barcelona (ICC-UB), Martí i Franquès 1, 08028 Barcelona, Spain
e-mail: htranin@icc.ub.edu</footnote>

<footnote>* Tables A.1 and A.2 are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/681/A16</footnote>
# Statistical study of a large and cleaned sample of ultraluminous and hyperluminous X-ray sources*

Hugo Tranin<fnref n="1" />,<fnref n="2" />,✉, Natalie Webb<fnref n="1" />, Olivier Godet<fnref n="1" />, and Erwan Quintin<fnref n="1" />

Received 12 September 2022 / Accepted 12 April 2023

## ABSTRACT

*Context.* Ultraluminous and hyperluminous X-ray (ULX and HLX) sources 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 an unknown nature, but they could be good candidates for accreting intermediate mass black holes (IMBHs). However, a significant portion of the sample of ULX and 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 to better assess their differences and understand their populations.

*Methods.* After a thorough classification of these X-ray sources, we were able to remove ~42% of sources with a signal-to-noise ratio ($S/N$) > 3 which were shown to be contaminants, to obtain the cleanest sample of ULXs and HLXs to date. From a sample of 1342 ULXs and 191 HLXs detected with a $S/N > 3\sigma$, we study the occupation fraction, hardness, variability, radial distribution, and preferred environment of the sources. We built their Malmquist-corrected X-ray luminosity functions (XLFs) and compared them with previous studies. Thanks to the unprecedented size of the sample, we were 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 ~$10^{40}$ erg s$^{-1}$. With our ULX sample, containing $\le 2\%$ of contaminants, we are able to confirm that ULXs are located preferentially in spiral galaxies and galaxies with higher star formation rates. 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 as well as lenticular and elliptical galaxies. We note that 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 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.** black hole physics – X-rays: galaxies – X-rays: binaries – methods: observational – accretion, accretion disks

## 1. Introduction

Supermassive black holes (SMBHs, in the mass range $10^6-10^{10} M_\odot$) have been found in the centres of the most massive galaxies observed so far and they 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 have yet to be fully understood. Some SMBHs observed at $z > 6$, which is less than 1 Gyr after the Big Bang, are more massive than $10^9 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 M_\odot$) mergers and sustained super-Eddington accretion onto these seeds (Volonteri et al. 2008; Haiman et al. 2013; Pacucci & Loeb 2022). However, compelling evidence for the existence of IMBHs 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 insight 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 \times 10^{38} (M/M_\odot) \text{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 (ULX) sources 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

<footnote n="1">¹ IRAP, Université de Toulouse, CNRS, CNES, 9 Avenue du Colonel Roche, 31028 Toulouse, France</footnote>

<footnote n="2">² Institut d'Estudis Espacials de Catalunya, Universitat de Barcelona (ICC-UB), Martí i Franquès 1, 08028 Barcelona, Spain
e-mail: htranin@icc.ub.edu</footnote>

<footnote>* Tables A.1 and A.2 are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/681/A16</footnote>
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 power law is steeper ($\Gamma = 2–4.5$) and cut off at much lower energies ($2–7$ keV compared to $\ge 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 (HLX) sources ($L_X > 10^{41}$ erg s$^{-1}$) may still remain excellent candidates to look for IMBHs, as evidenced with ESO 243-49 HLX-1 (Farrell et al. 2009), which is 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, 2012; Servillat et al. 2011). 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, and they are most often located at a 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 a galaxy catalogue based on Sloan Digital Sky Survey (SDSS) images. 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 with 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, 2011; Liu & Bregman 2005; Liu & Mirabel 2005; Walton et al. 2011; Wang et al. 2016) and recent (Earnshaw et al. 2019; Kovlakas et al. 2020; Inoue et al. 2021; Bernadich 2022; Walton et al. 2022) releases

of large X-ray catalogues. Focussing on the hardness (e.g., Earnshaw et al. 2019) and variability of ULXs (e.g., Sutton et al. 2013; Bernadich 2022), 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 power-law 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 M_{\odot}$ neutron star ($L_{Edd} \sim 2 \times 10^{38}$ erg s$^{-1}$) or a $100 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 and elliptical galaxies are found to be consistent with a dominant population of high-mass and low-mass X-ray binaries (HMXBs and LMXBs), respectively, which is 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 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 select number 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 spectral model fitting, which does not give a physical description of the data (e.g., Swartz et al. 2004), or from 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 (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 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).
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 power law is steeper ($\Gamma = 2–4.5$) and cut off at much lower energies ($2–7$ keV compared to $\ge 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 (HLX) sources ($L_X > 10^{41}$ erg s$^{-1}$) may still remain excellent candidates to look for IMBHs, as evidenced with ESO 243-49 HLX-1 (Farrell et al. 2009), which is 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, 2012; Servillat et al. 2011). 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, and they are most often located at a 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 a galaxy catalogue based on Sloan Digital Sky Survey (SDSS) images. 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 with 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, 2011; Liu & Bregman 2005; Liu & Mirabel 2005; Walton et al. 2011; Wang et al. 2016) and recent (Earnshaw et al. 2019; Kovlakas et al. 2020; Inoue et al. 2021; Bernadich 2022; Walton et al. 2022) releases

of large X-ray catalogues. Focussing on the hardness (e.g., Earnshaw et al. 2019) and variability of ULXs (e.g., Sutton et al. 2013; Bernadich 2022), 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 power-law 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 M_\odot$ neutron star ($L_{Edd} \sim 2 \times 10^{38}$ erg s$^{-1}$) or a $100 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 and elliptical galaxies are found to be consistent with a dominant population of high-mass and low-mass X-ray binaries (HMXBs and LMXBs), respectively, which is 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 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 select number 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 spectral model fitting, which does not give a physical description of the data (e.g., Swartz et al. 2004), or from 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 (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 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 mod-
erately high contamination fraction due to foreground stars
and (mainly) background AGN, often estimated to be ~20%
from the log N–log S diagram of cosmic X-ray back-
ground sources (e.g., Walton et al. 2011), increasing with
X-ray luminosity (e.g., Bernadich 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 has suggested different possible scenar-
ios 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 neu-
tron stars (e.g., Erkut et al. 2019; Kuranov et al. 2020); thermal-
timescale mass transfer, in particular from a Helium-burning sec-
ondary (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, 2017;
Israel et al. 2017; Koliopanos et al. 2017); and accretion of
clumpy stellar wind enhanced by X-ray ionisation (Krtička et al.
2022).

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

In Sect. 2 we explain our ULX selection method, involving
the use of a classification of X-ray sources, and the way we val-
idated the sample and built a complete sub-sample. The sam-
ple of contaminants and the selection of HLX candidates are
also described. In Sect. 3 we present a statistical study of our
clean sample, in terms of the X-ray luminosity function, ULX
rate evolution with the environment, and hardness and variabil-
ity properties. Samples of XRBs, ULXs, and HLXs are com-
pared. We discuss these results in Sect. 4 and compare them to
catalogues and results obtained in previous works. The nature of
HLXs is also investigated. We summarise our study in Sect. 5.
Unless stated otherwise, errors are quoted at the 1-σ level.

## 2. The sample

### 2.1. X-ray catalogues

To obtain the largest possible ULX sample, we consider the
three largest X-ray catalogues to date, generated from obser-
vations 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 ≤ 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⁴, 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 over-
lapping 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 visual inspection, we choose to
keep 1217 sources flagged as extended but not being in a con-
fused 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 detec-
tions with a conf_code lower than 256 (to remove sources with
unreliable parameters due to an overlapping extended source);
214 757 sources remain after this filtering process. As detailed in
Sect. 3.1, 199 sources not included at this stage (because flagged
as extended or confused by the *Chandra* pipeline) were added a
posteriori to the ULX catalogue, after visual inspection.

Likewise, from the *Swift*-XRT catalogue 2SXPS
(Evans et al. 2020) we select sources having detflag = 0
and fieldflag ≤ 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 ~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
130 162 clean sources.

### 2.2. X-ray matches

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 three 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-σ 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 two close *Chandra* sources are confused into a
single *XMM-Newton* or *Swift* one, since these instruments can-
not resolve them). We choose to remove these associations and
flag the corresponding source, unless they become ungrouped
when a crossmatch using the 1-σ 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 412 242 *XMM-Newton*, 214 757 *Chandra* and
130 162 *Swift* unique sources, the cross-correlation retrieves
687 291 unique sources including 65 366 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 asso-
ciated with a galaxy, this fraction of ambiguous associations
increases with galaxy distance and becomes ~10% in galaxies
at 40 Mpc (Fig. 1). Similarly, 0.7% of 4XMM-CSC2 associa-
tions are ambiguous, their median separation is 5.5 arcsec, and
~5% of 4XMM-CSC2 associations are ambiguous in galaxies
Nevertheless, ULX population studies suffer from a mod-
erately high contamination fraction due to foreground stars
and (mainly) background AGN, often estimated to be ~20%
from the log N–log S diagram of cosmic X-ray back-
ground sources (e.g., Walton et al. 2011), increasing with
X-ray luminosity (e.g., Bernadich 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 has suggested different possible scenar-
ios 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 neu-
tron stars (e.g., Erkut et al. 2019; Kuranov et al. 2020); thermal-
timescale mass transfer, in particular from a Helium-burning sec-
ondary (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, 2017;
Israel et al. 2017; Koliopanos et al. 2017); and accretion of
clumpy stellar wind enhanced by X-ray ionisation (Krtička et al.
2022).

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

In Sect. 2 we explain our ULX selection method, involving
the use of a classification of X-ray sources, and the way we val-
idated the sample and built a complete sub-sample. The sam-
ple of contaminants and the selection of HLX candidates are
also described. In Sect. 3 we present a statistical study of our
clean sample, in terms of the X-ray luminosity function, ULX
rate evolution with the environment, and hardness and variabil-
ity properties. Samples of XRBs, ULXs, and HLXs are com-
pared. We discuss these results in Sect. 4 and compare them to
catalogues and results obtained in previous works. The nature of
HLXs is also investigated. We summarise our study in Sect. 5.
Unless stated otherwise, errors are quoted at the 1-σ level.

## 2. The sample

### 2.1. X-ray catalogues

To obtain the largest possible ULX sample, we consider the
three largest X-ray catalogues to date, generated from obser-
vations 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 ≤ 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⁴, so we remove the 561 corresponding
sources. We also remove sources below the detection likelihood
SCDETNL = 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 over-
lapping 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 visual inspection, we choose to
keep 1217 sources flagged as extended but not being in a con-
fused 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 detec-
tions with a conf_code lower than 256 (to remove sources with
unreliable parameters due to an overlapping extended source);
214 757 sources remain after this filtering process. As detailed in
Sect. 3.1, 199 sources not included at this stage (because flagged
as extended or confused by the *Chandra* pipeline) were added a
posteriori to the ULX catalogue, after visual inspection.

Likewise, from the *Swift*-XRT catalogue 2SXPS
(Evans et al. 2020) we select sources having detflag = 0
and fieldflag ≤ 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 ~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
130 162 clean sources.

### 2.2. X-ray matches

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 three 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-σ 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 two close *Chandra* sources are confused into a
single *XMM-Newton* or *Swift* one, since these instruments can-
not resolve them). We choose to remove these associations and
flag the corresponding source, unless they become ungrouped
when a crossmatch using the 1-σ 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 412 242 *XMM-Newton*, 214 757 *Chandra* and
130 162 *Swift* unique sources, the cross-correlation retrieves
687 291 unique sources including 65 366 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 asso-
ciated with a galaxy, this fraction of ambiguous associations
increases with galaxy distance and becomes ~10% in galaxies
at 40 Mpc (Fig. 1). Similarly, 0.7% of 4XMM-CSC2 associa-
tions are ambiguous, their median separation is 5.5 arcsec, and
~5% of 4XMM-CSC2 associations are ambiguous in galaxies
<figure><img src="image_1.png" /><figcaption>Fig. 1. Fraction of 2SXPS-CSC2 and 4XMM-CSC2 ambiguous matches as a function of the host distance given by GLADE.</figcaption></figure>

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

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, 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 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 270 000 entries at $D < 200$ Mpc, ~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 visual 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 ~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 the corresponding 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 ~460 000 galaxies. Such information is necessary to distinguish the study of ULX populations in spiral and elliptical galaxies, which are significantly different as mentioned in Sect. 1. Similar to the definition used by Earnshaw et al. (2019) and Bernadich (2022), we define spiral and elliptical galaxies as those with $t \ge 0$ and $t < 0$, respectively. Irregular galaxies are thus included in the spiral sample in the same way as Walton et al. (2011, 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 ~1.5 million and ~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 SFR from the W3-band (12 µ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 TOP-CAT 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 18 506, 13 055 and 7243 sources for Chandra, XMM-Newton and Swift, respectively.

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

$$d > 3(POSERR + 0.5) \qquad (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
<figure><img src="image_2.png" /><figcaption>Fig. 1. Fraction of 2SXPS-CSC2 and 4XMM-CSC2 ambiguous matches as a function of the host distance given by GLADE.</figcaption></figure>

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

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, 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 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 270 000 entries at $D < 200$ Mpc, ~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 visual 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 ~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 the corresponding 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 ~460 000 galaxies. Such information is necessary to distinguish the study of ULX populations in spiral and elliptical galaxies, which are significantly different as mentioned in Sect. 1. Similar to the definition used by Earnshaw et al. (2019) and Bernadich (2022), we define spiral and elliptical galaxies as those with $t \ge 0$ and $t < 0$, respectively. Irregular galaxies are thus included in the spiral sample in the same way as Walton et al. (2011, 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 ~1.5 million and ~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 SFR from the W3-band (12 µ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 TOP-CAT 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 18 506, 13 055 and 7243 sources for Chandra, XMM-Newton and Swift, respectively.

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

$$d > 3(POSERR + 0.5) \qquad (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⁻¹ for the source to be considered as a ULX candidate. Sources brighter than $10^{41}$ erg s⁻¹ are kept in the same sample at this stage.

## 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 focussed on removing known AGN and stars (e.g., Kovlakas et al. 2020; Bernadich 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 (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 visually 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 Sect. 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 catalogues are detailed in Table 2, as well as the classification results in the training samples. Four known ULXs are contained in the training samples

<table><caption>**Table 1.** Classification parameters for 4XMM, CSC2 and 2SXPS.</caption><thead><tr><th></th><th colspan="4">4XMM</th></tr></thead><tbody><tr><td>Priors P<sub>AGN</sub>, P<sub>Soft</sub>, P<sub>XRB</sub></td><td colspan="4">67%, 8%, 25%</td></tr><tr><td>Weights α<sub>loc</sub>, α<sub>spe</sub>, α<sub>ctp</sub>, α<sub>var</sub></td><td colspan="4">4.8, 9.1, 1.0, 5.5</td></tr><tr><th></th><th colspan="4">CSC2</th></tr><tr><td>Priors P<sub>AGN</sub>, P<sub>Soft</sub>, P<sub>XRB</sub></td><td colspan="4">50%, 10%, 40%</td></tr><tr><td>Weights α<sub>loc</sub>, α<sub>spe</sub>, α<sub>ctp</sub>, α<sub>var</sub></td><td colspan="4">6.9, 5.3, 3.6, 4.2</td></tr><tr><th></th><th colspan="4">2SXPS</th></tr><tr><td>Priors P<sub>AGN</sub>, P<sub>Soft</sub>, P<sub>XRB</sub></td><td colspan="4">67%, 8%, 25%</td></tr><tr><td>Weights α<sub>loc</sub>, α<sub>spe</sub>, α<sub>ctp</sub>, α<sub>var</sub></td><td colspan="4">9.6, 9.1, 1.3, 4.5</td></tr></tbody></table>

**Notes.** For 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 Eq. (3) of Tranin et al. (2022) for further details.

<table><caption>**Table 2.** Confusion matrixes resulting from the classification of 4XMM-DR11, CSC2 and 2SXPS training samples.</caption><thead><tr><th></th><th colspan="4">4XMM</th></tr></thead><tbody><tr><td></td><td>AGN</td><td>Soft</td><td>XRB</td><td>Precision</td></tr><tr><td>→AGN</td><td>1482</td><td>26</td><td>10</td><td>98%</td></tr><tr><td>→Soft</td><td>15</td><td>33</td><td>18</td><td>50%</td></tr><tr><td>→XRB</td><td>29</td><td>20</td><td>320</td><td>87%</td></tr><tr><td>Recall</td><td>97%</td><td>42%</td><td>92%</td><td></td></tr><tr><th></th><th colspan="4">CSC2</th></tr><tr><td></td><td>AGN</td><td>Soft</td><td>XRB</td><td>Precision</td></tr><tr><td>→AGN</td><td>996</td><td>27</td><td>58</td><td>92%</td></tr><tr><td>→Soft</td><td>25</td><td>140</td><td>27</td><td>73%</td></tr><tr><td>→XRB</td><td>81</td><td>80</td><td>830</td><td>84%</td></tr><tr><td>Recall</td><td>90%</td><td>57%</td><td>91%</td><td></td></tr><tr><th></th><th colspan="4">2SXPS</th></tr><tr><td></td><td>AGN</td><td>Soft</td><td>XRB</td><td>Precision</td></tr><tr><td>→AGN</td><td>1568</td><td>31</td><td>14</td><td>97%</td></tr><tr><td>→Soft</td><td>8</td><td>43</td><td>13</td><td>67%</td></tr><tr><td>→XRB</td><td>65</td><td>36</td><td>363</td><td>78%</td></tr><tr><td>Recall</td><td>96%</td><td>39%</td><td>93%</td><td></td></tr></tbody></table>

of XRBs. The classifier is very efficient at retrieving AGN, and retrieves more than 80% of X-ray binaries. False positive rates are also low (≤15%) and can be further decreased by applying a probability threshold on the class under consideration. Figure 2 illustrates the results of the classification in the field of the spiral galaxy NGC 1365.

As a sanity check, we obtain that 99% of sources at $L_X > 10^{42}$ erg s⁻¹ and 94% of sources matching the galaxy centre (see Eq. (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.

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:

$$f_{\text{cont}} = 2 \frac{\pi a b n_{\text{cont}} \text{ sep}^2}{N_{\text{ULX}}} \qquad (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 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⁻¹ for the source to be considered as a ULX candidate. Sources brighter than $10^{41}$ erg s⁻¹ are kept in the same sample at this stage.

## 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 focussed on removing known AGN and stars (e.g., Kovlakas et al. 2020; Bernadich 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 (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 visually 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 Sect. 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 catalogues are detailed in Table 2, as well as the classification results in the training samples. Four known ULXs are contained in the training samples

<table><caption>**Table 1.** Classification parameters for 4XMM, CSC2 and 2SXPS.</caption><thead><tr><th colspan="4">4XMM</th></tr></thead><tbody><tr><td>Priors P<sub>AGN</sub>, P<sub>Soft</sub>, P<sub>XRB</sub></td><td>67%, 8%, 25%</td></tr><tr><td>Weights α<sub>loc</sub>, α<sub>spe</sub>, α<sub>ctp</sub>, α<sub>var</sub></td><td>4.8, 9.1, 1.0, 5.5</td></tr><tr><th colspan="4">CSC2</th></tr><tr><td>Priors P<sub>AGN</sub>, P<sub>Soft</sub>, P<sub>XRB</sub></td><td>50%, 10%, 40%</td></tr><tr><td>Weights α<sub>loc</sub>, α<sub>spe</sub>, α<sub>ctp</sub>, α<sub>var</sub></td><td>6.9, 5.3, 3.6, 4.2</td></tr><tr><th colspan="4">2SXPS</th></tr><tr><td>Priors P<sub>AGN</sub>, P<sub>Soft</sub>, P<sub>XRB</sub></td><td>67%, 8%, 25%</td></tr><tr><td>Weights α<sub>loc</sub>, α<sub>spe</sub>, α<sub>ctp</sub>, α<sub>var</sub></td><td>9.6, 9.1, 1.3, 4.5</td></tr></tbody></table>

**Notes.** For 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 Eq. (3) of Tranin et al. (2022) for further details.

<table><caption>**Table 2.** Confusion matrixes resulting from the classification of 4XMM-DR11, CSC2 and 2SXPS training samples.</caption><thead><tr><th colspan="5">4XMM</th></tr></thead><tbody><tr><td></td><td>AGN</td><td>Soft</td><td>XRB</td><td>Precision</td></tr><tr><td>→AGN</td><td>1482</td><td>26</td><td>10</td><td>98%</td></tr><tr><td>→Soft</td><td>15</td><td>33</td><td>18</td><td>50%</td></tr><tr><td>→XRB</td><td>29</td><td>20</td><td>320</td><td>87%</td></tr><tr><td>Recall</td><td>97%</td><td>42%</td><td>92%</td><td></td></tr><tr><th colspan="5">CSC2</th></tr><tr><td></td><td>AGN</td><td>Soft</td><td>XRB</td><td>Precision</td></tr><tr><td>→AGN</td><td>996</td><td>27</td><td>58</td><td>92%</td></tr><tr><td>→Soft</td><td>25</td><td>140</td><td>27</td><td>73%</td></tr><tr><td>→XRB</td><td>81</td><td>80</td><td>830</td><td>84%</td></tr><tr><td>Recall</td><td>90%</td><td>57%</td><td>91%</td><td></td></tr><tr><th colspan="5">2SXPS</th></tr><tr><td></td><td>AGN</td><td>Soft</td><td>XRB</td><td>Precision</td></tr><tr><td>→AGN</td><td>1568</td><td>31</td><td>14</td><td>97%</td></tr><tr><td>→Soft</td><td>8</td><td>43</td><td>13</td><td>67%</td></tr><tr><td>→XRB</td><td>65</td><td>36</td><td>363</td><td>78%</td></tr><tr><td>Recall</td><td>96%</td><td>39%</td><td>93%</td><td></td></tr></tbody></table>

of XRBs. The classifier is very efficient at retrieving AGN, and retrieves more than 80% of X-ray binaries. False positive rates are also low (≤15%) and can be further decreased by applying a probability threshold on the class under consideration. Figure 2 illustrates the results of the classification in the field of the spiral galaxy NGC 1365.

As a sanity check, we obtain that 99% of sources at $L_X > 10^{42}$ erg s⁻¹ and 94% of sources matching the galaxy centre (see Eq. (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.

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:

$$f_{\text{cont}} = 2 \frac{\pi a b n_{\text{cont}} \text{sep}^2}{N_{\text{ULX}}} \quad (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),
<figure><img src="image_2.png" /></figure>

**Fig. 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 http://xmm-ssc.irap.omp.eu/claxson/xray_analyzer.php?srcquery=53.4019%20-36.1406

<figure><img src="image_3.png" /></figure>

**Fig. 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 µ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.

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 ~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. From visual inspection of 150 sources, we estimate that more than 90% of reliable ULXs are retrieved (as detailed in Sect. 2.9) and that at most 15% of selected candidates are compatible with background contaminants. Removed candidates represent ~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 visually 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 Lite<fnref n="1" /> (Bonnarel et al. 2000; Boch et al. 2014), the Simbad database<fnref n="2" /> and the VizieR catalogue access tool<fnref n="3" /> (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 visually inspect 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% (Sect. 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 (Fig. 3, lower left panel).

## 2.7. HLX sample

Hyperluminous candidates were selected as described in Sect. 2.4, as off-nuclear sources with a mean observed X-ray luminosity $>10^{41}$ erg s⁻<fnref n="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 visually 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.

<footnote n="1">¹ https://aladin.u-strasbg.fr/AladinLite/</footnote>

<footnote n="2">² https://simbad.cds.unistra.fr/simbad/</footnote>

<footnote n="3">³ https://vizier.cds.unistra.fr/viz-bin/VizieR</footnote>
<figure><img src="image_3.png" /></figure>

**Fig. 2.** Illustration of ULX selection on the galaxy NGC 1365. Left: *XMM-Newton* X-ray image (right) DSS Optical image. The *D*₂₅ 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 http://xmm-ssc.irap.omp.eu/claxson/xray_analyzer.php?srcquery=53.4019%20-36.1406

<figure><img src="image_4.png" /></figure>

**Fig. 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*<sub>X</sub>/*F*<sub>W1</sub>): logarithm of the X-ray to infrared (W1-band, 3.4 µm) flux ratio. L<sub>X</sub>: 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.

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<sub>ULX</sub> is the number of ULX candidates
in the considered galaxy. Regardless of this quantity, sources
with a probability to be a contaminant higher than ~95% are also
classified as contaminant: it is notably the case of many fore-
ground sources, and many of the few spurious sources remain-
ing in the sample. Selected candidates are thus those following

$P_{\text{XRB}} > \max(0.05, f_{\text{cont}})$. Here, $P_{\text{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. From visual inspection of 150 sources, we estimate that more than 90% of reliable ULXs are retrieved (as detailed in Sect. 2.9) and that at most 15% of selected candidates are compatible with background contaminants. Removed candidates represent ~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 visually verified a large number of our selected ULX can-
didates, which greatly helped to develop and assess the filter-
ing pipeline described above. To this end, we used the virtual
observatory tools Aladin Lite<fnref n="1" /> (Bonnarel et al. 2000; Boch et al.
2014), the Simbad database<fnref n="2" /> and the VizieR catalogue access
tool<fnref n="3" /> (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 identi-
fied by the classification). We also visually inspect 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%
(Sect. 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 likeli-
hood of the AGN class is enhanced (Fig. 3, lower left panel).

2.7. HLX sample

Hyperluminous candidates were selected as described in Sect. 2.4, as off-nuclear sources with a mean observed X-ray luminosity >10<sup>41</sup> erg s<sup>-1</sup> 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 visually 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.

<footnote n="1">¹ https://aladin.u-strasbg.fr/AladinLite/</footnote>

<footnote n="2">² https://simbad.cds.unistra.fr/simbad/</footnote>

<footnote n="3">³ https://vizier.cds.unistra.fr/viz-bin/VizieR</footnote>
In particular, bright, reliable HLX candidates having a Gaia (Gaia Collaboration 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 billion photometric redshifts. In SDSS, this work was notably performed by Beck et al. (2016), inferring the redshift of 208 million sources up to z ~ 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 ~ 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 2023), using the low-resolution optical Gaia spectra. We use this set of redshift catalogues 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 zerr 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 sub-sample

As introduced for example in Walton et al. (2011), ULX properties and the 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⁻¹ to ensure that all ULXs are well-detected. For instance, Kovlakas et al. (2020; resp. Bernadich 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 Sect. 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 sensi-

<figure><img src="image_4.png" /><figcaption>Fig. 4. Cumulative distribution of GLADE galaxies X-ray $3\sigma$ sensitivities, computed from Eq. (3) as a function of flux. $\theta$ refers to the Chandra off-axis angle.</figcaption></figure>

tivity maps are not accessible for all 3 instruments. For each X-ray catalogue, we infer the flux – effective exposure time relation from subsets of ~$10^5$ detections with signal-to-noise close to 3. Chandra sources show the highest deviation from a single power-law model, because a third parameter, the off-axis angle $\theta$, is determinant in the sensitivity value. We empirically find that once we consider a Chandra effective exposure time $t' = 2t/\max(\theta, 2 \text{arcmin})$, the three X-ray subsets show a $\pm 0.5$ dex deviation from a single power-law model. The resulting $F_{lim}$ relations are as follows:

$$F_{lim} = \begin{cases} \left( \frac{1.2 \times 10^{-15}}{t'} \right)^{0.75} & \text{for CSC2} \\ \left( \frac{1 \times 10^{-13}}{t} \right)^{0.8} & \text{for 4XMM} \\ \left( \frac{5 \times 10^{-13}}{t} \right)^{0.8} & \text{for 2SXPS.} \end{cases} \quad (3)$$

The limiting X-ray sensitivities of Chandra, XMM-Newton and Swift are illustrated in Fig. 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–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.

## 2.9. Sample of contaminants

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(\text{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 visually 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.
In particular, bright, reliable HLX candidates having a Gaia (Gaia Collaboration 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 billion photometric redshifts. In SDSS, this work was notably performed by Beck et al. (2016), inferring the redshift of 208 million sources up to z ~ 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 ~ 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_{\text{norm}}) = 0.03$. Besides, most recently, identification as quasars or galaxies and redshifts were provided for 7.8 million Gaia DR3 sources (Gaia Collaboration 2023), using the low-resolution optical Gaia spectra. We use this set of redshift catalogues 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_{\text{host}} - z| > z_{\text{err}}$, with z and zerr 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_{\text{host}} - z| < z_{\text{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 sub-sample

As introduced for example in Walton et al. (2011), ULX properties and the 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⁻¹ to ensure that all ULXs are well-detected. For instance, Kovlakas et al. (2020; resp. Bernadich 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_{\text{lim}} = 4\pi F_{\text{lim}}D^2$. As further detailed in Sect. 3.1, the contribution of ULXs in each luminosity bin is later weighted by the number of galaxies having $L_{\text{lim}}$ below this luminosity.

The only prerequisite is thus to compute the sensitivity $F_{\text{lim}}$ for each galaxy in the sky coverage of each instrument. We use here a data-driven approach, for simplicity and because sensi-

<figure><img src="image_5.png" /><figcaption>Fig. 4. Cumulative distribution of GLADE galaxies X-ray $3\sigma$ sensitivities, computed from Eq. (3) as a function of flux. $\theta$ refers to the Chandra off-axis angle.</figcaption></figure>

tivity maps are not accessible for all 3 instruments. For each X-ray catalogue, we infer the flux – effective exposure time relation from subsets of ~$10^5$ detections with signal-to-noise close to 3. Chandra sources show the highest deviation from a single power-law model, because a third parameter, the off-axis angle $\theta$, is determinant in the sensitivity value. We empirically find that once we consider a Chandra effective exposure time $t' = 2t/\max(\theta, 2 \text{arcmin})$, the three X-ray subsets show a $\pm 0.5$ dex deviation from a single power-law model. The resulting $F_{\text{lim}}$ relations are as follow:

$$F_{\text{lim}} = \begin{cases} \left( \frac{1.2 \times 10^{-15}}{t'} \right)^{0.75} & \text{for CSC2} \\ \left( \frac{1 \times 10^{-13}}{t} \right)^{0.8} & \text{for 4XMM} \\ \left( \frac{5 \times 10^{-13}}{t} \right)^{0.8} & \text{for 2SXPS.} \end{cases} \quad (3)$$

The limiting X-ray sensitivities of Chandra, XMM-Newton and Swift are illustrated in Fig. 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_{\text{lim}}$ cut keeps 70–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.

## 2.9. Sample of contaminants

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(\text{POSERR} + 0.5)$ and $P_{\text{XRB}} < \max(0.05, f_{\text{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 visually 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.
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Normalised distributions of some properties of selected ULX candidates and candidates classified as contaminants, from the compilation of CSC2, 4XMM-DR11 and 2SXPS.</figcaption></figure>

The distributions of galactocentric distance, mean X-ray luminosities, GLADE distance and Hubble type are shown in Fig. 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 the selection of the HLX candidates, all but six sources above $10^{42} \text{ 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} \text{ 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.

## 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 Sect. 3.1. Table 3 summarises the number counts of ULXs, complete ULXs (i.e., with $L_X > L_{lim}$), HLXs, and galaxies of different morphologies

<figure><img src="image_6.png" /><figcaption>**Fig. 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).</figcaption></figure>

<table><caption>**Table 3.** Samples of ULXs, HLXs and galaxies under study.</caption><thead><tr><th></th><th>CSC2</th><th>4XMM</th><th>2SXPS</th><th>Total (unique)</th></tr></thead><tbody><tr><td>ULX</td><td>1234</td><td>667</td><td>304</td><td>2205 (1901)</td></tr><tr><td>cULX</td><td>925</td><td>470</td><td>266</td><td>1661 (1448)</td></tr><tr><td>HLX</td><td>75</td><td>115</td><td>13</td><td>203 (191)</td></tr><tr><td>cHLX</td><td>63</td><td>92</td><td>13</td><td>168 (158)</td></tr><tr><td>Galaxies</td><td>741</td><td>599</td><td>207</td><td>1547 (1303)</td></tr><tr><td>Spiral</td><td>377</td><td>320</td><td>171</td><td>868 (713)</td></tr><tr><td>Elliptical</td><td>206</td><td>112</td><td>19</td><td>337 (325)</td></tr></tbody></table>

**Notes.** Number counts of ULXs, complete ULXs (i.e., with $L_X > L_{lim}$, Sect. 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.

in each ULX sample. To our knowledge, our work provides both the cleanest (about 2% of contaminants, Sect. 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} \text{ erg s}^{-1}$ for each instrument. All the 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, Sect. 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} \text{ erg s}^{-1}$.

### 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} \text{ 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
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Normalised distributions of some properties of selected ULX candidates and candidates classified as contaminants, from the compilation of CSC2, 4XMM-DR11 and 2SXPS.</figcaption></figure>

The distributions of galactocentric distance, mean X-ray luminosities, GLADE distance and Hubble type are shown in Fig. 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 the selection of the HLX candidates, all but six sources above $10^{42} \text{ 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} \text{ 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.

## 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 Sect. 3.1. Table 3 summarises the number counts of ULXs, complete ULXs (i.e., with $L_X > L_{lim}$), HLXs, and galaxies of different morphologies

<figure><img src="image_7.png" /><figcaption>**Fig. 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).</figcaption></figure>

<table><caption>**Table 3.** Samples of ULXs, HLXs and galaxies under study.</caption><thead><tr><th></th><th>CSC2</th><th>4XMM</th><th>2SXPS</th><th>Total (unique)</th></tr></thead><tbody><tr><td>ULX</td><td>1234</td><td>667</td><td>304</td><td>2205 (1901)</td></tr><tr><td>cULX</td><td>925</td><td>470</td><td>266</td><td>1661 (1448)</td></tr><tr><td>HLX</td><td>75</td><td>115</td><td>13</td><td>203 (191)</td></tr><tr><td>cHLX</td><td>63</td><td>92</td><td>13</td><td>168 (158)</td></tr><tr><td>Galaxies</td><td>741</td><td>599</td><td>207</td><td>1547 (1303)</td></tr><tr><td>Spiral</td><td>377</td><td>320</td><td>171</td><td>868 (713)</td></tr><tr><td>Elliptical</td><td>206</td><td>112</td><td>19</td><td>337 (325)</td></tr></tbody></table>

**Notes.** Number counts of ULXs, complete ULXs (i.e., with $L_X > L_{lim}$, Sect. 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.

in each ULX sample. To our knowledge, our work provides both the cleanest (about 2% of contaminants, Sect. 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} \text{ erg s}^{-1}$ for each instrument. All the 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, Sect. 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} \text{ erg s}^{-1}$.

### 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} \text{ 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 Sect. 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, Sect. 4.3) after making a correction on the CSC2 sample. Indeed, we noted a deficit of *Chandra* sources at medium luminosities ($\ge 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 visual 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 power-law 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 Sect. 4.4.2.

Each XLF is then fitted with two models, a single power law

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

and a broken power law

$$n(> L) = n_{39}(L/10^{39})^{\alpha_1} \text{ if } L < L_b$$

$$n(> L) = n_{39}(L/10^{39})^{\alpha_2} \text{ 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 40 000 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 power law) and $n_{39} \in [0, 1]$,

<figure><img src="image_7.png" /><figcaption>Fig. 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) power-law fits detailed in Table 5. Bottom panel: deconvolved XLF of elliptical galaxies, and residuals from the single power-law fits detailed in Table 4. For comparison, we overlay the cumulative distributions of sources matching their host nucleus, and of ULX candidates classified as background contaminants.</figcaption></figure>

$\alpha_1 \in [0.3, 2]$, $\alpha_2 \in [0, 2]$, $\log(L_b) \in [39, 41]$ (broken power law). Results from the single power law and broken power law fits of the deconvolved cumulative XLF are detailed in Tables 4 and 5. Unlike the XLF of elliptical galaxies, in spiral galaxies, the broken power-law fit is always preferred over a single power-law fit. However, the different catalogues disagree on the exact
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 Sect. 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, Sect. 4.3) after making a correction on the CSC2 sample. Indeed, we noted a deficit of *Chandra* sources at medium luminosities ($\ge 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 visual 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 power-law 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 Sect. 4.4.2.

Each XLF is then fitted with two models, a single power law

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

and a broken power law

$$n(> L) = n_{39}(L/10^{39})^{\alpha_1} \text{ if } L < L_b$$

$$n(> L) = n_{39}(L/10^{39})^{\alpha_2} \text{ 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 40 000 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 power law) and $n_{39} \in [0, 1]$,

<figure><img src="image_8.png" /><figcaption>Fig. 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) power-law fits detailed in Table 5. Bottom panel: deconvolved XLF of elliptical galaxies, and residuals from the single power-law fits detailed in Table 4. For comparison, we overlay the cumulative distributions of sources matching their host nucleus, and of ULX candidates classified as background contaminants.</figcaption></figure>

$\alpha_1 \in [0.3, 2]$, $\alpha_2 \in [0, 2]$, $\log(L_b) \in [39, 41]$ (broken power law). Results from the single power law and broken power law fits of the deconvolved cumulative XLF are detailed in Tables 4 and 5. Unlike the XLF of elliptical galaxies, in spiral galaxies, the broken power-law fit is always preferred over a single power-law fit. However, the different catalogues disagree on the exact
<table><caption>**Table 4.** Parameters of the fits of the deconvolved XLF with a single power-law model.</caption><thead><tr><th>Sample</th><th>n<sub>39</sub></th><th>α</th><th>χ<sup>2</sup><sub>f</sub> (d.o.f.)</th></tr></thead><tbody><tr><td>CSC2, ETG</td><td>0.31<sup>+0.01</sup><sub>-0.02</sub></td><td>0.99<sup>+0.02</sup><sub>-0.02</sub></td><td>0.54 (28)</td></tr><tr><td>4XMM, ETG</td><td>0.16<sup>+0.02</sup><sub>-0.03</sub></td><td>0.99<sup>+0.05</sup><sub>-0.04</sub></td><td>0.09 (27)</td></tr><tr><td>2SXPS, ETG</td><td>0.03<sup>+0.02</sup><sub>-0.01</sub></td><td>1.02<sup>+0.30</sup><sub>-0.24</sub></td><td>0.68 (9)</td></tr><tr><td>Compiled, ETG</td><td>0.29<sup>+0.01</sup><sub>-0.01</sub></td><td>1.06<sup>+0.02</sup><sub>-0.02</sub></td><td>0.67 (28)</td></tr></tbody></table>

**Notes.** ‘ETG’ refer to the sample of complete ULXs in elliptical galaxies. $n_{39}$ and $\alpha$ are the amplitude and slope of the power-law model, respectively. See the text for more details.

location of the break, even at the three sigma level. This discrepancy is further discussed in Sect. 4.4. However, it is to be noted that all three datasets are well-fitted by a power law breaking at $L_{X,\text{break}} = 5 \times 10^{39}$ erg s$^{-1}$ (fixed parameter), with $\chi^2_f < 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 ~80 $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 Sect. 4.4.1.

## 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, SFR, 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. We find that ULX rates globally increase with the galaxy stellar mass, regardless of the SFR. Besides, ULX rates increase with the SFR, at least in spiral galaxies. The rate of ULX is higher in spiral galaxies, regardless of the SFR; and a significant ULX population does exist in elliptical galaxies.

Modelling this ULX rate – SFR relation with a power law, $\log(n_{ULX}) = \alpha \log(\text{SFR}) + \beta$ where $\alpha$ is the slope and $\beta$ the normalisation, we obtain $\alpha = 0.43 \pm 0.04$ for spiral galaxies shown in the middle panel of Fig. 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 SFR < 1 $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 Fig. 8, we can see an excess of ULXs in dwarf galaxies (hardly significant since it disappears with a different binning) 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).

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 (Fig. 9). We note that this trend was reported in Bernadich (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 Sect. 4.4.2.

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 Fig. 10, where the first bin is not shown since the ULX census in this inner part of the galaxy ($\le 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 sub-sample 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.

## 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 \le 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_{\text{XRB}} > 0.7$, $\text{sep} < 0.8$ and $L_{\text{X}} < 10^{39}$ erg s$^{-1}$ make up our sample of XRB candidates. The contamination rate among these selected candidates is 10–12% from visual 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
<table><caption>**Table 4.** Parameters of the fits of the deconvolved XLF with a single power-law model.</caption><thead><tr><th>Sample</th><th>n<sub>39</sub></th><th>α</th><th>χ<sup>2</sup><sub>f</sub> (d.o.f.)</th></tr></thead><tbody><tr><td>CSC2, ETG</td><td>0.31<sup>+0.01</sup><sub>-0.02</sub></td><td>0.99<sup>+0.02</sup><sub>-0.02</sub></td><td>0.54 (28)</td></tr><tr><td>4XMM, ETG</td><td>0.16<sup>+0.02</sup><sub>-0.03</sub></td><td>0.99<sup>+0.05</sup><sub>-0.04</sub></td><td>0.09 (27)</td></tr><tr><td>2SXPS, ETG</td><td>0.03<sup>+0.02</sup><sub>-0.01</sub></td><td>1.02<sup>+0.30</sup><sub>-0.24</sub></td><td>0.68 (9)</td></tr><tr><td>Compiled, ETG</td><td>0.29<sup>+0.01</sup><sub>-0.01</sub></td><td>1.06<sup>+0.02</sup><sub>-0.02</sub></td><td>0.67 (28)</td></tr></tbody></table>

**Notes.** ‘ETG’ refer to the sample of complete ULXs in elliptical galaxies. $n_{39}$ and $\alpha$ are the amplitude and slope of the power-law model, respectively. See the text for more details.

location of the break, even at the three sigma level. This discrepancy is further discussed in Sect. 4.4. However, it is to be noted that all three datasets are well-fitted by a power law breaking at $L_{X,\text{break}} = 5 \times 10^{39} \text{ erg s}^{-1}$ (fixed parameter), with $\chi^2_f < 1$.

Unlike the lower luminosity break seen in the XLF of elliptical galaxies at a few $10^{38} \text{ 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 ~80 $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 Sect. 4.4.1.

## 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, SFR, 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. We find that ULX rates globally increase with the galaxy stellar mass, regardless of the SFR. Besides, ULX rates increase with the SFR, at least in spiral galaxies. The rate of ULX is higher in spiral galaxies, regardless of the SFR; and a significant ULX population does exist in elliptical galaxies.

Modelling this ULX rate – SFR relation with a power law, $\log(n_{ULX}) = \alpha \log(\text{SFR}) + \beta$ where $\alpha$ is the slope and $\beta$ the normalisation, we obtain $\alpha = 0.43 \pm 0.04$ for spiral galaxies shown in the middle panel of Fig. 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 SFR < 1 $M_{\odot} \text{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 Fig. 8, we can see an excess of ULXs in dwarf galaxies (hardly significant since it disappears with a different binning) 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_{\ast}$ thus decreases with stellar mass, in agreement with Walton et al. (2011) and Kovlakas et al. (2020).

We find tentative evidence that ULXs and HXLs do not share exactly the same environment: the latter tend to be hosted in equal rates in spiral and elliptical galaxies (Fig. 9). We note that this trend was reported in Bernadich (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 Sect. 4.4.2.

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 \text{ arcsec}$, as a function of their galactocentric distance. The result is shown in Fig. 10, where the first bin is not shown since the ULX census in this inner part of the galaxy ($\le 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 sub-sample 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.

## 3.3. Comparison of XRB and ULX

It is now widely accepted that most ULXs below $\sim 10^{41} \text{ 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 \le 20 \text{ 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_{\text{XRB}} > 0.7$, $\text{sep} < 0.8$ and $L_{\text{X}} < 10^{39} \text{ erg s}^{-1}$ make up our sample of XRB candidates. The contamination rate among these selected candidates is 10–12% from visual 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
<table><caption>**Table 5.** Parameters of the fits of the deconvolved XLF of spiral galaxies with a broken power-law model.</caption><thead><tr><th>Sample</th><th>n<sub>39</sub></th><th>α<sub>1</sub></th><th>α<sub>2</sub></th><th>log(L<sub>break</sub>)</th><th>χ<sup>2</sup><sub>r</sub> (d.o.f.)</th></tr></thead><tbody><tr><td>CSC2, LTG</td><td>0.45<sup>+0.02</sup><sub>-0.02</sub></td><td>0.67<sup>+0.09</sup><sub>-0.10</sub></td><td>1.36<sup>+0.06</sup><sub>-0.05</sub></td><td>39.45<sup>+0.07</sup><sub>-0.06</sub></td><td>0.37 (27)</td></tr><tr><td>4XMM, LTG</td><td>0.35<sup>+0.02</sup><sub>-0.02</sub></td><td>0.81<sup>+0.06</sup><sub>-0.07</sub></td><td>1.24<sup>+0.05</sup><sub>-0.04</sub></td><td>39.73<sup>+0.11</sup><sub>-0.11</sub></td><td>0.19 (33)</td></tr><tr><td>2SXPS, LTG</td><td>0.29<sup>+0.02</sup><sub>-0.02</sub></td><td>0.80<sup>+0.04</sup><sub>-0.05</sub></td><td>1.73<sup>+0.12</sup><sub>-0.10</sub></td><td>39.97<sup>+0.07</sup><sub>-0.08</sub></td><td>0.33 (24)</td></tr><tr><td>Compiled, LTG</td><td>0.41<sup>+0.01</sup><sub>-0.01</sub></td><td>0.86<sup>+0.04</sup><sub>-0.06</sub></td><td>1.31<sup>+0.03</sup><sub>-0.03</sub></td><td>39.66<sup>+0.08</sup><sub>-0.09</sub></td><td>0.77 (33)</td></tr><tr><td>CSC2, LTG39</td><td>240<sup>+12</sup><sub>-12</sub></td><td>0.80<sup>+0.08</sup><sub>-0.10</sub></td><td>1.90<sup>+0.07</sup><sub>-0.14</sub></td><td>39.67<sup>+0.08</sup><sub>-0.11</sub></td><td>1.13 (12)</td></tr><tr><td>4XMM, LTG39</td><td>168<sup>+10</sup><sub>-10</sub></td><td>0.71<sup>+0.08</sup><sub>-0.09</sub></td><td>1.80<sup>+0.14</sup><sub>-0.17</sub></td><td>39.74<sup>+0.10</sup><sub>-0.10</sub></td><td>0.17 (14)</td></tr><tr><td>2SXPS, LTG39</td><td>150<sup>+10</sup><sub>-10</sub></td><td>0.65<sup>+0.10</sup><sub>-0.12</sub></td><td>1.39<sup>+0.25</sup><sub>-0.17</sub></td><td>39.70<sup>+0.20</sup><sub>-0.19</sub></td><td>0.13 (18)</td></tr><tr><td>Compiled, LTG39</td><td>424<sup>+15</sup><sub>-15</sub></td><td>0.89<sup>+0.05</sup><sub>-0.06</sub></td><td>1.97<sup>+0.03</sup><sub>-0.06</sub></td><td>39.71<sup>+0.05</sup><sub>-0.06</sub></td><td>0.60 (12)</td></tr></tbody></table>

**Notes.** ‘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 power-law model, respectively. See the text for more details.

<figure><img src="image_8.png" /><figcaption>**Fig. 8.** Rates of ULX as a function of different properties characterising galaxy environment.</figcaption></figure>

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

<figure><img src="image_9.png" /><figcaption>**Fig. 9.** Evolution with X-ray luminosity of the ratio between ULX rates $N(L_X)$ in spiral and elliptical galaxies.</figcaption></figure>

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. To combine the detections coming from all three instruments, we built a flux hardness ratio
<table><caption>**Table 5.** Parameters of the fits of the deconvolved XLF of spiral galaxies with a broken power-law model.</caption><thead><tr><th>Sample</th><th>n<sub>39</sub></th><th>α<sub>1</sub></th><th>α<sub>2</sub></th><th>log(L<sub>break</sub>)</th><th>χ<sup>2</sup><sub>r</sub> (d.o.f.)</th></tr></thead><tbody><tr><td>CSC2, LTG</td><td>0.45<sup>+0.02</sup><sub>-0.02</sub></td><td>0.67<sup>+0.09</sup><sub>-0.10</sub></td><td>1.36<sup>+0.06</sup><sub>-0.05</sub></td><td>39.45<sup>+0.07</sup><sub>-0.06</sub></td><td>0.37 (27)</td></tr><tr><td>4XMM, LTG</td><td>0.35<sup>+0.02</sup><sub>-0.02</sub></td><td>0.81<sup>+0.06</sup><sub>-0.07</sub></td><td>1.24<sup>+0.05</sup><sub>-0.04</sub></td><td>39.73<sup>+0.11</sup><sub>-0.11</sub></td><td>0.19 (33)</td></tr><tr><td>2SXPS, LTG</td><td>0.29<sup>+0.02</sup><sub>-0.02</sub></td><td>0.80<sup>+0.04</sup><sub>-0.05</sub></td><td>1.73<sup>+0.12</sup><sub>-0.10</sub></td><td>39.97<sup>+0.07</sup><sub>-0.08</sub></td><td>0.33 (24)</td></tr><tr><td>Compiled, LTG</td><td>0.41<sup>+0.01</sup><sub>-0.01</sub></td><td>0.86<sup>+0.04</sup><sub>-0.06</sub></td><td>1.31<sup>+0.03</sup><sub>-0.03</sub></td><td>39.66<sup>+0.08</sup><sub>-0.09</sub></td><td>0.77 (33)</td></tr><tr><td>CSC2, LTG39</td><td>240<sup>+12</sup><sub>-12</sub></td><td>0.80<sup>+0.08</sup><sub>-0.10</sub></td><td>1.90<sup>+0.07</sup><sub>-0.14</sub></td><td>39.67<sup>+0.08</sup><sub>-0.11</sub></td><td>1.13 (12)</td></tr><tr><td>4XMM, LTG39</td><td>168<sup>+10</sup><sub>-10</sub></td><td>0.71<sup>+0.08</sup><sub>-0.09</sub></td><td>1.80<sup>+0.14</sup><sub>-0.17</sub></td><td>39.74<sup>+0.10</sup><sub>-0.10</sub></td><td>0.17 (14)</td></tr><tr><td>2SXPS, LTG39</td><td>150<sup>+10</sup><sub>-10</sub></td><td>0.65<sup>+0.10</sup><sub>-0.12</sub></td><td>1.39<sup>+0.25</sup><sub>-0.17</sub></td><td>39.70<sup>+0.20</sup><sub>-0.19</sub></td><td>0.13 (18)</td></tr><tr><td>Compiled, LTG39</td><td>424<sup>+15</sup><sub>-15</sub></td><td>0.89<sup>+0.05</sup><sub>-0.06</sub></td><td>1.97<sup>+0.03</sup><sub>-0.06</sub></td><td>39.71<sup>+0.05</sup><sub>-0.06</sub></td><td>0.60 (12)</td></tr></tbody></table>

**Notes.** 'LTG' refer to the samples of complete ULXs in spiral galaxies. 'LTG39' stands for the samples of spiral galaxies having $L_{\text{lim}} < 10^{39}$ erg s$^{-1}$; they correspond to the classical, non-renormalised XLF. $n_{39}$, $\alpha_1$, $\alpha_2$ and $\log(L_{\text{break}})$ are the amplitude, faint-end slope, bright end slope and break luminosity of the broken power-law model, respectively. See the text for more details.

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Rates of ULX as a function of different properties characterising galaxy environment.</figcaption></figure>

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_{\text{AGN}} > P_{\text{Soft}}$) and

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Evolution with X-ray luminosity of the ratio between ULX rates $N(L_X)$ in spiral and elliptical galaxies.</figcaption></figure>

soft ($P_{\text{AGN}} < P_{\text{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. To combine the detections coming from all three instruments, we built a flux hardness ratio
<figure><img src="image_10.png" /><figcaption>Fig. 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.</figcaption></figure>

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}}. \qquad (4)$$

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 power-law model with (Γ = 1.7, nH = 3 × 10²⁰ cm⁻²) 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 ten 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³⁹ erg s⁻¹ are considered as XRB (resp. ULX). Figure 12 shows the detected luminosities of these 59 sources, with hardness being colour-coded.

While a substantial fraction (~50%) of sources do not present a significant hardness evolution (ΔHR < 0.2), most

<figure><img src="image_11.png" /><figcaption>Fig. 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: between 0.5–1 and 1–2 keV (SC_HR2), between 1–2 and 2–4.5 keV (SC_HR3), between 0.5–1.2 and 1.2–2 keV (hard_ms), between 1.2–2 and 2–7 keV (hard_hm), between 0.3–1 and 1–2 keV (HR1), and between 1–2 and 2–10 keV (HR2).</figcaption></figure>

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
<figure><img src="image_11.png" /><figcaption>Fig. 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.</figcaption></figure>

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}}. \qquad (4)$$

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 power-law model with (Γ = 1.7, nH = 3 × 10²⁰ cm⁻²) 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 ten 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³⁹ erg s⁻¹ are considered as XRB (resp. ULX). Figure 12 shows the detected luminosities of these 59 sources, with hardness being colour-coded.

While a substantial fraction (~50%) of sources do not present a significant hardness evolution (ΔHR < 0.2), most

<figure><img src="image_12.png" /><figcaption>Fig. 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: between 0.5–1 and 1–2 keV (SC_HR2), between 1–2 and 2–4.5 keV (SC_HR3), between 0.5–1.2 and 1.2–2 keV (hard_ms), between 1.2–2 and 2–7 keV (hard_hm), between 0.3–1 and 1–2 keV (HR1), and between 1–2 and 2–10 keV (HR2).</figcaption></figure>

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
<figure><img src="image_12.png" /><figcaption>**Fig. 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.</figcaption></figure>

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.

### 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 three 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 min) than the typical duration of our multi-mission observations (~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 (Fig. 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 min long bins, and recording the peak rate of each band. The result is shown in the bottom panel of Fig. 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 four galaxies observed for more than 200 ks: M 51, M 81, M 101 and NGC 300. After removal of these galaxies, the ratio between the
<figure><img src="image_13.png" /><figcaption>**Fig. 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.</figcaption></figure>

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.

### 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 three 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 min) than the typical duration of our multi-mission observations (~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 (Fig. 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 min long bins, and recording the peak rate of each band. The result is shown in the bottom panel of Fig. 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 four galaxies observed for more than 200 ks: M 51, M 81, M 101 and NGC 300. After removal of these galaxies, the ratio between the
<figure><img src="image_13.png" /><figcaption>**Fig. 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.</figcaption></figure>

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.

### 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 four 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 seven 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 sub-sample of initial and selected HLX candidates in Table 6.

The distributions of the distances and mean observed 0.5–10 keV X-ray luminosities of the 191 HLX candidates are shown in Fig. 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 visual inspection as lying visibly outside the extent of the host. M 82 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 ~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).

### 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 disc 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 (Fig. 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 (Fig. 16): the galaxy stellar mass, its SFR 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 \ge 3$ (equivalent to 1.5 major graduations in Fig. 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 SFR.
<figure><img src="image_14.png" /><figcaption>**Fig. 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.</figcaption></figure>

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.

## 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 four 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 seven 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 sub-sample of initial and selected HLX candidates in Table 6.

The distributions of the distances and mean observed 0.5–10 keV X-ray luminosities of the 191 HLX candidates are shown in Fig. 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 visual inspection as lying visibly outside the extent of the host. M 82 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 ~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).

## 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 disc 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 (Fig. 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 (Fig. 16): the galaxy stellar mass, its SFR 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 \ge 3$ (equivalent to 1.5 major graduations in Fig. 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 SFR.
<figure><img src="image_14.png" /><figcaption>**Fig. 14.** Images of a random sample of new HLX candidates. The images are structured as follows: Cols. (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σ X-ray position error circle is shown in cyan.</figcaption></figure>

### 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
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Images of a random sample of new HLX candidates. The images are structured as follows: Cols. (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σ X-ray position error circle is shown in cyan.</figcaption></figure>

### 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
<table><caption>**Table 6.** Sub-samples of HLXs under study.</caption><thead><tr><th></th><th>CSC2</th><th>4XMM</th><th>2SXPS</th><th>Total (unique)</th></tr></thead><tbody><tr><td>Initial candidates</td><td>195</td><td>360</td><td>110</td><td>665 (619)</td></tr><tr><td>With redshift</td><td>96</td><td>187</td><td>41</td><td>324 (298)</td></tr><tr><td>Background</td><td>82</td><td>165</td><td>35</td><td>282 (259)</td></tr><tr><td>Foreground</td><td>2</td><td>3</td><td>0</td><td>5 (5)</td></tr><tr><td>Distance match</td><td>12</td><td>19</td><td>6</td><td>37 (34)</td></tr><tr><td>Selected candidates</td><td>75</td><td>115</td><td>13</td><td>203 (191)</td></tr><tr><td>Robust</td><td>51</td><td>66</td><td>5</td><td>122 (115)</td></tr><tr><td>Weak</td><td>15</td><td>35</td><td>7</td><td>57 (54)</td></tr><tr><td>Galaxy pair</td><td>9</td><td>14</td><td>1</td><td>24 (22)</td></tr><tr><td>With opt. counterpart</td><td>29</td><td>63</td><td>8</td><td>100 (94)</td></tr><tr><td>Soft (HR67 &lt; -0.6)</td><td>2</td><td>5</td><td>1</td><td>8 (6)</td></tr><tr><td>Hard (HR67 &gt; -0.6)</td><td>67</td><td>110</td><td>12</td><td>189 (179)</td></tr></tbody></table>

**Notes.** Number counts of initial and selected HLXs in different sub-samples (see Sects. 2.7 and 3.4, respectively, for details on the sub-samples). The sample of 'soft HLXs' is defined in Sect. 3.5.1, and is the complementary of 'hard HLXs' (note that six Chandra candidates are not included in these samples, because their HR67 could not be calculated).

<figure><img src="image_15.png" /><figcaption>**Fig. 15.** Distribution of some properties of selected samples of ULXs and HLXs. Left: distance of the host. ULXs are detected up to 500 Mpc and HLXs up to 1100 Mpc. Right: 0.5–10 keV observed mean X-ray luminosity of the 191 HLXs candidates, including the 22 galaxy pair candidates.</figcaption></figure>

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 power-law spectrum, and in this case such values would only be produced by $\Gamma \le 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 (Soria et al., in prep.), and 4XMM J085253.8+180110, a weak candidate located in a cluster of galaxies. Another soft outlier (not visible in Fig. 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)

<figure><img src="image_16.png" /><figcaption>**Fig. 16.** Rates of ULX and HLX as a function of different properties characterising galaxy environment.</figcaption></figure>

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

### 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
<table><caption>**Table 6.** Sub-samples of HLXs under study.</caption><thead><tr><th></th><th>CSC2</th><th>4XMM</th><th>2SXPS</th><th>Total (unique)</th></tr></thead><tbody><tr><td>Initial candidates</td><td>195</td><td>360</td><td>110</td><td>665 (619)</td></tr><tr><td>With redshift</td><td>96</td><td>187</td><td>41</td><td>324 (298)</td></tr><tr><td>Background</td><td>82</td><td>165</td><td>35</td><td>282 (259)</td></tr><tr><td>Foreground</td><td>2</td><td>3</td><td>0</td><td>5 (5)</td></tr><tr><td>Distance match</td><td>12</td><td>19</td><td>6</td><td>37 (34)</td></tr><tr><td>Selected candidates</td><td>75</td><td>115</td><td>13</td><td>203 (191)</td></tr><tr><td>Robust</td><td>51</td><td>66</td><td>5</td><td>122 (115)</td></tr><tr><td>Weak</td><td>15</td><td>35</td><td>7</td><td>57 (54)</td></tr><tr><td>Galaxy pair</td><td>9</td><td>14</td><td>1</td><td>24 (22)</td></tr><tr><td>With opt. counterpart</td><td>29</td><td>63</td><td>8</td><td>100 (94)</td></tr><tr><td>Soft (HR67 &lt; -0.6)</td><td>2</td><td>5</td><td>1</td><td>8 (6)</td></tr><tr><td>Hard (HR67 &gt; -0.6)</td><td>67</td><td>110</td><td>12</td><td>189 (179)</td></tr></tbody></table>

**Notes.** Number counts of initial and selected HLXs in different sub-samples (see Sects. 2.7 and 3.4, respectively, for details on the sub-samples). The sample of 'soft HLXs' is defined in Sect. 3.5.1, and is the complementary of 'hard HLXs' (note that six Chandra candidates are not included in these samples, because their HR67 could not be calculated).

<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Distribution of some properties of selected samples of ULXs and HLXs. Left: distance of the host. ULXs are detected up to 500 Mpc and HLXs up to 1100 Mpc. Right: 0.5–10 keV observed mean X-ray luminosity of the 191 HLXs candidates, including the 22 galaxy pair candidates.</figcaption></figure>

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 power-law 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 (Soria et al., in prep.), and 4XMM J085253.8+180110, a weak candidate located in a cluster of galaxies. Another soft outlier (not visible in Fig. 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)

<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Rates of ULX and HLX as a function of different properties characterising galaxy environment.</figcaption></figure>

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

### 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
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. 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.</figcaption></figure>

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.

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

<figure><img src="image_20.png" /><figcaption>**Fig. 18.** Distribution of the hardness ratio HR67 between the 0.2–2 and 2–12 keV energy bands for selected HLX candidates.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 19.** Distribution of the variability ratio between observations ($\log(\text{Fratio}) = \log(F_{\text{max}}/F_{\text{min}})$), for ULXs and HLXs.</figcaption></figure>

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 visual 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 focussing on one catalogue decreases the survey coverage and misses some candidates as they may have been fainter than $10^{39} \text{ 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} \text{ erg s}^{-1}$, the evolution of the ULX rate with SFR and the decrease of specific ULX rate with galaxy stellar mass.

Bernadich (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 visual inspection of the sources
<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 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.</figcaption></figure>

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.

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

<figure><img src="image_21.png" /><figcaption>**Fig. 18.** Distribution of the hardness ratio HR67 between the 0.2–2 and 2–12 keV energy bands for selected HLX candidates.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. 19.** Distribution of the variability ratio between observations ($\log(\text{Fratio}) = \log(F_{\text{max}}/F_{\text{min}})$), for ULXs and HLXs.</figcaption></figure>

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 visual 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 focussing on one catalogue decreases the survey coverage and misses some candidates as they may have been fainter than $10^{39} \text{ 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} \text{ erg s}^{-1}$, the evolution of the ULX rate with SFR and the decrease of specific ULX rate with galaxy stellar mass.

Bernadich (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 visual 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 visual 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,\text{err}} > 10^{39} \text{ erg s}^{-1}$ in at least one detection, and only 632 of their candidates have $\langle L_X \rangle > 10^{39} \text{ erg s}^{-1}$.

Their work also focusses on the statistical study of bright ULXs, having $L_X + L_{X,\text{err}} > 5 \times 10^{40} \text{ erg cm}^{-2} \text{ 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} \text{ erg s}^{-1}$ (only 25 of their bright candidates have $\langle L_X \rangle > 10^{41} \text{ 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 \text{ 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 (Sect. 4.5). We do not find evidence in our sample for the softer average spectra for HLXs reported in Bernadich (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 (Fig. 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} \text{ 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 ~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). We obtain 27% of these sources discarded as contaminants, which is higher than their estimated background contamination rate (18–23%). Most of the ULXs of Bernadich (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 (~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} \text{ erg s}^{-1}$. Few of their galaxies are present in the GLADE catalogue: as a result, only ten of their candidates are present in our CSC2 sample. Two of them have $\langle L_X \rangle < 10^{41} \text{ 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 six candidates out of the possible ten 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 Sect. 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 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 power-law spectrum with $\Gamma = 1.69$ and $n_H \sim 3 \times 10^{20} \text{ 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 s}^{-1}}{4\pi D^2}$.

Because our centre exclusion radius $3(\text{POSERR}+0.5)$ depends on the source position error, we estimate the number of
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 visual 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,\text{err}} > 10^{39} \text{ erg s}^{-1}$ in at least one detection, and only 632 of their candidates have $\langle L_X \rangle > 10^{39} \text{ erg s}^{-1}$.

Their work also focusses on the statistical study of bright ULXs, having $L_X + L_{X,\text{err}} > 5 \times 10^{40} \text{ erg cm}^{-2} \text{ 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} \text{ erg s}^{-1}$ (only 25 of their bright candidates have $\langle L_X \rangle > 10^{41} \text{ 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 \text{ 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 (Sect. 4.5). We do not find evidence in our sample for the softer average spectra for HLXs reported in Bernadich (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 (Fig. 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} \text{ 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 ~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). We obtain 27% of these sources discarded as contaminants, which is higher than their estimated background contamination rate (18–23%). Most of the ULXs of Bernadich (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 (~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} \text{ erg s}^{-1}$. Few of their galaxies are present in the GLADE catalogue: as a result, only ten of their candidates are present in our CSC2 sample. Two of them have $\langle L_X \rangle < 10^{41} \text{ 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 six candidates out of the possible ten 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 Sect. 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 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 power-law spectrum with $\Gamma = 1.69$ and $n_H \sim 3 \times 10^{20} \text{ 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 s}^{-1}}{4\pi D^2}$.

Because our centre exclusion radius $3(\text{POSERR}+0.5)$ depends on the source position error, we estimate the number of