© ESO 2019

1. Introduction

Mass transfer rates in binary systems may vary in very broad limits, from very low to amounts vastly exceeding the Eddington limit of the accretor. Highly super-Eddington accretion rates are not surprising in binary systems containing a neutron star (NS) and a massive star filling its Roche lobe. In this case, the mass transfer runs on a relatively short thermal timescale of a massive star, which is hundreds of thousands of years. The mass transfer does not quench because the high mass ratio makes a binary system more likely to remain in contact or even to tighten, thus increasing the mass transfer rate (see e.g. Pavlovskii & Ivanova 2015; Pavlovskii et al. 2017).

The strong magnetic field makes it possible to transport the accreted matter deep into the gravitational well to the NS surface. Presumably, the matter falls to the NS surface in a thin curtain (see Fig. 1a in Basko & Sunyaev 1976). The high ratio of the radiating surface area to the volume, and the reduced scattering cross section in the strong magnetic field are invoked to explain the observed excess above the Eddington limit (see e.g. Mushtukov et al. 2015; Kawashima et al. 2016).

It has recently been realized that some of the ultra-luminous X-ray sources (ULX) in nearby galaxies are actually magnetized NSs accreting well above their Eddington limits. Bachetti et al. (2014) using NuSTAR data discovered coherent pulsations in the ULX X-2 in the galaxy M82. Later, two similar objects, NGC 7793 P13 (Israel et al. 2017a; Fürst et al. 2016) and NGC 5907 X-1 (Israel et al. 2017b), were found. NGC 5907 X-1 is exceptional in its luminosity, exceeding 10⁴¹ erg s⁻¹ during some of the observations. These are the prototypical members of the ULX pulsar (ULXP) family. More recently, other ULXPs have been found. Pulsations with a period of 20–30 s were discovered in NGC 300 ULX1, the supernova impostor SN2010da (Villar et al. 2016), with a very strong spin-up of ṗ = −1.75 × 10⁻⁷ s s⁻¹ (Bachetti et al. 2018). Also, M51 ULX-7 was identified as a ULXP (G. Israel, priv. comm.); and finally, the first ULXP in the Milky Way, Swift J0243.6+6124, was discovered (Kennea et al. 2017; Tsygankov et al. 2018; Wilson-Hodge et al. 2018). Though only a few persistent super-Eddington objects have been robustly identified as NSs to date, it is quite natural to expect more supercritically accreting NSs among the more general class of ULXs (see Kaaret et al. 2017).

Studies of super-Eddington accretion discs started simultaneously with the standard disc theory by Shakura & Sunyaev (1973). They suggested that, for the accretion rate Ṁ₀ exceeding a certain critical value Ṁ_cr, a wind emanates inside the spherization radius R_sph, driven by the radiation pressure, and removes just the right amount of matter to keep the disc at the local Eddington limit. The luminosity of the supercritical disc in this case exceeds the Eddington luminosity L_Edd by a logarithmic factor 1 + ln(Ṁ₀/Ṁ_cr), and the accretion rate decreases with radius as Ṁ(R) = Ṁ₀ R/R_sph. We refer to such a scenario as “classical mass loss”. Implications of this scenario to the discs around magnetized NSs were considered by Lipunov (1982), and recently by Grebenev (2017) and Mushtukov et al. (2019). Remarkably, because the properties of the inner disc in this model are independent of the outer boundary conditions, the size of the magnetosphere and the luminosity at high Ṁ₀ converge to universal values dependent on the NS magnetic moment μ alone.

Later it was realized that it is important to take into account other effects related to deviations from the thin disc approximation. The most significant departure from the standard disc model is the heat advection, which violates the locality assumption in the energy balance of the disc. The role of advection at low accretion rates was considered, for example, by Ichimaru (1977), Narayan & Yi (1994, 1995), and Abramowicz et al. (1995, 1996), and at high accretion rates by Begelman & Meier (1982), Abramowicz et al. (1988), and Beloborodov (1998). Numerical simulations of supercritical accretion onto a black hole were performed by Eggum et al. (1988), Ohsuga et al. (2005), Okuda et al. (2005), and more recently by Ohsuga & Mineshige (2011), McKinney et al. (2014), Sa̧dowski et al. (2014), Sa̧dowski & Narayan (2016), and Ogawa et al. (2017). Semi-analytical models for supercritical discs including heat advection and outflows were constructed by Lipunova (1999), Kitabatake et al. (2002), Fukue (2004), and Poutanen et al. (2007).

While most of this theoretical work on supercritical accretion was devoted to discs around black holes, the discovery of ULXPs has drawn attention to magnetized NSs accreting at high rates. Numerical magnetohydrodynamic (MHD) simulations of subcritical disc accretion onto magnetized neutron stars have been carried out in Romanova et al. (2002, 2003), Long et al. (2005), who studied extensively the properties of the magnetospheric MHD flow and the processes at the inner boundary of the disc, and the formation of funnel-shaped flows in particular. They consider magnetically driven winds, which are particularly favoured in the propeller regime (see review by Romanova & Owocki 2015). In our work we focus instead on radiatively driven outflows, which are expected to be launched by supercritical discs. It is important to note that even in highly supercritical ULXPs like M82 X-2 and NGC 7793 P13, the accretion disc outside the magnetosphere may remain in a subcritical regime. Chashkina et al. (2017; hereafter CAP) considered accretion onto a magnetized NS in a regime that can be characterized as intermediate, with an accretion rate in the range (10⁻⁹−10⁻⁶) M_⊙ yr⁻¹ for μ ∼ 10³¹−10³² G cm³ inferred for ULXPs (see Tsygankov et al. 2016). The disc outside the magnetosphere in this case remains geometrically thin, nearly Keplerian, and does not lose considerable amounts of matter. This allowed us to use certain results of the standard disc theory such as local radiation energy balance and simple conservation laws for mass and angular momentum. This approach works well if the formally calculated spherization radius is smaller than the Alfvén radius. However, super-Eddington NSs with lower magnetic moments or those accreting at higher rates should possess a supercritical advective disc, whose radiation pressure becomes sufficient to unbind part of the accreting material. Recently, simulations of supercritical accretion onto a non-magnetized NS have been performed by Abarca et al. (2018) in an attempt to explain non-pulsating ULXs.

In the present work, we develop a model of a supercritical accretion disc around a magnetized NS, applicable to a broader class of objects, including NSs with pulsar-scale magnetic fields, luminosities hundreds of times higher than the Eddington limit, and in particular extremely bright sources like NGC 5907 X-1. The structure of the paper is as follows. In Sect. 2 we describe our model and present the main equations describing disc accretion accounting for the effects of advection and the mass and momentum loss in a wind. Section 3 is devoted to the results; we present the dependences of the magnetospheric radius on the mass accretion rate and NS magnetic moment. We also consider the effects of the pulsar spin and the irradiation of the disc by the central source on the disc structure. We discuss our results and apply them to particular ULXPs in Sect. 4. We conclude in Sect. 5.

2. Model

2.1. Basic equations

We consider a NS of mass M, radius R_NS, and magnetic moment μ. We assume it is an aligned rotator with an unperturbed dipolar magnetic field inside the magnetospheric radius, which is equal to the inner radius of the disc R_in. Outside this radius, the magnetic field lines are opened by the ideally conducting accretion disc, as proposed by Parfrey et al. (2016). A simple sketch of the adopted structure of the disc is shown in Fig. 1. All the interactions between the disc and the magnetosphere are assumed to occur in a narrow strip at the edge of the disc and are reduced to two boundary conditions introduced in CAP and later in this section. One of the conditions, the pressure balance, as shown in 3D MHD simulations by Kulkarni & Romanova (2013) for lower accretion rates, corresponds to the radius from which the matter starts to flow along a funnel stream.

Fig. 1. Structure of an accretion disc around a ULXP. For very high mass accretion rates, the inner parts of the disc (inside the spherization radius Rsph shown here) enter the super-Eddington accretion regime. Inside Rsph the thin disc model is not applicable. Mass loss in a wind is shown by blue arrows. The accretion column, where most of the energy is released, is shown by yellow cones, and the red wavy lines refer to the radiation of the column that may affect the inner disc pressure balance. The red vertical line marks the effective boundary between the disc and the magnetosphere at Rin.

In our model, we include three effects important for an accretion disc at a near- or supercritical accretion rate: (i) advection, (ii) mass loss in the wind, and (iii) angular momentum loss in the wind. Advection is related to the increasing photon diffusion timescales at high mass accretion rates. It alters the energy equation. The matter loss in the wind makes the disc accretion rate dependent on radius, and also affects the angular momentum and the energy conservation equations because the wind carries some angular momentum and energy. Poutanen et al. (2007) assumed that the specific angular momentum of the wind is equal to that of the matter in the disc. However, if the wind is magnetized, it can efficiently remove the angular momentum from the disc, as in the centrifugally driven wind model of Blandford & Payne (1982).

We follow the basic framework of CAP and retain some of the important assumptions of the model. First of all, the boundary conditions at the inner disc edge remain the same. We take the torque balance in the form

(1)

and assume pressure balance, which yields for the viscous stress tensor the following relation:

(2)

Here, Ṁ_in is the mass accretion rate at the inner edge of the disc; Ω_NS is the angular velocity of the NS; Ω_in, H_in, and are the angular velocity, the half-thickness, and the vertically integrated rϕ-component of the viscous stress tensor at the inner boundary of the disc R_in; the dimensionless constant k_t parametrizes the efficiency of the angular momentum removal by the magnetic and viscous torques; and α is dimensionless viscosity parameter as defined by Shakura & Sunyaev (1973). The luminosity L = ηṀ_in c² is released close to the NS in the accretion column with efficiency η as measured by an imaginary observer in the equatorial plane at the inner edge of the disc. Because the accretion column emits anisotropically, the efficiency η may differ from the angle-averaged efficiency of an isotropic source and from the efficiency measured by an observer at infinity. The radial structure of the disc is described by the angular frequency Ω, thickness H, and the viscous stress W_rϕ. The angular frequency and the viscous stress should conform to the boundary conditions (1) and (2). Equation (1) is similar to the inner boundary condition in Spruit & Taam (1993) and Rappaport et al. (2004) in the case of a slowly rotating NS. The additional term in the RHS is radiation drag that can be important at high luminosities L >  L_Edd. We note that Eq. (1) is valid only for the case of the disc rotating faster than the magnetosphere. For the propeller case (Ω_in <  Ω_NS), the boundary condition would be different. Moreover, one could consider a radial distribution of the magnetic torque applied to the disc, as was done, for example, by Kluzniak & Rappaport (2007). On the other hand, Matt & Pudritz (2005b) argue that if Ω_in >  Ω_NS, the magnetic interaction can be limited to a narrow ring.

We take the radial component of the momentum equation in the form

(3)

where is the surface density and is the vertically integrated pressure. As in CAP, we ignore the dynamic term v_R dv_R/dR, which is suppressed by a factor of α²(H/R)² ≪ 1 with respect to the pressure gradient term.

In all the equations for the radial structure we use vertically integrated quantities. In the equations for the vertical structure we assume a fixed vertical effective polytropic index n, defined through the equation for the density ρ = ρ_c[1 − (z/H)²]ⁿ and for the pressure P = P_c[1 − (z/H)²]^n + 1 (Paczynski & Jaroszynski 1978). This allows us to compute analytically the connection factors between the midplane quantities such as ρ_c and P_c, and the corresponding vertically integrated values Σ and Π, as described in Appendix A. For the vertical temperature profile used in the advection term calculations we use the expression T  =  T_c[1 − (z/H)²]^{(n + 1)/4}, which is valid for the radiation-pressure-dominated case as T ∝ p^1/4; here T_c is the central temperature on the disc.

The angular momentum conservation equation is modified by an additional term corresponding to the angular momentum outflow in the wind,

(4)

where ψ ≥ 1 allows us to scale up the net angular momentum lost in the wind. If ψ = 1, the net angular momentum in the wind is equal to that in the disc (Lipunova 1999; Poutanen et al. 2007). Larger ψ can appear in magneto-centrifugal winds, where the sub-Alfvénic part of the transonic flow rotates approximately rigidly, thus increasing the angular momentum loss by a factor of ∼(R_wA/R)², where R_wA is the cylindrical radius of the Alfvén surface, and R is the launch point of the wind streamline (see also Matt & Pudritz 2005a who used this approach for magnetospheric outflows). In the centrifugally driven wind model of Blandford & Payne (1982), for instance, R_wA/R may be as high as 5. Simulations of the centrifugally driven winds (e.g. Ustyugova et al. 1999) suggest even greater lever lengths at low accretion rates, while at high accretion rates when the wind is radiatively driven the effect of magnetic stresses is smaller (see e.g. Proga 2000).

The mass accretion rate derivative dṀ/dR needs to be calculated using some additional equations describing the physics of the wind launching. We assume that some fraction ϵ_w ≤ 1 of the energy leaving the disc with radiation is spent to accelerate the outflow (Lipunova 1999; Poutanen et al. 2007),

(5)

where Ω_K is Keplerian angular velocity and T_eff is effective temperature. The model of optically thick, energy-driven wind developed by Poutanen et al. (2007) gives similar results. Physically, there could be wind losses everywhere in the disc, but we take into account only the continual radiation-driven wind which would operate within the spherization radius R_sph defined by the condition on the minimal relative thickness H/R >  (H/R)_cr. Because the details of wind launching are uncertain, we assume (H/R)_cr = 1. We note that for a disc around a black hole, the wind mass loss is determined only by the mass of the black hole and the outer accretion rate. For a disc around a magnetized neutron star, the mass lost in the wind also depends on the magnetic moment of the star and its spin period.

We use the α-viscosity prescription for the vertically integrated viscous stress:

(6)

The midplane pressure can be expressed as a sum of the radiation and gas pressures,

(7)

where m̃ is the mean particle mass and (Eq. (A.5)). Here we have taken into account the α-viscosity prescription and the relations between the central and vertically integrated quantities implied by the adopted vertical structure (Eqs. (A.3) and (A.4)). We do not use Eq. (7) directly, but as a supplement to the energy equation (see Sects. 2.2 and 2.3). Treatment of advection is a new part of the model that requires a separate consideration.

2.2. Advection

In a general case with advection, some fraction of energy is radially transported. The energy flux carried by radiation diffusion in the vertical direction is no longer equal to the local energy release,

(8)

where the total energy released in the disc at radius R per unit area is

(9)

and the radiation flux from both sides of the disc is (see Eq. (A.10))

(10)

Here, κ ≃ 0.34 cm² g⁻¹ is the Thomson scattering opacity. The advected flux Q_adv can be viewed (as done in e.g. Lipunova 1999) as the flux of heat carried with the flow, and thus can be expressed through the specific entropy per particle s

(11)

where the radial velocity (taken by its absolute value) is v_R = Ṁ/(2πRΣ). Taking into account Eq. (7) we can expand the expression above as (see Appendix B)

(12)

where the coefficients 𝒮, 𝒫, 𝒬, and ℛ are given by Eqs. (B.8)–(B.11). The energy balance can be written using Eqs. (8), (9), and (12) as

(13)

This relation can be re-written for the derivative of the angular velocity as

(14)

where the coefficients C_{Ω, Σ, w, T, free} are given by expressions (B.12)–(B.16).

2.3. Solving the disc equations

Unlike in our previous work (CAP), here there are five variables for which we need to solve differential equations: the angular frequency Ω, the tangential stress W_rϕ, the surface density Σ, the mass accretion rate Ṁ, and the midplane temperature T_c. There are also seven free global parameters: the viscosity parameter α, the wind efficiency parameter ϵ_w, the wind magnetization parameter ψ, the polytropic index n, the spin period p and magnetic moment of a NS μ, and the accretion efficiency η. We calculate the disc structure from its outer edge to the inner boundary.

We solve the five equations listed below:

1. The derivative of vertically integrated viscous stress that can be obtained from the radial Euler Eq. (3) and α-viscosity prescription (6):

   (15)

2. Whenever the outflow condition (H/R >  (H/R)_cr) is satisfied, the mass accretion rate changes with radius according to Eq. (5):

   (16)

3. The differential equation for the central temperature can be obtained by taking derivative of the pressure Eq. (7):

   (17)

   Here, we used the thickness of the disc H following from the hydrostatic equilibrium (see Eq. (A.2)) together with the gradient of W_rϕ substituted from Eq. (14). We also use the ratio of the gas pressure to the total pressure in the equatorial plane expressed using the adopted vertical structure (see Eq. (A.11)):

   (18)

4. To determine the radial dependence of the angular velocity, we substitute Eqs. (15) and (16) into angular momentum conservation Eq. (4):

   (19)

   The second term on the right-hand side is switched off for a thin disc (H/R <  (H/R)_cr) when there are no outflows.

5. The differential equation for the surface density was obtained from the advection Eq. (14) substituting all other derivatives from Eqs. (15)–(19):

   (20)

The disc structure equations are further converted to a more compact dimensionless form in Appendix C.

Our solution should satisfy a number of boundary conditions. The accretion rate at the outer boundary of the integration region should be equal to Ṁ₀, which is one of the global parameters of the simulation. Also two variables, Ω and W_rϕ, should satisfy boundary conditions (1) and (2). In order to achieve that we use two adjustable parameters. The first one is the tangential stress at the outer boundary . The second one is the magnetospheric radius R_in, which we parametrize through the ratio

(21)

where

(22)

is the Alfvén radius (see e.g. Elsner & Lamb 1977 and Sect. 6.3 in Frank et al. 2002). As a result of the iteration procedure we get the five sought for functions and the value for ξ. Because the accretion rate at the inner boundary of the disc is physically more relevant, we also use a differently normalized version of the relative magnetospheric radius:

(23)

3. Results

3.1. Global parameters

The relative magnetospheric size ξ is found as one of the two eigenvalues of the problem, the other being . Other global parameters obtained in the model include several quantities at the magnetospheric boundary: the relative thickness of the disc (H/R)_in, the fraction of mass reaching the magnetosphere Ṁ_in/Ṁ₀, and the advection fraction (Q_adv/Q⁺)_in. We also track the maximum thickness of the disc (H/R)_max.

All the simulations were made for the following values of the parameters: α = 0.1, k_t = 0.5, ϵ_w = 0.5, ψ = 1, n = 1, and η = 0 if not stated otherwise. We first consider a slowly rotating accretor with the spin period equal to ten equilibrium periods p = 10p_eq,¹ where the equilibrium period (Lipunov 1992; Illarionov & Sunyaev 1975) is defined as the spin period for which the Alfvén radius R_A equals the corotation radius ,

(24)

where μ₃₀ = μ/10³⁰ G cm³ is the dimensionless NS magnetic moment and m = M/1.4M_⊙ is the normalized NS mass. We also use the dimensionless mass accretion rate ṁ₀ = Ṁ₀/Ṁ_Edd normalized by the Eddington value Ṁ_Edd = 4πGM/cκ. In Fig. 2 we show the contours of ξ_eff in the ṁ₀−μ plane, covering two orders of magnitude in magnetic moment and five orders of magnitude in the mass accretion rate.

Fig. 2. Relative magnetospheric radius ξeff given by Eq. (23) is colour-coded and shown by contours on the ṁ0− μ plane. The spin period is p = 10 peq. The luminosity scale is calculated assuming efficiency η = 0.1, but irradiation effects on the disc structure were ignored.

The magnetospheric radius and the thickness of the disc are tightly related. There is good agreement with the results of CAP, as we can see from Fig. 3. The relative magnetospheric size ξ_eff remains a monotonic function of (H/R)_in and behaves in approximate accordance with Eq. (57) in CAP (within 5%). This long-period approximation works fine far from the equilibrium period, for instance, in outbursts.

Fig. 3. Relative magnetospheric radii ξ (dashed blue curve) and ξeff (green solid) as functions of the disc thickness at the magnetospheric boundary. The dotted red line corresponds to the long-period asymptotic given by Eq. (57) in CAP. Magnetic moment of the NS was set to μ30 = 100.

Not only mass but also angular momentum is lost in the wind. If there is no angular momentum flow due to stresses in the wind (ψ = 1), the fraction of angular momentum expelled from the disc depends only on the distribution of the mass loss over the radial coordinate. The total angular momentum flux through an annulus in the disc is composed of the angular momentum carried by the matter in the disc and of the viscous torque acting on the annulus,

(25)

where ω = Ω/Ω_K. Unless some of this angular momentum is removed in the wind, l_in = l(R_in) should be equal to l₀ = l(R_out). In Fig. 4 we show the ratio of the angular momentum fluxes l_in/l₀ as a function of the ratio of mass accretion rates Ṁ_in/Ṁ₀. Both ratios start at unity for conservative thin disc accretion and then decrease as Ṁ₀ increases. The slope of the curve in Fig. 4 shows the evolution of the mean net angular momentum in the wind. With increasing mass accretion rate, the outflow involves larger radii. Characteristic radii losing most of the angular momentum is approximately equal to the spherization radius that leads to the scaling (see Fig. 4):

(26)

Fig. 4. Fraction of the angular momentum flux retained in the super-Eddington disc lin/l0 as a function of the fraction of mass reaching the magnetosphere. The dotted black line corresponds to equal ratios (lin/l0 = Ṁin/Ṁ0, as expected if the net angular momentum is constant); the solid green line corresponds to the expected scaling of a growing spherization radius (Eq. (26)).

The amount of angular momentum lost in a centrifugal wind is enhanced approximately proportionally to ψ.

3.2. Dependence on the spin period

The transition to the propeller regime, when the accretion flow can no longer spin up the NS, can be traced using the fastness parameter ω_s defined as

(27)

When ω_s = 1, the inner rim of the disc rotates exactly with the same frequency as the magnetosphere, making our first boundary condition marginally satisfiable.

The main effect of increasing the fastness parameter (or decreasing the spin period of the accretor) on the properties of our solution is the increasing ratio ξ of the size of the magnetosphere to the Alfvén radius. The factor ξ_eff increases by about 40% between the slowly rotating NS case and the propeller limit (see Fig. 5).

Fig. 5. Parameter ξeff as a function of fastness parameter for magnetic dipole moments μ30 = 1 and μ30 = 10 and mass accretion rates ṁ0 = 104 and ṁ0 = 5 × 103.

Near the corotation, when R_co = R_in, the disc thickness approaches zero and the boundary condition for the viscous stresses (2) reduces to the zero-torque condition used in standard disc theory. This allows us to compare our results directly to some of the results obtained using the codes designed for black hole accretion by Poutanen et al. (2007), among others. Their model takes into account the advection effects and the outflows from the disc. We compare the spherization radius defined as the maximum distance from the NS where the condition for outflows H >  R is fulfilled.

In Fig. 6, we show the spherization radius normalized by the dimensionless mass accretion rate (green dots) and compare our results to the zero-torque case. For ϵ_w = 0.5 and ṁ₀ = 10⁴, Eq. (21) in Poutanen et al. (2007) predicts R_sph ≃ 0.575 ṁ₀(GM/c²), in reasonable agreement with our results close to equilibrium.

Fig. 6. Ratio of the spherization radius rsph (in units of gravitational radius Rg) to the dimensionless mass accretion rate ṁ0 = 104 as a function of the fastness parameter (green dots). The red horizontal line gives the result from Eq. (21) in Poutanen et al. (2007).

3.3. Effects of irradiation

All the previous results were calculated without irradiation effects by setting η = 0 (or L = 0 in Eqs. (1) and (2)). The real efficiency affecting the pressure balance is probably of the same order as the integrated accretion efficiency, though strong anisotropy of the radiation from the column is not excluded. Confirming the result of CAP, we find that radiation from the accretion column is an important factor affecting the structure of the disc, in particular, the radius of the magnetosphere. The magnetospheric radius increases by up to 30% for ṁ₀ ∼ 10⁴, as is shown in Fig. 7. The effect grows rapidly with the mass accretion rate and with the magnetic moment. The latter is a consequence of a rapidly growing ratio P_rad/P_mag with the disc inner radius in the dipole approximation (see e.g. CAP, Eq. (63)).

Fig. 7. Colour-coded contours of relative correction to the magnetospheric radius caused by irradiation ξ(η = 0.1)/ξ(η = 0)−1 shown in the ṁ0 − μ plane.

3.4. Dependence on other parameters

There are adjustable parameters that influence the structure of the disc. In Table 1 we show the global properties of the discs for different parameters. The fiducial model here is the model with α = 0.1, ϵ_w = 0.5, n = 1, η = 0, and ψ = 1. All the models are calculated for a NS with magnetic field μ = 10³⁰ G cm³ and accretion rate ṁ₀ = 3000. Other models differ from the fiducial model by a single parameter.

As in CAP, we find the viscosity parameter to be an important factor that alters the structure of the disc. For the parameters under consideration, changing α from 0.1 to 0.5 results in a twofold decrease in disc thickness and ξ, and quenches wind formation. Changing the vertical disc structure by increasing the effective polytropic index n also makes the disc thinner and less likely to form outflows.

3.5. Effects of advection and wind

Under the assumptions we use, including the adopted vertical structure and the minimum disc thickness for wind launching, advection starts to play a role rather early, when the entire disc is still subcritical. As a consequence, the disc thickness stabilizes at H ≃ R (see Fig. 8). A slimming effect of advection was noted earlier by Abramowicz et al. (1988), Beloborodov (1998), Lipunova (1999), and Lasota et al. (2016), among others. The local Eddington limit (H = R) is reached at the critical mass accretion rate of

(28)

i.e. 5–6 times higher than in CAP (see their Eq. (66)).

Advection starts to dominate in the energy balance already below this limit, at ṁ₀ ∼ 10³ (see Fig. 9), making the inner disc a huge reservoir of heat. For ṁ₀ above the limit given by Eq. (28), most of the gravitational energy released in the disc is stored as heat. To illustrate this, we calculated the cumulative luminosity of the disc integrated from some radius R to the outer radius . Similarly, by integrating Q_rad and Q_adv, we can define the cumulative radiative and advection powers. The total cumulative power is shown by a blue dashed line in Fig. 10; at the inner radius R_in it is in full agreement with theoretical prediction, L_theor = GMṀ/2R_in. We note that the ratio of the disc luminosity to the luminosity of the accretion column can be as low as ∼R_NS/R_in.

Fig. 8. Relative disc thickness as a function of the radius for a model with μ = 1030 G cm3, ṁ0 = 3 × 103, and p = 0.67 s. Our results are shown by the red dotted curve, whereas the solid green and blue curves correspond to the asymptotics for the zones B and A of the standard disc, respectively.

Fig. 9. Advection effects for the model shown in Fig. 8. The advection flux fraction Qadv/Q+ and the radiated fraction Qrad/Q+ are shown by the red solid and green dashed curves, respectively.

Fig. 10. Cumulative luminosities as functions of radius. The total power generated in the disc is shown by the dashed blue curve. The solid red and dotted black curves show the radiative and advected luminosities, respectively. The parameters of the model are μ30 = 1 and ṁ0 = 3000.

In Fig. 11, we plot the fraction of the initial mass accretion rate remaining in the disc for different ṁ₀ as a function of radius. At very high accretion rates, the supercritical wind blows away a considerable part of the accretion material and operates at all radii within R_sph. At some intermediate rates, 2000 ≲ ṁ₀ ≲ 3000, there is a prominent subcritical region near R_in where there is no wind. This is caused by the non-monotonic dependence of the disc height on radius. We note that the amount of the blown-away material depends on the condition used to switch on the wind. The wind launching relies on complex physical processes and can be described only very approximately in 1D models.

Fig. 11. Fraction of the mass accretion rate reaching radius R for a NS with magnetic moment μ30 = 1. The lines from top to bottom correspond to different ṁ0 in the interval 2000–3500 with steps of 500.

3.6. Magnetospheric radius for different accretion regimes

One of the most important outputs of our model is the radius of the magnetosphere. In many standard models (Ghosh et al. 1977; Koenigl 1991; Wang 1996; Kluzniak & Rappaport 2007), the magnetospheric radius is supposed to scale with the Alfvén radius, corresponding to the constant ξ in our notation. It is worth noticing that MHD calculations of magnetospheric flows with relatively weak accretion rates (Bessolaz et al. 2008; Kulkarni & Romanova 2013) demonstrate similar values of ξ ∼ 0.4 to those in Fig. 2. Chashkina et al. (2017) have demonstrated that the radius of the magnetosphere interacting with a thin radiation-pressure-dominated disc depends weakly on the mass accretion rate. Accounting for the effects of advection and wind losses makes the picture more complicated.

The dependence of the magnetospheric radius (in units of gravitational radius R_g = GM/c²) on the mass accretion rate is shown in Fig. 12 for a wide range of accretion rates. The inner disc regions of most X-ray pulsars are in the gas-pressure-dominated regime. As the accretion rate increases, radiation pressure becomes important. For pulsar-scale magnetic fields, μ ∼ 10³⁰ G cm³, this happens at luminosities of a few times L_Edd, which are quite reachable, for instance, in Be/X-ray binaries during strong outbursts like the recent super-Eddington outburst of SMC X-3 (Townsend et al. 2017; Tsygankov et al. 2017; Weng et al. 2017). As was shown in CAP (see Eq. (61)), the magnetospheric radius becomes almost independent of the accretion rate, if the radiation pressure dominates at the inner edge of a subcritical disc:

(29)

Fig. 12. Magnetospheric radius in units of Rg as a function of the accretion rate for a NS with magnetic moment μ30 = 1. Parts of the black solid curve with different slopes correspond to the different regimes of accretion near the magnetospheric boundary. Two standard solutions () are plotted with the grey dotted lines: ξ = 0.5 and ξ = 1 (spherically symmetric case).

Thus, provided there is a direct measurement of the magnetospheric radius, for example, from quasi-periodic oscillations, we can directly estimate the magnetic field of a NS, with a weak dependence on the viscosity parameter α.

The inner disc radius is defined mainly by the balance of pressures. The pressure inside the disc is related to its thickness. Thus, the dependence of H/R on ṁ₀ is crucial for the behaviour of r_in(ṁ₀). As the accretion rate increases, advection starts to play an important role. The relative thickness of the disc is no longer proportional to ṁ₀, and the magnetospheric radius again depends on ṁ₀. The interplay between wind losses and advection makes the radius dependence on mass accretion rate shallower than the ξ = const approximation historically proposed for spherical accretion, but stronger than r_in = const.

With the increasing magnetic field, the magnetospheric radius increases, and all the boundaries between different regimes shift to higher accretion rates, as shown in Fig. 13. In addition, the length of the plateau corresponding to the thin radiation-pressure-dominated inner disc gradually decreases with the magnetic field, becoming effectively zero at μ ∼ 10³² G cm³. Weakly magnetized objects, on the other hand, should have a prominent region of constant magnetosphere size. The plateau starts when radiation pressure begins to dominate over gas pressure at the radius of the magnetosphere. Position of the boundary between gas- and radiation-pressure-dominated regions of the standard disc scales with mass accretion rate as R_ab ∝ ṁ^16/21 (Shakura & Sunyaev 1973). Because R_in ∝ μ^4/9 in the radiation-pressure-dominated regime (see Eq. (29)), the left boundary of the plateau depends on the magnetic moment as ṁ_in,left ∝ μ^7/12. The right boundary of the plateau is determined by advection effects. Advection becomes important when Q_adv ∼ Q_rad, which implies H_in ∼ R_in. Hence the radius at which the inner disc becomes advective scales linearly with the mass accretion rate, and the accretion rate at the right end of the plateau is ṁ_in,right ∝ μ^4/9. Thus, the length of the plateau slowly decreases with magnetic moment as ṁ_in,right/ṁ_in,left ∝ μ^−5/36.

Fig. 13. Same as Fig. 12, but for different magnetic moments μ = 1030−1032 G cm3.

It is interesting to compare our results with the classical prescriptions (Ghosh et al. 1977; Wang 1996). Figure 12 demonstrates that the classical dependences are much steeper, having the slope of δ ≡ d log r_in/d log ṁ = −2/7. In the case of μ = 10³⁰ G cm³, the slope is δ ≈ −0.21 in the gas-pressure dominated case, δ ≈ −0.07 in the radiation-pressure-dominated regime, and δ ≈ −0.16 when advection dominates. The evolution of the local slope δ is plotted in Fig. 14. At low accretion rates, corresponding to the gas-pressure-dominated thin disc, we get δ ∼ −0.23 in accordance with the results of CAP (see their Sect. 5.1). The maximum value of δ depends on how prominent the thin radiation-pressure-dominated part of the disc is for a given magnetic field, changing from nearly zero for small μ to about −0.1 for magnetar-scale fields. The highest mass accretion rates tend to reproduce much steeper dependences, approaching δ ≃ −2/7 ≃ −0.29. We note that, for a comparatively narrow range of parameters (relevant for T Tauri star magnetospheres), Kulkarni & Romanova (2013) have obtained δ ≃ −0.2 from their 3D MHD simulations.

Fig. 14. Slope of the dependence of the magnetospheric radius on the accretion rate for the models shown in Fig. 13.

4. Application to ULXPs

4.1. NGC 5907 X-1 as a ULXP with a supercritical accretion disc

The ULXP NGC 5907 X-1 has a huge period derivative, even after averaging in time: its period has changed from 1.43 to 1.13 s over the 10 years of observations (Israel et al. 2017b). It is also remarkable that its luminosity exceeds 10⁴¹ erg s⁻¹ during some of the observations. The maximum detected period derivative (by absolute value) was ṗ = −5 × 10⁻⁹ s s⁻¹, about an order of magnitude higher than the average value for this object. The high value of |p̄| was put forward as an argument for this object being in a pure spin-up state, with the unconstrained braking term in the angular momentum equation being negligible.

The rate of change of the total angular momentum of a NS can be written as

(30)

where I ≃ (1−2) × 10⁴⁵ g cm² is the moment of inertia of the NS, and and K_sd are spin-up and spin-down torques. Ignoring the unknown spin-down contribution provides the constraint

(31)

resulting in an inequality for accretion rate

(32)

where I₄₅ = I/10⁴⁵ g cm², and ṗ₋₁₂ = ṗ/(10⁻¹² s s⁻¹). Alternatively, having an independent estimate for the accretion rate, we then can set a lower limit for the magnetic field of this object. For a maximum accretion rate of ∼6000, corresponding to isotropic X-ray emission, the normalized magnetic moment μ₃₀ cannot be less than ∼0.06.

In the opposite case, when a NS is close to equilibrium, its magnetospheric radius is about the size of the corotation radius. Thus, we can put an upper limit on the magnetic field, which suggests that the NS is still in the accretion regime. In terms of mass accretion rate, this condition may be written as

(33)

Figure 15 shows these two limits as applied to NGC 5907 X-1. The spin-up line shows the lower limit for Ṁ₀ set by inequality (32) using period derivative ṗ = −5 × 10⁻⁹ s s⁻¹. The area below this line is forbidden unless some additional spin-up process is present. The propeller line shows the boundary of the region where the magnetospheric radius is equal to the corotation radius. The region below this line is prohibited because the accretion disc rotates more slowly than the magnetosphere of the NS, and no stable accretion is possible. The maximum bolometric luminosity of this object is about 2 × 10⁴¹ erg s⁻¹ (Israel et al. 2017b), which leads to a lower limit on the mass accretion rate ṁ₀ > ṁ_in ∼ 6500 assuming efficiency η ∼ 0.15. This value of efficiency does not take into account any beaming effects that in principle can alter the observed efficiency value. This figure provides evidence for a truly high mass accretion rate in NGC 5907 X-1, as significantly low ṁ₀ ≲ 10³ are forbidden for any magnetic moments. Therefore, any beaming exceeding a factor of 5–6 is unlikely because with increasing beaming, the red line in Fig. 15 will move down and there will be no allowed region left. Another argument against strong beaming is the observed high pulse fractions and nearly sinusoidal profiles, inconsistent with collimation by a wind, and indicating that the emission site is seen during a large fraction of the spin period.

Fig. 15. Restrictions for the position of NGC 5907 X-1 in the ṁ0 − μ plane. The red horizontal line corresponds to the mass accretion rate of ṁ0 = 6500 (accretion efficiency η ≃ 0.15). The black solid line corresponds to peq = p = 1.137 s. The dotted spin-up line shows the lower limit for ṁin set by inequality (32) using period derivation |ṗ| = 5 × 10−9 s s−1. We use here M = 1.4 M⊙, R = 10 km, I = 1.5 × 1045 g cm2, and p = 1.137 s.

We can set an upper limit for the magnetic field of NGC 5907 X-1 as μ ≤ 7.5 × 10³¹ G cm³ and μ ≤ 5.45 × 10³¹ G cm³ if we take into account irradiation from the column (see Sect. 3.3). There is evidence for bimodal distribution in luminosities of NGC 5907 X-1 (see Fig. S2 in Israel et al. 2017b) that can be interpreted as a manifestation of the propeller effect (similarly to M82 X-2 in Tsygankov et al. 2016). This would mean that the source is close to the propeller line. Then instead of an upper limit on the magnetic moment, we get its accurate estimate. Beaming does not play a major role as magnetic field weakly depends on ṁ₀ (as , see Eq. (33)). The disc is expected to be supercritical (i.e. having winds) in its inner parts if μ₃₀ ≲ 14 (see Eq. (28)). At a pulsar magnetic field μ₃₀ ∼ 1, the expected outflow rate from the disc is about 20% of the inflow rate, or 5 × 10⁻⁶ M_⊙ yr⁻¹.

4.2. NGC 300 ULX1: a ULXP with a strong spin-up

The source NGC 300 ULX1 was identified as a ULXP only recently by Carpano et al. (2018). It is characterized by a moderate peak luminosity of L ≃ 4.7 × 10³⁹ erg s⁻¹, but the observed spin-up rate of this source is exceptionally high: its spin period has changed from 126 to 19 s in about 4 years and is still decreasing (see Vasilopoulos et al. 2018). Its highest spin-up rate is estimated as ṗ ≃ −5.5 × 10⁻⁷ s s⁻¹, which is the highest ṗ ever observed from an accreting NS. The long period makes NGC 300 ULX1 a promising candidate for a NS in a pure spin-up state. This can be seen in Fig. 16, where the observed properties of the object are in good agreement with a pulsar-scale μ ∼ 10³⁰ G cm³ of a NS that is currently being rapidly spun up. As in the case of NGC 5907 X-1 considered above, pure spin-up gives us a lower limit on the magnetic field while the propeller limit gives us an upper limit. Thus, the magnetic field of NGC 300 ULX1 is in the range μ = (1.5−85) × 10³⁰ G cm³, which is not affected much by the effects of irradiation. The range is very wide, but sets an upper limit for the beaming factor of 2.5 only, consistent with the constraints set by Binder et al. (2018) based on the number of He II-ionizing quanta. The lower boundary here corresponds to pure spin-up and implies an equilibrium period of about p_eq ∼ 0.2 s. To reach such a period, the object needs to accrete at the observed rate for about t ∼ p_s/ṗ ∼ 1 yr.

Fig. 16. Restrictions for the position of NGC 300 ULX1 in the ṁ0 − μ plane. The red horizontal line corresponds to the luminosity of 4.7 × 1039 erg s−1 assuming accretion efficiency η ≃ 0.15. The black solid line corresponds to peq = p = 20 s. The dotted spin-up line shows the lower limit for ṁ0 set by inequality (32) using period derivation |ṗ| = 5.5 × 10−7 s s−1 (Carpano et al. 2018). We use here M = 1.4 M⊙, R = 10 km, I = 1.5 × 1045 g cm2, and p = 20 s.

At the moment of the first observations of NGC 300 ULX1 in 2010, the NS should have had a period even longer than the 45 s measured in 2016. The question then arises why the NS had such a long period to start with. This period may have originated from a long episode of a very low rate of wind accretion if we assume that this NS was born rotating much faster. To estimate the minimum timescale of the spin-down to such a long period, we use Eq. (30), neglect any spin-up torque acting on the NS, and parametrize the spin-down torque as (Lipunov 1992). The time required to spin down to the observed period is about years. This means that, in order to explain a recently observed long spin period, the system with a normal pulsar-scale magnetic field should have spent tens of thousands of years accreting at a very low rate, constantly in propeller regime. The very high luminosity that we observe now is an exceptionally rare event in this system. This is consistent with the extraordinary properties of the optical B[e] transient SN2010da; this object is associated with the disappearance of a huge amount of hot dust observed before the outburst, very high peak optical luminosity, and a bright B[e] supergiant observed after the event (Villar et al. 2016). All this fits well with a catastrophic event in a previously faint Be/X-ray binary.

5. Conclusions

Our model provides a simple, physically motivated description of accretion onto a NS with a magnetosphere, where the interaction between the accretion disc and the magnetosphere of the NS is reduced to two boundary conditions. This allows us to reconstruct the structure of the disc and obtain a relative magnetospheric radius ξ, which is important for describing the spin evolution of magnetized NSs and for interpreting observational data on X-ray pulsars.

Depending on the dipole magnetic moment of a NS and the mass accretion rate, the inner parts of the accretion disc may appear in different regimes. For classical X-ray pulsars (μ ∼ 10³⁰ G cm³, ṁ ≲ 10), the accretion disc remains thin and gas-pressure-dominated, which implies a nearly classical scaling R_in ∝ ṁ − 0.22. As the mass accretion rate increases, a large portion of the disc can exist in a radiation-pressure-dominated regime. Unlike in the radiation-pressure-dominated disc without advection where the inner radius is independent of accretion rate (CAP), the present model with advection yields R_in ∝ ṁ^{−(0.05−0.1)}. For ULXPs, where the mass accretion rates reach ṁ₀ṁ₀∼ṁ₀10³−10⁴, the effects of advection and mass loss by the wind make the scalings similar to a spherical accretion case, resulting in a trend of about R_in ∝ ṁ^{−(0.2−0.3)}. However, if the strength of the magnetic field is one to two orders of magnitude higher than the usual pulsar values, the magnetospheric radius is larger, making the inner disc again geometrically thin and leading to a nearly flat dependence of the magnetospheric radius on ṁ.

Predictions for the magnetospheric radius can be tested through timing analysis of the stochastic component of the variability of X-ray sources, where quasi-periodic oscillations and breaks in the power-density spectra likely trace the dynamical timescales at the inner rim of the disc. Another test is spin period dynamics. We can compare the observed ṗ with theoretical predictions, as was done here for NGC 5907 X-1 and NGC 300 ULX1. The constraints we obtain from the observations of these two objects allow for a rather broad range of magnetic moments. However, the observational data set upper limits for beaming (not more than a factor of 5–6), confirming that ULXPs are intrinsically very luminous objects rather than X-ray sources whose luminosity is amplified by an order of magnitude or more as a result of anisotropy of their emission.

Acknowledgments

This research was supported by the grant 14.W03.31.0021 of the Ministry of Science and Higher Education of the Russian Federation. We also acknowledge the support from the Program of development of M.V. Lomonosov Moscow State University (Leading Scientific School “Physics of stars, relativistic objects and galaxies”) (GL, PA) and from the Väisälä foundation (AC). The authors thank Konstantin Malanchev for the fruitful discussion. AC and PA thank the lecturers and participants of the Astrosoma summer school (http://astrosoma.org) who provided an excellent environment for the development of this paper.

References

Appendix A: Disc vertical structure

For the vertical structure of the disc, we assume the form ρ = ρ_c(1 − x²)ⁿ, where x = z/H and n can be viewed as an effective vertical polytropic index. The vertical hydrostatic relation

(A.1)

implies a similar law for the vertical pressure profile P = P_c(1−x²)^n + 1, and a scaling relation for the disc thickness

(A.2)

where the surface density and vertically integrated pressure are related to the midplane quantities as

(A.3)

(A.4)

Here

(A.5)

A similar model for the vertical structure was considered by Kato et al. (2008). In the radiation-pressure-dominated regime T ∝ P^1/4, and therefore T ∝ (1 − x²)^{(n + 1)/4}. This approximation is not valid near the disc surface because T_eff ≠ 0. The vertical radiative energy flux is determined by the vertical radiation diffusion equation,

(A.6)

where ϵ = aT⁴ is the radiation energy density and D = c/(3κρ) is the diffusion coefficient. We denote the total energy release per unit surface area by Q⁺, and the total energy leaving the two sides of the disc by

(A.7)

The diffusion approximation allows us to connect the effective temperature to the temperature gradient inside the disc as

(A.8)

Taking into account Eqs. (A.3) and (A.4), this expression implies

(A.9)

and hence

(A.10)

The midplane gas pressure is P_g = k T_cρ_c/m̃. The gas-to-total pressure ratio then equals

(A.11)

where m̃ is the mean particle mass (for completely ionized gas of solar metallicity ṁ ≈ 0.6m_p).

Appendix B: Derivation of advection equations

The advective flux is given by expression (11), which contains the specific (per particle) dimensionless entropy

(B.1)

The radial derivative of the entropy, under our assumptions about the vertical structure, can be written as

(B.2)

where

(B.3)

(B.4)

(B.5)

and

(B.6)

with β being the gas-to-total pressure ratio introduced in Eq. (A.11). Vertical integration of Eq. (11) allows us to write the advective flux as

(B.7)

where the dimensionless constants are

(B.8)

(B.9)

(B.10)

(B.11)

To shorten the notations in Eq. (13), we also use the following combinations:

(B.12)

(B.13)

(B.14)

(B.15)

(B.16)

Appendix C: Dimensionless notations and equations

Here we list the dimensionless parameters and combinations used throughout the paper. Our notations here are identical to those in CAP. We normalize the NS mass as

(C.1)

The radius and the disc thickness r = R/R_g and h = H/R_g are measured in units of the gravitational radius R_g. The angular frequency is normalized by the local Keplerian frequency as

(C.2)

The characteristic magnetic moments of NSs lie in the range 10²⁸−10³² G cm³, hence we normalize μ as

(C.3)

The mass accretion rate is normalized by the Eddington value as

(C.4)

where

(C.5)

It is convenient to express the surface density in the units of the inverse opacity κ⁻¹. This quantity also has the physical meaning of the disc vertical optical depth

(C.6)

The dimensionless version of the vertically integrated tangential stress can be constructed as

(C.7)

For temperatures we use the normalization

(C.8)

where

(C.9)

The inner radius of the disc can be normalized either by the gravitational or by the Alfvén radius

(C.10)

where the dimensionless Alfvén radius is

(C.11)

and

(C.12)

We also introduce the natural time unit

(C.13)

that can be viewed as a Keplerian rotation period at R_g, and the dimensionless factor

(C.14)

The physical meaning of χ is the square of the dimensionless speed of sound (c_s/c)² corresponding to the characteristic temperature T_*.

Here we give all the equations in dimensionless form, as they were used to calculate the disc structure. The angular velocity at the inner boundary can be found from Eq. (1):

(C.15)

where p is the neutron star period in seconds. From Eq. (2) we can find the stress tensor at the boundary of the disc

(C.16)

The ratio of the gas pressure to the total pressure is, according to Eq. (A.11),

(C.17)

The differential Eqs. (15)–(20) in the dimensionless form are

(C.18)

(C.19)

(C.20)

(C.21)

Taking into account the sign of dΩ/dR <  0, we re-write Eq. (20) as

(C.22)

where the dimensionless versions of coefficients (B.12)–(B.15) are

(C.23)

(C.24)

(C.25)

(C.26)

All Tables

All Figures

Fig. 1. Structure of an accretion disc around a ULXP. For very high mass accretion rates, the inner parts of the disc (inside the spherization radius Rsph shown here) enter the super-Eddington accretion regime. Inside Rsph the thin disc model is not applicable. Mass loss in a wind is shown by blue arrows. The accretion column, where most of the energy is released, is shown by yellow cones, and the red wavy lines refer to the radiation of the column that may affect the inner disc pressure balance. The red vertical line marks the effective boundary between the disc and the magnetosphere at Rin.

In the text

	Fig. 2. Relative magnetospheric radius ξeff given by Eq. (23) is colour-coded and shown by contours on the ṁ0− μ plane. The spin period is p = 10 peq. The luminosity scale is calculated assuming efficiency η = 0.1, but irradiation effects on the disc structure were ignored.
In the text

	Fig. 3. Relative magnetospheric radii ξ (dashed blue curve) and ξeff (green solid) as functions of the disc thickness at the magnetospheric boundary. The dotted red line corresponds to the long-period asymptotic given by Eq. (57) in CAP. Magnetic moment of the NS was set to μ30 = 100.
In the text

	Fig. 4. Fraction of the angular momentum flux retained in the super-Eddington disc lin/l0 as a function of the fraction of mass reaching the magnetosphere. The dotted black line corresponds to equal ratios (lin/l0 = Ṁin/Ṁ0, as expected if the net angular momentum is constant); the solid green line corresponds to the expected scaling of a growing spherization radius (Eq. (26)).
In the text

	Fig. 5. Parameter ξeff as a function of fastness parameter for magnetic dipole moments μ30 = 1 and μ30 = 10 and mass accretion rates ṁ0 = 104 and ṁ0 = 5 × 103.
In the text

	Fig. 6. Ratio of the spherization radius rsph (in units of gravitational radius Rg) to the dimensionless mass accretion rate ṁ0 = 104 as a function of the fastness parameter (green dots). The red horizontal line gives the result from Eq. (21) in Poutanen et al. (2007).
In the text

	Fig. 7. Colour-coded contours of relative correction to the magnetospheric radius caused by irradiation ξ(η = 0.1)/ξ(η = 0)−1 shown in the ṁ0 − μ plane.
In the text

	Fig. 8. Relative disc thickness as a function of the radius for a model with μ = 1030 G cm3, ṁ0 = 3 × 103, and p = 0.67 s. Our results are shown by the red dotted curve, whereas the solid green and blue curves correspond to the asymptotics for the zones B and A of the standard disc, respectively.
In the text

	Fig. 9. Advection effects for the model shown in Fig. 8. The advection flux fraction Qadv/Q+ and the radiated fraction Qrad/Q+ are shown by the red solid and green dashed curves, respectively.
In the text

	Fig. 10. Cumulative luminosities as functions of radius. The total power generated in the disc is shown by the dashed blue curve. The solid red and dotted black curves show the radiative and advected luminosities, respectively. The parameters of the model are μ30 = 1 and ṁ0 = 3000.
In the text

	Fig. 11. Fraction of the mass accretion rate reaching radius R for a NS with magnetic moment μ30 = 1. The lines from top to bottom correspond to different ṁ0 in the interval 2000–3500 with steps of 500.
In the text

	Fig. 12. Magnetospheric radius in units of Rg as a function of the accretion rate for a NS with magnetic moment μ30 = 1. Parts of the black solid curve with different slopes correspond to the different regimes of accretion near the magnetospheric boundary. Two standard solutions () are plotted with the grey dotted lines: ξ = 0.5 and ξ = 1 (spherically symmetric case).
In the text

	Fig. 13. Same as Fig. 12, but for different magnetic moments μ = 1030−1032 G cm3.
In the text

	Fig. 14. Slope of the dependence of the magnetospheric radius on the accretion rate for the models shown in Fig. 13.
In the text

Fig. 15. Restrictions for the position of NGC 5907 X-1 in the ṁ0 − μ plane. The red horizontal line corresponds to the mass accretion rate of ṁ0 = 6500 (accretion efficiency η ≃ 0.15). The black solid line corresponds to peq = p = 1.137 s. The dotted spin-up line shows the lower limit for ṁin set by inequality (32) using period derivation |ṗ| = 5 × 10−9 s s−1. We use here M = 1.4 M⊙, R = 10 km, I = 1.5 × 1045 g cm2, and p = 1.137 s.

In the text

Fig. 16. Restrictions for the position of NGC 300 ULX1 in the ṁ0 − μ plane. The red horizontal line corresponds to the luminosity of 4.7 × 1039 erg s−1 assuming accretion efficiency η ≃ 0.15. The black solid line corresponds to peq = p = 20 s. The dotted spin-up line shows the lower limit for ṁ0 set by inequality (32) using period derivation |ṗ| = 5.5 × 10−7 s s−1 (Carpano et al. 2018). We use here M = 1.4 M⊙, R = 10 km, I = 1.5 × 1045 g cm2, and p = 20 s.

In the text