© ESO 2020

1. Introduction

For a non-magnetised accretor with a fluid or solid surface, non-relativistic 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⁸ 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 low-mass X-ray binaries (LMXBs). In particular, the BL apparently plays an important role in the so-called 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 highest-frequency, kilohertz quasi-periodic 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 low-frequency QPO types are quite similar in NS and black-hole 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 non-axisymmetric structures at the surface of the NS. This interpretation is known as a beat-frequency 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 beat-frequency models that do not involve a BL were proposed to explain the properties of kHz QPOs. The best known are the magnetospheric beat-frequency model (Psaltis et al. 1998), considering LMXBs as magnetised accretors with very compact magnetospheres, and the sonic-point 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 narrow-band 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).

Quasi-periodic 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 QPO-frequency 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 non-Keplerian 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 magneto-rotational 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 Kelvin-Helmholtz instability, oblique waves rather than vortices are generated. Moving in a shear velocity flow, the waves create Reynolds stress and thus provide effective non-local viscosity not only in the BL but also in the accretion disc. Most of the numerical studies mentioned above considered a two-dimensional 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 three-dimensional 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 wave-mediated 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 two-dimensional (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 so-called dwarf-nova 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.

Two-dimensional 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 (Jakob-Chien et al. 1995). Spectral methods provide much higher accuracy than finite-difference 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 open-source software ¹.

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

(1)

which is the approximate light-crossing or dynamical timescale at the event horizon. The corresponding radius is

(2)

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⁻²) for Σ and internally normalise the surface density by an arbitrary scale set either to 10⁴ or to 10⁸ g cm⁻² for the simulations presented in this paper. Evidence for a very long correlation timescale t_corr suggests a characteristic value of

(3)

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² g⁻¹.

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:

(4)

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¹⁵ − 10¹⁸ g s⁻¹ (10⁻¹¹–10⁻⁸ M_⊙ yr⁻¹) which implies effective temperatures of the order

(5)

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⁸ g cm⁻² 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,

(6)

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 radiation-pressure 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 ∝ β⁻¹). 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⁻⁴, 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³ g cm⁻².

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 thermal-timescale phenomena we intend to consider appear in fact in a quasi-stationary 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 pre-existing 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

(7)

Integration over r yields

(8)

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),

(9)

where α is the angular distance from the direction of the adopted disc rotation axis,

(10)

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, co-latitude θ, and longitude φ. Instead of θ, latitude λ = π/2 − θ is used below for visualisation. The system is aligned with the rotation of the star, but is itself non-rotating (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 quasi-stationary 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

(11)

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 steady-state 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 pre-existing 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,

(12)

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

(13)

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

(14)

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 gas-to-total pressure ratio equals

(15)

Proportionality of pressures also implies p ∝ ρT ∝ T⁴, which leads to p ∝ ρ^4/3, an effectively polytropic law. This implies a relation between the pressure p₀ and density ρ₀ at the bottom of the SL and the corresponding vertically integrated quantities, pressure Π = ∫pdr and surface density Σ = ∫ρdr,

(16)

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

(17)

which implies an identical relation for the vertically integrated quantities

(18)

The pressure ratio itself may be found as a function of E (or Π) and Σ. At the bottom of the layer,

(19)

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₀ = g_effΣ arising as a solution to Eq. (13) and Eq. (18) and substituting them into Eq. (19):

(20)

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

(21)

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

(22)

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 non-equilibrium stresses in the crust), which makes the polar-angle component , where Φ is gravitational potential, non-zero 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 re-write in terms of the two scalar quantities normally used in shallow-water approximation: vorticity,

(23)

and divergence,

(24)

Multiplying Eq. (22) by ρ and integrating over the total vertical extent of the SL yields

(25)

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

(26)

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 right-hand 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 δ

(27)

The term originates from the rotational deformation of the NS.

2.5. Energy conservation

In general form, energy conservation implies (Suleimanov & Poutanen 2006):

(28)

where the right-hand 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

(29)

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 pre-existing material. Vertically integrated pressure and energy are related by (18). Dissipation is calculated as kinetic energy lost by the flow:

(30)

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):

(31)

where 𝜘 is the Rosseland average opacity which we assume to be equal to the Thomson scattering opacity, 𝜘 = 𝜘_T ≃ 0.34 cm² g⁻¹. 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 non-zero, 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 pre-existing 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, hence

(32)

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 quasi-stationary, one-dimensional 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 co-rotating 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 rigid-body 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 rigid-body rotation is

(33)

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 under-density. 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 high-resolution 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.

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 steady-state accretion from a thin disc. To avoid strong velocity gradients causing high-frequency noise, the mass accretion rate approaches the steady-state value smoothly, following the exponential law

(34)

where the turn-on 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.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 dynamic-timescale 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⁻³ M_⊙ yr⁻¹ ≃ 6 × 10²⁰ g s⁻¹ (models 3LR and 3HR); secondly, we consider accretion on top of a thin layer with Σ₀ = 10⁴ g cm⁻². 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. Time-latitude 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 north-south symmetry breaking between t≃200 and 250 ms, and later dynamical-timescale evolution. In these simulation runs, β ≃ 0.99–1, close to 1 with an accuracy of about 10⁻⁵ outside the regions of intense mixing and heating instability.

Fig. 3. Time-latitude diagram for longitudinally averaged surface density (normalised by the surface density averaged over the sphere) in the realistic-accretion-rate 8LR simulation. Upper panels: three snapshots of surface density (104 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 gas-to-total 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 quasi-stationary 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⁻⁸ M_⊙ yr⁻¹) and low initial surface density (Σ₀ = 10⁴ g cm⁻²), 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 later-time evolution (starting at approximately t ∼ 0.2 s) is marked by a gradual and then dynamical-timescale 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 low-resolution 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 high-resolution runs, more fine-scale substructure starts to appear.

Fig. 5. Resolution effects after t = 0.03 s (ten spin periods) of evolution for the high-accretion-rate 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 high-accretion-rate 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 co-rotating with the star is a consequence of angular momentum exchange with the star. Whether or not the velocity drop exists in a quasi-stationary 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 non-linear regime, and forms a non-axisymmetric 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

(35)

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 quasi-stationary axisymmetric state on a timescale of about 1 s, higher latitudes still show variability and non-axisymmetry 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 high-accretion-rate 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 oblique-wave-mediated angular momentum transfer. Near the poles, Reynolds stress partially compensates the advective angular momentum flow caused by compression of the pre-existing polar cap material. Near the equator, at the same time, both advection and viscous transport of angular momentum are directed polewards.

For the realistic-accretion-rate 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 realistic-accretion-rate 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:

(36)

where

(37)

and α_obs is the angle at which the surface element is seen by an observer located at the co-latitude 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

(38)

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) power-density spectra for different i_obs. We plot the relative PDS multiplied by f² to decrease the contamination from the low-frequency noise component related to the overall shape of the simulated light curve. Several oscillation modes are visible, one of them for a pole-on observer. Their frequencies evidently correlate with the flux. Most of the non-axisymmetric 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 non-axisymmetric modes absent in the pole-on 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 (log10(f2PDS(f)) normalised to total power within a single time bin) for the “isotropic luminosity” calculated as described in Sect. 3.4 for the high-accretion-rate 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: iobs = π/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 quasi-periodic features around two-three 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.

Fig. 9. Same as Fig. 8, but for the model with a turned-off mass source, 3LRoff.

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 quasi-axisymmetric 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 pole-on observer is

(39)

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

Starting from t ≃ 0.15 s, the pole-on 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, non-axisymmetric 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 pre-existing 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 pole-on simulated light curve both in the original model with enhanced accretion, and in the switch-off simulation. The axisymmetric mode discussed above dominates for the pole-on case, for which a clear correlation between flux and frequency is observed. For an inclined observer, non-axisymmetric 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 iobs = 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. Pole-on 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 gas-pressure-dominated 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 pole-on observer iobs = 0, red diamonds to the intermediate inclination of iobs = π/4, and green triangles to iobs = π/2.

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

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

For the realistic mass-accretion-rate 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 iobs = 0, red diamonds iobs = 45°, green upward triangles iobs = 90°, and blue downward triangles iobs = 180°).

4. Discussion

In the dynamical spectra presented in Sect. 3.4, there are clearly at least two types of quasi-periodic variability signals: one disappears for a pole-on observer and is thus related to non-axisymmetric 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 dwarf-nova oscillations (DNOs). However, the width of the SL, independently from the assumptions, should increase with mass accretion rate. This is especially true for the radiation-pressure-supported 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 DNO-like sonic mode for a thin radiation-pressure-supported SL is

(40)

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 longitude-averaged energy dissipation (upper panel) and radiation flux (lower panel) for the high-accretion-rate 3LR simulation run. Mean values and root-mean-square 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 r-modes (i.e. Rossby waves). As we show in Appendix C, their frequencies at a given co-latitude θ form an equidistant spectrum with

(41)

where

(42)

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 rigid-body 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 non-rotating 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

(43)

which reproduces the characteristic values of the frequency in our simulations but somewhat over-estimates 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, over-predicting 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 three-dimensional (Babkovskaia et al. 2008).

The amplitudes of the oscillations observed in our simulations are rather small, of the order 10⁻⁴, 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 X-rays (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 kilometer-scale 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 time-dependent hydrodynamic SL on the surface of a NS. We used two-dimensional spectral modelling to resolve the evolution of the differentially rotating flow. We find that, though challenging due to the super-sonic 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, spun-down 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.

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

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

(A.1)

where g = −∇Φ is gravity without centrifugal terms, and then integrate them in radial direction. Gravitational potential has the form , where ΔΦ accounts for a non-spherical 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

(A.2)

Hence,

(A.3)

and the curl of the non-linear term (v∇)v in Eq. (A.1), after application of the triple vector product rule, takes the form

(A.4)

As ω has zero divergence, the curl of Eq. (A.1) becomes

(A.5)

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):

(A.6)

The right-hand 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 non-trivial to expand the advection term, ∇·((v · ∇)v). Therefore, let us first note that, according to Eq. (A.2),

(A.7)

Here, the last term on the right-hand side is identical to the advective left-hand-side term in the derivative of Eq. (A.1), hence

(A.8)

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

(A.9)

implying .

Appendix B: Numerical implementation and tests

B.1. Code description

The problem we consider is essentially a more physically elaborate version of shallow-water 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². A major challenge of our problem is that, unlike the classical shallow-water physics, it is far from the subsonic Rossby-approximation motions, and thus the time step is limited by several processes. One of the requirements for the time step is the Courant-Friedrichs-Lewy condition (Courant et al. 1928) of the form

(B.1)

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

(B.2)

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 high-frequency noise, diffusion-like 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 (low-pass filter). A usual form for this low-pass filter (see e.g. Parfrey et al. 2012) is hyperdiffusion

(B.3)

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 low-pass filter. If N_diss = 2, the filtering procedure is equivalent to a regular diffusion term added to the right-hand sides of all the main differential equations. The flow as a whole is negligibly affected if the dissipation factor for the lowest-order harmonic is indistinguishable from zero t_D = ≳|lne_M|Δt (Parfrey et al. 2012), where e_M is machine precision (e_M ≃ 2 × 10⁻¹⁶ 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 non-linear 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 large-scale flow is

(B.4)

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 high-frequency 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 high-level functionalities are operated from the more user-friendly 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 multi-core 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.

B.2.1. Zero-accretion-rate, rigid-body rotation case

As one of the tests, we try evolution of a layer with initial surface density Σ₀ = 10⁸ g cm⁻² and sound velocity c_s ≃ 1.7 × 10⁻³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⁻³ 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. Split-sphere tests

The purpose of this test set (twistLR, twistHR, stwistLR, and stwistHR) was to trace the development of sub- and super-sonic shear instabilities on a sphere. Rigid-body rotation (P_spin = 10 and 30 ms) was modified by a factor rapidly changing from −1 to 1 near the equator

(B.5)

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 Kelvin-Helmholtz instability. Instability at high wavenumbers is suppressed by the finite shear value (Ray 1982), hence the fastest-growing unstable modes are two- and three-armed, with the increment of about Ω_* (see Fig. B.2). The sharper the velocity gradient, the higher the fastest-growing 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 split-sphere test, sub- (left panel) and super-sonic (right panel) cases. Black lines correspond to the part of kinetic energy related to vθ, red to vφ. Dotted lines are used for lower-resolution models (s)twistLR, solid lines for high-resolution (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 split-sphere simulation, model twistHR.

For the mildly supersonic split-sphere 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 split-sphere 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 dynamical-timescale 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 small-amplitude shallow-water waves on a unit sphere. Perturbed quantities to be considered are density ρ = ρ₀ + δρ(θ, φ, 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_φφ)).

First-order perturbation of the continuity equation in such assumptions is

(C.1)

The two tangential Euler equations may, in general form, ignoring the terms containing radial velocities, be written as

(C.2)

and

(C.3)

For a super-sonic flow, contributions of the pressure variations on the right-hand 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,

(C.4)

and

(C.5)

where . Excluding the velocity components v_θ and δv_φ from Eqs. (C.5) and (C.4) yields a dispersion equation

(C.6)

where

(C.7)

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

(C.8)

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 co-rotating with the NS would produce an equidistant spectrum of eigenmodes

(C.9)

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. Time-latitude 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. Time-latitude diagram for longitudinally averaged surface density (normalised by the surface density averaged over the sphere) in the realistic-accretion-rate 8LR simulation. Upper panels: three snapshots of surface density (104 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 high-accretion-rate 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 high-accretion-rate 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 realistic-accretion-rate 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 (log10(f2PDS(f)) normalised to total power within a single time bin) for the “isotropic luminosity” calculated as described in Sect. 3.4 for the high-accretion-rate 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: iobs = π/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 turned-off 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 iobs = 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 pole-on observer iobs = 0, red diamonds to the intermediate inclination of iobs = π/4, and green triangles to iobs = π/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 iobs = 0, red diamonds iobs = 45°, green upward triangles iobs = 90°, and blue downward triangles iobs = 180°).
In the text

	Fig. 16. Time- and longitude-averaged energy dissipation (upper panel) and radiation flux (lower panel) for the high-accretion-rate 3LR simulation run. Mean values and root-mean-square 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 split-sphere test, sub- (left panel) and super-sonic (right panel) cases. Black lines correspond to the part of kinetic energy related to vθ, red to vφ. Dotted lines are used for lower-resolution models (s)twistLR, solid lines for high-resolution (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 split-sphere 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