# Ultraluminous X-ray sources

Andrew King, Jean-Pierre Lasota, Matthew Middleton

► To cite this version:

Andrew King, Jean-Pierre Lasota, Matthew Middleton. Ultraluminous X-ray sources. New As-
tron.Rev., 2023, 96, pp.101672. 10.1016/j.newar.2022.101672. hal-03958508

HAL Id: hal-03958508

https://hal.science/hal-03958508

Submitted on 27 Jan 2023

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

# Ultraluminous X-ray sources

Andrew King<sup>a,b,c</sup>, Jean-Pierre Lasota<sup>d,e</sup>, Matthew Middleton<sup>f,*</sup>

## ARTICLE INFO

**Keywords:**

Accretion

ULXs

Black holes

Neutron stars

## ABSTRACT

The study of ultraluminous X-ray sources (ULXs) has changed dramatically over the last decade. In this review we first describe the most important observations of ULXs in various wavebands, and across multiple scales in space and time. We discuss recent progress and current unanswered questions. We consider the range of current theories of ULX properties in the light of this observational progress. Applying these models to neutron-star ULXs offers particularly stringent tests, as this is the unique case where the mass of the accretor is effectively fixed.

## 1. Introduction

Ultraluminous X-ray sources (ULXs) are usually defined by fluxes which, if assumed isotropic, give luminosities

$$L_X > 10^{39} \text{ erg s}^{-1}. \qquad (1)$$

The implication is that ULXs appear to be above the Eddington critical luminosity

$$L_{\text{Edd}} = 1.3 \times 10^{38} \left( \frac{M}{M_{\odot}} \right) \text{erg s}^{-1}, \qquad (2)$$

for accretor masses $M \lesssim 10 M_{\odot}$ (see Eq. (20) for the exact definition of the luminosity $L_{\text{Edd}}$). Early systematic X-ray surveys found significant numbers of these objects: Fabbiano (1989) reported 16 sources with $L_X > 10^{39} \text{ erg s}^{-1}$. By now, more than about 1800 ULXs (Walton et al., 2022) are known, including several with $L_X > 10^{41}-10^{42} \text{ erg s}^{-1}$. ULXs are (almost) all in external galaxies, but not located in their nuclei, making them distinct from accreting supermassive black holes

From the beginning it was clear that the brightest ULXs required a significant shift in the standard paradigm for the then known accretion-powered sources, and this provoked a wide variety of speculations. An obvious possibility was that the accreting mass $M$ was unusually

large. Colbert and Mushotzky (1999) used a sample of ~20 ULXs with apparent X-ray luminosities in the range $10^{40}-10^{42} \text{ erg s}^{-1}$ to argue that the typical accretors were black holes with masses in the range $10^2 M_{\odot} \lesssim M \lesssim 10^4 M_{\odot}$. The suggested objects soon acquired the name ‘intermediate-mass black holes’, or IMBHs (Taniguchi et al., 2000). Models for IMBH$^1$ formation include rapid merging of stellar-mass black holes in globular clusters (Miller and Hamilton, 2002), or of stars in dense clusters (Portegies Zwart and McMillan, 2002). A problem for all these models was to ensure that the IMBH would find a companion star and then form a binary able to sustain mass transfer at the required rate (up to ~$10^{-5}-10^{-4} M_{\odot} \text{ yr}^{-1}$).

A quite different idea was that ULXs were fairly standard X-ray binaries with $M \lesssim 10 M_{\odot}$ in unusual states. Begelman (2002) suggested that the accretion discs in standard black hole X-ray binaries might be unstable to the photon-bubble instability, and simply radiate at super-Eddington luminosities. This is limited to luminosities $L_{\text{sph}} \lesssim 3 \times 10^{40} \text{ erg s}^{-1}$, and requires that mass loss from the disc surface is small. Körding et al. (2002) suggested that the high luminosities of ULXs could result from relativistic beaming in jets from microquasars oriented towards the observer. This is again limited, to $L_{\text{sph}} \lesssim 10^{40} \text{ erg s}^{-1}$. Similarly, it was known that magnetic fields with the strength ~$10^{12}-10^{13} \text{ G}$ found in accreting neutron stars reduced the electron scattering opacity

<footnote><sup>a</sup> Astrophysics Division, School of Physics & Astronomy, University of Leicester, Leicester LE1 7RH, UK</footnote>

<footnote><sup>b</sup> Astronomical Institute Anton Pannekoek, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, Netherlands</footnote>

<footnote><sup>c</sup> Leiden Observatory, Leiden University, Niels Bohrweg 2, NL-2333 CA Leiden, Netherlands</footnote>

<footnote><sup>d</sup> Institut d'Astrophysique de Paris, CNRS et Sorbonne Université, UMR 7095, 98bis Bd Arago, 75014 Paris, France</footnote>

<footnote><sup>e</sup> Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, ul. Bartycka 18, 00-716 Warsaw, Poland</footnote>

<footnote><sup>f</sup> Department of Physics and Astronomy, University of Southampton, Highfield, Southampton SO17 1BJ, UK</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: m.j.middleton@soton.ac.uk (M. Middleton).</footnote>

<footnote n="1">¹ Some dwarf galaxies have active nuclei powered by accretion on to black holes of masses ~$10^5 M_{\odot}$, with similar luminosities to some ULXs. These black holes appear to obey the same $M - \sigma$ relation between black hole mass and host spheroid velocity dispersion as universally holds for higher mass nuclear black holes and their host galaxies (see Baldassare et al. 2020, King and Nealon 2019 for a discussion). Unlike these objects, ULXs are not unique systems at the dynamical centres of their hosts. Accordingly, in this review we confine use of the term ‘IMBH’ to objects not in the nuclei of dwarf galaxies (see Section 2.7).</footnote>
for X-rays propagating perpendicular to the field (Canuto et al., 1971),<fnref n="2" /> but this works only for $L_{sph} \lesssim 3 \times 10^{40}$ erg s⁻¹, and even then requires a magnetar-strength field $B \sim 10^{14}$ G (Mushtukov et al., 2015b), although these do not seem to appear in binary systems (see Section 3.9). The limits for all three of these stellar-mass models are incompatible with the luminosities $L_X > 10^{41} - 10^{42}$ erg s⁻¹ observed in some ULXs. Narayan et al. (2017) suggested that radiatively inefficient accretion on to black holes could explain ULXs as the radiative output was tightly beamed. Here much of the accreting gas is swallowed by the hole without radiating.

King et al. (2001) noted that the non-magnetic neutron star in the low-mass X-ray binary Cyg X-2 (not a ULX) was known from evolutionary arguments (King and Ritter, 1999, see Section 3.11 below) to have survived a previous phase where it was subject to mass transfer at rates ~ $10^2 - 10^3$ times its Eddington accretion rate of ~$10^{-8}$ M☉ yr⁻¹. The neutron star had evidently ejected the vast bulk (~2–3 M☉) of the transferred mass, retaining no more than a few × 0.1 M☉. Ejection like this is expected in the simple picture of Shakura and Sunyaev (1973) of mass expulsion by radiation pressure from thin accretion discs fed at super-Eddington rates<fnref n="3" /> (see Section 3.3 below), so a plausible deduction was that the intense accretion disc wind might block the emitted X-rays from many lines-of-sight, perhaps leaving only narrow escape routes near the disc’s rotational poles. If these had total solid angle $4\pi b$, with $b \ll 1$, an observer viewing the source along these directions – but assuming that the emission was isotropic – would ascribe a luminosity

$$L_{\text{sph}} = \frac{1}{b} L \gg L \qquad (3)$$

to the source, where L is its true luminosity.<fnref n="4" />

Early discussion of this model was hampered by the lack of a clear relation between the beaming factor b and other physical variables, which left a spurious degree of freedom. This gap was later closed (see Section 3.6), but even before this, the model made interesting and testable predictions. Simple estimates showed that the characteristic velocity of the accretion disc wind was ~ 0.1c (King and Pounds, 2003, compare Section 3.14). If $b \lesssim 10^{-2}$, the X-ray fluxes of the early observed ULXs gave $L = bL_{sph} \lesssim 10^{39}$ erg s⁻¹, below the Eddington value for a 10 M☉ black hole. King et al. (2001) also noted that for $b \lesssim 10^{-3}$, the deduced L would be compatible with a neutron star accretor, as expected if Cyg X-2 is a survivor of a ULX phase, and with a white dwarf accretor in cases where the observed X-ray flux was soft (photon energies $\lesssim 0.1$ keV). The brightest ULXs would correspond to the same value $b = 10^{-3}$ of the beaming factor, but with a 10 M☉ black hole accretor.

This idea gave sensible answers for the likely populations and lifetimes of ULXs (King et al., 2001). It suggested that ULXs were an extremely common stage of high-mass X-ray binary evolution characterised by very large mass transfer rates. From the analysis of Cyg X-2's evolution, the best candidate for this stage was the preceding phase of super-Eddington mass transfer. This inevitably follows the standard wind-fed high-mass X-ray binary stage once the donor star fills its Roche lobe, as this is always more massive than the accretor. Then mass transfer shrinks the binary to conserve orbital angular momentum, maintaining very high mass rates limited only by the companion's ability to expand rapidly, usually on its thermal timescale (see Section 3.11). Begelman et al. (2006) noted that the extreme binary SS433

was probably a system of this type, but not viewed along the accretion disc axis, so was probably a ULX ‘viewed from the side’.

The current value of the mass transfer rate $M_{tr}$ in any system like this is (as usual) hard to determine from direct observations. In particular measuring the rate $\dot{P}$ of change of the orbital period P gives at best a wide upper limit $M_{tr} \lesssim |\dot{P}/P| M_2$ to the mass transfer rate (where $M_2$ is the mass of the donor star). This is clearly seen for high-mass X-ray binaries: Falanga et al. (2015) measure $\dot{P}/P \sim -10^{-6}$ yr⁻¹ for LMC X-4, Cen X-3, 4U1538-522 and SMC X-1, which would give maximum mass transfer rates ~ $10^{-6}$ M☉ yr⁻¹ and make all of these systems ULXs, but their observed luminosities are all sub-Eddington ($M_{tr} \lesssim 10^{-8}$ M☉ yr⁻¹).

We now know that a small group of high-mass X-ray binaries with Be-star companions can occasionally produce very high mass transfer rates for short intervals because of dynamical effects on the Be-star disc – see Section 3.15. These systems are unique in repeatedly making transitions from a ‘normal’ X-ray binary to a ULX and back. This property distinguishes very sharply between models for ULXs.

The short lifetimes of both types of high-mass X-ray binary (pre-ULX, or Be-star) showed that the disc-wind beaming model required an association of ULXs with star formation. This was already suspected, in particular following the discovery of 7 ULXs in Chandra observations of the Antennae, where the merger of two galaxies has provoked vigorous star formation (Fabbiano et al., 2001). Chandra observations of the Cartwheel galaxy (Gao et al., 2003) soon strongly confirmed it. They revealed more than 20 ULXs with $L_{sph} > 3 \times 10^{39}$ erg s⁻¹, four of them in the range 1–5.5 × $10^{41}$ erg s⁻¹, in a spreading ring of star formation triggered some $3 \times 10^8$ yr ago when a smaller intruder galaxy crossed the plane of the Cartwheel disc close to its nucleus. None of these features were easily compatible with the IMBH hypothesis for ULXs (King, 2004). The proof that at least some ULXs had stellar masses, and indeed neutron star accretors, as suggested by King et al. (2001) (see the discussion after Eq. (3)), finally arrived when Bachetti et al. (2014) found a coherent periodicity $P = 1.37$ s in the ULX M82 X-2, with $L_{sph} \approx 1.8 \times 10^{40}$ erg s⁻¹. The only plausible interpretation, as the spin of a magnetic neutron star, required $M \sim 1$ M☉, and showed that this object had an apparent X-ray luminosity $L_{sph}$ more than 100 times its Eddington luminosity. Note that this observation did not identify the cause of the neutron star’s super-Eddington luminosity – this is still disputed, essentially between geometrical beaming as suggested by King et al. (2001) and the effects of very strong magnetic fields – see the later part of the present review, and the reviews by Fabrika et al. (2021) and Mushtukov and Tsygankov (2022). Models of this type suggest that a significant, if not dominant, fraction of ULXs are accreting NSs (see abstract in Mushtukov et al. 2015a). It is interesting that only 6 out of the ~ 1800 non-Be ULXs have so far been observed to show pulsing, unlike what we might expect if this dominance holds. (We thank a referee for provoking this insight). A further problem here is that these very strong fields do not seem to be present in binaries which do not have high accretion rates. This picture then requires us to believe that strong-field neutron stars in binaries are unobservable unless the binary companion supplies them with mass at a super-Eddington rate (see Fig. 1).

The observation by Bachetti et al. (2014) removed the main argument motivating IMBH models, and discussion of ULXs is now tightly focussed. The finding of ULXs with $L_X > 10^{41} - 10^{42}$ erg s⁻¹, at least one of which contains a neutron-star accretor is a stringent test of theoretical models. It is compatible with disc-wind beaming, and this is also the only model naturally allowing standard Be-X-ray binaries to become ULXs during bright outbursts.

The main question in discussions of ULXs is whether they are all generically similar, or divide into subgroups with different physical models for each. We discuss this question in detail at the end of this review (Section 4), but note here that in the absence of observations

<footnote n="2">² The possibility of this effect is probably the reason why the occasional $L_X \approx 10^{39}$ erg s⁻¹ outbursts of the Be-star-magnetic neutron star system A0538-66 (White and Carpenter, 1978) did not attract more attention.</footnote>

<footnote n="3">³ Note that advection, neglected by Shakura and Sunyaev (1973), cannot change the accretion luminosity for a neutron star accretor.</footnote>

<footnote n="4">⁴ Although this is almost always called ‘beaming’, as in this review, it is important to note that it is simply collimation, rather than more complex processes such as relativistic beaming.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Soft (0.3–1.5 keV) X-ray image overlaid with the hard X-ray contours (1.577 keV) with sources labelled. Companion galaxies G1 and G2 are labelled as well. The third companion galaxy G3 is ~ 3'; in the northeast, outside the region shown in the figure (Gao et al., 2003).</figcaption></figure>

dividing ULXs into clearly distinct groups, Occam's razor<fnref n="5" /> leaves disc-wind beaming as the only model open to test against all observations of ULXs. For reviews taking differing points of view, see Fabrika et al. (2021), and Mushtukov and Tsygankov (2022), the latter for pulsed ULXs containing very strongly magnetic neutron star accretors.

## 2. Observations

Following the launch of NASA's Chandra (Weisskopf et al., 2000) and ESA's XMM-Newton (Jansen et al., 2001) satellites in 1999, the field of studying X-ray binaries in external galaxies was revolutionised. Both satellites carry X-ray CCD instruments sensitive in the 0.3–10 keV range, and have spectral, timing and imaging capabilities (with on-axis resolutions of < 5"), the latter being vital for successfully disentangling an X-ray bright object from the centre of its host galaxy.

### 2.1. Catalogues & X-ray surveys

After more than two decades of pointed observations, the number of individual sources now totals over 300,000 from Chandra (CSC 2.0: Evans et al. 2010) and over 500,000 from XMM-Newton (4XMM-DR9: Webb et al. 2020) respectively. These extensive catalogues (and their previous iterations) have led to the formation of valuable ULX catalogues, most recently by Kovlakas et al. (2020) utilising Chandra, Earnshaw et al. (2019) utilising XMM-Newton and Walton et al. (2022) using the combination of Chandra, XMM-Newton and Swift (with previous catalogues by Swartz et al. 2004, 2011, Walton et al. 2011, Wang et al. 2016).

In these catalogues, ULXs (numbering ~ 1800 in the most recent: Walton et al. 2022) are located by cross-matching against local

galaxy catalogues (e.g. HECATE: Kovlakas et al. 2021), with checks that sources lie within the host's optical D25 isophote but are not coincident with the central region, and estimates applied for foreground and background (AGN) contamination. These catalogues enable quantification and searches for variability (e.g. propeller states), targeted follow-up observations, and the study of correlations with host type, including star formation rate and metallicity, thereby testing the predictions of binary population synthesis models (e.g. Wiktorowicz et al. 2019, Kuranov et al. 2020).

Both of the most recent ULX catalogues find that spirals are the most common host galaxies – consistent with ULXs being associated with regions of high star formation. An excess of ULXs is found in lower metallicity galaxies (see also Soria et al. 2005, Mapelli et al. 2010, Prestwich et al. 2013, Brorby et al. 2014, Tzanavaris et al. 2016), the numbers scaling with star formation rate as (Kovlakas et al., 2020):

$$N_{ULX} = 0.45^{+0.06}_{-0.09} \times \text{SFR}_{M_{\odot}/\text{year}} + 3.3^{+3.8}_{-3.2} \times \frac{M_{\star}}{10^{12} M_{\odot}} \quad (4)$$

It is important to note that, regardless of the spatial coverage, the long timescale (months or longer) variability of ULXs (as discussed below) demands repeat observing to obtain a reliable census of the population. The all-sky survey by eROSITA (eRASS - Cappelluti et al. 2011) will give precisely these repeat scans, eight being taken over the space of four years (2020–2024) with a far greater sensitivity than the ROSAT all-sky survey (Voges et al., 1999) and coverage from 0.3–10 keV (matching the nominal range of Chandra and XMM-Newton).

The production of such extensive catalogues provides unique statistical insights, notably the X-ray luminosity function (XLF, such as that shown in Fig. 2) giving the combined numbers of low- and high-mass X-ray binary populations and how these depend on metallicity (Lehmer et al., 2021). In turn, these XLFs allow estimates for the impact of accreting binary populations on reionisation and galaxy evolution (Fragos et al., 2013; Kovlakas et al., 2022). XLFs can also be directly compared with population synthesis calculations which include

<footnote n="5">⁵ ‘Do not multiply hypotheses beyond necessity’, or more simply, ‘don’t invent two theories for the same thing.’</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 2. From Mineo et al. (2012) showing the point source X-ray luminosity functions (XLFs), normalised by star formation rate, created from Chandra observations of M101 (blue), NGC 3079 (red), M51 (cyan), Antennae (magenta), NGC 3310 (green).</figcaption></figure>

the evolution of systems through their various phases (including mass transfer, common envelope, SNe and gravitational inspiral). Recent attempts include those by Wiktorowicz et al. (2019) where the effect of geometrical beaming (following equation 55) is included, Khan et al. (2021) where the observational impact of precession is considered (see below) and Kuranov et al. (2020) who obtain the XLF for magnetised neutron stars in the absence of beaming and precession of the wind-cone (see Section 2.3.3).

## 2.2. Spectroscopy

Chandra and XMM-Newton have made it possible to study ULXs with (relatively) high signal-to-noise and energy resolution, and have provided a ‘canonical’ picture of ULX X-ray spectra. As data quality improved, the spectra from these satellites allowed single component spectral models to be excluded, instead indicating the presence of a soft excess and signs of a break to a steeper spectral index at the highest energies, where instrumental effective area typically diminishes (e.g. Kaaret et al. 2003, Miller et al. 2004b, Roberts et al. 2005, Stobartt et al. 2006).

If the soft excess is modelled as a thin disc (set to truncate at the ISCO), the characteristic temperatures are sufficiently low (0.1–0.2 keV) that IMBHs accreting at significantly sub-Eddington rates were proposed (Miller et al., 2004a,c). But these accretion rates imply an upscattered power-law component extending unbroken to very high energies (around the electron’s thermal energy, see e.g. Done et al. 2007). The presence of a break below 10 keV conflicts with this (as do theoretical arguments - see King 2004), supporting the alternative hypothesis of super-critical accretion (Shakura and Sunyaev, 1973). The first large-scale spectral analysis based on a picture of super-Eddington mass supply was performed by Gladstone et al. (2009). Here the harder emission was created by inverse Compton scattering of seed photons (presumed in this case to originate from the soft excess). Since then, the field has benefited from data from newer X-ray satellites, notably NuSTAR, whose ability to image in the 3–70 keV bandpass (Harrison et al., 2013) confirmed the <10 keV break with high signal-to-noise ratio (e.g. Bachetti et al. 2013, Mukherjee et al. 2015) and allowed for the most complete picture of the X-ray SEDs of ULXs (Fig. 3). This wide band-pass has led to important advances in

the field, including the discovery of the first ULX pulsar (Bachetti et al., 2014).

Current modelling of the X-ray spectra of ULXs remains to some degree phenomenological (although more physically motivated models have been developed: e.g. Vinokurov et al. 2013). Models used to describe the soft component tend to be either an advection dominated disc (typically DISKPBB – with a radial temperature profile $T \propto r^{-p}$ with $p < 0.75$, e.g. Mineshige et al. 1994), or a standard disk blackbody (typically DISKBB – Mitsuda et al. 1984). Each of these has been used to describe the wind, inclination and advection-modified emission from the super-critical disc or from an accretion curtain (depending on field strength and accretion rate – Mushtukov et al. 2017). The harder X-ray emission (typically > 1 keV) is variously modelled with Compton components (with empirical assumptions), broken power-laws and, more recently, the addition of emission from a neutron star accretion column (Walton et al., 2018c). Such two- or three-component spectral fits are motivated by spectral timing methods (e.g. the covariance spectrum - Middleton et al. 2015a) and pulse-resolved modelling (Brightman et al., 2016; Walton et al., 2018b). Although multiple components are often used, simple modelling has encouraged the wide adoption of three simple, empirical categories: soft ultraluminous, hard ultraluminous, and broadened dislike (Sutton et al. 2013 – Fig. 4). It is important to note that these definitions are merely descriptive. They do not imply precise statements about the relative contributions of the underlying physical components but provide a useful shorthand for the shape of the spectrum and its evolution with time.

The simple phenomenological modelling described above can be connected to the theoretical picture of super-Eddington mass supply. The characteristic temperatures obtained from modelling the spectra approximate the colour-corrected temperatures $T_{c,ph}$, $T_{c,sp}$ and $T_{c,max}$ of the quasi-blackbody emission from the disc at characteristic radii (Poutanen et al., 2007): the outer photosphere of the wind ($R_{ph}$), spherization radius ($R_{sph}$ - noting that the assumed position of this radius in Poutanen et al. 2007 slightly differs from Eq. (14)), and innermost edge of the disc ($R_{max}$):

$$T_{c,ph} = 0.8 f_{col} \left( \frac{\zeta \beta}{\epsilon_w} \right)^{1/2} m^{-1/4} \dot{m}^{-3/4} \text{ keV} \quad (5)$$

$$T_{c,sp} = 1.5 f_{col} m^{-1/4} \dot{m}^{-1/2} \text{ keV} \quad (6)$$

$$T_{c,max} = 1.6 f_{col} m^{-1/4} \text{ keV} \quad (7)$$

Associated values for the accretion rate then follow under the assumption of a fixed fraction of accretion energy deposited in the wind ($\epsilon_w$), colour temperature correction ($f_{col}$) as well as dynamical values associated with the wind (i.e. the wind velocity relative to the Keplerian velocity at $r_{sph}$, $\beta$, and the cotangent of the wind opening angle, $\xi$ – Poutanen et al. 2007). As an example, typical values of the soft X-ray component’s temperature lie between 0.1 and 0.4 keV (Middleton et al., 2015a); for values of $f_{col} \sim 2$ this indicates Eddington-scaled mass supply rates in the range $\sim 20–300$ for a 10 M$_{\odot}$ black hole and $50–800$ for a 1.4 M$_{\odot}$ neutron star.

One must be cautious when modelling ULX spectra, as descriptions are often degenerate (see Section 2.3.4), and the presence of a neutron star can have a major impact on the interpretation of characteristic temperatures (see Section 2.2 for more detail). As an example, strong dipole fields could truncate the thin disc (i.e. before the local Eddington limit is reached) which comes with the natural corollary that such a disc emits below the Eddington limit for a neutron star mass (see Middleton et al. 2019a). Additional uncertainty in the spectral fitting (regardless of the nature of the compact object) arises as Compton down-scattering within the wind itself can distort the temperature of the photons produced at various radii. It remains likely that the soft component is inclination-dependent, so that a more edge-on system appears increasingly softer until the temperature at the face of the wind dominates (and from
<figure><img src="image_5.png" /><figcaption>Fig. 3. From Pinto et al. (2017) illustrating the variety of observed spectra of ULXs. Each dataset has been unfolded through the instrumental response only (XMM-Newton and/or NuSTAR) and show two spectral states for each ULX (one hard and one soft). The fluxes on the y-axis are arbitrary only.</figcaption></figure>

which we expect a luminosity ~ $L_{Edd}$ to be emitted: Poutanen et al. 2007, Zhou et al. 2019, Qiu and Feng 2021). Together with changes in accretion rate, the pre-dominance along the line-of-sight of the wind photosphere provides a natural explanation for the Soft (and super-soft, e.g. Feng et al. 2016) Ultraluminous Sources which fit well into the simplified unified model of ULX spectra (Poutanen et al., 2007; Middleton et al., 2015a; Soria and Kong, 2016; Urquhart and Soria, 2016a; Pinto et al., 2017).

A notable absence from ULX spectra (with the notable exception of Swift J0243.6+6124: Bykov et al. 2022, see also the claim of a feature in NGC 7456 ULX-1: Mondal et al. 2021) is that of Fe Kα emission, with Walton et al. (2012) placing strong constraints on its presence in Ho IX X-1. This absence is remarkable, given the strong fluorescent

yield of Fe and its typically high abundance. Importantly this may point towards an extremely highly ionized gas environment, consistent with an irradiated conical outflow (see also Medvedev and Fabrika 2010, Pinto et al. 2021). It is quite likely that weak Fe K in reflection would be produced far out in the wind cone where the material is not fully ionized, and also in the disc beyond $R_{\text{sph}}$ provided that this radius is small enough that the outer disc is irradiated by X-rays with energies > 7 keV.

More elaborate models of the emission from slim discs already exist (e.g. Straub et al. 2013), and models for fitting to data may be obtained directly from 3D radiative MHD simulations for super-Eddington mass supply on to black holes (e.g. Narayan et al. 2017) and magnetised neutron stars (Kawashima and Ohsuga, 2020). Still more elaborate spectral modelling may include self-consistent photoionization.
<figure><img src="image_6.png" /><figcaption>Fig. 4. From Kaaret et al. (2017) showing the variety of phenomenological states of ULXs (SSUL: supersoft ultraluminous; BD: broadened disc; H/SUL hard and soft ultraluminous) as defined by Sutton et al. 2013 as compared to the standard canonical soft (thermal) and hard state spectra of XRBs.</figcaption></figure>

### 2.2.1. High energy-resolution X-ray spectroscopy

Phenomenological models can reasonably describe ULX spectra as observed by current instrumentation, although clear residuals to these models are found between 0.5 and 2 keV (Stobbart et al., 2006) in multiple ULXs (Middleton et al., 2015c), and seen in both XMM-Newton and Chandra spectra (ruling them out as instrumental artefacts - Roberts et al. 2006, Pintore et al. 2015). Initially thought to originate in diffuse star-forming plasma (Strohmayer et al., 2007; Strohmayer and Mushotzky, 2009; Miller et al., 2013), the luminosities are too high for this interpretation (Middleton et al., 2014). Instead, both via high angular resolution studies with Chandra (Sutton et al., 2015), and by studying the correlated variability of the features with spectral hardness (Middleton et al., 2015c; Pinto et al., 2020b), it was determined that they must be associated with the ULX (and not with the surrounding nebula, e.g. Sathyaprakash et al. 2019a).

Instruments with higher energy resolution (e.g. XMM-Newton's Reflection Grating Spectrometer, RGS) have shown that the features detected at low energy resolution in the CCD-quality spectra are composed of blueshifted absorption lines and rest-frame emission lines (Pinto et al. 2016, 2017 – but see Kosec et al. 2018a for apparently blueshifted emission lines in NGC 5204 X-1). In NGC 1313 X-1, counterpart Fe absorption lines have also been found in NuSTAR CCD-quality data, with a blueshift consistent with that inferred from the lower energy resolved lines (Walton et al., 2016b). The limit implied on the outflow velocity is > 0.2 c, considerably faster than observations of thermal (or magneto-centrifugal) winds from sub-Eddington XRBs (e.g. Ponti et al. 2012), but consistent with radiatively driven winds from a super-critical inflow (e.g. King and Pounds 2003, Jiang et al. 2014, Kitaki et al. 2021).

Ultrafast outflows have now been located in several ULXs via individual searches of RGS spectra (NGC 1313 X-1 and NGC 5408 X-1: Pinto et al. 2016, NGC 55 ULX: Pinto et al. 2017, NGC 247 ULX-1: Pinto et al. 2021) and notably in the transient NS ULXs, NGC 300 ULX-1 (Kosec et al., 2018b) and Swift J0243.6+6124 (van den Eijnden et al., 2019). Intriguingly this latter source simultaneously launches both a wind and a jet (see Section 2.4.1). In addition to these studies, a larger scale search using a faster algorithm of 19 ULXs by Kosec et al. (2021) indicates that rest-frame emission lines and blue-shifted absorption lines are indeed a ubiquitous feature.

The velocity of the wind, $v_{wind}$, a measure of the ionisation state, $\xi$ (both obtained by using consistent models for photoionised plasma –

Pinto et al. 2016), and an assumption for the ionising luminosity, allows the kinetic luminosity of the wind to be estimated from:

$$P_{\text{wind}} = \frac{1}{2} \dot{m} v_{\text{wind}}^2 = 2\pi L_{\text{ion}} \frac{\Omega}{\xi} v_{\text{wind}}^3 m_{\text{p}} \quad (8)$$

The filling and covering factor of the wind ($\Omega$) in the above formula remain unknown, but the high projected velocities (which are therefore lower limits on the true velocity) imply that the kinetic luminosities for these winds are substantial, and may even dominate over the observed radiated power from the source (Pinto et al., 2016). It is worth noting that such observations are potentially in conflict with MHD models which predict lower fractions of the radiative power spent in launching the wind (Jiang et al., 2014). Questions remain over the nature of the rest-frame emission lines – conceivably these may be due to collisional ionization as a result of interaction with the wind of the companion star (e.g. Oskinova 2005, Mauerhan et al. 2010) or may be associated with lower velocity winds driven from an outer, irradiated disc (Middleton et al., 2021a).

Studies have recently shifted to determining how the wind appears to change as a function of spectral hardness. As first reported in Middleton et al. (2015c), the strength of the CCD residuals in NGC 1313 X-1 appear to anti-correlate with increasing spectral hardness (i.e. the features are weaker when the source appears spectrally harder). The interpretation was that the wind or system was precessing, or the mass accretion rate changing (the optical depth along our line-of-sight through the wind changing with either/both: Poutanen et al. 2007). A more recent in-depth study of NGC 1313 X-1 based on a series of long XMM-Newton observations, has resolved the atomic features in emission and absorption (using the RGS) and made a direct comparison to changes in the X-ray spectrum (Pinto et al., 2020b). These data have allowed for important progress to be made, including the discovery of a slower, more neutral component of the wind (moving with a projected velocity of ~ 0.08c rather than the fast wind of ~ 0.2c), pointing towards a complex outflow (as one might expect from simulations - e.g. Takeuchi et al. 2013). These high energy-resolution results broadly confirm the picture from the CCD-quality residuals (as does the larger sample result of Kosec et al. 2021), with the observation of an anti-correlation between spectral hardness and line strength, possibly driven by changes in accretion rate and/or inclination angle via precession (Middleton et al., 2015b,d). Importantly the correlated changes further confirm that the winds are associated with the ULX rather than some surrounding, more distant material (which would not respond on such timescales). The high quality data made available via long RGS exposures has allowed, for the first time, photoionization modelling of a ULX wind (Pinto et al. 2020a, see also Pinto et al. 2021 for a comparison to other ULXs), indicating that the outflowing material in the innermost regions is potentially thermally unstable. In future, the abilities of XRISM and Athena (with its X-ray Integral Field Unit) will allow for far deeper studies of ULX winds.

### 2.2.2. Spectral studies of neutron star ULXs

The discovery of several pulsating neutron star ULXs (Bachetti et al. 2014, Israel et al. 2017a,b, Fürst et al. 2016, Tsygankov et al. 2017, Carpano et al. 2018, Sathyaprakash et al. 2019b, Doroshenko et al. 2018 – see Table 2) has led to scrutiny of the X-ray spectrum in light of the expected differences from the standard super-critical model as a consequence of the presence of the magnetosphere – accretion curtain and column – and surface. The truncation of the disc at the magnetospheric radius and free-fall of material onto the magnetic poles, leads to two likely configurations depending on the accretion rate. At very high accretion rates or low dipole field strengths, the spherization radius is reached at radii greater than the magnetospheric radius – if the wind is optically thick (see Vasilopoulos et al. 2019, Middleton et al. 2019b) it will act to collimate some of the radiation within. Regardless, we expect thermal emission from $R_{sph}$ down to $R_{mag}$, which tends increasingly towards a blackbody as $R_{sph}$ tends towards $R_{mag}$. At radii
<figure><img src="image_7.png" /><figcaption>Fig. 5. Pulse-on versus pulse-off spectra for NGC 7793 P13 from Walton et al. (2018b) using data from NuSTAR.</figcaption></figure>

larger than R<sub>sph</sub>, the emission is from an unmodified, approximately thin disc (thin below ~ 30% Eddington, e.g. McClintock et al. 2006), although this emission is probably modified and scattered by the wind covering those regions between R<sub>sph</sub> and R<sub>ph</sub> (Poutanen et al. 2007). For strong magnetic fields or relatively low accretion rates, R<sub>sph</sub> < R<sub>mag</sub> and the disc emission beyond R<sub>mag</sub> originates from an Eddington-limited thin disc. In either accretion-rate regime, within R<sub>mag</sub>, the material falls freely along magnetic field lines until meeting the standing shock front in the accretion column. The key characteristic of the magnetospheric curtain is the optical depth τ which depends on the accretion rate (and therefore dipole field strength) and angle from the accretion disc plane. According to Mushtukov et al. (2017) this gives:

$$ \tau \approx \frac{70 L_{39}^{6/7} B_{12}^{2/7}}{\beta(\theta)} \left( \frac{\cos \theta_0}{\cos \theta} \right)^3 \quad (9) $$

where θ<sub>0</sub> ≈ π - (R/R<sub>M</sub>)<sup>1/2</sup> is the angle to the accretion column base and β(θ) = v(θ)/c is the local dimensionless velocity along the magnetic field lines. This optical depth determines the shape of the thermal spectrum from the curtain (becoming optically thick for L<sub>39</sub> > B<sub>12</sub><sup>1/4</sup>) as well as the effect of dilution of the variability (Mushtukov et al. 2019), with a range of emergent photon temperatures due to the inclination dependence of the optical depth. A rough estimate for the temperature of the curtain is provided by Mushtukov et al. (2017) as T<sub>out</sub> ~ 0.5L<sub>39</sub><sup>11/28</sup>B<sub>12</sub><sup>2/7</sup>m<sup>-1/14</sup>R<sub>6</sub><sup>-5/7</sup> keV (where R<sub>6</sub> = R/10<sup>6</sup> cm and m = M/M<sub>☉</sub>) which can be ≥ 1 keV for typical ULX luminosities, and can therefore explain one of the hard components in the spectra of ULXs.

The final set of likely spectral components is then: at high accretion rates/low field strength – thin disc, thick disc and wind, curtain and column, and for low accretion rates/high dipole field strength – thin disc, curtain and column. Using the latter description of the inflow, Koliopanos et al. (2017) explored the possible range of dipole magnetic field strengths for 18 ULXs (containing both pulsing ULXs and those without) under the assumption that R<sub>sph</sub> < R<sub>mag</sub>, finding values ranging from ~ 10<sup>11</sup>–10<sup>13</sup> G (although the assumption of a thin disc appears to be at odds with the observed super-Eddington luminosity from the soft component).

The spectrum of the pulsed emission can also be extracted via selecting ‘off’ versus ‘on’ pulse times and obtaining the difference spectrum. The latter has been performed by Brightman et al. (2016)

in the case of M82 X-2 and Walton et al. (2018b) in the case of NGC 7793 P13 (Fig. 5), who found the accretion column to be well-modelled by a power-law with an exponential cutoff (see also Potts Schmidt et al. 2016 for a similar finding for SMC X-3). The unambiguously spectrally hard nature of the pulsed emission has naturally led to suggestions that other, similarly hard ULXs (e.g. Ho IX X-1: Walton et al. 2017 and IC 342 X-1: Middleton et al. 2015a) may also harbour neutron stars (Pintore et al., 2017; Koliopanos et al., 2017; Walton et al., 2018c; Gúrpide et al., 2021). Showing whether the hard excess is indeed due to the presence of an accretion column will be valuable for identifying neutron star ULXs without the need for deep pulsation searches. This in turn will help reveal ULX demographics and test our grasp of binary evolution.

### 2.2.3. Cyclotron resonance scattering features

In the presence of strong magnetic fields, electron and proton orbits are quantised into Landau levels; such particles can resonantly scatter incident photons of sufficient energy depending on the quantum mechanical cross section (see Schwarm et al. 2017) and leave cyclotron resonance scattering features (CRSFs) imprinted on the spectrum (although any subsequent scattering, for instance in an accretion curtain, can reduce their detectability - see Mushtukov et al. 2017). The transition energy of such lines is related to the magnetic field strength and particle type according to:

$$ \Delta E_e^- = \frac{11.6}{(1+z)} \left( \frac{B}{10^{12} G} \right) \text{ keV} \quad (10) $$

$$ \Delta E_{p^+} = \frac{6.3}{(1+z)} \left( \frac{B}{10^{15} G} \right) \text{ keV} \quad (11) $$

where z is the gravitational redshift, given by:

$$ z = \left( 1 - \frac{2GM}{r_{\text{cyc}}c^2} \right)^{-1/2} - 1 \quad (12) $$

where r<sub>cyc</sub> is the distance from the surface of the neutron star to where the line is formed. Studies of CRSFs in Galactic systems (e.g. Tsygankov et al. 2006, Jaisawal and Naik 2016) have been useful for studying the regions in which these lines are formed (see the review of Staubert et al. 2019).

Brightman et al. (2018) reported the first CRSF in a ULX, at an energy of ≈ 4.5 keV in a Chandra observation of M51 ULX-8 (see Fig. 6); as this energy does not correspond to known instrumental features and cannot be readily explained as a highly blueshifted absorption line, a CRSF is a possible explanation. Given the line energy, the implications are either electrons orbiting in a field with a strength of ~ 5 × 10<sup>11</sup> G (for z = 0.25) or protons orbiting in a field around 9×10<sup>14</sup> G. Given the narrow width of the line, it was concluded by Brightman et al. (2018) that the feature is a proton CRSF. Combined with constraints from spectral fitting, it was concluded by Middleton et al. (2019a) that the dipole field is likely to be weak with the CRSF the result of a stronger, higher order multipolar (e.g. quadrupole) component closer to the neutron star (e.g. Israel et al. 2017a,b, Brice et al. 2021), a situation now confirmed in the case of Swift J0243.6+6124 (Kong et al., 2022).

At the time of writing, there is only a single other claim of a CRSF in a ULX, in the transient system, NGC 300 ULX-1 which was identified as a neutron star from its 31.6 s period (Carpano et al., 2018). The presence of a CRSF was inferred from spectral modelling of the pulsed component (Walton et al., 2018a). Although the line energy (≈ 13 keV) appears consistent with the estimate of the dipole field strength from spin-up, doubts have been raised by Koliopanos et al. (2019) as to the presence of the line, as applying a spectral model based on emission from the accretion curtain and additional hard power-law tail, do not argue for its inclusion.
<figure><img src="image_8.png" /><figcaption>Fig. 6. Figure adapted from Brightman et al. (2018) showing the location of ULX-8 in M51 (left) and the CRSF identified at ~4.5 keV in a Chandra observation (right).</figcaption></figure>

### 2.2.4. Spectral evolution

ULXs are well known to change in both brightness and in spectral shape over timescales of months or less. Most ULXs tend to have a fairly predictable behaviour, with a given source becoming brighter as the hard and soft spectral components both increase in brightness, with the harder component tending to increase more than the soft (see e.g. Middleton et al. 2015a) although this depends on how the spectrum is modelled (as this delineates the hard and soft components). As a result, for some of the brightest ULXs, there tends to be a clear correlation between the source becoming spectrally harder (with more flux emerging above ~1 keV compared to that below) and brighter (Weng and Feng, 2018) although this is not universal (see Gúrpide et al. 2021).

The coupled changes in X-ray spectrum and observed luminosity can be explained by a given ULX undergoing a change in accretion rate, inclination angle (due to precession) or a combination of the two (Middleton et al., 2015a). Under the condition that the neutron star dipole field is weak, so that $r_{sph} < r_M$, or the compact object is a black hole, increasing the accretion rate increases the spherization radius and the optical depth of the wind (as it is an integrated quantity - see Poutanen et al. 2007) and decreases the temperature at both the spherization radius and outer photospheric radius.

Accordingly, we should expect to see an anticorrelation between the soft component's luminosity (which may be driven by radii interior to $R_{sph}$) and the temperature at the spherization radius. Identifying and interpreting a correlation (or lack of one) is not necessarily straightforward, as increasing $R_{sph}$ or moving to a more face-on orientation leads to steadily increased beaming of the photons interior to $R_{sph}$ (Walton et al., 2017, 2020) which can make spectral modelling difficult.

Conversely, the hard X-ray emission is not expected to be particularly sensitive to accretion rate (see Eq. (2)) but is the region deepest within the wind-cone and therefore experiences the largest amount of collimation and beaming (this is an important point when considering the effects on the local environment - see Section 2.6). An increase in accretion rate is only likely to affect the luminosity of the hard component if the wind opening angle is sensitive to accretion rate (c.f. Jiang et al. 2014, 2019). For a fixed inclination where sightlines graze or intercept the wind, this picture is further complicated by obscuration and reprocessing (Middleton et al. 2015a, although the effect of absolute collimation at any angle can be inferred numerically - Dauser et al. 2017).

Given the observed precession of SS433 (known to be highly supercritical and established to be a ULX seen edge-on - see Section 2.8.3),

it is reasonable to expect precession to occur in other ULXs as well (for a variety of reasons, and with some support from observation e.g. Pasham and Strohmayer 2013, Middleton et al. 2015c, Weng and Feng 2018). It is therefore important to consider the effect such a changing effective inclination of the super-critical flow might have (see Middleton et al. 2015a for details). Also, whilst changes in accretion rate and/or inclination lead to relatively predictable coupled evolution in spectrum and brightness, we also expect an impact on the timing signatures (see Middleton et al. 2015a), and detected strength and speed of winds (see Pinto et al. 2020a, Pinto et al. 2021). Such additional considerations are vital as they are independent of the spectral model used to describe the emission which can often be degenerate (see Section 2.3.4).

The majority of the spectral evolution observed in ULXs can be fitted into the above picture. Notably the model can broadly explain observations of super soft ULXs, which have extremely low temperature soft components ($kT_{bb} \approx 50–150$ eV) and very little emission above 1 keV (e.g. Kong and Di Stefano 2003, Urquhart and Soria 2016a, Soria and Kong 2016) and which can be explained by an edge-on orientation and/or very high accretion rates (where the wind is highly obscuring). Conversely, the brightest and spectrally hardest ULXs, tend to become slightly softer as they brighten, consistent with increased beaming of emission within $r_{sph}$ (Walton et al., 2017, 2020) and where the ULX is viewed at lower inclinations. In spite of the success of this simple model, a number of important open questions remain.

One of the most obvious questions facing the unified inclination/accretion rate model of ULXs is how the presence of a strong dipole field affects the flow. For the most part, and as long as $R_M \ll R_{sph}$, the expectation is that the spectrum will evolve along similar paths, as the position of the magnetospheric radius is expected to be independent of the mass-transfer rate for classical supercritical flows (in the absence of advection, see Chashkina et al. 2017, 2019 for the effects of including advection). Rather more drastic deviations from the model above occur for those NS ULXs (note not necessarily ULX pulsars as there are many reasons why a pulsation may be absent - e.g. Mushtukov et al. 2020) which are close to spin equilibrium, as a propeller state can switch off the accretion curtain and column, leading to a quenching of the hardest spectral components (while the other components remain ‘on’). There are already hints of propeller states occurring in one ULX pulsar (Tsygankov et al., 2016b) and such states might account for some of the more highly variable ULXs (e.g. Earnshaw et al. 2018, 2020, Song et al. 2020, although a propeller state is ruled out for recent large-scale
changes in X-ray flux observed in NGC 300 ULX: Vasilopoulos et al. 2019, and NGC 7793 P13: Fürst et al. 2021).

One of the most intriguing observations is that the shape and normalisation of the emission above 10 keV tends to stay remarkably constant in a number of bright ULXs as their spectra evolve (see Walton et al. 2020). This can be explained so long as the innermost regions are maximally beamed and any changes in accretion rate do not change the opening angle of the wind cone appreciably (although see Jiang et al. 2019 for indications that the cone's opening angle may well respond to larger changes in accretion rate). An alternative possibility is that the emission at these energies originates from a stable accretion column, and applying a phenomenological model for the column can indeed provide a good description of the hard emission (Pintore et al., 2017; Walton et al., 2018c).

There are important questions remaining over how the individual spectral components fit into the above picture, typically explored through use of the temperature-luminosity plane. Such analyses have a long history in the study of ULXs, with indications of a strong negative L-T correlation ($L \propto T^{-3.5}$) reported across a sample of ULXs (e.g. Kajava and Poutanen 2008, 2009). Given the possibility that the ULX population is heterogeneous, it is important to also study sources individually. This is becoming increasingly feasible with intense observing campaigns.

Walton et al. (2020) examined nine observations of NGC 1313 X-1 (combining XMM-Newton with NuSTAR data) and fitted the spectra with a two-component thermal model of a disc black-body at low energies and an advection dominated disc at high energies (and an additional component to account for the hard excess). From this they found a positive (and relatively steep) correlation between the temperature and luminosity of the hotter, advection dominated disc (the soft component has a flatter, possibly negative correlation). This correlation is further subdivided into high and low luminosity branches and may provide an intriguing test of what processes drive the changing appearance of ULXs. A similar analysis has been performed by Robba et al. (2021) for the ULX pulsar NGC 1313 X-2, although a second luminosity-temperature branch is not as evident as for NGC 1313 X-1. As shown in several works (Luangtip et al., 2016; Gúrpide et al., 2021), the presence and nature of any correlation (at low or high temperatures) can be affected by the choice of spectral model, reinforcing the important role other methods can play in understanding the nature of the flow and the components in the X-ray spectrum (see Section 2.3.4).

### 2.2.5. Open questions and future directions

Important questions remain regarding the modelling of ULX continua, including:

* the apparent stability at energies above the spectral break (Walton et al., 2020)

* the presence of fainter spectral components which may correspond to bremsstrahlung emission, as seen from SS433's jets (e.g. Walton et al. 2015)

* the trends of luminosity versus temperature for the individually modelled components (Kajava and Poutanen, 2009; Walton et al., 2020).

* the impact of spin and advection on the hard spectral emission

* the prevalence (or lack thereof) of CRSFs

Many of these issues may be resolved as self-consistent spectral models derived directly from 3D simulations of black hole and neutron star ULXs become available, (e.g. Ohsuga and Mineshige 2011, Jiang et al. 2014, Sądowski et al. 2015, Sądowski and Narayan 2016, Narayan et al. 2017, Takahashi et al. 2018, Kawashima and Ohsuga 2020).

## 2.3. Timing analysis

X-ray spectroscopy has allowed for a broad categorisation of ULXs, but as with other accreting systems, timing properties (over short and long timescales) are just as important for testing models and for revealing new behaviour.

### 2.3.1. Short timescales

Because of its large effective area, XMM-Newton has been the favoured instrument for studying the short (intra-observation) timescale variability of ULXs. The simplest analysis uses the Fast Fourier Transform (FFT) of the lightcurve to obtain the periodogram (or, if averaged, the power spectral density – PSD). As this is typically normalised to be in (rms/mean)$^2$/Hz units, the integral of the PSD gives the fractional root mean square (rms), which, when Poisson noise corrected and plotted against energy, yields the rms spectrum (see Uttley et al. 2014 for a comprehensive review).

The PSDs extracted from single observations of ULXs are generally described by either a red-noise power law (power scaling as $v^{-1}$ → $v^{-2}$) or a broken power-law (e.g. Strohmayer and Mushotzky 2003, Dewangan et al. 2006, Strohmayer et al. 2007, Heil et al. 2009, Atapin et al. 2019), designated as the broadband noise. Binning the rms reveals a tentative linear relationship with X-ray flux (Heil et al., 2009; Heil and Vaughan, 2010; Hernández-García et al., 2015) as expected from multiplicative propagation of fluctuations in accretion rate (Lyubarskii, 1997; Churazov et al., 2001; Arévalo and Uttley, 2006). Heil et al. (2009) observed that the variability of some ULXs appeared very low, while conversely others appeared extremely high (up to 30% fractional rms – Middleton et al. 2015a). This can be explained via either the dilution of intrinsic variability from a neutron-star accretion column as photons scatter inside the accretion curtain before eventually escaping (Mushtukov et al., 2019) or by mass loss in the wind, which may remove variability from the inflow (Middleton et al., 2015a). The substantial amounts of short timescale (<10 ks) variability in those spectrally softer sources – in which there is strong evidence for mass-loaded winds (e.g. NGC 5408 X-1 and NGC 6946 X-1) – could arise from stochastic obscuration by clumps in the wind (Middleton et al., 2011; Sutton et al., 2013), produced via radiative hydrodynamic (e.g. Takeuchi et al. 2013) or Kelvin Helmholtz instabilities.

In addition to the broadband noise, there have been several claims of quasi-periodic oscillations (QPOs – with power concentrated over a small frequency range) at 10 s of mHz (NGC 5408 X-1: Strohmayer and Mushotzky 2009, Pasham and Strohmayer 2012, NGC 6946 X-1: Rao et al. 2010, NGC 1313 X-1: Pasham et al. 2015, M82 X-1: Strohmayer and Mushotzky 2003, IC 342 X-1: Agrawal and Nandi 2015) as well as at higher (~Hz) frequencies in M82 X-1 (Pasham et al., 2014). There have been various attempts to use QPOs to obtain the accretor mass, based upon a simple frequency scaling (e.g. Strohmayer and Mushotzky 2009, Pasham et al. 2015), which tend to give masses typical of IMBHs. However, while QPOs are detected from supercritically accreting systems such as TDEs (e.g. Pasham et al. 2019) and in GRS 1915+105 (the well known 67–68 Hz QPO: Belloni and Altamirano 2013), it is not clear that there are well-studied objects with well-understood signals available for direct comparison. At present, the origin of ULX QPOs remains unclear, although models invoking modulations of accretion rate (Okuda et al., 2009) or precession of the radiation-pressure supported disc – which can also explain the lag seen around the QPO frequency (De Marco et al. 2013 – see Section 2.3.4) have both been put forward (Middleton et al., 2019b). With the high throughput of forthcoming instruments (namely ESA's Athena), studies of ULX QPOs will become considerably more revealing.

### 2.3.2. Pulsations

Following the discovery of pulsations from M82 X-2 (Bachetti et al., 2014), one of the most important tools now regularly applied to ULX data is the accelerated pulsation search (see e.g. Ransom 2001, Ransom et al. 2002, Andersen and Ransom 2018). This can account for the orbit of the neutron star about the system barycenter, which acts to smear the power in the pulsations over neighbouring Fourier frequency bins. Accelerated searches within both new and archival data have led to the discovery of several more pulsating ULXs (NGC 5907 ULX-1: Israel et al. 2017a; NGC 7793 P13: Fürst et al. 2016, Israel et al. 2017b; SMC X-3: Tsygankov et al. 2017; NGC 300 ULX-1: Carpano et al.
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** From Bachetti et al. (2020) showing the positive glitch (with a change in frequency $\Delta\nu \sim 8 \times 10^{-11}$ Hz/s) in M82 X-2, revealed by extensive monitoring with NuSTAR.</figcaption></figure>

2018 – first identified as a supernova by Monard 2010; NGC 1313 X-2: Sathyaprakash et al. 2019b; M51 ULX-7: Rodríguez Castillo et al. 2020; RX J0209.6-7427: Vasilopoulos et al. 2020b, Chandra et al. 2020) and Swift J0243.6+6124: Kennea et al. 2017.

The pulsing ULXs typically have periods ranging from ~1–30 s (see Table 2 and Fig. 8). As is apparent from Table 2, the secular $\dot{P}$ in these systems (i.e. $\Delta P/\Delta T$ where $\Delta T$ is typically months or longer – except in the case of NGC 300 ULX-1) reveals a collection of objects which tend to spin up at extremely high rates (notably NGC 300 ULX: Carpano et al. 2018) when compared to Galactic X-ray pulsars (although M82 X-2 is observed to spin down between 2014 and 2016: Bachetti et al. 2020). The secular spin up/down rates may differ greatly from the ‘instantaneous’ $\dot{P}$ obtained over the course of an observation due to the influence of orbital dynamics and changes in accretion rate between observations (Israel et al., 2017a; Bachetti et al., 2020). In addition to these long timescale trends, glitches (perhaps associated with the coupling of differentially rotating inner crust – and possibly core – to outer crust: Anderson and Itoh 1975) have now been observed in M82 X-2 (Bachetti et al. 2020 – Fig. 7). These imply a change in spin frequency $\delta\nu/\nu \approx 10^{-5}$, similar in magnitude to glitches seen in magnetars and a small number of accreting neutron stars (Galloway et al., 2004; Dib et al., 2009; Dib and Kaspi, 2014; Archibald et al., 2016; Serim et al., 2017). Interestingly, there are two (or possibly three) anti-glitches in NGC 300 ULX-1 (Figure 8, Ray et al. 2019), where the neutron star undergoes a dramatic spin-down $\delta\nu/\nu > -10^{-4}$. This may be the result of the enormous spin-up rate creating a sizeable lag between the rotation of the crust and superfluid interior (Ray et al., 2019), an idea which is supported by the lack of corresponding radiative changes during the glitch. Only spin-up glitches have been observed in non-accreting pulsars, but accreting pulsars (Klochkov et al., 2009) and magnetars show both spin-up and spin-down glitches (e.g. Dib and Kaspi 2014). The anti-glitches seen in NGC 300 ULX-1 are the largest yet-observed in any kind of pulsar system.

The known ULX pulsars have pulse fractions (defined as $[\text{Pulse max} - \text{Pulse min}]/[\text{Pulse max} + \text{Pulse min}]$) which increase with photon energy (as seen in Galactic X-ray pulsars: Lutovinov and Tsygankov 2009), and take values from <5%–100% depending on the observation,

source and energy band (see Table 2 for more details). Notably, Fig. 9 shows that several of the ULX pulsars have an a highly sinusoidal pulse profile (although both NGC 300 ULX-1 and SMC X-3 are notable exceptions), unlike that expected from a pencil beam geometry. This suggests a large emitting region, either a fan-beam geometry (Basko and Sunyaev, 1976) or accretion curtain, the latter reprocessing (and therefore potentially diluting) the pulsed emission (Mushtukov et al., 2017, 2019).

We can obtain a rough indication of the statistical significance of a given pulsation (with a width $\Delta\nu$, fractional rms, $r$ – related to the pulse fraction – in an observation of length $T$ from a source of count rate $I$, with background count rate $B$) from $n_{\sigma} = 1/2 [I^2/(I+B)]\sqrt{T/\Delta\nu}$ (Lewin et al. 1997. (A true measure of the significance relies on considering both the noise and the number of independent free trials). Of the ~500 known ULXs (e.g. Earnshaw et al. 2019), the vast majority do not have sufficient data quality for deep pulse searches (down to small $r$). In those with sufficient data quality and not already identified as PULXs, pulsations appear to be absent down to the detection threshold (e.g. Doroshenko et al. 2015, Israel et al. 2017b).

There are several reasons why pulsations may be absent in a given ULX or in a given observation. Pulsations are observed to be transient in the known ULX pulsars (and indeed, in Galactic accreting X-ray pulsars as well, e.g. Kretschmar et al. 2000) for reasons not yet well understood (see Bachetti et al. 2014, Israel et al. 2017a, Sathyaprakash et al. 2019b, Rodríguez Castillo et al. 2020, Bachetti et al. 2020), and conditions during a given observation may not be favourable for detecting the pulsed emission from the column (and/or curtain). As an example, we may be prevented from observing the pulsations by dilution by a stronger stable component in the spectrum (Walton et al. 2018b, e.g. if the disc/wind has precessed or accretion rate has changed Pasham and Strohmayer 2013, Middleton et al. 2015a,c, 2018, 2019b). This can also happen if the neutron star has free-body precessed (Fürst et al., 2017) so that the pulses are considerably weaker or beamed away from us. It may also be that we will never observe pulsations in some neutron star ULXs as we may never observe at suitable angles to view the accretion column. Alternatively the field may have been suppressed by accretion (e.g. Urpin and Geppert 1995), the neutron star and inner disc may be aligned (King and Lasota, 2020), scattering in the wind cone may have suppressed pulsations on timescales shorter than the light-crossing time (Mushtukov et al., 2020), or some ULXs may simply harbour black holes (see Middleton and King 2017). Binary evolution would seem to make the latter inevitable to some degree - (Wiktorowicz et al., 2019).

### 2.3.3. Long timescales

It is well accepted that LMXBs undergoing recurrent outbursts (see Hameury and Lasota 2020) may appear as ULXs for relatively short periods of time (e.g. Middleton et al. 2013, Carpano et al. 2018, van den Eijnden et al. 2018). Many bright ULXs are persistent (although we note the discovery of a new transient in M51 reaching $10^{40}$ erg/s: Brightman et al. 2020a), but their luminosities can fluctuate by orders of magnitude, periodically/quasi-periodically (on timescales of 10 s of days: Kaaret et al. 2006, Kong et al. 2016, Walton et al. 2016a) or seemingly at random. There are several possible explanations for such variations: changes in accretion rate driven by effects in the outer disc (see Fig. 19, Section 3.16 and Hameury and Lasota 2020, Middleton et al. 2021a), an accreting neutron star entering a propeller phase, orbital variability as a neutron star passes through the decretion disc of a Be star, and precession of the disc/wind under external torques. In a number of the pulsating ULXs, the orbital period can be deduced from the Doppler modulation of the pulsations. This gives periods of a few days (see Table 2), with NGC 7793 P13 having a proposed period of 64 days (Motch et al. 2014 although see Fuerst et al. 2021 for details regarding the periodicities in this ULX). Other than NGC 7793 P13, the ~ day orbital periods are considerably shorter than the ‘super-orbital’ periods typically revealed by Swift monitoring but
<figure><img src="image_10.png" /><figcaption>Fig. 8. Data from Ray et al. (2019), highlighting the glitches in the ULX pulsar NGC 300 ULX-1. The upper (orange) points show offset residuals from the best-fitting model of the frequency evolution without glitches. These anti-glitches (with $\delta v/v > -10^{-4}$) are the largest yet observed in any pulsar system, and may be a result of the lag between the rotation of the crust and superfluid interior, driven by the enormous spin-up rate (Table 2).</figcaption></figure>

in some cases observed in the optical/UV (e.g. NGC 7793 P13: Fuerst et al. 2021). Both sets of periods are reported in Table 2.

Precession has been proposed as one means by which to explain the ~10 s of days super-orbital variability of ULXs by a number of authors. Pasham and Strohmayer (2013) suggested that the changes in the spectrum of M82 X-1 associated with a 62 day period (Kaaret et al., 2006; Kaaret and Feng, 2007) could result from changes in the projected areas of the inner disc (perhaps driven by a radiative warp: Pringle 1996). Middleton et al. (2015c) suggested that changes in the strength of atomic features, now known to be associated with the wind (Pinto et al., 2016, 2020b), could result from precession of the disc and wind in NGC 1313 X-1, and Luangtip et al. (2016) appeal to a combination of precession and changes in accretion rate in order to explain the changing appearance of Ho IX X-1.

There are a number of possible mechanisms for driving a changing view of the inner regions of ULXs including precession of a neutron star's magnetic dipole (Lipunov and Shakura, 1980):

$$t_{\text{mag,prec}} \approx 1.5 \times 10^4 \frac{\mu_{30}^{-2} I_{45} R_{M.8}^3 P^{-1}}{\cos \phi (3 \cos \zeta - 1)} \text{ years} \quad (13)$$

where $\mu_{30} = B_{12}R_{6}^{3}$, $B_{12}$ is the surface field strength in units of $10^{12}$G, $R_{6}$ is the neutron star radius in units of $10^{6}$cm, $I_{45}$ is the neutron star moment-of-inertia in units of $10^{45}$ gcm$^2$, $\phi$ is the angle between the normal to the disc plane and the spin axis, and $\zeta$ is the angle between the magnetic axis and the spin axis. To match observations of ULX superorbital periods of months or less, field strengths > $10^{14}$G are required (Mushtukov et al., 2017; Vasilopoulos et al., 2020a).

Should the accretion flow be misaligned with the neutron star or black hole spin axis (which seems quite plausible: King and Lasota 2020) then frame-dragging of the fluid induces a torque. Similar to the hot inner flow of X-ray binaries, a thick super-critical disc/wind is expected to have an effective viscosity parameter $\alpha < H/R$. The frame-dragging induced torque is then communicated as a bending wave from the inner edge of the flow out to the photospheric radius (Middleton et al., 2018, 2019b) which may then induce solid body precession (Fragile et al., 2007) on an estimated timescale of:

$$t_{\text{LT,prec}} \sim \frac{GM\pi}{3c^3 a_*} r_{\text{sph}}^3 \left[ \frac{1 - \left(\frac{r_{\text{in}}}{r_{\text{sph}}}\right)^3}{\ln\left(\frac{r_{\text{sph}}}{r_{\text{in}}}\right)} \right] \text{s} \quad (14)$$

for the disc and

$$t_{\text{LT,prec}} \sim \frac{GM\pi}{3c^3 a_*} r_{\text{ph}}^3 \left[ \frac{1 - \left(\frac{r_{\text{in}}}{r_{\text{ph}}}\right)^3}{\ln\left(\frac{r_{\text{ph}}}{r_{\text{in}}}\right)} \right] \text{s} \quad (15)$$

for the wind (here r is in units of gravitational radius ($GM/c^2$) and $a_*$ is the dimensionless spin parameter). Both periods are sensitive to the mass accretion rate. This can be related to properties of the X-ray spectrum and may tell us the nature of the accretor (neutron star or black hole, Middleton et al. 2019b) without the requirement for pulsations which are themselves transient or may be diluted beyond detection (Bachetti et al., 2014; Mushtukov et al., 2020).

Besides these torques, other effects may warp and precess the disc, e.g tidal torques from the secondary star (although these are unlikely to be strong in cases with known ~ day orbital periods, e.g. Bachetti et al. 2014), magnetic warping (Lai, 2003; Pfeiffer and Lai, 2004), and radiative warping (Pringle, 1996). In the latter case, where super-critical flows – specifically a large scale-height wind and optically thick wind – are present, it is not clear how efficiently the outer disc is irradiated (see also Middleton et al. 2021a). But if the dipole field is strong and/or the accretion rate low, irradiation may well be possible, and similar warps and precession to that seen in Her X-1 is then quite plausible (see Petterson 1977, Kotze and Charles 2012). Regardless of the mechanism driving the precession, it is possible to build a simplified model of the motion of the cone and then numerically extract the lightcurve. Dauser et al. (2017) developed a model which can produce a variety of lightcurves (available within xspec: Arnaud et al. 1999 as ulxlc), under the assumption that the X-ray source lies at the apex of a cone of outflowing, totally opaque, achromatically scattering plasma. Applying this model to the lightcurve of NGC 5907 ULX-1 implies a small total opening angle (~10–15 deg) and substantial beaming (amplification by a factor 60–90) although it is as yet unclear whether this is in conflict with scattering reducing the observed pulse fraction (Mushtukov et al., 2020).

If the neutron star is asymmetric about its rotation axis, free-body precession may also occur on a timescale of (Jones and Andersson 2001):

$$t_{\text{free,prec}} \approx \frac{2.32 \times 10^{-5} P_{\text{NS}}}{\epsilon} \frac{0.5}{1s} \cos \chi \text{ days} \quad (16)$$

where $\chi$ is the angle of the misaligned distortion of the neutron star relative to the rotation axis, and $\epsilon$ is the measure of distortion, related to the surface magnetic field strength through (Lander and Jones, 2009):

$$\epsilon = 2 \times 10^{-11} \left( \frac{B_{\text{surf}}}{10^{12}G} \right)^2 \quad (17)$$

which, notably, Vasilopoulos et al. (2020a) used to estimate the dipole field strength in M51 ULX-7.

It is interesting to consider how the different physical origins for superorbital periods can be distinguished. A natural method would be
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Pulse profiles for the PULXs shown in Table 2, extracted from data taken by *XMM-Newton* (NGC 5907 ULX-1, M51 ULX-7, NGC 7793 P13, NGC 300 ULX-1, NGC 1313 X-1) and *NuSTAR* (M82 X-2 and SMC X-3). Values for the count rates are not provided due to bandpass and detector differences, however, errors are provided for the purposes of comparing the level of constraints within each phase bin.</figcaption></figure>

to consider how the timescales respond to the accretion rate and spin-up/down of the neutron star. Lense-thirring precession (Eqns (14) & (15)) is potentially highly sensitive to the accretion rate (far more so than the spin) and so, where the accretion rate is variable (on long timescales), we would expect variability in the superorbital period (with correlated changes in the X-ray spectrum). Conversely, magnetic and free-body precession is more sensitive to the spin than accretion rate, and the superorbital period should therefore track changes in the neutron star's spin over time.

It is also worth noting that plotting the superorbital versus the orbital period for Galactic systems and those ULXs with known periods (see Table 2, implies a potential correlation (see Townsend and Charles 2020), although care must be taken in interpreting any correlation between logarithmically scaled values which include systems known to be of very different types. It is plausible that the superorbital periods are in some cases responding to tidal interactions, while the mass supply to the compact object may correlate with the binary period, leading to a positive correlation through e.g. Eq. (15). This currently

remains an open question but we can expect progress as the detected number of constrained superorbital periods increases.

In addition to the smooth, periodic/quasi-periodic variations seen in the X-ray and optical bands, there are a variety of changes in the lightcurves of ULXs which can help reveal the nature of the system. Dips are found in NGC 55 ULX-1 (Stobbart et al., 2004); NGC 5408 X-1 (Grisé et al., 2013), NGC 247 ULX-1 (Alston et al., 2021), M51 ULX-7 (a known ULX pulsar, Vasilopoulos et al. 2021), and J125048.6+410743 (Lin et al., 2013). Significantly, these ULXs tend to have relatively soft spectra (M51 ULX-7 is somewhat harder), possibly implying a correlation with sight-lines intercepting regions in the wind subject to Rayleigh Taylor instabilities (e.g. Middleton et al. 2011, 2015a, Takeuchi et al. 2013).

Whilst dips are seen in the lightcurves of a small number of ULXs, eclipses seem to be mostly absent from the currently known ULX population (Middleton and King, 2016). This would make sense if edge-on ULXs emit around the Eddington limit (Poutanen et al., 2007) implying sub-ULX luminosities in many cases (Zhou et al., 2019). However, there
have been reports of eclipses in two ULXs in M51 (Urquhart and Soria, 2016b; Wang et al., 2018). These give constraints on the angle between the companion star's orbit and the X-ray bright inner regions (noting that these do not have to be aligned in a single plane). If there is an estimate of the mass of the donor star, this then constrains the compact object mass (Wang et al., 2018).

Neutron stars probably account for the majority of observed ULXs (King and Lasota, 2016; Middleton and King, 2017; King et al., 2017; Walton et al., 2018c; Wiktorowicz et al., 2019). Such systems are likely to approach spin equilibrium during their lifetimes, where the co-rotation radius approaches the magnetospheric radius. Whilst the effective torque vanishes when the two are close (Dai and Li 2006), should some process push the magnetospheric radius outside the co-rotation radius, a propeller state ensues, where the emission from the accretion column and curtain – but importantly not from radii down to the magnetospheric radius – is switched off. In systems where the magnetospheric radius is large, the result would be a severe drop in flux (although centrifugally driven winds from around the magnetospheric radius would remain). Such propeller states, widely observed in Galactic systems (e.g. Cui 1997, Campana et al. 2008, Tsygankov et al. 2016a) have been invoked to explain the large changes in luminosity seen in M82 X-2 (Tsygankov et al., 2016b). However, in NGC 7793 P13 (Fuerst et al., 2021) and NGC 300 ULX-1 (Vasilopoulos et al., 2019), the spin has continued increasing at the same rate when the source is considerably fainter in the X-rays, implying that some other mechanism – possibly obscuration – drives the change, rather than a propeller state. (In M51 ULX-7, it is as yet unclear whether a similar X-ray off state implies a propeller state or is associated with precession: Vasilopoulos et al. 2021). The predicted large-scale changes in X-ray flux in the propeller state have motivated searches for pulsar candidates among known ULXs. Earnshaw et al. (2018) searched within those ULXs in the 3XMM DR4 catalogue (Rosen et al., 2016), finding five objects which underwent a factor of 10 or more change in flux (with one propeller candidate: M51 ULX-4). Song et al. (2020) compiled XMM-Newton, Swift and Chandra flux lightcurves of 296 ULXs from the (Earnshaw et al., 2019) catalogue, and found that 25 out of 296 ULXs changed flux by more than a factor 10 (based on the available sampling), 17 of which show a bi-modal flux distribution, potentially consistent with sources at spin equilibrium.

### 2.3.4. Spectral-timing studies

As well as standalone studies focussing on either the spectral evolution or timing characteristics of ULXs, there are an increasing number of studies which bring the two together. This spectral-timing analysis can be performed in the time or Fourier domain in a given energy band (e.g. the fractional rms or rms spectrum: Revnivtsev et al. 1999, Gilfanov et al. 2000, and the covariance spectrum: Wilkinson and Uttley 2009), or via lags between bands (via the CCF or Fourier lag spectrum) as a function of energy or frequency (see Uttley et al. 2014 for a comprehensive review). Spectral-timing studies require high signal-to-noise or large amounts of variability (such that the Poisson noise is not dominant), making such analyses challenging for ULXs. However, there have been a number of successful attempts which have permitted otherwise inaccessible insights.

Covariance spectra (Wilkinson and Uttley, 2009) have an advantage over rms spectra, as the Poisson noise between bands is uncorrelated; when the variability is weighted by the linear coherence between bands (see Vaughan and Nowak 1997, Uttley et al. 2014), the constraints on the variability in each energy bin are improved. Middleton et al. (2015a) presented the covariance spectra for a sample of variable ULXs with relatively high data quality, and showed that the variability on relatively short (≤ ks) timescales matched well to the shape of the hard component (seen also in the studies of the covariance spectra of NGC 55 ULX-1: Pinto et al. 2017 and NGC 1313 X-1: Kara et al. 2020). This can be readily explained if the variability emerges from only this component on the timescales being studied. However, it is less

clear why spectrally harder ULXs (where this component is dominant) appear to be less variable on similar timescales (e.g. Heil et al. 2009). One possible explanation is that the variability is extrinsic and being driven by obscuration by optically thick clumps of wind (Middleton et al., 2011), probably created due to fluid instabilities (Takeuchi et al., 2013). Clumps may also then be responsible for dips seen in the softer ULXs (e.g. NGC 55 ULX-1: Stobart et al. 2004 and NGC 247: Alston et al. 2021) which, in the standard paradigm, are viewed at higher inclinations (Poutanen et al., 2007; Middleton et al., 2015a).

The covariance spectrum can also be used to isolate the non-variable component in ULXs. This can be modelled and better compared to theories for emission at soft energies (typically either a thin disc beyond the magnetospheric radius or emission from close to the spherisation radius). This can be complex, depending on the role of variable absorption, but it is at least possible to compare luminosities, and rule out the presence of a thin disc. This was demonstrated in Middleton et al. (2019a) (Fig. 10), and allowed for independent constraints to be placed on the dipole field strength in M51 ULX-7 (see also discussions in Brightman et al. 2018, Christodoulou et al. 2019).

Time (or phase) lags between bands may be extracted in the time domain via the CCF (as shown to great effect when studying a variety of features in XRBs, e.g. lags corresponding to propagation in jets: Gandhi et al. 2010) or via the Fourier cross-spectrum (see Uttley et al. 2014 for a detailed step-by-step guide). The lags are related to the nature and geometry of the accretion flow via the transfer function and are imprinted by the light-travel time through reverberation or transmission (the latter affected by the amount and state of intervening gas). Lags were first identified in a 2006 observation of the ULX, NGC 5408 X-1 (Heil and Vaughan, 2010) between a soft (0.3-1 keV) and hard (1-10 keV) band. This 'soft lag' is ~ a few seconds in NGC 5408 X-1 (De Marco et al., 2013) at frequencies of ~ 10 mHz, is ~ 150 s at 0.1 mHz in NGC 1313 X-1 (Kara et al., 2020), and is a few hundred seconds at around the same frequency in NGC 55 ULX-1 (Pinto et al., 2017). Soft lags have been interpreted as reprocessing in distant material (presumably the wind) and intriguingly appear around the QPO frequencies seen in NGC 5408 X-1 (De Marco et al., 2013). Should the lag be imprinted as a consequence of resonant line scattering or Compton down-scattering, the magnitude can provide an estimate of the optical depth and thickness of the wind. An alternative explanation is that hard photons created in the accretion column are Comptonised by a surrounding envelope (Mushtukov et al., 2017) and a lag is imprinted due to the light travel time across the curtain and the time taken to Compton scatter (Mushtukov et al., 2019). This picture would imply very large or dense curtains in the case of ULXs such as NGC 55 ULX-1.

### 2.4. Multi-wavelength studies

ULXs are defined (albeit empirically) by their X-ray luminosity, although much of the gravitational energy dissipated in the accretion flow may emerge in other wavebands through reprocessing, or as kinetic luminosity in the wind and/or jet. It is therefore crucial to obtain as broad an SED as possible and below we discuss insights obtained in several key bands (note that we discuss observations of bubble nebulae surrounding ULXs separately in Section 2.6.1).

#### 2.4.1. Radio

Radio monitoring campaigns (Körding et al., 2005) and cross-matching of X-ray catalogues with radio surveys (Pérez-Ramírez et al., 2011) have returned limited numbers of ULX radio counterparts, perhaps suggesting that any jets are inherently weak or transient in nature. Targeted observations have tended to focus on the brightest ULXs (HLXs at ~ $10^{41}$ erg/s) as these remain the most promising candidates for X-ray bright hosts of IMBHs other than AGN in dwarf galaxies (Farrell et al., 2009). Their suggested sub-Eddington accretion rates would naturally lead to the expectation of flat spectrum radio emission associated with
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** From Middleton et al. (2019a) showing (left-hand panel) the covariance spectrum of M51 ULX-8 (in green) with the time-averaged Chandra data (in black). The covariance picks out the linearly correlated variability (see Wilkinson and Uttley 2009) and allows the spectrum to be decomposed. In this case, the additional soft component (shown in red in the middle and right-hand panels) is above the Eddington limit for a neutron star and so probably inconsistent with a thin disc. This allowed Middleton et al. (2019a) to place indirect constraints on the dipole magnetic field strength (a higher-order multipole field is not excluded).</figcaption></figure>

a compact jet (if in an analogous hard X-ray state as seen in X-ray binaries and AGN). Placing such HLXs on the radio-X-ray fundamental plane (Merloni et al., 2003; Falcke et al., 2004; Gültekin et al., 2009) opens the possibility of constraining the compact object mass, and this approach has been adopted for a number of ULXs (NGC 4088-X1: Mezcua and Lobanov 2011, NGC2276-3c: Mezcua et al. 2015; N5457-X9: Mezcua et al. 2013a, IC 342 X-1: Cseh et al. 2012).

Confirming the presence of a jet in ULXs can be more difficult given the presence of a radio bright nebula in some cases. These can be up to 100 s of pc across: Mezcua et al. 2013b), potentially analogous to the W50 nebula surrounding SS433 (e.g. Pakull et al. 2006, Pakull and Grisé 2008, Russell et al. 2011, Cseh et al. 2012), and energised by termination shocks between the surrounding gas and outflows. In only a single bright ULX – Ho II X-1 – has the nebula emission been resolved by using VLBI into discrete core and lobe emission (Cseh et al. 2014, 2015 – although there are also twin lobes either side of a bright radio core in the IMBH candidate NGC 2276-3c, see Mezcua et al. 2015). The fading of Ho II X-1's lobes over a period of around 2 years indicates that the jets are sporadic and may have more in common with the discrete ejections in near-Eddington LMXBs, placing a limit on the mass of the compact object in Ho II X-1 of > 25 $M_{\odot}$ (Cseh et al., 2015).

Conversely, jet-inflated nebulae may provide one of the few ways to locate edge-on ULXs. A key example is the radio-optical nebula in M83 (30 times more luminous than W50 at 5 GHz) and collimated jets in NGC 7793 (S26: Pakull et al. 2010, Soria et al. 2010). In both cases, the estimated jet power lies above $10^{40}$ erg/s, but the observed X-ray luminosity to-date is 3–4 orders of magnitude lower, and the sources would therefore not qualify as ULXs. These systems are in many ways similar to SS433 (albeit with even higher jet kinetic luminosities) and provide ample evidence for similar edge-on ULXs being located in nearby galaxies.

As well as jets associated with highly super-Eddington rates of mass transfer (definitely super-Eddington in SS433), transient jet ejections have been observed in some Galactic XRBs moving to less extreme, near-Eddington accretion rates (Fender et al., 2009). These have been successfully detected in nearby galaxies and can be associated with the appearance of a low luminosity ULX (e.g. M31: Middleton et al. 2013). It is plausible that the transient jets detected from HLX-1 (Webb et al.,

2012) suggest the source was moving to near-Eddington rates, which, at a luminosity of > $10^{42}$ erg/s would imply a BH mass of ~10,000 $M_{\odot}$, consistent with the mass estimate from modelling the X-ray spectrum with relativistic disc models (Davis et al., 2011). In future, with the advent of sensitive all-sky radio monitoring (SKA), the detection of such radio events should become commonplace.

Finally, it is intriguing to note that a jet event has also been detected from the transient X-ray pulsar, Swift J0243.6+6124 – which contains a neutron star with $B \approx 10^{11}$ G (King et al., 2017) – when accreting above its Eddington limit (van den Eijnden et al., 2018). At peak radio brightness, the flux density from the jet was < 100$\mu$Jy at 6 GHz (at an assumed distance ≈ 7 kpc – van den Eijnden et al. 2018). Similar jets are therefore unlikely to be discovered in nearby galaxies (although there is the theoretical potential to drive more powerful jets through tapping the neutron star’s angular momentum: Parfrey et al. 2016).

## 2.5. Infra-red to ultraviolet

Generally speaking, observations of accreting systems at energies intermediate to those of X-rays and radio allow us to study the impact of reprocessing, the highest energy populations of electrons in jets and the nature of the companion star. In ULXs, such studies are of vital importance for constraining the true energy budget and age of the system, for understanding the impact on the surroundings and how such high rates of mass transfer actually occur.

### 2.5.1. Infrared

Campaigns studying ULXs in the IR have focused mostly on the nature of the donor star, using a combination of imaging and spectroscopy from ground-based instruments. Initial imaging studies by Heida et al. (2014) identified a number of possible red supergiant (RSG) companion stars (11/62 ULXs in 37 galaxies within 10 Mpc) from a combination of data from WHT/LIRIS, MMT/SWIRC and VLT/ISAAC. Heida et al. (2016) later performed NIR spectroscopy using Keck/MOSFIRE (Multi-Object Spectrometer for Infra-Red Exploration), identifying candidate counterparts for two ULXs (one in NGC 925 and one in NGC 4136) with RSGs. Building on the work of Heida et al. (2016), a large study by López et al. 2017, 2020 made a systematic search for NIR
counterparts to 113 ULXs, with the use of a combination of instruments (imaging via LIRIS/WHT, SOFI/NTT and WIRC on the Palomar Hale 5-m telescope, and spectroscopy with Keck/MOSFIRE). This study yielded candidate counterparts for 38 ULXs, 5 of which are identified as RSGs. The numbers obtained from ground-based studies match well to mid-IR studies performed using *Spitzer*, where Lau et al. (2019) located 12 candidate RSGs in a sample of 96 ULXs within 10 Mpc. These results point towards a highly probable association of RSGs with some ULXs, and an additional – somewhat less considered – form of mass transfer, i.e. that of wind-fed Roche-lobe overflow (e.g. Copperwheat et al. 2007, El Mellah et al. 2019) which has now been incorporated into models for binary population synthesis (Wiktorowicz et al., 2021).

There have been a number of focused campaigns to explore the brightest ULXs, and specifically the PULXs (determining the companion star in the latter being of obvious value when considering the evolution of these systems). Unfortunately, for most PULXs the IR fields are extremely crowded, precluding the identification of a unique counterpart (Heida et al., 2019a). However the companion of the SN impostor, NGC 300 ULX-1 has been identified as an RSG (Heida et al., 2019b).

As well as emission from the companion star, it is quite possible that synchrotron emission from jets and from circumbinary dust could both contribute to observed flux in the IR band. Lau et al. (2019) identify a mid-IR excess and a probable circumbinary disc in several ULXs (including Ho II X-1), and associate variable mid-IR emission from the X-ray spectrally hard ULX, Ho IX X-1, with a jet (as previously suggested by Dudik et al. 2016, Sathyaprakash et al. 2022 in the case of the pulsating ULX, NGC 1313 X-2).

Indirect constraints on the SED which photoionizes gas surrounding the ULX (see Section 2.6) can be extracted through modelling IR lines. Berghea et al. (2010a) identified Ho II X-1's nebula (thought to be inflated via jets) in *Spitzer* observations and, in Berghea et al. (2010b), later performed photoionization modelling to suggest that the lines are consistent with irradiation by the observed soft X-ray and UV emission, requiring little to no beaming (of these components). This is completely consistent with the picture where the soft emission is less beamed than the hard X-rays. The MF16 nebula surrounding NGC 6946 X-1 has also been detected by *Spitzer*, with photo-ionisation modelling by Berghea and Dudik (2012) indicating the possible presence of an O-type supergiant donor star.

### 2.5.2. Optical

Identifying a unique optical counterpart to a ULX potentially gives a powerful probe of the system age via binary evolution arguments. If atomic lines can be isolated and associated with the companion, then subsequent studies may allow the mass of the compact object itself to be constrained by standard radial velocity techniques. There are a number of issues which make such investigations challenging. Given the distance to many ULXs, there may be several optical counterparts associated with the X-ray astrometric position (e.g. Gladstone et al. 2013). The distance can make lower mass companion stars undetectable, perhaps misleadingly favouring an association with any high mass star within the neighbourhood of the ULX. Further, the photometric signal of the star may be contaminated by direct emission from the wind photosphere depending on the accretion rate (Poutanen et al., 2007), and emission from the uncovered thin disc which may be increased by irradiation (should the geometry, irradiation pattern and wind optical depth permit this – Sutton et al. 2014, Yao and Feng 2019). Finally the spectroscopic signatures (e.g. He II emission lines) may also be contaminated or dominated by the wind (Fabrika et al., 2015).

Despite these challenges, there have been a number of attempts to find and study the optical counterparts of ULXs (besides the results discussed in the IR above). Many studies (too many to recount in this review) over the last 20 years have targeted individual ULXs (see e.g. Roberts et al. 2001), whilst others have aimed to describe the broader population. As important examples of the latter, the sensitivity

and angular resolution of HST allowed Roberts et al. (2008) to locate counterparts to four out of six ULXs, whilst Tao et al. (2011) studied the optical counterparts of 13 ULXs, finding the optical flux and colour to vary on the timescale of days to years, ruling out changes in extinction as the origin. Gladstone et al. (2013) utilised HST photometry to study 33 ULXs with *Chandra* positions, locating 13 ± 5 counterparts and finding O-star companions to be rare.

Where unique counterparts have been identified, there have been attempts to constrain the mass and species of compact object via dynamical arguments (most using the broadened He II λ4686 emission line):

* Using *Gemini* observations, Roberts et al. (2011) attempted to constrain the mass of the compact object powering NGC 1313 X-2 (which is now known to be a neutron star: Sathyaprakash et al. 2019b) but were unable to do so due to the absence of sinusoidal motion in the line velocities.

* Cseh et al. (2013) were unable to locate an orbital period, but placed a limit on the mass (<510$M_\odot$) of the compact object in NGC 5408 X-1 from a limit on the semi-amplitude in the absence of periodic modulations, and under the assumption that the He II line is produced in the irradiated outer accretion disc.

* Liu et al. (2013) used the He II λ4686 line to place a lower limit on the mass of the compact object in M101 ULX-1 of $5M_\odot$.

* Motch et al. (2014) placed an upper limit on the mass of NGC 7793 P13 (<15 $M_\odot$) by modelling the ~ 64 d period in the V and U band lightcurves as irradiation of the well-identified B9Ia companion star; this result is of course consistent with the discovery that a neutron star powers this ULX (Fürst et al., 2016; Israel et al., 2017b).

Use of the He II λ4686 emission line to determine the compact object mass via dynamical means is problematic, as Fabrika et al. (2015) show that this line in ULXs probably originates in the outer regions of an outflowing wind, as in SS433. There remains a question about the width of the He II lines in NGC 7793 P13 – these cannot originate from the companion star (Fabrika et al., 2015), but also appear to be modulated on the same 64 d timescale. This may point towards the wind precessing, as it does in SS433 (see Section 2.3.3 for a discussion of super-orbital periods).

### 2.5.3. UV

Models which consider the temperature profile of the super-critical inflow (Poutanen et al., 2007) already imply that, for high mass accretion rates, the outer photosphere of the wind should emit in the UV around the Eddington limit for the compact object. Compton downscattering of photons from the innermost regions could also produce extreme luminosities in the optical/UV band. This may explain the inference that SS433 emits at over $10^{40}$ erg/s in the UV (Dolan et al., 1997; Waisberg et al., 2019).

UV observations of ULXs are limited by the difficulties inherent in observing in much of this band, but a particularly notable result is that of Kaaret et al. (2010), who observed NGC 6946 X-1 with HST, finding a spectrally soft ULX showing strong variability and QPOs, surrounded by a photoionized and shock-heated optical nebula (MF16 - see below). Spectral fitting carried out by these authors implied a UV luminosity in excess of $10^{39}$ erg/s, consistent with the luminosity needed to produce the observed He II λ4686 emission-line luminosity (Abolmasov et al., 2008).

More recently, monitoring of ULXs with Swift has started to provide hints of how the X-rays (via the XRT) and UV (via UVOT) may be correlated (or anti-correlated). As shown in Fuerst et al. (2021), the ~ 64 d optical period observed in NGC 7793 P13 (Motch et al., 2014), is also seen in the X-rays (with a precise value of 65.21 ± 0.15 d, see also Hu et al. 2017), in the UV (with a period of 63.75$^{+0.17}_{-0.12}$ d) and in the fits to the changing pulse frequency, yielding an orbital period of 64.86 ± 0.19 d. It is clear from Fig. 11 that the UV brightness is not in phase with the X-rays (and indeed, the UV modulation is considerably stronger when the X-rays are dim).
<figure><img src="image_13.png" /><figcaption>Fig. 11. From Fuerst et al. (2021) showing the long term Swift X-ray (XRT) and UV (UVOT, U band) lightcurve of the ULX pulsar, NGC 7793 P13, with modulations around 65 days. The drop in X-ray flux in 2019 does not appear to result in a substantial decrease in the spin-up rate as indicated by the lower panels (see Fuerst et al. 2021 for details). It is also intriguing that the UV and X-ray modulations appear to be – to some extent – anti-correlated.</figcaption></figure>

## 2.6. The environment of ULXs

ULXs are a product of their host environment, and interact with it through kinetic and radiative feedback. In this section we discuss the small and larger scale environments of ULXs and the insights these provide.

### 2.6.1. ULX nebulae

As first reported by Pakull and Mirioni (2002), a substantial number of ULXs are observed to be surrounded by a bubble of warm gas (although the phenomenology is in reality somewhat varied), visible in the optical via emission lines (and in some cases detected in the radio as well). As indicated in Table 1, these ‘ULX nebulae’ are typically ~10-

100 s pc across (the largest presently known being ~500 pc in diameter) and tend to be considerably larger than typical SNRs (e.g. Asvarov 2014) although there is sometimes disagreement on the actual size (e.g. Moon et al. 2011 versus Kaaret et al. 2004). Bearing this point in mind, the dimensions provided in Table 1 do not necessarily account for the entire region influenced by the ULX (e.g. there is a reported 800 pc zone of weakly ionised gas surrounding NGC 1313 X-1: Pakull and Mirioni 2002).

Analogous to the W50 nebula of SS433, we expect gas from the original SNR to be inflated by jets and winds from the inner accretion flow (e.g. Feng and Kaaret 2008, Pakull and Grisé 2008) and photoionized by the UV-X-ray radiation field (which we broadly expect to be anisotropic). Indeed, ULX nebulae show a mixture of photoionized
<caption>Table 1</caption>

The properties of ULX nebulae, with non-ULX (i.e. hidden ULX) systems included at the bottom in a detached table. Objects and values from Pakull and Mirioni (2002) unless indicated otherwise. Dimension is diameter unless additional dimension is given. Corresponding references are: a: Ramsey et al. (2006); b: Cseh et al. (2012), c: Roberts et al. (2003); d: Kaaret et al. (2004); e: Miller (1995); f: Moon et al. (2011); g: Pakull and Grisé (2008); h: Grisé et al. (2012); i: Soria et al. (2006); j: Soria et al. (2021); k: Abolmasov et al. (2007a); l: Dopita et al. (2012) m: Pakull et al. (2010); n: Fabrika (2004); o: Soria et al. (2014).

<table><thead><tr><th>ULX</th><th>Dimension(s)</th><th>Inferred age (years)</th><th>Radio/IR (jet) counterpart?</th><th>L<sub>kin</sub> (10<sup>39</sup> erg/s)</th></tr></thead><tbody><tr><td>NGC 1313 X-1</td><td>240 pc</td><td></td><td></td><td></td></tr><tr><td>NGC 1313 X-2</td><td>350 × 500 pc<sup>a</sup></td><td>~1 × 10<sup>6</sup> (b)</td><td></td><td>~10</td></tr><tr><td>IC 342 X-1</td><td>110 pc<sup>c</sup></td><td>~6 × 10<sup>5</sup></td><td>Y</td><td>~0.7</td></tr><tr><td>Ho II X-1</td><td>45 pc<sup>d</sup></td><td></td><td>Y</td><td></td></tr><tr><td>Ho IX X-1</td><td>250 pc<sup>e</sup></td><td>~1 × 10<sup>6</sup></td><td>Y</td><td>~10</td></tr><tr><td>M81 X-6</td><td>115 × 42 pc<sup>f</sup></td><td></td><td></td><td></td></tr><tr><td>IC 2574 X-1</td><td></td><td></td><td></td><td></td></tr><tr><td>NGC 2403 X-1</td><td>300 pc<sup>g</sup></td><td></td><td></td><td></td></tr><tr><td>NGC 4559 X-7</td><td></td><td></td><td></td><td></td></tr><tr><td>NGC 4631 H7</td><td></td><td></td><td></td><td></td></tr><tr><td>NGC 4861 X-1</td><td></td><td></td><td></td><td></td></tr><tr><td>NGC 4861 X-2</td><td></td><td></td><td></td><td></td></tr><tr><td>NGC 5204 X-1</td><td>360 pc</td><td></td><td></td><td></td></tr><tr><td>NGC 5408 X-1</td><td>60 pc<sup>h</sup></td><td>~1 × 10<sup>5</sup> (i)</td><td>Y</td><td>0.7 (i)</td></tr><tr><td>NGC 5885 ULX</td><td>200 × 300 pc<sup>g</sup></td><td>6 × 10<sup>5</sup> (j)</td><td>Y</td><td>~20 (j)</td></tr><tr><td>NGC 6946 X-1</td><td>20 × 34 pc<sup>k</sup></td><td></td><td>Y</td><td></td></tr><tr><td>NGC 7793 S26</td><td>185 × 350 pc (l)</td><td>4 × 10<sup>5</sup> (l)</td><td>Y</td><td>~10–100 (l,m)</td></tr><tr><td>SS433 (W50)</td><td>50–70 pc<sup>a</sup></td><td>~2 × 10<sup>5</sup> (b)</td><td>Y</td><td>>0.1</td></tr><tr><td>M83 MQ1</td><td>13 pc<sup>a</sup></td><td>1.3–2.0 × 10<sup>4</sup> (o)</td><td>Y</td><td>~30 (n)</td></tr></tbody></table>

(notably He II 4686λ) and shock ionized lines ([OI] 6300λ and [OIII] 5007λ) – some nebulae appearing to be dominated by one process over the other (although both probably occur in any given source e.g. Abolmasov et al. 2007a). ULX nebulae have conclusively been shown to not be inflated by O star associations (i.e. a superbubble) and are unlikely to be due to a hypernova explosion (Pakull and Grisé, 2008).

Applying a model where the ULX nebula is inflated via the mechanical action of outflows (i.e. a shock dominated rather than photo-ionised dominated nebula) gives interesting estimates for the age of the system. This in turn provides clues as to how the binary may have evolved, and for how long super-critical accretion can be maintained. From Weaver et al. (1977) the radius of the bubble is (compare Eq. (117))

$$R = \left(\frac{125}{154\pi}\right)^{1/5} \left(\frac{L_{\text{tot}}}{\rho_0}\right)^{1/5} t^{3/5} \qquad (18)$$

where R is the radius of the bubble at time t, $L_{tot}$ is the mechanical luminosity of the wind and/or jet, and the ISM density (assumed to be uniform) is $\rho_0 = \mu m_p n$, where $\mu$ is the mean atomic weight of the gas, $m_p$ is the proton mass and $n$ is the hydrogen number density. The characteristic age of the bubble is then given by $\tau = 3R/5v_{\text{exp}}$ where $v_{\text{exp}}$ is expansion velocity. These equations allowed Cseh et al. (2012) to determine the energy content and age of several ULX nebulae, the results of which are shown in Table 1.

Assuming $L_k = 1/2 \dot{m}_{\text{wind}} v_{\text{wind}}^2$, and that the wind provides the entire energy reservoir of the bubble (i.e. provides $L_{tot}$), then equation (18) may alternatively be written in terms of the wind properties (e.g. Pinto et al. 2020a):

$$R = 0.76 \left( \frac{\dot{m}_{\text{wind}} v_{\text{wind}}^2}{2\rho_0} \right)^{1/5} t^{3/5} \qquad (19)$$

where $\dot{m}_{\text{wind}}$ is the mass loss rate in the wind and $v_{\text{wind}}$ is the wind velocity. For measured wind velocities in ULXs (e.g. Pinto et al. 2016, Kosec et al. 2018b), a kinetic luminosity of a few × 10³⁹ erg/s is plausible which gives bubble radii ~100 pc, consistent with observation (Table 1). We note however that (a) the kinetic energy of the impacting wind ~ $L_{Edd}$, and (b) the radius depends on the driving luminosity (radiative and mechanical) to the power 1/5, so that even a difference of a factor 100 alters the radius only by a factor 2.5.

The He II 4686λ recombination line is more easily created via photoionization than by shock-excitation (e.g. Berghea et al. 2010b), so

using this line allows us to treat the nebula as a UV-X-ray bolometer, and infer the irradiating luminosity by using photoionization codes such as CLOUDY (Ferland et al., 1998). This has been performed for a number of ULXs, notably Ho II X-1 (Pakull and Mirioni 2002, Berghea et al. 2010b who used the OIV line in the IR band) and NGC 5408 X-1 (Kaaret and Corbel, 2009), with the finding that the required luminosities are a good match to the X-ray luminosities observed from the source. Such consistency is used as an argument against strong geometrical beaming in these ULXs as the nebula should see the isotropic emission (Pakull and Mirioni, 2002). However, the ULX wind photosphere emits a luminosity ~ $L_{Edd}$ which – depending on the mass transfer rate – can have a temperature ~ 10⁶ K, implying strong near-isotropic soft X-rays and EUV (see the discussion of Eqs. (110), (111) below). In agreement with this, both Ho II X-1 and NGC 5408 X-1 are spectrally soft and potentially viewed at larger inclinations than some of the spectrally harder ULXs. Should this be true then it is not surprising that the irradiating luminosity inferred from photoionisation modeling is consistent with the observed luminosity, as we would naturally expect the more obscured (at such angles) hard X-ray emission to be the most geometrically beamed (rather than the more isotropic soft emission).

In principle, the number of ULX bubble nebulae can also be used to explore the role collimation/beaming plays in the wider population based on the reasoning that, for every ULX observed with a given beaming factor, there should be many more which we do not observe as a ULX. As the nebula’s photoionisation responds mostly to unbeamed soft X-rays and EUV photons, such a premise would hold only if energy was also imparted through shocks, which would indeed seem to be the case: (Abolmasov et al., 2007a). This then implies we should detect bubble nebulae with X-ray faint (non-ULX) counterparts; indeed there is already support for such observations both in other galaxies (e.g. Soria et al. 2014) and within our own (SS433: see Section 2.8.3). Searches for X-ray quiet bubble nebulae are made somewhat harder, as the ULX nebulae with strong photoionisation signatures (NGC 5408 X-1 and Ho II X-1) appear to be smaller (10s of pc rather than 100s of pc: Table 1), precluding detection of such relatively small structures in the absence of high angular resolution, narrow filter instruments (Russell et al., 2011). However, Pakull and Grisé (2008) searched VLT images of NGC 1313 for additional bubble nebulae without ULX counterparts and found only one possible candidate. As shown by Wiktorowicz et al. (2019), should geometric beaming be strongly related to accretion
<caption>**Table 2**</caption>

Key parameters of interest for the ULX pulsars with values taken from those papers indicated by the reference attached to the ULX name. The period over which the spin up rate is calculated is given in parentheses as is the energy range for the measured pulse fraction. a: *Bachetti et al.* (2014), b: *Bachetti et al.* (2020). It is worth noting that, in the case of M82 X-2, there are periods of both spin up and spin down – see *Bachetti et al.* (2020) for details. c: *Brightman et al.* (2019); d: *Israel et al.* (2017a); e: *Walton et al.* (2016a); f: *Fürst et al.* (2016); g: *Israel et al.* (2017b); h: *Fuerst et al.* (2021); i: *Carpano et al.* (2018); j: *Heida et al.* (2019b), k: *Edge et al.* (2004); l: *Townsend et al.* (2017); m: *Tsygankov et al.* (2017); n: *Corbet et al.* (2004); o: *Rodríguez Castillo et al.* (2020); p: *Brightman et al.* (2020a); q: *Sathyaprakash et al.* (2019b); r: *Trudolyubov et al.* (2007); s: *Doroshenko et al.* (2018); t: *Vasilopoulos et al.* (2020b); u: *Brightman et al.* (2018); v: *Trudolyubov* (2008).

<table><thead><tr><th>ULX</th><th>P<sub>pulse</sub> (s)</th><th>Spin up rate (ss<sup>-1</sup>)<br>(baseline)</th><th>Pulse fraction (E range)</th><th>P<sub>orb</sub></th><th>P<sub>s-orb</sub> (d)</th><th>L<sub>max</sub> (×10<sup>39</sup> erg s<sup>-1</sup>)</th><th>M<sub>2</sub> (M<sub>☉</sub>)</th></tr></thead><tbody><tr><td>M82 X-2<sup>a,b,c</sup></td><td>1.4</td><td>1 × 10<sup>-11</sup> (2.5 years)</td><td>5%–13% (3–30 keV)<br>8%–23% (10–30 keV)</td><td>2.5</td><td>~60</td><td>20</td><td>≥5.2</td></tr><tr><td>NGC 5907 ULX-1<sup>d,e</sup></td><td>1.1</td><td>~8 × 10<sup>-10</sup> (11 years)</td><td>12% (0.2–2.5 keV)<br>20% (7–30 keV)</td><td>5.3</td><td>78</td><td>200</td><td>?</td></tr><tr><td>NGC 7793 P13<sup>f,g,h</sup></td><td>0.4</td><td>4 × 10<sup>-11</sup> (4 years)</td><td>18%–22% (0.1–12 keV)</td><td>64</td><td>64</td><td>~10</td><td>18–23 (B9I)</td></tr><tr><td>NGC 300 ULX-1<sup>i,j</sup></td><td>31.6</td><td>~6 × 10<sup>-7</sup> (3.6 days)</td><td>~60% (0.2–10 keV)<br>~70% (3–20 keV)</td><td>>290</td><td>?</td><td>5</td><td>8–10 (RSB?)</td></tr><tr><td>SMC X-3<sup>k,l,m,n</sup></td><td>7.8</td><td>~4 × 10<sup>-9</sup> (110 days)</td><td>30%–100% (3–70 keV)</td><td>~45</td><td>?</td><td>~3</td><td>>3.7 (Be?)</td></tr><tr><td>M51 ULX-7<sup>o,p</sup></td><td>2.8</td><td>~1 × 10<sup>-9</sup> (13 years)</td><td>&lt;5%–20% (0.1–12 keV)</td><td>2</td><td>38</td><td>~10</td><td>>8</td></tr><tr><td>NGC 1313 X-2<sup>q</sup></td><td>1.5</td><td>1.2 × 10<sup>-12</sup> (98 days)</td><td>~5% (0.3–10 keV)</td><td>&lt;4</td><td>?</td><td>20</td><td>?</td></tr><tr><td>NGC 2403 ULX<sup>r</sup></td><td>18</td><td>3.4 × 10<sup>-10</sup> (?)</td><td>?</td><td>60–100</td><td>?</td><td>1.2</td><td>Be?</td></tr><tr><td>Swift J0243.6+6124<sup>s</sup></td><td>9.86</td><td>2.2 × 10<sup>-10</sup> (?)</td><td>?</td><td>28.3</td><td>?</td><td>≥1.5</td><td>Be?</td></tr><tr><td>RXJ0209.6-7427<sup>t</sup></td><td>9.3</td><td>1.165 × 10<sup>-10</sup> (?)</td><td>?</td><td>> 50</td><td>?</td><td>1–2</td><td>Be</td></tr><tr><td>M51 ULX-8<sup>u</sup></td><td>?</td><td>?</td><td>?</td><td>8–400</td><td>?</td><td>2</td><td>40(?)</td></tr><tr><td>NGC 1313 PULX<sup>v</sup></td><td>765.6</td><td>?</td><td>38% (0.3–7 keV)</td><td></td><td>?</td><td>1.6</td><td>Be?</td></tr></tbody></table>

rate (King, 2009), then binary population synthesis would imply that (in the absence of precession) the difference between the observed and underlying ULX population should be a factor of 5–15 unless black holes are the dominant population (in which case the observed is very similar to the underlying population). This would imply ~10-30 ULXs should be present in NGC 1313 (assuming neutron stars dominate the underlying population), around 1/3 of which would have large bubble nebulae (Pakull and Grisé, 2008). This is not significantly inconsistent with observation and any tension is relaxed if the population demographic is not 100% neutron stars or if the beaming is not as extreme a function of accretion rate as assumed

A number of ULXs have radio nebulae (see Table 1) of a similar size to their optical nebulae (Lang et al., 2007; Cseh et al., 2012). In Ho II X-1 (which shows a clear photo-ionisation dominated optical nebula), the radio nebula has been resolved into a discrete lobe-core-lobe structure (Cseh et al., 2014). The presence and variability of the radio structure in Ho II X-1 implies that repeat ejections also energise the surrounding gas (Cseh et al., 2015) in a similar manner to SS433 (W50), the inflated bubble nebula in M83 S2 (Soria et al., 2020), M82 Q1 (Soria et al., 2014) and the collimated jet source in NGC 7793 (S26: Soria et al., 2010).

### 2.6.2. ULX host environment

There is a trend favouring increased numbers of ULXs with lower host metallicity (e.g. Mapelli et al. 2010) and higher host star formation rate (see Section 2.1). This makes it important to explore the sub-galactic environment where ULXs are located, as this can provide powerful constraints on formation.

Besides the ionised environment (which can be extremely large, e.g. the 800 pc region around NGC 1313 X-1: Pakull and Mirioni 2002), it has been widely observed that ULX counterparts tend to be coincident with young stellar clusters or OB associations (e.g. Zezas et al. 2002, Goad et al. 2002, Gao et al. 2003, Grisé et al. 2006, 2008, 2012, Swartz et al. 2009, Liu et al. 2007, Abolmasov et al. 2007b). Although ionizing, the stellar environment provides insufficient flux to explain the emission line luminosities seen from the bubble nebulae (and the expected kinetic input due to supernovae is also not enough to inflate them: Ramsey et al. 2006) – indeed, in the case of IC 342 X-1, there are no O stars within 300 pc of the ULX (Feng and Kaaret, 2008).

ULXs are also found (albeit rarely) in globular clusters. Relatively few (17 in total) are known at the time of writing: three in NGC 1399 (Shih et al., 2010; Irwin et al., 2010; Dage et al., 2019), two in NGC 4649 (Roberts et al., 2012; Dage et al., 2019), five in NGC 4472 (Maccarone et al., 2007, 2011) and a further seven in M87 (Dage et al., 2020).

ULXs located in globular clusters are expected to differ markedly from those located in the wider population. Given the dense environment of a globular cluster, such ULXs likely have a dynamical origin Ivanova et al. (2010) – distinctly unlike those ULXs located in star forming regions (e.g. Brightman et al. 2020b), and the companion stars are also considerably older (and therefore less massive, e.g. the ULX may be feeding from a white dwarf: Steele et al. 2014, potentially in a short period, ultracompact binary King 2011). In terms of brightness, the typical X-ray luminosity is ~10<sup>39</sup> erg s<sup>-1</sup> with the only globular cluster ULX found to exceed 4 × 10<sup>39</sup> erg s<sup>-1</sup> being GCU8 in NGC 1399 (Dage et al., 2019).

### 2.7. ‘Hyperluminous’ X-ray sources?

The term ‘hyperluminous’ has sometimes been used to refer to sources with inferred (isotropic) luminosities in excess of 10<sup>41</sup> erg s<sup>-1</sup>. However there is little evidence to support the implicit suggestion that these bright ULXs form a physically distinct class. They certainly do not need to contain unusually high mass black holes or IMBH: the ULX pulsar NGC 5907 ULX-1 reaches such luminosities (e.g. Israel et al. 2017a) and evidently contains a neutron star. The well known extreme source SS433 is probably a ULX viewed ‘from the side’ (i.e. from outside the X-ray beam: see Section 3.14.1 and inferred on evolutionary grounds (see King et al. 2000) to be supplied with mass at a rate ~10<sup>-5</sup> M<sub>☉</sub> yr<sup>-1</sup>. From the beaming formula (55), if viewed from within this beam, SS433 would have an inferred isotropic luminosity L<sub>sph</sub> ~10<sup>41</sup> erg s<sup>-1</sup> if the accretor is a 10 M<sub>☉</sub> black hole. We discuss the source HLX-1 in Section 3.16.2 below.

### 2.8. Galactic ULXs

A number of Galactic sources are known to reach or exceed their Eddington limit for extended periods of time (i.e. discounting type I bursts etc.). Their proximity allows for more detailed tests of how super-Eddington mass supply operates. Here we report those sources which would qualify as ULXs, with observed (or inferred) X-ray luminosities ≥ 10<sup>39</sup> erg s<sup>-1</sup>.

#### 2.8.1. Low mass X-ray binaries

Several low mass X-ray binaries (LMXBs) containing black holes are known to exceed their Eddington limit during outbursts mediated by the classical disc instability (Lasota, 2001), recently modified to account for large discs powering super-Eddington outbursts (Hameury and Lasota, 2020). GRS 1915+105 reaches rates close to its Eddington
limit (and has been in outburst for at least 30 years: Castro-Tirado et al. 1992), GRO J1655-40 may have entered a super-Eddington state in its 2005 outburst (Neilsen et al., 2016), V404 Cygni is thought to have exceeded its Eddington limit in its recent 2015 outburst (Motta et al. 2017, Miller-Jones et al. 2019, and potentially in the earlier outburst if not all of the radiation was visible to us) in a similar fashion to V4641 Sgr in 1999 (Revnivtsev et al., 2002). In addition to these confirmed black hole systems, neutron star LMXBs certainly exceed their own Eddington limit and may even approach the empirically defined ULX threshold of ~ $10^{39}$ erg/s (for instance the bursting pulsar, GRO 1744-28: Sazonov et al. 1997). It is notable that all of these systems have fairly long orbital periods (2.6–31 d), which results in a large outer disc radius and the creation of a large reservoir of material for accretion (see Hameury and Lasota 2020 for details).

In addition to large X-ray and optical luminosities, and powerful mass-loaded outflows, optically thin, flaring jet emission is observed to coincide with outbursts of those sources approaching their Eddington limits (Fender et al., 1997). Such radio-bright, highly variable events provide another way of detecting extra-galactic systems entering a near-Eddington/super-Eddington state (and may favour the detection of black hole over neutron star systems: Middleton et al. 2013, van den Eijnden et al. 2018).

### 2.8.2. High mass X-ray binaries

A number of accretion powered pulsars fuelled by high mass counterparts (i.e. high mass X-ray binaries: HMXBs) are known to reach ~ $10^{39}$ erg/s, including the X-ray millisecond pulsar AO538-66 (Skinner et al., 1982), SMC X-1 (e.g. Bonnet-Bidaud and van der Klis 1981) and SMC X-3 (mentioned already in this review).

A recent important observation of Galactic HMXBs in the context of ULXs is that of Swift J0243.6+6124 (Kennea et al., 2017). This magnetised neutron star (with an inferred dipole field of ≥ $10^{12}$G: Doroshenko et al. 2018) accretes from the decretion disc of an O9.5Ve star (Reig et al., 2020), reaches a (0.5–10 keV) luminosity of $10^{39}$ erg s⁻¹, shows pulsations on a timescale of ~ 10s (Jenke and Wilson-Hodge, 2017), shows evidence for fast outflows (van den Eijnden et al., 2019; Tao et al., 2019), and produces radio jets (van den Eijnden et al., 2018). These jets are relatively weak radio emitters, with a maximum flux density of 100 µJy at 6 GHz (at an assumed distance of around 7 kpc: van den Eijnden et al. 2018), and are therefore unlikely to be detected in other galaxies in the absence of substantial Doppler boosting. However, detections of jet events from Galactic HMXBs may become more widespread with the advent of sensitive all sky radio monitoring (SKA) and will permit tests of neutron star jet models (e.g. Parfrey et al. 2016).

### 2.8.3. SS433

Discovered by Stephenson and Sanduleak (1977), the A type supergiant companion star in SS433 (Gies et al., 2002) is transferring mass to the compact object at a rate ~ $10^{-5}$ M☉ yr⁻¹ King et al. (2000). This would imply $\dot{m} \sim 100$ and 1000 for a 10 M☉ black hole or 1.4 M☉ neutron star respectively. This strongly super-critical mass transfer rate suggests that SS433 is ‘intrinsically’ a Galactic ULX similar to those with luminosities ≥ $10^{40}$ erg s⁻¹ seen in other galaxies (Begelman et al., 2006). However it does not qualify as a ULX on the grounds of observed luminosity – see Section 3.14.1.

One of the most dramatic features of SS433 is the baryon-loaded jet (Fabian and Rees, 1979; Margon et al., 1979; Kotani et al., 1994; Marshall et al., 2013) which is ejected at ~0.26c and is twisted into a helical shape by the 162 day precession of the surrounding disc and wind (driven by a slaved disc – van den Heuvel et al. 1980).

The baryon-loading, low velocity ($v/c \approx 0.25$) and 162.4 d precession period of the jet probably all result because it collides with the precessing outer disc (Begelman et al. 2006. The most like cause of the precession is radiation-driving (Pringle, 1997: MN 292, 136).

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** From Pinto et al. (2020b), showing simulations which highlight the spectral quality of observations of ULXs (in this case NGC 1313 X-1) which will be made possible by the launch of *Athena* with its X-ray integral field unit (X-IFU). The simulations utilise absorption and emission models derived from existing observations (see Pinto et al. 2020b for more details).</figcaption></figure>

This picture also explains why the 162.4d jet precession has a small modulation at the 13d orbital period.

We limit ourselves here to a discussion of the similarities and differences between SS433 and other ULXs, and direct the reader to Fabrika (2004) for a comprehensive review.

SS433 can be regarded as a key example of a ULX in our own Galaxy, merely inclined such that we do not see directly into the wind-cone (Begelman et al., 2006; Middleton et al., 2021b). As indicated in the previous sections, besides the similarities of the broad He II recombination lines (Fabrika et al., 2015) to those seen in other ULXs, there are a number of sources which look distinctly like SS433, viewed at similar edge-on inclinations. Notably the inflated nebula in M83 and collimated jet source in NGC 7793 (Pakull et al., 2010; Soria et al., 2010; Dopita et al., 2012). In both of these cases, the jet power appears to be greater than in SS433 for similar or slightly higher apparent X-ray luminosities. Indeed, very many ULXs must also be oriented in such a manner such that we only observe a relative proportion of the entire ULX population (see Middleton and King 2017, Wiktorowicz et al. 2019). Although it is rather hard to study the innermost regions of SS433, it provides a unique opportunity to explore how such systems are fed, its feedback into the local ISM, the nature of the binary (e.g. probes of circumbinary structure: Waisberg et al. 2019), and the structure of the X-ray obscuring regions (Middleton et al., 2021b).

### 2.8.4. Future ULX observations

It is apparent that the study of ULXs has made dramatic progress as instrumentation has improved (e.g. the discovery of pulsations, winds and jets). In the near future, a number of developments will see the field evolve in new, exciting ways.
## 2.9. X-rays

Much of our understanding of ULXs derives from the X-ray band (even if this only actually grants us a restricted view of the full energetics of the system). Several new X-ray satellites are anticipated in the next two decades, including the Astro-H recovery mission, XRISM (dedicated to high energy-resolution X-ray imaging: XRISM Science Team 2020), and ESA's flagship Athena mission (scheduled to launch in the 2030s; Nandra et al. 2013), comprising a wide field imager (WFI) and X-ray integral field unit (X-IFU). Both instruments onboard Athena represent a significant improvement in effective area for their respective specialisms: the WFI is around an order of magnitude more sensitive than XMM-Newton's EPIC pn detector (from ~0.5–2 keV), and the X-IFU is between 45 (at 1 keV) and 6 (at 7 keV) times more sensitive even than XRISM. This will allow considerably deeper searches for pulsations in known ULXs, placing important constraints on the underlying population demographic. The sensitivity will also allow new searches for QPOs and reveal new details in the lag spectra (Heil and Vaughan, 2010; De Marco et al., 2013; Kara et al., 2020; Middleton et al., 2021b). Finally, the impact on studying ULX outflows and how they vary with spectral state will be dramatic, allowing results previously obtained in ~day exposures to be obtained in only a couple of hours (see Fig. 12 which shows a simulation of the atomic lines detected in NGC 1313 X-1 as seen by Athena: Pinto et al. 2020b).

## 2.10. Multi-wavelength

As we have discussed in the preceding sections, multiwavelength observations of ULXs have been highly revealing, allowing the discovery of jets, studies of optical nebulae and likely donor stars, and providing insights into the nature of the accretion flow (irradiation and emission from the outer photosphere).

The future of optical observing of ULXs – as with other variable astronomical objects of interest – is poised to change dramatically with the introduction of all-sky, high angular resolution, high-throughput observing via the Vera Rubin LSST, which is due to start operations in 2023. LSST will provide highly sensitive snapshots (e.g. a limiting r-band magnitude of 24.7 in a single 30 s exposure at SNR = 5) of ~18,000 deg² of the Southern sky (DEC = -65 → +5). Combined with a median seeing-limited angular resolution of 0.7", this will allow deeper probes of possible optical counterparts and explore optical variability for signs of irradiation or precession (e.g. Motch et al. 2014, Middleton et al. 2019b). LSST is not equipped with narrow filters (nor any spectroscopic instrumentation), precluding locating new bubble nebulae, but the MUSE IFU on the VLT is set to transform our understanding of ULX feedback by spatially resolving different regions of ULX bubble nebulae. Finally, the introduction of the ELT (sometime after 2027), will also provide a step change in signal-to-noise, allowing deep spectroscopic studies of ULX counterparts to be performed. By identifying lines from the companion star, mass functions will no doubt follow.

At longer wavelengths, there will be a great deal of interest in the JWST view of ULXs, as this will provide opportunities to locate further red super-giant companions (e.g. Heida et al. 2016), study the surrounding nebulae (e.g. using the OIV line: Berghea et al. 2010b) at a resolution even higher than can be achieved by MUSE, search directly for IR bright ejecta or indirectly search for the presence of jets via IR variability (e.g. Dudik et al. 2016).

At the longest wavelengths, SKA will transform our view of the variable radio sky. In relation to ULXs, the enormous sensitivity of SKA will allow us to locate LMXB ULXs entering outburst in nearby galaxies (Middleton et al., 2013), search for variability in radio bubble nebulae (i.e. unresolved jets), and potentially resolve jet structures (Cseh et al., 2014, 2015). In terms of the impact on studies of Galactic ULXs, SKA will allow the deepest searches for jet ejections associated with accreting pulsars (van den Eijnden et al., 2018) and allow models for their production to be tested (e.g. Parfrey et al. 2016).

## 3. Theory

### 3.1. Definitions and notation

The Eddington luminosity is

$$L_{\text{Edd}} = \frac{4\pi G M m_{\rho} c}{\sigma_{T}} = 1.26 \times 10^{38} \left(\frac{M}{M_{\odot}}\right) \text{erg s}^{-1}, \quad (20)$$

where $m_\rho$ is the proton mass, $\sigma_T$ the Thomson scattering cross-section. In this review $L_{Edd}$ defined by Eq. (20) is a unit of luminosity depending on accretor's mass only.

The Eddington accretion rate is defined as

$$\dot{M}_{\text{Edd}} = 1.4 \times 10^{18} \eta_{0.1}^{-1} \left(\frac{M}{M_{\odot}}\right) \text{g s}^{-1} \quad (21)$$

$$= 2.2 \times 10^{-8} \eta_{0.1}^{-1} \left(\frac{M}{M_{\odot}}\right) M_{\odot} \text{yr}^{-1}, \quad (22)$$

where $\eta = 0.1\eta_{0.1}$ is the radiative efficiency of accretion, $\kappa_{es}$ the electron scattering (Thomson) opacity. We will often use accretion rate measured in units of Eddington accretion rate and stellar mass in solar masses:

$$\dot{m} = \frac{\dot{M}}{M_{\text{Edd}}} \quad \text{and} \quad m = \frac{M}{M_{\odot}}. \quad (23)$$

### 3.2. Discs with super-Eddington mass supply rates

At high local accretion rates $\dot{m} \gtrsim 1$, the radial structure of a stationary ($\dot{M} = \text{const.}$) disc around a compact object can be divided into three parts (Shakura and Sunyaev, 1973). In the outermost regions the pressure is dominated by gas and the electron scattering contribution to the opacity can be neglected. For radii

$$R_{es} < 5 \times 10^3 \dot{m}^{2/3} R_g, \quad (24)$$

electron scattering is the main source of opacity, but gas pressure is still dominant. Finally, for radii

$$R_{pr} < 1.9 \times 10^2 \dot{m}^{16/21} m^{-8/21} R_g, \quad (25)$$

the pressure is dominated by radiation, while the opacity is mainly due to electron scattering. $R_g \equiv GM/c^2$ is the gravitational radius.

The standard accretion disc model (Shakura and Sunyaev, 1973; Frank et al., 2002) applies only to geometrically thin discs, i.e., discs whose height satisfies the condition $H \ll R$. When the luminosity approaches the Eddington value this approximation breaks down.

The approximate vertical dynamical equilibrium equation

$$\frac{H}{R} \approx \frac{c_s}{v_K}, \quad (26)$$

where $c_s = \sqrt{P/\rho}$, is the isothermal sound speed, $P$ and $\rho$ are the (total, gas + radiation) pressure and density respectively, and $v_K = \sqrt{GM/R}$, can be written as

$$\frac{H}{R} \approx \frac{L(R)}{L_{\text{Edd}}}. \quad (27)$$

To derive Eq. (27), we have used the expression for the radiative pressure $P_r = 4\sigma T_c^4/3c$ and the solution for the local vertical radiative transfer equation $T_c^4 = (3/8)\kappa T_{\text{eff}}^4$ (Lasota, 2016).

Therefore already at $L \sim 0.1 L_{Edd}$ the thin-disc approximation ceases to be valid and for $L \gtrsim L_{Edd}$ one has $H \gtrsim R$ so that the accretion flow is more spherical than disc-like.

The point at which the local flux takes the Eddington value is called the spherization radius (Shakura and Sunyaev, 1973) and can be seen as the boundary separating the thin disc from a more spherical inner accretion flow.
### 3.2.1. Spherization radius

In a stationary disc ($M = \text{const}$), the local radiative flux generated by the disc's (anomalous) viscosity becomes equal to the local Eddington flux at the spherization radius $R_{\text{sph}}$. The viscous flux (Frank et al., 2002) can be written as

$$F_{\text{visc}} = \frac{3}{8\pi} \frac{R_g}{R^3} \dot{M} c^2 \left[ 1 - \left( \frac{R_{\text{in}}}{R} \right)^{1/2} \right], \quad (28)$$

and the Eddington flux is

$$F_{\text{Edd}} = \frac{L_{\text{Edd}}}{4\pi R^2}. \quad (29)$$

From $F_{\text{visc}} = F_{\text{Edd}}$, one obtains

$$R_{\text{sph}} = 15 \dot{m} R_g, \quad (30)$$

We have neglected the boundary term in square brackets since $R_{\text{sph}} \gg R_{\text{in}}$ for $\dot{m} \gg 1$ and $R_{\text{in}} \approx 6R_g$. The spherization radius defined in Eq. (30) is 1.1 times larger than the spherization radius in Shakura and Sunyaev (1973). Ohsuga and Mineshige (2014) assume $R_{\text{sph}} \equiv (\dot{M}c^2/L_{\text{Edd}}) R_S$, where $R_S = 2R_g$,<fnref n="6" /> which gives a value of $R_{\text{sph}}$ 1.3 times larger than in Eq. (30).

The outer regions of accretion discs around compact objects can be strongly X-ray self-irradiated. Then the condition for $R_{\text{sph}}$ is

$$F_{\text{Edd}} = F_{\text{irr}} + F_{\text{visc}}, \quad (31)$$

where we put (see Dubus et al., 1999):

$$F_{\text{irr}} = c \frac{\eta \dot{M}_{\text{in}} c^2}{4\pi R^2}, \quad (32)$$

where $c$ is a constant containing all the information about the disc irradiation process and $\dot{M}_{\text{in}}$ is the accretion rate onto the accretor. The spherization radius becomes

$$R_{\text{sph}} = 15 \left(1 - C \frac{\dot{m}}{2}\right)^{-1} \dot{m} R_g. \quad (33)$$

For “standard” values of $c \sim 10^{-2}$ (Dubus et al., 2001) disc irradiation can be neglected in the definition of the spherization radius for $\dot{m} \ll 200(0.01/C)$. One should notice, however, that at very high accretion rates it is nor clear how self-irradiation can proceed. Even at sub-Eddington accretion rates when applying irradiation models to observed systems, one concludes that scattering in the wind is still not sufficient to produce the observed disc heating, even in combination with direct illumination (Tetarenko et al., 2020). On the other hand, until now, the disc considered in this context were rather smaller than those supposed to be present in many ULXs so it this latter case a component of order $L_{\text{Edd}}$ emerging from the photosphere of the spherical wind, might irradiate the outer parts of the disc. Fortunately, as argued above, this effect is presumably not important at most luminosities of interest.

High accretion rates implying large optical depth and high radial velocities in the inner regions of the flow can result in photons being trapped in the accretion flow.

### 3.2.2. Trapping radius

The trapping radius is defined as the place where the photon diffusion (escape) time $H\tau/c$ equals the viscous infall time $R/v_r$

$$R_{\text{trapp}} = \frac{H \tau v_r}{c} = \frac{H \kappa \Sigma}{c} \frac{\dot{M}}{2\pi R \Sigma} = 20 \frac{H}{R} \dot{m} R_g, \quad (34)$$

where $\Sigma = 2 \int \rho dz$ is the disc surface density. For $R < R_{\text{trapp}}$ photons trapped in the flow are entrained (advected) faster than they can escape, so that if the accreting body is a black hole, most of the accretion energy is accreted by the black hole and adds to its mass. In this state there is in principle no upper limit to the accretion rate

$\dot{M}$ but the accretion luminosity grows only logarithmically with $\dot{M}$ (see Eq. (45)). In other words, very super-Eddington mass supply rates do not imply (very) super-Eddington luminosities.

Eq. (34) assumes that the accretion flow moves only in the radial direction, and that there are no outflows. Shakura and Sunyaev (1973) proposed an alternative to such an advection-dominated solution (which they considered “improbable”) where for $R \le R_{\text{sph}}$ the local flux is everywhere equal to its Eddington value.

### 3.3. Shakura-Sunyaev (1973) second (windy) solution

To keep $F_{\text{vis}} = F_{\text{Edd}}$, everywhere for $R \le R_{\text{sph}}$, we must arrange that $\dot{M} \sim R$ (see Eqs. (28) and (29)). Hence

$$\dot{m}(R) \approx \dot{m}_0 \frac{R}{R_{\text{sph}}} \quad (35)$$

Using $\sigma T_{\text{eff}}^4 \equiv F_{\text{visc}} = F_{\text{Edd}}$ the luminosity of disc described by the windy SS73 solution is

$$L_{\text{thick}} = 2 \int_{R_{\text{in}}}^{R_{\text{sph}}} \sigma T_{\text{eff}}^4 2\pi R \, dR \approx L_{\text{Edd}} \ln \dot{m}. \quad (36)$$

The total disc luminosity is then

$$L_{\text{total}} = L_{\text{thin}} + L_{\text{thick}} = 4\pi \left( \int_{R_{\text{in}}}^{R_{\text{sph}}} \sigma T_{\text{eff}}^4 R \, dR + \int_{R_{\text{sph}}}^{R_{\infty}} \sigma T_{\text{eff}}^4 R \, dR \right) \approx L_{\text{Edd}} (1 + \ln \dot{m}). \quad (37)$$

In this SS73 ‘windy’ solution, maintaining the radiation flux at the local Eddington value requires the accretion rate to decrease with radius as $\dot{M} \sim R$, i.e. that matter is expelled from the disc flow as a wind. But the trapping radius is close to the spherization radius, so one can also legitimately consider discs in which the “excess” radiation is trapped in the inflow and advected on to the accretor, while the accretion rate remains constant. We examine these advective solutions in the next two subsections. They have played a very important role in the understanding of high-rate accretion onto compact objects but, as we will see, they do not give suitable models for ULX accretion.

### 3.4. Slim discs

For a black hole accretor, the advected energy disappears behind the event horizon (growing the hole’s mass) and we can regard advection as an additional cooling mechanism. In contrast, for a neutron star accretor, the advected energy must be radiated from the stellar surface or the magnetosphere. Then the accretion luminosity violates the Eddington limit and the assumption of steady radial infall.

The energy conservation equation in a stationary disc is

$$F_{\text{visc}} = F_r + F_{\text{adv}}, \quad (38)$$

where $F_r$ is the radiative flux and $F_{\text{adv}}$ is the advective term (Lasota, 2016)

$$F^{\text{adv}} = \frac{\Sigma T v_r}{R} \frac{ds}{d \ln R} = \frac{\dot{M}}{2\pi R^2} c_s^2 \xi_a. \quad (39)$$

Here $T$ is the midplane temperature, $v_r$ the radial velocity, $c_s = \sqrt{P/\rho}$ the speed of sound, $s$ the specific entropy and $\xi_a \sim 1$ is a slowly varying function related to the entropy gradient. From Eqs. (28), (26) and (39) we get

$$\frac{F_{\text{adv}}}{F_{\text{visc}}} \sim \frac{c_s^2}{\Omega_K^2 R^2} \approx \left(\frac{H}{R}\right)^2, \quad (40)$$

so that for thin ($H/R \ll 1$) discs the advective term can be neglected. But at luminosities approaching $L_{\text{Edd}}$ this term becomes dominant, leading to

$$F_{\text{visc}} \approx F_{\text{adv}}. \quad (41)$$

<footnote n="6">⁶ In the literature both units of length, $R_S$ and $R_g$ are used.</footnote>
Solutions of this equation represent advection dominated accretion flows (ADAFs), which for high accretion rates are called “slim disc” (Abramowicz et al., 1988; Sądowski, 2009; Sadowski, 2011; Sądowski et al., 2011).

Solutions representing vertically-integrated slim-disc structure have the form

$$ \dot{m} \sim \kappa_{es} r^{1/2} \alpha \Sigma \qquad (42) $$

on the $\dot{M}(\Sigma)$ plane. $\kappa_{es}$ is the electron-scattering opacity and $\alpha$ the disc viscosity coefficient.

For these pure advection solutions $H/R \approx 1$, independently of accretion rate and radius (see Eqs. (40), (41) and Lasota, 2016). But since $F_{adv} \sim \dot{m}^2$ while $F_r \sim \dot{m}^{1/2}$ (Lasota et al., 2016), radiative cooling overtakes advection below $\dot{m} \le 1$, and the disc structure switches to the standard radiation-pressure-dominated SS73 disc:

$$ \dot{m} \sim \kappa_{es}^{-1} r^{3/2} (\alpha \Sigma)^{-1}. \qquad (43) $$

Using Eq. (27) the radiative flux of the advection-dominated flow becomes

$$ F_r = \sigma T_{\text{eff}}^4 \approx \frac{L_{\text{Edd}}}{4\pi R^2}, \qquad (44) $$

so that the luminosity of this part of the accretion flow is

$$ L_{\text{slim}} = 2 \int_{R_{\text{in}}}^{R_{\text{trans}}} \sigma T_{\text{eff}}^4 2\pi R dR \approx L_{\text{Edd}} \cdot \ln \frac{R_{\text{trans}}}{R_{\text{in}}} \approx L_{\text{Edd}} \ln \dot{m}. \qquad (45) $$

Eqs. (36) and (45) both result from the mechanical equilibrium equation of a radiation-pressure dominated disc

$$ F_r = \frac{c}{\kappa_{es}} g_z \approx \frac{c}{\kappa_{es}} \frac{GM}{R^2}, \qquad (46) $$

for $H/R \approx 1$, and $g_z \approx (GM/R^2)H/R$ (Paczyński, 1982; Poutanen et al., 2007). But although the luminosities are similar, in the first case only $\dot{m} \approx 1$ is accreted on to the compact object while in an advection dominated flow the whole of $\dot{m} \gg 1$ arrives at the accretor's surface.

Since $H/R$ in slim discs is independent of the accretion rate, one cannot get so-called “thick discs” or “Polish doughnuts” by increasing the rate at which matter is supplied to the compact object. Contrary to assertions found in the literature (e.g., Poutanen et al., 2007) slim discs are not the same objects as Polish doughnuts. As we have seen, they do not significantly increase the true emitted luminosity above $L_{Edd}$, and do not collimate or beam it such that the apparent luminosity may be super-Eddington in a restricted solid angle. Accordingly they do not offer an explanation for the defining feature of ULXs. However for a black hole ULX they may describe the fate of some of the super-Eddington mass supply.

### 3.5. Polish doughnuts

Unlike slim discs, for which $H/R \le 1$, Polish doughnuts (Kozłowski et al., 1978; Abramowicz et al., 1978; Jaroszynski et al., 1980; Paczynsky and Wiita, 1980) have by construction $H/R \gg 1$, with long central funnels along which most of the radiation is emitted. They were devised in the late 1970s and early 1980s to try to explain the high luminosities and jets of QSOs by assuming that these objects are fed at highly super-Eddington rates (see, e.g., Sikora, 1981). They are often confused with slim discs and incorrectly assumed to be advection dominated accretion flows. They have been invoked as solutions to the ULX problem, but we will see that this is not correct.

Polish doughnuts are models of stationary and axially symmetric accretion structures around black holes. All their properties are obtained from a single function $\ell(R)$ describing the specific (per unit mass) angular momentum distribution at the doughnut's photosphere. There is no explicit assumption about the doughnut's interior. One assumes only that (a.) the photosphere coincides with an equipressure surface. b.) the specific angular momentum at the photosphere is assumed to be given by some function $\ell(R)$ (one often assumes $\ell(R) = \text{const.}$), c.)

radiation is emitted at the photosphere at the local Eddington flux, that is, the local flux $F_r$ is given by

$$ F_r = \frac{c}{\kappa} g_{\text{eff}}, \qquad (47) $$

Paczyński (1982), where boldface symbols denote vectors. The specific angular momentum is assumed to be Keplerian at the inner and outer boundaries: $\ell(R_{in}) = \ell_K(R_{in})$ and $\ell(R_{out}) = \ell_K(R_{out})$, and the inner radius lies between the innermost stable circular orbit (ISCO) and the innermost bound orbit (IBCO), and $R_{IBCO} < R_{in} < R_{ISCO}$. Because of the strong radial pressure gradients (negligible for thin discs) needed to have $H/R \gg 1$, the inner flow edge is pushed inside the ISCO.

The doughnut shape is obtained from

$$ \frac{dH}{dR} = - \left( \frac{g_{\text{eff}}^R}{g_{\text{eff}}^Z} \right) \Bigg|_{Z=H} \equiv f(R, H), \qquad (48) $$

where $R, z$ are the cylindrical coordinates, and the radial component of the effective gravitational acceleration $g_{\text{eff}}^R$ is a function of $\ell(R)$.

The inner part of the doughnut forms a funnel with an opening angle obeying $\tan \alpha = 1/\chi \sim \sqrt{8\rho}$, where $\rho = (R_{in} - R_{IBCO})/R_S$ ($0 \le \rho \le 1$) and $\chi = (H/R)_{\text{max}}$ (Wielgus et al., 2016).

The luminosity is mainly emitted from the funnel, i.e.,

$$ L_{\text{apparent}} \approx \frac{L_{\text{funnel}}}{1 - \cos \alpha} \approx \frac{1}{b} L_{\text{funnel}}. \qquad (49) $$

With $b = (1 - \cos \alpha) \sim \rho$, the funnel can be very narrow, but this does not lead to an increase of the apparent luminosity. Indeed,

$$ L_{\text{apparent}} \approx \frac{1}{b} L_{\text{funnel}} = \frac{\epsilon_{\text{rad}}}{b} M c^2 \approx \text{const.} M c^2, \qquad (50) $$

because the accretion radiation efficiency $\epsilon_{\text{rad}} \sim \rho$, since by definition the binding energy at $R_{IBCO}$ is equal to zero.

So although Eq. (49) at first sight resembles Eq. (3), in Polish doughnuts one does not have $L_{\text{apparent}} \gg L_{\text{true}}$. This means that, contrary to initial appearances, Polish doughnuts are not relevant to BH ULXs (Wielgus et al., 2016), and we already know that they are not applicable to NS ULXs. Together with the similar conclusions we arrived at in Sub Section 3.4 for slim discs, this severely narrows the range of potential explanations for ULXs.

We are left with only one possibility for making an accretion flow appear super-Eddington while simultaneously ensuring that the mass supply rate very close to the accretor is not significantly super-Eddington. In the next Subsection we ask if the disc winds required by the Shakura-Sunyaev 'windy' solution can collimate the modestly super-Eddington intrinsic luminosity $L_{Edd}(1 + \ln \dot{m})$ (cf (37)) given by a super-Eddington mass supply. We note that Eq. (37) is an approximation since the spherization radius is not a uniquely defined quantity, which we should not see as a rigid limit defining the boundary of the wind. Eq. (37) is valid in both the SS73 and advection-dominated cases, but when these effects are simultaneously present one gets a slightly different formula (see e.g., Poutanen et al. 2007). In the following we therefore use $L = L_{Edd}(1 + \ln \dot{m})$ as a "universal" reference value for unbeamed luminosity.

### 3.6. Beaming by accretion disc winds

We have seen that in the SS73 'windy' solution, the accretion flow carries only the local Eddington mass rate at each accretion disc radius, while the excess is blown away in winds at each radius near the disc centre. These winds must be quasispherical, but leave open funnels around the central disc axis (since they cannot achieve zero angular momentum). The funnels offer a suitable way of collimating the accretion luminosity, as all other photon escape routes have high optical depth. We will see in Section 3.8 that numerical solutions of such flows are still some way from providing a picture which is easy to apply to observed ULXs.

Fortunately, soft X-ray observations of ULXs give a remarkable insight into beaming by disc winds. The soft components of ULXs can be
fitted with blackbody spectra, giving the temperature $T$ and apparent (assumed isotropic) luminosity $L_{\text{soft}}$. Kajava and Poutanen (2009) show that the brightest soft X-ray components (above $3 \times 10^{39}$ erg s$^{-1}$) in 9 ULXs cluster around the relation

$$L_{\text{soft}} \approx 7 \times 10^{40} T_{0.1 \text{ keV}}^{-4} \text{ erg s}^{-1}, \quad (51)$$

where $T_{0.1 \text{ keV}}$ is the blackbody temperature expressed in units of 0.1 keV/k, where $k$ is Boltzmann's constant. However, not all ULXs show such a relationship (cf Walton et al. 2020, Gúrpide et al. 2021) and – even in the absence of additional spectral components in the presence of a neutron star accretor – the combination of anisotropy and/or beaming, precession, and the lack of a suitable model for extracting $L_{\text{soft}}$ unambiguously (see, e.g. Robba et al. 2021), must result in complexity and a range of observed slopes.

At first sight (51) seems paradoxical, as one would expect a black-body source of fixed size to vary as $L \propto T^4$ instead.<fnref n="7" /> But the relation (51) instead implies that the radius of the blackbody source varies along with the luminosity and temperature. This is of course what we should expect for emission from the central part of accretion discs fed at super-Eddington rates, progressively ejecting the excess accreting matter from inside the spherization radius, as the latter varies $\propto \dot{m}$ (see Section 3.1). We can see this explicitly by parametrizing the true emitted luminosity (pre-collimation) as

$$L = 4\pi p R^2 \sigma T^4 \quad (52)$$

where $R = rR_g \propto M$ with $r \approx 15\dot{m}$ (see Eq. (30)), and $\sigma$ is the Stefan-Boltzmann constant. Setting $L = L_{\text{Edd}} \propto M$, where $l \sim 1$, and $p$ allows for geometrical projection, we can write

$$L \propto R^2 T^4 p \propto M^2 T^4 r^2 p \propto L^2 T^4 \frac{r^2 p}{l^2}, \quad (53)$$

where we use $R \propto rR_g$ at the first step and $M \propto L/l$ at the second. The final form implies $L \propto T^{-4}$. An observer assuming that the flux is isotropic with the observed value, rather than collimated, now infers a total blackbody luminosity $L_{\text{sph}} = b^{-1}L$, so that, inserting the constants,

$$L_{\text{sph}} = \frac{4\pi c^6}{\sigma \kappa^2} \frac{1}{T^4} \frac{l^2}{p b r^2} = \frac{2.8 \times 10^{44}}{T_{0.1 \text{ keV}}^4} \left( \frac{l^2}{p b r^2} \right) \text{ erg s}^{-1}. \quad (54)$$

The fact that observation gives $L \propto T^{-4}$ (cf Eq. (51)) means that the factor $l^2/pbr^2$ in this equation must be a constant. Since $r \approx (15/2)\dot{m}$ (Eq. (30)), consistency with Eq. (51) requires

$$b \approx \frac{73}{\dot{m}^2}, \quad (55)$$

where we have assumed that the slowly-varying factors $l \sim p \sim 1$. (Tight beaming means that we observe the source along the disc axis, so we view the central disc plane orthogonally, making $p \approx 1$.)

The physical origin of the simple form (55) is straightforward. The beaming factor is the total solid angle of the two open cones around the disc axis, i.e

$$b = (1 - \cos \theta) \quad (56)$$

where $\theta$ is the opening half-angle. Near the spherization radius, where the wind outflow is maximal (see Section 3.3), the flow geometry scales with $\dot{m}$. But all flows, regardless of $\dot{m}$, reduce locally to the same Eddington inflow near the accretor, so we expect $\theta \sim \text{constant} \times \dot{m}^{-1}$. For sufficiently large $\dot{m}$ we have $\theta \ll 1$ and so $b \sim \theta^2/2 \propto \dot{m}^{-2}$.

We see that the simple form (55) formally applies only for $\dot{m} \gtrsim 8$, since, at such rates, $b < 1$. Evidently it is legitimate to adopt a form of

$b$ interpolating smoothly between $b = 1$ for $\dot{m} = 0$, and the form (55) for $\dot{m} \gtrsim \sqrt{73} \approx 8$. For example, Hameury and Lasota (2020) use

$$b \approx \frac{73}{73 + \dot{m}^2}. \quad (57)$$

Importantly, although we have derived this relation by considering only the soft X-ray emission of ULXs, the fact that it is purely geometrical and involves only electron scattering means that it must apply to all forms of radiation from ULXs.

The true emission from the central accretion disc is $L \sim L_{\text{Edd}}(1 + \ln \dot{m})$, so the apparent luminosity is

$$L_{\text{sph}} = \frac{L_{\text{Edd}}}{b} \sim M \dot{m}^2 \sim \frac{\dot{M}^2}{M}, \quad (58)$$

where $\dot{M}$ is the mass supply rate, since $\dot{m} = \dot{M}/\dot{M}_{\text{Edd}} \propto \dot{M}/M^2$. So for a given mass supply rate, accretors of lower mass have higher apparent luminosities, as they are more super-Eddington, and the consequent increase in beaming outweighs the lower Eddington luminosity. This means that neutron stars are favoured as ULXs over black hole accretors in binaries where the mass transfer rate is insensitive to the accretor mass. This condition holds for mass transfer driven by the evolutionary expansion of the donor star, (see Sections 3.11 and 3.13), and also for cases where the increased mass transfer is a transient effect caused by a disc instability, e.g. in Be-X-ray binaries, and in soft X-ray transients, where the disc is subject to the thermal-viscous instability (cf Hameury and Lasota 2020).

The beaming formula (55) gives a simple explanation of the observed soft X-ray correlation (51). We can easily see that other explanations are problematic.

If we assume there is no beaming, i.e. $b = 1$, then (54) requires

$$r = 62 \left( \frac{l}{p^{1/2}} \right). \quad (59)$$

The true accretion luminosity is $L = l L_{\text{Edd}}$, and so (54) would require an accreting mass

$$M_{\text{accretor}} = \frac{500 M_{\odot}}{T_{0.1 \text{ keV}}^4} \left( \frac{1}{l} \right) \quad (60)$$

for each object.

Without beaming, the typical size, $r$, of the soft X-ray emitting region cannot self-consistently be much larger than $O(1)$ for an accretion disc spectrum. Then since $p \sim 1$, (59) requires $l \sim 1/62$, we find $M_{\text{accretor}} \gtrsim 3 \times 10^4 M_{\odot}$ from (60), i.e. a significantly massive IMBH.

A magnetar model would require a neutron star mass, so (60) requires $l \gtrsim 350$. But then from (59) there must be an implausibly large soft X-ray region $r \sim 1.6 \times 10^4$, which is $\sim 500$ neutron star radii.<fnref n="8" /> In addition, the implied luminosity is $L \sim 5 \times 10^{40}$ erg s$^{-1}$, beyond the domain of validity for magnetar models for ULXs (Canuto et al., 1971).

In the context of the model described above, the observed correlation (51) provides strong evidence in favour of accretion-disc wind beaming (and indeed, yields sensible results: King and Lasota 2020). We also note that the innermost (hottest) regions will naturally be the most collimated and therefore beamed, whilst one would expect radiation produced from around $r_{\text{sph}}$ to be more isotropic. The final luminosity we observe may therefore be a fairly complicated function and will naturally be inclination dependent. But regardless of the final functional form of the beaming, it is inevitable that where the wind is optically thick and subtends a large scaleheight (as seen in every MHD simulation of supercritical accretion), collimation and beaming must result.

<footnote n="7">⁷ Indeed one does find $L \propto T^4$ for sources which are permanently below $3 \times 10^{39}$ erg s$^{-1}$ (Kajava and Poutanen, 2009). It is likely that these sources are black holes with masses large enough that their luminosities are sub-Eddington.</footnote>

<footnote n="8">⁸ We will see in Section 3.14 that a large photosphere is possible in ULXs producing strong outflows, i.e. if a strong accretion-disc wind produces beaming.</footnote>
<figure><img src="image_15.png" /><figcaption>Fig. 13. The $L_{39} - \dot{v}_{-10}$ diagram for XRPs and PULXs. Red dots: the ten PULXs with known spin-up rates (values from KLK17). Blue diamonds: selected (for comparison) X-ray pulsars (see Table A.8 - in the Appendix).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 14. Schematic picture of a PULX according to Mushtukov et al. (2017, 2019). $R_m$ and $R_{sp}$ are respectively the magnetospheric and spherization radius. $\lambda$ is the angle between the disc plane and the magnetic axis. Accreting matter not lost in the wind is accreted on to the central object, forming an envelope which is optically thick at accretion rates typical of PULXs. In the KLK17 model, the configuration is similar, but $R_m \lesssim R_{sp}$ and the wind does not end up at the magnetospheric radius. (see Fig. 15). Note that the Figure of Mushtukov et al. (2017, 2019) shown here assumes that the neutron star spin is aligned with the disc axis, although this is not explicitly stated. We argue in Section 3.7.4 that this assumption is in general not justified. In addition the configuration shown appears to be axisymmetric, i.e. has the magnetic axis parallel to the spin (and disc) axes. This would of course preclude pulsing, so there must in reality be some deviation from axisymmetry.</figcaption></figure>

## 3.7. Models of pulsing ULXs (PULXs)

The discovery (beginning in 2014) of a small group of ULXs showing coherent pulsing has had a significant effect on the study of ULXs in general. First, it subjected existing models of ULXs to a stringent new test, as the mass of the accretor must be close to $1.4 M_\odot$. As remarked above, it removed the main motivation for the IMBH model, and we will see in Section 3.7.1 below, that it fitted naturally and self-consistently into the disc-wind beaming picture, even though this was not developed with PULXs in mind.

But the removal of one candidate model (IMBH) did not simplify the discussion, as the discovery of PULXs naturally stimulated suggestions that magnetic effects on the accretion flow might be significant.

### 3.7.1. The KLK17 model of PULXs

The discovery of PULXs has shown that ULX apparent luminosities can be highly super-Eddington. The luminosities of the classic (non Be-X) PULXs range from ~20 to ~385 times the Eddington value for a neutron star (see Table 2). For normal (pulsar-like) magnetic fields ~$10^{11}-10^{13}$ G, such values exclude the possibility that they correspond to the actual luminosities since this would require implausibly high accretion rates (see Eqs. (36) and (45)). We will show in Section 3.9 that the hypothesis that PULXs contain magnetars is extremely unlikely. Then the only physically reasonable option left is that the luminosity of magnetic ULXs is beamed according to Eq. (58). Given this, King and Lasota (2016) proposed a beamed-emission model for the first detected PULX, M82 ULX-2. This model was successfully used by King et al. (2017) (hereafter KLK17) to explain the properties of the three first-discovered PULXs, and later King and Lasota (2019) applied this model to the growing number of observed PULXs, including the Be-PULXs. King and Lasota (2020) explained how one obtains significant pulsed X-ray luminosity fraction in beamed radiation (see Section 3.7.5). It is assumed throughout that for neutron star magnetic fields of the standard X-ray binary strength (i.e. $B \lesssim 10^{13}$ G) the surface emission is locally isotropic.

Kluźniak and Lasota (2015) noted that PULXs are sharply distinct from other X-ray pulsars, not simply in their luminosities (by definition $L > 10^{39}$ erg s$^{-1}$) but also in their very large spinup rates ($\dot{v} > 10^{-10}$ s$^{-2}$ – up to two orders of magnitude larger than normal XRP). Between superoutbursts, transient Be-PULXs are normal spinning-down X-ray pulsars, but when their luminosity reaches ~$10^{39}$ erg s$^{-1}$ they become rapidly spinning-up sources with $\dot{v} \gtrsim 10^{-10}$ s$^{-2}$ (see Table 2), demonstrating that such high spin-up rates are a generic property of PULXs.<fnref n="9" />

The luminosities and spin-up rates of PULXs are tightly correlated (Fig. 13). This correlation strongly implies that the dominating torque in the system is provided by accretion, as assumed in KLK17 and Vasilopoulos et al. (2018), since

$$ \dot{v} = \frac{\dot{J}(R_M)}{2\pi I} = \frac{\dot{M}(GM R_M)^{1/2}}{2\pi I} \propto \dot{M}^{6/7} \qquad (61) $$

<footnote n="9">⁹ One might insist that it is a coincidence that Be-X systems show the same spin-up as other PULXs when they reach the same luminosities, but here one rather expects that “coincidences mean you’re on the right path” (Van Booy, 2009).</footnote>
where $R_M \propto \dot{M}^{-2/7}$ (e.g. Frank et al. 2002) is the magnetospheric radius and $I$ the neutron star’s moment of inertia (see Fig. 14).

If we define PULXs by the twin properties of luminosities larger than $10^{39}$ erg s$^{-1}$ and spin-up rates larger than ~ $10^{-10}$ s$^{-2}$, these two quantities should form the basis for models of these sources. At the very least, any cogent model must be able to reproduce these values. KLK17 chose the first option, using the observed quantities $L_X$, $\dot{v}$ as input for their PULX model. Although there are no direct mass measurements for neutron stars in PULXs, the likely mass range is limited. KLK17 make no assumptions about the magnetic field strengths, and these are outputs of the model. The resulting values lie in the typical range for standard pulsing X-ray binaries (see Table 3), and compare favourably with the dipole magnetic field strengths measured from CRSFs (one of them for a non-pulsing neutron star ULX – see also Section 2.2.3).

For the mass inflow, KLK17 use the super-critical (‘windy’) solution of Shakura and Sunyaev (1973) described in Section 3.3. In this picture, the super-Eddington part of the mass supply below $R_{\text{sph}}$ is eventually expelled as a wind. This simultaneously provides the beaming, making PULXs appear super-Eddington and explains the strong outflows observed from PULXs. The neutron star gains mass only at its effective Eddington rate (Eq. (36)).

The spherization radius defined in Eq. (30) can be written as

$$R_{\text{sph}} \approx 2.3 \times 10^6 \dot{m}_0 m_1 \text{ cm}, \qquad (62)$$

where $\dot{M}_0 = \dot{m}_0 \dot{M}_{\text{Edd}}$ is the accretion rate at $R_{\text{sph}}$, assumed equal to the mass transfer rate. For $R > R_{\text{sph}}$ the disc is assumed to be a standard SS73 accretion disc. From Eq. (63) we have

$$\dot{M}(R) \approx \dot{m}_0 \dot{M}_{\text{Edd}} \frac{R}{R_{\text{sph}}}. \qquad (63)$$

for $R < R_{\text{sph}}$.

KLK17 assume that the Shakura and Sunyaev (1973) ‘windy’ model describes the accretion flow between $R_{\text{sph}}$ and the magnetospheric radius $R_M$ defined by the equation (Frank et al., 2002)

$$R_M = 1.2 \times 10^8 q \dot{m}^{-2/7} m_1^{-3/7} \mu_{30}^{4/7} \text{ cm}, \qquad (64)$$

where $q \sim 1$ is a factor taking into account the geometry of the accretion flow at the magnetosphere. Then from Eq. (63)

$$\dot{M}(R_M) \approx \dot{m}_0 \frac{R_M}{R_{\text{sph}}}. \qquad (65)$$

Applying the model to PULXs with measured spin-up gives $R_M < R_{\text{sph}}$ in all cases. This implies that in these objects the accretion flow inside the magnetosphere is highly super-Eddington. Since it is hypersonic but forced to follow the magnetic fieldlines it is highly dissipative, and one expects it to generate an outflow similar to that of Shakura and Sunyaev (1973), limiting the local luminosity to its Eddington value.

The total luminosity is then given by Eq. (36)

$$L \approx L_{\text{Edd}} [1 + \ln \dot{m}_0]. \qquad (66)$$

The luminosity from both parts of the super-Eddington outflow is assumed to be beamed. Outside the magnetosphere the beaming factor is taken to be

$$b \approx \frac{73}{\dot{m}^2}, \qquad (67)$$

as in Eq. (55).

For accretion rates such that radiation is geometrically beamed as described by Eqs. (66) and (67), one can deduce $\dot{m}_0$ from the observed X-ray luminosity $L = L_X$ by combining these two equations into

$$L_X \approx 2.0 \times 10^{36} \dot{m}_0^2 [1 + \ln \dot{m}_0] m_1 \text{ erg s}^{-1}. \qquad (68)$$

Having $\dot{m}_0$ one obtains $R_{\text{sph}}$ from Eq. (62).

The second observed quantity, the spinup, follows from Eqs. (61) and (64) as

$$\dot{v} = 3.3 \times 10^{-11} q^{1/2} \dot{m}^{6/7} m_1^{6/7} \mu_{30}^{2/7} I_{45}^{-1} \text{ s}^{-2}. \qquad (69)$$

This gives the values of the magnetic moment $\mu = BR^3$ (where B is the field and R the neutron-star radius) which in turn allows one to calculate the values of $R_M$ and $\dot{m}(R_M)$ from Eqs (64) and (65):

$$\mu_{30} = 0.04 q^{-7/4} \dot{v}_{-10}^{3/2} m_1^{1/2} I_{45}^{3/2} \text{ G cm}^3 \qquad (70)$$

$$R_M = 1.0 \times 10^7 \dot{v}_{-10}^{2/3} m_1^{1/3} I_{45}^{2/3} \text{ cm} \qquad (71)$$

$$\dot{m}(R_M) = 9.9 \dot{v}_{-10}^{2/3} m_1^{-2/3} I_{45}^{2/3}. \qquad (72)$$

where $10^{-10} \dot{v}_{-10} = \dot{v}$.

The results are presented in Table 3, with the values of the equilibrium period $P_{\text{eq}} = 0.23 q^{7/6} m_1^{-1/3} \mu_{30}^{2/3}$ s, and the lower limit on the time to reach equilibrium at the present spin-up rate

$$t_{\text{eq}} \equiv \frac{1}{\dot{v}} \left( \frac{1}{P_{\text{eq}}} - \frac{1}{P} \right). \qquad (73)$$

(The corotation radius $R_{\text{co}} \equiv (GM P_s^2 / 4\pi^2)^{1/3}$ is always larger than the magnetospheric and spherization radii. Thus only possibly NGC 5907 ULX1 is close to its equilibrium spin, and none of the observed PULXs is a propeller.)

For NGC300 ULX-1, the KLK17 model predicts a neutron-star magnetic field $B > 5 \times 10^{12}$G (since $q < 1$ and $m_1 > 1$, the entries provided in the 4th column of Table 3 are a factor of few lower than the predicted magnetic field values). A magnetic field of very similar strength is deduced by Walton et al. (2018a) from the CRSF inferred to be present in the X-ray spectrum of this PULX. According to the KLK17 model, all seven PULXs have magnetic fields in the range $10^{11} - 10^{13}$G, i.e. always in the standard pulsing X-ray binary range and below the value defining magnetars.

The predicted values of the beaming factor $b$ are in the range ~ 0.03 – 0.5, except for the very luminous source in NGC 5709, for which $b \approx 0.01$. For NGC 300 ULX-1, the model gives $b \approx 0.13$: Binder et al. (2018) use the KLK17 model and obtain $b \approx 0.25$, because they deduce $\dot{m}_0$ from the average rather than the maximum luminosity. The beaming factor $b \approx 0.12$ for P13 in NGC7793 is consistent with observations of the X-ray irradiation of the neutron-star companion (Motch, 2018). According to the KLK17 model, the luminosity of most of the known PULXs is only mildly beamed. This is worth stressing since critics of the KLK17 model often claim that there are problems with ‘strong beaming’.

We see that for all seven PULXs with known spin-up rates $\dot{v}$, one has $R_M \lesssim R_{\text{sph}}$, which is probably the condition for observing pulses at all if the mass transfer is super-Eddington. The small difference between $R_{\text{sph}}$ and $R_M$ means that the flow is strongly super-Eddington on reaching $R_M$. Most of this cannot land on the neutron star (let alone its polecaps only) and so must be ejected. It is worth stressing that nothing in the assumptions of the model guarantees its self-consistency, i.e. that $R_{\text{sph}} > R_M$. But in every case this is satisfied by the deduced values. Further, exactly the same equations describe both the Be-star magnetic accretors and those PULXs which do not have Be-star companions: the first group cannot of course be modelled at all by assuming very strong dipole magnetic fields.

The spin-up time-scales $t_{\text{eq}}$ are probably much shorter than the lifetimes of the individual PULXs (see Table for NGC 1313 X-2) so it seems very unlikely that we observe these systems during their only approach to spin equilibrium. Instead, although most of them are probably rather far from their $P_{\text{eq}}$ they have nevertheless the time to alternate spin-up and spin-down phases. Evidently, we can only see these systems during spin-up phases (so that $\dot{v}$ has its maximum value) because centrifugal repulsion during spin-down presumably reduces the accretion rate and so the luminosity. In addition the magnetic fields of the PULXs appear to be significantly lower than is usual for a new-born neutron star ($10^{12}-10^{13}$G) which is consistent with hypothesis that accretion of even a relatively small mass severely reduces the surface fields of neutron stars. This is this is central to the concept of
<table><caption>**Table 3**

KLK17 model: derived properties of neutron star ULXs.</caption><thead><tr><th>Name</th><th>m&#x0307;<sub>0</sub></th><th>b</th><th>B<sub>0</sub> q<sup>7/4</sup>m<sub>1</sub><sup>-1/2</sup>I<sub>45</sub><sup>-3/2</sup>R<sub>6</sub><sup>3</sup> [G]</th><th>R<sub>sph</sub>m<sub>1</sub><sup>-1</sup> [cm]</th><th>R<sub>M</sub>m<sub>1</sub><sup>-1/3</sup>I<sub>45</sub><sup>-2/3</sup> [cm]</th><th>P<sub>eq</sub>q<sup>-7/6</sup>m<sub>1</sub><sup>1/3</sup> [s]</th><th>t<sub>eq</sub> [yr]<sup>1</sup></th></tr></thead><tbody><tr><td>M82 ULX2</td><td>46</td><td>0.03</td><td>4.0 &times; 10<sup>10</sup></td><td>1.1 &times; 10<sup>8</sup></td><td>1.0 &times; 10<sup>7</sup></td><td>0.02</td><td>15600</td></tr><tr><td>NGC 7793 P13</td><td>25</td><td>0.12</td><td>1.1 &times; 10<sup>11</sup></td><td>5.8 &times; 10<sup>7</sup></td><td>1.6 &times; 10<sup>7</sup></td><td>0.09</td><td>1386</td></tr><tr><td>NGC5907 ULX1</td><td>95</td><td>0.01</td><td>9.4 &times; 10<sup>12</sup></td><td>2.2 &times; 10<sup>8</sup></td><td>1.1 &times; 10<sup>8</sup></td><td>1.86</td><td>0</td></tr><tr><td>NGC300 ULX1</td><td>24</td><td>0.13</td><td>5.3 &times; 10<sup>11</sup>&#x2019;</td><td>5.5 &times; 10<sup>7</sup></td><td>3.2 &times; 10<sup>7</sup></td><td>0.19</td><td>297</td></tr><tr><td>M51 ULX7<sup>a</sup></td><td>28</td><td>0.09</td><td>1.9 &times; 10<sup>11</sup></td><td>6.4 &times; 10<sup>7</sup></td><td>2.0 &times; 10<sup>7</sup></td><td>0.08</td><td>1337</td></tr><tr><td>M51 ULX7<sup>b</sup></td><td>28</td><td>0.09</td><td>6.9 &times; 10<sup>9</sup></td><td>6.4 &times; 10<sup>7</sup></td><td>4.6 &times; 10<sup>6</sup></td><td>0.01</td><td>~10<sup>5</sup></td></tr><tr><td>NGC 1313 X-2</td><td>46</td><td>0.03</td><td>5.3 &times; 10<sup>10</sup></td><td>1.1 &times; 10<sup>8</sup></td><td>1.8 &times; 10<sup>6</sup></td><td>0.03</td><td>8641</td></tr><tr><td>SMC X-3<sup>Be</sup></td><td>18</td><td>0.23</td><td>2.3 &times; 10<sup>10</sup></td><td>4.1 &times; 10<sup>7</sup></td><td>7.1 &times; 10<sup>6</sup></td><td>0.006</td><td>76621</td></tr><tr><td>NGC 2403 ULX<sup>Be</sup></td><td>13</td><td>0.43</td><td>2.5 &times; 10<sup>11</sup></td><td>3.0 &times; 10<sup>7</sup></td><td>2.3 &times; 10<sup>7</sup></td><td>0.16</td><td>578</td></tr><tr><td>Swift J0243.6+6124<sup>Be</sup></td><td>14</td><td>0.37</td><td>1.3 &times; 10<sup>11</sup></td><td>3.2 &times; 10<sup>7</sup></td><td>1.7 &times; 10<sup>7</sup></td><td>0.07</td><td>2047</td></tr><tr><td>RXJ0209.6-7427<sup>Be</sup></td><td>17</td><td>0.25</td><td>5.3 &times; 10<sup>10</sup></td><td>3.2 &times; 10<sup>7</sup></td><td>1.8 &times; 10<sup>6</sup></td><td>0.03</td><td>8665</td></tr><tr><td>NGC 1313 ULX<sup>Becc</sup></td><td>15</td><td>0.32</td><td>?</td><td>3.5 &times; 10<sup>7</sup></td><td>?</td><td>?</td><td>?</td></tr><tr><td>M51 ULX8</td><td>17</td><td>0.25</td><td>~3 &times; 10<sup>11</sup>&#x2019;</td><td>3.2 &times; 10<sup>7</sup></td><td>2.7 &times; 10<sup>7</sup>&#x2019;</td><td>Non-pulsing</td><td></td></tr></tbody></table>

– Systems with $^{Be}$ superscript have Be–star companions.

¹ – calculated using the value of $P_{eq}q^{-7/6}m_1^{1/3}$ from the previous column.

† – consistent with observations Walton et al. (2018c).

• – from observations (Brightman et al., 2018; Middleton et al., 2021b).

• – for a ~10¹² G magnetic field.

a – for $\dot{v} = 2.8 \times 10^{-10}$; b – for $\dot{v} = 3.1 \times 10^{-11}$ (Vasilopoulos et al., 2019).

c – unknown $\dot{v}$.

pulsar recycling, which is implicated in the production of millisecond pulsars (see Bhattacharya and van den Heuvel, 1991, and references therein). Things are even more complicated by the fact that the vast majority of ULXs do not pulse, despite containing (magnetized) neutron stars, which suggests that alignment of spin and central disc axes is rapidly suppressing (perhaps temporarily?) the pulsations, which are observed to be transient anyway. As expected from the KLK17 model, Be–star systems are prominent among the PULXs, since accretion is relatively weak and transient.

Parametrizing $R_M = f R_{sph}$, with $f \approx 0.3-0.9$, then from Eq. (65), $\dot{m}(R_M) = f \dot{m}_0$ (KLK17) and from Eqs (69) and (72) we find

$$ \dot{v} = 7.8 \times 10^{-10} f^{7/6} q^{7/6} \mu_{30}^{2/3} m_1^{-1/7} I_{45}^{-1} \quad (74) $$

in agreement with the observed spin-up values of PULXs (cf Fig. 13).

It has been claimed (e.g. Mushtukov et al., 2020) that the KLK17 picture suffers from fundamental weaknesses since ‘strong beaming’ would be incompatible with the high observed pulsed fraction of the ULX radiation. As we noted (two paragraphs above), the beaming is not ‘strong’. Further, in Sect. (3.7.5) using the results of King and Lasota (2020), we will show that the KLK17 model is fully compatible with the observed X-ray pulse fraction (which varies in time and frequency), and also explains the very small number of observed PULX.

### 3.7.2. The g úrpide et al. PULX models

Gúrpide et al. (2021) suggest an alternative scenario in which the hardest ULXs are powered by strongly magnetised neutron stars, so that the high-energy emission is dominated by the hard direct emission from the accretion column. It is assumed that one can explain softer sources as weakly magnetised neutron stars or black holes, in which the presence of outflows naturally explains their softer spectra through Compton down–scattering and their spectral transitions. Outflows are also invoked to explain the dilution of the pulsed emission in sources containing neutron stars. According to these authors NGC 7793 P13 and NGC 1313 X–2 would be strongly magnetized with $R_M > R_{sph}$ in contradiction with the prediction of the KLK17 model. However, they find that NGC 5907 ULX–1 appears harder when dimmer. This is difficult to reconcile with their scenario and they have to conclude that for this source $R_M \approx R_{sph}$, in agreement with the KLK17 model prediction. However, Gúrpide et al. (2021) do not take into account the possibility that the flow geometry is as described in Section 3.7.4. In this case, under appropriate conditions one sees the “magnetic column” directly when $R_M > R_{sph}$. In any case, one would expect that the flux–hardness relation would also depend on the intrinsic properties of the emitter and not on the source geometry only. The

model geometry of King and Lasota (2020) is used by Pintore et al. (2021) with additional ingredients to interpret the spectral variation during the flares of NGC 4559 X7.

Although M51 ULX–8 is not a PULX, we know that it is a magnetic system. Assuming that its X-ray radiation is beamed, one can obtain its $\dot{m}_0$ and so both $R_M$ and $R_{sph}$. As in PULXs this gives $R_M \lesssim R_{sph}$ (see Table 3). King and Lasota (2019) (see also Section 3.7.5) speculate that M51 ULX–8 might be a PULX, and suggest that pulsations may be seen in future observations of this system. Brightman et al. (2018) point out that rather long XMM-Newton observations would be needed to detect ~1s pulsations if the pulsed fraction is $\lesssim 45\%$, as in most PULXs. (Interestingly the exception is NGC 300 ULX1, in which Walton et al. (2018a) inferred the presence of a cyclotron feature).

Although simple, the KLK17 model correctly reproduces the main observed properties of PULXs and has passed this test with each newly discovered PULX, in particular with the Be–PULX systems (see Table 3). In the KLK17 picture, PULXs are “normal” (non–magnetar) X-ray pulsars caught in a high accretion–rate episode of their binary evolution (non–Be systems: see Sections 3.11, 3.12) or orbital history (Be systems: see Section 3.15).

### 3.7.3. The Erkut et al. (2020) PULX models

Erkut et al. (2020) (hereafter ETEA20) base their models of PULX on the well–known approach to the disc–magnetosphere interaction best known from the paper of Ghosh and Lamb (1979). They assume that the luminosity in the magnetosphere is at most critical, using, as the relation between the observed X-ray luminosity and the accretion rate:

$$ L_X = \frac{\eta}{b} \dot{M}_* c^2, \quad (75) $$

where $\dot{M}_*$ is the accretion rate at the neutron-star surface and $b$ is a beaming factor, and assumes $\dot{M}_* = \dot{M}_{in}$, where $\dot{M}_{in}$ is the accretion rate at the innermost disc radius $R_{in}$. (This is separated from the magnetosphere by a boundary region of width $\Delta R \equiv \delta R \cdot b$.)

In other words

$$ L_X \le L_c \approx \left[ 1 + 311 \left( \frac{B}{B_c} \right)^{4/3} \right] L_E, \quad (76) $$

where $B_c \equiv m_e^2 c^3 / \hbar e = 4.4 \times 10^{13} \text{G}$ (see Eq. (94)).

The beaming considered by ETEA20 is assumed *magnetic* (i.e. not due to an outflow in origin, with beaming factor

$$ b_M = \frac{A_c}{\gamma A}, \quad (77) $$
where $A_c$ is the polar cap area, $A = 4\pi R_*^2$ is the total surface area of the neutron star, and $\gamma < 1$ “is a normalization constant corresponding to the maximum fractional area of the polar cap” (Erkut et al., 2020). Since they are physically different, and act in different parts of the accretion flow in PULXs, direct comparison of the values of ETEA20’s $b_M$ with those of KLK17’s outflow-generated $b$ is not of much interest.

The accretion flow between $R_{in}$ and $R_{sph}$ (defined as in Eq. (30), with a factor 27/2, instead of 15) is allowed to be super-Eddington and is assumed to be described by the ‘windy’ Shakura and Sunyaev (1973) solution (Eq. (63)).

The inner radius of the disk is determined by the balance between magnetic and material stresses

$$ \frac{d}{dR} (\dot{M} R^2 \Omega) = -R^2 B_{\phi}^{+} B_z, \qquad (78) $$

where $\Omega$ is the angular velocity of the innermost disc matter within the boundary region, $B_z$ is the poloidal magnetic field, $B_{\phi}^{+} = \gamma_{\phi} B_z$ is the toroidal magnetic field at the surface of the disc, and $\gamma_{\phi}$ is the azimuthal pitch that can be expressed as $\gamma_{\phi} \approx \omega_{\star} - 1$, where $\omega_{\star}$ is the so-called “fastness parameter”

$$ \omega_{*} \equiv \frac{\Omega_{*}}{\Omega_{K} (R_{in})} = \left( \frac{R_{in}}{R_{co}} \right)^{3/2}. \qquad (79) $$

At spin equilibrium, the fastness parameter is assumed to be equal to a (model-dependent) critical value $\omega_{\star} = \omega_c$. ETEA20 use $\omega_c = 0.75$.

Integrating Eq. (78), using Eq. (75) and defining poloidal field $B_z \approx -\mu/R^3$ one obtains an expression for the inner disc radius

$$ R_{in} = \left( \frac{\mu^2 \sqrt{GM_*} \delta}{b_M R_* L_X} \right)^{2/7}. \qquad (80) $$

Since one can write for a strong enough magnetic field

$$ \frac{A_c}{A} \approx \frac{R_*}{4R_{in}} \cos^2 \alpha, \qquad (81) $$

where $\alpha$ is the angle between the rotation and magnetic axes (Frank et al., 2002, the star’s spin and the disc rotation axes are assumed to be aligned), one gets $b = (\cos^2 \alpha/4\gamma) R_* / R_{in}$.

Using $\mu = B_z R_*^3$ one can then write

$$ b_M = \left( \frac{\cos^7 \alpha L_X}{32\gamma^{7/2} \sqrt{GM_*} B^2 R_*^{3/2} \delta} \right)^{2/5}, \qquad (82) $$

for the beaming factor and express the fastness parameter as

$$ \omega_{*} = \left( \frac{\gamma \sqrt{GM_{*}} B^{2} R_{*}^{4} \delta}{\cos^{2} \alpha L_{X}} \right)^{3/5} \frac{2\pi}{P \sqrt{GM_{*}}}, \qquad (83) $$

where $P = 2\pi/\Omega_*$.$^3$ In addition to observed quantities, $L_X$ (ETEA20 use the observed flux $F_X$, but since, with one exception, PULXs are extragalactic, this precaution does not seem to be necessary), and $P$, eqns (82) and (83) contain the unobservable (at least in PULXs) quantities $\alpha$ and $\gamma$. They can be eliminated by the use of the torque equation that ETEA20 write under the form

$$ 2\pi I \dot{v} = \omega_{*}^{1/3} n_0 \sqrt{GM_* R_{co}} \dot{M}_*, \qquad (84) $$

where $n_0$ is a constant of order unity corresponding to a dimensionless torque. Given the neutron-star mass (ETEA20 try $M_* = 2 M_\odot$ for M51 ULX-7), moment of inertia and radius, one then obtains two equations for the four quantities defining the physics of the system: $B$, $b_M$, $\omega_*$ and $\delta$ as a function of the observables $L_X$, $P$ and $\dot{P}$ (or $\dot{v}$):

$$ B = \frac{\omega_{*}}{R_{*}^{3}} \sqrt{\frac{2GM_{*}I|\dot{P}|}{\pi n_{0} \delta}}, \qquad (85) $$

and

$$ \omega_{*} = 16\pi^4 GM_{*} \left( \frac{I |\dot{P}|}{n_0 L_X P^{7/3} R_{*} b_M} \right)^3. \qquad (86) $$

These two equations allow only us only to limit the range of magnetic-field strength and fastness-parameter values compatible with observation. The boundary region parameter is assumed to lie in the range $0.01 \le \delta \le 0.3$. The determined range of physical parameters can be made narrow by assuming that $L_X = L_c$, or $L_X < L_c$ and the fastness parameter can be set equal to the equilibrium value $\omega_* = \omega_c = 0.75$.

In their paper, ETEA20 consider several alternative scenarios: (i) the PULXs are far from spin equilibrium and an efficient standard spin-up torque is used to account for the observed spin-up rates, (ii) the PULXs are so close to spin equilibrium that the fastness parameter is given by its critical value, (iii) the X-ray luminosities of the PULXs can be well represented by the maximum critical luminosity, (iv) the conditions described in (i) and (iii) both apply, (v) the conditions described in (ii) and (iii) both apply, and (vi) the X-ray luminosities of the PULXs are subcritical and either the condition described in (i) or the condition described in (ii) applies. They found that the narrowest ranges for $B$, $b_M$, and $\omega_*$ are found when they use the critical luminosity condition along with either the observed spin-up rates or the spin-equilibrium condition. Scenario (iv), based on the observed spin-up rates at critical luminosity, is the only one that works well for all the PULXs, implying $B$ in the $10^{11}-10^{13}$ G range. The results for scenario (iv) are presented in Table 4 (Table 5 in Erkut et al., 2020). The ETEA20 models do not allow for a prediction of the pulsed fraction of the emitted radiation, nor the pulse shapes, and the authors refer to the “accretion curtain” of Mushtukov et al. (2017) as a possible source of the pulsed luminosity.

The conclusions of the Erkut et al. (2020) paper, i.e. that PULXs have sub-critical magnetic fields and their luminosity is subject to modest beaming, are the same as those obtained by King and Lasota (2016), King et al. (2017), King and Lasota (2019, 2020), but they differ in their physical motivations. First, the beaming mechanisms are different in the two approaches. Second, while ETEA20 require the accretion flow to be always at most critical, and assume that the accretion rate is constant inside the magnetosphere, KLK17 allow the accretion flow to be super-Eddington and losing mass also when following the magnetic fieldlines. It is perhaps worth noting that these two very different models both exclude magnetar-strength fields.

### 3.7.4. The origin of “sinusoidal” pulses in PULXs

X-ray pulse light curves for PULXs are generally described as being ‘sinusoidal’, i.e. without obvious eclipses, and continuously modulated (see Section 2.3.2). In other words, the X-ray light curves of PULXs are never flat and never zero.

Mushtukov et al. (2017, 2019) argued that this was understandable if, perhaps because of the large optical depth, the photospheric surface emitting the X-rays occupied a large fraction of the magnetosphere. Then this region would neither be permanently in view, nor periodically occulted. Accordingly, Mushtukov et al. (2017, 2019) proposed that the X-rays are emitted by an optically thick envelope defined by the accretion flow over the neutron-star magnetosphere. This is consistent with the accretion along fieldlines being a bending hypersonic flow, which must therefore shock as suggested by King and Lasota (2019) (see also Section 3.7.1). The spin modulation occurs because the magnetic axis is inclined to the spin axis.

Mushtukov et al. (2017, 2019) do not attempt to reproduce the observed PULX lightcurves and do not explain why magnetized systems such as ULX-8 in M51 do not pulse. The simplest explanation would be that the spin and magnetic axes are aligned, but then one would have to explain why the few pulsing systems are different in this respect from the probable majority of NULXs.

### 3.7.5. Beamed pulsed X-ray radiation

An argument often made against the beaming model of PULXs is the belief that the observed high ($\le 45\%$) pulsed fraction of the X-ray emission cannot appear in radiation passing through a funnel because reflections from the funnel walls would destroy the signal coherence.
<table><caption>**Table 4**
ETEA20 models of PULXs [scenario (iv)].</caption><thead><tr><th>Source</th><th>L<sub>X</sub> (erg s<sup>-1</sup>)</th><th>ψ (s<sup>-2</sup>)</th><th>b<sub>M</sub></th><th>ω<sub>*</sub></th><th>B (G)</th></tr></thead><tbody><tr><td>M82 X-2</td><td>6.5 × 10<sup>39</sup></td><td>1.1 × 10<sup>-10</sup></td><td>0.16–0.25</td><td>0.060–0.22</td><td>2.0–2.9 × 10<sup>12</sup></td></tr><tr><td>ULX NGC 7793 P13</td><td>9.4 × 10<sup>39</sup></td><td>1.2 × 10<sup>-11</sup></td><td>0.021–0.027</td><td>0.065–0.14</td><td>1.2–3.2 × 10<sup>11</sup></td></tr><tr><td>ULX NGC 5907</td><td>1.5 × 10<sup>41</sup></td><td>3.9 × 10<sup>-9</sup></td><td>0.19–0.30</td><td>0.16–0.61</td><td>2.7–3.8 × 10<sup>13</sup></td></tr><tr><td>NGC 300 ULX1</td><td>4.3 × 10<sup>39</sup></td><td>5.6 × 10<sup>-10</sup></td><td>~1</td><td>0.015</td><td>6.8 × 10<sup>12</sup></td></tr><tr><td>M51 ULX-7</td><td>7.0 × 10<sup>39</sup></td><td>3.1 × 10<sup>-11</sup></td><td>0.052–0.0</td><td>0.019–0.063</td><td>6.2–10 × 10<sup>12</sup></td></tr><tr><td>Swift J0243.6+6124</td><td>1.7 × 10<sup>39</sup></td><td>2.4 × 10<sup>-10</sup></td><td>~1</td><td>0.041</td><td>3.5× 10<sup>12</sup></td></tr></tbody></table>

Note. In this scenario the magnetic field, beaming fraction, and fastness parameter for Swift J0243.6+6124 are only marginally obtained at $b_M \sim 1$ for $M_* = 1.2M_\odot$ and $R_* = 10$ km, and no common solution is found for NGC 1313 X-2 (even a mass as low as $0.9M_\odot$ does not help). For all other sources $M_* = 1.4M_\odot$ and $R_* = 10$ km are used.

<figure><img src="image_17.png" /><figcaption>**Fig. 15.** Effect of spin orientation on pulsing. Left: the neutron-star spin and the accretion disc beaming axes are strongly misaligned, so that a significant part of the pulsed emission is periodically occulted by the neutron-star body at certain spin phases. This gives a large pulse fraction. Right: the neutron-star spin and central disc axes are assumed to be substantially aligned (this special condition is assumed without any explicit statement in papers asserting that beaming cannot produce pulsed emission). Much of the primary X-ray emission is scattered by the walls of the beaming ‘funnel’ before escaping. In this kind of configuration the pulse fraction is reduced or suppressed. King and Lasota (2020).</figcaption></figure>

We note first that although many (if not most) ULXs probably contain magnetised neutron stars, only a small fraction show pulses. Second, King and Lasota (2020) point out that the argument above assumes (without any explicit statement) that the spin axis of the neutron star is aligned with the accretion disc axis, and so with the beam axis.<fnref n="10" /> In fact observations of high-mass X-ray binaries show that the neutron star spin is generally *not* aligned with the disc axis, probably because of the asymmetric nature of the supernova giving birth to the neutron star.

The beamed emission will appear sinusoidal, (i.e. the light curve has no constant intervals) if three conditions hold simultaneously:

(a) there is significant misalignment of spin and beam (disc) axes (cf Fig. 15), and

(b) the emitting region is not essentially symmetrical around the spin axis, and

(c) its linear scale is comparable to the neutron star radius, so that a significant fraction of its area, but not the whole of it, is occluded by this star as it spins

Two processes can align the spin and disc (beam) axes. One is warping of the central disc into the neutron-star spin plane through differential torques (both magnetic and precessional), and the other is

direct accretion of angular momentum (characterized by $R_M$) from an unwarped disc lying in the binary orbital plane. The latter seems more likely, as any interruption in the accretion flow means that disc warping has to ‘start again’. Neither process will be very efficient in the Be-star PULXs, where accretion rate is less than $\lesssim 10 M_{Edd}$ and confined to short transient episodes.

The outcomes of the two spin orientations are very different.

* If there is strong misalignment of the spin and beaming axes, a significant part of the pulsed emission can escape without scattering, giving a large pulse fraction (Fig. 15, left). This must be maximal at the highest X-ray energies, as scattering makes the X-rays both softer and less pulsed. This correlation of pulse fraction and X-ray energy is well known for observed PULXs (Kaaret et al., 2017).

* If instead there is substantial alignment of the spin and central disc axes, much of the primary X-ray emission is scattered by the walls of the beaming ‘funnel’ before escaping (Fig. 15, right). Since the light-travel time across the funnel is usually comparable with the pulse duration, the pulse amplitude can be severely reduced. This would result in an inability to detect many highly collimated PULXs through pulse searches (Mushtukov et al., 2019); as data accumulate and improve in quality, this picture will be tested more thoroughly. The pulse fraction can of course also be reduced or entirely removed if enough matter accretes to weaken the neutron-star magnetic field. This presumably happened in the case of Cyg X-2, which is a survivor of a phase of strongly super-Eddington accretion (King and Ritter, 1999). Here the neutron star is not noticeably magnetic, and has probably gained a few tenths of $M_\odot$ during the super-Eddington phase.

These outcomes appear to be consistent with observations of ULXs. Very few ULXs show pulsing. But if most ULXs are the direct descendants of high-mass X-ray binaries (HMXBs, once the companion star fills its Roche lobe), or Be-star HMXBs, the vast (unpulsed) majority must contain neutron stars. First, almost all HMXBs contain neutron stars rather than black holes (King and Lasota, 2016; King et al., 2017). Second, neutron-star systems are more super-Eddington and more beamed than black-hole systems with the same mass transfer rate, and therefore have *higher* apparent ULX luminosities (see Section 3.6). In addition, there is at least one ULX whose spectrum shows a CRSF corresponding to a pulsar-strength magnetic field (Middleton et al., 2021b), but which does not show pulsing, strongly suggesting that is it an aligned accretor.

The vast majority of ULXs do not have detected pulses, despite containing neutron stars. Unless this is caused by a lack of data for ULX pulse searches it implies that alignment of spin and central disc axes is rapid, and possibly that field suppression through accretion occurs also. As expected in this picture, Be-star systems are prominent among the PULXs since accretion rates are relatively low and the episodes of increased luminosity are transient.

<footnote n="10">¹⁰ This unstated assumption has a similar status to the implicit belief that ULX nebulae must be Strömgren spheres powered by the direct radiation of the ULX, rather than wind bubbles made by its powerful near-spherical wind – see Sections 2.6.1 and 3.14. In both cases the effect is to ‘rule out’ beaming.</footnote>
## 3.8. Numerical simulation of super-Eddington accretion flows

Most of the numerical simulations of super-Eddington accretion on to compact bodies concern black holes. The purely absorbing inner boundary conditions make the problem simpler than for neutron stars, which have hard surfaces and, in the most interesting case, strong and rotating magnetic fields. We first review the simulations of black-hole super-Eddington accretion.

### 3.8.1. Numerical simulations of black-hole accretion

Analytical or quasi-analytical solutions of super-Eddington accretion on to black holes represent two extreme possibilities: either the local emission is kept close to the Eddington value by blowing out the “excess” matter (Shakura and Sunyaev, 1973), or by advecting the “excess” energy towards the black hole, where it disappears from the Universe (Abramowicz et al., 1988). In both cases the total luminosity radiated by the accretion flow is $1+\ln \dot{m}_0$ times the Eddington luminosity, where $\dot{m}_0$ corresponds either to a constant accretion rate (slim disc) or is equal to the rate at which matter is brought to the accretion flow, the accretion rate decreasing to keep the local flux at the Eddington value (wind/outflow solution). Although the accretion rate reaching the black hole can differ between the two cases by several orders of magnitude, in neither case can the emitted luminosity be larger than ~10 $L_{Edd}$.<fnref n="11" /> Then for a ULX with $L_X > 10 L_{Edd}$, this apparent luminosity cannot be equal to the physical luminosity, whatever the accretion flow model: the observed super-Eddington luminosity must be beamed. Both extreme models imply beaming depending on the accretion rate.

For slim discs, Wiktorowicz et al. (2017) fitted the relative semi-height of the scattering photosphere calculated analytically and numerically by Lasota et al. (2016) by the formula:

$$ \frac{H_{\text{phot}}}{R} \approx \frac{1.6}{1 + \frac{4}{\dot{m}}}, \qquad (87) $$

and from

$$ \left( \frac{H_{\text{phot}}}{R} \right)^{-1} = \tan \frac{\theta}{2}, \qquad (88) $$

obtained the beaming factor

$$ b = 1 - \cos \frac{\theta}{2}. \qquad (89) $$

This beaming becomes nearly constant for $\dot{m} \gtrsim 100$ which gives the maximum beaming for $b \approx 0.15$. For the windy Shakura and Sunyaev (1973) model with beaming given by Eq. (55), $b = 0.15$ corresponds to $\dot{m} = 22$. One can expect that in reality the solution describing super-Eddington accretion on to a black hole may have both outflows and advection, but the relative proportions can only emerge from numerical calculations and comparison with observations. Unfortunately both these methods encounter serious difficulties.

Before we discuss numerical solutions describing super-Eddington mass supply to black holes, we should note that, without estimates of the black-hole masses in ULXs, these models remain untested, since we do not know the true Eddington luminosities. We now know that the masses of black holes descending from stellar evolution can differ by an order of magnitude (Abbott et al., 2020), so this is a real handicap. Currently there is only one upper limit on the compact-object mass in a ULX (Motch et al., 2014), and here the accretor was later found to be a neutron star (Fürst et al., 2016).

In discussing numerical simulations of super-Eddington accretion flows, we emphasize that these cannot so far represent the conditions expected in ULXs, i.e. a disc of at least $\gtrsim 10^5 R_g$ in size, geometrically thin for $R \gtrsim 15 \dot{m} R_g$. This is not a surprising result: MHD codes

attempting to model the disc viscosity from a first-principles treatment of the magnetorotational instability currently find values of the dimensionless viscosity parameter $\alpha$ significantly smaller ($\lesssim 0.01$) than observational estimates from time variations, which strongly suggest $\alpha \gtrsim 0.3$ (Tetarenko et al., 2018). The reason for this discrepancy is likely to be numerical diffusivity, as already recognized in the solar and MHD-fluids literature. One should add, however, that the viscosity parameter can be reliably determined in the outer disc's regions only (see, e.g., Kotko and Lasota 2012).

Most, if not all, numerical simulations of such flows start with a torus, with an assumed distribution of angular momentum that after relaxation forms an accretion flow. In most cases the “circularization radius” of a torus-generated flow (the Keplerian orbit corresponding to the torus angular momentum) is kept rather small (~30–50 $R_g$). This is because the simulations aim to get quasi-stationary solutions, but this requires run times significantly longer than the viscous time of the accretion flow.

Since the viscous time at radius $R$ is

$$ t_{\text{vis}} \sim \alpha^{-1} \left(\frac{H}{R}\right)^{-2} \left(\frac{R^3}{GM}\right)^{1/2}, \qquad (90) $$

one can write

$$ t_{\text{vis}} \sim 1.1 \times 10^5 \left( \frac{R}{100 R_g} \right)^{3/2} \left( \frac{\alpha}{0.1} \right)^{-1} \left( \frac{H/R}{0.5} \right)^{-2} t_g, \qquad (91) $$

where the dynamical time $t_g \equiv R_g/c$, is often used to express the duration of numerical calculations ($t_g = 5 \times 10^{-5}$ s, for a 10 M$_{\odot}$ black hole). Typical 3D GRRMHD (general-relativistic radiation magnetohydrodynamic) calculations by Sądowski and Narayan (2016) require runs of 20000 $t_g$, so it is not surprising that their quasi-stationary solutions require rather small discs, with $R \lesssim 40 R_g$.

The pioneering simulations of super-Eddington accretion by Eggum et al. (1988) were able to represent only 0.6 s of physical evolution of the system, while the improved scheme of Okuda (2002) increased this duration to 1.6 s, both shorter than the corresponding viscous time. Ohsuga et al. (2005), using 2D radiation hydrodynamic (RHD) simulations were the first to obtain a quasi-steady structure of a supercritical accretion flow on to a black hole. It took another 10 years for other groups to join the Japanese team of Ken Ohsuga in trying to model supercritical accretion on to compact bodies. Sądowski and Narayan (2015, 2016) used a general-relativistic RMHD code first in two, then in three dimensions, while Jiang et al. (2014) performed 3D simulations in the framework of the so-called pseudo-Newtonian potential of Paczynsky and Wiita (1980).<fnref n="12" /> Except for Takahashi et al. (2016), who used the general-relativistic description of the black-hole accretor, other simulations of the Ohsuga group (Ohsuga et al., 2005; Ohsuga and Minesshige, 2011; Hashizume et al., 2015; Kitaki et al., 2018, 2021) used the pseudo-Newtonian description of particle motion in the gravitational field of the accreting body. All these simulations produce rather similar results: strong outflows and regions of the inflow dominated by advection. They differ mostly in details, except for Jiang et al. (2014) who observe an effect not seen by the other teams: vertical advection of radiation caused by magnetic buoyancy that transports energy faster than photon diffusion. This allows a significant fraction of the photons to escape from the surface of the disk before being advected into the black hole; vertical “advection” reduces horizontal advection. The reason for this discrepancy between simulation results is not understood. The main suspect is the treatment of radiation in the code. While Jiang et al. (2014) solve the full radiative diffusion equations, other authors use approximations, such as flux limited diffusion (FLD)

<footnote n="11">¹¹ For “hyper-Eddington” accretion on to an intermediate-mass or super-massive black hole with $\dot{m} = 5000$, the luminosity is still ~30 $L_{Edd}$ (Sakurai et al., 2016).</footnote>

<footnote n="12">¹² Since the Paczyński-Wiita “potential” mimics the general-relativistic description of massive particle orbits in the black-hole’s equatorial plane, it should actually be called “pseudo-relativistic”, but it is too late to change the widely accepted nomenclature. (The name “Paczynski” was misspelled in the original publication, hence the spelling used in the reference list.)</footnote>