Issue 
A&A
Volume 638, June 2020



Article Number  A142  
Number of page(s)  21  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/201936958  
Published online  26 June 2020 
Kilohertz quasiperiodic oscillations from neutron star spreading layers
^{1}
Department of Physics and Astronomy, 20014 University of Turku, Finland
email: pavel.abolmasov@gmail.com
^{2}
Sternberg Astronomical Institute, Moscow State University, Universitetsky pr. 13, 119234 Moscow, Russia
^{3}
Physics Department and Columbia Astrophysics Laboratory, Columbia University, 538 West 120th Street, New York, NY 10027, USA
^{4}
Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA
^{5}
Nordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, 10691 Stockholm, Sweden
^{6}
Space Research Institute of the Russian Academy of Sciences, Profsoyuznaya str. 84/32, 117997 Moscow, Russia
Received:
20
October
2019
Accepted:
4
May
2020
When the accretion disc around a weakly magnetised neutron star (NS) meets the stellar surface, it should brake down to match the rotation of the NS, forming a boundary layer. As the mechanisms potentially responsible for this braking are apparently inefficient, it is reasonable to consider this layer as a spreading layer (SL) with negligible radial extent and structure. We perform hydrodynamical 2D spectral simulations of an SL, considering the disc as a source of matter and angular momentum. Interaction of new, rapidly rotating matter with the preexisting, relatively slow material corotating with the star leads to instabilities capable of transferring angular momentum and creating variability on dynamical timescales. For small accretion rates, we find that the SL is unstable for heating instability that disrupts the initial latitudinal symmetry and produces large deviations between the two hemispheres. This instability also results in breaking of the axial symmetry as coherent flow structures are formed and escape from the SL intermittently. At enhanced accretion rates, the SL is prone to shearing instability and acts as a source of oblique waves that propagate towards the poles, leading to patterns that again break the axial symmetry. We compute artificial light curves of an SL viewed at different inclination angles. Most of the simulated light curves show oscillations at frequencies close to 1 kHz. We interpret these oscillations as inertial modes excited by shear instabilities near the boundary of the SL. Their frequencies, dependence on flux, and amplitude variations can explain the highfrequency pair quasiperiodic oscillations observed in many lowmass Xray binaries.
Key words: accretion / accretion disks / stars: neutron / methods: numerical
© ESO 2020
1. Introduction
For a nonmagnetised accretor with a fluid or solid surface, nonrelativistic disc accretion releases only about half of the gravitational binding energy. The other half is stored as kinetic energy of the flow (see e.g. Sibgatullin & Sunyaev 2000). As the accretor is unlikely to rotate close to its breakup limit, a greater proportion of this energy will still dissipate close to the surface of the accretor in what is called a boundary layer (BL). There is no commonly accepted view of a BL; even its basic geometry is uncertain.
Boundary layers are expected to appear during disc accretion onto stars and compact objects with relatively weak magnetic fields that are incapable of creating a magnetosphere. For neutron stars (NS), this happens for a surface magnetic field B ≲ 10^{8} G if the accretion rate is approximately Eddington. Lower mass accretion rates set a stronger limit for the magnetic field that is proportional to the square root of the accretion rate (see e.g. Lamb et al. 1973). This case is relevant for old NSs in lowmass Xray binaries (LMXBs). In particular, the BL apparently plays an important role in the socalled Z and atoll sources (as classified by Hasinger & van der Klis 1989). Combined spectral and timing observations of LMXBs allow the contribution of the BL to be separated from that of the accretion disc (Gilfanov et al. 2003; Revnivtsev & Gilfanov 2006). Emission of the boundary layer is hotter, with a colour temperature of about 2.5 keV. The BL spectral component is more variable than the disc on short, dynamical timescales. In particular, the highestfrequency, kilohertz quasiperiodic oscillations (kHz QPOs; see e.g. van der Klis 2000) can be interpreted as some type of BL activity. This is supported by the fact that, while the properties of all the lowfrequency QPO types are quite similar in NS and blackhole LMXBs, kHz QPOs behave in a profoundly different way for NS sources (Motta et al. 2017). It appears natural to attribute this difference to the existence of a solid surface and a BL.
These QPOs are often observed in pairs, and the distance between peaks is either close to the frequency of burst oscillations (itself well consistent with the spin frequency of the NS; see van der Klis 2000) or to one half of this value (Méndez et al. 1998; Wijnands et al. 2003). As the observational data accumulate, the picture becomes more complicated, favouring rather a universal, sometimes variable frequency difference of Δf ∼ 300 Hz (Méndez & Belloni 2007), close to but not equal or proportional to the spin frequency. Normally, there are only two kHz peaks present. One explanation is a bright spot rotating at a frequency of the order of Keplerian frequency (for a conventional NS with a mass of 1.5 M_{⊙}, and radius 12 km, linear Keplerian frequency is f_{K} ≃ 1.5 kHz) and producing one peak due to visibility effects and the other due to interaction with nonaxisymmetric structures at the surface of the NS. This interpretation is known as a beatfrequency model (first apparently proposed by Lamb et al. 1985 in the context of a different QPO type) and implies strict equality between the difference in the frequencies of the two QPOs and the spin frequency of the NS. Also, different types of beatfrequency models that do not involve a BL were proposed to explain the properties of kHz QPOs. The best known are the magnetospheric beatfrequency model (Psaltis et al. 1998), considering LMXBs as magnetised accretors with very compact magnetospheres, and the sonicpoint beat frequency model (Miller et al. 1998). The latter relies on the strong gravity effects that for a conventional NS place the last stable Keplerian orbit at a radius somewhat larger than the radius of the star. Both models require some mechanism generating narrowband variability at the local Keplerian frequency. Resonances between Keplerian, epicyclic, and vertical epicyclic frequencies emerging in general relativistic solutions are apparently a good explanation for the kHz QPOs in black hole systems (Kluzniak & Abramowicz 2001; Kluźniak et al. 2004), but predict a fixed ratio of 2:3 for the peak frequencies. In NS systems, the frequency ratio varies in rather broad limits, and 2:3 is only a crude approximation. The caveats of the conventional resonance model were considered by Rebusco (2008).
Quasiperiodic oscillation frequencies shift with time, remaining correlated with the observed flux on short timescales (hours and less) and uncorrelated on longer timescales. This creates parallel tracks in the flux versus QPOfrequency plot (Méndez et al. 1999; van der Klis 2001), suggesting that a BL possesses a characteristic correlation timescales much longer than even the viscous timescales in the disc (on which the mass accretion rate varies). It is difficult to suggest a way in which the accretion disc could produce the parallel tracks, as its variability is governed essentially by a single parameter: the mass accretion rate. If the BL is weakly coupled with the surface of the star, it is a good candidate for such an “integrator”. The rich phenomenology of kHz QPOs is a potentially important source of information about the NS itself and the physics of the flows close to its surface. However, besides the numerous observational clues and the variety of existing models, little is known about the mechanisms and exact relations between the quantities.
The classical approach to the BL considers the flow as some part of the disc (Pringle 1977; Popham & Narayan 1995) where the rotation velocity deviates strongly from Keplerian rotation and matches the rotation rate of the star at the inner edge. Strongly nonKeplerian rotation means that the approximations used as the basis for the thin disc approach are no longer valid, and the radial structure of the flow is strongly affected by the radial pressure gradient. Another problem is the efficiency of angular momentum exchange in the BL. In a hot (ionized) accretion disc with a rotation law close to Keplerian, there is an outward angular momentum transfer provided by the magnetohydrodynamic (MHD) turbulence generated by magnetorotational instability (MRI; introduced in the astrophysical context by Balbus & Hawley 1991), which is only operational when the angular velocity decreases with the cylindrical radial coordinate. For a BL, the rotation profile does not in general fulfil the necessary condition for MRI. It is unclear which physical mechanisms are responsible for the angular momentum exchange between the accreted matter and the surface of the star. In practically any possible BL model, the Rayleigh stability criterion (Rincon et al. 2007) is satisfied. Hydrodynamic turbulence is still produced for large enough Reynolds numbers (see Zhuravlev & Razdoburdin 2018 for a detailed analysis), but the amplitudes of turbulent motions are apparently insufficient for efficient angular momentum transfer. The very existence of the extremely long correlation timescale mentioned above suggests that irrespective of the mechanism of angular momentum transfer, it is extremely slow and inefficient.
A good candidate for such a mechanism is supersonic shear instability at the interface between the star and the BL (Belyaev et al. 2013; Hertfelder & Kley 2015; Philippov et al. 2016; Belyaev & Quataert 2018). Unlike the classical subsonic KelvinHelmholtz instability, oblique waves rather than vortices are generated. Moving in a shear velocity flow, the waves create Reynolds stress and thus provide effective nonlocal viscosity not only in the BL but also in the accretion disc. Most of the numerical studies mentioned above considered a twodimensional problem either in the equatorial plane or at a fixed latitude on a conical surface, as did Philippov et al. (2016). Belyaev & Quataert (2018) considered a threedimensional local MHD problem with a fixed equation of state, mainly addressing the structure of the flow in the equatorial plane. Depending on the particular simulation setup, the contribution of the wavemediated angular momentum transfer varies from negligibly small to somewhat comparable to that of the MHD turbulence generated in the disc.
If the angular momentum exchange rate between the accreting matter and the material of the star is smaller than the angular momentum supply from the disc, rapidly rotating matter would accumulate on the surface, pushed to higher latitudes by pressure gradient. The radial dimension of such a flow, as well as of a conventional BL, is second order in relative disc thickness, much smaller than its vertical (along the polar angle) extent (Papaloizou & Stanley 1986). Consequently, one can treat the flow as twodimensional (2D) on the surface of the NS fed by matter and angular momentum injection from the disc. This approach, known as the spreading layer (SL) approach, was introduced by Inogamov & Sunyaev (1999) and further developed in Suleimanov & Poutanen (2006) and Inogamov & Sunyaev (2010). For the case of accreting white dwarfs, this model was considered by Piro & Bildsten (2004a) and used to explain a certain type of QPO observed in cataclysmic variables, the socalled dwarfnova oscillations (DNOs; Piro & Bildsten 2004b). The angular momentum transfer within the SL depends on the dynamics of the flow itself as well as existing oscillation and instability modes. It is quite probable that certain hydrodynamical phenomena will provide an efficient way to transfer momentum within the SL and thus define its internal dynamics. In this paper, we try to address the issue of angular momentum transfer within the layer using numerical hydrodynamical simulations.
Twodimensional hydrodynamics on a rotating sphere is an important subject in geophysics and astrophysics, as many spherical bodies, including planets, stars, and NSs, have fluid atmospheres. Vertically integrated equations of hydrodynamics lead to the system of shallow water equations (see e.g. Vreugdenhil 1990), normally used in geophysics for weather forecasts in combination with spectral methods (JakobChien et al. 1995). Spectral methods provide much higher accuracy than finitedifference methods on an equally fine grid, and are quite robust and stable for subsonic flows. However, rotation in a SL is almost always supersonic, which makes the flow compressible and its simulations potentially prone to numerical instabilities. This makes numerical simulations of spreading layers technically challenging. On the other hand, using spherical harmonics is natural in spherical coordinates and allows the singularity of the spherical grid near its axis to be avoided. We provide our full simulation code SLAYER as opensource software ^{1}.
The paper is organised as follows. In Sect. 2, we formulate the physical problem and derive all the basic equations. Results of the simulations are given in Sect. 3. Applications and limitations of the model are discussed in Sect. 4. We conclude in Sect. 5. A detailed description of the numerical techniques is given in Appendix B.
2. Physical model
2.1. Scales and dimensionless quantities
Natural timescale in the vicinity of a relativistic object of mass M is
which is the approximate lightcrossing or dynamical timescale at the event horizon. The corresponding radius is
The radius of the NS is taken as R_{*} ≃ 12 km ≃ 6R_{g} assuming a mass of 1.4 M_{⊙} (see e.g. the estimates of masses and radii in Miller & Lamb 2016; Nättilä et al. 2017).
All the geometrical and kinematic quantities are naturally normalised by combinations of these spatial and temporal scales. In particular, velocities in units of the speed of light c are used. We use physical quantities (g cm^{−2}) for Σ and internally normalise the surface density by an arbitrary scale set either to 10^{4} or to 10^{8} g cm^{−2} for the simulations presented in this paper. Evidence for a very long correlation timescale t_{corr} suggests a characteristic value of
The physical meaning of this value is the mean surface density if t_{corr} is a characteristic time of mass depletion or replenishment in the SL. The vertical optical depth of the layer is simply 𝜘Σ, where 𝜘 is opacity. Here, we set 𝜘 to Thomson electron scattering opacity for Solar metallicity, 𝜘_{T} ≃ 0.34 cm^{2} g^{−1}.
The problem has a complex hierarchy of timescales, starting with the dynamical scale, which is normally the smallest. The characteristic dynamical timescale is the Keplerian period near the surface of the NS:
The local thermal timescale depends on the effective and internal temperatures. For LMXBs, in Z and brighter atoll states, mass accretion rates span the range 10^{15} − 10^{18} g s^{−1} (10^{−11}–10^{−8} M_{⊙} yr^{−1}) which implies effective temperatures of the order
As half of the accretion power is emitted by the disc, we use 8π in the denominator for this estimate. Internal temperature, T_{in}, is about a factor (𝜘Σ)^{1/4} ∼ 100 times larger if Σ ∼ 10^{8} g cm^{−2} is taken as a representative value. The thermal timescale is then easy to estimate, separating contributions of the radiation and gas pressure. If β is the gas pressure fraction,
where Q^{−} represents the energy lost via radiation (introduced more rigorously in Sect. 2.5). Dependence on the gas pressure fraction follows from the relations between the vertically integrated and local quantities derived below in Sect. 2.3. As the gas becomes radiationpressure dominated, β approaches zero, and the thermal timescale becomes longer. As we show below (Eq. (15)), the gas pressure ratio is related to the Eddington factor 𝜘Q^{−}/cg_{eff} = 1 − β. Violation of the local Eddington limit in our model is impossible as it makes the effective gravity negative, and the material of the layer unbound. On the other hand, the amount of internal energy stored in the layer is essentially unlimited (E ∝ β^{−1}). However, for very small β, the vertical thickness of the layer (see Eq. (21)) becomes comparable to R_{*} and thus limits the thermal energy stored in the SL. For T_{in} = 100 keV, the minimal possible gas pressure fraction is about β_{min} ≃ 10^{−4}, which corresponds to a thermal timescale longer than a day, meaning that very rapid accretion effectively runs in a radiatively inefficient regime. If energy release is close to equilibrium with radiation losses, the effective temperature does not change significantly, and most of the variations of the thermal timescale are related to the energy stored by gas and trapped radiation. Thermal and dynamical timescales become comparable if the surface density is small, Σ ≲ 10^{3} g cm^{−2}.
The challenge of SL simulations for the case of LMXBs is in the relatively long accretion timescale. While phenomena like kHz QPOs manifest themselves on dynamical timescales of milliseconds, the accretion rate is relatively stable at the scales of minutes to hours (viscous times of the inner disc), and the putative time of mass growth and depletion in a SL is apparently of the same order. However, these timescales are much longer than both the thermal and dynamical timescales.
Hence, all the dynamical and thermaltimescale phenomena we intend to consider appear in fact in a quasistationary SL where mass accreted during the considered timescales is negligibly small. However, to reach such an equilibrium state, one needs to simulate either an episode of much more rapid accretion or, alternatively, accretion atop a much thinner atmosphere, where the surface densities of preexisting and newly accreted material would be comparable on a reasonable timescale of the simulation run. We try both approaches.
2.2. Geometry and mass conservation
Continuity equation for density ρ and velocity v in the most general Newtonian form is
Integration over r yields
where S^{±} are the source and sink terms for surface density. The source term is set explicitly as an inclined belt (to account for the case where the disc is inclined with respect to the rotation axis of the NS),
where α is the angular distance from the direction of the adopted disc rotation axis,
where disc inclination i with respect to the spin axis of the star sets the direction of the symmetry and rotation axis of the source (see Fig. 1). Here, we use a spherical coordinate system consisting of radial coordinate r, colatitude θ, and longitude φ. Instead of θ, latitude λ = π/2 − θ is used below for visualisation. The system is aligned with the rotation of the star, but is itself nonrotating (inertial). The SL is considered thin in the radial direction, making r, in a sense, a vertical coordinate, along which hydrostatic equilibrium is assumed to hold. See Sect. 2.3 for more details.
Fig. 1.
Illustration of the model geometry. The tilted blue arc near the equator shows the source of mass and momentum. The spin axis of the star, marked with Ω_{*}, is inclined with respect to the disc axis (Ω_{d}), by an angle i. 
Adding a sink allows to limit the growth of the mass of the SL and approach a quasistationary state. Matter existing long enough on the surface of the NS should adopt its velocity and physical properties. As our model adopts a simplified vertical structure, adding a sink allows to draw an effective boundary between the SL and the material of the NS. In this paper, we ignore the sink term, but include it in the equations for future use in the form
where t_{depl} is a depletion timescale set explicitly. This form is valuable for its simplicity and allows us to study the role of the mass of the SL by considering a gradual depletion regime (mass accretion is switched off but the sink is not) and the steadystate spreading regime by turning simultaneously on both the source and the sink.
We also use a tracer quantity a to separate the contributions of the preexisting and newly accreted material. It is assumed to be a passive scalar quantity transferred with the flow and initially equal to zero. The source of this quantity is designed in such a way as to reproduce the evolving fraction of newly accreted matter,
2.3. Vertical structure
In the vertical (radial) direction, SL is supported by thermal (gas and radiation) pressure together with the relevant centrifugal force component. As we consider the vertical extent of the SL to be infinitely small, the timescale of vertical dynamical relaxation should also be small, and therefore we neglect the effect of radial velocity in momentum and energy equations. This allows us to write down a hydrostatic balance equation
which needs to be supplemented by another equation to calculate the vertical profiles of pressure and density simultaneously. Let us assume, following Inogamov & Sunyaev (1999) and Suleimanov & Poutanen (2006), that the heat is released at the bottom of the SL, and the optical depth is high enough to use the radiation diffusion approximation. Hence, radiation flux F is constant with height
where p_{rad} is radiation pressure. Total pressure is assumed to be contributed by p_{rad} and gas pressure p_{gas}. Together, Eqs. (13) and (14) imply a constant pressure ratio β = p_{gas}/p as long as opacity is constant with height. Hence, gas, radiation, and total pressure scale with each other, and the gastototal pressure ratio equals
Proportionality of pressures also implies p ∝ ρT ∝ T^{4}, which leads to p ∝ ρ^{4/3}, an effectively polytropic law. This implies a relation between the pressure p_{0} and density ρ_{0} at the bottom of the SL and the corresponding vertically integrated quantities, pressure Π = ∫pdr and surface density Σ = ∫ρdr,
Constancy of β allows also to link surface energy density E = ∫εdr (where ε is the volumetric energy density) and integrated pressure Π with each other. Local energy density may be expressed as
which implies an identical relation for the vertically integrated quantities
The pressure ratio itself may be found as a function of E (or Π) and Σ. At the bottom of the layer,
where m ≃ 0.6m_{p} is the mean mass of a massive particle, which allows us to solve implicitly for β, taking into account the expression for p_{0} = g_{eff}Σ arising as a solution to Eq. (13) and Eq. (18) and substituting them into Eq. (19):
The thickness of the layer may be found as the solution of hydrostatic equation, taking into account the constancy of β (see also Eq. (32) in Suleimanov & Poutanen 2006) as
Strictly speaking, the adopted dissipation at the bottom of the layer is a very simplified picture, as the energy dissipation is in general distributed along the vertical coordinate in a way that is dependent on the unknown mechanisms involved. As long as diffusion approximation is valid, vertical distribution of the dissipation processes leads to two systematic effects: a decrease in the effective optical depth, and deviation from the effective vertical polytrope used in this section. Both are likely to introduce correction factors of the order unity, which would be easy to calculate once our model is extended to three dimensions with a reasonable set of boundary conditions on the surface of the NS.
2.4. Momentum equations
We start with Euler equations with additional source and sink terms related to the momentum of the matter being accreted and to the friction between the SL and NS surface. Their general vector form is
where g is gravity without the contribution of centrifugal force and is assumed to be directed along the radius vector. At the same time, the surface of the NS is close to being equipotential and thus is deformed due to rotation (we neglect all the other sources of deformation, such as magnetic fields and nonequilibrium stresses in the crust), which makes the polarangle component , where Φ is gravitational potential, nonzero even after vertical integration.
The radial component of the momentum equation reduces to the hydrostatic equation considered in Sect. 2.3. The two tangential components of the equation are convenient to rewrite in terms of the two scalar quantities normally used in shallowwater approximation: vorticity,
and divergence,
Multiplying Eq. (22) by ρ and integrating over the total vertical extent of the SL yields
A detailed derivation of the equations for vorticity and divergence is given in Appendix A.
Taking the radial curl component of Eq. (25) results in an equation for vorticity
where Ω_{*} is the angular frequency of the NS. The velocity field of the accreting matter v_{d} is assumed to be uniform rotation with the angular frequency c_{K}Ω_{K}, where c_{K} is the deviation from the Keplerian rotation law in the disc, set to 0.9 in all the simulations. Vorticity of this velocity field is ω_{d} = 2c_{K}Ω_{K} cos α, where α is given by Eq. (10). The last term in (26) describes viscous coupling between the SL and the surface of the NS. Our lack of knowledge about the nature and strength of this coupling is included in the unknown friction timescale t_{fric}. Preceding terms containing S^{+} appear due to accretion of matter with a given vorticity. In addition to this, the righthand side of the equation includes a baroclinic term equal to zero if a fixed equation of state is adopted, or if the distributions of pressure and density are exactly axisymmetric. This term can create vorticity through entropy variations.
Taking divergence of Eq. (25) provides an equation for δ
The term originates from the rotational deformation of the NS.
2.5. Energy conservation
In general form, energy conservation implies (Suleimanov & Poutanen 2006):
where the righthand side accounts for heat exchange with the NS (q_{NS}), heat released within the layer (q^{+}), and radiation losses q^{−}. After integration, all the q quantities result in corresponding capital Q quantities: fluxes through the surface and energy release per unit area.
We treat the hydrodynamics of the SL as ideal, though the numerical solution techniques used (described later in Appendix B) provide dissipation on small scales close to the spatial resolution used. If momentum transfer is dominated by turbulent motions forming a direct cascade similar to Kolmogorov cascade (Monin et al. 2007) where energy is transferred from larger to smaller scales, the exact nature and properties of the viscous dissipation at the small scales are irrelevant. However, conservation of energy implies that viscous dissipation should act as an additional source of internal energy. All the kinetic energy lost by the flow should reappear as heat. We assume this heat to appear at the bottom of the flow and to diffuse upwards leaving the SL from its upper surface. As already mentioned in Sect. 2.3, for dissipation taking place throughout the volume, this introduces a systematic uncertainty of the order unity in the equations of vertical structure.
Taking into account momentum conservation, friction, and viscous dissipation, and integrating the energy equation vertically, we end up with the equation
where Q^{+} is the heat released in the spreading layer, Q_{NS} is the heat received from the neutron star, Q^{−} is the radiation flux lost from the surface, and the additional term Q_{acc} corresponds to the thermal energy introduced with the accreting matter and released during its mixing with the preexisting material. Vertically integrated pressure and energy are related by (18). Dissipation is calculated as kinetic energy lost by the flow:
where the dissipation in velocity is related both to the friction term in Eqs. (26) and (27) and with numerical dissipation on the grid scale, which we discuss in detail in Appendix B. The energy radiated away from the surface is set by radiation energy diffusion (see Sect. 2.3):
where 𝜘 is the Rosseland average opacity which we assume to be equal to the Thomson scattering opacity, 𝜘 = 𝜘_{T} ≃ 0.34 cm^{2} g^{−1}. If β approaches zero, the SL becomes a “levitating layer” supported mainly by radiation pressure (Inogamov & Sunyaev 1999). For β ≪ 1, the energy loss term is nearly independent of the physical conditions inside the layer.
Additional terms related to the accreting matter are more difficult to constrain from a physical point of view. It is natural to assume that the initial temperature of the newly introduced material is nonzero, and hence the increase in surface density is accompanied with an increase in surface energy density as well. In addition, as the velocities of the accreting and the preexisting matter are in general different, some of the kinetic energy is dissipated during the process of mixing (the exact physical mechanism could be kinetic, hydrodynamic, or MHD). Energy and momentum conservation laws predict that the amount of dissipated energy per unit of accreted mass is (v_{d} − v)^{2}/2, hence
where again the “d” index corresponds to the properties of the mass source.
Our approach to the vertical structure, as well as to momentum and energy conservation, is mostly in line with the works of Inogamov & Sunyaev (1999, 2010), and Suleimanov & Poutanen (2006). However, momentum and energy equations introduced below assume neither axisymmetry nor stationarity, both of which were crucial for the quasistationary, onedimensional SL model. Relaxing these assumptions requires that we specify certain physical mechanisms. In particular, the stationary nature of the flow allows us to relate surface density and latitudinal velocity in algebraic form independently of the mechanism of angular momentum loss. In addition, Inogamov & Sunyaev (1999) do not consider rotational deformation of the NS, meaning that a corotating atmosphere in their assumptions should strongly concentrate near the equator. However, we take into account small equilibrium deformation of the star. Apart from making the simulations more realistic, this equilibrium deformation allows a simple set of initial conditions with uniform surface density and pressure.
2.6. Initial conditions
The simplest possible initial conditions are constant surface density and pressure in combination with rigidbody rotation. This may be achieved if the rotation rate is exactly equal to the rotation frequency of the NS, Ω_{*}, and the deformation of the NS makes its surface equipotential. Configuration is stable and may survive for simulation times vastly exceeding the durations of the runs used in this work. In Appendix B.2.1, we use this initial condition configuration to check the numerical stability and dissipation of our numerical scheme.
Vorticity of a rigidbody rotation is
As the motions are limited to pure rotation, initial divergence is strictly zero. To the basic initial condition set, a small (5%) perturbation was added in the form of an over or underdensity. The perturbation is designed as an entropy variation not affecting the pressure distribution.
3. Results
3.1. Model setup
The equations derived in Sect. 2, including the equations of mass (Eq. (8)), momentum (Eqs. (26) and (27)), and energy (Eq. (29)) conservation, were solved using our 2D spectral modelling code. We refer to Appendix B for a detailed description of the numerical techniques. There, we also describe the tests for numeric performance, stability, and accuracy. We list all the SL models calculated for this paper in Table 1. Letters “LR” and “HR” in a simulation ID always refer to “low” (128×256) or “high” (256 × 512) resolution. Consistency between the corresponding low and highresolution runs is an important test for numerical effects (noise and diffusion). Below, we describe the setups of all these models, while the description of the results is given in the following sections.
Spreading layer simulations.
All the models include a NS rotating with a spin period P_{spin} = 3 ms. Initially, all the matter on the surface rotates together with the star as a rigid body. As we so far do not include any friction with the NS surface, rotation of the star affects the results only through the initial conditions and the deformation of the stellar surface (potential term in Eq. (27)).
Most models include a mass source corresponding to a steadystate accretion from a thin disc. To avoid strong velocity gradients causing highfrequency noise, the mass accretion rate approaches the steadystate value smoothly, following the exponential law
where the turnon timescale was set to 10P_{spin} for all the simulations. The spatial distribution of the source of mass is always a Gaussian function of cos α as given by Eq. (9) with the standard deviation of cos α_{0} = 0.1, corresponding to the width of about 6°. Rotation rate of the newly accreted material is set to c_{K} = 0.9 in Keplerian units.
As the timescales observed in real LMXBs differ by six orders of magnitude, it is difficult to perform a single realistic simulation resolving the dynamictimescale variability over a timescale that is sufficient to see the changes in the SL structure. We use two approaches to avoid this difficulty: first, we consider “enhanced accretion” with Ṁ = 10^{−3} M_{⊙} yr^{−1} ≃ 6 × 10^{20} g s^{−1} (models 3LR and 3HR); secondly, we consider accretion on top of a thin layer with Σ_{0} = 10^{4} g cm^{−2}. Both tricks shorten the effective evolution timescales by several orders of magnitude. We expect the hydrodynamics to be equally effective in both configurations even though the radiation timescales are vastly different.
As an alternative to both these approaches, we calculate a model reproducing the evolution of a spreading layer after switching off the mass source (model 3LRoff). It starts with the final snapshot of 3LR and then gradually cools down for another 0.5 s. Finally, in the model 3LRinc, we consider the case where the source is inclined with respect to the initial rotation plane.
3.2. General properties
Each simulation covers several tenths of a second, which is significantly longer than the dynamical timescales, and for the parameters of our simulations is comparable to the timescale on which the initial mass content of the layer is replaced by newly accreted matter. Accretion is concentrated (everywhere except 3LRinc) at low latitudes. As the surface mass and energy density increase near the equator, matter is pushed toward the poles. Surface density in the equatorial region tends to become larger due to accumulation of the accreted material, but often becomes lower, as the mixing between the streams moving at different velocities leads to high levels of dissipation (see Eq. (32)). Equatorial regions tend to spin faster in all the models (see Fig. 2), but the excess angular momentum is redistributed by oblique waves (for rapid accretion modes) or by the heating instability (for 8LH/8HR). The main difference between the models with different mass accretion rate is the role of cooling: for rapid accretion, heating occurs much faster and the layer effectively accumulates heat, while for the models 8LR/8HR, radiation efficiently cools down some portions of the flow.
Fig. 2.
Timelatitude diagrams for longitudinally averaged angular frequency normalised to the Keplerian rotation frequency at the radius of the star R_{*}. The upper and lower panels correspond to the rapid accretion model 3LR and 8LR, respectively. 
For the low mass accretion rate (model 8LR), the evolution of the SL is illustrated by Fig. 3. We note the decrease in surface density at latitudes of about ±20° (t≃50–100 ms) due to heating, the gradual northsouth symmetry breaking between t≃200 and 250 ms, and later dynamicaltimescale evolution. In these simulation runs, β ≃ 0.99–1, close to 1 with an accuracy of about 10^{−5} outside the regions of intense mixing and heating instability.
Fig. 3.
Timelatitude diagram for longitudinally averaged surface density (normalised by the surface density averaged over the sphere) in the realisticaccretionrate 8LR simulation. Upper panels: three snapshots of surface density (10^{4} g cm^{−2} units) during the later stages of the heating instability development. 
For the high mass accretion rate models, 3LR and 3HR, the main process driving the subsequent evolution is the velocity contrast between the new and old matter, leading to a shear instability. A sequence of snapshots of the radiation flux and velocity field for model 3HR is shown in Fig. 4. The gastototal pressure ratio β for enhanced accretion models varies between about 0.1 in the accretion region and ∼0.9 near the poles.
Fig. 4.
Radiation flux (in Eddington units, ) snapshots for the model 3HR. White streamlines show the velocity field and the black contour corresponds to the accretion tracer a = 0.5. For all but the first snapshot, newly accreted matter dominates between the black contours. 
As we do not include the sinks in this paper, no quasistationary picture is expected to be reached. However, heating and velocity gradients created by accretion lead to at least two important dynamical effects relevant for the dynamics of SLs. In the two groups of simulations, with ‘enhanced’ and ‘normal’ mass accretion rates, two different instabilities emerge. For the realistically low mass accretion rate (Ṁ = 10^{−8} M_{⊙} yr^{−1}) and low initial surface density (Σ_{0} = 10^{4} g cm^{−2}), the equatorial belt forming out of the accreting matter during the first milliseconds of accretion is cooled efficiently, and most of the subsequent dissipation takes place at higher latitudes where the newly added material mixes with the old, slowly rotating NS atmosphere. This results in a heating instability: local displacement of the cool equatorial belt material leads to increased dissipation on the opposite side of the equator, which results in a pressure gradient increasing the initial displacement. This is easily seen in Fig. 3, where the latertime evolution (starting at approximately t ∼ 0.2 s) is marked by a gradual and then dynamicaltimescale development of a strong mass asymmetry between the southern and the northern hemispheres. Development of the instability also breaks the axial symmetry, especially during the period of rapid evolution. As a large amount of matter migrates between the polar and equatorial regions, some part of the flow acquires very large, and some very slow rotational frequency (even smaller than the rotation velocity of the NS).
At the same time, for the case of enhanced accretion, cooling and heating timescales are much longer, and all the observed dynamical effects are purely hydrodynamical. They are reasonably reproduced by the lowresolution simulations, though increasing resolution reveals more details and adds regularity to the observed turbulent patterns (see Fig. 5). Spectral simulations are able to capture the oblique wave patterns and “curling” of the equatorial spreading layer belt even with the low numerical resolution runs. For the highresolution runs, more finescale substructure starts to appear.
Fig. 5.
Resolution effects after t = 0.03 s (ten spin periods) of evolution for the highaccretionrate models 3LR (low numerical resolution; left panels) and 3HR (high numerical resolution; right panels). Top: vorticity maps. Bottom: emitted radiation flux Q^{−} (normalised by the local Eddington value) map. 
3.3. Angular momentum transport within the layer
The highaccretionrate models 3LR and 3HR demonstrate a system of oblique waves generated by the velocity discontinuity. The discontinuity itself in our simulations is a consequence of the initial conditions. In real sources, the existence of old material corotating with the star is a consequence of angular momentum exchange with the star. Whether or not the velocity drop exists in a quasistationary case with a sink and friction is an open question that may be resolved by running a simulation with sink terms for a sufficiently long time. In the new, rapidly rotating part of the flow, a system of standing waves rapidly evolves into a nonlinear regime, and forms a nonaxisymmetric wiggle structure seen in Fig. 5 that is subsequently smeared off. At large latitudes, where the old, slowly rotating matter dominates, the waves create a correlation between orthogonal velocity components. Velocity correlation provides a Reynolds stress component
where angular brackets ⟨…⟩ denote averaging in time and longitude. In Fig. 6, we show the value of T_{θφ} calculated for the model 3LR for the period of time 0.1–0.3 s after the start. Apparently, Reynolds stress is small in comparison to pressure but surprisingly stable in its sign, showing a clear poleward angular momentum transfer. While the SL itself appears to approach a quasistationary axisymmetric state on a timescale of about 1 s, higher latitudes still show variability and nonaxisymmetry up to the end of the simulation. In real astrophysical sources, the mass accretion rate is variable on subsecond timescales, meaning that even if the Reynolds stress is a reaction to the variations in mass accretion rate, it should always be present.
Fig. 6.
Azimuthally and temporally averaged dynamical quantities for the highaccretionrate model 3LR. Visualised quantities are Reynolds stress (black), mean velocity product v_{θ}v_{φ} (blue), and sound velocity (red dotted) as functions of latitude. All the quantities were averaged over the period of time between 0.1 and 0.3 s, and over the azimuthal angle. Solid and dashed lines correspond to positive (southward motion or momentum transfer) and negative quantities (northward). Dotted black curves show the 1σ standard deviation interval for the Reynolds stress. 
The existence of a small but significant hydrodynamical stress suggests a long local viscous timescale corresponding to obliquewavemediated angular momentum transfer. Near the poles, Reynolds stress partially compensates the advective angular momentum flow caused by compression of the preexisting polar cap material. Near the equator, at the same time, both advection and viscous transport of angular momentum are directed polewards.
For the realisticaccretionrate models, 8LR and 8HR, heating instability creates a strong flow of mass toward one of the poles. At the later stage of this process, when most of the accumulated mass flips to one side, axial symmetry is broken, which effectively creates a very large Reynolds stress spreading the angular momentum of the rapidly rotating matter over latitudes (see Fig. 7). Reynolds stress rapidly removes the angular momentum from the equatorial stream in both directions (we note the sign change near −10° latitude).
Fig. 7.
Azimuthally and temporally averaged dynamical quantities for the realisticaccretionrate model 8LR. Symbols and quantities are the same as in Fig. 6. Averaging was performed from t = 0.2 to 0.28 s, when the initial hemisphere asymmetry is developed due to heating instability. 
3.4. Artificial light curves
To produce artificial light curves, we use a simplified approach that ignores all the relativistic effects. We choose an inclination of the observer i_{obs} and integrate the bolometric flux Q^{−} emitted from the surface:
where
and α_{obs} is the angle at which the surface element is seen by an observer located at the colatitude of i_{obs} and at the azimuthal angle of φ = 0 in the spherical coordinate system. Such an approach allows us to reproduce the effects of visibility of any moving features on the surface of the star. The factor four in Eq. (36) allows us to interpret L_{obs} as isotropic luminosity, equal to the actual luminosity for an isotropic source.
Power density spectra (PDS) were calculated using the standard FFT algorithm (Cooley & Tukey 1965) with fractional normalisation. If the light curve is set as a series of observed luminosities, L_{k}, at the equidistant instances of time, t_{k}, the Fourier power density (in Miyamoto normalisation, see e.g. Miyamoto et al. 1991; Nowak et al. 1999) is found as a function of frequency, f, as
The frequency grid on which the PDS is calculated is equally spaced with Δf = 1/T, where T is the time span. Spectral power defined this way is a measure of the relative amplitude of a variability mode. For a broad spectral peak, variability amplitude may be estimated as , where ΔPDS stands for the excess spectral power associated with the particular spectral detail, and summation is done over the relevant spectral interval.
In Fig. 8 we show the dynamic (calculated inside 20 separate time bins) powerdensity spectra for different i_{obs}. We plot the relative PDS multiplied by f^{2} to decrease the contamination from the lowfrequency noise component related to the overall shape of the simulated light curve. Several oscillation modes are visible, one of them for a poleon observer. Their frequencies evidently correlate with the flux. Most of the nonaxisymmetric structures in this simulation are moving slightly faster than the star itself: their contribution is visible just above Ω_{*}. There is also power at about double spin frequency and in the very beginning near the third harmonic. The perturbation seen during the first ∼0.1 s of the simulation is mostly related to the initial perturbation rotating at the spin frequency. However, very little variability is seen by a polar observer, except for a single peak initially close to 1.5Ω_{*} and then gradually increasing its frequency towards 700–800 Hz. This signal is visible for all the inclinations but in general is weaker than the nonaxisymmetric modes absent in the poleon dynamic spectra. The properties of this mode fit well into the concept of a latitudinally propagating surface wave moving in a waveguide, similar to the modes considered by Piro & Bildsten (2004b). We discuss the implications of such an interpretation later in Sect. 4.
Fig. 8.
Dynamical PDSs (log_{10}(f^{2}PDS(f)) normalised to total power within a single time bin) for the “isotropic luminosity” calculated as described in Sect. 3.4 for the highaccretionrate model 3LR. The three upper panels show dynamical spectra for the observer’s inclinations of π/2, π/4, and 0, respectively. White horizontal lines show the spin (lower) and Keplerian frequencies. Lowermost panel: corresponding light curves: i_{obs} = π/2, π/4, and 0 cases are shown with blue dotted, green dashed, and black solid lines. 
In Fig. 9 we show how the dynamic PDS changes when rapid accretion stops (model 3LRoff, that starts with the end of simulation 3LR). The quasiperiodic features around twothree spin frequencies retain at least for the period when the layer remains hot (before t ≃ 0.95 s). The axisymmetric mode is evidently split into two QPO features. A hint of such a split is visible also at t∼0.3–0.4 s in Fig. 8. All the characteristic features, as in the original simulation with accretion, correlate with flux.
For the inclined simulation, 3LRinc (Fig. 10), there is an early stage (0 − 0.2 s) of the collision between the two flows inclined to each other. Subsequently, a quasiaxisymmetric configuration forms, and the variability pattern becomes similar to those of the aligned models discussed above. Unlike the aligned case, the observed luminosity changes in very narrow limits, probably because the size of the SL is now determined by the geometry of the inflow rather than by angular momentum transfer. Dissipation is smoothly distributed over the whole latitude range between −i and i, and saturates at a level Q^{−} ≃ cg_{eff}/𝜘_{T}. The apparent luminosity seen by a poleon observer is
Starting from t ≃ 0.15 s, the poleon PDS shows one stable peak at about 1 kHz and a hint of another peak at about 1.5 kHz, sometimes split in two (see Fig. 14 showing the PDS integrated over the integral 0.4–0.65 s; a similar picture is seen for t ∼ 0.2 − 0.3 s). Unlike the aligned case, nonaxisymmetric modes at later stages appear slower than the axisymmetric. This is probably related to the overall change in angular momentum of the layer that is affected by the preexisting matter rotating in a different direction.
Peak frequencies in the dynamical PDSs are shown in Fig. 11 as functions of flux. There is an evident signal for the poleon simulated light curve both in the original model with enhanced accretion, and in the switchoff simulation. The axisymmetric mode discussed above dominates for the poleon case, for which a clear correlation between flux and frequency is observed. For an inclined observer, nonaxisymmetric modes are stronger. Interestingly, an inclined source is capable of exciting variability modes with frequencies lower than the spin frequency.
Fig. 11.
Peak frequency (calculated as the position of the maximum of f × PDS) as a function of observed luminosity (calculated using expression (36)) for i_{obs} = 0 (black circles), π/4 (red diamonds), and π/2 (blue triangles). Simulation runs 3LR (left) and 3LRoff (right). Error bars show flux dispersion within the time bin and the size of the frequency bin where the maximum was detected. Horizontal dotted green lines show spin frequency. 
We also consider PDSs integrated over time intervals where the shapes of dynamic PDSs remain relatively stable and/or show hints of additional spectral features. In Fig. 12, we show such a spectrum for the 3LR simulation, computed for t = 0.2–0.5 s. Poleon PDS is dominated by a single narrow peak at about 800 Hz. At large inclinations, this single QPO transforms into two, and an additional third peak emerges at about one spin frequency. For 3LRoff, splitting of the main peak for i_{obs} = 0 is visible at t ∼ 0.5–0.7 s (shown in Fig. 13) and later at t ∼ 0.8 s. Later, the structure of the layer starts rapidly changing due to rapid cooling and switching to the gaspressuredominated regime.
Fig. 12.
Integral PDS (time span 0.2 to 0.5 s) of the simulation 3LR for different inclinations. Black dots correspond to a poleon observer i_{obs} = 0, red diamonds to the intermediate inclination of i_{obs} = π/4, and green triangles to i_{obs} = π/2. 
For the realistic massaccretionrate configuration, oscillations appear simultaneously with the development of the heating instability, and their contribution is only clearly visible at particular inclinations (see Fig. 15). The two poles behave in a profoundly different way because of the asymmetry formed by heating instability. Further simulations of the 8LR case are difficult because of the very strong density contrasts formed during the instability.
Fig. 15.
Integral PDS (time span from 0.27 to 0.32 s) for the simulation 8LR for different inclinations (black circles i_{obs} = 0, red diamonds i_{obs} = 45°, green upward triangles i_{obs} = 90°, and blue downward triangles i_{obs} = 180°). 
4. Discussion
In the dynamical spectra presented in Sect. 3.4, there are clearly at least two types of quasiperiodic variability signals: one disappears for a poleon observer and is thus related to nonaxisymmetric structures (waves and vortices produced by shear instabilities); and the other is present at all inclinations. The frequency of this mode clearly increases with the flux, approximately as (see Fig. 11). We leave a detailed study of the properties of the predicted QPOs to a separate paper.
It is natural to interpret this oscillation mode as a surface mode existing within the SL, as was done by Piro & Bildsten (2004b) for dwarfnova oscillations (DNOs). However, the width of the SL, independently from the assumptions, should increase with mass accretion rate. This is especially true for the radiationpressuresupported case where the radiative flux is fixed by equilibrium with effective gravity F = cg_{eff}/𝜘, and the growth of the area over which the dissipation is spread should reflect the growth of the mass accretion rate. Approximately, the width of the SL grows linearly with the mass accretion rate. As the speed of sound depends weakly on the mass accretion rate, the frequency of a DNOlike sonic mode for a thin radiationpressuresupported SL is
where is the speed of sound. However, this approach assumes that the observed oscillations are produced near the equator. The equatorial belt is indeed responsible for most of the energy dissipation, but the radiation flux is broadly distributed over the surface due to the importance of radiation pressure, and most of the variability comes from high latitudes (see Fig. 16). This is visible in Fig. 16, where most of the dissipation in the system takes place directly at the equator; most of the radiation flux comes from intermediate latitudes (30 − 40°), while the most variable regions are at higher latitudes. In addition, the SL itself cannot work as a proper waveguide because of the very strong velocity shear that exceeds the Keplerian frequency.
Fig. 16.
Time and longitudeaveraged energy dissipation (upper panel) and radiation flux (lower panel) for the highaccretionrate 3LR simulation run. Mean values and rootmeansquare deviations are shown, respectively, with black solid and red dotted curves. In each panel, the relevant quantity is normalised by the maximal averaged value. 
We can only conclude that the oscillations present in the observational data and in simulations are not the resonance frequencies for sonic waves but rather correspond to a different type of oscillation. The best candidate for these oscillation modes are rmodes (i.e. Rossby waves). As we show in Appendix C, their frequencies at a given colatitude θ form an equidistant spectrum with
where
is the epicyclic frequency (the frequency at which a portion of matter conserving its angular momentum and only affected by gravity and inertial forces would oscillate in a latitudinal direction), Ω = Ω(θ) is the rotation frequency, and m is a whole number. For rigidbody rotation, Ω_{e} ≃ 2Ω cos θ. If the variability is excited in a slowly rotating region outside the SL itself, and the epicyclic frequency changes slowly throughout this region, we see one peak corresponding to the nonrotating m = 0 mode at f ≲ 2f_{spin} and aliases at frequencies differing by Δf ≃ f_{spin}. This is similar to the spectra obtained in our simulations (see Sect. 3.4), and at the same time similar to the pair of QPOs in LMXBs. As seen in Fig. 17, there is a maximum of epicyclic frequency roughly in the interaction region between the SL and the slowly rotating matter. In addition, the epicyclic frequency in this region is very close to the local rotation frequency, meaning that the perturbations are in resonance with rotation. As most variability comes from higher latitudes, we propose that the oscillations are excited at intermediate latitudes (30 − 50° for 3LR), probably by shear instabilities in the interaction region, and then propagate towards the poles.
Fig. 17.
Epicyclic (black solid) and local rotation (red dotted) frequencies for the model 3LR, averaged in time between 0.3 and 0.5 s and over the azimuthal angle. Green dashed horizontal lines correspond to the rotation rate of the NS and its second harmonic. 
Let us assume that the oscillations are always generated at the latitude of the rim of the SL, all the energy is dissipated within the layer, and the local flux is equal to the Eddington flux cg_{eff}/𝜘. Flux scaling with the Eddington limit means that the luminosity should grow approximately linearly with the surface area of the layer, L ≃ L_{Edd} cos θ_{SL}. The epicyclic frequency should therefore scale as
which reproduces the characteristic values of the frequency in our simulations but somewhat overestimates the dependence on flux. As the radiating region does not exactly coincide with the SL, and the width of the SL also depends on the parameters of the inflow (the extreme case being the case of a strongly inclined source; see Eq. (39)), and g_{eff} is generally smaller than gravity, we expect the linear scaling to be a very crude approximation, overpredicting the slope of the actual (seen in simulations) frequency dependence on flux.
A similar type of QPO spectrum consisting of the local epicyclic frequency and its aliases with the rotation frequency was obtained by Erkut et al. (2008) and Belyaev (2017) who considered a BL as a part of the accretion disc. However, in the case considered in these papers, both Ω_{e} and Ω are close to the local Keplerian frequency and are not supposed to be sensitive to the spin of the NS. Belyaev (2017) also studied the conditions for the excitation of the oscillation modes, which are to some degree also applicable to our results, as the existence of a strong velocity shear in combination with differential rotation is a universal feature of any BL model. Shear instabilities may be excited without an initial velocity discontinuity, but the spatial scales of velocity variations should be smaller than the size of the simulation domain (Belyaev & Rafikov 2012), and the velocity profile should have an inflection point (Hertfelder & Kley 2015). However, neither this model nor our simulations are so far capable of explaining clearly why only two peaks are observed in the PDSs of real LMXBs (though see above, Sect. 3.4), and why, in a large number of LMXBs, the distance between the peaks is actually half of the spin frequency. Explaining and predicting the details of QPO features in the PDS requires more profound studies, both analytical and numerical. Both classical and spreading layer approaches have their limitations, as the real motions are likely to be threedimensional (Babkovskaia et al. 2008).
The amplitudes of the oscillations observed in our simulations are rather small, of the order 10^{−4}, but nevertheless the peaks themselves are clearly significant. This is much smaller than the observed ∼10% (Méndez et al. 2001), probably due to several reasons. First, the observational data provide us with spectral variability that does not always follow the variations of the bolometric flux. It appears that the large amplitudes of the kHz QPOs in harder Xrays (E ≳ 5 keV) are related to the temperature variations of the radiating surface, converted in the blackbody approximation to exponentially strong monochromatic flux variations. Still, the temperature variations required to reproduce the observed flux variability is several per cent. Another reason could be related to the algorithm we use to calculate the observables: it does not take into account either relativistic effects or the shape of the photosphere of the SL. In some of our models, the vertical thickness of the layer reaches hundreds of meters. The kilometerscale variations of the shape of the photosphere potentially have important implications for the observed variability of BLs. Last but probably most important, we cannot exclude that the oscillations generated in the SL resonate or are amplified elsewhere, such as for example in an optically thin hot corona or in the accretion disc. This suggestion, though speculative, can help us to understand why only two harmonics are normally seen: structures more extended than BLs are unlikely to have resonance frequencies higher than ∼1.5 kHz.
5. Conclusions
In this paper, we consider a timedependent hydrodynamic SL on the surface of a NS. We used twodimensional spectral modelling to resolve the evolution of the differentially rotating flow. We find that, though challenging due to the supersonic compressible nature of the flow, spectral simulations of a SL on the surface of a NS may be quite productive. We mainly consider the interaction of a new material rotating close to Keplerian velocity with the old, spundown atmosphere of the NS, and this interaction produces a set of hydrodynamical phenomena that have a huge impact on the dynamics of the system. In particular, the velocity shear is susceptible to shear instability modes that provide angular momentum transfer within the layer and excite inertial oscillation modes closer to the poles where they produce variability patterns closely resembling kHz QPOs in real LMXB systems.
We also used a wrapper class spharmt (https://gist.github.com/jswhit/3845307) for shtns quantities and operators written by Jeffrey Whitaker.
Acknowledgments
We would like to thank Alexander Philippov for his help with the spectral filtering, and Roman Rafikov for discussions about different oscillation modes. We also thank Yuri Levin, Andrei Beloborodov, Andrei Gruzinov, Lorenzo Sironi, Joe Patterson, and Anatoly Spitkovsky for discussions on the different aspects of boundary layer physics and the referee for numerous valuable comments. PA acknowledges the support from the Program of development of M. V. Lomonosov Moscow State University (Leading Scientific School “Physics of stars, relativistic objects and galaxies”). JP was supported by the grant 14.W03.31.0021 of the Ministry of Science and Higher Education of the Russian Federation. We acknowledge the use of the Finnish Grid and Cloud Infrastructure and the Finnish IT Center for Science (CSC).
References
 Babkovskaia, N., Brandenburg, A., & Poutanen, J. 2008, MNRAS, 386, 1038 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Belyaev, M. A. 2017, ApJ, 835, 238 [NASA ADS] [CrossRef] [Google Scholar]
 Belyaev, M. A., & Quataert, E. 2018, MNRAS, 479, 1528 [CrossRef] [Google Scholar]
 Belyaev, M. A., & Rafikov, R. R. 2012, ApJ, 752, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Belyaev, M. A., Rafikov, R. R., & Stone, J. M. 2013, ApJ, 770, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Cooley, J. W., & Tukey, J. W. 1965, Math. Comput., 19, 297 [CrossRef] [MathSciNet] [Google Scholar]
 Courant, R., Friedrichs, K., & Lewy, H. 1928, Math. Ann., 100, 32 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Erkut, M. H., Psaltis, D., & Alpar, M. A. 2008, ApJ, 687, 1220 [CrossRef] [Google Scholar]
 Gilfanov, M., Revnivtsev, M., & Molkov, S. 2003, A&A, 410, 217 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gottlieb, D., & Shu, C.W. 1997, SIAM Rev., 39, 644 [CrossRef] [MathSciNet] [Google Scholar]
 Hasinger, G., & van der Klis, M. 1989, A&A, 225, 79 [NASA ADS] [Google Scholar]
 Hertfelder, M., & Kley, W. 2015, A&A, 579, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Inogamov, N. A., & Sunyaev, R. A. 1999, Astron. Lett., 25, 269 [NASA ADS] [Google Scholar]
 Inogamov, N. A., & Sunyaev, R. A. 2010, Astron. Lett., 36, 848 [NASA ADS] [CrossRef] [Google Scholar]
 JakobChien, R., Hack, J., & Williamson, D. 1995, J. Comput. Phys., 119, 164 [CrossRef] [Google Scholar]
 Kluzniak, W., & Abramowicz, M. A. 2001, Acta Phys. Polonica B, 32, 3605 [Google Scholar]
 Kluźniak, W., Abramowicz, M. A., Kato, S., Lee, W. H., & Stergioulas, N. 2004, ApJ, 603, L89 [Google Scholar]
 Lamb, F. K., Pethick, C. J., & Pines, D. 1973, ApJ, 184, 271 [NASA ADS] [CrossRef] [Google Scholar]
 Lamb, F. K., Shibazaki, N., Alpar, M. A., & Shaham, J. 1985, Nature, 317, 681 [NASA ADS] [CrossRef] [Google Scholar]
 Méndez, M., & Belloni, T. 2007, MNRAS, 381, 790 [NASA ADS] [CrossRef] [Google Scholar]
 Méndez, M., van der Klis, M., & van Paradijs, J. 1998, ApJ, 506, L117 [NASA ADS] [CrossRef] [Google Scholar]
 Méndez, M., van der Klis, M., Ford, E. C., Wijnands, R., & van Paradijs, J. 1999, ApJ, 511, L49 [NASA ADS] [CrossRef] [Google Scholar]
 Méndez, M., van der Klis, M., & Ford, E. C. 2001, ApJ, 561, 1016 [NASA ADS] [CrossRef] [Google Scholar]
 Miller, M. C., & Lamb, F. K. 2016, Eur. Phys. J. A, 52, 63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Miller, M. C., Lamb, F. K., & Psaltis, D. 1998, ApJ, 508, 791 [NASA ADS] [CrossRef] [Google Scholar]
 Miyamoto, S., Kimura, K., Kitamoto, S., Dotani, T., & Ebisawa, K. 1991, ApJ, 383, 784 [NASA ADS] [CrossRef] [Google Scholar]
 Monin, A., Yaglom, A., & Lumley, J. 2007, Statistical Fluid Mechanics: Mechanics of Turbulence, Dover books on physics No. v. 1 (Dover Publications) [Google Scholar]
 Motta, S. E., Rouco Escorial, A., Kuulkers, E., MuñozDarias, T., & Sanna, A. 2017, MNRAS, 468, 2311 [CrossRef] [Google Scholar]
 Nättilä, J., Miller, M. C., Steiner, A. W., et al. 2017, A&A, 608, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nowak, M. A., Vaughan, B. A., Wilms, J., Dove, J. B., & Begelman, M. C. 1999, ApJ, 510, 874 [NASA ADS] [CrossRef] [Google Scholar]
 Papaloizou, J. C. B., & Stanley, G. Q. G. 1986, MNRAS, 220, 593 [NASA ADS] [CrossRef] [Google Scholar]
 Parfrey, K., Beloborodov, A. M., & Hui, L. 2012, MNRAS, 423, 1416 [NASA ADS] [CrossRef] [Google Scholar]
 Philippov, A. A., Rafikov, R. R., & Stone, J. M. 2016, ApJ, 817, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Piro, A. L., & Bildsten, L. 2004a, ApJ, 610, 977 [NASA ADS] [CrossRef] [Google Scholar]
 Piro, A. L., & Bildsten, L. 2004b, ApJ, 616, L155 [NASA ADS] [CrossRef] [Google Scholar]
 Popham, R., & Narayan, R. 1995, ApJ, 442, 337 [NASA ADS] [CrossRef] [Google Scholar]
 Pringle, J. E. 1977, MNRAS, 178, 195 [NASA ADS] [CrossRef] [Google Scholar]
 Psaltis, D., Méndez, M., Wijnands, R., et al. 1998, ApJ, 501, L95 [NASA ADS] [CrossRef] [Google Scholar]
 Ray, T. P. 1982, MNRAS, 198, 617 [NASA ADS] [CrossRef] [Google Scholar]
 Rebusco, P. 2008, New A Rev., 51, 855 [NASA ADS] [CrossRef] [Google Scholar]
 Revnivtsev, M. G., & Gilfanov, M. R. 2006, A&A, 453, 253 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rincon, F., Ogilvie, G. I., & Cossu, C. 2007, A&A, 463, 817 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schaeffer, N. 2013, Geochem. Geophys. Geosyst., 14, 751 [NASA ADS] [CrossRef] [Google Scholar]
 Sibgatullin, N. R., & Sunyaev, R. A. 2000, Astron. Lett., 26, 699 [NASA ADS] [CrossRef] [Google Scholar]
 Suleimanov, V., & Poutanen, J. 2006, MNRAS, 369, 2036 [NASA ADS] [CrossRef] [Google Scholar]
 Tikhonov, A., & Samarskii, A. 2013, Equations of Mathematical Physics, Dover Books on Physics (Dover Publications) [Google Scholar]
 van der Klis, M. 2000, ARA&A, 38, 717 [NASA ADS] [CrossRef] [Google Scholar]
 van der Klis, M. 2001, ApJ, 561, 943 [NASA ADS] [CrossRef] [Google Scholar]
 Vreugdenhil, C. B. 1990, VKI, Computational Fluid Dynamics, 1, 48 (SEE N90–27989 22–34), 1 [Google Scholar]
 Wijnands, R., van der Klis, M., Homan, J., et al. 2003, Nature, 424, 44 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Zhuravlev, V. V., & Razdoburdin, D. N. 2018, A&A, 619, A44 [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Derivation of the equations for divergence and vorticity
To get the equations for δ = ∇ · v and ω = (∇ × v)_{r}, we need to take divergence and curl of the system of dynamical equations containing sources and sinks
where g = −∇Φ is gravity without centrifugal terms, and then integrate them in radial direction. Gravitational potential has the form , where ΔΦ accounts for a nonspherical shape of the star and depends on the latitude and longitude. For the friction force f_{fric}, we use the form . Here, v_{NS} = Ω_{∗}R_{∗} sin_{θ} is the rotation velocity field of the NS itself and is directed azimuthally.
Before deriving the actual equations for vorticity and divergence, we note that, for any velocity field, its vector product with its curl ω = ∇ × v is
Hence,
and the curl of the nonlinear term (v∇)v in Eq. (A.1), after application of the triple vector product rule, takes the form
As ω has zero divergence, the curl of Eq. (A.1) becomes
where ω_{d} = 2c_{K}Ω_{K} cos α is the vorticity in the source of matter (see Sect. 2.3). Taking radial integral of the radial component of Eq. (A.5) is straightforward, as the equation does not contain radial derivatives. Finally, we get Eq. (26):
The righthand side of this equation contains a baroclinic term capable of creating vorticity out of density and pressure variations, and three terms related to the vorticity of the accreted matter ω_{source} and friction with the surface.
Another equation describing the time evolution of δ comes from taking divergence of dynamical equation. It is rather nontrivial to expand the advection term, ∇·((v · ∇)v). Therefore, let us first note that, according to Eq. (A.2),
Here, the last term on the righthand side is identical to the advective lefthandside term in the derivative of Eq. (A.1), hence
The potential term ΔΦ appears because of rotational deformation of the star. As the surface of the star should be an equipotential surface, total gravitational and centrifugal potential
implying .
Appendix B: Numerical implementation and tests
B.1. Code description
The problem we consider is essentially a more physically elaborate version of shallowwater hydrodynamics, supplemented with energy transfer and source and sink terms. As the problem is formulated for the surface of a sphere, it is natural to use a spectral code working with spherical harmonics. We used the shtns library (Schaeffer 2013) designed for hydrodynamical and geophysical applications^{2}. A major challenge of our problem is that, unlike the classical shallowwater physics, it is far from the subsonic Rossbyapproximation motions, and thus the time step is limited by several processes. One of the requirements for the time step is the CourantFriedrichsLewy condition (Courant et al. 1928) of the form
where Δx_{min} is the minimal physical size of a cell in the simulation domain, u is the fastest relevant signal propagation velocity, and C ≲ 1 is a constant related to the particular solver used in simulations. The existence of sources and sinks for several physical quantities sets additional upper limits for the time step and requires us to adjust the time step with the physical conditions. We compute the actual time step as a harmonic sum of several time steps
where C_{sound, adv, thermal, accr} ≳ 1 are dimensionless adjustable parameters regulating the role of each time step (note that a larger coefficient here means a smaller time step). Thermal, Δt_{thermal} = minE/(Q^{+} + Q^{−}), and accretion, Δt_{accr} = minΣ/(S^{+} + S^{−}), time steps are estimated as the time steps required to resolve temporally thermal and accretion processes, respectively. Such a choice for a variable time step ensures that all the relevant physical processes can be resolved: sonic wave propagation, dissipation, radiation losses, and accretion.
As spectral methods tend to produce highfrequency noise, diffusionlike dissipation terms were introduced for all the five principal quantities ω, δ, Σ, E, and a. This is done for the numerical implementation of Eqs. (8), (26), (27), and (29). In spectral space, an additional dissipation term may be viewed as a multiplier cutting off high frequencies (lowpass filter). A usual form for this lowpass filter (see e.g. Parfrey et al. 2012) is hyperdiffusion
where is the Laplacian acting on the spherical harmonic of degree l (see Tikhonov & Samarskii 2013), L_{max} is its maximal value (corresponding to the grid resolution), t_{D} is the characteristic dissipation timescale on the size of the grid, and N_{diss} ≥ 2 is a free parameter describing the shape of the lowpass filter. If N_{diss} = 2, the filtering procedure is equivalent to a regular diffusion term added to the righthand sides of all the main differential equations. The flow as a whole is negligibly affected if the dissipation factor for the lowestorder harmonic is indistinguishable from zero t_{D} = ≳lne_{M}Δt (Parfrey et al. 2012), where e_{M} is machine precision (e_{M} ≃ 2 × 10^{−16} in all the simulations we run).
Effectively, filtering with N_{diss} > 2, while more efficient than normal diffusion, cuts all the high frequencies very sharply. Such spectral truncation leads to its own noise, especially close to discontinuities. This is known as Gibbs phenomenon (Gottlieb & Shu 1997) and leads to oscillations that may be amplified by the nonlinear nature of our system of equations. For Σ and E that can vary sharply and span several orders of magnitude, but become unphysical if negative, the Gibbs phenomenon may become disastrous. A reasonable solution introducing the little Gibbs effect and providing an extremely accurate treatment to the largescale flow is
which we use for all the simulations in this paper.
The energy lost by the flow is added as a source of internal energy, as described in Sect. 2.5. By adding dissipation as a source of energy, we are likely to introduce highfrequency noise to the energy field, and therefore the dissipation field was smoothed in the same way as the basic quantities (see Eq. (B.3)), but with a shorter diffusion timescale. The exact value of the dissipation smoothing parameter affects the thermal stability of the simulation but does not change the overall dynamics. From the physical point of view, it only ensures that the dissipation does not significantly vary within a single resolution element.
The code itself is written as a hybrid PYTHON3/C++ program. All the numerically heavy (spectral) calculations are solved with the C++ SHTNS library, whereas the main loop and the related highlevel functionalities are operated from the more userfriendly PYTHON3 driver. In our experience, this provides a good balance between numerical efficiency and ease of use. The spherical harmonic calculations are parallelised using the shared memory paradigm with OPENMP pragmas. This enables us to take advantage of multicore platforms ranging from powerful desktop computers to occupying one complete node in computing clusters. The Hierarchical Data Format (HDF5) is used to save and store the simulation results.
B.2. Tests
In Table B.1, we list the test models we calculated with their basic parameters.
Test simulations.
B.2.1. Zeroaccretionrate, rigidbody rotation case
As one of the tests, we try evolution of a layer with initial surface density Σ_{0} = 10^{8} g cm^{−2} and sound velocity c_{s} ≃ 1.7 × 10^{−3}c without accretion or depletion (test models NDLR/NDHR). As in all the other models, an initial perturbance of 5% was introduced. The NS spin period was set to 3 ms. The Mach number of this flow is about 50. For this test, we also turn off dissipation heating and radiation losses. Without thermal effects, rotation profile in such a model should not change with time, and the perturbation proceeds rotating with the surface of the star. To check the accuracy of this solution, it is sufficient to correct for the rotation angle, interpolate from one grid to another, and estimate the standard deviation or the maximal deviation between the map calculated by the code and the interpolated initial map.
Figure B.1 shows how the maximal relative difference in surface density evolves with time. The errors around 10^{−3} are interpolation errors. As we can see, supersonic rotation is reasonably well tracked for multiple rotation periods, and the accuracy is better for a finer grid.
Fig. B.1.
Maximal relative error in surface density for the test models NDLR (black) and NDHR (red). The blue horizontal segment in the lower panel has the length of one spin period. The dotted green horizontal line corresponds to the amplitude of the initial perturbation, 0.05. 
B.2.2. Splitsphere tests
The purpose of this test set (twistLR, twistHR, stwistLR, and stwistHR) was to trace the development of sub and supersonic shear instabilities on a sphere. Rigidbody rotation (P_{spin} = 10 and 30 ms) was modified by a factor rapidly changing from −1 to 1 near the equator
where Δθ was set to 0.1 for all the models. The choice of the effective temperature (set by Q_{NS}, see Sect. 2.5) makes some simulations subsonic and others supersonic. For Δθ ≪ π, subsonic configuration is unstable to KelvinHelmholtz instability. Instability at high wavenumbers is suppressed by the finite shear value (Ray 1982), hence the fastestgrowing unstable modes are two and threearmed, with the increment of about Ω_{*} (see Fig. B.2). The sharper the velocity gradient, the higher the fastestgrowing mode. Conservation of angular momentum prevents the formation of a single vortex. The primary instability mode changes the overall velocity field into a set of vortices centered in the initial equatorial region. Vorticity evolution during the instability development phase is shown in Fig. B.3.
Fig. B.2.
Energy evolution and relaxation for the splitsphere test, sub (left panel) and supersonic (right panel) cases. Black lines correspond to the part of kinetic energy related to v_{θ}, red to v_{φ}. Dotted lines are used for lowerresolution models (s)twistLR, solid lines for highresolution (s)twistHR. Blue dashed lines show an exponential law ∝e^{Ω*t}. 
Fig. B.3.
Four vorticity snapshots (t = 0.04, 0.05, 0.06, and 0.08 s) of the Kelvin–Helmholtz instability development in the splitsphere simulation, model twistHR. 
For the mildly supersonic splitsphere test (P_{spin} = 10 ms), a supersonic shear instability develops on similar timescales close to the rotation period is used as the basis for the splitsphere rotation. However, instead of vortices, a system of standing shock waves is formed (see Fig. B.4).
Fig. B.4.
Snapshot of the stwistLR simulation at t = 0.04 s, when the instability is fully developed and evolves into a system of shock waves. Vorticity is shown in the left and the surface density in the right panel. 
As we can see, development of shear instabilities conforms well with the expectations based on analytical and numerical studies of the subject. First, there is dynamicaltimescale exponential growth of the instability, subsequently evolving into an equipartition turbulent stage. For the subsonic case, numerical resolution does not affect the results considerably during the time span of the calculations.
Appendix C: Frequencies of inertial modes on a differentially rotating sphere
Assuming an axisymmetric, differentially rotating velocity background, we linearise the set of dynamic equations and derive a dispersion relation for smallamplitude shallowwater waves on a unit sphere. Perturbed quantities to be considered are density ρ = ρ_{0} + δρ(θ, φ, t), longitudinal velocity v_{φ} = Ω(θ) sin θ + δv_{φ}(θ, φ, t), and latitudinal velocity v_{θ} = δv_{θ}(θ, φ, t), where the terms with δ are small perturbations. The background flow is assumed to be pure differential rotation parametrized by angular velocity distribution Ω(θ). All the perturbations are expressed in exponential form ∝exp(i(ωt − k_{θ}θ − k_{φ}φ)).
Firstorder perturbation of the continuity equation in such assumptions is
The two tangential Euler equations may, in general form, ignoring the terms containing radial velocities, be written as
and
For a supersonic flow, contributions of the pressure variations on the righthand side of the equations are of secondary importance, though in reality they are responsible for pressure and gravity oscillation modes. If we ignore the pressure variations, the two Euler equations become, respectively,
and
where . Excluding the velocity components v_{θ} and δv_{φ} from Eqs. (C.5) and (C.4) yields a dispersion equation
where
is the square of the local epicyclic frequency in the sense that a particle with a conserved angular momentum, confined to the surface of the star and being a subject of gravity and centrifugal force, will oscillate in latitudinal direction at this frequency. The possible values of k_{φ} are restricted by the longitudinal periodic boundary conditions to be k_{φ} = m, where m is a whole number. Thus, the spectrum of possible inertial oscillation frequencies takes the form
Without any loss of generality, we choose the sign in Eq. (C.8) to be plus. In the case of m = 0 and Ω ≃ Ω_{*}, the only axisymmetric inertial mode has ω_{inertial, 0} = Ω_{e} ≃ 2Ω cos θ reproducing the Coriolis oscillation regime. Variability occurring in the regions corotating with the NS would produce an equidistant spectrum of eigenmodes
with the frequencies differing by the rotation frequency of the star.
All Tables
All Figures
Fig. 1.
Illustration of the model geometry. The tilted blue arc near the equator shows the source of mass and momentum. The spin axis of the star, marked with Ω_{*}, is inclined with respect to the disc axis (Ω_{d}), by an angle i. 

In the text 
Fig. 2.
Timelatitude diagrams for longitudinally averaged angular frequency normalised to the Keplerian rotation frequency at the radius of the star R_{*}. The upper and lower panels correspond to the rapid accretion model 3LR and 8LR, respectively. 

In the text 
Fig. 3.
Timelatitude diagram for longitudinally averaged surface density (normalised by the surface density averaged over the sphere) in the realisticaccretionrate 8LR simulation. Upper panels: three snapshots of surface density (10^{4} g cm^{−2} units) during the later stages of the heating instability development. 

In the text 
Fig. 4.
Radiation flux (in Eddington units, ) snapshots for the model 3HR. White streamlines show the velocity field and the black contour corresponds to the accretion tracer a = 0.5. For all but the first snapshot, newly accreted matter dominates between the black contours. 

In the text 
Fig. 5.
Resolution effects after t = 0.03 s (ten spin periods) of evolution for the highaccretionrate models 3LR (low numerical resolution; left panels) and 3HR (high numerical resolution; right panels). Top: vorticity maps. Bottom: emitted radiation flux Q^{−} (normalised by the local Eddington value) map. 

In the text 
Fig. 6.
Azimuthally and temporally averaged dynamical quantities for the highaccretionrate model 3LR. Visualised quantities are Reynolds stress (black), mean velocity product v_{θ}v_{φ} (blue), and sound velocity (red dotted) as functions of latitude. All the quantities were averaged over the period of time between 0.1 and 0.3 s, and over the azimuthal angle. Solid and dashed lines correspond to positive (southward motion or momentum transfer) and negative quantities (northward). Dotted black curves show the 1σ standard deviation interval for the Reynolds stress. 

In the text 
Fig. 7.
Azimuthally and temporally averaged dynamical quantities for the realisticaccretionrate model 8LR. Symbols and quantities are the same as in Fig. 6. Averaging was performed from t = 0.2 to 0.28 s, when the initial hemisphere asymmetry is developed due to heating instability. 

In the text 
Fig. 8.
Dynamical PDSs (log_{10}(f^{2}PDS(f)) normalised to total power within a single time bin) for the “isotropic luminosity” calculated as described in Sect. 3.4 for the highaccretionrate model 3LR. The three upper panels show dynamical spectra for the observer’s inclinations of π/2, π/4, and 0, respectively. White horizontal lines show the spin (lower) and Keplerian frequencies. Lowermost panel: corresponding light curves: i_{obs} = π/2, π/4, and 0 cases are shown with blue dotted, green dashed, and black solid lines. 

In the text 
Fig. 9.
Same as Fig. 8, but for the model with a turnedoff mass source, 3LRoff. 

In the text 
Fig. 10.
Same as Fig. 8, but for the model with an inclined source, 3LRinc. 

In the text 
Fig. 11.
Peak frequency (calculated as the position of the maximum of f × PDS) as a function of observed luminosity (calculated using expression (36)) for i_{obs} = 0 (black circles), π/4 (red diamonds), and π/2 (blue triangles). Simulation runs 3LR (left) and 3LRoff (right). Error bars show flux dispersion within the time bin and the size of the frequency bin where the maximum was detected. Horizontal dotted green lines show spin frequency. 

In the text 
Fig. 12.
Integral PDS (time span 0.2 to 0.5 s) of the simulation 3LR for different inclinations. Black dots correspond to a poleon observer i_{obs} = 0, red diamonds to the intermediate inclination of i_{obs} = π/4, and green triangles to i_{obs} = π/2. 

In the text 
Fig. 13.
Same as Fig. 12, but for the simulation 3LRoff and between t = 0.58 and 0.67 s. 

In the text 
Fig. 14.
Same as Fig. 12, but for the simulation 3LRinc and for the time span t = 0.35 − 0.5 s. 

In the text 
Fig. 15.
Integral PDS (time span from 0.27 to 0.32 s) for the simulation 8LR for different inclinations (black circles i_{obs} = 0, red diamonds i_{obs} = 45°, green upward triangles i_{obs} = 90°, and blue downward triangles i_{obs} = 180°). 

In the text 
Fig. 16.
Time and longitudeaveraged energy dissipation (upper panel) and radiation flux (lower panel) for the highaccretionrate 3LR simulation run. Mean values and rootmeansquare deviations are shown, respectively, with black solid and red dotted curves. In each panel, the relevant quantity is normalised by the maximal averaged value. 

In the text 
Fig. 17.
Epicyclic (black solid) and local rotation (red dotted) frequencies for the model 3LR, averaged in time between 0.3 and 0.5 s and over the azimuthal angle. Green dashed horizontal lines correspond to the rotation rate of the NS and its second harmonic. 

In the text 
Fig. B.1.
Maximal relative error in surface density for the test models NDLR (black) and NDHR (red). The blue horizontal segment in the lower panel has the length of one spin period. The dotted green horizontal line corresponds to the amplitude of the initial perturbation, 0.05. 

In the text 
Fig. B.2.
Energy evolution and relaxation for the splitsphere test, sub (left panel) and supersonic (right panel) cases. Black lines correspond to the part of kinetic energy related to v_{θ}, red to v_{φ}. Dotted lines are used for lowerresolution models (s)twistLR, solid lines for highresolution (s)twistHR. Blue dashed lines show an exponential law ∝e^{Ω*t}. 

In the text 
Fig. B.3.
Four vorticity snapshots (t = 0.04, 0.05, 0.06, and 0.08 s) of the Kelvin–Helmholtz instability development in the splitsphere simulation, model twistHR. 

In the text 
Fig. B.4.
Snapshot of the stwistLR simulation at t = 0.04 s, when the instability is fully developed and evolves into a system of shock waves. Vorticity is shown in the left and the surface density in the right panel. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.