The MHD model of kHz QPOs

# The Magneto Hydro Dynamical Model of KHz Quasi Periodic Oscillations in Neutron Star Low Mass X-ray Binaries (II)

Chang-Sheng Shi11affiliation: National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100012, China 44affiliation: College of Material Science and Chemical Engineering, Hainan University, Hainan 570228, China 66affiliation: Key Laboratory of Modern Astronomy and Astrophysics (Nanjing University), Ministry of Education, Nanjing 210093, China , Shuang-Nan Zhang11affiliation: National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100012, China 22affiliation: Key Laboratory of Particle Astrophysics, Institute of High Energy Physics, Chinese Academy of Sciences, Beijing 100049, China; zhangsn@ihep.ac.cn 5 ⋆5 ⋆affiliationmark: , Xiang-Dong Li33affiliation: Department of Astronomy, Nanjing University, Nanjing 210093, China 66affiliation: Key Laboratory of Modern Astronomy and Astrophysics (Nanjing University), Ministry of Education, Nanjing 210093, China
###### Abstract

We study the kilohertz quasi-periodic oscillations (kHz QPOs) in neutron star low mass X-ray binaries (LMXBs) with a new magnetohydrodynamics (MHD) model, in which the compressed magnetosphere is considered. The previous MHD model (Shi & Li 2009) is re-examined and the relation between the frequencies of the kHz QPOs and the accretion rate in LMXBs is obtained. Our result agrees with the observations of six sources (4U 0614+09, 4U 1636–53, 4U 1608–52, 4U 1915–15, 4U 1728–34, XTE 1807–294) with measured spins. In this model the kHz QPOs originate from the MHD waves in the compressed magnetosphere. The single kHz QPOs and twin kHz QPOs are produced in two different parts of the accretion disk and the boundary is close to the corotation radius. The lower QPO frequency in a frequency-accretion rate diagram is cut off at low accretion rate and the twin kHz QPOs encounter a top ceiling at high accretion rate due to the restriction of innermost stable circular orbit.

accretion, accretion disks — magnetohydrodynamics — X-rays: binaries — stars: neutron
55affiliationtext: Physics Department, University of Alabama in Huntsville, Huntsville, AL 35899, USA

## 1 Introduction

The fastest variability, i.e., the high frequency quasi-periodic oscillations (QPOs), has been observed in both neutron star low mass X-ray binaries (NS-LMXBs) and black hole LMXBs, which gives us an important channel to understand the physics of accretion process in the accretion disks of LMXBs. Especially it gives us a clue to find the parameters of the compact stars, e.g. their equation of state (van der Klis et al. 2006). In NS-LMXBs the frequencies of the high frequency QPOs approach the Kepler frequency of the matter at the surface of the NS and often exceed 1000 Hz, so named as kilohertz (kHz) QPOs. The kHz QPOs always appear as broad peaks in their Fourier power spectra and were often discovered in pairs; they are labeled as upper QPOs () and lower QPOs () according to their frequencies. The kHz QPOs change with the source states and the X-ray count rates (Méndez et al. 1999) (probably the mass accretion rate). In addition there is an interesting phenomenon that the X-ray intensity and kHz QPOs do not exhibit a one-to-one relation in 4U 1608–52 and 4U 1636–53, i.e. “the parallel track” (Méndez 2003). Barret et al. (2005, 2006) discovered that the lower QPO frequency from 4U 1636–53 in the frequency-count rate diagram has a maximum value (around 920 Hz), above which frequency no QPO was detected for any count rate; they named this phenomenon the ceiling of the frequency. In this work, we simply call this maximum value the ceiling frequency.

Many models have been suggested to explain the physics of kHz QPOs in NS-LMXBs (see e.g. van der klis 2006, Shi & Li 2009). Lin et al. (2011) discussed the frequency relationship of kHz QPOs for 4U 1636–53 and Sco X–1 by comparing the observations to some theoretical models. Generally the models of the kHz QPOs mainly include several types, i.e. beat-frequency models (Miller et al. 1998; Cui 2000; Campana 2000; Lamb & Miller 2001, 2003; ), rotation, precession and epicyclic frequency models (Stella & Vietri 1999; Romanova & Kulkarni 2009; Bachetti et al. 2010), disk-oscillation and resonance models (Osherovich & Titarchuk 1999; Abramowicz & Kluzńiak W. 2001, 2003; Kato 2001, 2004; Urbanec et al. 2010), wave models (Zhang 2004; Li & Zhang 2005; Rezania & Samson 2005; Shi & Li 2009, 2010). The above classification is not strict, because some factors are overlapped in those studies. In most of the models in NS-LMXBs, the characteristic radius is very important, because the changing QPOs frequencies are very dependent on the radius in those models. The location of their production region of the kHz QPOs is discussed and several kinds of places have been suggested, e.g. the surface of the star (such as Romanova & Kulkarni 2009, Bachetti et al. 2010), the accretion disk (such as Shi & Li 2009, 2010, Sriram et al. 2011) or the magnetic field lines (such as Zhang 2004, Li & Zhang 2005).

Zhang (2004) considered the Alfvén waves as the source of the kHz QPOs of the accreting X-ray binaries. Rezania and Samson (2005) also discussed the MHD turbulence effect coming from the accretion process when the plasma hits the magnetosphere, i.e. the excited resonant shear Alfvén waves in a region of enhanced density gradients from collision in the magnetosphere led to the kHz QPOs. Shi & Li (2009) have considered the two magnetohydrodynamics (MHD) oscillation modes in NS magnetospheres including the effect of gravity and the rotation of a NS as the source of the twin kHz QPOs. A linear relation about the frequencies of the upper QPOs and the lower QPOs was obtained and the model fitted the observation about the change of the upper kHz QPO frequencies and the frequency difference with lower kHz QPOs frequencies well.

In this work, we re-examine the previous MHD model (Shi & Li 2009) and find the relation between the kHz QPOs and the accretion rate in NS-LMXBs, based on the interpretation of the MHD waves of Shi & Li (2009). Unlike the result of others, we find that the Alfvén-like transverse waves only exist under special conditions. We start in section 2 with the new MHD model and give the solutions of the dispersion equation on the MHD wave frequencies. Then, we compare our results with observations in section 3. Finally, we make discussion in section 4 and summarize our result in section 5.

## 2 The MHD model

In this section, we consider that a steady standard -disk of Shakura & Sunyaev (1973) in a NS-LMXB is truncated by the stellar magnetosphere, and the plasma is accreted to the polar cap along the magnetic field (e.g., Ghosh et al. 1977, Elsner & Lamb 1977). As shown in Fig. 1, in the accretion process, the plasma hits the magnetic field lines and compresses the primary polar magnetic field which might lead to a deformation of the magnetic field and some instability (Elsner & Lamb 1977). The MHD waves produced at the magnetosphere radius from a small perturbation lead to the kHz QPOs.

### 2.1 The MHD Waves

We start from a balance by gravity, barometric pressure, magnetic pressure for the plasma in the border between the magnetosphere and the steady thin accretion disk rotating around the NS; this gives a definition of the magnetosphere radius of the LMXBs. The balance equation can be expressed as follows,

 ρ0∂u0∂t+ρ0u0⋅∇u0=−∇P0+1μ(∇×B0)×B0+2ρ0u0×Ω+ρ0Ω2r0−ρ0(Ω⋅r0)Ω−ρ0GMr30r0, (1)

where is the plasma velocity, the vacuum magnetic conductivity, the barometric pressure, the magnetic field, the displacement from the NS, the plasma density, the gravitational constant, the mass of the NS and the angular velocity of the NS, respectively. The subscript “0” denotes variables in the equilibrium state and the bold italic expresses vectors. It is very different from Shi & Li (2009) that now we consider that the magnetic field and the density of the plasma in the magnetosphere are not uniform, but we consider the same initial conditions of Shi & Li (2009) that the balance is steady at first, i.e., , , and . We consider that the plasma is ideal conductor, so the vacuum electroconductivity , and then we can obtain the equation according to the Faraday principle of electromagnet induction as follows,

 ∂B0∂t=(B0⋅∇)u0−(u0⋅∇)B0−(∇⋅u0)B0. (2)

The continuity and the adiabatic condition are always used in the accretion process and they can be expressed as follows,

 ∂ρ0∂t+∇⋅(ρ0u0)=0, (3)
 P0ρ−γ0=const, (4)

where is the adiabatic index. The plasma at the magnetosphere radius is always disturbed by strong excitation (Rezania & Samson 2005), so now we discuss the MHD waves that are produced from a small disturbance. The disturbed MHD equations are written as,

 ρdudt+ρu⋅∇u=−∇P+J×B+2ρu×Ω+ρΩ×(r×Ω)−ρGMr3r, (5)
 ∂B∂t=∇×(u×B)=(B⋅∇)u−(u⋅∇)B−(∇⋅u)B, (6)
 ∂ρ∂t+∇⋅(ρu)=0, (7)
 Pρ−γ=const, (8)

where , , , , and the subscript “s” denotes the variation of a physical quantity due to the perturbation, except in the sound velocity () below. We consider that a small perturbation lead to the MHD waves, so , , , , . Combining Equations (1)(8), we can obtain the equations about the variation of the perturbed physical quantities in the first-order approximation in the corotation reference frame (so ):

 ρ0∂us∂t=−∇Ps+1μ[(∇×B0)×Bs+(∇×Bs)×B0]+2ρ0us×Ω+ρ0Ω2rs+ρsΩ2r0−ρ0(Ω⋅rs)Ω−ρs(Ω⋅r0)Ω−ρsGMr30r0−ρ0GM(rsr30−3r0⋅rsr0r50), (9)
 ∂Bs∂t=(B0⋅∇)us−(∇⋅us)B0−(us⋅∇)B0, (10)
 ∂ρs∂t=−∇ρ0⋅us−ρ0∇⋅us, (11)
 Ps=γP0ρ0ρs. (12)

Unlike Equation (13) in Shi & Li (2009), Equation (9) includes the previously neglected term that comes from , since sometimes ; in addition we also consider the nonuniform distribution of the density and magnetic field for the plasma in the magnetosphere in those equations.

Differentiating Equation (9) and substituting (12) into it, we obtain

 ρ0∂2us∂t2≃−∂∂t(∇γP0ρ0ρs)+1μ∂∂t[(∇×B0)×Bs+(∇×Bs)×B0]+2ρ0∂∂t(us×Ω)+ρ0Ω2us+∂∂tρsΩ2r0−ρ0(Ω⋅us)Ω−∂∂tρs(Ω⋅r0)Ω−∂∂tρsGMr30r0−ρ0GMr30(us−3r0⋅usr0r20). (13)

Now we discuss the physical process in the rectangular coordinate system () in the corotation reference frame (see Fig. 1). We assume that (i) the accretion disk does not warp, i.e. the compressed magnetic field lines are normal to the disk in the equatorial plane, the magnetic field and the density are only functions of the longitudinal displacement () after being compressed, i.e. , , ; (ii) the MHD waves propagate along the magnetic field lines or they would be dissipated in the disk easily (Shi & Li 2009), i.e. the wave vector can be expressed as, , ( is the wavenumber). The balance equation of the plasma at the magnetosphere radius, Equation (1), can be simplified as,

 0=−∇P0+1μ(∇×B0)×B0+ρ0Ω2r0−ρ0GMr30r0, (14)

and the further simplified form is,

 Ω2r0−GMr30r0=(c2sρ0∂ρ0∂r+1μ0ρ0B0∂B0∂r)∣r=r0, (15)

where is the sound velocity.

After substituting Equations (10)(12) and (15) into Equation (13), we can obtain,

 ρ0∂2us∂t2≃γP0ρ0[∇ρ0×(∇×us)+(∇ρ0⋅∇)us+(us⋅∇)(∇ρ0)+∇ρ0(∇⋅us)+ρ0∇(∇⋅us)]+1μ{∇B0×[(e⋅∇)us]+B0∇×[(e⋅∇)us]−∇(∇⋅us)×B0−(∇⋅us)∇×B0−∇×[(us⋅∇)B0]}×B0+1μ(∇×B0)×(B0⋅∇)us−1μ(∇×B0)×(∇⋅us)B0−1μ(∇×B0)×(us⋅∇)B0+ρ0Ω2us−ρ0(Ω⋅us)Ω−ρ0GMr30(us−3r0⋅usr0r20)+2ρ0∂∂t(us×Ω)−(∇ρ0⋅us+ρ0∇⋅us)£[Ω2r0−(Ω⋅r0)Ω−GMr30r0−∇(γρ0P0)], (16)

where is the unit vector that has the same direction with . We then carry out Fourier transformation () for Equation (16) and simplify the result, and obtain the following dispersion equations:

 (−ω2−Ω2−2ω2k + k2V2A+m)usx=−ik1μρ0B0∂B0∂xusz+ik(γ−1)c2sρ0∂ρ0∂xusz−i2ωΩusy, (17)
 (−ω2−Ω2+ω2k + k2V2A)usy=−i2ωΩusx, (18)
 (−ω2+ω2k + k2c2s)usz=(c2sρ0∂ρ0∂x+1μ0ρ0B0∂B0∂x)ikusx, (19)
 m=−[∂∂r(c2sρ0∂ρ0∂r+1μ0ρ0B0∂B0∂r)]∣r=r0, (20)

where is the Alfvén velocity and the Kepler angular frequency. If the magnetic field and the plasma are uniform, and (so ), this magnetosphere radius is also the corotation radius of the NS from Equation (15), and the frequencies of the MHD Alfvén-like waves are,

 ω2=k2V2A+12Ω2±12√Ω4+16k2V2AΩ2. (21)

Since the coefficients of the highest order derivatives in Eq. (16) are not small, we cannot make expansions for small coefficients nor assume that these coefficients are equal to 0. Therefore our approximation is a classical method when the total energy is greater than the potential energy in all space, compared with the WKB approximation method.

In the case , several solutions are obtained from Equations (17) (20).

1. If and , it is a sound-like wave solution, on condition that .

2. If and , it is a Alfvén-like wave solution, , on condition that the Coriolis force is neglected.

3. If and , i.e. the magnetosphere radius is also the corotation radius, we can obtain the MHD waves as, , that revert to Equation (21).

4. If and , we can obtain the dispersion equation from Equations (17) (20) as follows,

 (−ω2−Ω2−2ω2k + k2V2A+m)(−ω2−Ω2+ω2k + k2V2A)(−ω2+ω2k + k2c2s)−4Ω2ω2(−ω2+ω2k + k2c2s)=k2r20(Ω2−ω2k)2(−ω2−Ω2+ω2k + k2V2A)−(Ω2−ω2k)(−ω2−Ω2+ω2k + k2V2A)γk2c2sr0ρ0(∂ρ0∂r∣r=r0) (22)

The above solutions of the dispersion equations are special solutions under special conditions except the ones of Equation (22), so only the solutions of Equation (22) are ordinary ones and we consider them as the source of the kHz QPOs in NS-LMXBs. Because there is not a concise analytic solution for Equation (22), we will discuss the numerical solutions below.

The magnetosphere radius is the characteristic radius of the kHz QPOs in this model. Since Lamb et al. (1973) gave a definition about the magnetosphere radius in the globular accretion process onto a compact star, many authors have explored the outer boundary of the magnetosphere of a pulsar (e.g. Cui 1997) and they obtained not identical but analogous conclusions. The magnetosphere radius obtained by Lamb et al. (1973) is written as,

 rm≃2.29∗106μ4/726˙M−2/716m−1/71.4⊙ cm, (23)

where is the magnetosphere radius in the accretion process, the magnetic moment of the star, the accretion rate, the mass of the NS and the subscripts “26”, “16”, “” express the quantities in units of , , 1.4 times the mass of the sun, respectively.

McCray & Lamb (1976) considered that the magnetic pressure could resist the falling spherical plasma layer and they found a magnetosphere radius ranging from to . Elsner & Lamb (1977) continued to discuss the spherical accretion process and they considered that the central star rotates sufficiently slowly. In addition a lot of results about the magnetosphere radius close to were calculated if the quantities (, ) were adopted when many conditions, such as the radiation pressure, the different structures of the magnetic field were considered (Davidson & Ostriker 1973; Burnard et al. 1983; Mitra 1992; Li & Wang 1995; Wang 1996; Cui 1997; Weng & Zhang 2011; Shakura et al. 2012).

Long et al. (2005) considered the disk accretion onto the magnetized star with a dipole magnetic field and found by simulation that the magnetosphere radius is half the one from Elsner & Lamb (1977). Recently Kulkarni & Romanova (2013) made a three-dimensional MHD simulations of magnetospheric accretion at a quasi-equilibrium state, in which the gravitational, centrifugal and pressure gradient forces are in balance. Then a different dependence of the magnetosphere radius on the accretion rate of the LMXBs, the magnetic moment, the mass and the radius of the NS is found as,

 rm≈2.50×106μ2/526˙M−1/516m−1/101.4⊙R3/106 cm, (24)

where denotes the radius of the NS in units of . The more gradual changing trend with and comes from the compression of the magnetosphere by the disk matter, which leads to the non-dipole magnetic field of the external magnetosphere.

We also consider that a central NS with a dipole magnetic field, whose mass is , accretes plasma from the standard thin disk, in which the density of the plasma is,

 ρ=3.1×10−8α−7/10˙M11/2016m5/8⊙R−15/810f11/5 g cm−3, (25)

where , is the viscosity parameter (Frank et al. 2002).

We can obtain the barometric pressure and is temperature in the standard thin disk, where is Boltzmann constant and is mass of the proton. We consider the balance of the plasma by magnetic pressure, barometric pressure and collision, and the magnetosphere radius can be estimated according to, (Eslner & lamb 1977; Romanova et al. 2002). Due to in the standard thin disk when the radius ranges from to , the balance condition can be simplified as, . The magnetosphere radius can be estimated when the dipole magnetic field is adopted,

 rm≃6.64×105α4/15μ16/2726˙M−34/13516m−7/271.4⊙f−136/135 cm. (26)

We have also made the strict computation on the magnetosphere radius by the equation, , and the result is very close to the above magnetosphere radius.

If we select the characteristic value, i.e., , , and (King et al. 2007), all the magnetosphere radius can be simplified. The magnetosphere radius in the spherical accretion process can be simplified from Equation (23) as,

 rms≃2.29×106˙M−2/716B4/7∗8cm, (27)

where denotes the magnetic flux density at the surface of the NS in the unit of . After the parameters being substituted into Equation (26) the magnetosphere radius is simplified as,

 rmA≃3.60×105˙M−34/13516B16/27∗8f−136/135cm, (28)

which we call “model A” when the radius is used. We also simplify Equation (24) as,

 rmB≃2.50×106˙M−1/516B2/5∗8 cm, (29)

which we call “model B” when the compressed magnetosphere radius () is considered.

As shown in Figure 2, the magnetosphere radius in the spherical accretion process () is much larger and covers a larger range than the others. The magnetosphere radius in model A approaches the radius of a NS when is high enough or the magnetic field at the surface of NS is weak enough. The compressed magnetosphere radius () changes slower than the other two radii ( and ), which reveals the clear differences between the compressed magnetic field and the dipolar magnetic field. In this study we only discuss the accretion process of the plasma in accretion disk of LMXBs, so the magnetosphere radius of Equation (27) for spherical accretion will not be used hereafter. The innermost stable circular orbit (ISCO) is the smallest stable circular orbit of accretion plasma and the plasma will fall to the central NS after it passes through the orbit due to a dynamical instability for circular geodesics in general relativity. The ISCO of the accretion plasma around a NS with the mass () is marked as the dash lines in Figure 2 according to its radius (). Inside ISCO, the assumption of the standard disk model breaks down, so we assume that the disk is truncated at ISCO at high .

### 2.3 The general solutions in two kinds of magnetic configurations

We suppose that the MHD waves at the magnetosphere radius are the origin of the kHz QPOs, so numerical solutions of Equation (22) should be obtained according to the different parameters, , and for the characteristic value of NS-LMXBs, i.e., , , and (King et al. 2007). In the two models (i.e. model A and model B), two types of magnetic field configurations are adopted, distinguished by that the magnetic field is compressed or not.

In model A, the standard -disk contacts with the dipolar magnetic field at the magnetosphere radius and an equilibrium state is reached. We substitute the balance equation, (i.e. Equation (15)), into Equation (22) and obtain the equation when the dipolar magnetic field is considered,

 (−ω2−Ω2−2ω2k + k2V2A+m)(−ω2−Ω2+ω2k + k2V2A)(−ω2+ω2k + k2c2s)−4Ω2ω2(−ω2+ω2k + k2c2s)=(1−γ)k2r20(Ω2−ω2k)2(−ω2−Ω2+ω2k + k2V2A)−3γk2V2A(Ω2−ω2k)(−ω2−Ω2+ω2k + k2V2A). (30)

In Equation (30) the Kepler angular frequency, Alfvén velocity, sound velocity can be obtained as follows, , , , where .

In model B, the standard -disk compresses the magnetic field untill an equilibrium state is reached and the relevant physical quantities for Equation (22) are obtained as, , , where . Because the dipole magnetic field is compressed within the magnetosphere (), we adopt an approximation about the magnetic field structure that all the magnetic field lines, which are compressed between the magnetosphere radius () and the magnetosphere radius of the uncompressed magnetic field (), are uniform. During the compressing process the magnetic flux remains unchanged, so the magnetic field can be estimated from the equation,

 ∫rmBrmAB′02πrdr=∫∞rmABp2πrdr, (31)

where is the magnetic flux density of the compressed magnetic field and is the one of the uncompressed dipole magnetic field in the equatorial plane. The Alfvén velocity can be obtained as,

 VA=1.51∗108B11/8∗8˙M−37/8016(r2mB,6−r2mA,6)−1r−1mA,6f−11/10B cm/s. (32)

In order to derive the general solutions we should derive the secondary derivative of the density and the magnetic field (i.e. in Equation (20)) in model A and model B, we assume that the condition on the balance of plasma (i.e. Equation (15)) is extended to the vicinity of the magnetosphere radius, so . Besides that, we use the distribution of the density of the -disk () in Equation (22) in model B or the distribution of the dipolar magnetic field () in Equation (30) in model A. Finally we obtain the general numerical solutions after the two kinds of magnetosphere radii in Equations (28) and (29) are substituted into Equations (22) and (30), respectively (see Figure 3).

As shown in Figure 3, we can find that there are twin solutions (the upper solution and the lower solution ) and single solutions both in the two models; the solutions within can be considered as twin solutions and the solutions outside are single solutions in the left panel of Figure 3. In the left panel, the four characteristics of the solutions of Equation (30) in model A are listed as follows. (i) The solutions are divided into two parts, i.e. the left solutions corresponding to low and the right ones corresponding to high . (ii) A transition accretion rate corresponding to the transition radius () can be found by the vertical dash-dotted line; the transition radius is the border between the twin solutions and the single solution for an accretion rate and it is near the corotation radius. (iii) In the right part there are only twin solutions for an accretion rate. (iv) and increase with until they reach their ceiling frequencies at  Hz and  Hz, due to the restriction of ISCO.

In the right panel of Figure 3, the general real solutions of Equation (22) are obtained and their five characteristics are listed as follows. (i) The single solution described by the bold solid line is also the upper solution and supposed as the origin of the single kHz QPOs. (ii) and increase with until they reach their ceiling frequencies at  Hz and  Hz when the accretion disk is truncated by ISCO. (iii) There is a transition point (corresponding to a transition radius ) that splits the solutions into the twin solutions and the single one; the lower real solutions will disappear when is lower than the transition accretion rate. (iv) The turning point of the single solutions (corresponding to a transition radius ) separates the decreasing trend from the increasing trend of and means the changing of the key factor dominating the balance of the plasma. (v) All the two transition radii ( and ) are close to the corotation radius at which the Kepler rotation frequency of the plasma equals to the spin of a NS.

The transition accretion rates in the two models (corresponding to and ) all separate the twin solutions from the single solutions and they originate from the same reason: with the decrease of the quantity () decreases from a positive value into a negative value and meanwhile the twin solutions change into a single solution, i.e. the lower real solutions disappear. According to the transition condition, , and the balance Equation (15), we can obtain a positive expression . In model A, of the dipolar magnetic field leads to , i.e., a compressed dense plasma layer. Conversely in model B, of the standard disk density leads to , i.e., a compressed stronger magnetic field topology.

The changing trend of the single solutions in model B is different from the one in model A; the main factor can be found from Figure 4 in model B and Equation (22), as follows. (i) The Kepler angular frequency is very high ( & ) and it is the key factor to determine the solutions of Equation (22). (ii) With the increase of , decreases and the other variables (, ) in Equation (22) increase. (iii) The expression is smaller than except a small section near the turning point (see panel E of Figure 4), so Equation (22) can be simplified as,

 (ω2+2Ω2+4ω2k)(−ω2−Ω2+ω2k)(−ω2+ω2k)+4Ω2ω2(−ω2+ω2k)≃k2r20(Ω2−ω2k)2(ω2+Ω2−ω2k). (33)

(iv) In Equation (33), it can be inferred that the expression dominates the changing trend because all the other expressions keep the changing trend with increasing . The changing trend of is shown in panel F of Figure 4 and its transition accretion rate is also close to that of the single solution.

We can draw an conclusion that the frequencies of the solutions of Equation (22) will change with the increase of the accretion rate in the same trend with .

The negative factor () in Equation (30) in model A from the balance Equation (15) leads to a different changing trend of the frequencies of the single solutions with the accretion rate. The configuration of the magnetic field is related to the distribution of the density by the balance equation. Finally we can conclude that the differences of the configuration of the magnetic field and the distribution of the density lead to the different changing trends of the frequencies of the solutions in model A and B.

## 3 Comparison with observations

### 3.1 Twin kHz QPOs

We take several steps to select the suitable parameters, such as the wavenumber, the surface magnetic field of NSs to match the observations.

1. Choose the initial value of the magnetic flux density (e.g. ) and the accretion rate (e.g. ), and then substitute them into or , and the detailed expressions on , .

2. The solutions of Equation (22) or (30) can be obtained when different wavenumbers, , are considered.

3. Select the value of for which the model predicted frequency difference () is closest to the discovered average value of (such as 322 Hz for 4U 0614+09).

4. The solutions of Equation (22) or (30) can be obtained again with the above , the same accretion rate (e.g. ) and different values of .

5. Select the value of , for which the twin solutions are closest to the observed frequencies of the twin kHz QPOs.

6. After substituting the above , and different values of into Equation (22) or (30), we can obtain new solutions.

7. Change the value of slightly and and repeat step 6 until the solutions match the frequencies of most observed twin kHz QPOs.

8. Change the value of slightly and repeat step 6 until the solutions match the frequencies of most observed twin kHz QPOs.

9. Repeat the above steps (7) and (8) in turn until the numerical solutions match the maximum number of the observed twin kHz QPOs.

10. With the last and , we can obtain both and changing with in the two models.

11. With all the parameters determined above (and listed in Table 1), We compute by requiring .

12. For each , we numerically calculate .

13. Finally in Figure 5, we plot the numerically calculated and observed twin kHz QPOs as functions of .

Because the analytic solutions are too complex and their expressions are too long to be used, we just choose the “best” parameters that match the most data as much as possible by eye-balling. As shown in Table 1, most of the wavenumbers () are close to in Table 1. Rezania and Samson (2005) discussed that the wavelength of the MHD wave in LMXBs was in the order of magnitude of the radius of NS.

The data points describe the observational data and the solid lines come from our numerical solutions in Figure 5. In the left panels of Figure 5, we consider the dipole magnetic field and use the magnetosphere radius in model A from Equation (28). In , suitable values of for several groups of observed QPOs are not found with model A and thus those QPOs are not plotted. The right panels describe the result while the compressed magnetic field is considered and the magnetosphere radius from Equation (29) is used.

In the left panels of Figure 5 we could not find a group of parameters for 4U 0614+09, 4U 1608–52, 4U 1636–53 in order to match the observations and we can only adopt not one but two different () for 4U 1636–53 for the best result to match the observations including all the twin kHz QPOs. The relation between and is reproduced well in the right panels of Figure 5 and the result from model B is much better than that from model A.

In the right panels of Figure 5, our numerical solutions in 4U 1608–52 and 4U 1636–53 deviate the observation slightly in the tails of the curves maybe due to ISCO. Barret et al. (2005, 2006) discovered the ceiling of the lower QPO frequency in 4U 1636–53 and 4U 1608–52 in a frequency-count rate diagram. As shown in Figure 5, however the ceiling of in 4U 1636–53 and 4U 1608–52 seem to be clearer than . This is is different from Barret et al. (2005) and the main reason is that in plotting our results in the figure, is required to determine and so deviations can only exist in . Then the magnetosphere radius can be determined with Equation (29).

We can estimate the masses of the two NSs in 4U 1636–53 and 4U 1608–52 if the magnetosphere radius corresponding to the ceiling frequency is ISCO (see Figure 3), i.e., , where is the ISCO of a NS with and in Equation (29) when the QPO frequency reach the ceiling frequency. Then the ceiling frequencies and their errors can be found as follows (see Figure 6):

1. The histogram of the lower (or upper) QPO frequencies is plotted with the bin size of 20 Hz (slightly larger than the error (about 17 Hz) of the observational data);

2. We look for a peak in the histogram near the high frequency end;

3. We pick all the QPO frequencies if their absolute differences between the peak frequency are less than 40 Hz;

4. The average of those QPO frequencies is considered as the ceiling frequency and their root-mean-square (RMS) is the error of the ceiling frequency.

With the ceiling frequency of either or , we can obtain from the numerical solution for model B. Then with (see Table 1) determined from the fitting (see section 3.1 for details), can be determined. In Figure 7, we show the relations between the ceiling frequencies as functions of .

For 4U 1636–53 we have and , corresponding to and , respectively. In comparison, was obtained by Kaaret et al. (1997) by assuming the maximum maximum as the Kepler frequency at ISCO in a Kerr spacetime. Similarly, in 4U 1608–52 or for or , respectively.

### 3.2 Single kHz QPOs

The corresponding for single kHz QPOs cannot be identified by the above method for the twin kHz QPOs and can be considered as our model prediction in the right panels of Figure 5. In the left panels of Figure 5, the single numerical solutions in model A are too small to match the observed single kHz QPOs so we will discuss the single kHz QPOs only in model B.

Generally we use a Lorentzian function to describe the finite-width peak in the power spectrum, i.e. QPOs, and the QPOs are confirmed only if the qualify factor is larger than 2. As shown in Figure 8, the single kHz QPOs with low count rate outside the box are related to low qualify factor and the other one of the twin kHz QPOs may also be related to a low qualify factor. The other kHz QPOs with very low qualify factor may be missed and some reported single kHz QPOs may be one of the twin kHz QPOs.

Barret et al. (2006) believed that the origin of was different from the one of due to the different changing trends of their qualify factors. It is also possible that the qualify factors is related to the different origins of twin kHz QPOs and single kHz QPOs. As shown in the right panel of Figure 8 for 4U 1728–34, three reported “single” kHz QPOs clustered around the twin kHz QPOs probably belong to twin kHz QPOs, but with one of the twin kHz QPOs missed from detection due to possibly lower signal to noise ratio. On the other hand, the reported single kHz QPOs in the box of each of the two panels of Figure 8 are located quite distinctively separated from all others, and are thus considered true single kHz QPOs; their frequencies decrease with the increase of the count rate (probably an indicator of accretion rate), as predicted by model B shown in Figure 5. The real single kHz QPOs in the other sources are identified by the same method as above and the unselected single kHz QPOs listed in the references in the caption of Figure 5 are omitted from Figure 5.

## 4 Discussion

In this study we only consider that the rotation axis of NSs coincides with the magnetic axis, i.e. the inclination angle is zero. Méndez et al. (2001) concluded that the properties of the kHz QPOs are determined only by the mass accretion rate through the disk. Lamb et al. (2009) considered that the magnetic inclination is likely to be very small for accretion-powered millisecond pulsars. Romanova & Kulkarni (2009) found that the moving spots in the magnetic boundary layer regime might produce QPO features in some cases. The further simulation by Kulkarni & Romanova (2013) showed that the magnetosphere radius mainly depends on the accretion rate in Equation (24) but not on the misalignment angle of the dipole magnetic field (also see Figure 4 in paper of Kulkarni & Romanova 2013). It seems that there is little influence on the kHz QPOs from the magnetic inclination because the central frequencies of kHz QPOs in our model are mainly determined by the magnetosphere radius.

With the decrease of the accretion rate the frequencies of the single kHz QPOs increase outside and will exceed the ceiling of the twin kHz QPOs in model B when the accretion rate of LMXBs is very low (about for the parameters in the right panel of Figure 3), however the expected high frequency cannot be detected due to the following reasons. (i) A NS-LMXB in such a very low accretion rate cannot be detected with sufficiently high signal to noise ratio. (ii) we can confirm kHz QPOs only when the fluctuation of the signal is small enough and signal-noise ratio is high enough. (iii) The magnetosphere radius was not simulated by Kulkarni et al. (2013) when the the accretion rate is very low and maybe it can be considered as a uncompressed one, i.e the modes of the kHz QPOs in the compressed magnetic field may convert to the ones in the dipolar magnetic field.

In our model the surface magnetic field of the NS is considered as an invariant and so the accretion rate of the LMXBs that determines the magnetosphere radius is the key parameter. In Figures 3 and 5 the accretion rate is not an observed result because it is determined by our selection from comparing the central frequencies of the observational lower kHz QPOs to our lower numerical solution. Generally the energy spectrum from each observation needs to be fitted to calculate the corresponding physical accretion rate. However the quality of the available data from RXTE is not good enough to do that. This is why we compare in the twin kHz QPOs with the observation by means of model derived accretion rate and the detailed relation between the frequencies of kHz QPOs and the measured accretion rate needs to be tested with future observations of better X-ray instruments.

## 5 Summary

In this study we re-examined the MHD model of kHz QPOs and the relation between the frequencies of the kHz QPOs and the accretion rate is predicted in the two models. In model A, the magnetic field of a NS keeps its dipolar topology and the accretion disk is compressed due to the magnetic pressure of its dipolar field. In model B, the accretion disk keeps the standard -disk and the magnetosphere of a NS is compressed due to the gas pressure of the standard disk. We find that the results of model B match the observations much better. Our main results are summarized as follows.

1. The Alfvén-like wave and the sound-like wave at the magnetosphere radius only exist under special conditions.

2. We predict that the accretion process of NS-LMXBs may happen in two different areas: sometimes only single kHz QPOs are produced and sometimes twin kHz QPOs are produced. There is a transition radius for the changing kHz QPOs from the single kHz QPOs to the twin kHz QPOs in frequency-accretion rate diagrams, i.e. in model A or in model B as shown in Fig. 3.

3. In model B, the frequency of the single kHz QPOs decreases first and then increases with increasing accretion rate; the transition radius is , as shown in Fig. 3. All the frequencies of the twin kHz QPOs increase with the increase of the accretion rate.

4. The transition radii (, , ) are all near the corotation radius.

5. In model B, the lower QPO frequency in a frequency-accretion rate diagram is cut off at low accretion rate and the twin kHz QPOs encounter a top ceiling at high accretion rate due to the restriction of ISCO.

6. In model B, the mass of the NS in 4U 1636–53 is estimated as and the mass of the NS in 4U 1608–52 is estimated as about , provided that the observed ceiling frequencies of each NS originate from the magnetosphere radius truncated at the ISCO of the NS.

This work is supported by 973 Program of China under grant 2014CB845802, by the National Natural Science Foundation of China under grant Nos. 11133002, 11373036, 11203009, 11133001 and 11333004, by the Strategic Priority Research Program ”The Emergence of Cosmological Structures” of the Chinese Academy of Sciences, Grant No. XDB09000000, by the Qianren start-up grant 292012312D1117210, and by the 54th postdoctoral Science Foundation funding of China.

## References

• Abramowicz & Kluzńiak (2001) Abramowicz, M. A., & Kluzńiak, W., 2001, A&A, 374, 19
• Abramowicz et al. (2003) Abramowicz, M. A., Bulik, T., Bursa, M., Kluzńiak, W., 2003, A&A, 404, L21
• Altamirano et al. (2008) Altamirano, D., van der Klis, M., Méndez, M., Jonker, P. G., Klein-Wolt, M., Lewin, W. H. G., 2008, ApJ, 685, 436
• Arons (1993) Arons, J., 1993, ApJ, 408, 160
• Bachetti et al. (2010) Bachetti, M., Romanova, M. M., Kulkarni, A., Burderi, L., di Salvo, T., 2010, MNRAS, 403, 1193
• Barret et al. (2005) Barret, D., Olive, J. F., Miller, M. C., 2005, AN, 326, 808
• Barret et al. (2006) Barret, D., Olive, J. F., Miller, M. C., 2006, MNRAS, 370, 1140
• Boirin et al. (2000) Boirin, L., Barret, D., Olive, J. F., Bloser, P. F., Grindlay, J. E., 2000, A&A, 361, 121
• Burnard et al. (1983) Burnard, D. J., Arons, J., Lea, S. M., 1983, ApJ, 266, 175
• Campana (2000) Campana, S. 2000, ApJ, 534, L79
• Cui (1997) Cui, W., 1997, ApJ, 482, L163
• Cui (2000) Cui, W., 2000, ApJ, 534, L31
• Cumming et al. (2001) Cumming, A., Zweibel, E., Bildsten, L., 2001, ApJ, 557, 958
• Davidson & Ostriker (1973) Davidson, K., & Ostriker, J. P., 1973, ApJ, 179, 585
• Di Salvo et al. (2001) Di Salvo T., Méndez M., van der Klis M., Ford E., Robba N. R., 2001, ApJ, 546, 1107
• Di Salvo et al. (2003) Di Salvo T., Méndez M., van der Klis M., 2003, A&A, 406, 177
• Elsner & Lamb (1977) Elsner, R. F., & Lamb F. k., ApJ, 1977, 215, 897
• Frank et al. (2002) Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics: Third Edition (Cambridge University Press)
• Ghosh et al. (1977) Ghosh, P., Lamb, F. K., Pethick, C. J., 1977, ApJ, 217, 578
• Ghosh & Lamb (1978) Ghosh, P., & Lamb, F. K. 1978, ApJ, 223, L83
• Jonker et al. (2000) Jonker, P. G., Méndez, M., & van der Klis M., 2000, ApJ, 540,29
• Jonker et al. (2002) Jonker, P. G., Méndez, M., & van der Klis M., 2002, MNRAS, 336, L1
• Kaaret (1997) Kaaret, P., Ford, E., Chen, K., 1997, ApJ, 480, L27
• Kato (2001) Kato, S., 2001, PASJ, 53, 37
• Kato (2004) Kato, S., 2004, PASJ, 56, 905
• Kulkarni & Romanova (2013) Kulkarni, A. K., & Romanova, M. M., 2013, MNRAS, 433, 3048
• Lamb et al. (1973) Lamb, F. K., Pethick, C. J., Pines, D. 1973, ApJ, 184, 271
• Lamb & Miller (2001) Lamb, F. K., & Miller, M. C., 2001, ApJ, 554, 1210
• Lamb et al. (2009) Lamb, F. K., Boutloukos, S., Van Wassenhove, S., Chamberlain, R. T., Lo, K. H., Clare, A., Yu, W. F., Miller, M. C., 2008, ApJ, 706, 417
• Lovelace et al. (2005) Lovelace, R. V. E., Romanova, M. M., Bisnovatyi-Kogan, G. S., 2005, ApJ, 625, 957
• Li & Wang (1995) Li, X. D., & Wang, Z. R., 1995, AP&SS, 225, 47
• Li & Zhang (2005) Li, X. D., & Zhang, C. M., 2005, ApJ, 635, L57
• Lin et al. (2011) Lin, Y. F., Boutelier, M., Barret, D., Zhang, S. N., 2011, ApJ, 726, 74
• Linares et al. (2005) Linares, M., van der Klis, M., Altamirano, D., Markwardt, C. B., 2005, ApJ, 634, 1250
• McCray & Lamb (1976) McCray, R., & Lamb, F. K., 1976, ApJ, 204, 115
• Méndez et al. (2001) Méndez, M., van der Klis, M., Ford, E. C., 2001, ApJ, 561, 1016.
• Méndez (2003) Méndez, M., 2003, ASPC, 308, 289
• Migliari et al. (2003) Migliari S., van der Klis M., Fender R. P., 2003, MNRAS, 345, L35
• Miller et al. (1998) Miller, M. C., Lamb, F. K., Psaltis, D., 1998, ApJ, 508, 791
• Mitra (1992) Mitra, A., 1992, A&A, 257, 807
• Osherovich & Titarchuk (1999) Osherovich, V. & Titarchuk, L. 1999, ApJ, 522, L113
• Pringle & Rees (1972) Pringle, J. E. & Rees, M. J., 1972, A&A ,21, 1
• Rappaportet al. (2004) Rappaport, S. A., Fregeau, J. M., Spruit, H., 2004, ApJ, 2004, 606, 436
• Rezania & Samson (2005) Rezania, V., & Samson, J. C., 2005, A&A, 436, 999
• Romanova et al. (2002) Romanova, M. M., Ustyugova, G. V., Koldoba, A. V., & Lovelace, R. V. E., 2002, ApJ, 578, 420
• Romanova & Kulkarni (2009) Romanova, M. M., & Kulkarni, A. K., 2009, MNRAS, 398, 1105
• Spruit & Taam (1993) Spruit, H. C., & Taam, R. E., 1993, Apj, 402, 593
• Stella & Vietri (1998) Stella, L., & Vietri, M. 1998, ApJ, 492, L59
• Strohmayer et al. (1996) Strohmayer, T. E., Zhang, W., Swank, J. H., Smale, A., Titarchuk, L., Day, C., Lee, U., 1996, ApJ , 469, 9
• Shakura et al. (2012) Shakura, N., Postnov, K., Kochetkova, A., Hjalmarsdotter, L., 2012, MNRAS, 420, 216
• Shi & Li (2009) Shi, C. S., & Li, X. D., 2009, MNRAS, 392, 264
• Shi & Li (2010) Shi, C. S., & Li, X. D., 2010, ApJ, 714, 1227
• Sriram et al. (2011) Sriram, K., Rao, A. R., Choi, C. S., 2011, ApJ, 743, 31
• Urbanec et al. (2010) Urbanec, M., Török, G., Šrámková, E., Čech, P., Stuchlík, Z., Bakala, P., 2010, A&A, 522, 72
• van der Klis et al. (2006) van der Klis M., 2006, in Lewin W. H. G., van der Klis M., eds, Compact Stellar X-ray Sources. Cambridge Univ. Press, Cambridge, p. 39
• van Straaten et al. (2000) van Straaten, S., Ford, E. C., van der Klis, M., Méndez, M., Kaaret, P., 2000, ApJ, 540, 1049
• van Straaten et al. (2002) van Straaten, S., van der Klis, M., di Salvo, T., Belloni, T., Psaltis, D., 2002, ApJ, 568, 912
• van Straaten et al. (2003) van Straaten, S., van der Klis, M., Méndez, M., 2003, ApJ, 596, 1155
• Wang (1996) Wang, Y. M., 1996, ApJ, 465, 111
• Weng & Zhang (2011) Weng, S. S., & Zhang, S. N., 2011, ApJ, 739, 42
• Wijnands et al. (1997) Wijnands, R. A. D., van der Klis, M., van paradijs, J., Lewin, W. H. G., Lamb, F. K., Vaughan, B., Kuulkers, E., 1997, ApJ, 479, L141
• Zhang (2004) Zhang, C. M. , 2004, A&A, 423, 401.
• Zhang et al. (2006) Zhang, F., Qu, J., Zhang, C. M., Chen, W., Li, T. P., 2006, ApJ, 646, 1116