On the Bimodal Spin-period Distribution of Be/X-Ray Pulsars

Be/X-ray binaries (BeXBs) are a subclass of high-mass X-ray binaries, usually composed of a neutron star (NS) and an O/Be type star (Reig 2011). The O/Be stars in BeXBs are generally fast-rotating main-sequence stars with a circumstellar disk, featured by Balmer emission lines. These observational properties can be well reproduced by the decretion disk model (Lee et al. 1991), which is similar to an accretion disk except for the inverse direction of mass flow. Note that the circumstellar disk is not a stable structure, undergoing dissipation and recurrence (Haigh et al. 1999), the size and the structure of which are strongly affected by the gravity and the orbital motion of the NS (Negueruela & Okazaki 2001; Negueruela et al. 2001; Okazaki et al. 2002). Part of the disk material is captured by the NS, preferentially near periastron. The accreted material is then channeled onto the NS surface by its magnetic field lines, so the NS usually appears as an Be/X-ray pulsar (BeXP; Lamb et al. 1973).

BeXBs are further divided into persistent and transient sources. A small population of BeXBs are persistent sources, characterized by nearly circular orbits (with eccentricities e ≲ 0.2) and low X-ray luminosities (L_X ∼ 10³⁴–10³⁵ erg s⁻¹), while other BeXBs are transient sources with eccentric orbits (with e ≳ 0.3; Reig 2011). They experience outbursts from time to time, which are separated by long, quiescent intervals. There are two types of outbursts in transient BeXBs. Type I or normal outbursts are characterized by (quasi-)periodicity with peak luminosities ranging from ∼10³⁶ to ∼10³⁷ erg s⁻¹. They are thought to be caused by rapid accretion of the NS when it passes periastron. Type II or giant outbursts show nonperiodicity, and usually have peak luminosities ≳10³⁷ erg s⁻¹. The origin of type II outbursts is still open to debate.

Knigge et al. (2011) found that the distribution of the spin periods (P) of BeXPs in the Milky Way (MW) and the Small/Large Magellanic Clouds (SMC/LMC) show a bimodal feature, peaked at ∼10 s and ∼200 s, respectively. These authors proposed that the two populations may be related to two types of supernovae, i.e., iron core-collapse supernovae and electron-capture supernovae, respectively. Cheng et al. (2014) suggested that the two populations of BeXBs are likely associated with different accretion modes in the two types of outbursts. Because the peak luminosity and the accretion disk structure are different in type I and II outbursts, NSs undergoing different types of outbursts are subject to different accretion torques, which lead to bimodal distribution of the spin periods.

The accretion process in BeXBs is very complicated and uncertain. Although it is generally believed that the NS undergoes disk accretion during outbursts, scenarios involving wind accretion have also been proposed (Ikhsanov 2007; Shakura et al. 2012). Even for disk accretion, there were arguments that the viscous timescale of the standard thin disk (Shakura & Sunyaev 1973) fails to explain the observed outburst duration of type I outbursts, suggesting an advection-dominated accretion flow (ADAF) instead (Okazaki et al. 2013).

Several authors (Chashkina & Popov 2012; Klus et al. 2014; Ho et al. 2014) investigated the long-term spin evolution of BeXPs in the SMC and arrived at the conclusion that a large population of the BeXPs are magnetars, i.e., NSs with surface magnetic fields stronger than 4 × 10¹³ G. But they are in conflict with the observations of cyclotron lines in the Galactic X-ray pulsars (Coburn et al. 2002; Caballero & Wilms 2012; Walter et al. 2015, and references therein).

Considering the fact that most BeXBs are transient sources, the transition between outbursts and quiescence should play an important role in the spin evolution of the NSs. In most of the previous studies, only the long-term average evolution was taken into account, which actually deviates from the evolution in the case of transient accretion. In this work, we attempt to investigate the spin evolution in BeXPs in a more self-consistent way by considering accretion in both outbursting and quiescent stages. We take into account the possible influence of various types of outbursts and the structure of the accretion flows. In Section 2, we describe our theoretical considerations. In Section 3, we describe the data sample and obtain possible constraints on the model parameters through comparison with observations. We summarize our work in Section 4.

2.1. Spin Evolution

We assume that the NS in a BeXB is interacting with an accretion disk originating from the Be star's wind. We first introduce the light-cylinder radius R_L and the inner radius of the accretion disk R₀. Here, R_L is defined to be

$$\begin{eqnarray}&&{R}_{{\rm{L}}}=\displaystyle \frac{c}{{{\rm{\Omega }}}_{{\rm{S}}}},\end{eqnarray} \tag{ 1 }$$

where c is the speed of light and Ω_S is the angular velocity of NS; R₀ is evaluated by

$$\begin{eqnarray}&&{R}_{0}=\phi {R}_{{\rm{A}}},\end{eqnarray} \tag{ 2 }$$

where ϕ is a factor around unity (Ghosh & Lamb 1979a, 1979b; Long et al. 2005), taken to be 0.5 in our work, and R_A is the Alfvén radius for spherical accretion

$$\begin{eqnarray}&&{R}_{{\rm{A}}}={\left(\displaystyle \frac{{\mu }^{2}}{\dot{M}\sqrt{2{GM}}}\right)}^{2/7},\end{eqnarray} \tag{ 3 }$$

where μ is the magnetic moment of the NS, $\dot{M}$ is the accretion rate, G is the gravitational constant, and M is the mass of the NS, taken to be 1.4 M_☉. According to the relation between R_L and R₀, the spin evolution is divided into two states: the rotation-powered state (R_L < R₀) and the accretion-powered state (R_L > R₀).³

In the rotation-powered state, the spin evolution of the NS is determined by the torque N_dip due to magnetic dipole radiation (Shapiro & Teukolsky 1983), i.e.,

$$\begin{eqnarray}&&{N}_{\mathrm{dip}}=-\displaystyle \frac{{\mu }^{2}{{\rm{\Omega }}}_{{\rm{S}}}^{3}}{9{c}^{3}}.\end{eqnarray} \tag{ 4 }$$

Considering the condition R_L = R₀, we obtain the critical spin period that separates the rotation- and accretion-powered states,

$$\begin{eqnarray}&&{P}_{\mathrm{cr}}\approx 0.23{\left(\displaystyle \frac{\mu }{{10}^{30}{\rm{G}}{\mathrm{cm}}^{3}}\right)}^{4/7}{\left(\displaystyle \frac{\dot{M}}{{10}^{14}{\rm{g}}{{\rm{s}}}^{-1}}\right)}^{-2/7}\,{\rm{s}}.\end{eqnarray} \tag{ 5 }$$

The critical spin period is much smaller than the current spin periods of most BeXPs with typical magnetic fields and accretion rates. Therefore, the rotation-powered state is usually ignorable.

In the accretion-powered state, we assume that the spin evolution of the NS is dominated by disk accretion and the interaction between the disk and the NS magnetic field, according to the magnetically threaded accretion disk model (Ghosh & Lamb 1979a, 1979b). The accretion torque N_d exerted by the disk consists of two terms,

$$\begin{eqnarray}&&{N}_{{\rm{d}}}={N}_{0}+{N}_{{\rm{M}}}.\end{eqnarray} \tag{ 6 }$$

Here N₀ is the material torque due to mass accretion,

$$\begin{eqnarray}&&{N}_{0}=\dot{M}{{\rm{\Omega }}}_{{\rm{K}}}({R}_{0}){R}_{0}^{2},\end{eqnarray} \tag{ 7 }$$

where Ω_K is the Keplerian angular velocity in the disk, and N_M is the torque resulting from the magnetic field-disk interaction, due to the shearing motion between the magnetic field lines of the NS and the disk material. Equation (6) can be expressed in the following form (Ghosh & Lamb 1979a, 1979b),

$$\begin{eqnarray}&&{N}_{{\rm{d}}}={N}_{0}n(\omega ),\end{eqnarray} \tag{ 8 }$$

where ω is the fastness parameter defined by

$$\begin{eqnarray}&&\omega \equiv \displaystyle \frac{{{\rm{\Omega }}}_{{\rm{S}}}}{{{\rm{\Omega }}}_{{\rm{K}}}({R}_{0})},\end{eqnarray} \tag{ 9 }$$

and n(ω) is a dimensionless function, which is, according to Ghosh & Lamb (1979a, 1979b),

$$\begin{eqnarray}&&n\approx 1.39\left\{1-\omega \left[4.03{\left(1-\omega \right)}^{0.173}-0.878\right]\right\}{\left(1-\omega \right)}^{-1}.\end{eqnarray} \tag{ 10 }$$

This equation was derived based on the assumption that the NS magnetic field penetrates the disk via the Kelvin–Helmholtz instability, turbulent diffusion, and reconnection with small-scale fields within the disk. However, Wang (1987) pointed out that this will lead to the pressure of the wound field to exceed the thermal pressure at large radius and disrupt the disk. There are various forms of the n(ω) function proposed in the literature (e.g., Wang 1995; Li & Wang 1996; Kluźniak & Rappaport 2007). Here we adopt the following simplified equation

$$\begin{eqnarray}&&n\approx 1-\displaystyle \frac{\omega }{{\omega }_{{\rm{c}}}},\end{eqnarray} \tag{ 11 }$$

which is similar to the torque equation in Lipunov (1982) and used in previous studies (e.g., Li 1999; Hartmann 2002; Ho et al. 2014; Ekşi et al. 2015). In Equation (11), ω_c is a critical value of the fastness parameter at which n = 0. This means that N_d > 0 when ω < ω_c and N_d < 0 when ω > ω_c. Also note that Equation (11) can be used to cover both accretor and propeller states: when $\omega \to 0$, ${N}_{{\rm{d}}}=\dot{M}{R}_{0}^{2}{{\rm{\Omega }}}_{{\rm{K}}}({R}_{0});$ when ω ≫ 1, ${N}_{{\rm{d}}}\sim -\dot{M}{R}_{0}^{2}{{\rm{\Omega }}}_{{\rm{S}}}$.

It should be pointed out that Equation (8) is only applicable to steady accretion onto NSs, while the accretion processes in BeXBs are much more complicated, so we consider further modifications of the accretion torque by taking into the following factors.

(1) Due to the shear between the magnetosphere of the NS and the accretion disk, there is an open-field region on the magnetosphere, where the accreting material can potentially escape as wind (Lovelace et al. 1995; Romanova et al. 2003). The coupling between the magnetic field of the NS and the wind material generates a spindown torque on the NS. In our model, we assume that a fraction η of the accreting material leaves the system as wind, i.e., ${\dot{M}}_{{\rm{w}}}=\eta \dot{M}$. Then, the torque exerted by the wind N_w is given by

$$\begin{eqnarray}&&{N}_{{\rm{w}}}=-{\eta }_{{\rm{w}}}{\dot{M}}_{{\rm{w}}}\sqrt{{{GMR}}_{{\rm{A}},{\rm{w}}}},\end{eqnarray} \tag{ 12 }$$

where ${R}_{{\rm{A}},{\rm{w}}}$ is the Alfvén radius of the wind, and η_w is an efficiency coefficient for the magnetic field–wind interaction. In our work, we take η_w = 1. Then, the total torque becomes

$$\begin{eqnarray}&&N={N}_{{\rm{d}}}+{N}_{{\rm{w}}}={N}_{0}n(\omega )+{N}_{{\rm{w}}}.\end{eqnarray} \tag{ 13 }$$

(2) Most BeXBs are transient sources subject to (quasi-)periodic type I outbursts (Reig 2011). The spin evolution during type I outbursts is described by

$$\begin{eqnarray}&&{N}_{{\rm{I}}}={N}_{0,{\rm{I}}}n({\omega }_{{\rm{I}}})+{N}_{{\rm{w}},{\rm{I}}},\end{eqnarray} \tag{ 14 }$$

where the subscripts I represent the parameters evaluated for type I outbursts. Assuming that the detection of X-ray pulsations marks the occurrence of the outburst phase, we derive that the duty cycle x of type I outbursts roughly ranges from 0.005 to 0.05 from the data sample compiled by Yang et al. (2017) .

(3) Type II outbursts are rare but more violent compared with type I outbursts. The origin and accretion physics of type II outbursts are still open to debate. While Okazaki et al. (2013) argued that the Bondi–Hoyle–Lyttleton accretion is responsible for type II outbursts, disk accretion seems to be observationally preferred (e.g., Haigh et al. 1999; Sugizaki et al. 2017). Here we adopt the suggestion by Sugizaki et al. (2017) and write the torque N_II during type II outbursts to be

$$\begin{eqnarray}&&{N}_{\mathrm{II}}={N}_{0,\mathrm{II}}n({\omega }_{\mathrm{II}})+{N}_{{\rm{w}},\mathrm{II}},\end{eqnarray} \tag{ 15 }$$

where the subscripts II represent that the parameters are evaluated for type II outbursts. We introduce the duty cycle y to describe the occurrence of type II outbursts, given by

$$\begin{eqnarray}&&y=\mathrm{outburst}\,\mathrm{rate}\times \mathrm{outburst}\,\mathrm{duration}.\end{eqnarray} \tag{ 16 }$$

The observations by Sugizaki et al. (2017) suggest that y is about 0.1 with the outburst rate $\sim 1/1000\,{\mathrm{day}}^{-1}$ and the outburst duration ∼100 days.

(4) During the quiescent phase BeXBs are very faint with X-ray luminosities ${L}_{{\rm{X}},{\rm{l}}}$ as low as $\sim {10}^{33}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ (Yang et al. 2017), so the accretion disk may evolve from an optically thick, geometrically thin disk to an ADAF. Consequently, the accretion disk becomes geometrically thick with sub-Keplerian rotation, so we add a parameter A to modify the disk rotation for an ADAF,

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{D}}}=A{{\rm{\Omega }}}_{{\rm{K}}}.\end{eqnarray} \tag{ 17 }$$

The typical value of A is 0.2–0.3 (Narayan & Yi 1995). Then, the torque N_q in the quiescent phase can be rewritten to be

$$\begin{eqnarray}&&{N}_{{\rm{q}}}={{AN}}_{0,{\rm{q}}}n({\omega }_{{\rm{q}}}/A)+{N}_{{\rm{w}},{\rm{q}}},\end{eqnarray} \tag{ 18 }$$

where the subscript q represents the parameters evaluated in the quiescent phases.

With the above factors taken into account, the averaged total torque exerted by an accretion disk is written to be

$$\begin{eqnarray}&&N\approx {{xN}}_{{\rm{I}}}+{{yN}}_{\mathrm{II}}+(1-x-y){N}_{{\rm{q}}}.\end{eqnarray} \tag{ 19 }$$

Combining Equations (6)–(19) and setting N = 0, we obtain an approximate estimate of the equilibrium spin period P_eq for the NS in a transient BeXB (see the Appendix for the derivation),

$$\begin{eqnarray}\begin{array}{rcl}{P}_{\mathrm{eq}} & \simeq & \left(\displaystyle \frac{2\pi {\phi }^{3/2}}{{2}^{3/14}}\right)\left[\displaystyle \frac{{\mu }^{6/7}}{{\omega }_{{\rm{c}}}{\left({GM}\right)}^{5/7}}\right]\\ & & \times \,\left[\displaystyle \frac{x{\dot{M}}_{{\rm{I}}}^{3/7}+(1-x-y){\dot{M}}_{{\rm{q}}}^{3/7}+y{\dot{M}}_{\mathrm{II}}^{3/7}}{x{\dot{M}}_{{\rm{I}}}^{6/7}{C}_{1,1}^{-1}+(1-x-y){\dot{M}}_{{\rm{q}}}^{6/7}{C}_{1,{\rm{A}}}^{-1}+y{\dot{M}}_{\mathrm{II}}^{6/7}{C}_{1,1}^{-1}}\right],\end{array}\end{eqnarray} \tag{ 20 }$$

where

$$\begin{eqnarray}&&{C}_{\mathrm{1,1}}={\left(1-\displaystyle \frac{{\eta }^{6/7}}{{\phi }^{1/2}}\right)}^{-1},\ \mathrm{and}\ {C}_{1,A}={\left(A-\displaystyle \frac{{\eta }^{6/7}}{{\phi }^{1/2}}\right)}^{-1}.\end{eqnarray} \tag{ 21 }$$

In addition, we require that the magnitude of the spindown torque generated by the wind is always weaker than that of the material torque, that is ${C}_{\mathrm{1,1}}\,\mathrm{and}\,{C}_{1,{\rm{A}}}\gt 0$.

2.2. Spin-change Timescale and Spin Equilibrium

In the accretion-powered state, the spin evolution of the NS is determined by the total torque exerted on the NS by the accretion disk

$$\begin{eqnarray}&&-2\pi I\displaystyle \frac{\dot{P}}{{P}^{2}}=N,\end{eqnarray} \tag{ 22 }$$

where I is the moment of inertia of the NS and $\dot{P}$ is the derivative of the spin period. This equation can be solved analytically, and according to its solution, we obtain an estimate of the typical spin-change timescale τ_spin (see the Appendix for the derivation),

$$\begin{eqnarray}&&{\tau }_{\mathrm{spin}}=\left[\displaystyle \frac{I{\omega }_{{\rm{c}}}}{{\phi }^{2}}\displaystyle \frac{{\left(2{GM}\right)}^{2/7}}{{\mu }^{8/7}}\right]{\left(\displaystyle \frac{{GM}}{{\bar{L}}_{{\rm{X}}}R}\right)}^{3/7},\end{eqnarray} \tag{ 23 }$$

where $\dot{M}$ is converted into L_X via

$$\begin{eqnarray}&&{L}_{{\rm{X}}}=\displaystyle \frac{{GM}\dot{M}}{R},\end{eqnarray} \tag{ 24 }$$

R is the radius of the NS, and ${\bar{L}}_{{\rm{X}}}$ is the averaged X-ray luminosity, defined by

$$\begin{eqnarray}&&{\bar{L}}_{{\rm{X}}}^{3/7}={{xL}}_{{\rm{X}},{\rm{I}}}^{3/7}+(1-x-y){L}_{{\rm{X}},{\rm{q}}}^{3/7}+{{yL}}_{{\rm{X}},\mathrm{II}}^{3/7}.\end{eqnarray} \tag{ 25 }$$

It is interesting to note that τ_spin is independent of the structure of the accretion disk and the wind loss, though these two factors affect the value of P_eq. Equation (23) shows that ${\tau }_{\mathrm{spin}}\propto {\bar{L}}_{{\rm{X}}}^{-3/7}{\mu }^{-8/7}$, i.e., systems with higher average accretion rate and stronger magnetic field can reach spin equilibrium with shorter time. Assuming steady accretion and adopting the following parameters: x = 0, y = 0, ${L}_{{\rm{X}},{\rm{q}}}={10}^{34}$ erg s⁻¹, μ = 10³⁰ G cm³, R = 10⁶ cm, M = 1.4 M_☉, I = 10⁴⁵ g cm², ϕ = 0.5, and ω_c = 0.8, we obtain τ_spin ≲ 10⁶ yr, which implies that a typical NS is able to reach the spin equilibrium within 10⁶ yr even if we only take into account accretion during quiescence. Obviously τ_spin become shorter if we consider the contribution from type I/II outbursts. Considering that the ages of most BeXBs are roughly 10 Myr (Antoniou et al. 2010), it is safe to assume that most of the observed BeXPs have evolved to be near the spin equilibrium.

We then perform a numerical test of the above analysis. We adopt the following parameters: P₀ = 1 s, $({L}_{{\rm{X}},{\rm{I}}},{L}_{\mathrm{II}},{L}_{{\rm{X}},{\rm{q}}})\,=({10}^{37},{10}^{38},{10}^{34})\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, μ = 10³⁰ G cm³, R = 10⁶ cm, A = 0.1 or 1, x = 0.05 or 0.005, y = 0 or 0.1, P_orb = 100 days, η = 0 or 0.1, and ω_c = 0.8. Our results are shown in Figure 1. In this figure, the upper and lower panels depict the results with x = 0.05 and 0.005, respectively. In the first, second, third, and fourth columns, we set $(A,y,\eta )$ to be (1, 0, 0), (1, 0, 0.1), (0.1, 0, 0.1), and (0.1, 0.1, 0.1), respectively. In each panel, the blue line represents the evolutionary track of P, the black vertical line the age of 10⁶ yr, and the black horizontal line the equilibrium spin period P_eq given by Equation (20). The numerical calculations confirm our analytical result that current BeXPs should be close to spin equilibrium. A remarkable feature is that type II outbursts can significantly change the overall distribution of the equilibrium spin periods. Cheng et al. (2014) showed that the BeXPs with P < 40 s are more likely to experience type II outbursts compared with those with P > 40 s. The calculated results presented in Figure 1 demonstrates that the NSs can be significantly spun up during type II outbursts.

Zoom In Zoom Out Reset image size

In the next section we perform statistical calculations to examine whether the equilibrium spin periods in transient accretion can reproduce the bimodal distribution of BeXPs.

We first describe the data sample that will be used to compare with the model prediction. There are 103 BeXPs in the sample, including the sources from MW, SMC, and LMC (data taken from Raguzova & Popov 2005; Reig 2011; Townsend et al. 2011; Sturm et al. 2012; Vasilopoulos et al. 2014; Coe & Kirk 2015; Vasilopoulos et al. 2016, 2017; Yang et al. 2017). The distribution of the pulse period P is presented in Figure 2, where the black vertical line represents P = 40 s. The upper panel shows the histogram of the real P distribution, and the lower panel depicts the normalized distribution of P (blue solid line) and the fitted curve (green solid line)⁴ with a double-Gaussian function, which is defined as

$$\begin{eqnarray}\begin{array}{rcl}F(x) & = & \displaystyle \frac{A}{\sqrt{2\pi {\sigma }_{1}^{2}}}\exp \left[-\displaystyle \frac{{\left(x-{\mu }_{1}\right)}^{2}}{2{\sigma }_{1}^{2}}\right]\\ & & +\,\displaystyle \frac{1-A}{\sqrt{2\pi {\sigma }_{2}^{2}}}\exp \left[-\displaystyle \frac{{\left(x-{\mu }_{2}\right)}^{2}}{2{\sigma }_{2}^{2}}\right].\end{array}\end{eqnarray} \tag{ 26 }$$

The fitting parameters are as follows

• A = 0.321 ± 0.057;
• ${\mu }_{1}=0.953\pm 0.056$, μ₂ = 2.298 ± 0.067;
• ${\sigma }_{1}=0.296\pm 0.063$, σ₂ = 0.524 ± 0.071;

where x, μ_1,2, and σ_1,2 are in units of seconds, and the error for each parameter was calculated by the fitting method.⁵ We apply the Hartigans' dip test (Hartigan & Hartigan 1985) to our sample. The P-value is 0.4308 for the distribution in logarithmic scale, and 0.9918 in linear scale.⁶

Zoom In Zoom Out Reset image size

Then, we introduce the properties of the input parameters for our statistical calculation:

• 1.  
  The duty cycle x of type I outbursts: based on the data compiled in Yang et al. (2017) for the occurrence of X-ray outbursts and detection of X-ray pulsations, we adopt a log-uniform distribution in the range of 0.005–0.05.
• 2.  
  The X-ray luminosity ${L}_{{\rm{X}},{\rm{I}}}$ during type I outbursts: we adopt a log-uniform distribution in the range of 10³⁶–10³⁷ erg s⁻¹ (Reig 2011).
• 3.  
  The X-ray luminosity L_X,II during type II outbursts: Cheng et al. (2014) showed that the peak X-ray luminosity of the BeXPs is truncated around $2\,\times {10}^{38}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, and the distribution seems to prefer the high-luminosity end. We accordingly adopt a uniform distribution in the range of 10³⁷–2 × 10³⁸ erg s⁻¹.
• 4.  
  The X-ray luminosity in the quiescent phase ${L}_{{\rm{X}},{\rm{q}}}$: we adopt a log-uniform distribution in the range of 10³⁴–10³⁶ erg s⁻¹ (Yang et al. 2017).
• 5.  
  The critical value of the fastness parameter ω_c: various studies show that it is in the range of 0.35–0.95 (Ghosh & Lamb 1979a, 1979b; Wang 1995; Li & Wang 1996; Long et al. 2005; Zanni & Ferreira 2009, 2013). In our work, we adopt a uniform distribution in the range of 0.5–1.
• 6.  
  The ADAF parameter A: we adopt a uniform distribution in the the range 0.2–0.3 (Narayan & Yi 1995).

In the next section we will vary the distribution forms and examine their influence. Other (free) parameters are listed as follows:

• 1.  
  The disk wind fraction η.
• 2.  
  The NS magnetic fields B: we adopt a log-normal distribution for B, i.e.,

  $$\begin{eqnarray}&&F(\mathrm{lg}B| {\mu }_{\mathrm{lg}B},{\sigma }_{\mathrm{lg}B})=\displaystyle \frac{1}{\sqrt{2\pi {\sigma }_{\mathrm{lg}B}^{2}}}\exp \left[-\displaystyle \frac{{\left(\mathrm{lg}B-{\mu }_{\mathrm{lg}B}\right)}^{2}}{2{\sigma }_{\mathrm{lg}B}^{2}}\right],\end{eqnarray} \tag{ 27 }$$

  where the average value ${\mu }_{\mathrm{lg}B}$ and the variance ${\sigma }_{\mathrm{lg}B}$ are left as free parameters.
• 3.  
  The fraction f of BeXPs that undergo type II outbursts.
• 4.  
  The duty cycle y of type II outbursts.

With the above parameters and the equilibrium spin period assumption we model the bimodal distribution of BeXPs with a Markov Chain Monte Carlo method.⁷ We consider three types of solutions in our model: the no-wind solution with η = 0, the weak-wind solution with η = 0.1, and the strong-wind solution with η = 0.5 (hereafter Solution-1, Solution-2, and Solution-3, respectively). The calculated results are shown in Figure 3.⁸ There are three rows and two columns in each figure. The first, second, and third rows present the results of Solution-1, Solution-2, and Solution-3. The left column presents the posterior probability distribution of ${\mu }_{\mathrm{lg}{\rm{B}}}$, ${\sigma }_{\mathrm{lg}{\rm{B}}}$, f, and y, where the blue solid lines represent the best fitting points, and the right column presents the comparison between the observed distribution (blue solid lines) and the model distribution (green solid lines) with the best fitting parameters, which are summarized in Table 1.

Zoom In Zoom Out Reset image size

Table 1.  Fitting Parameters for the Period Distribution^a

η	${\mu }_{\mathrm{lg}B}$	${\sigma }_{\mathrm{lg}B}$	f	y
0	${13.34}_{-0.11}^{+0.09}$	${0.45}_{-0.07}^{+0.10}$	${0.41}_{-0.06}^{+0.40}$	${0.33}_{-0.17}^{+0.62}$
0.1	${12.97}_{-0.11}^{+0.12}$	${0.39}_{-0.09}^{+0.12}$	${0.39}_{-0.07}^{+0.07}$	${0.08}_{-0.04}^{+0.34}$
0.5	${12.15}_{-0.12}^{+0.12}$	${0.28}_{-0.09}^{+0.12}$	${0.36}_{-0.06}^{+0.05}$	${0.04}_{-0.01}^{+0.04}$

Note.

^aThe errors are given under the 90% confidence level.

Download table as:  ASCIITypeset image

It is seen from Figure 3 and Table 1 that the three solutions can well reproduce the bimodal distribution of the spin periods. The short-period subpopulation mainly consists of BeXBs, which undergo type II outbursts. There is a tendency that μ_lgB, f, and y all become smaller with increasing η. This is not difficult to understand. As the wind torque acts to spin down the NS, a weaker magnetic field is required with stronger wind, given a specific equilibrium spin period. When η increases from 0 to 0.5, the average value of μ_lgB decreases from 13.34 (B ∼ 2.2 × 10¹³ G) to 12.15 (B ∼ 1.4 × 10¹² G), which are still in the reasonable range for X-ray pulsars. The fraction f of BeXBs with type II outbursts varies from 0.41 to 0.36. According to the data in Cheng et al. (2014) we infer that the observed value of f for transient BeXBs is roughly 0.2. But this should be taken as a lower limit because type II outbursts must have been missed because of limited monitoring by X-ray telescopes. The duty cycle y of type II outbursts changes from 0.33 to 0.04, which is roughly consistent with the observed value ∼0.1.

We mention that the derived parameter distributions depend on the values of the input parameters, which are subject to some uncertainties. For example, in our calculation, we use a constant luminosity during outbursts and ignore its evolution during the rising and decay phases for simplicity. The adopted distributions of ${L}_{{\rm{X}},{\rm{I}}}$, ${L}_{{\rm{X}},\mathrm{II}}$, and ${L}_{{\rm{X}},{\rm{q}}}$ are somewhat speculative and need to be refined by more dense observations. In addition, we neglect the possible influence of the eccentricities in BeXBs and the possible dependence of f on both the orbital period and the eccentricity. To see how the specific distribution functions of the input parameters potentially affect our results, we compare the results by taking into account different combinations of the distribution functions, with LU for log-uniform and U for uniform distributions. In the cases of Solution-2 and Solution-3, and we consider the following combinations for $(x,{L}_{{\rm{X}},{\rm{I}}},{L}_{{\rm{X}},\mathrm{II}},{L}_{{\rm{X}},{\rm{q}}})$,: (LU, LU, U, LU), (U, LU, U, LU), (LU, U, U, LU, U), (LU, LU, U, U), and (LU, LU, LU, LU). Our results are summarized in Table 2. We see that in these different combinations we obtain results that are compatible with each other. Therefore, our results are not sensitively dependent on the adopted distribution functions.

Table 2.  Fitting Parameters for Different Combinations of the Distribution Functions

η	Distributions of Parametersa				${\mu }_{\mathrm{lg}B}$	${\sigma }_{\mathrm{lg}B}$	f	y
	x	${L}_{{\rm{X}},{\rm{I}}}$	${L}_{{\rm{X}},\mathrm{II}}$	${L}_{{\rm{X}},{\rm{q}}}$
	LU	LU	U	LU	${12.97}_{-0.11}^{+0.12}$	${0.39}_{-0.09}^{+0.12}$	${0.39}_{-0.07}^{+0.07}$	${0.08}_{-0.04}^{+0.34}$
	U	LU	U	LU	13.12${}_{-0.10}^{+0.11}$	${0.40}_{-0.08}^{+0.10}$	${0.37}_{-0.05}^{+0.08}$	${0.15}_{-0.04}^{+0.77}$
0.1	LU	U	U	LU	13.10${}_{-0.10}^{+0.11}$	${0.40}_{-0.09}^{+0.10}$	${0.37}_{-0.05}^{+0.07}$	${0.15}_{-0.06}^{+0.74}$
	LU	LU	U	U	${12.87}_{-0.11}^{+0.11}$	${0.41}_{-0.08}^{+0.12}$	${0.40}_{-0.06}^{+0.41}$	${0.07}_{-0.02}^{+0.09}$
	LU	LU	LU	LU	${12.96}_{-0.09}^{+0.11}$	${0.38}_{-0.08}^{+0.12}$	${0.37}_{-0.05}^{+0.08}$	${0.18}_{-0.06}^{+0.72}$
	LU	LU	U	LU	${12.15}_{-0.12}^{+0.12}$	${0.28}_{-0.09}^{+0.12}$	${0.36}_{-0.06}^{+0.05}$	${0.04}_{-0.01}^{+0.04}$
	U	LU	U	LU	12.39${}_{-0.14}^{+0.10}$	${0.35}_{-0.11}^{+0.11}$	${0.36}_{-0.06}^{+0.06}$	${0.08}_{-0.02}^{+0.57}$
0.5	LU	U	U	LU	12.35${}_{-0.13}^{+0.11}$	${0.32}_{-0.10}^{+0.13}$	${0.36}_{-0.06}^{+0.05}$	${0.06}_{-0.01}^{+0.43}$
	LU	LU	U	U	${11.90}_{-0.15}^{+0.08}$	${0.28}_{-0.12}^{+0.15}$	${0.37}_{-0.07}^{+0.05}$	${0.03}_{-0.01}^{+0.02}$
	LU	LU	LU	LU	${12.15}_{-0.11}^{+0.12}$	${0.27}_{-0.17}^{+0.14}$	${0.36}_{-0.06}^{+0.05}$	${0.09}_{-0.03}^{+0.26}$

Note.

^aLU and U for log-uniform distribution and uniform distribution respectively.

Download table as:  ASCIITypeset image

We summarize our work as follows:

• (1)  
  NSs in most BeXBs undergo spin-up during outbursts and spin down during quiescence, so their spin evolution is controlled by the competition between the spin-up and spindown torques. Therefore, it is inappropriate to use long-term, average steady accretion to calculate the spin evolution and derive the magnetic fields of the accreting NSs.
• (2)  
  Most NSs in BeXBs have reached the equilibrium spin periods because the spin-change timescale is sufficiently short compared to the ages of BeXBs.
• (3)  
  Considering the transient behavior of BeXBs, it seems able to reproduce the bimodal spin period distribution with the assumption of the equilibrium spin period for NSs with typical magnetic fields ($\sim {10}^{12}\mbox{--}{10}^{13}$ G). The critical factors that determine the spin period distribution are the properties of type I and II outbursts.

This work was supported by the National Key Research and Development Program of China (2016YFA0400803), the Natural Science Foundation of China under grant Nos. 11333004, 11773015, 11573016, and Project U1838201 supported by NSFC and CAS.

Software: emcee (Foreman-Mackey et al. 2013), corner.py (Foreman-Mackey 2016), scipy (Jones et al. 2001), R (R Core Team 2018), dip test (Maechler 2016).

$$\begin{eqnarray}&&\begin{array}{l}I\displaystyle \frac{2\pi }{{P}^{2}}(-\dot{P})=\left[{{xN}}_{0,{\rm{I}}}\left(1-\displaystyle \frac{{\eta }^{6/7}}{{\phi }^{1/2}}\right)\right.\\ \qquad +\,\left.(1-x-y){N}_{0,{\rm{q}}}\left(A-\displaystyle \frac{{\eta }^{6/7}}{{\phi }^{1/2}}\right)+{{yN}}_{0,\mathrm{II}}\left(1-\displaystyle \frac{{\eta }^{6/7}}{{\phi }^{1/2}}\right)\right]\\ \qquad -\,\left[{{xN}}_{0,{\rm{I}}}\displaystyle \frac{{\omega }_{{\rm{I}}}}{{\omega }_{{\rm{c}}}}+(1-x-y){N}_{0,{\rm{q}}}\displaystyle \frac{{\omega }_{{\rm{q}}}}{{\omega }_{{\rm{c}}}}+{{yN}}_{0,\mathrm{II}}\displaystyle \frac{{\omega }_{\mathrm{II}}}{{\omega }_{{\rm{c}}}}\right]\\ \quad =\,\displaystyle \frac{{\phi }^{1/2}{\left({GM}\right)}^{3/7}{\mu }^{2/7}}{{2}^{1/14}}\left[x{\dot{M}}_{{\rm{I}}}^{6/7}\left(1-\displaystyle \frac{{\eta }^{6/7}}{{\phi }^{1/2}}\right)\right.\\ \qquad +\,\left(1-x-y\right){\dot{M}}_{{\rm{q}}}^{6/7}\left(A-\displaystyle \frac{{\eta }^{6/7}}{{\phi }^{1/2}}\right)\\ \qquad +\,\left.y{\dot{M}}_{\mathrm{II}}^{6/7}\left(1-\displaystyle \frac{{\eta }^{6/7}}{{\phi }^{1/2}}\right)\right]-\displaystyle \frac{2\pi {\phi }^{2}}{{\omega }_{{\rm{c}}}}\displaystyle \frac{{\mu }^{8/7}}{{\left(2{GM}\right)}^{2/7}}\\ \qquad \times \,\left[x{\dot{M}}_{{\rm{I}}}^{3/7}+(1-x-y){\dot{M}}_{{\rm{q}}}^{3/7}+y{\dot{M}}_{\mathrm{II}}^{3/7}\right]{P}^{-1}.\end{array}\end{eqnarray} \tag{ 28 }$$

Here we define two parameters C₁ and C₂, which are, respectively, given by

$$\begin{eqnarray}&&\begin{array}{l}{C}_{1}=\displaystyle \frac{{\phi }^{1/2}{\left({GM}\right)}^{3/7}{\mu }^{2/7}}{{2}^{1/14}}\left[x{\dot{M}}_{{\rm{I}}}^{6/7}\left(1-\displaystyle \frac{{\eta }^{6/7}}{{\phi }^{1/2}}\right)+\left(1-x-y\right)\right.\\ \quad \times \,\left.{\dot{M}}_{{\rm{q}}}^{6/7}\left(A-\displaystyle \frac{{\eta }^{6/7}}{{\phi }^{1/2}}\right)+y{\dot{M}}_{\mathrm{II}}^{6/7}\left(1-\displaystyle \frac{{\eta }^{6/7}}{{\phi }^{1/2}}\right)\right]\displaystyle \frac{1}{2\pi I}\end{array}\end{eqnarray} \tag{ 29 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}{C}_{2} & = & \displaystyle \frac{2\pi {\phi }^{2}}{{\omega }_{{\rm{c}}}}\displaystyle \frac{{\mu }^{8/7}}{{\left(2{GM}\right)}^{2/7}}\\ & & \times \,\left[x{\dot{M}}_{{\rm{I}}}^{3/7}+(1-x-y){\dot{M}}_{{\rm{q}}}^{3/7}+y{\dot{M}}_{\mathrm{II}}^{3/7}\right]\displaystyle \frac{1}{2\pi I}.\end{array}\end{eqnarray} \tag{ 30 }$$

Then, Equation (28) can be rewritten as

$$\begin{eqnarray}&&\dot{P}=-({C}_{1}P-{C}_{2})P.\end{eqnarray} \tag{ 31 }$$

The equilibrium spin P_eq is given by setting $\dot{P}=0$, i.e.,

$$\begin{eqnarray}&&{P}_{\mathrm{eq}}=\displaystyle \frac{{C}_{2}}{{C}_{1}}.\end{eqnarray} \tag{ 32 }$$

Assuming that the torque generated by the wind is smaller than the material torque, mathematically $(1-{\eta }^{6/7}/{\phi }^{1/2})\gt 0$ and $(A-{\eta }^{6/7}/{\phi }^{1/2})\gt 0$, we have ${C}_{1}\ne 0$. Then, Equation (31) has the following solution

$$\begin{eqnarray}&&P=\displaystyle \frac{{C}_{2}/{C}_{1}}{1-(1-{C}_{2}/{C}_{1}/{P}_{0}){e}^{-{C}_{2}T}},\end{eqnarray} \tag{ 33 }$$

where T is the evolutionary time, and P₀ is the spin period at T = 0. Note that, according to Equation (32), C₂/C₁ can be substituted by P_eq. Then, the spin evolution is given by

$$\begin{eqnarray}&&P=\displaystyle \frac{{P}_{\mathrm{eq}}}{1-(1-{P}_{\mathrm{eq}}/{P}_{0}){e}^{-{C}_{2}T}}.\end{eqnarray} \tag{ 34 }$$

According to the exponential term in Equation (34), we can define the spin-evolutionary timescale τ_spin,

$$\begin{eqnarray}&&{\tau }_{\mathrm{spin}}={C}_{2}^{-1}.\end{eqnarray} \tag{ 35 }$$