Broadband aperiodic variability in X-ray pulsars: accretion rate fluctuations propagating under the influence of viscous diffusion

TL;DR

This paper models aperiodic X-ray variability in accreting neutron stars by analyzing viscous diffusion in truncated accretion discs, linking observed spectral features to disc dynamics and magnetic field influences.

Contribution

It introduces the most general analytical solutions for viscous diffusion in truncated accretion discs, enhancing understanding of variability origins in X-ray pulsars.

Findings

Break in power density spectra linked to dynamo process timescale.

Initial variability timescale is a few times longer than local Keplerian time.

Model fits well with observational data from A 0535+26.

Abstract

We investigate aperiodic X-ray flux variability in accreting highly magnetized neutron stars - X-ray pulsars (XRPs). The X-ray variability is largely determined by mass accretion rate fluctuations at the NS surface, which replicate accretion rate fluctuations at the inner radius of the accretion disc. The variability at the inner radius is due to fluctuations arising all over the disc and propagating inwards under the influence of viscous diffusion. The inner radius varies with mean mass accretion rate and can be estimated from the known magnetic field strength and accretion luminosity of XRPs. Observations of transient XRPs covering several orders of magnitude in luminosity give a unique opportunity to study effects arising due to the changes of the inner disc radius. We investigate the process of viscous diffusion in XRP accretion discs and construct new analytical solutions of the diffusion equation applicable for thin accretion discs truncated both from inside and outside. Our solutions are the most general ones derived in the approximation of Newtonian mechanics. We argue that the break observed at high frequencies in the power density spectra of XRPs corresponds to the minimal time scale of the dynamo process, which is responsible for the initial fluctuations. Comparing data from the bright X-ray transient A 0535+26 with our model, we conclude that the time scale of initial variability in the accretion disc is a few times longer than the local Keplerian time scale.

Equations

$$R_{\mathrm{m}}=2.4\times 10^8\Lambda B_{12}^{4/7}{L}_{37}^{-2/7}{m}^{1/7}{R}_{\mathrm{NS},6}^{10/7}\mathrm{cm},$$

$$t_{\mathrm{K}}(R_{\mathrm{m}})=\frac{2\pi }{{\Omega }_{\mathrm{K}}}\simeq 1.9{\Lambda }^{3/2}{B}_{12}^{6/7}{L}_{37}^{-3/7}{M}_{1.4}^{-2/7}{R}_{\mathrm{NS},6}^{15/7}\mathrm{s}.$$

$$v_{\mathrm{r}}\simeq \alpha v_{\phi }{\left(\frac{H}{R}\right)}^2\ll v_{\phi },$$

$$t_{\mathrm{v}}(R)=\frac{t_{\mathrm{K}}(R)}{3\pi \alpha }{\left(\frac{H(R)}{R}\right)}^{-2}\gg t_{\mathrm{K}}(R).$$

$$\frac{H(R)}{R}\approx 0.03{\alpha }^{-1/10}{L}_{37}^{3/20}{m}^{-21/40}{R}_{\mathrm{NS},6}^{3/20}{R}_8^{1/8},$$

$$t_{\mathrm{d}}\approx k_{\mathrm{d}}{t}_{\mathrm{K}},$$

$$f_{\mathrm{K}}=t_{\mathrm{K}}^{-1}\simeq 0.53{\Lambda }^{-3/2}{B}_{12}^{-6/7}{L}_{37}^{3/7}{M}_{1.4}^{2/7}{R}_{\mathrm{NS},6}^{-15/7}\mathrm{Hz}$$

$$f_{\mathrm{v}}=t_{\mathrm{v}}^{-1}\simeq 0.94{\alpha }_{0.1}{f}_{\mathrm{K}}(R){\left(\frac{H(R)}{R}\right)}^2\ll f_{\mathrm{K}}(R),$$

$$f_{\mathrm{K}}\gtrsim f_{\mathrm{d}}\gg f_{\mathrm{v}}.$$

$$\frac{\partial \Sigma (R,t)}{\partial t}=\frac{1}{R}\frac{\partial }{\partial R}\left[R^{1/2}\frac{\partial }{\partial R}\left(3\nu \Sigma R^{1/2}\right)\right],$$

$$\Sigma (R,t)=\underset{R_{\mathrm{in}}}{\overset{R_{\mathrm{out}}}{\int }}G(R,R^{\prime },t-t_0)\Sigma (R^{\prime },t_0)\mathrm{d}{R}^{\prime },$$

$$\nu (R)={\nu }_0{\left(\frac{R}{R_0}\right)}^n.$$

$$\nu =\alpha {\Omega }_{\mathrm{K}}{R}^2{\left(\frac{H}{R}\right)}^2\propto R^{3/4},$$

$$\dot{M}(R,t)=6\pi R^{1/2}\frac{\partial }{\partial R}\left(\nu \Sigma (R,t)R^{1/2}\right).$$

$$\dot{M}(R,t)=\underset{R_{\mathrm{in}}}{\overset{R_{\mathrm{out}}}{\int }}{G}_{\dot{M}}(R,R^{\prime },t){\otimes }_t\frac{\partial {\Sigma }^{\ast }(R^{\prime },t)}{\partial t}\mathrm{d}{R}^{\prime },$$

$$G_{\dot{M}}(R,R^{\prime },t)=6\pi R^{1/2}\frac{\partial }{\partial R}\left(\nu G(R,R^{\prime },t)R^{1/2}\right).$$

$${\overline{G}}_{\dot{M}}(R,R^{\prime },f)=\underset{-\infty }{\overset{\infty }{\int }}\mathrm{d}x{G}_{\dot{M}}(R,R^{\prime },x)e^{-2\pi ifx}$$

$$S_{\dot{M}}(R,f)\simeq \underset{R_{\mathrm{in}}}{\overset{R_{\mathrm{out}}}{\int }}\frac{\mathrm{d}{R}^{\prime }}{(R^{\prime }{)}^2}\Delta R(R^{\prime })|{\overline{G}}_{\dot{M}}(R,R^{\prime },f){|}^2{S}_a(R^{\prime },f),$$

$$G(R,R^{\prime },t)=(2-n)R^{-n-1/4}{R}^{\prime 5/4}{R}_{\mathrm{out}}^{n-2}$$

$$\times \sum _i\exp \left[-2{\left(1-\frac{n}{2}\right)}^2{k}_i^2\frac{t}{t_{\mathrm{v}}}\right]\frac{V_l(k_i{x}^{\prime },k_i{x}_{\mathrm{in}})V_l(k_{i}x,k_i{x}_{\mathrm{in}})}{V_l^2(k_i{x}_{\mathrm{out}},k_i{x}_{\mathrm{in}})},$$

$$V_l(u,\nu )=J_l(u)J_{-l}(\nu )-J_{-l}(u)J_l(\nu ),$$

$$k_i[J_{l-1}(k_i{x}_{\mathrm{out}})J_{-l}(k_i{x}_{\mathrm{in}})-J_{-l-1}(k_i{x}_{\mathrm{out}})J_l(k_i{x}_{\mathrm{in}})]$$

$$-\frac{2l}{x_{\mathrm{out}}}{J}_{-l}(k_i{x}_{\mathrm{out}})J_l(k_i{x}_{\mathrm{in}})=0.$$

$${\overline{G}}_{\dot{M}}(R,R^{\prime },f)=6\pi R^{1/2}{\left(\frac{R_{\mathrm{out}}}{R_{\mathrm{in}}}\right)}^{2-n}$$

$$\times \frac{\partial }{\partial R}\{R^{n+0.75}\sum _i\frac{V_l(k_i{x}^{\prime },k_i{x}_{\mathrm{in}})V_l(k_{i}x,k_i{x}_{\mathrm{in}})}{V_l^2(k_i{x}_{\mathrm{out}},k_i{x}_{\mathrm{in}})}$$

$$\times \frac{1}{4\pi if/f_{\mathrm{in}}+k_i^2(n-2{)}^2}\}.$$

$${\overline{G}}_{\dot{M}}(R,R^{\prime },f)=\frac{6\pi }{4}(2-n)R^{-1/4}{\left(\frac{R_{\mathrm{out}}}{R_{\mathrm{in}}}\right)}^{2-n}$$

$$\times \sum _i\frac{V_l(k_i{x}^{\prime },k_i{x}_{\mathrm{in}})W_l(k_{i}x,k_i{x}_{\mathrm{in}})}{V_l^2(k_i{x}_{\mathrm{out}},k_i{x}_{\mathrm{in}})}\frac{1}{4\pi if/f_{\mathrm{in}}+k_i^2(n-2{)}^2},$$

$$W_l(k_{i}x,k_i{x}_{\mathrm{in}})$$

$$F(R,t)=3\pi h{\nu }_0\Sigma (R,t){\left(\frac{R}{R_0}\right)}^n,$$

$$\dot{M}(R,t)=\frac{\partial F(R,t)}{\partial h}.$$

$$F(R)=F_{\mathrm{in}}+\underset{h_{\mathrm{in}}}{\overset{h}{\int }}\mathrm{d}h\dot{M}(R),$$

$$F(R,t)=F_{\mathrm{in}}+\underset{h_{\mathrm{in}}}{\overset{h}{\int }}\mathrm{d}h{\dot{M}}_0(R)+F_2(R,t),$$

Full text

††pubyear: 2019††pagerange: Broadband aperiodic variability in X-ray pulsars: accretion rate fluctuations propagating under the influence of viscous diffusion–B.2

1 Introduction

X-ray pulsars (XRPs) are highly magnetized (typical magnetic field strength at the surface $\gtrsim 10^{12}\,{\rm G}$) neutron stars (NSs) in close binary systems, whose luminosity in X-rays is caused by accretion from the companion star (Walter et al., 2015). The accretion flow, in the form of a stellar wind or accretion disc, is interrupted by the strong NS magnetic field ($B$-field) at the magnetospheric radius, $R_{\rm m}$. This radius depends on the accretion flow geometry, mass accretion rate $\dot{M}$ and $B$-field strength and structure (Pringle & Rees, 1972; Davidson & Ostriker, 1973; Lipunov, 1978; Aly, 1980). If the mass accretion rate is sufficiently high for matter to go through the centrifugal barrier caused by rotation of the strongly magnetized NS, the accretion flow penetrates into the magnetosphere and follows magnetic field lines to reach the NS surface in regions located close to the magnetic poles. There the accretion flow loses its kinetic energy, which is radiated mostly in the X-ray energy band. Misalignment between the rotational and magnetic axes results in the phenomenon of XRPs. If the matter cannot penetrate through the centrifugal barrier, the accretion process stops, leading to the so-called ”propeller” effect (Illarionov & Sunyaev, 1975; Syunyaev & Shakura, 1977; Wang & Robertson, 1985; Romanova et al., 2004; Tsygankov et al., 2016a; Tsygankov et al., 2016b).

The accretion luminosity of known XRPs covers several orders of magnitude from $\sim 10^{33}\,{\rm erg\ \rm s^{-1}}$ up to $\sim 10^{40}\,{\rm erg\ \rm s^{-1}}$. Many sources are transients and show significant variability of luminosity. The brightest XRPs belong to the recently discovered class of pulsating ultra-luminous X-ray sources (ULXs, see e.g. Bachetti et al. 2014; Israel et al. 2017). The observed X-ray luminosity in pulsating ULXs exceeds the Eddington value, which is $L_{\rm Edd}\approx 2\times 10^{38}\,{\rm erg\ \rm s^{-1}}$ for a NS, by a factor of hundreds (Israel et al., 2017).

The X-ray flux originates from regions located close to the NS magnetic poles. The geometry of the radiating regions can be affected by radiation pressure and is determined by the mass accretion rate. At relatively low mass accretion rates ($\lesssim 10^{17}\,{\rm g\,s^{-1}}$), the accretion flow is stopped by Coulomb collisions and/or plasma oscillations at surface layers of the NS (Zel’dovich & Shakura, 1969), while at high mass accretion rates ($\gtrsim 10^{17}\,{\rm g\,s^{-1}}$), the radiation pressure becomes strong enough to stop the flow at some height above the stellar surface (Basko & Sunyaev, 1976; Mushtukov et al., 2015b). In the latter case, an accretion column forms. The column can provide luminosity well above the Eddington value because it is confined by the strong magnetic field, and thus radiation pressure can be largely reduced due to the reduction of scattering cross-sections in a strong B-field (Wang & Frank, 1981; Mushtukov et al., 2015a), plus photon bubbles can come into play reducing effective radiation pressure (Arons, 1992). The radiation of the accretion column is likely beamed towards the NS surface (Kaminker et al., 1976; Lyubarskii & Syunyaev, 1988) and, thus, a large fraction of X-ray flux is reprocessed/reflected by the atmosphere (Poutanen et al., 2013).

XRPs show strong aperiodic variability of X-ray flux over a very broad frequency range similar (modulo mass scaling) to what is detected in accreting black holes (BHs, see e.g. Revnivtsev et al. 2000) and active galactic nuclei (AGN, see e.g. McHardy et al. 2004). The time scale of the observed variability extends down to milliseconds. Comparing accreting NSs with BHs we note that the power density spectra of weakly magnetized NSs contains much stronger variability at frequencies close of one kHz (Sunyaev & Revnivtsev, 2000). The power density spectrum (PDS) typically includes a broad component, which can be approximated by a broken (or double broken) power-law (Hoshino & Takeshima, 1993), and narrow features that are classified as quasi-periodic oscillations (QPOs, see Takeshima et al. 1994). Both components are detected to evolve with long term trends in the observed luminosity (Revnivtsev et al. 2009).

The aperiodic variability in X-ray binaries and AGNs is naturally explained by the propagating fluctuations model (Lyubarskii, 1997; Churazov et al., 2001; Titarchuk et al., 2007). According to the model, the initial fluctuations of the mass accretion rate arise at different radial coordinates in the accretion disc and then propagate inwards and outwards, modulating the fluctuations arising at other radial coordinates. In this scenario, different time scales are injected into the accretion flow at different distances from the central object, while the observed variability of X-ray flux reflects the variability of the mass accretion rate at the inner parts of accretion disc (Kotov et al., 2001; Ingram & van der Klis, 2013; Ingram, 2016; Mushtukov et al., 2018).

There is a fundamental difference between the X-ray flux variability in accreting BHs and XRPs: in the case of BH binaries, an observer detects X-ray photons originating from the accretion disc itself, while in the case of XRPs, photons (seed photons at least, see below) mostly originate from the NS surface. We can imagine a few sources of aperiodic variability in XRPs: (i) the variability of X-ray flux due to the mass accretion rate fluctuations at the inner disc radius, which determines directly the variability of the mass accretion rate at the NS surface; (ii) the variability caused by changes of the beam pattern due to changes in the geometry of the radiating region (Mushtukov et al., 2018); (iii) the variability due to reprocessing of seed X-ray photons by the accretion flow (in the form of a hot optically thin accretion flow at low mass accretion rates, see e.g. Narayan & Yi 1995, or an optically thick accretion flow at high mass accretion rates, see e.g. Mushtukov et al. 2017); (iv) production of photon bubbles in the accretion column at high mass accretion rates (Arons, 1992; Begelman, 2006). The time scales of photon emission and diffusion in radiating regions also affect the observed variability properties, but only on time scales much shorter than those of the processes discussed above. In this paper, we focus on the aperiodic variability due to propagating fluctuations of the mass accretion rate in a disc, and assume that the fluctuations in X-rays replicate fluctuations of the mass accretion rate at the inner disc radius.

The changes of PDS with luminosity of the XRP are largely caused by changes of the inner disc radius, which depends on the mass accretion rate (see Fig.1, see e.g. Revnivtsev et al. 2009). In a number of cases, the magnetic field strength of the XRP is known (it can be measured from cyclotron line scattering features in the X-ray spectrum), and thus the disc inner radius $R_{\rm m}$, and the Keplerian frequency there, can be estimated. The similar behavior of the observed break frequency and estimated Keplerian frequency at the magnetospheric radius reveals that the changes in PDS break frequency with X-ray luminosity are largely caused by changes in $R_{\rm m}$ (Revnivtsev et al., 2009). The power-law dependence of break frequency on mass accretion rate is consistent with the break frequency being proportional to the Keplerian frequency at $R_{\rm m}$. Constraint of the constant of proportionality obtained through detailed physical modelling offers the opportunity to use the measured PDS break frequency to estimate $R_{\rm m}$ for XRPs with poorly constrained magnetic field strength. Moreover, such constraints will feed into our understanding of accretion disc variability in general, providing diagnostics of the accretion flow geometry of XRBs and AGN through timing properties. Similarly, the low frequency breaks detected in some XRPs contain information about the outer disc radius (Gilfanov & Arefiev, 2005).

In this paper, we focus on XRPs where accretion takes place through a geometrically thin accretion disc (this assumption puts limitations on the range of mass accretion rate under consideration). We focus on variability caused by processes in the accretion disc: mass accretion rate fluctuations arising and propagating in the disc due to the process of viscous diffusion. We analyze the effects arising from truncation of the accretion disc at certain inner $R_{\rm in}$ and outer $R_{\rm out}$ radii, assuming that the disc loses its mass through the inner radius only. We consider the general case of non-zero torque at the inner disc radius and provide a theoretical background for calculations of the PDS in accreting strongly magnetized NSs.

The paper consists of five sections. In Section 2 we discuss the basic features of accretion discs in XRPs at different luminosity states, accretion disc geometry and typical time scales in the accretion flow. The analytical theory of propagating fluctuations of the mass accretion rate is discussed in Section 3, where we introduce a new analytical solution of the viscous diffusion equation accounting for disruption of the accretion disc both from the inside and outside (see Section 3.2). Section 4 presents our numerical results based on the theory and novel Green functions developed in Section 3. Summary and conclusions are given in Section 5.

2 Accretion discs in X-ray pulsars

In this paper, we will define and apply a model for the aperiodic X-ray variability of X-ray pulsars. In this Section we first explore the expected disc geometry (Section 2.1) and accretion regime (Section 2.2) before considering typical timescales of the system (Section 2.3), and finally commenting on specific features of the accretion flow in some classes of XRPs that we do not consider here, but must be included in advanced models of aperiodic variability (Section 2.4)

2.1 Accretion disc geometry

The accretion disc in XRPs is truncated at the magnetospheric radius, $R_{\rm m}$, which can be roughly estimated by comparing the $B$-field pressure with the ram pressure of the accreting material (Lipunov, 1987; Frank et al., 2002):

[TABLE]

where $\Lambda<1$ is a constant which depends on the accretion flow geometry, with $\Lambda=0.5$ being a commonly used value for the case of accretion through the disc (Ghosh & Lamb, 1979; Lai, 2014), $B_{12}$ is the magnetic field strength $B$ at the NS surface in units of $10^{12}\,{\rm G}$, $L_{37}$ is the accretion luminosity $L$ in units of $10^{37}\,{\rm erg\ \rm s^{-1}}$, $m$ is the NS mass in units of solar masses $M_{\odot}$, and $R_{\rm NS,6}$ is the NS radius in units of $10^{6}\,{\rm cm}$. Apart from the fact that $\Lambda$ depends on a specific model of the $B$-field and disc interaction, the inner radius of the accretion disc is also affected by the magnetic dipole inclination (see e.g. Lipunov 1978; Scharlemann 1978; Aly 1980). For the typical $B$-field strengths and mass accretion rates in XRPs, the inner disc radius is so large that the disc cannot produce a noticeable fraction of luminosity in the X-ray energy band. The boundary layer of the accretion disc, where material penetrates into the magnetosphere, is formed due to instabilities developing at the inner disc radius (magnetic Kelvin-Helmholtz and Rayleigh-Taylor, and reconnection, see e.g. Spruit & Taam 1990; Lai 2014). The accretion flow, penetrating into the magnetosphere, follows magnetic field lines and reaches the NS magnetic poles to provide most of the X-ray flux at the surface, where the kinetic energy is converted into heat.

The outer radius of the accretion disc is determined by tidal torques in the binary system, which removes the angular momentum from the matter near the disc edge and prevents disc spreading. The tidal radius $R_{\rm tid}$ is likely close to the size of the Roche lobe: $R_{\rm tid}\approx(0.8-0.9)R_{\rm L}$ (Paczynski, 1977; Papaloizou & Pringle, 1977).

2.2 Accretion regime

An optically thick accretion disc can be divided into three zones according to the dominating pressure and opacity sources (see e.g. Shakura & Sunyaev 1973; Suleimanov et al. 2007). Gas pressure and Kramers opacity dominate in the outer regions (C-zone), gas pressure and electron scattering dominate in the intermediate zone (B-zone) and radiation pressure dominates in the inner zone (A-zone). Fig. 2 illustrates the boundaries between these three regimes. As marked by the dot-dashed line, the effective temperature in the outer parts of the disc may drop below 6500K depending on mass accretion rate. In this case, the hydrogen recombines and a ‘cooling’ wave propagates inward (e.g. Lasota 2001) and can reach $R_{\rm in}$ making all the hydrogen in the disc neutral. Stable accretion from such a ”cold disc” state has recently been discovered in a few XRPs (Tsygankov et al., 2017a; Tsygankov et al., 2017b). The grey region in Fig. 2 denotes an advection dominated accretion flow (ADAF) regime (Narayan & Yi, 1995), which is beyond the scope of this paper.

The dashed lines in Fig. 2 show the magnetospheric radius for a selection of typical values of XRP magnetic field strength. We see that, for all but the highest accretion rates, the inner disc is expected to be in the C-zone regime. Our model for the aperiodic variability XRPs will therefore assume that the disc is entirely in the C-zone. We will apply our model to A 0535+26, whose disc extent is also plotted in Fig. 2 for two accretion states. The tidal radius of this source is not well known, but that the disc does not appear to be in the cold disc regime implies that $R\lesssim 6\times 10^{9}$ cm. Note that known orbital period and an estimated mass of a companion star in A 0535+26 ($P_{\rm orb}\simeq 110$ days and $M_{\rm comp}\sim 14\,M_{\sun}$, see Giangrande et al. 1980) imply the tidal radius $R_{\rm tid}\gg 6\times 10^{9}\,{\rm cm}$, but the outer regions of accretion disc are expected to be in a cold state and only slightly affect the process of viscous diffusion in internal hot part of accretion flow.

2.3 Typical time scales

2.3.1 Time scales in an accretion disc

There are a few time scales which arise naturally from the physics of disc accretion onto a magnetized NS: (i) the dynamical time scale $t_{\rm K}$ given by Keplerian rotation, (ii) the viscous time scale $t_{\rm v}$, (iii) the time scale of the dynamo process, which provides the mechanism of viscosity in the accretion disc, and (iv) the time scale of instabilities $t_{\rm ins}$ at the inner disc radius, which are responsible for matter penetration into the NS magnetosphere.

The velocity of Keplerian rotation in accretion disc at radius $r$ is $v_{\varphi}=(GM/r)^{1/2}$ which corresponds to angular velocity $\Omega_{\rm K}=(GM/r^{3})^{1/2}$. Thus, the dynamical timescale at the disc inner radius is given by

[TABLE]

The viscous time scale $t_{\rm v}$ is determined by the radial velocity in the accretion disc $v_{\rm r}$, which is a result of viscous diffusion of accreting matter in the disc. The diffusion process is described by the viscous diffusion equation and depends on the kinematic viscosity $\nu$ (see Section 3). In the simplified case of an $\alpha$-disc (Shakura, 1972; Shakura & Sunyaev, 1973), the radial velocity is given by

[TABLE]

where $0<\alpha<1$ is the dimensionless viscosity parameter, $H$ is the disc scale height and $c_{\rm s}$ is the isothermal sound speed. The viscous time scale can be roughly estimated as

[TABLE]

In the C-zone of an accretion disc, the relative disc scale height can be estimated as (Suleimanov et al., 2007)

[TABLE]

where $R_{8}$ is the radial coordinate $R$ in accretion disc in units of $10^{8}\,{\rm cm}$.

Another important time scale of the problem is the typical timescale of initial fluctuations of the mass accretion rate in the disc $t_{\rm d}$. The fluctuations are likely caused by a magnetic dynamo that generates a poloidal field component in a random fashion and serves as the source of viscosity in the accretion flow, in the form of correlated fluctuations in magnetic stress (Balbus & Hawley, 1991; Hawley et al., 1995; Brandenburg et al., 1995; Liska et al., 2018). The dynamo timescale has been shown to be close to local Keplerian time-scale (Tout & Pringle, 1992; Stone et al., 1996): $t_{\rm d}\gtrsim t_{\rm K}$. Here, we assume

[TABLE]

where the exact value of $k_{\rm d}>1$ is not known a priori, requiring detailed numerical simulations or/and detailed observational diagnostics.

These typical time scales determine the typical frequencies in the accretion disc. The Keplerian frequency:

[TABLE]

the viscous frequency:

[TABLE]

where $\alpha_{0.1}\equiv\alpha/0.1$, and the frequency of dynamo processes $f_{\rm d}=t_{\rm d}^{-1}$. The introduced frequencies are related as

[TABLE]

2.3.2 The time scale of emitting processes in the vicinity of the neutron star surface

The geometry of the emitting region in the vicinity of the NS surface depends on the mass accretion rate and $B$-field strength (Mushtukov et al., 2015b). At low mass accretion rates (sub-critical, $\lesssim 10^{17}\,{\rm g\,s^{-1}}$), the accretion flow is braked by Coulomb collisions in the NS atmosphere (Zel’dovich & Shakura, 1969). If the mass accretion rate is high enough (super-critical, $\gtrsim 10^{17}\,{\rm g\,s^{-1}}$), the accretion flow is stopped above NS surface by a radiation dominated shock (Basko & Sunyaev, 1976; Mushtukov et al., 2015a). Below the shock region, the accretion flow slowly settles to the stellar surface emitting energy in X-rays. Radiation from the accretion column can be beamed towards the NS surface (Lyubarskii & Syunyaev, 1988; Poutanen et al., 2013). In this case, the photon energy flux is partly reprocessed/reflected by the NS atmosphere. The time scales of photon diffusion from a hot spot/accretion column and reprocessing of X-ray flux by the atmosphere are much shorter than the Keplerian time scale at the inner disc radius. Thus, these processes do not influence the PDS in the frequency range of our interest.

At high mass accretion rates and relatively low $B$-field strength ($B<10^{12}\,{\rm G}$) at the NS surface, one would expect instabilities in the accretion column to result in in X-ray flares (Basko & Sunyaev, 1976). The time scale of the flares is expected to be in the range $10^{-4}<t_{\rm flares}<10^{-1}\,{\rm s}$.

2.4 Accretion disc features in some classes of X-ray pulsars

2.4.1 Accretion discs in transients

Outbursts of X-ray transients are likely driven by the development of a disc instability triggered by gradual accumulation of matter in the disc or episodic material capture from a companion star (like it happens in Be/X-ray transients, see e.g. Reig 2011) and caused by the strong dependence of viscosity on temperature for temperatures in the range $T\sim 6500\,{\rm K}$ (see e.g. Lasota 2001). The development of this instability results in cooling and heating waves propagating in the accretion disc. The physical conditions (particularly, the viscosity) in the disc are significantly different at opposite sides of the cooling/heating waves. It is likely that only the ”hot” part of the accretion disc effectively contributes to the propagating fluctuations in mass accretion rate (the initial fluctuations are smaller in weakly ionized parts of the disc and suppression of variability at high frequencies is stronger in a cold part of the accretion disc because of the much longer viscous time scale). We expect that the evolution of the outburst and dynamics of cooling/heating waves in accretion disc can affect the PDS of aperiodic variability and even cause QPOs at low frequencies. However, these effects are beyond the scope of the paper and will be considered separately.

2.4.2 Accretion discs in ULX pulsars

It has been recently discovered that the luminosity of XRPs can exceed hundreds of Eddington luminosities (Bachetti et al., 2014; Israel et al., 2017) and at least a fraction of ULXs host accreting NSs. The exact mechanism of such an extreme mass accretion rate is still under debate (Basko & Sunyaev, 1976; Paczynski, 1992; Mushtukov et al., 2015a, 2017; King et al., 2017), but it is clear that the conditions of accretion discs in these systems are different from those in normal XRPs. In particular, one would expect a geometrically thick inner part of the accretion disc, where the pressure is dominated by radiation (see A-zone in Fig. 2). The thick inner part of the disc can affect the timing properties of aperiodic variability due to different radial dependence of the kinematic viscosity $\nu$ and, therefore, different transfer properties of propagating fluctuations. A radiation pressure dominated region and possible mass losses from the accretion disc (Lipunova, 1999; Poutanen et al., 2007) can also result in a dependence of the inner disc radius on the mass accretion rate different from the one given by equation 1 (see Tab. 1 in Psaltis & Chakrabarty 1999, see also Chashkina et al. 2017, 2019; Mushtukov et al. 2019).

Another feature of magnetized NSs at extremely high mass accretion rates is an optically thick envelope formed by the accretion flow moving from magnetospheric radius to the central object (Mushtukov et al., 2017). The envelope hides the NS from a distant observer and affects the observational manifestation of ULX pulsars including their fast aperiodic variability. Particularly, the envelope suppresses high-frequency variability of the X-ray energy flux (Mushtukov et al., 2019).

3 Propagating fluctuations of the mass accretion rate

In this Section, we develop a model for propagating accretion rate fluctuations in an XRP disc with inner radius at the magnetospheric radius, $R_{\rm in}=R_{m}$, and outer radius at the tidal radius, $R_{\rm out}=R_{\rm tid}$. We consider small amplitude local fluctuations of mass accretion rate/surface density arising at each radius in the accretion disc and propagating due to the process of viscous diffusion (Mushtukov et al., 2018). The fluctuations propagate both inwards and outwards (Mushtukov et al. 2018) contributing to the total variability of the mass accretion rate at each radial coordinate. The propagation of fluctuations is governed by the viscous diffusion equation (Section 3.1). In Section 3.2, we present our solution to the diffusion equation that, for the first time, accounts for a finite inner and outer radius which we will use to model the PDS of XRPs. In Section 3.3, we explore the properties of our new Green function and in Section 3.4, we define our model for the power spectrum of initial fluctuations assumed to be injected into the disc.

3.1 The viscous diffusion equation

Propagation of the mass accretion rate fluctuations can be accurately described by the solutions of the equation of viscous diffusion (Lynden-Bell & Pringle, 1974):

[TABLE]

where $\Sigma$ is the local surface mass density, $R$ is the radial coordinate, $\nu$ is the kinematic viscosity, and $t$ is time. In the particular case whereby kinematic viscosity is not dependent on the local surface density, the equation of viscous diffusion is linear. Then the solutions of the equation can then be found using the Green functions:

[TABLE]

where $G(R,R^{\prime},t)$ is a Green function describing evolution of the surface density, $\Sigma(R,t_{0})$ is the surface density distribution over the radial coordinate $R$ at $t=t_{0}$, and $R_{\rm in}$ and $R_{\rm out}$ are inner and outer disc radius respectively.

The particular Green functions are determined by viscosity dependence on the radial coordinate $\nu(R)$ and boundary conditions at the inner $R_{\rm in}$ and outer $R_{\rm out}$ radii of the disc. The solution to the equation of viscous diffusion simplifies if we assume a power-law dependence of kinematic viscosity on radial coordinate:

[TABLE]

In zone-C of the Shakura & Sunyaev (1973) disc model, which is relevant to the XRP discs we consider here, we have:

[TABLE]

i.e. $n=3/4$ (assuming constant $\alpha$). The assumption (10) is applicable for a wide range of mass accretion rates and radii if we consider relatively small oscillations of accretion rate around an average value. Under assumption (10), Green functions have been derived analytically for a few particular cases: (1) for the case of $R_{\rm in}=0$, $R_{\rm out}=\infty$ the Green functions was derived by Lynden-Bell & Pringle (1974); (2) for the case of $R_{\rm in}>0$, $R_{\rm out}=\infty$ the Green function was derived by Tanaka (2011); (3) for the case of $R_{\rm in}=0$, $R_{\rm out}<\infty$ and zero mass accretion rate at $R_{\rm out}$ the Green function was derived by Lipunova (2015). The zero torque/zero mass accretion rate condition at $R_{\rm in}$ was adopted in each of the mentioned works. Green functions accounting for the effects of general relativity in the Kerr geometry have been derived by Balbus (2017). In this paper, we present a Green function solution to equation (8) that accounts for $R_{\rm in}>0$ and $R_{\rm out}<\infty$, and can therefore be used to describe XRP discs.

Variability of the surface density is accompanied by variability of local mass accretion rate:

[TABLE]

Thus, the mass accretion rate variability can be represented as

[TABLE]

where $\otimes_{x}$ denotes the convolution in $x$-variable, $\Sigma^{*}$ is the initial fluctuation of the surface density due to local MHD processes, and $G_{\dot{M}}(R,R^{\prime},t)$ is the Green function for the mass accretion rate:

[TABLE]

Equation (14) gives the Green function in the time domain, but in the analysis of timing properties of mass accretion rate fluctuations it is convenient to use the Green function in the frequency domain.

The Green function in the frequency domain is given by the Fourier transform:

[TABLE]

and has the physical meaning of a transfer functions that describes how the variability at radial coordinate $R^{\prime}$ affects variability at the the radial coordinate $R$ (see e.g. Mushtukov et al. 2018). The Green function in the frequency domain belongs to the complex plane and contains information both about the amplitude and the time delay of the transferred fluctuations. The PDS of mass accretion rate fluctuations at radius $R$ can be calculated as

[TABLE]

where $S_{a}(R^{\prime},f)$ is the PDS of initial fluctuations, $\Delta R(R)$ is a typical radial scale on which the initial fluctuations can be considered coherent, and the integral is taken over the whole range of radial coordinates in the accretion disc. It is likely that $\Delta R(R)\sim H(R)$ (Hogg & Reynolds, 2016) and we use the estimation given by (5) in our numerical calculations. Thus, the PDS of the mass accretion rate at any radius (including the inner disc radius) is entirely determined by the PDS of initial fluctuations (see Section 3.4) and the transfer functions given by suitable Green functions in the frequency domain (see Sections A and 3.2). Note that the expressions introduced in this section do not account for non-linear effects arising because of interaction of propagating fluctuations with each other. This approximation is good in the case of a small fractional rms (Mushtukov et al., 2018).

Propagation of fluctuations of the mass accretion rate under the influence of viscous diffusion suppresses high-frequency ($f>f_{\rm v}$) variability (Kotov et al., 2001; Mushtukov et al., 2018). The initial variability at a given radius in the accretion disc affects mass accretion rate variability both in the inner and outer part of the disc (Mushtukov et al., 2018).

3.2 New Green function for an accretion disc with finite inner and outer radii

The Tanaka (2011) Green function solution to equation (8) is reproduced in Appendix A. This solution accounts for a finite disc inner radius but assumes $R_{\rm out}=\infty$ and is therefore not appropriate for modelling XRP accretion discs, for which the disc’s inner radius is not much smaller than its outer radius. We therefore derive a new Green function solution to equation (8) with finite $R_{\rm in}$ and $R_{\rm out}$ (Appendix B), which we will use for our subsequent analysis of XRP PDSs. We assume zero mass accretion rate at the outer radius, $\dot{M}(R_{\rm out})=0$, implying that the disc loses mass from the inner radius only. This is appropriate for $R_{\rm out}$ set by the tidal radius or by the boundary between the hot and cold parts of the disc (Lipunova & Malanchev, 2017). Note, that mass inflow from the companion star is possible at any radial coordinate within the outer disc radius. The exact geometry of the mass inflow is determined by the geometry of the binary system and the mechanism of mass transfer (wind accretion, accretion from the Lagrangian point or matter capture during the periastron passage). We will first assume zero torque at the inner radius (Section 3.2.1), before considering general torque at the inner radius (Section 3.2.2). Remarkably, we find that our Green function derived assuming zero torque at $R_{\rm m}$ can be generally applied to XRP discs.

3.2.1 The case of a zero-torque inner boundary

In this limit, our Green function is (see Appendix B and Fig. 3a):

[TABLE]

where $x=(R/R_{\rm out})^{1-n/2}$,

[TABLE]

$J_{l}(x)$ are the Bessel functions of the first order, and $k_{i}$ are roots of the transcendental equation

[TABLE]

Using (14) and (15) we get the Green function for the mass accretion rate in the frequency domain (see Fig.5 and 6):

[TABLE]

The expression (3.2.1) can be reduced to

[TABLE]

where

[TABLE]

and $f_{\rm in}=f_{\rm v}(R_{\rm in})$ is the viscous frequency at the inner disc radius. Equation (3.2.1) gives an analytic expression for the Green function in the frequency domain. In the limiting case of $R_{\rm out}\rightarrow\infty$, the expression (3.2.1) turns to (A), as expected.

3.2.2 The case of a non-zero-torque inner boundary

Whilst a zero-torque inner boundary condition may be typical for accretion discs around black holes, in the case of accretion onto magnetised NSs, the torque at the inner disc radius may have a finite value, which is determined by the the spin-up rate of the rotating magnetized NS. Stable accretion onto a magnetised NS is possible if the the mass accretion rate is high enough for the accretion flow to penetrate through the centrifugal barrier set up by the rotating magnetosphere (Illarionov & Sunyaev, 1975). Non-zero torque at the inner disc radius affects the torque distribution all over the disc.

Viscous torques in accretion disc are determined by the local surface density and viscosity. If the viscosity depends on radial coordinate only, the linear dependence between the local viscous torque $F$ and surface density $\Sigma$ holds:

[TABLE]

where $h=(GMR)^{1/2}$ is the specific angular momentum. The local mass accretion rate is determined by the derivative

[TABLE]

In the case of stationary accretion in a disc with a non-zero torque at the inner radius, the the viscous torques in the disc $F$ are given by

[TABLE]

where $F_{\rm in}$ is the torque at the inner disc radius, i.e. non-zero torque at the inner disc radius results in an increase of viscous torques all over the disc by a constant value. The time dependent torque can be represented as the sum of the time independent torque corresponding to stable accretion (23) and a fluctuating viscous torque $F_{2}$ on top of it:

[TABLE]

where $\dot{M}_{0}$ is the average mass accretion rate. The first two terms on the right hand side of (24) represent the time independent solution of the equation of viscous diffusion. Because the equation of viscous diffusion is considered to be linear, the fluctuating part of the solution is a solution of the viscous diffusion equation by itself. Note that $F_{2}(R,t)|_{R=R_{\rm in}}=0$ and, therefore, it satisfies the equation of viscous diffusion with zero torque boundary condition at $R_{\rm in}$. According to (22) the time dependent mass accretion rate can be obtained from (24):

[TABLE]

where the second term represents local fluctuations of the mass accretion rate on top of the average $\dot{M}_{0}$. Thus, the fluctuations of the mass accretion rate $\dot{m}(R,t)={\partial F_{2}(R,t)}/{\partial h}$ are described by a Green function derived for zero torque at the inner disc radius (see Section 3.2.1).

3.3 Properties of the transfer functions of propagating fluctuations

Fig. 3 shows an example of our Green function for the surface density (top; Equation 17) and mass accretion rate (bottom; Equation 3.2.1). We see that an initial perturbation in the surface density at $R^{\prime}=10$ spreads out over time, with more material propagating inwards than outwards. This creates a perturbation in the accretion rate that is initially positive inside $R^{\prime}=10$ and negative outside of $R^{\prime}=10$ before the accretion rate slowly settles back to its equilibrium value (zero in the Figure). In the frequency domain, the Green functions of the mass accretion rate $\overline{G}_{\dot{M}}(R,R^{\prime},f)$ play the role of transfer functions that describe how variability at Fourier frequency $f$ and radial coordinate $R^{\prime}$ in the disc affects variability at the coordinate $R$, see equation (16). In particular, the absolute value of the transfer function demonstrates the survival level of initial variability: $|\overline{G}_{\dot{M}}(R,R^{\prime},f)|=1$ corresponds to the case of fluctuations propagating without suppression, while $|\overline{G}_{\dot{M}}(R,R^{\prime},f)|\ll 1$ corresponds to significant suppression of initial variability. In this paper, we are interested in mass accretion rate variability at the inner disc radius and focus on the inward propagation of fluctuations (i.e. $R<R^{\prime}$).

The process of viscous diffusion effectively suppresses variability at frequencies $f\gtrsim f_{\rm v}(R^{\prime})$, where $f_{\rm v}(R^{\prime})$ is a local viscous frequency corresponding to the radial coordinate of initial fluctuations, unless the radial coordinate of initial fluctuations is very close to the inner disc radius (see Fig. 4). Fluctuations arising close to the inner disc radius can propagate inwards without significant suppression. Fig. 5 shows that the effect of finite $R_{\rm out}$ becomes important for the propagation of fluctuations that originated at large radii, and Fig. 6 shows that the effect is more pronounced for smaller $R_{\rm out}$ (as one would expect). Fig. 7 compares the fluctuations for three different Green function solutions: (a) $R_{\rm in}=0$ and $R_{\rm out}=\infty$, (b) $R_{\rm in}>0$ and $R_{\rm out}=\infty$, and (c) $R_{\rm in}>0$ and $R_{\rm out}<\infty$. We see that the boundary conditions affect the effectiveness of propagation of mass accretion rate fluctuations in the following ways (see Fig.7):

(i) Variability originating near the non-zero inner radius of the disc is weakly supressed and the inner disc regions contribute more high-frequency variability comparing to the case with $R_{\rm in}=0$ (compare Fig.7a and Fig.7b for $f/f_{\rm v}>10^{-1}$), and thus the inner disc regions contribute more high-frequency variability at the inner radius in the former case. This happens because there is no flow of angular momentum to $R_{\rm in}$ from $R<R_{\rm in}$ and the time scale of radial transfer becomes smaller.

(ii) The accretion discs with outer radius $R_{\rm out}<\infty$ with $\dot{M}(R_{\rm out})=0$ preserve some extra variability originating from the outer parts of accretion disc (see Fig.6, compare Fig.7b and Fig.7c for $f/f_{\rm v}<10^{-2}$).

3.4 Power spectra of initial perturbations

An essential ingredient of the model predicting timing properties of broadband aperiodic variability is the properties of initial fluctuations produced in the accretion disc. Starting from given properties of the initial fluctuations, we can describe their transfer properties using suitable Green functions of the viscous diffusion equation (see equation 16) to obtain predictions on the aperiodic variability in X-rays. However, the timing properties of the initial perturbations of the surface density and mass accretion rate are not known precisely. As discussed in Section 2.3.1, they may be caused by a magnetic dynamo process (Balbus & Hawley, 1991; Hawley et al., 1995; Brandenburg et al., 1995), which has typical timescale $t_{\rm d}\approx k_{\rm d}f_{\rm K}^{-1}$, where $k_{\rm d}\sim few$ (Tout & Pringle, 1992). Here, we specify the PDS of initial fluctuations with a Lorentzian function 111 The full information about the timing properties of initial perturbations is given by the cross-spectrum of initial variability at various radial coordinates in the accretion disc (Mushtukov et al., 2018). However, one can use the PDS of initial fluctuations $S_{a}(f)$ instead of the cross-spectrum in order to get an approximate solution (see equation 16). The approximation is good if different rings in the accretion disc produce uncorrelated initial perturbations of the mass accretion rate.

[TABLE]

where $F_{\rm var}$ is the fractional variability amplitude per radial decade generated by turbulence in the disc (following e.g. Arévalo & Uttley 2006; Ingram & Done 2011, 2012; Ingram & van der Klis 2013). This means that the integrated power of the perturbations $P(R)=\int_{0}^{\infty}S_{a}(f)df={F_{\rm var}}/{(\ln 10~{}R)}\propto 1/R$ for a constant $F_{\rm var}$, which is expected from MRI turbulence in a constant $h/r$ disc (e.g. Churazov et al. 2001; Arévalo & Uttley 2006). Since we are ignoring the influence of non-linear effects on the shape of the power spectrum, the $F_{\rm var}$ model parameter acts as a simple normalisation of the predicted XRP power spectrum.

The frequencies $f_{0}$ and $f_{\rm br}$ are model parameters, and we assume $f_{\rm br}(r)=f_{\rm d}(r)\propto f_{\rm K}(r)/k_{\rm d}$, where the dynamo coefficient $k_{\rm d}$ is a model parameter. The choice of Lorentzians is motivated by the fact that an exponentially decaying periodic signal

[TABLE]

is described by

[TABLE]

in the frequency domain, which gives the Lorentzian when squared.

4 Modeling the Power Density Spectra

Fluctuations of X-ray energy flux in XRPs are determined by fluctuations of the mass accretion rate at the inner edge of the accretion disc. Therefore, we model the PDS of fluctuations of the mass accretion rate at the inner disc radius and compare it with the PDS of observed fluctuations of X-ray energy flux. We use dimensionless radial coordinates and dimensionless frequencies in our calculations, i.e. the radial coordinate and the frequency are measured in units that set the physical scales in the system. In particular, the frequency is measured in units of the viscous frequency at the disc inner radius, which is related to viscous and geometrical properties of the accretion disc (see Section 2.3.1).

4.1 The major physical parameters and examples of PDS

The PDS of the mass accretion rate at the inner disc radius is determined by a set of parameters: (i) the inner and outer radii of the accretion disc and boundary conditions there (Fig. 8); (ii) the PDS of initial fluctuations (Fig. 9); (iii) the dependence of the kinematic viscosity on the radial coordinate in the disc (Fig. 10); (iv) the dependence of total power of initial fluctuation on the radial coordinate in the disc (Fig. 11). We see that the broadband variability covers a few orders of magnitude in the frequency domain. The typical theoretical PDS has two breaks. The break at lower frequencies corresponds to the viscous frequency $f_{\rm v}$ at the outer disc radius, while the break at high frequencies corresponds to the typical frequencies of initial fluctuations at the inner disc radius (see Fig. 8). Both breaks are smooth and can cover an order of magnitude in frequency, which poses difficulties for interpretation of observational results. The exact shape of the break at high frequency is strongly affected by the shape of the PDS of initial fluctuations (see Fig. 9).

The decrease of the inner disc radius (which can be caused by increasing mass accretion rate in XRPs, see equation 1) results in a shift of the high-frequency break to even higher frequencies and a decrease of the PDS at frequencies below the high-frequency break (compare black solid and blue dashed-dotted lines in Fig. 8). The shift of the break frequency is caused by the inclusion of a new inner part of the disc that is able to produce fast variability at time scales not achievable in an accretion disc with larger inner radius. The same inner part of the accretion disc additionally suppresses variability from the outer parts. It leads to a decrease of the PDS at low frequencies (compare black and blue lines of Fig. 8 in the frequency range $10^{-2}\lesssim f/f_{\rm v}(R_{\rm m})\lesssim 10^{1}$). The PDS above the high-frequency break deviates from a power-law.

4.2 Comparison with data

We compare the theoretical PDS of X-ray flux variability with the observed PDS in two luminosity states of the X-ray transient A 0535+26 (see Fig. 12), 222 The data for the plot are from RXTE/PCA, which observed A 0535+262 during its 2009 outburst. From each observation, light curves were extracted using the heasoft software package XRONOS. Power spectra were calculated from light curve segments of 512 s duration using the powspec tool with Miyamoto normalization and averaged together from observations within a luminosity interval. Poisson noise level was subtracted from the power spectra, which were then multiplied by frequency. The power spectra presented here correspond to observations with average luminosities of $1.7\times 10^{35}\,{\rm erg\ \rm s^{-1}}$ and $3.8\times 10^{36}\,{\rm erg\ \rm s^{-1}}$. which is one of the brightest XRPs on the X-ray sky located at distance around 2 kpc from the Sun (Steele et al., 1998). A cyclotron line scattering feature at $E_{\rm cyc,0}\sim 45\,{\rm keV}$ and its first harmonic at $E_{\rm cyc,1}\sim 100\,{\rm keV}$ have been detected (Caballero, 2007). Thus, A 0535+26 is a source for which we can estimate magnetic field strength at the NS surface: $B\sim 4\times 10^{12}\,{\rm G}$ and, therefore, the inner disc radius at any given mass accretion rate (see equation 1).

Using our model, we calculate the theoretical PDS of mass accretion rate fluctuations at the inner disc radius and compare it with the observed PDS of X-ray flux variability in two luminosity states of A 0535+26. The input parameters of our model are $n$, $R_{\rm m}$ [in cm], $R_{\rm out}$ [in cm], $f_{\rm br}(R_{\rm m})$ [in Hz], $f_{\rm br}(R_{\rm m})/f_{\rm v}(R_{\rm m})$, $f_{0}$ and $F_{\rm var}$. We use $n=3/4$ (typical for a gas pressure dominated $\alpha$-disc) and calculate the magnetospheric radius of the two luminosity states by combining the B-field measurement with the observed luminosities of the two states

[TABLE]

To obtain theoretical PDSs we choose the inner disc radii to be

[TABLE]

according to equation (1). We fix the outher disc radius at $R_{\rm out}=10^{10}\,{\rm cm}$, but note that the model in the frequency range explored is almost insensitive to this parameter. We set the the break frequencies at the inner radius to

[TABLE]

and $f_{\rm br}/f_{\rm v}=100$ at the inner disc radius (see Fig. 12). We fix $f_{0}=0$ for both observations and set $F_{\rm var}$, which acts simply as a normalisation constant of the PDS, separately for the two observations.

For the above values of $R_{\rm m}$, and assuming $M_{1.4}=1$, the Keplerian frequencies at the inner disc radii are expected to be

[TABLE]

In both cases this is a few times higher than the break frequency used in our modeling (see the gaps between the breaks and Keplerian frequencies in Fig. 12), giving $k_{d}\sim 9$ and $\sim 5$ for the low and high flux state respectively. The possibility of an underestimate of the inner disc radius (and corresponding overestimation of the Keplerian frequency) from relation (1) due to uncertainty in $\Lambda$ can be largely excluded because the inner disc radius in A 0535+26 was verified by measurements of NS spin evolution during a few outburst in 2009-2015 (Sugizaki et al., 2017), which have shown that the magnetic field is likely dominated by its dipole component and the inner disc radius is well described by the relation (1) with $\Lambda=0.5$. The estimations of the inner radius show that the disc in these particular cases is geometrically thin and gas pressure dominated in the considered time intervals (see Fig. 2). Therefore, the exponent $n$ in (10) can be taken to be $3/4$.

Thus, we have to conclude that the initial variability is generated at frequencies which are noticeably lower than the local Keplerian frequency. The inference that initial fluctuations generated at frequencies a few times below the local Keplerian frequency is consistent with the idea that initial fluctuations are driven by a dynamo process (King et al., 2004) resulting from the magneto-rotational instability (Balbus & Hawley, 1991). Thus, the modeling of XRP PDSs provides an opportunity to measure the time scales of the dynamo process in accretion discs observationally. However, these measurements require much more detailed analysis of PDS variability with X-ray luminosity, which is beyond the scope of this paper.

5 Summary and discussion

We have considered the physical processes responsible for the broadband aperiodic variability of X-ray energy flux in accreating highly magnetized NSs - XRPs. In the case of XRPs, emission in the X-ray energy band originates from small regions located close to the surface of the NS, while the accretion disc is truncated at a large distance by a strong magnetic field and thus does not produce an appreciable fraction of X-ray energy flux. 333With the exception of a narrow energy band, where the accretion disc produces an iron K$\alpha$ line due to the reflection of X-ray flux from the accretion disc (George & Fabian, 1991). Note that, the light crossing time of the inner radius is much shorter than the Keplerian time scale, and the variability of the reflected emission should replicate the variability of the accretion rate at the NS surface.

The observed aperiodic variability of X-ray flux in XRPs is caused by mass accretion rate variability in the vicinity of the NS surface, which replicates the variability of the mass accretion rate at the inner radius of the accretion disc. Because the majority of X-ray photons originate from a small region at the NS surface, the timing properties of aperiodic variability in X-rays are expected to be independent of the energy band. Possible dependencies of the timing properties on the photon energy are expected to be caused by further reprocessing of X-ray photons, which beyond the scope of this paper. Fluctuations of the mass accretion rate at the inner disc radius are shaped by fluctuations arising all over the disc and propagating inwards due to the process of viscous diffusion (Lyubarskii, 1997; Churazov et al., 2001; Kotov et al., 2001). The process of viscous diffusion modifies the timing properties of propagating fluctuations, suppressing the variability at high frequencies (Mushtukov et al., 2018). The diffusion process is determined by properties of viscosity in the accretion disc and boundary conditions: displacement of the inner and outer radii of the disc and conditions there (local torques and mass accretion rates).

We have developed a theoretical base for calculations of mass accretion rate variability at the inner radius of the accretion disc, which largely shapes the timing properties of XRPs. Using known analytical solutions of the equation of viscous diffusion (8), we have investigated the transfer properties of the discs. A new Green function of the viscous diffusion equation accounting for fixed inner and outer radii of the accretion disc has been derived (see Section 3.2). Our Green function generalizes the previous solutions obtained for accretion discs in the Newtonian approximation (Lynden-Bell & Pringle, 1974; Tanaka, 2011; Lipunova, 2015). The obtained Green function provides a more complete description of the propagation of mass accretion rate fluctuations under the influence of viscous diffusion in geometrically thin accretion discs.

XRPs are unique objects because they provide the possibility to probe timing properties of mass accretion rate variability within a small range of radial coordinates in the accretion disc (in contrast to accreting BHs, where the observer detects X-ray photons originating from the extended inner region of the disc, see e.g. Ingram & van der Klis 2013; Mushtukov et al. 2018). The geometry of the accretion disc (particularly, the inner disc radius) depends on the mass accretion rate. The timing properties of aperiodic variability, and particularly the PDS of the variability, depend on accretion disc geometry and, therefore, on the mass accretion rate. As a result, transient XRPs and analyses of their variability can be used to probe the geometry (its inner radius) of a disc and the physical conditions there: the time scale of the dynamo process, which is assumed to be responsible for the initial fluctuation of viscosity.

In the case of known strength and structure of NS magnetic field and, therefore, known inner radius of the disc at given luminosity (see Section 2.1), we can test physical conditions in the accretion disc. Particularly, one can measure the time scale of initial fluctuations of the mass accretion rate. Using the observed PDS of the transient X-ray pulsar A 0535+26 and modeling it with our theory, we conclude that the typical frequency of initial fluctuations is lower (by a factor of $\sim 5-9$) than the local Keplerian frequency (see Section 4.2). Note, that this conclusion directly follows from the assumption that the break in the PDS corresponds to the characteristic time scale at the inner boundary of the disk, which is not suppressed by the process of viscous diffusion. This statement is in agreement with the model, where the initial fluctuations are caused by the dynamo process, driven by the magneto-rotational instability in the accretion disc (King et al., 2004). The coefficient of proportionality between the Keplerian time scale and dynamo time scale is a matter of first-principal numerical simulations. According to our results, this coefficient can be obtained from observations of X-ray transients with known magnetic fields. This, however, requires detailed analysis of PDS variability with mass accretion rate and is beyond the scope of the paper.

It is worth noting that there are still a few open issues in the problem: (i) the exact timing properties of initial fluctuations of the mass accretion rate, and (ii) the timing properties of instabilities developing at the inner disc radius. These problems have to be addressed by numerical MHD simulations.

Acknowledgements

This research was supported by the Netherlands Organization for Scientific Research (AAM), the grant 14.W03.31.0021 of the Ministry of Education and Science of the Russian Federation (AAM and SST), RFBR grant 18-502-12025 (GVL), the Royal Society (AI), the Väisälä Foundation (SST). The authors would like to acknowledge networking support by the COST Actions CA16214 and CA16104. We are also grateful to Victor Doroshenko and Nikolai Shakura for discussion and a number of useful comments.

Appendix A Accretion disc with a finite inner radius

There are a few known analytical solutions of the equation of viscous diffusion. In order to describe the viscous evolution of the accretion disc in an XRP, one should account for a non-zero inner radius of a disc. Green functions for the case of the non-zero inner disc radius and infinite outer radius ($R_{\rm out}=\infty$) have been derived by Tanaka (2011). In the case of zero torque at $R_{\rm in}$ the Green functions of the surface density are given by

[TABLE]

where $x=(R/R_{\rm in})^{1-n/2}$, $t_{\rm v,\,in}=\frac{2}{3}R_{\rm in}^{2}/\nu(R_{\rm in})$ is the local viscous time at the inner disc radius,

[TABLE]

[TABLE]

and $J_{l}(x)$ and $Y_{l}(x)$ are the Bessel functions of the first and second kind. The corresponding Green function for the mass accretion rate $G_{\dot{M}}(R,R^{\prime},t)$ can be found according to equation (14). The Green function for the mass accretion rate in the frequency domain takes the following form:

[TABLE]

The Green functions given by (26) and (A) ignore the effects arising from the disc truncation at larger radii.

Appendix B Green functions for the disc with finite inner and outer radii

The equation of viscous diffusion in an accretion disc (8) can be rewritten in terms of the viscous torque

[TABLE]

and the specific angular momentum $h=(GMR)^{1/2}$ as follows:

[TABLE]

If the kinematic viscosity $\nu$ is a function of radial coordinate only, as it is determined by equation (10), the equation (28) is a linear diffusion equation. Its solution can be found in terms of Green functions, specified by the boundary conditions in the disc (Lynden-Bell & Pringle, 1974; Tanaka, 2011; Lipunova, 2015).

B.1 Zero torque at the inner disc radius

Tanaka (2011) found the form of basis functions which should be used to fulfill a boundary condition at a non-zero $R_{\rm in}$ (see his equation 26). Extending a particular solution of Lipunova (2015) to a case of $R_{\mathrm{in}}\neq 0$, we can construct a Green function $G_{F}(x,x_{1},t)$ for the torque $F$, valid for a disc with the finite inner and outer radii:

[TABLE]

where for $0<l<1$ and $l=1/(4-2\,n)$ we use the following combination of the Bessel functions:

[TABLE]

the viscous time at the outer disc radius $t_{\mathrm{v}}=2{R_{\rm out}^{2}}/{(3\nu(R_{\rm out}))}\,$, $x\equiv(R/R_{\rm out})^{(2-n)/2}$, and $k_{i}$ are the roots of the transcendent equation (notice that $h/h_{\rm out}=\xi=x^{2l}$)

[TABLE]

which is equivalent to

[TABLE]

and then to

[TABLE]

Equation (B.1) has infinite countable number of roots, which can be found numerically.

The solution of the viscous diffusion equation is given by

[TABLE]

Condition (31) expresses the homogeneous outer boundary condition on the accretion rate since $\dot{M}=\partial F/\partial h$. Notice that, by setting $x_{\mathrm{in}}=0$, expression (29) is reduced to a Green function found in in Lipunova (2015).

The Green functions for the surface density $G(R,R^{\prime},t)$ satisfy equation

[TABLE]

Equation (34) can be rewritten as

[TABLE]

or

[TABLE]

Thus, the Green functions describing the evolution of the surface density are given by

[TABLE]

which results in

[TABLE]

The Green function for the accretion disc of fixed inner and outer radii is expressed by the infinite series (29), where each term contains a root of equations (B.1). Numerically, we use a limited number of terms in (29). The smaller the ratio $(R_{\rm out}/R_{\rm in})$, the smaller the number of terms in a series, which are necessary to reach a certain accuracy.

B.2 Zero mass accretion rate at the inner disc radius

In the case of the zero mass accretion rate at the inner disc radius $R_{\rm in}$ and, therefore, on the NS surface, which is appropriate for the case of so-called ”dead discs” (Syunyaev & Shakura, 1977) in the propeller state of accretion, the Green function for the torque is modified:

[TABLE]

where the function $V_{*}(u,v)$ is defined as

[TABLE]

and $k_{i}$ are roots of transcendent equation similar to (31), but written for the function $V_{*}(u,v)$.