Super-Eddington accretion discs with advection and outflows around magnetized neutron stars

TL;DR

This paper introduces a semi-analytical model for super-Eddington accretion discs around magnetized neutron stars, accounting for heat advection and wind-driven mass loss, predicting physical profiles and magnetospheric radius behavior.

Contribution

It provides a novel semi-analytical framework for super-Eddington accretion discs with advection and outflows, including the magnetospheric radius as an eigenvalue, applicable to various neutron star systems.

Findings

Magnetospheric radius depends weakly on accretion rate when advection is mild.

Approaching the Eddington limit, the disc becomes advection-dominated with classical Alfven scaling.

Mass loss via wind increases the magnetospheric radius.

Abstract

We present a model for a super-Eddington accretion disc around a magnetized neutron star taking into account advection of heat and the mass loss by the wind. The model is semi-analytical and predicts radial profiles of all basic physical characteristics of the accretion disc. The magnetospheric radius is found as an eigenvalue of the problem. When the inner disc is in radiation-pressure-dominated regime but does not reach its local Eddington limit, advection is mild, and the radius of the magnetosphere depends weakly on the accretion rate. Once approaching the local Eddington limit, the disc becomes advection-dominated, and the scaling for the magnetospheric radius with the mass accretion rate is similar to the classical Alfven relation. Allowing for the mass loss in a wind leads to an increase of the magnetospheric radius. Our model may be applied to a large variety of magnetized neutron stars accreting close to or above their Eddington limits: ultra-luminous X-ray pulsars, Be/X-ray binaries in outbursts, and other systems. In the context of our model we discuss the observational properties of NGC 5907~X-1, the brightest ultra-luminous pulsar known so far, and NGC 300~ULX-1 which is apparently a Be/X-ray binary experiencing a very bright super-Eddington outburst.

Figures

a in Basko & Sunyaev [1976](#bib.bib7)). Large ratio of the radiating surface area to the volume as well as 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](#bib.bib39); Kawashima et al., [2016](#bib.bib27)).

Tables

Table 1. Table 1: Disc properties.

parameters	$\xi$	${\left(H/R\right)}_{\mathrm{in}}$	${\left(H/R\right)}_{\max }$	${\dot{M}}_{\mathrm{in}}/{\dot{M}}_0$	$l_{\mathrm{in}}/l_0$
fiducial	$1.16$	$0.97$	$1.04$	$0.88$	$0.82$
$\eta =0.1$	$1.31$	$1.08$	$1.09$	$0.84$	$0.73$
${ϵ}_w=1$	$1.18$	$0.93$	$1.02$	$0.91$	$0.88$
${ϵ}_w=0.1$	$1.13$	$0.99$	$1.06$	$0.97$	$0.95$
$n=1.5$	$1.10$	$0.92$	$1.00$	$0.99$	$0.99$
$n=3$	$1.05$	$0.81$	$0.91$	$1$	$1$
$\psi =1.5$	$1.15$	$0.96$	$1.02$	$0.91$	$0.81$
$\alpha =0.5$	$0.68$	$0.61$	$0.64$	$1$	$1$

Equations

$${\dot{M}}_{\mathrm{in}}({\Omega }_{\mathrm{in}}-{\Omega }_{\mathrm{NS}})R_{\mathrm{in}}^2=k_{\mathrm{t}}\frac{{\mu }^2{H}_{\mathrm{in}}}{R_{\mathrm{in}}^4}+L\frac{{\Omega }_{\mathrm{in}}}{c^2}{H}_{\mathrm{in}}{R}_{\mathrm{in}},$$

$$W_{r\varphi }^{\mathrm{in}}=2\alpha H_{\mathrm{in}}\left(\frac{{\mu }^2}{8\text{uppi}{R}_{\mathrm{in}}^6}+\frac{L}{4\text{uppi}{R}_{\mathrm{in}}^{2}c}\right).$$

$${\Omega }^{2}R=\frac{1}{\Sigma }\frac{\text{rotatebox[origin=t]20.0scalebox0.95[1.0]∂}\Pi }{\text{rotatebox[origin=t]20.0scalebox0.95[1.0]∂}R}+\frac{GM}{R^2},$$

$$\frac{\mathrm{d}\left(\dot{M}(R)\Omega R^2\right)}{\mathrm{d}R}=\frac{\mathrm{d}}{\mathrm{d}R}\left(2\text{uppi}{R}^2{W}_{r\varphi }\right)+\frac{\mathrm{d}\dot{M}(R)}{\mathrm{d}R}\Omega R^2\psi ,$$

$${ϵ}_{\mathrm{w}}{Q}_{\mathrm{rad}}={ϵ}_{\mathrm{w}}2{\sigma }_{\mathrm{SB}}{T}_{\mathrm{eff}}^4=\frac{{\Omega }_{\mathrm{K}}^{2}R}{4\pi }\frac{\mathrm{d}\dot{M}(R)}{\mathrm{d}R},$$

$$W_{r\varphi }=\alpha \Pi .$$

$$\frac{1}{G_{n+1}}\frac{W_{r\varphi }}{\alpha H}=\frac{a{T}_{\mathrm{c}}^4}{3}+\frac{1}{G_n}\frac{\Sigma k{T}_{\mathrm{c}}}{H\tilde{m}},$$

$$Q^{+}=Q_{\mathrm{rad}}+Q_{\mathrm{adv}},$$

$$Q^{+}=W_{r\varphi }R\left|\frac{\mathrm{d}\Omega }{\mathrm{d}R}\right|,$$

$$Q_{\mathrm{rad}}=2{\sigma }_{\mathrm{SB}}{T}_{\mathrm{eff}}^4=\frac{16}{3\kappa \Sigma }(n+1)G_n{\sigma }_{\mathrm{SB}}{T}_{\mathrm{c}}^4.$$

$$Q_{\mathrm{adv}}={\int }_{-H}^H\rho v_R\frac{kT}{\tilde{m}}\frac{\mathrm{d}s}{\mathrm{d}R}\mathrm{d}z,$$

$$Q_{\mathrm{adv}}$$

$$W_{r\varphi }R\left|\frac{\mathrm{d}\Omega }{\mathrm{d}R}\right|=\frac{16}{3\kappa \Sigma }(n+1)G_n{\sigma }_{\mathrm{SB}}{T}_{\mathrm{c}}^4-\frac{1}{2\pi (n+1)}\frac{1}{G_n}\frac{\dot{M}{W}_{r\varphi }}{\alpha R\Sigma }$$

$$\left|\frac{\mathrm{d}\Omega }{\mathrm{d}R}\right|=C_{\Omega }-\frac{\mathrm{d}\ln \Sigma }{\mathrm{d}R}{C}_{\Sigma }-\frac{\mathrm{d}\ln W_{\mathrm{r}\varphi }}{\mathrm{d}R}{C}_{\mathrm{wrf}}-\frac{\mathrm{d}\ln T_{\mathrm{c}}}{\mathrm{d}R}{C}_{\mathrm{T}}-C_{\mathrm{free}},$$

$$\frac{\mathrm{d}{W}_{r\varphi }}{\mathrm{d}R}=\alpha \Sigma \left({\Omega }^{2}R-\frac{GM}{R^2}\right).$$

$$\frac{\mathrm{d}\dot{M}}{\mathrm{d}R}={ϵ}_{\mathrm{w}}\frac{8\pi {\sigma }_{\mathrm{SB}}{T}_{\mathrm{eff}}^4}{{\Omega }_{\mathrm{K}}^{2}R}={ϵ}_{\mathrm{w}}\frac{64\pi {\sigma }_{\mathrm{SB}}}{3\kappa \Sigma {\Omega }_{\mathrm{K}}^{2}R}(n+1)G_n{T}_{\mathrm{c}}^4.$$

$$\frac{\mathrm{d}\ln T_{\mathrm{c}}}{\mathrm{d}R}$$

$$\beta =\frac{P_{\mathrm{g}}}{P_{\mathrm{tot}}}=\frac{2(n+1)}{2n+3}\frac{\alpha \Sigma k{T}_{\mathrm{c}}}{\tilde{m}\Pi }.$$

$$\frac{\mathrm{d}\Omega }{\mathrm{d}R}$$

$$\frac{\mathrm{d}\ln \Sigma }{\mathrm{d}R}$$

$$\xi =\frac{R_{\mathrm{in}}}{R_{\mathrm{A}}},$$

$$R_{\mathrm{A}}={\left(\frac{{\mu }^2}{2{\dot{M}}_0\sqrt{2GM}}\right)}^{2/7}$$

$${\xi }_{\mathrm{eff}}=R_{\mathrm{in}}{\left(\frac{{\mu }^2}{2{\dot{M}}_{\mathrm{in}}\sqrt{2GM}}\right)}^{-2/7}=\xi {\left(\frac{{\dot{M}}_{\mathrm{in}}}{{\dot{M}}_0}\right)}^{2/7}.$$

$$p_{\mathrm{eq}}=2\text{uppi}(GM{)}^{-5/7}{\left(\frac{{\mu }^2}{2\sqrt{2}\dot{M}}\right)}^{3/7}\simeq 1.3m^{-5/7}{\dot{m}}_0^{-3/7}{\mu }_{30}^{6/7}\mathrm{s},$$

$$l=\dot{M}\omega \sqrt{GMR}-2\pi R^2{W}_{r\varphi },$$

$$\frac{l_0-l_{\mathrm{in}}}{l_{\mathrm{in}}}\simeq \sqrt{\frac{R_{\mathrm{sph}}}{R_{\mathrm{in}}}}\frac{{\dot{M}}_0-{\dot{M}}_{\mathrm{in}}}{{\dot{M}}_{\mathrm{in}}}.$$

$${\omega }_{\mathrm{s}}=\frac{{\Omega }_{\mathrm{NS}}}{\Omega (R_{\mathrm{in}})}.$$

$${\dot{m}}_{\mathrm{cr}}\simeq 2000{\left(\frac{\alpha }{0.1}\right)}^{2/9}{\mu }_{30}^{4/9},$$

$$\frac{R_{\mathrm{in}}}{R_{\mathrm{g}}}\approx 170(\alpha /0.1{)}^{2/9}{\mu }_{30}^{4/9}{m}^{-10/9}.$$

$$\frac{\mathrm{d}(I{\Omega }_{\mathrm{NS}})}{\mathrm{d}t}=K_{\mathrm{su}}-K_{\mathrm{sd}},$$

Full text

1 Introduction

Mass transfer rates in binary systems may vary in very broad limits, from very low to the 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 time scale, being the thermal time scale of a tens-solar-mass star, which is millions of years. The mass transfer does not quench because the large 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). Large ratio of the radiating surface area to the volume as well as 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).

Recently it was realised that some of ultraluminous X-ray sources (ULX) in the nearby galaxies are actually accreting magnetized NSs. Bachetti et al. (2014) using NuSTAR data discovered coherent pulsations in the ULX X-2 in the galaxy M82. Later, two more similar objects, NGC 7793 P13 (Israel et al., 2017b; Fürst et al., 2016) and NGC 5907 X-1 (Israel et al., 2017a) were found. NGC 5907 X-1 is exceptional in its luminosity, exceeding $10^{41}\rm erg\,s^{-1}$ 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 20–30 s period were discovered in NGC 300 ULX1, the supernova impostor SN2010da (Villar et al., 2016), with a very strong spin up of $\dot{p}=-1.75\times 10^{-7}$ s s*-1* (Bachetti et al., 2018). Also, M51 ULX-7 was identified with a ULXP (G. Israel, priv. communication), 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 few persistent super-Eddington objects are robustly identified so far as NSs, 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 $\dot{M}_{0}$ exceeding a certain critical value $\dot{M}_{\rm cr}$, a wind emanates inside the spherization radius $R_{\mathrm{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_{\rm Edd}$ by a logarithmic factor $1+\ln(\dot{M}_{0}/\dot{M}_{\rm cr})$, and the accretion rate decreases with radius as $\dot{M}(R)=\dot{M}_{0}\,R/R_{\rm sph}$. We will 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). 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 $\dot{M}_{0}$ converge to universal values dependent on the NS magnetic moment $\mu$ only.

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, that violates the locality assumption in the energy balance of the disc. The role of advection for the case of optically thin discs was considered, for example, by Ichimaru (1977), Narayan & Yi (1994, 1995), Abramowicz et al. (1995, 1996), and for the optically thick case by Begelman & Meier (1982) and Abramowicz et al. (1988). Numerical simulations of supercritical accretion onto a black hole were performed by Eggum et al. (1988), Beloborodov (1998), 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. Recently, simulations of supercritical accretion onto a non-magnetized NS were performed by Abarca et al. (2018). 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 sub-critical regime. Chashkina et al. (2017) (hereafter CAP) considered accretion onto a magnetized NS in a regime that may be characterized as intermediate, with the accretion rate in the range $(10^{-9}-10^{-6}){\rm\,M_{\odot}\,yr^{-1}}$ for $\mu\sim 10^{31}-10^{32}{\rm\,G\,cm^{3}}$ 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 accreting at larger rates should possess a supercritical advective disc, whose radiation pressure becomes sufficient to unbind part of the accreted material.

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 and luminosities of hundreds of 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. Sect. 3 is devoted to the results, where we present the dependencies 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 will consider a NS of mass $M$, radius $R_{\rm NS}$ and magnetic moment $\mu$. We assume it is an aligned rotator with the unperturbed dipolar magnetic field inside the magnetospheric radius, which is equal to the inner radius of the disc $R_{\rm in}$. Outside this radius, the magnetic field lines are opened by the ideally conducting accretion disc in the way it was proposed by Parfrey et al. (2016). 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 described by only two boundary conditions introduced in CAP and later in this section. A simple sketch of the adopted structure of the disc is shown in Fig. 1.

In our model, we include three effects important for an accretion disc at near- or super-critical accretion rate: (i) advection, (ii) mass loss in the wind, and (iii) angular momentum losses in the wind. Advection is related to the increasing photon diffusion time scales at large 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 remove efficiently the angular momentum from the disc, as in a 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

[TABLE]

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

[TABLE]

Here, $\dot{M}_{\rm in}$ is the mass accretion rate at the inner edge of the disc, $\Omega_{\rm NS}$ is the angular velocity of the NS, $\Omega_{\rm in}$, $H_{\rm in}$, and $W_{r\phi}^{\rm in}$ are the angular velocity, the half-thickness and the vertically integrated $r\phi$- component of the viscous stress tensor at the inner boundary of the disc $R_{\rm in}$, the dimensionless constant $k_{\rm t}$ parametrizes the efficiency of the angular momentum removal by the magnetic and viscous torques, and $\alpha$ is dimensionless viscosity parameter as defined by Shakura & Sunyaev (1973). The luminosity $L=\eta\,\dot{M}_{\rm in}\,c^{2}$ is released close to the NS in the accretion column with efficiency $\eta$ 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 $\eta$ may differ from the angle-averaged efficiency of an isotropic source as well as from the efficiency measured by an observer at infinity. The radial structure of the disc is described by the angular frequency $\Omega$, thickness $H$, and the viscous stress $W_{r\phi}$. 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 large luminosities $L>L_{\rm Edd}$. Note that equation (1) is valid only for the case of the disc rotating faster than the magnetosphere. For the propeller case ($\Omega_{\rm in}<\Omega_{\rm NS}$), the boundary condition would be different. Moreover, one could consider a radial distribution of the magnetic torque applied to the disc, as done, for example, by Kluzniak & Rappaport (2007). On the other hand, Matt & Pudritz (2005b) argue that if $\Omega_{\rm in}>\Omega_{\rm NS}$, the magnetic interaction can be limited to a narrow ring.

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

[TABLE]

where $\Sigma=\displaystyle\int^{H}_{-H}\rho\,{\rm d}z$ is the surface density and $\Pi=\displaystyle\int^{H}_{-H}P\,{\rm d}z$ is the vertically-integrated pressure. As in CAP, we ignore the dynamic term $v_{R}\,{\rm d}v_{R}/{\rm d}R$ which is suppressed by a factor of $\alpha^{2}(H/R)^{2}\ll 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 $\rho=\rho_{\rm c}[1-(z/H)^{2}]^{n}$ and for the pressure $P=P_{\rm c}[1-(z/H)^{2}]^{n+1}$ (Paczynski & Jaroszynski, 1978). This allows us to compute analytically the connection factors between the midplane quantities such as $\rho_{\rm c}$, $P_{\rm c}$, and the corresponding vertically-integrated values $\Sigma$, $\Pi$ as described in Appendix A. For the vertical temperature profile used in the advection term calculations, we use the expression $T=T_{\rm c}[1-(z/H)^{2}]^{(n+1)/4}$ valid for the radiation-pressure dominated case, as $T\propto p^{1/4}$, here $T_{\rm 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:

[TABLE]

where $\psi\geq 1$ allows us to scale up the net angular momentum lost in the wind. If $\psi=1$, the net angular momentum in the wind is equal to that in the disc (Lipunova, 1999; Poutanen et al., 2007). Larger $\psi$ can appear in magneto-centrifugal winds, where the sub-Alfvenic part of the transonic flow rotates approximately rigidly, thus increasing the angular momentum loss by a factor of $\sim(R_{\rm wA}/R)^{2}$, where $R_{\rm 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_{\rm wA}/R$ may be as large as 5. Simulations of the centrifugally driven winds (e.g. Ustyugova et al., 1999), Lovelace1995 suggest even larger lever lengths at low accretion rate, 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 ${\rm d}\dot{M}/{\rm d}R$ needs to be calculated using some additional equations describing the physics of the wind launching. We assume that some fraction $\epsilon_{\rm w}\leq 1$ of the energy leaving the disc with radiation is spent to accelerate the outflow (Lipunova, 1999; Poutanen et al., 2007):

[TABLE]

where $\Omega_{\rm K}$ is Keplerian angular velocity and $T_{\rm 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_{\rm sph}$ defined by the condition on the minimal relative thickness $H/R>\left(H/R\right)_{\rm cr}$. Because the details of wind launching are uncertain, we assume $\left(H/R\right)_{\rm cr}=1$. 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 depends also on the magnetic moment of the star and its spin period.

We use the $\alpha$-viscosity prescription for the vertically integrated viscous stress

[TABLE]

The midplane pressure may be expressed as a sum of the radiation and gas pressures:

[TABLE]

where $\tilde{m}$ is the mean particle mass and $G_{n}=\displaystyle\int^{1}_{-1}(1-x^{2})^{n}{\rm d}x$, see equation (38). Here we have taken into account the $\alpha$-viscosity prescription and the relations between the central and vertically-integrated quantities implied by the adopted vertical structure (equations 36 and 37). We will not use equation (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 general case with advection, some fraction of energy is radially transported. The energy flux carried by radiation diffusion in the vertical direction is no more equal to the local energy release:

[TABLE]

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

[TABLE]

and the radiation flux from both sides of the disc is (see equation 43):

[TABLE]

Here $\kappa\simeq 0.34\,$cm2g*-1* is the Thomson scattering opacity. The advected flux $Q_{\rm adv}$ may be viewed (as it is done, for instance, in Lipunova 1999) as the flux of heat carried with the flow, and thus may be expressed through the specific entropy per particle $s$

[TABLE]

where the radial velocity (taken by its absolute value) is $v_{R}=\dot{M}/(2\pi R\Sigma)$. Taking into account equation (7) we can expand the expression above as (see Appendix B)

[TABLE]

where the coefficients ${\cal S}$, ${\cal P}$, ${\cal Q}$, and ${\cal R}$ are given by equations (B)–(55). The energy balance can be written using equations (8), (9) and (12) as:

[TABLE]

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

[TABLE]

where the coefficients $C_{\Omega,\Sigma,w,{\rm T,free}}$ are given by expressions (56)–(60).

2.3 Solving the disc equations

Unlike the previous work (CAP), the number of variables we need to solve differential equations for is five: the angular frequency $\Omega$, the tangential stress $W_{r\phi}$, the surface density $\Sigma$, the mass accretion rate $\dot{M}$, and the midplane temperature $T_{\rm c}$. There are also several adjustable free parameters: viscosity parameter $\alpha$, wind efficiency parameter $\epsilon_{\rm w}$, wind magnetization parameter $\psi$, polytropic index $n$, spin period $p$ and the magnetic moment of a NS $\mu$, and the accretion efficiency $\eta$. We calculate the disc structure from its outer edge to the inner boundary. We solve the five equations listed below:

The derivative of vertically-integrated viscous stress that can be obtained from the radial Euler equation (3) and $\alpha$-viscosity prescription (6):

[TABLE] 2. 2.

Whenever the outflow condition ($H/R>(H/R)_{\rm cr}$) is satisfied, the mass accretion rate changes with radius according to equation (5)

[TABLE] 3. 3.

The differential equation for the central temperature may be obtained by taking derivative of the pressure equation (7)

[TABLE]

Here, we used the thickness of the disc $H$ following from the hydrostatic equilibrium (see equation 35) together with the gradient of $W_{r\phi}$ substituted from equation (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 equation 44)

[TABLE] 4. 4.

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

[TABLE]

The second term on the right-hand side is switched off for a thin disc ($H/R<(H/R)_{\rm cr}$) when there no outflows. 5. 5.

The differential equation for the surface density was obtained from the advection equation (14) substituting all other derivatives from equations (15)–(4):

[TABLE]

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

For two of the variables, $\Omega$ and $W_{r\phi}$, the boundary conditions exist that need to be satisfied. All the quantities at the outer boundary of the integration region should conform to the thin disc model with $\dot{M}=\dot{M}_{0}$ being one of the global parameters of the simulation and $W_{r\phi}=W_{r\phi}^{\rm out}$, an adjustable parameter varied independently of the relative magnetospheric radius $\xi$ to satisfy the inner boundary conditions:

[TABLE]

where

[TABLE]

is the Alfvén radius (see, for instance, Elsner & Lamb 1977 and section 6.3 in Frank et al. 2002).

Because the accretion rate at the inner boundary of the disc is physically more relevant, we will also use a differently normalized version of the relative magnetospheric radius:

[TABLE]

3 Results

3.1 Global parameters

The relative magnetospheric size $\xi$ is found as one of the two eigenvalues of the problem, the other being $W_{r\phi}^{\rm out}$. Other global parameters obtained in the model include several quantities at the magnetospheric boundary: the relative thickness of the disc $\left(H/R\right)_{\rm in}$, the fraction of mass reaching the magnetosphere $\dot{M}_{\rm in}/\dot{M}_{0}$, and the advection fraction $\left(Q_{\rm adv}/Q^{+}\right)_{\rm in}$. We also track the maximal thickness of the disc $\left(H/R\right)_{\rm max}$.

All the simulations were made for the parameters: $\alpha=0.1$, $k_{\rm t}=0.5$, $\epsilon_{\rm w}=0.5$, $\psi=1$, $n=1$ and $\eta=0$ if not stated otherwise. We first consider a slowly rotating accretor with the spin period equal to ten equilibrium periods $p=10p_{\rm eq}$,111When $p\gg p_{\rm eq}$, the exact value of the spin period barely affects $\xi_{\rm eff}$, see Sect. 3.2. where the equilibrium period ( Lipunov 1992, Illarionov & Sunyaev 1975) is defined as the spin period for which the Alfvén radius $R_{\rm A}$ equals to the corotation radius $R_{\rm co}=\left(\Omega_{\rm NS}^{2}/GM\right)^{1/3}$ :

[TABLE]

where $\mu_{30}=\mu/10^{30}$ G cm3 is dimensionless NS magnetic moment, $m=M/1.4M_{\odot}$ is normalized NS mass. We also use dimensionless mass accretion rate $\dot{m}_{0}=\dot{M}_{0}/\dot{M}_{\rm Edd}$ normalized by the Eddington value $\dot{M}_{\rm Edd}=4\uppi GM/c\kappa$. In Fig. 2, we show the contours of $\xi_{\rm eff}$ in the $\dot{m}_{0}-\mu$ plane, covering two orders of magnitude in magnetic moment and five orders of magnitude in the mass accretion rate.

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

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 ($\psi=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,

[TABLE]

where $\omega=\Omega/\Omega_{\rm K}$. Unless some of this angular momentum is removed in the wind, $l_{\rm in}=l(R_{\rm in})$ should be equal to $l_{0}=l(R_{\rm out})$. In Fig. 4, we show the ratio of the angular momentum fluxes $l_{\rm in}/l_{0}$ as a function of the ratio of mass accretion rates $\dot{M}_{\rm in}/\dot{M}_{0}$. Both ratios start at unity for conservative thin disc accretion and then decrease as $\dot{M}_{0}$ 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 $R_{\rm sph}=\displaystyle\frac{3}{2}\dot{m}_{0}\frac{GM}{c^{2}}$ that leads to the scaling (see Fig. 4):

[TABLE]

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

3.2 Dependence on the spin period

The transition to the propeller regime, when the accretion flow cannot spin up the NS anymore, may be traced using the fastness parameter $\omega_{\rm s}$ defined as

[TABLE]

When $\omega_{\rm 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 fastness parameter (or decreasing spin period of the accretor) on the properties of our solution is in the increasing ratio $\xi$ of the size of the magnetosphere to the Alfven radius. Factor $\xi_{\rm eff}$ increases by about 40% between the slowly rotating NS case and the propeller limit (see Fig. 5).

Near the corotation, when $R_{\rm co}=R_{\rm in}$, 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, e.g. by Poutanen et al. (2007). Their model takes into account the advection effects and the outflows from the disc. We compare spherization radius defined as the maximal distance from the NS where the condition for outflows, $H>R$, is fulfilled.

In Fig. 6, we show the spherization radius normalized by dimensionless mass accretion rate (green dots) and compare our results to the zero-torque case. For $\epsilon_{\rm w}=0.5$ and $\dot{m}_{0}=10^{4}$, equation (21) by Poutanen et al. (2007) predicts $R_{\rm sph}\simeq 0.575\dot{m}_{0}({GM}/{c^{2}})$, in reasonable agreement with our results close to equilibrium.

3.3 Effects of irradiation

All the previous results were calculated without irradiation effects by setting $\eta=0$ (or $L=0$ in equations 1 and 2). The real efficiency affecting the pressure balance is probably of the same order with 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 and the radius of the magnetosphere, in particular. The magnetospheric radius increases by up to 30% for $\dot{m}_{0}\sim 10^{4}$, as it is shown in Fig. 7. The effect grows rapidly with the mass accretion rate, as well as with the magnetic moment. The latter is a consequence of a rapidly growing ratio $P_{\mathrm{rad}}/P_{\mathrm{mag}}$ with the disc inner radius in the dipole approximation (see for instance CAP, equation 63).

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 $\alpha=0.1$, $\epsilon_{\rm w}=0.5$, $n=1$, $\eta=0$, and $\psi=1$. All the models are calculated for a NS with magnetic field $\mu=10^{30}$ G cm3 and accretion rate $\dot{m}_{0}=3000$. Other models differ from the fiducial one in a single parameter.

As in CAP, we find viscosity parameter to be an important factor altering the structure of the disc. For the parameters under consideration, changing $\alpha$ from 0.1 to 0.5 results in a twofold decrease in disc thickness and $\xi$ and quenches wind formation. Changing 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 minimal disc thickness for wind launching, advection starts to play a role rather early, when all the disc is still sub-critical. As a consequence, the disc thickness stabilizes at $H\simeq R$ (see Fig. 8). Slimming effect of advection was noted earlier by, for example, Abramowicz et al. (1988), Beloborodov (1998), Lipunova (1999), and Lasota et al. (2016). The local Eddington limit ($H=R$) is reached at the critical mass accretion rate of

[TABLE]

that is 5–6 times higher than in CAP (equation 66).

Advection starts to dominate in the energy balance already below this limit, at $\dot{m}_{0}\sim 10^{3}$ (see Fig. 9), making the inner disc a huge reservoir of heat. For $\dot{m}_{0}$ above the limit given by equation (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 $L_{\rm tot}=\int_{R}^{R_{\rm out}}2\pi Q^{+}RdR$. Similarly, by integrating $Q_{\rm rad}$ and $Q_{\rm 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_{\rm in}$ it is in full agreement with theoretical prediction, $L_{\rm theor}=GM\dot{M}/2R_{\rm in}$. We note that the ratio of the disc luminosity to the luminosity of the accretion column can be as low as $\sim R_{\rm NS}/R_{\rm in}$.

In Fig. 11, we plot the fraction of the initial mass accretion rate remaining in the disc for different $\dot{m}_{0}$ as a function of radius. At very high accretion rates, supercritical wind blows away a considerable part of the accretion material and operates at all radii within $R_{\rm sph}$. At some intermediate rates, $2000\lesssim\dot{m}_{0}\lesssim 3000$, there is a prominent sub-critical region near $R_{\rm in}$, where there is no wind. This is caused by the non-monotonic dependence of the disc height on radius. We note that amount of the blown-away material depends on the condition which is used to switch on the wind. The wind launching relies on complex physical process and can be described only very approximately in 1-D models.

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), magnetospheric radius is supposed to scale with the Alfvén radius, corresponding to constant $\xi$ in our notation. One of the important results of CAP was understanding that the physical radius of the magnetosphere interacting with a thin radiation-pressure-dominated disc barely changes with mass accretion rate. Accounting for the effects of advection and wind losses makes the picture more complicated.

Magnetospheric radius (in units of gravitational radius $R_{\rm g}=GM/c^{2}$) dependence 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, $\mu\sim 10^{30}$ G cm3, this happens at luminosities of a few times $L_{\mathrm{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). As we have 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 sub-critical disc:

[TABLE]

Thus, provided with 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 $\alpha$.

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 $\dot{m}_{0}$ is crucial for the behaviour of $r_{\mathrm{in}}(\dot{m}_{0})$. As the accretion rate increases, advection starts to play important role. The relative thickness of the disc is no more proportional to $\dot{m}_{0}$, and the magnetospheric radius again depends on $\dot{m}_{0}$. The interplay between wind losses and advection makes the radius dependence on mass accretion rate shallower than the $\xi=const$ approximation historically proposed for spherical accretion but stronger than $r_{\rm 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 $\mu\sim 10^{32}$G cm3. 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_{\rm ab}\propto\dot{m}^{16/21}$ (Shakura & Sunyaev, 1973). Because $R_{\rm in}\propto\mu^{4/9}$ in the radiation-pressure-dominated regime (see Eq. 29), the left boundary of the plateau depends on the magnetic moment as $\dot{m}_{\rm in,left}\propto\mu^{7/12}$. The right boundary of the plateau is determined by advection effects. Advection becomes important when $Q_{\rm adv}\sim Q_{\rm rad}$, that implies $H_{\rm in}\sim R_{\rm 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 $\dot{m}_{\rm in,right}\propto\mu^{4/9}$. Thus the length of the plateau slowly decreases with magnetic moment as $\dot{m}_{\rm in,right}/\dot{m}_{\rm in,left}\propto\mu^{-5/36}$.

It is interesting to compare our results with the classical prescriptions (Ghosh et al., 1977; Wang, 1996). Fig. 12 demonstrates that the classical dependencies are much steeper, having the slope of $\delta\equiv{\rm d}\log r_{\rm in}/{\rm d}\log{\dot{m}}={-2/7}$. In the case of $\mu=10^{30}$ G cm3, the slope is $\delta\approx-0.21$ in the gas-pressure dominated case, $\delta\approx-0.07$ in the radiation-pressure-dominated regime, and $\delta\approx-0.16$ when advection dominates. Evolution of the local slope $\delta$ is traced in Fig. 14. Low accretion rate asymptotic corresponding to gas-pressure-dominated thin disc stably reproduces $\delta\sim-0.23$ in accordance with the results of CAP (see their Sect. 5.1). The maximal value of $\delta$ depends on how prominent is the thin radiation-pressure-dominated part of the disc for given magnetic field, changing from nearly zero for small $\mu$ to about $-0.1$ for magnetar-scale fields. Largest mass accretion rates tend to reproduce much steeper dependencies, approaching $\delta\simeq-2/7\simeq-0.29$.

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 during the 10 years of observations (Israel et al., 2017a). It is also remarkable that its luminosity exceeds $10^{41}\rm erg\,s^{-1}$ during some of the observations. The maximal detected period derivative (by absolute value) was $\dot{p}=-5\times 10^{-9}{\rm s\,s^{-1}}$, about an order of magnitude larger than the average value for this object. The large value of $|\dot{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

[TABLE]

where $I\simeq(1-2)\times 10^{45}{\rm\,g\,cm^{2}}$ is the moment of inertia of the NS, $K_{\rm su}=\dot{M}\sqrt{GMR_{\rm in}}$ and $K_{\rm sd}$ are spin-up and spin-down torques. Ignoring the unknown spin-down contribution allows us to get a constraint

[TABLE]

resulting in an inequality for accretion rate

[TABLE]

here $I_{45}=I/10^{45}{\rm\,g\,cm^{2}}$, and $\dot{p}_{-12}=\dot{p}/(10^{-12}$ s s*-1*). Alternatively, having an independent estimate for the accretion rate, we then can set a lower limit for the magnetic field of this object. For maximum accretion rate of $\sim 6000$, corresponding to isotropic X-ray emission, the normalised magnetic moment $\mu_{30}$ cannot be less than $\sim 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 magnetic field, suggesting the NS is still in the accretion regime. In terms of mass accretion rate, this condition may be written as

[TABLE]

Fig. 15 shows these two limits as applied to NGC 5907 X-1. The spin-up line shows the lower limit for $\dot{m}_{0}$ set by inequality (32) using period derivative $\dot{p}=-5\times 10^{-9}\,{\rm s\,s^{-1}}$. 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 corotational radius. The region below this line is prohibited because the accretion disc rotates slower than the magnetosphere of the NS, and no stable accretion is possible. Maximal bolometric luminosity of this object is about $2\times 10^{41}{\rm erg\,s^{-1}}$ (Israel et al., 2017a), that leads to a lower limit on the mass accretion rate $\dot{m}_{0}>\dot{m}_{\rm in}\sim 6500$ assuming efficiency $\eta\sim 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 an evidence for a truly high mass accretion rate in NGC 5907 X-1, as significantly low $\dot{m}_{0}\lesssim 10^{3}$ are forbidden for any magnetic moments. Therefore, any beaming exceeding a factor of several is unlikely because with increasing beaming, the red line in Fig. 15 will move down, so if beaming is too large there would be no allowed region. 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.

One can set an upper limit for the magnetic field of NGC 5907 X-1 as $\mu\leq 7.5\times 10^{31}$ G cm3 and $\mu\leq 5.45\times 10^{31}$ G cm3 if we take into account irradiation from the column (see Sect. 3.3). There is an evidence for bimodal distribution in luminosities, see Fig. S2 in Israel et al. (2017a) 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 $\dot{m}_{0}$ (as $\mu\propto\dot{m}_{0}^{1/2}$). The disc is expected to be supercritical (i.e. having winds) in its inner parts if $\mu_{30}\lesssim 14$ (see Eq. 28). At a pulsar magnetic field $\mu_{30}\sim 1$, the expected outflow rate from the disc is about 20% of the inflow rate, or $5\times 10^{-6}M_{\odot}\,{\rm yr}^{-1}$.

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\simeq 4.7\times 10^{39}\rm erg\,s^{-1}$, but the observed spin-up rate of this source is exceptionally high: its spin period has changed from 45 to 20 s in less than 2 years. Its spin period derivative was estimated as $\dot{p}\simeq-5.5\times 10^{-7}$s s*-1*, that is the highest $\dot{p}$ ever observed from an accreting NS. Large period makes NGC 300 ULX1 a promising candidate for a NS in a pure spin-up state. This may be seen in Fig. 16 where the observed properties of the object are in good agreement with a pulsar-scale $\mu\sim 10^{30}$G cm3 of a NS which 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 $\mu=(1.5-85)\times 10^{30}$ G cm3, 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_{\rm eq}\sim 0.2$ s. To reach such a period, the object needs to accrete at the observed rate for about $t\sim p_{\rm s}/\dot{p}\sim{\rm 1yr}$.

At the moment of the first observations of NGC 300 ULX1 in 2010, the NS should have had a period even larger than $45$ s measured in 2016. The question then arises why the source had such a large period to start with. This period may have originated from a long episode of a very low-rate wind accretion if we assume that this NS was born rotating much faster. To estimate the minimal time scale of the spin-down to such a long period, we use equation (30), neglect any spin-up torque acting on the NS and parametrize the spin-down torque as $K_{\rm sd}=\kappa_{t}{\mu^{2}}/R^{3}_{\rm co}$ (Lipunov, 1992). The time required to spin down to the observed period is about $40\,000\,\mu_{30}^{-2}$ years. It 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 (Villar et al., 2016): 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. All this fits well with a catastrophic event in an initially broad and 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 a couple of boundary conditions. This allows us to reconstruct the structure of the disc and obtain a relative magnetospheric radius $\xi$ which is important for describing the spin evolution of magnetized NSs as well as 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 ($\mu\sim 10^{30}{\rm\,G\,cm^{3}}$, $\dot{m}\lesssim 10$), the accretion disc remains thin and gas-pressure-dominated that implies a nearly classical scaling $R_{\rm in}\propto\dot{m}^{-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_{\rm in}\propto\dot{m}^{-(0.05-0.1)}$. For ULXPs, where the mass accretion rates reach $\dot{m}\sim 10^{3}-10^{4}$, 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_{\rm in}\propto\dot{m}^{-(0.2-0.3)}$. However, if the magnetic field is one-two orders of magnitude larger 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 $\dot{m}$.

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 time scales at the inner rim of the disc. Another test is spin period dynamics. We can compare the observed $\dot{p}$ with theoretical predictions as it was done here for NGC 5907 X-1 and NGC 300 ULX1. The constraints we get from 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 several), confirming that ULXPs are intrinsically very luminous objects rather than X-ray sources whose luminosity is amplified by an order of magnitude or more by anisotropy of their emission.

Acknowledgements

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 support from the Russian Science Foundation grant 14-12-00146 (AC, GL, PA) and from the Väisälä foundation (AC). AC and PA thank lecturers and participants of the Astrosoma summer school333http://astrosoma.ru who have provided an excellent environment for the development of this paper.

Appendix A Disc vertical structure

For the vertical structure of the disc, we assume the form $\rho=\rho_{\rm c}(1-x^{2})^{n}$, where $x=z/H$, and $n$ may be viewed as an effective vertical polytropic index. The vertical hydrostatic relation

[TABLE]

implies a similar law for the vertical pressure profile $P=P_{\rm c}\left(1-x^{2}\right)^{n+1}$, and a scaling relation for the disc thickness

[TABLE]

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

[TABLE]

[TABLE]

Here

[TABLE]

A similar model for the vertical structure was considered by Kato et al. (2008). In the radiation-pressure-dominated regime $T\propto P^{1/4}$, therefore $T\propto(1-x^{2})^{(n+1)/4}$. This approximation is not valid near the disc surface, because $T_{\rm eff}\neq 0$. The vertical radiative energy flux is determined by the vertical radiation diffusion equation:

[TABLE]

where $\epsilon=aT^{4}$ is the radiation energy density and $D=c/(3\kappa\rho)$ is the diffusion coefficient. We will denote the total energy release per unit surface area by $Q^{+}$, and the total energy leaving the two sides of the disc by

[TABLE]

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

[TABLE]

Taking into account equations (36) and (37), this expression implies

[TABLE]

and hence

[TABLE]

The midplane gas pressure is $P_{\rm g}=k\,T_{\rm c}\,\rho_{\rm c}/\tilde{m}$. The gas-to-total pressure ratio then equals to

[TABLE]

where $\tilde{m}$ is the mean particle mass (for completely ionized gas of solar metallicity $\tilde{m}\approx 0.6m_{\rm p}$).

Appendix B Derivation of advection equations

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

[TABLE]

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

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

with $\beta$ being the gas-to-total pressure ratio introduced in equation (44). Vertical integration of equation (11) allows us to write the advective flux as

[TABLE]

where the dimensionless constants are

[TABLE]

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

[TABLE]

Appendix C Dimensionless notation and equations

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

[TABLE]

The radius and the disc thickness $r=R/R_{\rm g}$ and $h=H/R_{\rm g}$ are measured in units of the gravitational radius $R_{\rm g}$. The angular frequency is normalized by the local Keplerian frequency as

[TABLE]

The characteristic magnetic moments of NSs lie in the range $10^{28}-10^{32}$ G cm3, hence we normalize $\mu$ as

[TABLE]

The mass accretion rate is normalized by the Eddington value as

[TABLE]

where

[TABLE]

It is convenient to express the surface density in the units of the inverse opacity $\kappa^{-1}$. This quantity has also the physical meaning of the disc vertical optical depth

[TABLE]

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

[TABLE]

For temperatures we use the following normalization

[TABLE]

where

[TABLE]

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

[TABLE]

where the dimensionless Alfvén radius is

[TABLE]

and

[TABLE]

We also introduce the natural time unit

[TABLE]

that may be viewed as Keplerian rotation period at $R_{\rm g}$, and the dimensionless factor

[TABLE]

The physical meaning of $\chi$ is the square of the dimensionless speed of sound $(c_{\rm s}/c)^{2}$ 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 equation (1):

[TABLE]

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

[TABLE]

The ratio of the gas pressure to the total pressure is, according to equation (44),

[TABLE]

The differential equations (15)–(20) in the dimensionless form are:

[TABLE]

Taking into account the sign of $d\Omega/dR<0$, we re-write equation (20) as

[TABLE]

where the dimensionless versions of coefficients (56)–(59) are

[TABLE]