Super-Eddington accretion onto a magnetized neutron star

TL;DR

This paper models super-Eddington accretion onto magnetized neutron stars, revealing how disc thickness and radiation influence magnetosphere size and magnetic field estimates, with implications for ultraluminous X-ray pulsars.

Contribution

It introduces a nearly-standard accretion disc model interacting with the magnetosphere, accounting for disc thickness and irradiation effects, refining magnetic field estimates in neutron star accretors.

Findings

Magnetosphere size scales with disc thickness and radiation effects.

Radiation-pressure-dominated discs have magnetosphere size independent of accretion rate.

Accounting for disc and irradiation effects reduces neutron star magnetic field estimates.

Abstract

Most of ultraluminous X-ray sources are thought to be objects accreting above their Eddington limits. In the recently identified class of ultraluminous X-ray pulsars, accretor is a neutron star and thus has a fairly small mass with a small Eddington limit. The accretion disc structure around such an object affects important observables such as equilibrium period, period derivative and the size of the magnetosphere. We propose a model of a nearly-standard accretion disc interacting with the magnetosphere only in a thin layer near the inner disc rim. Our calculations show that the size of the magnetosphere may be represented as the classical Alfv\'en radius times a dimensionless factor $\xi$ which depends on the disc thickness only. In the case of radiation-pressure-dominated disc, the size of the magnetosphere does not depend on the mass accretion rate. In general, increasing the disc thickness leads to a larger magnetosphere size in units of the Alfv\'en radius. For large enough mass accretion rates and magnetic moments, it is important to take into account not only the pressure of the magnetic field and the radiation pressure inside the disc, but also the pressure of the radiation produced close to the surface of the neutron star in accretion column. The magnetospheric size may increase by up to factor of two as a result of the effects related to the disc thickness and the irradiation from the central source. Accounting for these effects reduces the estimate of the neutron star magnetic moment by a factor of several.

Tables

Table 1. Table 1: Magnetic moments μ 30 subscript 𝜇 30 \mu_{30} for different values of η 𝜂 \eta and α 𝛼 \alpha . All the calculations were made for m ˙ = 500 ˙ 𝑚 500 \dot{m}=500 ( L = 10 40 𝐿 superscript 10 40 L=10^{40}  erg  s -1 ) and p s = 1.37 subscript 𝑝 s 1.37 p_{\rm s}=1.37  s, aimed to reproduce the properties of ULX-pulsar M82 X-2.

	$0.01$	$0.05$	$0.1$	$0.5$	$1$
$0$	161	72	51	22	16
$0.1$	157	62	36	–	–
$0.2$	153	52	10	–	–

Equations

$$\zeta \Omega R_{\mathrm{G}}^2>\sqrt{GM{R}_{\mathrm{in}}},$$

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

$$H\simeq \frac{3}{8\text{uppi}}\frac{\kappa \dot{M}}{c},$$

$$\frac{\text{rotatebox[origin=t]20.0scalebox0.95[1.0]∂}\text{v}}{\text{rotatebox[origin=t]20.0scalebox0.95[1.0]∂}t}+(\text{v}\cdot \nabla )\text{v}=-\frac{1}{\rho }\nabla P-\nabla \phi +\text{N},$$

$$v_R\frac{\text{rotatebox[origin=t]20.0scalebox0.95[1.0]∂}{v}_R}{\text{rotatebox[origin=t]20.0scalebox0.95[1.0]∂}R}+v_z\frac{\text{rotatebox[origin=t]20.0scalebox0.95[1.0]∂}{v}_R}{\text{rotatebox[origin=t]20.0scalebox0.95[1.0]∂}z}-{\Omega }^{2}R=-\frac{1}{\rho }\frac{\text{rotatebox[origin=t]20.0scalebox0.95[1.0]∂}P}{\text{rotatebox[origin=t]20.0scalebox0.95[1.0]∂}R}-\frac{GM}{R^2}.$$

$$\frac{v_{R,z}^2}{R}\sim {\left(\alpha {\left(\frac{H}{R}\right)}^2\right)}^2{\Omega }_{\mathrm{K}}^{2}R\ll \frac{1}{\rho }\frac{\mathrm{d}p}{\mathrm{d}r}\sim {\left(\frac{H}{R}\right)}^2{\Omega }_{\mathrm{K}}^{2}R.$$

$${\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}.$$

$$P_{\mathrm{c}}={\rho }_{\mathrm{c}}\frac{GM}{R^3}\frac{H^2}{4},$$

$$H=\sqrt{\frac{5\Pi }{\Sigma }\frac{R^3}{GM}}.$$

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

$$\dot{M}\frac{\mathrm{d}(\Omega R^2)}{\mathrm{d}R}=\frac{\mathrm{d}}{\mathrm{d}R}(2\text{uppi}{R}^2{W}_{r\varphi }).$$

$$F=-D{\nabla }_zϵ=-D\frac{\mathrm{d}ϵ}{\mathrm{d}z},$$

$$D=c/(3\kappa \rho )$$

$$\frac{\mathrm{d}F}{\mathrm{d}z}=\alpha PR\left|\frac{\mathrm{d}\Omega }{\mathrm{d}R}\right|,$$

$$P=nkT+\frac{a{T}^4}{3},$$

$$T_{\mathrm{s}}^4=T_{\mathrm{c}}^4-\frac{73}{120}\frac{1}{a}\frac{\alpha \kappa {\rho }_{\mathrm{c}}{P}_{\mathrm{c}}R{H}^2}{c},$$

$$2{\sigma }_{\mathrm{SB}}{T}_{\mathrm{eff}}^4=-\alpha \Pi R\frac{\mathrm{d}\Omega }{\mathrm{d}R},$$

$$2{\sigma }_{\mathrm{SB}}{T}_{\mathrm{c}}^4=-R\frac{\mathrm{d}\Omega }{\mathrm{d}R}{W}_{r\varphi }\left[\frac{219}{1024}\kappa \Sigma +1\right].$$

$$\frac{15}{16}\frac{W_{r\varphi }}{\alpha H}=\frac{a{T}_{\mathrm{c}}^4}{3}+\frac{3}{2}\frac{\Sigma k{T}_{\mathrm{c}}}{H{m}_{\mathrm{p}}},$$

$$\dot{M}({\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).$$

$$W_{r\varphi }=\frac{\dot{M}}{2\text{uppi}{R}^2}(\Omega R^2-{\Omega }_{\mathrm{in}}{R}_{\mathrm{in}}^2)+\frac{R_{\mathrm{in}}^2}{R^2}{W}_{r\varphi }^{\mathrm{in}}.$$

$$\begin{array}[]{l}\dfrac{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}\Pi}{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}R}=\displaystyle\frac{1}{\alpha}\displaystyle\frac{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}W_{r\phi}}{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}R}\displaystyle=\frac{1}{\alpha}\left(\displaystyle\frac{\dot{M}\Omega_{\rm in}R^{2}_{\rm in}}{\uppi R^{3}}+\displaystyle\frac{\dot{M}}{2\uppi}\displaystyle\frac{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}\Omega}{\text{\rotatebox[origin={t}]{20.0}{\scalebox{0.95}[1.0]{$\partial$}}}R}-2\frac{R^{2}_{\rm in}}{R^{3}}W^{\rm in}_{r\phi}\right).\\ \end{array}$$

$${\Omega }^{2}R=\frac{1}{\Sigma \alpha }\left(\frac{\dot{M}{\Omega }_{\mathrm{in}}{R}_{\mathrm{in}}^2}{\text{uppi}{R}^3}+\frac{\dot{M}}{2\text{uppi}}\frac{\text{rotatebox[origin=t]20.0scalebox0.95[1.0]∂}\Omega }{\text{rotatebox[origin=t]20.0scalebox0.95[1.0]∂}R}-2\frac{R_{\mathrm{in}}^2}{R^3}{W}_{r\varphi }^{\mathrm{in}}\right)+\frac{GM}{R^2}.$$

$$\frac{\mathrm{d}\Omega }{\mathrm{d}R}=-\frac{T_{\mathrm{c}}^4}{R}\frac{2{\sigma }_{\mathrm{SB}}}{W_{r\varphi }}{\left(\frac{219}{1024}\kappa \Sigma +1\right)}^{-1}.$$

$$T_{\mathrm{c}}^4+\frac{9}{2}\frac{\Sigma k}{a{m}_{\mathrm{p}}H}{T}_{\mathrm{c}}-\frac{45}{16}\frac{W_{r\varphi }}{\alpha aH}=0.$$

$$m=\frac{M}{1.4{\mathrm{M}}_{\odot }}.$$

$$R_{\mathrm{g}}=\frac{GM}{c^2}.$$

$$\omega =\frac{\Omega }{\sqrt{GM/R^3}}.$$

$${\mu }_{30}=\frac{\mu }{{\mu }_0}=\frac{\mu }{10^{30}\mathrm{G}\mathrm{c}{\mathrm{m}}^3}.$$

Full text

††pagerange: Super-Eddington accretion onto a magnetized neutron star–A††pubyear: 2017

1 Introduction

Magnetospheric accretion is a very important topic in astrophysics, especially for magnetized neutron stars. Interaction of falling matter with magnetosphere determines observational properties of neutron stars such as their equilibrium period and the period derivative $\dot{p}$. At low mass accretion rates, matter enters the magnetosphere through a quasi-spherical envelope (Shakura et al., 2012). When the angular momentum of the matter is large enough, accretion disc is likely to be formed. The disc forms when the angular momentum is larger than the Keplerian angular momentum at the magnetospheric radius $R_{\rm in}$ (Illarionov & Sunyaev, 1975):

[TABLE]

where $\zeta$ is a numerical factor that depends on the structure of the wind, $R_{\rm G}$ is the Bondi radius, $M$ is the neutron star mass. In the disc case, the neutron star luminosity is determined by the mass loss of the star. For the case of accretion via wind, the mass accretion rate is lower than the mass loss of the companion star by a factor of $(R_{G}/a)^{2}\ll 1$, where $a$ is the binary separation. The Roche-lobe overflow in high-mass X-ray binaries (HMXB) is a relatively rare phenomenon (the only persistent Galactic source of this kind is probably SS 433, see Fabrika 2004). HMXB with Roche-lobe overflow have higher luminosities due to a larger accretion rate and can be easily observed in other galaxies.

There are two approaches to magnetospheric accretion. The first one considers diamagnetic, finite-conduction disc threaded by stellar magnetic fields over a wide range of radii. The pioneers of this approach were Ghosh et al. (1977); other works based on this idea by Kluzniak & Rappaport (2007) and Wang (1987) used slightly modified magnetic field distribution in the disc as well as non-Keplerianity near the inner edge of the disc. The second approach is to consider a disc interacting with the stellar magnetic field only in a very narrow region near the inner boundary of the disc (see, for example, recent simulations by Parfrey et al. 2016a, b). One of the first works by Scharlemann (1978) assumes diamagnetic currents that are modulated by a single ring current located at the edge of the magnetosphere. Further developing this approach Aly (1980) studied the influence of non-radial magnetic forces. The instabilities due to velocity differences at the disc-magnetosphere boundary were studied by Anzer & Boerner (1980). Another effect of the angular velocity difference is opening of some magnetic lines considered by Aly & Kuijpers (1990), that can produce outflows along the open field lines (Shu et al., 1994; Lovelace et al., 1995).

Infalling matter creates a hot polar cap or an accretion column (Basko & Sunyaev, 1976) that radiates X-ray emission at several keV. The structure of the accretion column, however, depends on how the matter enters the magnetosphere and, specifically, upon the radius of the magnetosphere $R_{\rm in}$. Magnetospheric radius can be estimated from equality of the disc ram pressure and the magnetic pressure:

[TABLE]

where $R_{\rm A}$ is known as Alfvén, or Alfvenic, radius, $\mu$ is the neutron star magnetic moment and $\dot{M}$ is the accretion rate. The dimensionless coefficient $\xi$ is determined by the accretion flow structure. For instance, for spherical accretion, $\xi=1$. Analytical models of disc accretion onto magnetized neutron star predict $\xi$ from $\simeq 0.5$ (Ghosh et al., 1977; Kluzniak & Rappaport, 2007; Koenigl, 1991) to $1$ (Wang, 1996). MHD simulations for small magnetospheres of young stars and white dwarfs give $\xi=0.4-0.5$ (Long et al., 2005; Bessolaz et al., 2008). Recent simulations by Kulkarni & Romanova (2013) show variable $\xi\simeq 0.77({\mu^{2}}/{\dot{M}})^{-0.086}$ that may be understood as different scaling of the magnetosphere size with the basic parameter set, roughly as $R_{\rm in}\propto(\mu^{2}/\dot{M})^{1/5}$. This is similar to the predictions of some analytical models like that of Spruit & Taam (1990).

Until recently, the discs in all the X-ray pulsars were well understood as thin, gas-pressure dominated, and nearly Keplerian, and trimmed from the inside by the neutron star magnetosphere. Even for high mass accretion rates the size of the magnetosphere is usually large enough for the accretion disc to remain very far from the local Eddington limit. The accretion disc thickness at high enough mass accretion rates is determined by radiation pressure (the so-called zone A, see Shakura & Sunyaev 1973) and equals to

[TABLE]

where $\kappa=0.34$ cm2g*-1* is the (Thomson scattering) opacity. At some radius known as spherization radius (Shakura & Sunyaev, 1973), the disc thickness becomes larger than the disc radius. Even if the Eddington luminosity $L_{\rm Edd}={4\uppi GMc}/{\kappa}$ is exceeded by a considerable amount, energy release in the accretion disc is likely to remain insufficient for super-Eddington accretion unless spherization radius becomes larger than the radius of the magnetosphere. For this to take place, the total luminosity must exceed the Eddington value by a factor of about $R_{\rm in}/R_{*}$ (where $R_{*}$ is the neutron star radius) that is of the order of hundreds or thousands for X-ray pulsars. Therefore, here we will concentrate on the quite likely case when the total luminosity is above the Eddington limit but the disc is everywhere sub-critical.

Column luminosity is not strictly limited by the classical Eddington limit because of geometrical and opacity effects (Basko & Sunyaev, 1976; Mushtukov et al., 2015). It is quite possible for a neutron star to exceed the Eddington limit by a factor of $\sim 10^{2}-10^{3}$. Such objects can shine very bright in the X-ray range and explain at least some part of the population of ultraluminous X-ray sources (ULXs) in nearby galaxies.

Recently several ULXs have been shown to have coherent pulsations at periods $\sim 1$ s (Bachetti et al., 2014; Israel et al., 2017a, b), that identifies these objects with neutron stars rather than black holes. The high luminosities ($\sim 10^{40}-10^{41}{\rm erg\,s^{-1}}$) of these sources is well in excess of the Eddington limit for a neutron star that clearly points to the importance of radiation pressure leading to a high disc thickness. Most of the radiation in super-Eddington pulsars is released in the accretion column. The accretion disc, at the same time, becomes illuminated by the radiation flux exceeding the local Eddington limit, and thus its structure as well as the position of the disc-magnetosphere interface should be affected significantly. Both radiation from the central source (accretion column) and the magnetic field of the neutron star create a radial pressure gradient within the disc, violating one of the basic assumptions of the standard disc theory.

In this paper we develop a model of the accretion disc with deviations from the Keplerian rotation caused by the disc thickness and by the influence of magnetic field pressure at its inner edge. We assume that the disc-magnetosphere interaction is confined to a thin layer at the inner boundary of the disc. We also study the influence of radiation pressure from accretion column on the disc structure and the magnetospheric radius for ULX-pulsars.

In Section 2 we describe our model. Techniques used to solve the equations are given in Section 3. Our main results are presented in Section 4 and discussed in Section 5.

2 The model

In this section we give the whole system of equations that describes the structure of the disc. We start from the equations of hydrodynamics. The momentum equation in the vector notation takes the form:

[TABLE]

where $\rho$ is density, $P$ is pressure, and v is the velocity field in the disc and $\varphi$ is gravitational potential. The last viscous term N does not affect the radial and vertical components of the Euler equation. The azimuthal component of the equation is equivalent to the equation of angular momentum transfer that will be considered separately. In steady state and axial symmetry, the radial component of this equation in cylindrical coordinate system can be rewritten as:

[TABLE]

Here $P$ is the pressure, $v_{R}$ is the radial and $v_{z}$ is the vertical velocity and $\Omega$ is the angular velocity in the disc. Radiation pressure of the central source shifts the radial force equilibrium of the Keplerian disc ($\Omega^{2}R=GM/R^{2}$), and the radial pressure gradient can no more be neglected in the inner parts of the disc. This leads to two effects: the angular velocity of the disc becomes non-Keplerian, and the inner parts of the disc are additionally decelerated in radial direction. Thus in this case, the radial velocity is even smaller than in the standard disc, allowing us to neglect the first term, quadratic in $v_{R}$, on the left-hand side of equation (5). Because we expect $v_{z}$ to be also small, the second term can be neglected as well. In general, in thin accretion discs, first two terms in equation  (5) are of the order (see e.g. Urpin, 1984)

[TABLE]

Integrating equation (5) over the vertical coordinate $z$, or, more precisely, over $\rho{\rm d}z$, yields:

[TABLE]

Here $\Pi=\int^{H}_{-H}{p{\rm d}z}$ is the vertically integrated pressure, $\Sigma=\int^{H}_{-H}\rho{\rm d}z$ is the surface density, and $H$ is the thickness of the disc.

To relate the vertically-integrated and equatorial-plane parameters of the disc, it is important to make assumptions about the vertical structure of the disc, which is discussed in detail in Appendix A. From hydrostatic equilibrium together with the assumptions about the vertical structure, we have the central pressure in the disc:

[TABLE]

where $\rho_{\rm c}$ is the central (mid-plane) density. The vertical structure together with equation (8) provide us with the thickness of the disc:

[TABLE]

Following the standard accretion disc theory we assume alpha-prescription (Shakura & Sunyaev, 1973):

[TABLE]

where $W_{r\phi}$ is the $r\phi$-component of the vertically integrated viscous stress tensor. Angular momentum conservation equation is:

[TABLE]

We assume $\dot{M}$ constant with radius that allows straightforward integration of this equation. Following Shakura & Sunyaev (1973), we assume that the energy released in the accretion disc is radiated locally from its surface. In this case, the local energy release will be coupled to the radiation emitted from the disc surface and to the angular velocity gradient in the disc allowing to link the rotation law of the disc $\omega(R)$ with the vertically-integrated pressure. The pressure is itself linked to the energy release rate by the vertical radiation diffusion equation

[TABLE]

where $F$ is the vertical energy flux, $\epsilon=aT^{4}$ is the radiation energy density and

[TABLE]

is the diffusion coefficient. Local energy dissipation rate is

[TABLE]

where $P$ is the total (gas+radiation) pressure:

[TABLE]

where $n\simeq\dfrac{2\rho}{m_{\rm p}}$ is the proton or electron number density for proton-electron plasma and $a$ is the radiation constant. Solving equations (12) and (14) one obtains the following relation between the surface temperature $T_{\rm s}$ and the central temperature $T_{\rm c}$ :

[TABLE]

where the dimensionless coefficient is derived in Appendix A. Assuming that the surface temperature $T_{\rm s}$ equals the effective temperature determined by the overall energy release in the disc,

[TABLE]

we get an expression for the central disc temperature:

[TABLE]

This expression differs from the similar equation for the standard disc (see equation 2.24 in Shakura & Sunyaev 1973) because we account differently for the vertical structure of the disc.

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

[TABLE]

where we have taken into account $\alpha$-prescription and the relations between the central and vertically-integrated quantities implied by the adopted vertical structure (see Appendix A).

Close to the disc boundary, the infalling matter has angular velocity that differs from the angular velocity of the neutron star. Entering the magnetosphere, the accreting material should somehow lose the excess of angular momentum. Following Lai (2014) and Aly & Kuijpers (1990), we assume that this excess angular momentum is removed by the magnetic and radiation stresses at the magnetospheric boundary:

[TABLE]

where $k_{\rm t}={B_{\phi}}/{B_{z}}$ is the ratio of magnetic field components and $L=\eta\dot{M}c^{2}$ is the luminosity released close to the neutron star, in the accretion column with efficiency $\eta$. Parameter $k_{\rm t}$ can vary from $k_{\rm t}=0$ (purely poloidal magnetic field in the region near corotation radius) to $k_{\rm t}\simeq 1$ (Lai, 2014). Some authors tend to use $k_{\rm t}=1/3$ (Lipunov et al., 1992). In this work we use $k_{\rm t}=0.5$. In Section 4.3, the effects of different $k_{\rm t}$ will be considered.

At the boundary, internal disc pressure is balanced with the external pressure of the radiation source $P_{\rm rad}={L}/{4\uppi R^{2}c}$ and the magnetic field pressure $P_{\rm mag}={\mu^{2}}/{8\uppi R^{6}}$. Note that the radiation pressure can be even higher by a factor up to two due to reflection of photons, depending on the scattering albedo of the disc. Thus vertically integrated pressure at $R_{\rm in}$ is $\Pi=(P_{\rm rad}+P_{\rm mag})H_{\rm in}$. The $\alpha$-prescription (10) gives us the vertically integrated stress at the boundary:

[TABLE]

3 Solving the disc structure

3.1 The basic system of equations

Integrating equation (11), we obtain a general expression for the vertically-integrated viscous stress:

[TABLE]

Using equations (22) and (10) we obtain the radial gradient of vertically-integrated pressure

[TABLE]

Substituting expression (23) into the radial Euler equation (7), we have the following expression linking deviations from Kepler’s law and the pressure gradient expanded from the previous equation:

[TABLE]

This equation is solved for $\Sigma$ at each radius. We still need one more equation to calculate the derivative of the rotation frequency and the unknown $\Omega(R)$. Expressing ${\rm d}\Omega/{\rm d}R$ from equation (18) yields:

[TABLE]

This is the only differential equation we need to solve. All the quantities in its right-hand side are found locally by solving algebraic equations: $W_{r\phi}$ is found from equation (22), $\Sigma$ from equation (24), and the central temperature can be calculated from the equation of state in the equatorial plane following from equation (19):

[TABLE]

We solve this quartic equation for $T_{\rm c}$ using Ferrari’s formula, see Hirose et al. (2009) for details.

3.2 Dimensionless equations

In this section we give all the equations in dimensionless form, as they were used to calculate the disc structure. Here we list the dimensionless parameters and combinations we use throughout the paper. We normalize the neutron star mass as

[TABLE]

The radii and disc thicknesses $r=\displaystyle\frac{R}{R_{\rm g}}$ and $h=\displaystyle\frac{H}{R_{\rm g}}$ are measured in units of the gravitational radius

[TABLE]

The angular frequency is normalized by the local Keplerian frequency as

[TABLE]

The characteristic magnetic moments of neutron stars 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 vertical optical depth of the disc

[TABLE]

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

[TABLE]

For the temperature, 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]

and the dimensionless factor

[TABLE]

From equation (21) we can find the stress tensor at the boundary of the disc

[TABLE]

The angular velocity at the inner boundary can be found from equation (20):

[TABLE]

where $p_{\rm s}$ is the neutron star period in seconds. A general expression for the vertically-integrated viscous stress is then obtained from equation (22):

[TABLE]

The relative thickness of the disc follows from equation (9):

[TABLE]

The relative thickness of the disc at the inner boundary can be obtained from equations (42) and (45)

[TABLE]

From equations (24) and (25) we find the equation for the dimensionless angular velocity

[TABLE]

which is solved numerically for $\tau$. The temperature is obtained using equation (26) which in dimensionless form can be written as

[TABLE]

The only differential equation (25) in dimensionless form becomes

[TABLE]

3.3 Boundary conditions

Two boundary conditions on the inner boundary of the disc are given by equations (42) and (43). The first one arises from equation (20) that says that the excess of angular momentum is removed by magnetic and radiation stresses acting upon the magnetospheric boundary. If the neutron star and the disc are in corotation, $\Omega_{\rm ns}=\Omega_{\rm in}$, the left-hand side of equation (20) is zero, and the thickness of this disc should be zero as in the standard disc with a zero-stress boundary condition. If the neutron star rotates faster than the disc, the left-hand side of equation (20) becomes negative while the right-hand side is always positive. Therefore, stationary accretion becomes impossible. This situation is known as propeller regime.

Boundary conditions can be solved together with equations (47) and (48). There are 5 unknowns in this system of equations: angular velocity $\omega_{\rm in}$, magnetospheric radius $r_{\rm in}$, disc thickness $h_{\rm in}$, surface density (optical depth) $\tau_{\rm in}$ and the central disc temperature $t_{\rm in}$ for a given accretion rate, magnetic field and the spin period. This means that for any $\xi=r_{\rm in}/r_{\rm A}$, the inner radius is specified and, thus, we can calculate the remaining 4 parameters at the inner boundary of the disc.

There are two terms in the right-hand side of equation (20): the first one depends on the magnetic field and accretion rate and the second one depends on the spin period of the neutron star. The relative importance of the two terms is set by the ratio of the neutron star period and the equilibrium period

[TABLE]

This period is determined from the equality of the Alfvén and the corotation radii.111Note that this definition of the equilibrium period does not take into account factor $\xi$ in equation (2) for the Alfvenic radius. If the spin period of a neutron star is less than $\xi^{3/2}p_{\rm eq}$, magnetosphere rotates faster than the disc and the matter from the disc cannot accrete onto the neutron star. Accretion is possible only if the spin period of the neutron star is longer than equilibrium multiplied by $\xi^{3/2}$. If $p_{\rm s}\leq p_{\rm eq}\xi^{3/2}$, condition (43) cannot be satisfied, and the neutron star enters propeller regime. It is thus natural to expect the structure of the boundary conditions and the disc itself to depend primarily upon the ratio $p_{\rm s}/p_{\rm eq}$, or on the ratio $R_{\rm co}/R_{\rm A}=\left(p_{\rm s}/p_{\rm eq}\right)^{2/3}$. At the same time, it is interesting to consider the role of the magnetic field and of the mass accretion rate at a fixed value of $p_{\rm s}/p_{\rm eq}$.

The properties of our boundary conditions are illustrated by Figs 1–3. In Fig. 1, we consider the influence of the mass accretion rate on the relative thickness at the inner disc radius. For this we fixed the magnetic moment $\mu$ and the spin period in $p_{\rm eq}$ units, $p_{\rm s}=10p_{\rm eq}$, and considered different accretion rate values. The value of $\xi$ is still unknown, but most of the curves have a single maximum of $(H/R)_{\rm in}$ that is the upper limit for the possible disc thickness at the inner boundary.

The relative thickness at $r_{\rm in}$ calculated using our set of boundary conditions is plotted for different accretion rates as a function of $\xi$ in Fig. 1 for $\mu_{30}=0.1$ (left panel) and $\mu_{30}=100$ (right panel). The accretion rate changes in the range $\dot{m}=0.1-1000$ from the blue (bottom) to the red (top) curve. In case of a small magnetic field, the thickness of the disc rapidly grows with the accretion rate and for high accretion rates it becomes comparable with the radius. When magnetic field is high, the disc stays thin even for high accretion rates.

Fig. 2 shows similar curves for constant spin period in equilibrium units $p_{\rm s}=10p_{\rm eq}$ and accretion rate $\dot{m}=10$ for different magnetic moments $\mu_{30}=0.1-100$. High magnetic fields keep the disc thin. For all the parameter values, disc thickness tends to be strictly limited at small $\xi$. Magnetic stress at the boundary is proportional to $r_{\rm in}^{-6}$ and increases very rapidly at small $\xi$. Thus, for small $\xi$, even a small thickness of the disc is sufficient to remove the excess angular momentum. That is why the disc becomes thin and stays very close to Keplerian rotation.

Finally we study the influence of the spin period on the conditions at the disc boundary. The disc thickness at the inner disc edge is shown as a function of $\xi$ for constant magnetic moment $\mu_{30}=1$ and mass accretion rate $\dot{m}=10$ and different spin periods $p_{\rm s}=(0.2-20)p_{\rm eq}$ in Fig. 3. Smallest spin periods produce very thin accretion discs as $\Omega_{\rm in}-\Omega_{\rm ns}$ tends to be very small. For large $\xi$ and small spin periods, the boundary conditions cannot be satisfied. Some of the curves disappear at large $\xi$, where equation (47) does not have solutions any more. As this happens at small periods, it is logical to identify the effect with entrance into propeller regime. An alternative explanation is that, for such parameter sets, there is no stationary solution, and the accretion disc will undergo some kind of non-stationary accretion.

3.4 Searching for the solution

The inner boundary of the disc determined by the value $\xi=R_{\rm in}/R_{\rm A}$ is, in fact, unknown but may be found by matching the boundary conditions to the solution of equation (49). Equation (49) is stable if integrated inwards. We find surface density from pressure gradient, hence we need to assume a small non-Keplerianity at the outer boundary of the disc. In our calculation we use $\omega_{\rm out}=0.99$, but this value does not influence the disc structure anywhere but in a very small region near the outer disc radius. The disc structure rapidly evolves toward a solution independent of the starting value of $\omega$.

We start our calculations by solving the boundary conditions. For a given $\xi$ we solve equations (42), (43), (47) and (48). After that we start integrating the radial structure of the disc from the outer boundary. For a given point we calculate viscous stresses $w_{r\phi}$ using equation (44). Then we substitute it to equations (46), (47), and (48) and solve these equations together to find the surface density $\tau$ of the disc, disc thickness $h$ and central temperature $t_{\rm c}$. We substitute these quantities into equation (49) and find the increment of angular velocity. Then we solve differential equation (49) with a leapfrog integration scheme, repeating all the steps listed above in this paragraph at each radius.

To find the disc boundary we use the shooting method: first we assume some $\xi$, using it we calculate $\omega_{\rm in}$ with the boundary condition (43) and then calculate the whole structure of the disc from the outside. We then vary $\xi$ until the angular velocity at the boundary is equal to $\omega_{\rm in}$ calculated from the boundary condition. We also require the surface density and disc thickness calculated solving the disc structure equations to coincide with the corresponding quantities at the boundary. Sometimes, this additional requirement allows to exclude spurious solutions. This intersection is found using bisection method. We search for $\xi$ in the range $0.2-2$ that usually contains one solution. The value of $\xi$ is found with accuracy $\sim 10^{-3}$. In some cases, for example, for $p\lesssim p_{\rm eq}$ or for very thick, supercritical discs (see next section) equation (47) does not have any solutions.

4 Results

In this Section, we consider the results of our simulations for different parameters $\mu_{30}$, $\dot{m}$, $p_{\rm s}$, $\eta$, $\alpha$ and $k_{\rm t}$. While the first three are likely to vary from one object to another and even for a single object, the other three are poorly known mainly due to unknown details of the physics involved. The effective value of $\eta$ is affected not only by the efficiency of accretion but also by the beaming and collimation of the radiation of the accretion column. The last two parameters, $\alpha$ and $k_{\rm t}$, depend on the complicated physics of magnetic field amplification and reconnection in the disc and the magnetosphere. Below, in Sections 4.1 and 4.2, we will assume $\eta=0$ (i.e. no irradiation from the accretion column), $\alpha=0.1$, and $k_{\rm t}=0.5$. We will also mainly focus on the case on slowly rotating magnetosphere, that proves to be a good approximation for large mass accretion rates. In Section 4.3, the effects of spin period and irradiation will be considered.

4.1 Disc structure

Depending on the accretion rate and the size of the magnetosphere, the disc around a neutron star can be thin and gas-pressure dominated, or it can come closer to the neutron star and be thicker and radiation-pressure-dominated in its inner parts.

The structure of a thin disc for $\mu_{\rm 30}=1$, $\dot{m}=10$ and $p_{\rm s}=10p_{\rm eq}=4.87$ s is presented in Fig. 4. The upper panel shows deviation from Keplerian rotation $1-\omega$ as a function of the radial coordinate $r=R/R_{\rm g}$. The middle panel shows the relative thickness $H/R$ of the disc. Our results are shown with the solid red line. The thickness predicted by the standard disc model within zone A (radiation-pressure-dominated zone) or zone B (gas-pressure-dominated zone) is plotted with the dashed green or dotted blue lines, respectively. We use the following expressions for the disc thickness in the radiation-dominated zone:

[TABLE]

or in a dimensionless form:

[TABLE]

In the gas-dominated zone the disc thickness is

[TABLE]

or

[TABLE]

rewriting formulas from Shakura & Sunyaev (1973) using our dimensionless variables. There is an additional factor $\sqrt{5}$ in both equations that accounts for the vertical structure we use. The correction factor $1-(r_{\rm in}/{r})^{1/2}$ arises from the boundary condition $W^{\rm in}_{r\phi}=0$ used in the standard disc model.

The bottom panel of Fig. 4 shows the optical depth as a function of radius. For the standard disc, the optical depth in the radiation-dominated zone is

[TABLE]

and in the gas-dominated zone:

[TABLE]

Far from the inner boundary, the structure of our disc is very close to the standard gas-pressure-dominated thin disc. Near the inner boundary, a different boundary condition makes the disc thicker.

At larger mass accretion rates and smaller magnetic moments, the inner disc becomes radiation-dominated. In Fig. 5, we show an example of a moderately-thick disc structure for $\mu_{\rm 30}=10$, $\dot{m}=500$ and $p_{\rm s}=10p_{\rm eq}=48.65$ s, giving $\xi=0.78$. Maximal deviations from the Keplerian rotation in this case are about $1.5$ per cent, while the maximal relative thickness is about $0.5$. Because of the boundary conditions, the disc has a non-zero thickness at $r_{\rm in}$ but $(H/R)_{\rm max}$ is smaller than that for the standard disc because the matter is accumulated near the inner boundary and the disc becomes denser and thinner because while $\Pi$ stays the same, surface density becomes higher, hence disc thickness decreases as we can see from equation (9). For such a disc it is important to account possible outflows (Poutanen et al., 2007) together with the whole complex of effects connected to the disc thickness such as radial advection and non-local energy release (Abolmasov & Chashkina, 2015).

For the parameters characteristic for normal bright X-ray pulsars ($\mu_{30}\sim 1-10$, $\dot{m}\sim 0.1-1$, see e.g. Walter et al. 2015), the standard thin disc model is well reproduced. As for the standard radiation-pressure-dominated disc, the maximal $H/R$ is reached at the radius $r\simeq(1.5-2)r_{\rm in}$. The exact position of the maximum and its value depend on parameters $k_{\rm t}$ and $\alpha$ (see Section 5).

Factor $\xi$ grows with the increasing disc thickness and becomes close to 1 when the relative thickness of the disc becomes large. As we will see below in Section 4.2, $\xi$ mainly depends on the disc thickness and the magnetic pitch parameter $k_{\rm t}$ and is relatively insensitive to other parameters.

4.2 Magnetospheric radius and disc thickness

The magnetospheric radius is found as the solution of the Cauchy problem that we obtain using the shooting method. In this section, we will consider the size of the magnetosphere and the thickness of the disc in its inner parts for different $\mu$ and $\dot{m}$. We will focus on the case of very slow rotation $p\gg p_{\rm eq}$ and $\eta=0$. The effects of spin period and irradiation will be considered in the next section.

Figs 6 and 7 show the contours of $\xi$ and relative thickness $(H/R)_{\rm in}$, respectively, in the $\mu_{30}$–$\dot{m}$ plane. All calculations are made for $p_{\rm s}=10p_{\rm eq}$, when the corotation radius is considerably larger than the radius of the magnetosphere. As long as $p_{\rm s}\gg p_{\rm eq}$, the spin period does not affect the size and the structure of the magnetosphere. The white lines in Fig. 6 correspond to the constant maximal relative thicknesses of the disc: $\left(H/R\right)_{\rm max}=0.03$, $0.1$ and $0.3$.

As we said in the previous section, the maximal relative disc thickness can be reached not at the boundary but at some distance of about $(1-2)r_{\rm in}$. This is similar to the case of the radiation-pressure-dominated standard disc where zero-stress boundary condition leads to a non-monotonic thickness behavior with a maximum at $\sim 2r_{\rm in}$. For smaller $\dot{m}\lesssim 10$, all the disc remains in gas-pressure-dominated state and the maximal thickness is reached at its outer rim. However, the value of $H/R$ in this case is practically independent of radius, except in its innermost parts (see Fig. 4). Fig. 8 shows the contours of $(H/R)_{\rm max}$ in $\dot{m}$-$\mu$ parameter space with colors.

Figs 6–8 suggest that all the three quantities: $\xi$, $(H/R)_{\rm in}$ and $(H/R)_{\rm max}$, depend on each other. This dependence is especially tight for $\xi$ and $(H/R)_{\rm in}$. It appears that, for any values of $\dot{m}$ and $\mu$, $\xi$ depends only on $(H/R)_{\rm in}$. This is illustrated by Fig. 9, showing $\xi$ as a function of $(H/R)_{\rm in}$ for the same data that were used in Figs 6 and  7. Here, the scatter around a universal dependence $\xi(H/R)$ is less than $1$ per cent. The results of the calculations are shown with black circles and the solid red line is the $\xi=2^{-3/7}[k_{\rm t}(H/R)_{\rm in}]^{2/7}$ slow-rotation asymptotic we derive below (see equation 57).

All the results discussed above were obtained for small fastness parameter $\omega_{\rm s}=\Omega_{\rm ns}/\Omega_{\rm in}\simeq(R_{\rm in}/R_{\rm c})^{3/2}$. This situation can be realized in transient objects during outbursts when the accretion rate is much higher than the average. For such objects we can neglect the second term in the boundary condition (43). The boundary condition may be then solved for $\xi$ under assumption $\omega_{\rm in}=1$ that is accurate up to the terms of the order $(H/R)_{\rm in}^{2}$:

[TABLE]

where $h_{\rm in}$ in general case should be calculated numerically. For purely gas- and radiation-dominated standard discs this limit can be obtained analytically.

For the gas-dominated disc it is safe to use $\xi$ in the general form given by equation (57) and the thickness follows equation (53) without the correction factor, because in this case the disc thickness is almost constant everywhere. Expression for the magnetospheric size becomes then

[TABLE]

In the radiation-pressure-dominated disc, the relative disc thickness tends to zero more rapidly at the boundary and we cannot neglect the correction due to boundary conditions. The relative disc thickness can be obtained by substituting equations (48), (42) and (46) to equation (47):

[TABLE]

The first factor in the right-hand side of this equation is $\sim r_{\rm A}/\alpha\gg 1$, therefore the two terms in the square brackets should nearly cancel each other because $(h/r)_{\rm in}\lesssim 1$. This gives us an estimate for $\xi$ in the radiation-pressure-dominated disc:

[TABLE]

The two expressions for $\xi$ given by equations (58) and (60) reproduce the scaling with $\mu_{30}$ and $\dot{m}$ (see Fig. 9) but cannot be used as universal approximations because the structure of the disc is modified near the inner boundary. Still they work surprisingly well for the case of a slowly-rotating magnetosphere ($p_{\rm s}\gg p_{\rm eq}$) and no irradiation ($\eta=0$). These two approximations are shown with green open circles and blue crosses in Fig. 9. Each point corresponds to a unique combination of $\mu_{30}$ and $\dot{m}$, therefore in general the abscissa and ordinate of the points are not expected to follow any single curve. However, in their applicability regions, the two approximations show very small scatter.

In radiation-dominated disc $\xi\propto\dot{m}^{2/7}$ and $r_{\rm A}\propto\dot{m}^{-2/7}$, therefore $r_{\rm in}=\xi r_{\rm A}$ does not depend on $\dot{m}$. Increasing accretion rate leads to an increase in the disc thickness in a way that the disc pressure remains constant. Hence, the pressure balance condition can be satisfied at the same radius for different $\dot{m}$.

4.3 Irradiation, rotation and free parameters

In Sections 4.1 and 4.2 we ignored the effects of irradiation by setting $\eta=0$. However, the accretion efficiency of a neutron star can be $\eta=0.1-0.2$. Radiation of the central X-ray source affects the pressure balance condition (42) as well as, through radiation drag, the other boundary condition (43). As a result of these two effects, the inner radius becomes larger, especially for thicker discs, and the disc thickness at the inner boundary also increases. The difference in $\xi$ for $\eta=0$ and $\eta=0.1$ is shown in Fig. 10. Thick discs are more sensitive to accretion efficiency because the intercepted luminosity fraction is proportional to the solid angle subtended by the disc $\propto(H/R)_{\rm in}$. As expected, the radiation does not affect the magnetospheric radius in the thin gas-dominated disc regime, but for a thick disc, radiation can increase magnetospheric radius by about $30$ per cent. Irradiation increases the inner relative thickness while $(H/R)_{\rm max}$ stays about the same for a disc with and without irradiation.

Though the slow-rotation estimate (57) was derived under the assumption $\eta=0$, it holds to a good accuracy (of about $\eta(H/R)_{\rm in}$) even when both $\xi$ and the disc thickness are strongly affected by irradiation. In Fig. 11, we show the effects of irradiation and spin period of the neutron star accreting at a very high rate ($\dot{m}=10^{3}$, $\mu_{30}=10$). The dependence of $\xi$ on the relative neutron star spin period $p_{\rm s}/p_{\rm eq}$ is shown with and without irradiation for a large mass accretion rate, when the structure of the disc is especially vulnerable to irradiation. The black circles show $\xi$ for disc without irradiation ($\eta=0$), and the red open circles are for the disc irradiated from inside ($\eta=0.1$). The blue dotted line shows $\xi_{\infty}$ calculated according to equation  (57). The slow-rotation approximation given by expression (57) works well whenever the period is much larger than the equilibrium and there is no irradiation, and gives the lowest possible value of $\xi$ for given $\mu$ and $\dot{m}$. In the radiation-dominated regime with $\eta=0$, the physical size of the magnetosphere practically does not depend on the spin period of neutron star, as it does not depend on the mass accretion rate.

In general case, in the fast-rotation limit, near corotation, disc becomes thinner, and $\xi$ increases. In Fig. 12, we show how variations in neutron star spin frequency affects the structure of the disc. Solutions exist above the poorly constrained minimal period that we interpret as the propeller limit and discuss in more detail in Section 5.3.

There are two additional free parameters in the model: viscosity of a disc $\alpha$ and the ratio of magnetic field components $k_{\rm t}=B_{\phi}/B_{z}$. Their ratio $\alpha/k_{\rm t}$ has an important impact on accretion disc structure, being the ratio of tangential components of the viscous stress inside the disc $W_{r\phi}=2\alpha H_{\rm in}\displaystyle\frac{\mu^{2}}{8\uppi R^{6}_{\rm in}}$ to the stress created by the magnetic fields in the magnetosphere $W_{r\phi}=2k_{\rm t}H_{\rm in}\displaystyle\frac{\mu^{2}}{4\uppi R^{6}_{\rm in}}$. In the limiting case when these stresses are equal, $\alpha=2k_{\rm t}$, position of the relative thickness $(H/R)_{\rm max}$ maximum in the disc reaches the inner boundary. Effects of variations in $\eta$, $\alpha$ and $k_{\rm t}$ are shown in Fig. 13 for two different magnetic moments. As it is easy to see from this figure, the maximal disc thickness is insensitive to $\eta$, a large $\alpha/k_{\rm t}\gtrsim 1$ make the inner disc thicker (maximal thickness is reached at $r_{\rm in}$), and larger viscosity generally makes the magnetospheric size larger and the disc thicker.

Though the radiative efficiency of neutron star accretion seems to lie somewhere around $\eta\simeq 0.15-0.2$, irradiation is also affected by anisotropy of the radiation emitted by the accretion column, that can add a factor of two uncertainty in the effective value of $\eta$. Besides, if some of the accreting material is expelled from the magnetosphere (that is likely if the radiation pressure is large, see Section 5.2), the mass accretion rate in the accretion column becomes smaller, this proportionally decreases the effective value of $\eta$. This justifies a large range of possible efficiencies considered in Fig. 13.

4.4 Deviations from Kepler’s law

All the disc structures discussed above show very small deviations from the Keplerian rotation because the discs are quite thin. Non-Keplerianity increases with disc thickness as $1-\omega\sim(H/R)^{2}$. According to the boundary condition equation (20), the disc thickness at the inner boundary increases when parameter $k_{\rm t}$ becomes smaller. For small $k_{\rm t}\lesssim\alpha/2$, the disc thickness strongly peaks at the inner boundary (see Fig. 13), and non-Keplerianity can become significant.

The structure of a thin disc with $k_{\rm t}=0.05$ for a neutron star with $\mu_{\rm 30}=1$, $\dot{m}=100$ and $p_{\rm s}=10p_{\rm eq}=0.9$ s is presented in Fig. 14. The difference from Keplerian rotation $1-\omega$ is shown as a function of the radial coordinate $r$. Here non-Keplerianity can reach $\sim 10$ per cent and the thickness $(H/R)_{\rm in}\simeq 1$. For even thicker discs, deviations from Kepler’s law can formally reach even larger values, but the model itself becomes inapplicable.

5 Discussion

5.1 Magnetosphere in radiation-dominated regime

A neutron star accreting at a high enough rate will have the relative magnetosphere size given by equation (60) that implies the physical size of magnetosphere

[TABLE]

independent of the mass accretion rate.

The rate of changes of the total angular momentum of a NS can be written down as follows:

[TABLE]

where $N_{\rm up}=\dot{M}\sqrt{GMR_{\rm in}}$ is the spin-up torque due to accretion, $N_{\rm down}$ is the spin-down torque. If accretion rate is large enough we can neglect the second term. In the radiation-pressure-dominated disc, $R_{\rm in}$ does not depend on the accretion rate and $|\dot{p}_{\rm s}|\propto\dot{M}$. For the gas-pressure-dominated disc, $|\dot{p}_{\rm s}|\propto\dot{M}^{427/483}$ ($427/483\simeq 6/7$ with the accuracy of about 3 per cent). These are the two limiting cases, and for real objects $|\dot{p}_{\rm s}|\propto\dot{M}^{6/7}-\dot{M}$, that is in a good agreement with observations of X-ray pulsars (Bildsten et al., 1997).

5.2 Eddington limit(s)

The challenge in understanding the super-Eddington accretion is in the large radiation fields that seemingly prevent the accretion that causes them. In the case of the accretion discs of ULX-pulsars, one should take into account two radiation fields: one from the disc itself and one from the central X-ray source (accretion column) emitting most of the luminosity. Correspondingly, we will talk about local and non-local Eddington limits. Impact of the second, non-local radiation field becomes important when its pressure becomes comparable to the magnetic field pressure

[TABLE]

If we now substitute expression (60) for $\xi$, the pressure ratio becomes

[TABLE]

For large mass accretion rates, and also large magnetic field and viscosity, irradiation from accretion column is important. Usually, the adopted picture of magnetospheric accretion assumes that matter is captured by magnetic field at the boundary of the disc. If the radiation pressure from the accretion column is larger than the magnetic pressure, optically thin infalling matter may be blown away by the radiation pressure force escaping the magnetosphere (Mushtukov et al., 2017). Optically thick accretion flow, on the other hand, can remain bound, but its dynamics will still be affected by the radiation pressure.

Let us estimate the critical (non-local Eddington) mass accretion rate when magnetic pressure equals radiation pressure. When the two pressures are equal, equation (64) implies:

[TABLE]

When accretion rate is larger than this value, our model is incomplete as some part of the accreting material is probably ejected from the magnetosphere. This limitation is plotted on Fig. 10 by the black solid line.

The local Eddington limit is reached when $(h/r)_{\rm max}=1$. In our model relative thickness is maximal at $r\simeq 2r_{\rm in}$. Let us consider a standard radiation-dominated disc far from the corotation. Using the disc thickness from equation  (51) and estimating $\xi$ from equation (60), one can obtain the maximum possible accretion rate:

[TABLE]

This limit is plotted in Figs 6–8 and 10 by the black dashed lines. For moderate magnetic moments $\mu_{30}\sim 1-10$, both limits are comparable $\dot{m}_{1}\sim\dot{m}_{2}\sim 500$, but for larger magnetic fields there is a region in parameter space where only non-local Eddington limit is violated. Unlike the local limit, $\dot{m}_{1}$ depends strongly on viscosity and practically does not depend on the magnetic field strength.

5.3 Propeller regime and ULX

The first ULX-pulsar, M82 X-2 (Bachetti et al., 2014), has a high enough luminosity of $L_{\rm X}\sim 10^{40}$ erg s*-1* for its inner parts of the disc to be radiation-pressure-dominated. Tsygankov et al. (2016) assumed that M82 X-2 is very close to corotation (see also Bachetti et al. 2014, Lyutikov 2014, Dall’Osso et al. 2015 for more discussion on this topic), and we observe it in two states: the low luminosity corresponds to the propeller regime and the high luminosity to the accretion regime. Therefore, the observed period of $p_{\rm s}=1.37$ s is the minimal possible period of a neutron star with a given magnetic field when accretion is still possible.

In the propeller regime, the angular velocity at the inner boundary of the disc is equal to the angular velocity of the neutron star. In this case, condition (20) implies $h_{\rm in}=0$. In the radiation-pressure-dominated regime, it leads to expression (59) that is approximately valid whenever the disc is radiation-pressure-dominated and irradiation is negligible. Usually, the size of the magnetosphere is not known directly, but we may assume $\omega_{\rm in}=1$ and thus $p_{\rm s}=p_{*}r_{\rm in}^{3/2}$. This allows to estimate the magnetic field as a function of the (accurately measured) spin period and (less constrained) mass accretion rate and efficiency

[TABLE]

If the accretion column radiates isotropically, $\eta\dot{m}=L/L_{\rm Edd}$. In a more general anisotropic case, $\eta$ is not bound to be equal to the apparent radiative efficiency of neutron star accretion. Isotropic efficiency is well approximated by

[TABLE]

At some luminosity, expression in the brackets in equation (67) becomes zero. This luminosity coincides with the non-local Eddington limit (65) if the disc is in co-rotation $p_{\rm s}=p_{*}r_{\rm in}^{3/2}$. For case of M82 X-2 ($p_{\rm s}=1.37$ s, $\dot{m}=500$ corresponding to the isotropic luminosity $L=10^{40}$ erg s*-1*), $\alpha=0.1$ and $k_{\rm t}=0.5$, equation (67) yields $\mu\simeq 5.1\times 10^{31}$G cm3 if irradiation is turned off. The magnetospheric radius in this case corresponds to $\xi\simeq 0.7$. Setting $\eta=0.1$ decreases the required magnetic field to $\mu\simeq 3.7\times 10^{31}$ G cm3 and $\xi\simeq 0.8$. These magnetic moment values correspond, respectively, to magnetic fields at the pole of $B=10^{14}$ G and $B=7.4\times 10^{13}$ G for the neutron star radius of $R_{*}=10$ km. The magnetic field value estimated without irradiation is in a good agreement with that of Tsygankov et al. (2016). As we see, irradiation from the accretion column decreases magnetic field by a factor of about 1.4, if a standard set of auxiliary parameters is used. However, the limit (67) has strong dependence on viscosity, and may be used to get an upper limit for $\alpha$. It gives $\alpha\lesssim 0.21$ for efficiency $\eta=0.1$ and $\alpha\lesssim 0.1$ for $\eta=0.2$. The estimated magnetic moments for different $\eta$ and $\alpha$ are presented in Table 1.

Another ULX-pulsar, NGC 7793 P13 (Israel et al., 2017b), has a period of $p_{\rm s}=0.42$ s and a luminosity of $L_{\rm X}=5\times 10^{39}$ erg s*-1*. Assuming that the observed period is the minimal period the magnetic field is $\mu_{30}=7\times 10^{30}$ G cm3 and $\xi\simeq 0.7$ in the case without irradiation and $\mu_{30}=8.8\times 10^{30}$ G cm3 and $\xi\simeq 0.7$ for the case with irradiation. The maximal relative thickness of the disk is $(H/R)_{\rm max}=0.2$, the disc is nearly standard. We see that for this object irradiation is not so important.

The third ULX-pulsar in NGC 5907 (Israel et al., 2017a) has a very high apparent luminosity of $L_{\rm X}=1.5\times 10^{40}-10^{41}$ erg s*-1* and its period was changing very rapidly from $p_{1}=1.43$ s to $p_{2}=1.13$ s during past 11 years. From equation (65) we can obtain the minimal magnetic field that can explain such a large luminosity, $\mu=6\times 10^{34}$ G cm3 which is rather unphysical. The most probable cause for this discrepancy is the unusually large luminosity of this object, probably affected by geometrical beaming. A beaming factor of several can reduce the luminosity and leads to results similar to that for M82 X-2. Still, for such objects with extreme mass accretion rates, the disc may become super-critical, and we need to consider effects like outflows and advection to explain the properties of the object.

Let us estimate period derivative $\dot{p}_{\rm s}$ expected for the three ULX-pulsars. Neglecting the unknown spin-down torque component in equation  (62) provides an upper limit for the period derivative. For the radiation-dominated disc, the magnetospheric radius is given by equation  (61), that gives us the period derivative:

[TABLE]

where $I=I_{45}10^{45}$g cm2 is the moment of inertia of a neutron star. For M82 X-2 we have $|\dot{p}_{\rm{spin-up}}|=(6.3-6.7)\times 10^{-9}$ s s*-1*. The observed period derivative (see observation 007 in Table 3 of Bachetti et al. 2014) is $|\dot{p}_{\rm obs}|=8\times 10^{-12}$ s s*-1* that is much smaller than the expected spin-up rate. There are episodes of spin-up and spin-down mentioned in the paper, so we can conclude that the object is close to the equilibrium (see also Ekşi et al. 2015, Lyutikov 2014 and Dall’Osso et al. 2015). The most interesting ULX-pulsar in NGC 5907 has the maximal observational spin-up rate of $|\dot{p}_{\rm obs}|=9.6\times 10^{-9}$ s s*-1* (Israel et al., 2017a) which is two times smaller than we obtain from equation (69) $|\dot{p}_{\rm s}|=2\times 10^{-8}$ s s*-1*. If we trust the luminosity estimate, this object should be in strong spin-up.

5.4 The role of advection

In our model we assume that the energy released in the disc is radiated locally from its surface. However, radiation transfer may be altered by radial advection if the mass accretion rate is high enough. Advection is important when the timescale for radial motion becomes comparable to or smaller than the timescale of vertical radiation diffusion. The characteristic timescale for radial drift is:

[TABLE]

where

[TABLE]

is the radial velocity. The diffusion time scale is:

[TABLE]

where the factor of 5 accounts for the fact that the energy release is distributed in height, and the effective optical depth is smaller by the same factor (see optical depth multiplier in equation 18). Substituting expression for $D$ from equation (13) we get the ratio of the two timescales:

[TABLE]

where $r_{\rm sp}\simeq\displaystyle\frac{3}{2}\dot{m}$ (Shakura & Sunyaev, 1973) is the spherization radius, inside which the standard disc model is expected to break down due to large disc thickness. As long as the local Eddington limit is not reached in the disc, or, equivalently, the radius of the magnetosphere is $r_{\rm in}\gtrsim r_{\rm sp}$ and $h_{\rm in}\lesssim r_{\rm in}$, accretion disc remains sub-critical, and the effects of advection, as well as other effects related to disc thickness, remain of minor importance. We will consider all these effects, together with the effects of outflows, in a separate paper.

6 Conclusions

We have developed a new model of the accretion disc around a strongly magnetized neutron star accreting at a high, possibly super-Eddington, rate. Calculating the radial structure of the disc allowed us to find more accurately the position of the disc-magnetosphere interface. Further outside, the disc very slightly deviates from the standard model. The correction factor to the classical Alfvenic radius is of the order unity, but becomes larger with increasing disc thickness reaching values of about two.

At high mass accretion rates and high magnetic moments, the pressure of the radiation field produced by the central source (accretion column) becomes stronger than the magnetic field pressure at the magnetospheric boundary meaning that the dynamics of the inner parts of the disc becomes strongly affected by the irradiation from the inner source. For magnetic moments $\mu\gtrsim 10^{30-32}$ G cm3 (position of this boundary strongly depends on accretion efficiency and viscosity parameter), this effect becomes important well below the accretion rate when the local Eddington limit is reached in the accretion disc.

The properties of ULX-pulsars can be generally explained in our model with magnetic moments of the order $10^{31}$ G cm3. The magnetar-scale fields proposed by Tsygankov et al. (2016) for M82 X-2 may be over-estimated by a factor of $1.5-5$ for $\eta=0.1-0.2$ because the large disc thickness and the irradiation by the central source increase the size of the magnetosphere. However, the corrections due to irradiation cannot be very high. To surpass the Eddington limit in the accretion column and allow accretion at $\dot{m}\gtrsim 500$, the local magnetic field should be of the order $\sim 10^{14}$ G (Mushtukov et al., 2015).

We also showed that for large enough disc thicknesses of $H/R\gtrsim 0.1$, when the inner disc becomes radiation-pressure-dominated, the size of the magnetosphere does not depend on the mass accretion rate. This fact should be reflected in the variability properties of the objects with high enough accretion rates ($\dot{m}\gtrsim 10$) and should also affect the spin-up rate dependence on the luminosity. However, the existing observational data do not allow yet to check these predictions.

Acknowledgements

This research was supported by the Center for International Mobility, the Väisälä Foundation, Nordita Visiting PhD Fellowship (AC), the Foundations’ Professor Pool, the Finnish Cultural Foundation, the National Science Foundation grant PHY-1125915 (JP) and the Academy of Finland grant 268740 (PA). AC and PA acknowledge support from RSF grant 14-12-00146 (vertical disc structure). The authors are grateful to Phil Armitage, Andrei Beloborodov, Wlodek Kluzniak, Dong Lai, Galina Lipunova, Nikolai Shakura, Sasha Tchekhovskoy and Sergey Tsygankov for useful discussions.

Appendix A Vertical disc structure

The vertical structure of a thin accretion disc depends on the opacity and pressure sources and is normally calculated by solving simultaneously the equations of vertical hydrostatic balance and radiation release and transfer. We cannot completely ignore the vertical structure because transition between vertically-integrated quantities like $\Sigma$ and equatorial-plane quantities like $\rho_{\rm c}$ requires some assumptions about the dependence of local quantities on $z$. To avoid calculating vertical structure for every point, we adopt a unified density dependence on relative vertical coordinate

[TABLE]

where $x=z/H$. The hydrostatic equilibrium equation

[TABLE]

allows to calculate the vertical pressure profile knowing that of the density as

[TABLE]

where $s\simeq 1$ is some coefficient. According to Shakura et al. (1978), $s=1.17$ for the radiation-dominated disc with convection. We will hereafter use $s=1$ unless otherwise stated. For $s=1$, central pressure and density are related as

[TABLE]

Thus, vertically integrated pressure and surface density are:

[TABLE]

[TABLE]

Now let us show the derivation of relation (16) between the surface temperature and the central temperature in the disc. The local energy dissipation rate:

[TABLE]

Integrating this expression along vertical direction from 0 to $z$ and using equation (76) we get the energy flux:

[TABLE]

Then we substitute this flux into the vertical diffusion equation $F=-D\nabla_{z}\epsilon=-\displaystyle\frac{c}{3\kappa\rho}{\rm d}(aT^{4})/{\rm d}z$, taking $\rho$ from equation (74) we get:

[TABLE]

where $T_{\rm s}$ and $T_{\rm c}$ are the surface and central temperatures, respectively. Assuming that the surface temperature equals the effective temperature

[TABLE]

we get an expression for the central disc temperature:

[TABLE]

In the paper we were mainly interested in the objects with high luminosities. The inner parts of their discs are in radiation-pressure-dominated regime, that justifies the usage of $s=1$. However, if we consider gas-pressure-dominated disc, $s$ can be larger than 1, about $2-3$ (Ketsaris & Shakura, 1998). For $s=2$ and $3$ the thickness of the disc is expressed as

[TABLE]

and

[TABLE]

respectively. Equation (84) derived for different values of $s$ from 1 to 3 has a very stable form, with the expression in the brackets always equal to $1+\left(0.21\pm 0.01\right)\kappa\Sigma$. Energy release in the disc is very much the same for any vertical structure. The resulting $\xi$ for different vertical structures differs by no more than $10$ per cent, mainly as a result of varying disc thickness multiplier (compare equations 9, 85 and 86), that changes by less than $25$ per cent.