MHD simulations of oscillating cusp-filling tori around neutron stars -- missing upper kHz QPO

TL;DR

This study uses axisymmetric MHD simulations to analyze oscillating cusp-filling tori around neutron stars, revealing dominant radial epicyclic modes and their impact on accretion and luminosity, with implications for understanding kHz QPOs.

Contribution

It presents the first detailed MHD simulation of oscillating cusp-filling tori around neutron stars, highlighting the suppression of vertical oscillations and their effect on accretion modulation.

Findings

Radial epicyclic mode dominates oscillations.

Vertical oscillations are suppressed by the cusp.

Mass accretion rate modulates boundary layer luminosity.

Abstract

We performed axisymmetric, grid-based, ideal magnetohydrodynamic (MHD) simulations of oscillating cusp-filling tori orbiting a non-rotating neutron star. A pseudo-Newtonian potential was used to construct the constant angular momentum tori in equilibrium. The inner edge of the torus is terminated by a "cusp" in the effective potential. The initial motion of the model tori was perturbed with uniform sub-sonic vertical and diagonal velocity fields. As the configuration evolved in time, we measured the mass accretion rate on the neutron star surface and obtained the power spectrum. The prominent mode of oscillation in the cusp torus is the radial epicyclic mode. It would appear that vertical oscillations are suppressed by the presence of the cusp. From our analysis it follows that the mass accretion rate carries a modulation imprint of the oscillating torus, and hence so does the boundary layer luminosity.

Figures

Tables

Table 1. Table 1: A list of models with model name, resolution, strength of the magnetic field and the initial profile of the magnetic field. The type of velocity perturbation is also listed.

Model	Resolution	$B_{\mathrm{NS}}$	$B_{\mathrm{NS}}$ profile	Perturbation
IHDV	420 $\times$ 300	$-$	$-$	vertical
IHDD	420 $\times$ 300	$-$	$-$	diagonal
IMHDD1	420 $\times$ 300	900 $G$	dipole	diagonal
IMHDD2	140 $\times$ 100	9000 $G$	dipole	diagonal

Equations

$${\Phi }_{\mathrm{KL}}=-\frac{GM}{3r_{\mathrm{s}}}\left(e^{3r_{\mathrm{s}}/r}-1\right),$$

$${\omega }_{\mathrm{R}}={\left(1-\frac{3r_{\mathrm{s}}}{r}\right)}^{1/2}{\Omega }_{\mathrm{K}},$$

Full text

††pubyear: 2017††pagerange: MHD simulations of oscillating cusp-filling tori around neutron stars – missing upper kHz QPO–MHD simulations of oscillating cusp-filling tori around neutron stars – missing upper kHz QPO

1 Introduction

Accreting neutron stars (NS) in low mass X-ray binaries (LMXBs) manifest quasi-periodic oscillations (QPOs), pairs of which appear at kHz frequencies (van der Klis et al., 1996; Strohmayer et al., 1996; van der Klis, 2000). There are different regions that produce NS X-ray emissions: the optically thick disk, the hot corona above the disk, and the equatorial belt of accretion on the stellar surface, the boundary layer (BL), if the NS magnetic dipole is not sufficiently strong to channel accretion onto poles. Phase resolved spectroscopy reveals that the X-ray luminosity in QPOs originate from the boundary layer on the neutron star surface (Gilfanov et al., 2003; Gilfanov & Revnivtsev, 2005). On theoretical grounds, low-frequency QPOs have also been attributed to the BL (Titarchuk & Osherovich, 1999), but in this contribution we only discuss the twin kHz QPOs.

Several neutron star models have radii smaller than the innermost stable circular orbit (ISCO), such that the disk is decoupled from the neutron star surface (Kluzniak & Wagoner, 1985). The accreted fluid releases a significant amount of energy in the boundary layer. It was shown that between 69$\%$ and 86$\%$ of the total luminosity is emitted from the boundary layer when matter falls from the ISCO onto the neutron star (Syunyaev & Shakura, 1986). The variable flow rate of the matter leaving the disk and crossing the “relativistic accretion gap” between the ISCO and the neutron star (Kluzniak & Wagoner, 1985) will lead to an intensified luminosity modulation as the matter enters the boundary layer (Paczynski, 1987). The frequencies of inner disk oscillations should thus lead to the variability of the boundary layer X-ray flux (Kluźniak & Abramowicz, 2005; Abramowicz et al., 2007).

The oscillations of neutron star accretion disks may modulate the X-ray flux by a general relativistic effect, the overflow of “Roche lobe” (or a cusp-filling torus) near the ISCO. The overflow is achieved by terminating the inner edge of the torus by a “cusp”, which is the point of self-crossing equipotential surface for any stationary rotating fluid with given angular momentum ($\ell_{\rm{c}}$). A cusp-filling torus accretes onto the star. The accretion is driven by gravity and pressure gradients, without need of viscous phenomena to transport the angular momentum (Paczynski, 1977; Abramowicz et al., 1978). An analytic model of mass accretion rate modulation as the fluid flows through the cusp of a torus oscillating in two axisymmetric eigenmodes has been constructed in a pseudo-Newtonian formulation (Abramowicz et al., 2007), and elaborated in general relativity for realistic equations of state to explain the twin peak QPOs observed in low mass X-ray neutron star binaries (Török et al., 2016).

The recently launched astronomy satellite ASTROSAT, equipped with the Large X-ray Proportional Counters (LAXPC), will provide better opportunities to investigate the QPOs in LMXBs with accreting neutron stars. The LAXPC has the largest effective area in the hard X-ray band and is capable of high time resolution X-ray measurements in the 2$-$80 keV band (Paul, 2013).

We have used the PLUTO code 111Freely available at http://plutocode.ph.unito.it/ (Mignone et al., 2007) to perform axisymmetric, grid-based ideal magnetohydrodynamic (MHD) simulations of oscillating cusp-filling tori around neutron stars. In order to mimic the effects of general relativity we used a pseudo-potential obtained by Kluźniak & Lee (2002, hereafter KL), which reproduces the ratio of orbital and epicyclic frequencies for the Schwarzschild metric. KL potential has been used by Parthasarathy et al. (2016) in the recent numerical investigation of QPOs in perturbed tori orbiting a non-rotating black hole.

The Letter is structured as follows. In Sect. 2 we explain the numerical set-up and present the results of our simulations. In Sect. 3 we discuss the results and summarize in Sect. 4.

2 Simulations and Results

We consider a torus accreting onto a canonical neutron star of mass 1.4 $M_{\odot}$ and radius 10 km. The numerical set-up and velocity perturbations in spherical coordinates ($r,\theta$) are adapted from Parthasarathy et al. (2016). The surface of the neutron star is at 2.4 $r_{\rm{s}}$, where the Schwarzschild radius is $r_{\rm{s}}=$ 2 $GM/c^{2}$. The outer radial domain is at 5 $r_{\rm{s}}$. The cusp of the torus is located at $r_{\rm{cusp}}=2.62r_{\rm{s}}$ and the radius of pressure maximum ($r_{\rm{c}}$) is at 3.47 $r_{\rm{s}}$ (see Fig. 1). The inner and outer co-latitudinal boundaries are located at $(\pi/2-0.3)$ and $(\pi/2+0.3)$, respectively. A description of the models is given in Tab. 1.

A radially increasing logarithmic grid is used, which preserves the cell aspect ratio. In such a grid the intrinsic characteristic spatial scale grows linearly with the distance from the stellar surface, exactly as the grid spacing of latitude. The grid resolution introduces about the same amount of numerical dissipation at all radii. Uniform grid is used in the latitudinal direction. The inner radial boundary is prescribed with a pseudo-reflective condition for modelling the neutron star surface, in which the components of the velocity and magnetic fields normal to the cell interface do not flip the signs.

In order to imitate the effects of general relativity we use the KL potential,

[TABLE]

where $r$ is the distance from the center of the neutron star. The KL potential reproduces the ratio of orbital and epicyclic frequencies for the Schwarzschild metric

[TABLE]

where $\omega_{\rm{R}}$ is the radial epicyclic frequency and $\Omega_{\rm{K}}$ = 2$\pi\nu_{\rm{K}}$ is the orbital frequency (see also Giussani et al. (2014); Parthasarathy & Kluźniak (2014)). The eigenfrequency of the radial epicyclic mode as a function of the equilibrium radial position of the torus ($r_{\rm{c}}$), is equal to 0.36 in units of Keplerian frequency ($\nu_{\rm{K}}$ at $r_{\rm{c}}$). In a spherically symmetric potential, the vertical epicyclic frequency is equal to $\Omega_{\rm{K}}$.

We set-up a stationary hydrodynamic cusp-filling torus around a non-rotating neutron star in the limit of vanishing magnetic dipole fields in the manner described hereafter. A constant angular momentum torus in equilibrium is initialised around the central compact body prescribed with KL potential (see Eq. 1). The inner radius of the torus coincides with the point of self-crossing equipotential surface ($r_{\rm cusp}$). The initial motion of the cusp-filling torus is unperturbed. There is a steady outflow of matter through the cusp of the stationary torus onto the star. The evolved quasi-equilibrium state of the accreting cusp-filling torus is used as the initial configuration for the hydrodynamic models (see Tab. 1). The motion of the evolved torus is triggered with uniform sub-sonic vertical and diagonal velocity fields to achieve the models IHDV and IHDD, respectively. The MHD models (see Tab. 1) include a weak NS dipole, and are perturbed with uniform sub-sonic diagonal velocity fields from the beginning.

We measured the mass accretion rate ($\dot{M}$) on the neutron star surface and obtained the power spectrum of $\dot{M}$ using the fast Fourier transform (FFT) technique. The power spectrum of $\dot{M}$ for the equilibrium model of the stationary overflow of Roche lobe does not manifest peaks that correspond to a periodicity. The results of models IHDV and IHDD are shown in Fig. 2, and of models IMHDD1 and IMHDD2 in Fig. 3. The most prominent mode excited with the vertical and diagonal perturbations (see Tab. 1) is the radial epicyclic mode of the cusp torus. The first harmonic of the radial mode is clearly visible in the power spectrum of the models IHDV (Fig. 2), IMHDD1 and IMHDD2 (Fig. 3). Additional frequencies are seen in the simulation of the model IHDV (Fig. 2).

3 Discussion

We numerically investigate the modulation of the neutron star boundary layer luminosity by constructing a torus with a constant distribution of angular momentum in the KL potential (Fig. 1). The fluid outflows from the torus through the cusp. The models of cusp-filling tori were initially perturbed with uniform sub-sonic vertical and diagonal velocity fields (Tab. 1). The results of our simulations are discussed below.

The model IHDV has been perturbed with uniform sub-sonic vertical velocity perturbation. The restoring force of a vertical oscillation is the vertical gradient of the gravitational potential. Since the cusp corresponds to the self-crossing of the equipotential surface, the gradient of the effective potential vanishes at the cusp. This seems to be suppressing vertical motion in the inner parts of the torus. The prominent mode of oscillation excited by the vertical motion of the model IHDV is the radial epicyclic mode of the torus (Fig. 2). In the model IHDV we observe the first harmonic of the radial mode and additional frequencies.

Abramowicz et al. (2007) analytically showed that the mass accretion rate retains the modulation imprint of the vertically oscillating disk, at double the mode frequency. Parthasarathy et al. (2016) numerically studied the frequencies of oscillation of a torus with an initial vertical motion, and demonstrated that the power of vertical motion in the density $L_{2}$ norm is concentrated in the first harmonic at twice the orbital frequency. In both the cases this is related to the fact that the equatorial plane is crossed twice per every vertical oscillation cycle. We expected to see a prominent vertical oscillation in the model IHDV. However, it was not observed in the $\dot{M}$ measured on the neutron star surface and the corresponding power spectrum of $\dot{M}$ (see model IHDV in Fig. 2). The uniform initial velocity perturbation does not excite the vertical eigenmode of the cusp-filling torus, since there is no restoring force at the cusp.

The models IHDD, IMHDD1 and IMHDD2 have been perturbed with uniform sub-sonic diagonal velocity perturbations. The diagonal perturbation may excite the radial and vertical eigenmodes (Parthasarathy et al., 2016; Mishra et al., 2017). We infer the radial epicyclic mode of the torus from the diagonal perturbations of the models IHDD (Fig. 2), IMHDD1 and IMHDD2 (Fig. 3), respectively. The vertical eigenmode was not excited in the simulations of models IHDD, IMHDD1 and IMHDD2. The first harmonic of the radial mode is present in the power spectrum of the models IMHDD1 and IMHDD2 (Fig. 3).

The results obtained from all the cusp-filling tori models are the same, although the initial conditions were not identical in the HD and MHD regimes (see Sec. 2). The most prominent peak in the power spectrum of vertically and diagonally perturbed tori (see Figs. 2 and  3) is the radial epicyclic mode frequency. The conclusion of our numerical analysis is robust in showing the modulation of the boundary layer luminosity by disk oscillations and a strong presence of radial epicyclic frequency in all models. However, the vertical epicyclic frequency is suppressed.

The variation of the quality factor and amplitude of the upper and lower kHz QPOs in a low-mass X-ray binary source was studied in Barret et al. (2005a); Barret et al. (2005b). They analysed the archival data from the Rossi X-ray Timing Explorer (RXTE). The quality factor gives a measure of the coherence of the oscillator. It was shown that the upper and lower QPOs follow different paths in the quality factor versus frequency diagram. The lower kHz QPOs exhibited much higher coherence (higher Q-factor) than the upper kHz QPOs. If we assume that in low-mass X-ray binary sources, the vertical eigenmode correspond to the upper kHz QPOs and the lower kHz QPOs to the radial epicyclic mode, our results seem to be in qualitative agreement with this finding. The upper kHz QPO is hard to discern in the boundary layer luminosity of our simulation. In contrast, the radial epicyclic frequency is a robust feature in the power spectrum of the mass outflow from the cusp-filling torus, appearing even when the initial perturbation is purely vertical.

4 Conclusions

Using the PLUTO code we performed axisymmetric ideal MHD simulations of oscillating cusp-filling tori around non-rotating neutron stars in the KL potential. For the first time we numerically study the modulation of the mass accretion rate on the neutron star surface. We draw the inference that, both the vertical and diagonal excitations of cusp-filling tori excite the radial epicyclic mode in such tori. Based on Parthasarathy et al. (2016), we expected that the vertical epicyclic mode may be excited, however the uniform velocity perturbations fail to excite the vertical eigenmode of the cusp-filling torus. Our simulations confirm the modulation mechanism of the neutron star boundary layer luminosity by disk oscillations (Abramowicz et al., 2007). However, the vertical eigenmode is suppressed as there is no restoring force at the cusp. The strong presence of the radial epicyclic frequency in the boundary layer luminosity (presumed to be the lower kHz QPO) and the absence of the vertical epicyclic frequency (presumed to be the upper kHz QPO) may be interpreted as corresponding to the high Q-factor of the lower kHz QPO (Barret et al., 2005a; Barret et al., 2005b) and the low Q-factor of the upper kHz QPO (Barret et al., 2005b).

The results of our numerical analysis are relevant to the ASTROSAT mission equipped with LAXPC (Paul, 2013) which will investigate the quasi-periodic oscillations in low mass X-ray binaries with low magnetic field neutron stars.

Acknowledgements

The authors are grateful to the anonymous referee for the helpful comments. This research was supported by the Polish NCN grant 2013/08/A/ST9/00795, and in part by the NSF under grant No. NSF PHY11-25915.