Magnetospheric Flows in X-ray Pulsars I : Instability at super-Eddington regime of accretion

Within the magnetospheric radius, the geometry of accretion flow in X-ray pulsars is shaped by a strong magnetic field of a neutron star. Starting at the magnetospheric radius, accretion flow follows field lines and reaches the stellar surface in small regions located close to the magnetic poles of a star. At low mass accretion rates, the dynamic of the flow is determined by gravitational attraction and rotation of the magnetosphere due to the centrifugal force. At the luminosity range close to the Eddington limit and above it, the flow is additionally affected by the radiative force. We construct a model simulating accretion flow dynamics over the magnetosphere, assuming that the flow strictly follows field lines and is affected by gravity, radiative and centrifugal forces only. The magnetic field of a NS is taken to be dominated by the dipole component of arbitrary inclination with respect to the accretion disc plane. We show that accretion flow becomes unstable at high mass accretion rates and tends to fluctuate quasi-periodically with a typical period comparable to the free-fall time from the inner disc radius. The inclination of a magnetic dipole with respect to the disc plane and strong anisotropy of X-ray radiation stabilise the mass accretion rate at the poles of a star, but the surface density of material covering the magnetosphere fluctuates even in this case.


INTRODUCTION
X-ray pulsars (XRPs) are accreting strongly magnetised neutron stars (NSs) in compact binary systems (see Mushtukov & Tsygankov 2022 for recent review).Magnetic field at the NS surface in XRPs confirmed by detection of the cyclotron lines (Staubert et al. 2019) is known to be of the order of 10 12 G or even stronger.Large scale magnetic field is expected to be dominated by the dipole component (that can be distorted by the accretion process, see Lai 2014 for review), but in some specific cases there was reported possibility of the field heavily contributed by non-dipole components (see, e.g., Tsygankov et al. 2017;Israel et al. 2017a;Mönkkönen et al. 2022).Observed luminosity of XRPs covers more then seven orders of magnitude from 10 32 erg s −1 up to 10 41 erg s −1 .The lower edge of accretion luminosity range is partly related to the onset of the "propeller" state of accretion when the cen-⋆ E-mail: alexander.mushtukov@physics.ox.ac.uk (AAM) trifugal barrier due to the rotating NS magnetosphere prevents accretion onto the stellar surface (Illarionov & Sunyaev 1975;Raguzova & Lipunov 1998; Ustyugova et al. 2006).The brightest XRPs belong to the class of recently discovered pulsating ultra-luminous X-ray sources (ULXs, Bachetti et al. 2014;Israel et al. 2017a;Fürst et al. 2016;Israel et al. 2017b;Carpano et al. 2018 andFabrika et al. 2021;King et al. 2023 for review).The apparent luminosity of the brightest ULX pulsar -NGC5907 X-1 -exceeds 10 41 erg s −1 (Israel et al. 2017a).The relation between the actual and apparent luminosity of bright XRPs is still debated by the community.Specific geometry of X-ray emission region in close proximity to a NS surface does not affect significantly the apparent luminosity of bright XRPs (see, e.g., Markozov & Mushtukov 2024).However, there are arguments in favor of a large difference between actual and apparent luminosity due to the geometrical beaming at high mass accretion rates (King 2009;King et al. 2017;King & Lasota 2020) expected because of possibly strong radiation driven outflows launched from accretion discs (Shakura & Sunyaev 1973;Poutanen et al. 2007;King 2009;Jiang et al. 2014).At the same time it has been shown that strong beaming of X-ray radiation is inconsistent with large pulsed fraction observed in six ULX pulsars known up to date, which says in favour of relatively small difference between the actual and apparent luminosity (Mushtukov et al. 2021;Mushtukov & Portegies Zwart 2023).
Accretion flow in the systems hosting strongly magnetised NSs is truncated by magnetic field of a NS at the magnetospheric radius Rm, which is determined by the NS magnetic field strength, mass accretion rate at Rm and geometry of the accretion flow (i.e., disc/wind accretion, disc thickness and physical conditions in it, see Chashkina et al. 2019 for detailed discussion).At the magnetospheric radius, the flow penetrates into the magnetosphere due to the instabilities, including magnetic Rayleigh-Taylor (Arons & Lea 1976;Kulkarni & Romanova 2008) and Kelvin-Helmholtz instabilities (Burnard et al. 1983).Development of the instabilities results in formation of a boundary layer between the disc and magnetosphere of a NS (Lai 2014).In the boundary layer, the angular velocity of plasma experiences a transition from the Keplerian rate to the NS spin rate.
Within the magnetospheric radius, the flow follows magnetic field lines and finally reaches the surface of a NS within small regions (area ∼ 10 10 cm 2 , see Section 4.3 in Mushtukov & Tsygankov 2022) located close the magnetic poles of a star.Depending on the mass accretion rate and corresponding energy release, the flow is stopped due to the Coulomb collisions in the atmosphere of a NS (expected for Ṁ ≲ 10 17 g s −1 , see, e.g., Zel'dovich & Shakura 1969;Nelson et al. 1995) or at the radiation pressure dominated shock above NS surface (expected for the case of Ṁ ≳ 10 17 g s −1 , see Basko & Sunyaev 1976;Wang & Frank 1981;Zhang et al. 2022).
The dynamics of accretion flow in between the magnetospheric radius and NS surface is determined by the geometry of magnetic field lines, gravitational force, centrifugal force and radiative force, which comes into play at sufficiently high mass accretion rates.Stable accretion flow covering the magnetosphere of a NS is expected to be optically thick at high mass accretion rates Ṁ ≳ few × 10 18 g s −1 typical for ULXs pulsars and bright Be XRPs (Reig 2011) at the peaks of their outbursts (see, e.g., Mushtukov et al. 2017Mushtukov et al. , 2019a)).The optically thick envelope might strongly affect the key observational properties of bright XRPs, including their energy spectra, X-ray polarisation, pulse profiles and timing properties of aperiodic variability.
This paper presents the first simulations of accretion flow dynamics between the inner disc radius and NS surface.We account for the geometry of accretion flow shaped by a magnetic field of a NS, and the influence of gravitational, centrifugal and radiative forces.The stellar magnetic field is assumed to be purely dipole with a given arbitrary inclination of the magnetic dipole with respect to the accretion disc plane.

Geometry of accretion flow
In this paper, we consider a particular case of accretion from the disc, i.e. the material starts its motion towards the poles of a NS from a specific regions, where the disc is interrupted by the B-field.We assume that the rotational axis of a NS is orthogonal to accretion disc plane, which is likely a case in strongly magnetised accreting NSs.The disc in XRPs is interrupted at the magnetospheric radius Rm, which is dependent on the mass accretion rate, NS magnetic field strength and structure.In the case of the field dominated by the dipole component, Rm can be estimated as where Λ < 1 is a factor depending on the accretion flow geometry (see, e.g., Chapter 6.3 in Frank et al. 2002, andChashkina et al. 2019 for recent discussion related to the case of high mass accretion rates), B12 is a surface magnetic field strength in units of 10 12 G, Ṁ17 is mass accretion rate in units of 10 17 g s −1 , m is a mass of a NS M in units of M⊙, and R6 is the NS radius in units of 10 6 cm.The effective temperature in accretion disc can be estimated as where σSB is the Stefan-Boltzmann constant and r is a distance from a central compact object, and β ∈ [0; 1] is a coefficient determined by a distance where stress disappears (β = 1 in the case of stress disappearance at Rm, while β = 0 corresponds to stress disappearance at r ≪ Rm).Combining ( 1) and ( 2) we can estimate roughly the effective temperature at the inner disc radius as The temperature of plasma settling magnetic field lines at the inner disc radius is probably higher because (2) and (3) do not account for the interaction between the magnetosphere of a NS and accretion disc.We assume that the accretion flow within the magnetospheric radius cannot move across magnetic field lines.In the spherical coordinates (r, θ, φ), the geometry of a dipole field line is given by where Rmax is a linear scale describing dipole magnetic field line, λ is the latitude, and θ = π/2 − λ is corresponding colatitude in the reference frame related to the B-field.The linear scale of a field line Rmax in ( 4) is typically assumed to be equal to Rm, but it is a case of magnetic dipole aligned with the disc only.In the case of inclined magnetic dipole where is the co-latitude, where the disc plane crosses the dipole surface, α is the magnetic obliquity, i.e. the angle between the rotational axis of the NS and the magnetic dipole axis (see Fig. 1).Assuming the inner disc radius independent of the azimuthal coordinate, we get Rmax dependent on φ and magnetic obliquity α.
The area on the dipole surface S is related to the latitude λ and the azimuthal angle φ as  while the length along the field lines x is related to the latitude as The angle between the position vector of the point at the dipole surface and the tangent line to the dipole magnetic field line (see Fig. 2) is Further in this paper we present maps describing the distribution of different quantities (accelerations, surface density, etc.) over the dipole surface of a NS magnetosphere.The position of a point on the dipole surface (Rmax is fixed) is uniquely given by two angles: λ (or θ) and φ.To represent magnetospheric maps, we apply Aitoff transformation to the coordinates λ and φ, as is standardly done to represent maps of spherical surfaces: where ϱ = arccos cos λ 2 cos φ .

Dynamics of accretion flow
We consider the accretion flow moving under the influence of three forces: the gravitational, the radiative and the centrifugal force.We assume that the flow covering the magnetosphere of a NS is supersonic and do not account for the influence of internal gas pressure gradients.The gravitational force is directed towards the compact object, and the corresponding acceleration along the dipole field line is given by The gravitational acceleration along the field lines turns to zero at the equator.At λ > 0 (λ < 0), the acceleration a grav,|| is positive (negative).Its absolute value of the acceleration increases towards to the central object.
The centrifugal force is directed perpendicularly to the rotational axis of a NS and dependent on the NS spin period and specific point on the magnetospheric surface.In the particular case of magnetic dipole aligned with the rotational axis, the corresponding centrifugal acceleration along the field lines is given by where ω is the angular velocity of a NS.In a general case of inclined magnetic dipole, (12) has to be calculated as where nB is the unit vector directed along the magnetic field line, "•" and "×" denote the scalar and the vector productions respectively.In the particular case of the orthogonal rotator (i.e., α = π/2), equation ( 13) is reduced to The total acceleration along the field lines due to the gravitational and centrifugal forces is dependent on the inclination of the magnetic dipole with respect to the rotation axis and NS spin period.In general, the acceleration is dependent both on the latitude and the azimuth angle (see Fig. 3).Sufficiently small spin periods result in the appearance of regions with acceleration directed towards the disc plane, i.e. total acceleration due to the gravitational and centrifugal forces tends to prevent accretion flow from its motion towards the poles of a NS.It is the centrifugal barrier, which can lead to the propeller regime of accretion in XRPs (Illarionov & Sunyaev 1975;Ustyugova et al. 2006).It is expected that the accretion flow is able to penetrate through the centrifugal barrier in the case of the sufficiently high initial velocity of the flow or in the case of a sufficient geometrical thickness of the disc (some aspects of the influence of accretion disc geometrical thickness on the propeller regime were discussed by Çıkıntoglu & Ekşi 2022).In the case of a magnetic dipole inclined to the rotation axis, the acceleration becomes dependent on the azimuthal coordinate on the dipole surface φ (see Fig. 3).At a sufficiently small spin period and sufficiently large inclination of the magnetic dipole with respect to the disc plane, there are regions at the dipole surface where material can stay stably, i.e. the acceleration due to the gravitational and centrifugal forces is directed towards these regions (appearance of these regions was discussed by Abolmasov &Biryukov 2020 andLyutikov 2023).The inclination of the magnetic dipole, initial accretion flow velocity, NS spin period and thickness of the disc determine the regions where the material is collected.
Interaction of magnetospheric flow with X-ray photons emitted by the central objects results in appearance of radiative force.The radiative force is determined by local Xray energy flux at the magnetospheric surface, the fraction of intercepted radiation and local inclination of magnetospheric surface in respect to the predominant direction of photons momentum.In the first approximation, the radiative force is directed in oppositely to the gravitational force.The fraction of intercepted X-ray photons is determined by the surface density of magnetospheric flow and mechanism of opacity.In our simulations, we assume that the opacity is dominated by the non-magnetic Compton scattering and taken to be κe = 0.34 cm 2 g −1 .Non-magnetic opacity is a valid approximation because magnetic field strength decreases rapidly with a distance from a NS and the cyclotron energy Ecyc ≃ 11.6 B12 keV on the most of magnetospheric surface is much smaller then the typical energy of photons.In the limiting case of total absorption of radiation, the radiative force along the dipole field lines acting on a unit area can be estimated as where F (λ) is the photon energy flux.In the real situation, only a fraction of photons are absorbed/scattered by the envelope and contribute to be radiative force.Thus, (15) has to be rewritten as where Σ is the local surface density of the magnetospheric flow.The estimation ( 16) assumes that the scattering is isotropic and, therefore, the photons effectively transfer their momentum to the accretion flow due to the first scattering event only.The equations ( 15) and ( 16) can be rewritten via the radiation force acceleration to keep the uniform with the previous description of the forces.The radiation acceleration inside the accretion flow is where m is Lagrangian coordinate of the point inside the flow.It is determined as dm = ρ dz, where z is a coordinate along the local normal to the magnetospheric surface, and ρ is the plasma mass density.The radiation acceleration averaged across the accretion flow is where Σ is a total surface density integrated across the accretion flow at the given position.The last equation coinsides with ( 16), if we take that a rad (λ, Σ) Σ = dfR/dS.In average, the photon scatterings following the first one do not affect the dynamics of the flow.
In the case of optically thick flow (i.e., Σκe > 1), the photon energy flux varies over the optical depth: some photons are scattered and reflected back by the surface layers, and only a fraction of photons penetrates deeper inside the flow covering magnetic dipole.Therefore, the acceleration due to the radiative force is dependent on the optical depth, and the layers located closer to the central source experience stronger accelerations.It might result in a gradient of velocity across the accretion channel.In this paper, we do not account for the influence of the radiative force gradient and possible effects of viscosity inside the accretion channel.

Radiation from the neutron star surface
NS radius is expected to be much smaller than its magnetosphere (see Eq. 1), which is a good approximation in the case of XRPs.Thus, accounting for the influence of radiative force on the accretion flow dynamics, we consider a NS as a point source of X-ray radiation.Because accretion flow is directed towards NS magnetic poles by the stellar magnetic field, the emission from the central object can be strongly anisotropic.The angular distribution of X-ray luminosity is not known and described in our model by a special parameter asp ∈ [−1; +1] in our model: where index k ∈ {1, 2} denotes the pole of a NS, and θB,i ∈ [0, π] is the colatitude in the reference frame related to the magnetic dipole.Parameter asp = 0 corresponds to the case of isotropic radiation from the central object, asp = −1 corresponds to the fan beam diagram, and asp > 0 corresponds to pencil beam diagram.The beam pattern is shaped by the geometry of the emitting region at the NS surface (Gnedin & Sunyaev 1973) and can be strongly influenced by the gravitational bending in the vicinity of a NS (Riffert & Meszaros 1988;Kraus 2001;Mushtukov et al. 2018).The total luminosity of the central object Ltot = L1 + L2 can be variable and dependent on the current mass accretion rate onto the poles of a NS.The mass accretion rate on one pole can be different from the mass accretion rate on another pole, which results in a difference between L1 and L2 and the corresponding X-ray energy flux at the magnetospheric surface: where r is determined by (4).

Mass accretion rate at the inner disc radius
The average mass accretion rate at the inner disc radius is a parameter of our simulations.It can be considered to be constant or variable due to fluctuating mass accretion rate in the disc (see, e.g., Mushtukov et al. 2019b).Stochastic variability of the mass accretion rate at the inner disc radius is described by the power density spectra (PDS) and rootmean square (rms) variability.The broad-band component of the PDS in X-ray binaries is generally well described by a twice-broken power law (Hoshino & Takeshima 1993).In this paper, we restrict ourselves to investigating the effect of the rms on the magnetospheric accretion and, for simplicity, assume that the PDS is given by a power law in a given range of Fourier frequencies f , i.e.
where q is fixed at (−1).Because the disc can be geometrically thin or thick depending on the mass accretion rate (see estimations of the disc thickness in Suleimanov et al. 2007), the particles participating in our simulations start their motion along dipole magnetic field lines at λ0 ̸ = λ d = π/2 − θ d .

NUMERICAL MODEL
We consider a one-dimentional motion of accretion flow along dipole field lines given by (4).In the case of the rotational axis orthogonal to the disc plane (i.e.α = 0), the accretion flow starts its motion at the latitude where δλ is determined by the semi-thickness of the accretion disc at the magnetospheric radius and reaches the NS surface at coordinate The accretion channel in between coordinates λ0 and λNS is divided into N λ intervals in λ-coordinate and Nφ intervals in the azimuthal angle.As a result, we construct a twodimensional grid composed of N λ × Nφ cells covering the magnetospheric surface.For the simulations represented in this paper, we have taken N λ = 630 and Nφ = 24.Each simulation covers time interval of 30 s, is composed of a series of time steps and traces the motion of quasi-particles along the field lines.
At each time step of the simulation, we have information about the coordinates of all particles, their velocities and the total mass within each cell of a grid.At each time step, we perform the following procedures: (i) We start by determining the mass accretion rate at the inner disc radius based on the average mass accretion rate and rms of mass accretion rate variability.The mean mass accretion rate and its rms are parameters in our simulations.We follow the algorithm proposed by Timmer & Koenig (1995) to get a time series of the mass accretion rate and thus the mass accretion rate at a given moment.In our calculations, we take q = −1 and simulate the mass accretion rate time series on the base of assumed rms and PDS in the Fourier frequency interval 0.1 Hz < f < 100 Hz.To model the time series, we take 5 × 10 3 frequencies in the considered frequency range (for details see Timmer & Koenig 1995).(ii) We add Np particles of equal mass to each of Nφ grid cells located near the edge of accretion disc, i.e. we add NpNφ particles.The mass of particles added to the grid at a given time step ti is related to the mass accretion rate Ṁ at the inner disc radius as where ∆ti is a time interval of the simulation and index j enumerates particles participating in a simulation.In the simulations represented in this paper, we use Nφ = 24 and Np = 2.The particles are injected to the accretion flow at a random position within the cell located at the accretion disc plane, i.e. the initial coordinates are where Xj ∈ (0; 1) is a random number.The initial velocity of each particle along magnetic field lines vini is a parameter of the simulation.It is expected to be close to the thermal velocity of protons at the inner disc radius: where T keV is the plasma temperature at Rm in units of keV.In our simulations we use vini ∼ 10 7 cm s −1 .(iii) On the base of known coordinates of particles and their masses, we get the total mass in each cell of the grid and local surface density of the magnetospheric accretion flow Σ(λ, φ, ti).(iv) We calculate the time interval towards the new time step.The time step of the simulation is variable and dependent on the typical particle velocities in the cells.
It is taken to be sufficiently small to prevent particle transitions to the cells further than the nearest ones ∆ti = min min where vj is a velocity of a particle j, and the minimum is taken over all particles currently participating in simulation.
(v) We get the acceleration of particles along magnetic field lines.The particle motion is calculated by accounting for gravitational, radiation and centrifugal forces (see Section 2.2).We calculate the vector sum of three forces for each cell of the grid and get acceleration along field lines for each particle currently participating in the simulation.The gravitational and centrifugal accelerations are not dependent on time and can be pre-calculated.
The radiative force, on the contrary, can be variable due the variations of accretion luminosity and local surface density of accretion flow covering NS magnetosphere.(vi) On the basis of calculated accelerations for particles, we get their velocities and coordinates λj(ti+1) = λj(ti) + ∆λj(ti), ( 29) (vii) We update the list of particles participating in the simulation.In particular, we account for particles accreted to the NS magnetic poles: particles that have reached one of the magnetic poles are taken into account when calculating the mass accretion rate and are removed from the list of particles.The mass accretion rate onto a NS magnetic pole is calculaterd as where k ∈ {1, 2} denotes one of the magnetic poles and the summation is performed over the particles that have reached the surface in the time interval ∆ti.We also account for particles that have returned to the accretion disc.These particles remain in the simulation but are restarted from a random azimuthal coordinate at the next time step.The mass accretion rate calculated in (30) and final velocity of particles are used to get accretion luminosity where γj = [1 − (vj/c) 2 ] −1/2 is the final Lorentz factor of a particle.The luminosity ( 31) is used to calculate photon energy flux (20) at the magnetospheric surface and then the radiative force at the next time step of simulation.(viii) We return to step (i).
Particle motion changes the mass density in each grid cell, and, as a result, the influence of three forces on the particles in the cells is recalculated at each time step.
To investigate features of the mass accretion rate variability in the frequency domain, we redact the time series obtained by (30), where the time intervals between element of a time series are variable and short (see equation 27), which results in a Poisson noise.In particular, we combine nearby elements of the simulated time series into one and get an updated time series with larger and homogeneous time ∆t ≃ 2 × 10 −3 s.Using the updated time series of the mass accretion rate at the inner disc radius and at the NS surface, we calculate its Fourier transform: where f is the Fourier frequency and Nt is a number of simulated points in the time series.The power density spectrum is calculated using normalisation: where Macc = N t j=1 Ṁk (tj)∆t is the total mass accreted onto a NS during the time interval covered by a time series.Calculating the Fourier transform (32) and PDS (33) we exclude the first second of the simulation which is required for the accretion flow to reach the NS surface.The Nyquist frequency for the simulated time series is fN = 0.5/∆t ≃ 250 Hz and further we discuss PDS shape below this frequency only.In the case of normalisation given by ( 33), the noise caused by finite number of particles participating in a simulation is expected at the level where Ntot is a total number of particles accreted onto the stellar surface during the entire simulation.
To investigate stability of magnetospheric coverage by the accretion flow, we calculate the local average surface density, local averaged squared surface density and its local standard deviation

RESULTS OF NUMERICAL SIMULATIONS
In our simulations of accretion flow dynamics, we initialize the flow at t = 0 from the edges of the accretion disc, which is fixed radius Rm = 10 8 cm.It takes a fraction of a second (approximately 0.1 s for the specified Rm) for the accretion flow to reach the NS surface.This duration is slightly longer than the free-fall time from the inner disc radius due to the flow following curved magnetic field lines.The radiative force comes into play as soon as the flow reaches the stellar surface.Before that, the flow dynamic is shaped by the gravitational and centrifugal forces only.In this section, we first discuss the results of simulations performed for a low mass accretion rates from the disc ( Ṁ ≲ 10 17 g s −1 ), when the influence of the radiative force can be neglected (Section 4.1).In the following we discuss the effects introduced by strong radiative force at high accretion rates ( Ṁ ≳ ṀEdd ∼ 10 18 g s −1 , Section 4.2).

The influence of gravitational and centrifugal forces
The accelerations along B-field lines due to the NS gravity and the centrifugal force are dependent on the mass of a NS, spin period, size and inclination of the magnetic dipole with respect to the disc plane.The maps of this acceleration are .The mass accretion rate onto the NS surface calculated for different constant mass accretion rates at the inner disc radius: Ṁ (Rm) = 10 17 g s −1 (∼ 0.1 ṀEdd , black line), 10 18 g s −1 (∼ ṀEdd , blue line), and 10 19 g s −1 (∼ 10 ṀEdd , red line).In the case of low mass accretion rates (black line), the mass accretion rate at the NS surface replicates the mass accretion rate at the inner disc radius.Sufficiently high mass accretion rate from the disc results in appearance of instability of accretion process over the magnetospheric surface, which makes mass accretion rate at the NS surface strongly variable (blue and red lines).At the mass accretion rates ∼ 10 19 g s −1 from the disc (red line), the mass accretion rate at the stellar surface shows quasi-periodic oscillations.Parameters: m = 1.4,Rm = 10 8 cm, P spin = 10 s, isotropic central source, and α = 0.
pre-calculated (see Fig. 3).Sufficiently fast rotation of a NS magnetosphere results naturally in arising of the centrifugal barrier, i.e. the regions at the magnetospheric surface, where the local acceleration is directed back to the disc plane.The centrifugal barrier is stronger in the case of faster rotation of a NS.It is expected that the material from the disc can penetrate through the centrifugal barrier in the case of the sufficiently high initial velocity of the material (in this case, the barrier slows down the flow but cannot stop it) or sufficiently large geometrical thickness of accretion disc, when the matter is collected from the upper layers of the disc (the latest effect was discussed by Çıkıntoglu & Ekşi 2022).
In the case of an inclined magnetic dipole, there are regions at the inner disc radius where the centrifugal and gravitational forces prevent material motion towards one of the magnetic poles of a NS.It results in a partial covering of NS magnetosphere by the flow from the disc (see Fig. 4).Additionally, in the case of the sufficiently fast rotation of inclined magnetic dipole, there are regions at the NS magnetosphere where the projections of gravitational and centrifugal forces cancel each other.The appearance of these regions of stable equilibrium may result in local accumulation of accreting material (see Fig. 4).The fate of the matter locally accumulated at the magnetospheric surface might be affected by the cooling and heating processes.In particular, sufficiently fast cooling of matter accumulated at the magnetosphere can result in events of episodic accretion (Shakura et al. 2012).Analyses of heating and cooling processes, however, require specific analyses, which are out of the scope of this paper.
According to our simulations, the coverage of the NS magnetosphere is dependent both on the inclination of the magnetic dipole and NS spin period: the larger the inclination and the faster the rotation of a NS, the smaller the area of magnetosphere covered by accretion flow (see Fig. 4).We note that the arising of the "propeller" effect in our simulations is dependent on the inclination of magnetic dipole: the spin pe-riod of 0.5 s is small enough to stop accretion in the case of α = 0 • , but still allows accretion at α = 30 • and α = 60 • (see Fig. 4).

Radiative force coming into play
Influence of the radiative force depends on the mass accretion rate and corresponding energy release at the NS surface.At low, sub-Eddington ( Ṁ ≪ 10 18 g s −1 ), mass accretion rates, the radiative force does not play a significant role and stable mass accretion rate at the inner disc radius is replicated (with a time delay) at the stellar surface (see black line corresponding to Ṁ = 10 17 g s −1 in Fig. 5).High mass accretion rates, however, result in high luminosity and radiative force comparable to the gravitational and centrifugal forces.It causes instability of accretion flow at the mass accretion rates above the Eddington value ( Ṁ ≳ 10 18 g s −1 , see blue and red curves in Fig. 5 corresponding to the mass accretion rates 10 18 g s −1 and 10 19 g s −1 from the disc respectively) and sharp increase in rms of the mass accretion rate variability at the NS surface.Sufficiently intensive mass inflow at the inner disc radius ( Ṁ ≳ 10 19 g s −1 ) leads to quasi-periodic oscillations (QPOs) of the mass accretion rate at the stellar surface (see red line in Fig. 5).The time scale of the QPOs is close to the free-fall time scale: The mechanism standing behind the QPOs is related to the influence of radiative force on the accretion flow and can be explained as follows: (i) accretion flow reaches NS surface and results in high luminosity of the central object; (ii) high luminosity leads to the radiative force, which is large enough to stop the flow on its way towards a star, (iii) mass accretion rate onto a NS drops as well as the accretion luminosity, which leads to a drop of the radiative force affecting the Different panels show the PDS for different mass accretion rates from the disc: 10 19 g s −1 , 5 × 10 18 g s −1 , 2 × 10 18 g s −1 , and 10 18 g s −1 (from top to bottom).Parameters: Rm = 10 8 cm, P spin = 10 s, isotropic source of emission.
magnetospheic flow motion, (iv) under condition of small radiative force, accretion flow continues its motion towards a NS and reaches stellar surface within the free-fall timescale, (v) the cycle returns to the first step and continues quasiperiodically.The mass accretion rate onto one magnetic pole can differ from the mass accretion rate onto the other one.However, there is still a correlation between the two mass accretion rates.
In a high luminosity state, the PDS of the simulated time series of the mass accretion rate at the NS surface shows QPO peaks and its harmonics (see Fig. 6abc).In the range of mass accretion rate 2 × 10 18 g s −1 ≲ Ṁ ≲ 5 × 10 18 g s −1 , QPO appears at a frequency ∼ 7 Hz and is followed by a series of harmonics (see vertical dashed lines in Fig. 6bc).The frequency is almost independent on the mass accretion rate.The frequency of this QPO corresponds to the time interval that is slightly shorter than the free-fall time scale (37).This dis-crepancy may be attributed to the non-zero initial velocity of the flow and the non-zero geometric thickness of the accretion disc.At mass accretion rates exceeding ∼ 5 × 10 18 g s −1 ), we see a transition where the QPOs with a frequency of ∼ 7 Hz and its harmonics are gradually replaced by QPOs with a frequency of ∼ 5 Hz and its associated harmonics (see grey vertical shaded lines in Fig. 6a).This phenomenon is likely attributable to the diminishing radiative force, which decreases but does not entirely vanish, thereby continuing to influence the flow dynamics even during the intervals between quasiperiodic flares.The quality factors of the most distinct QPOs at high mass accretion rates reach a value of ≳ 40.At relatively low mass accretion rate (∼ 10 18 g s −1 , see Fig. 6d), PDS of the mass accretion rate shows a broad hump at high frequencies (≳ 20 Hz).
The strength of QPOs of the mass accretion rate due to the influence of the radiative force is strongly affected by the inclination of NS magnetic dipole in respect to the accretion disc plane and beam pattern of X-ray radiation produces at the NS surface.Both inclination of magnetic dipole (see Fig. 7) and beaming of radiative pattern (see Fig. 9) stabilise the mass accretion rate at the stellar surface.In the case of accretion onto the inclined magnetic dipole, stabilisation of accretion results in desapearance of QPOs from the PDS, but PDS shows broadband high-frequency component (see Fig. 8).
Stabilisation of accretion flow for the inclined magnetic dipole is attributed to the partial coverage of NS magnetosphere by the flow (see Fig. 4).Under this condition, a fraction of X-ray photons can freely leave a system without interaction with the magnetospheric flow.In addition, the flow covering NS magnetosphere has larger surface density Σ in comparison to the case of magnetic dipole aligned with accretion disc.The higher surface density results in a reduction of the average acceleration due to the radiative force, as described by equation (see equation 18).
Analyzing the relative fluctuations in surface density across the magnetospheric surface (i.e., its local average value given by equation 35 and standard deviation given by equation 36), it becomes evident that these fluctuations intensify toward the central object, as illustrated in Fig. 10.In the case of an inclined magnetic dipole, the strength of these fluctuations also exhibits azimuthal dependence.Specifically, surface density fluctuations are more pronounced towards the edges of the accretion curtain enveloping the magnetosphere, as highlighted in the lower panel of Fig. 10.It is essential to emphasize that these fluctuations in surface density hold the potential to significantly influence the process of pulse formation.A thorough exploration of this effect will be the focus of an upcoming publication.
Fluctuations of the mass accretion rate at the inner disc radius (caused by, in particular, the stochastic nature of viscous diffusion) do not affect qualitatively the process of magnetospheric accretion.The fluctuations at Rm increase rms of mass accretion rate variability at the NS surface.In the case of the average mass accretion rate close to the Eddington value, the fluctuating mass inflow at the magnetospheric boundary might result in variable rms of the mass accreiton rate at the NS surface (see Fig. 11).It happens because the current mass accretion rate at Rm can be occasionally above the Eddington limit due to the fluctuations and then the magnetospheric accretion process turns into the unstable mode.At the mass accretion rates much larger than the Eddington limit (we have tested the average mass accretion rate 10 19 g s −1 ), the QPOs are preserved when fluctuations of the accretion rate at the inner disc radius are taken into account (see Fig. 12).A detailed analysis of accretion timing and calculation of mass accretion rate PDS requires, however, the simulation of long time series and will be discussed in a separate publication.
The initial velocity of the accretion flow at the inner disc radius (within the range 10 6 cm s −1 ≲ vini ≲ 10 7 cm s −1 ) has only a negligible impact on the simulated time series of the mass accretion rate and its PDS.Note that the strength of the correlation between the mass accretion rate onto one and another magnetic pole can be dependent on the details of beam pattern formation and affected by light bending in the vicinity of a NS.We are planning to investigate these effects in a separate publication.

DISCUSSION AND SUMMARY
We conducted numerical simulations to explore the dynamics of accretion flow over the magnetosphere of bright X-ray pulsars.We assumed magnetic field geometry dominated by the dipole component and have taken into account the influence of the gravitational, centrifugal and radiative forces.The magnetic dipole was taken to be inclined with respect to the plane of the accretion disc, while the rotation axis of a NS was assumed to be orthogonal to the disc, which is a good approximation for the case of accretion onto strongly magnetised NSs, where the surface field strength B0 ≳ 10 11 G.These assumptions find support in the determination of geometrical parameters from several XRPs using polarimetric data from the IXPE observatory (see e.g.Doroshenko et al. 2022;Tsygankov et al. 2023;Mushtukov et al. 2023;Malacaria et al. 2023;Suleimanov et al. 2023).The motion of accretion flow along the magnetospheric surface is constrained by the local direction of the magnetic field lines.Our calculations were performed under the assumption of opacity being dominated by non-magnetic Compton scattering.Deviations of the magnetic field geometry from the dipole one due to the accretion process (Lai 2014) can affect the dynamics of the magnetospheric accretion flow. 1 The calculation of magnetic field geometry perturbed by accretion requires the solution of the RMHD equations, which is beyond the scope of this paper but will be done in the future publications.
Our examination of the accretion process under low mass accretion rates-where the influence of the radiative force is negligibly small-reaffirms findings previously reported in the ).The mass accretion rate at the inner disc radius is constant and fixed at Ṁ = 10 19 g s −1 .The rotation axis is taken to be orthogonal to the accretion disc plane.Parameters: m = 1.4,Rm = 10 8 cm, P spin = 10 s. literature: (i) For sufficiently small spin periods and inclined magnetic dipoles, specific regions on the magnetospheric surface exhibit material entrapment due to the interplay between gravitational and centrifugal forces (see Fig. 4, see also Abolmasov & Biryukov 2020;Lyutikov 2023).(ii) the conditions required for the arising of the "propeller" effect are dependent on the geometrical thickness of accretion discs and inclination of a magnetic dipole with respect to the disc plane ( Çıkıntoglu & Ekşi 2022).
We have demonstrated that accretion flow over the NS magnetosphere is feasible but exhibits instability at super-Eddington mass accretion rates, particularly when the radiative force becomes influential in shaping the flow dynamics (see Fig. 5).This instability tends to be quasi-periodic, with a typical period closely aligning with the dynamical timescale at the magnetospheric radius.Interestingly, the inclination of the magnetic dipole concerning the disc plane (see Fig. 7) and radiation beaming toward the axis of the magnetic dipole (see Fig. 9) act as stabilizing factors.Notably, quasi-periodic variability transforms into broadband variability at high Fourier frequencies (see Fig. 8).An unstable mass accretion rate onto the NS surface manifests as pronounced fluctuations in the surface density of matter covering the NS magnetosphere (see Fig. ,10).Intriguingly, even under conditions of an inclined magnetic dipole and a relatively stable mass accretion rate onto the stellar surface, significant surface density fluctuations persist.These fluctuations atop a variable mass accretion rate at the NS surface can profoundly impact the pulse profile formation process, rendering pulsations unstable and challenging to detect.However, it's essential to note that this assumption requires validation through further numerical simulations.
Appearance of high-frequency QPOs and broadband variability at high mass accretion rates was reported in bright X-ray transient GRO J1744-28, where the QPO frequency is fQPO ≃ 40 Hz (see, e.g., Mönkkönen et al. 2019).If this QPO is associated with the instability of magnetospheric accretion flow, the inner disc radius is expected to be ∼ 2×10 7 cm.The appearance of high-frequency component (f > 1 Hz) of PDS at high mass accretion rates was reported in the other Galactic bright X-ray transient -Swift J0243.6+6124(Doroshenko et al. 2020).
Possibly strong fluctuations of the mass accretion rate can affect the geometry of the emitting region in close proximity to a NS surface: the height of accretion columns (Basko & Sunyaev 1976;Lyubarskii & Syunyaev 1988;Mushtukov et al. 2015) is expected to be dependent on the mass accretion rate.Moreover, according to some theoretical models, variation of the mass accretion rate lead to fluctuating column height on the way to its stable state (Abolmasov & Lipunova 2022).Fluctuating geometry of the emitting region influences the beam pattern (Gnedin & Sunyaev 1973;Kraus 2001;Mushtukov et al. 2018) and formation of pulsation even without accounting for the influence of the magnetospheric accretion flow.
Detection of high-frequency QPOs can be difficult in the case of systems with a strong radiation driven outflows from the disc.The outflows should collimate X-ray photons from the central object.Due to the collimation process, the photons might experience a number of reflections/reprocessings by the walls of accretion cavity.As a result, one would expect suppression of any variability on time scales smaller than the typical time scale of photon travel inside the cavity (Mushtukov et al. 2021).The mass accretion rate onto NS surface (red curve) calculated for the case of fluctuating mass accretion rate at the inner disc radius (grey curve).The average mass accretion rate at Rm is Ṁ = 7 × 10 17 g s −1 (given by horizontal black dashed line), rms i = 0.05.One can see that fluctuations of the mass accretion rate at the NS surface are sensitive to the mass accretion rate at the inner disc radius.It is especially a case when the average mass accretion rate is close to the Eddington limit (as it is here).The rotation axis is taken to be orthogonal to the accretion disc plane.Parameters: m = 1.4,Rm = 10 8 cm, P spin = 10 s, isotropic central source.

Figure 1 .
Figure1.Schematic illustration of accretion geometry onto strongly magnetised NS.Magnetic field is considered to be dominated by the dipole component, while the rotation axis is assumed to be orthogonal to accretion disc plane.

Figure 2 .
Figure 2. The schematic geometry of dipole magnetic field line, where Rmax is the linear scale describing magnetic field line, λ (θ d ) is the latitude (co-latitude) of a point, and χ is the angle between the position point and the tangent line to the field line (9).

Figure 3 .
Figure 3.The colour maps represent acceleration along dipole magnetic field lines due to the gravitational and centrifugal forces at the magnetospheric surface.The red (blue) colour corresponds to the case of the acceleration directed towards the larger (smaller) latitudes, i.e. towards the upper (lower) pole of a star.The solid black line represents coordinates where the disc plane crosses the dipole surface.The magnetic obliquity is taken to be α = 0 (the first line), 30 • (the second line) and 60 • (the third line).Left, middle and right columns correspond to different NS spin periods: P spin = 10 s, 0.5 s and 0.1 s respectively.Parameters: m = 1.4,Rm = 10 8 cm.

Figure 4 .
Figure4.The maps illustrating the distribution of accretion flow surface density over the magnetospheric surface.Different maps correspond to different NS spin periods (P = 10 s, 0.5 s and 0.1 s) and various magnetic obliquity α.The inclination of the magnetic axis with respect to the rotation axis of a NS results in a fractional covering of the magnetosphere.In panels corresponding to α = 0 • and P = 0.5 s, α = 0 • and P = 0.1 s, α = 30 • and P = 0.1 s, and α = 60 • and P = 0.1 s the accretion flow does not reach the NS surface due to the centrifugal force and we see the propeller effect in action.In the case of large obliquity, sufficiently fast rotation of a NS results in the appearance of regions where material tends to accumulate but does not move towards a NS.Parameters: Ṁm = 10 17 g s −1 , m = 1.4,Rm = 10 8 cm.

Figure 6 .
Figure6.The modeled PDS of the mass accretion rate fluctuations at the NS surface.Different panels show the PDS for different mass accretion rates from the disc: 10 19 g s −1 , 5 × 10 18 g s −1 , 2 × 10 18 g s −1 , and 10 18 g s −1 (from top to bottom).Parameters: Rm = 10 8 cm, P spin = 10 s, isotropic source of emission.

Figure 9 .
Figure9.The mass accretion rate onto NS surface calculated for different beam patterns: the red line illustrates the case of isotropic central source, while the black line corresponds to the case of beamed emission from the poles of a NS (parameter asp = 2 in 19).The mass accretion rate at the inner disc radius is constant and fixed at Ṁ = 10 19 g s −1 .The rotation axis is taken to be orthogonal to the accretion disc plane.Parameters: m = 1.4,Rm = 10 8 cm, P spin = 10 s.

Figure 10 .
Figure10.Maps of the relative standard deviation of the surface density (see equaions 35 and 36) on the magnetospheric surface for the mass accretion rate fixed at Ṁm = 10 19 g s −1 and different magnetic obliquity α = 0 • , 15 • and 30 • (from top to bottom).In the case of aligned rotation and magnetic axis, the relative variability of the surface density is stronger at larger latitudes.In the case of inclined magnetic dipole, only a fraction of magnetospheric surface is covered by accretion flow and the variability becomes dependent on the azimuthal angle: the regions located closer to the edges of accretion curtain tend to be more variable.Parameters: m = 1.4,Rm = 10 8 cm, P = 10 s.
Figure11.The mass accretion rate onto NS surface (red curve) calculated for the case of fluctuating mass accretion rate at the inner disc radius (grey curve).The average mass accretion rate at Rm is Ṁ = 7 × 10 17 g s −1 (given by horizontal black dashed line), rms i = 0.05.One can see that fluctuations of the mass accretion rate at the NS surface are sensitive to the mass accretion rate at the inner disc radius.It is especially a case when the average mass accretion rate is close to the Eddington limit (as it is here).The rotation axis is taken to be orthogonal to the accretion disc plane.Parameters: m = 1.4,Rm = 10 8 cm, P spin = 10 s, isotropic central source.

Figure 12 .
Figure12.The red (blue) curve at the upper (lower) panel shows the mass accretion rate at the NS surface calculated for the case of constant (fluctuating) mass accretion rate at the inner disc radius.The latter is shown with a black line at the upper (lower) panel.The average mass accretion rate at Rm is 10 19 g s −1 .In the lower panel, the rms of mass accretion rate fluctuations at Rm is taken to be 0.1.Parameters: m = 1.4,Rm = 10 8 cm, P spin = 10 s, isotropic central source, α = 0 • .