Wind-accreting Symbiotic X-ray Binaries

TL;DR

This paper introduces an improved model for symbiotic X-ray binaries that accounts for non-stationary accretion processes, better matching observed distributions of system properties and suggesting low matter expulsion efficiency during binary evolution.

Contribution

The model incorporates non-stationary quasi-spherical accretion and stricter criteria for plasma entry and settling, enhancing agreement with observed system distributions.

Findings

Model reproduces observed distributions of orbital and spin periods.

Galactic number of systems aligns with observations.

Low matter expulsion efficiency is necessary for system formation.

Abstract

We present a new model of the population of symbiotic X-ray binaries (SyXBs) that takes into account non-stationary character of quasi-spherical sub-sonic accretion of the red giant's stellar wind onto slowly rotating neutron stars. Updates of the earlier models are given, which include more strict criteria of slow NS rotation for plasma entry into the NS magnetosphere via Rayleigh-Taylor instability, as well as more strict conditions for settling accretion for slow stellar winds, with an account of variations in the specific angular momentum of captured stellar wind in eccentric binaries. These modifications enabled a more adequate description of the distributions of observed systems over binary orbital periods, NS spin periods and their X-ray luminosity in the $\sim 10^{32}-10^{36}$~erg s$^{-1}$ range and brought their model Galactic number into reasonable agreement with the observed one. Reconciliation of the model and observed orbital periods of SyXBs requires a low efficiency of matter expulsion from common envelopes during the evolution that results in the formation of NS-components of symbiotic X-ray systems.

Figures

Tables

Table 1. Table 1: Parameters of observed and suspected SyXBs. Listed are the source name, the NS spin period P ∗ superscript 𝑃 P^{*} (s), the binary orbital period P orb subscript 𝑃 orb P_{\rm orb} , the X-ray luminosity L X subscript 𝐿 X L_{\rm X} , the source distance d 𝑑 d .

SyXB	$P^{\ast }$	$P_{\mathrm{orb}}$	$L_X$	d
	(s)	(day)	(erg s-1)	(Kpc)
GX 1+4	$\simeq 140$(1)	1161(2,3)	${10}^{35}$—${10}^{36}$(4)	4.3(2)
		295$\pm$70(5)
		304(6,7)
4U 1954+319	$\sim 18300$(8)	$\gtrsim 400$?(9)	$4\times {10}^{32}$(9)	1.7(9)
4U 1700+24	?	404$\pm$20(6)	$2\times {10}^{32}-{10}^{34}$(10)	$0.42\pm 0.4$(10)
		4391(23)
Sct X-1	113(11)	?	$2\times {10}^{34}$(11)	$\ge 4^{\left(11\right)}$
IGR J16194-2810	?	?	$\le 7\times {10}^{34}$(12)	$\le {3.7}^{\left(12\right)}$
IGR J16358-4726	5850(13)	?	$3\times {10}^{32}-2\times {10}^{36}$(14)	5-6; 12-13(15)
CGCS 5926	?	$\sim {151}^{\left(16\right)}$	$\le 3\times {10}^{32\left(16\right)}$	${5.2}^{\left(16\right)}$
CXOGBS J173620.2-293338	?	?	$\sim 9\times {10}^{32}$(17)	?
XTE J1743-363	?	?	?	$\sim 5^{\left(18\right)}$
XMMU J174445.5-295044	?	?	$\gtrsim 4\times {10}^{34}$(19)	${3.1}_{-1.1}^{+1.8}$(19)
3XMM J181923.7-170616IGR	$407.9$(20)	?	$2.78\times {10}^{34}{d}_{10}^2$(20)	?
IGR J17329-2731	$6680\pm 3$(21)	?	?	${2.7}_{-1.2}^{+3.4}$(21)
IGR J17197-3010			$\lesssim 1.6\times {10}^{35}$(22)	6.3-16.6(22)

Equations

$$P^{\ast }>27[\mathrm{s}]{\dot{M}}_{16}^{1/5}{\mu }_{30}^{33/35}(M_{\mathrm{x}}/M_{\odot }{)}^{-97/70}.$$

$$f(u)\approx {\left[\frac{t_{\mathrm{ff}}}{t_{\mathrm{cool}}}\right]}^{1/3},$$

$${\dot{M}}_{\mathrm{B}}\approx \frac{1}{4}{\dot{M}}_{\mathrm{o}}\frac{v_{\mathrm{rel}}}{v_{\mathrm{w}}}{\left(\frac{R_{\mathrm{B}}}{a}\right)}^2,$$

$$R_{\mathrm{B}}=\delta \frac{2G{M}_{\mathrm{x}}}{v_{\mathrm{rel}}^2+c_{\mathrm{s}}^2}.$$

$$R_{\mathrm{A}}\sim {\left(\frac{f(u){\mu }^2}{{\dot{M}}_{\mathrm{x}}\sqrt{G{M}_{\mathrm{x}}}}\right)}^{2/7}.$$

$$f(u{)}_{\mathrm{Comp}}\approx 0.22{\zeta }^{7/11}{\dot{M}}_{x,16}^{4/11}{\mu }_{30}^{-1/11}$$

$$f(u{)}_{\mathrm{rad}}\approx 0.1{\zeta }^{14/81}{\dot{M}}_{x,16}^{6/27}{\mu }_{30}^{2/27}$$

$${\dot{M}}_{x,16}^{\mathrm{Comp}}\simeq 0.1\zeta {\dot{M}}_{B,16}^{11/7}{\mu }_{30}^{-1/7}$$

$${\dot{M}}_{x,16}^{\mathrm{rad}}\simeq 0.05{\zeta }^{2/9}{\dot{M}}_{B,16}^{9/7}{\mu }_{30}^{2/21}$$

$${\dot{M}}_{B,16}\gtrsim 0.1{\zeta }^{-2/9}{\mu }_{30}^{5/6}.$$

$$\frac{dI{\omega }^{\ast }}{dt}=\eta Z{\dot{M}}_{\mathrm{x}}{\omega }_{\mathrm{B}}{R}_{\mathrm{B}}^2-Z(1-z/Z){\dot{M}}_{\mathrm{x}}{\omega }^{\ast }{R}_{\mathrm{A}}^2,$$

$$\delta {\dot{M}}_{\mathrm{B}}/{\dot{M}}_{\mathrm{B}}=\delta {\rho }_{\mathrm{w}}/{\rho }_{\mathrm{w}}-3\delta v_{\mathrm{w}}/v_{\mathrm{w}}(v_{\mathrm{w}}/v_{rel}{)}^2.$$

$${\alpha }_{\mathrm{CE}}\left(-\frac{G{M}_1{M}_2}{2a_i}+\frac{G{M}_c{M}_2}{2a_f}\right)=-\frac{G(M_1-M_c)M_1}{\lambda R_1}.$$

$$\frac{\Delta J}{J}=\gamma \frac{\Delta M}{M},$$

$$\frac{\mathrm{SFR}(t)}{M_{\odot}\,\mathrm{yr}^{-1}}=\left\{\begin{array}[]{lr}11e^{-\frac{t-t_{0}}{\tau}}+0.12(t-t_{0}),&t\geq t_{0}\\ 0,&t<t_{0}\end{array}\right.$$

$$r=\frac{p}{1+e\cos \theta },$$

$$v_{\mathrm{w}}(r)=v_{\mathrm{esc}}f(r).$$

$$f(r)=\alpha {\left(1-\frac{R_{\mathrm{o}}}{r}\right)}^{\beta },$$

$$f(r)=\left(1-\frac{5\text{mbox}km/s}{v_{\mathrm{esc}}}\right)\exp \left(-{\left[\frac{1}{5}\left(\frac{r}{R_{\mathrm{o}}}-1\right)\right]}^2\right)+\frac{5\text{mbox}km/s}{v_{\mathrm{esc}}}.$$

$$v_n=\sqrt{\frac{GM}{p}}(1+e\cos \theta ),$$

$$v(r)=\sqrt{(v_{\mathrm{w}}(r)-v_r{)}^2+v_n^2}.$$

$$v_r=\sqrt{\frac{GM}{p}}(e\sin \theta ).$$

$$K_{sd}(r)=Z{\dot{M}}_{\mathrm{x}}{v}_n{R}_{\mathrm{B}}(R_{\mathrm{B}}/r)$$

$$P_{\mathrm{eq}}^e=2\text{piup}\frac{⟨{\dot{M}}_{\mathrm{B}}{R}_{\mathrm{A}}^2⟩}{⟨{\dot{M}}_{\mathrm{B}}{v}_n{R}_{\mathrm{B}}^2/r⟩}.$$

$${\dot{M}}_{\mathrm{B}}(r)=\frac{1}{4}{\dot{M}}_{\mathrm{o}}\left(\frac{v(r)}{v_{\mathrm{w}}(r)}\right){\left(\frac{R_{\mathrm{B}}(r)}{r}\right)}^2,$$

$$F(e)=\frac{\sqrt{\frac{GM}{a}}}{v_0}\frac{⟨{\left(\frac{1+e\cos \theta }{1-e^2}\right)}^{6/7}{\left(\frac{v}{v_{\mathrm{w}}}\right)}^{3/7}{\left(\frac{v_0}{v}\right)}^{12/7}⟩{\left(\frac{v_0}{v_{w,0}}\right)}^{4/7}}{⟨{\left(\frac{1+e\cos \theta }{1-e^2}\right)}^3\left(\frac{v_n}{v_{\mathrm{w}}}\right){\left(\frac{v_0}{v}\right)}^7⟩}\le 1.$$

$$v_0=\sqrt{\frac{GM}{a}}{\left[\frac{2M_{\mathrm{o}}}{M}\frac{a}{R_{\mathrm{o}}}f(a{)}^2+1\right]}^{1/2}$$

Full text

††pubyear: 2018††pagerange: Wind-accreting Symbiotic X-ray Binaries–A

1 Introduction

Symbiotic X-ray binaries (SyXBs) represent a small group of binary systems hosting an accreting neutron star (NS) and a late-type giant (K1-M8). Candidate systems CGCS 5926 (Masetti et al., 2011) and CXOGBS J173620.2-29333 (Hynes et al., 2014) possibly host carbon stars. Despite the fact that the first SyXB was identified more than 40 yr ago (Davidsen et al., 1977), currently, only about a dozen of SyXBs are detected (see Table 1 listing the basic parameters of objects known at the time of writing).

In SyXBs, a NS can accrete matter lost by the low-mass red giant companion as a relatively slow stellar wind ($v_{\mathrm{w}}\sim$ several 10 km s*-1*) and during the Roche lobe overflow. SyXBs are characterized by the relatively low mean X-ray luminosity $L_{\mathrm{x}}\sim 10^{32}-10^{36}$ erg s*-1* and manifest outburst activity typical for wind-accreting NSs. The relative faintness and non-stationarity of the sources make it difficult to detect and analyse the NS rotation, therefore the NS spin periods $P^{*}$ have been measured only for a handful of sources (see Table 1). Spin periods range from several hundred to several thousand seconds, which is related to the character of accretion onto magnetized NSs (see below). The estimated orbital periods of SyXBs $P_{\mathrm{b}}$ may be as high as several years (e.g., in GX 1+4). Ellipticity of orbits is detected for GX 1+4 ($e$=0.10, Hinkle et al. (2006)) and 4U 1700+24 (e=$0.26\pm 0.15$, Galloway et al. (2002)).

SyXB studies provide an independent probe of stellar winds from the red giant components and the physics of wind accretion onto NSs. The observed characteristics of SyXBs can be reasonably reproduced in the model of settling accretion onto slowly rotating magnetized NS developed by Shakura et al. (2012, 2018), with an account of the NS spin evolution in binary systems (Lipunov et al., 2009). The Galactic population of SyXBs was simulated in our earlier papers (Lü et al., 2012; Kuranov & Postnov, 2015), in which the algorithm of the calculations and the model are presented in more detail.

The present study is motivated by further development of the settling accretion model (Shakura & Postnov, 2017; Shakura et al., 2018) and availability of new observations of SyXBs. The main modifications of the accretion theory can be summarised as follows.

1. The model of stellar winds from low-mass post-main-sequence stars is updated. The stellar wind velocity $v_{\mathrm{w}}$ of low-mass giants, which is crucial for the gravitational capture of the matter by a NS in a binary system is now calculated according to the prescription suggested in Lü et al. (2006).
2. The criteria and parameters of the subsonic settling accretion are improved in accordance with the latest understanding of this regime.
3. An account of the transitions between the Compton and radiative plasma cooling with changing X-ray luminosity and the photoionization heating of the accreting plasma at low relative wind velocities is included in the model.
4. The reduction of NS equilibrium periods in eccentric binaries is taken into account. These amendments enabled us to reach a better description of the observed location of SyXBs in the $P^{*}-L_{\mathrm{x}}$ diagram and to improve the agreement of predicted and observed Galactic numbers of SyXBs.

The paper is organised as follows. The improvements in the settling accretion theory are briefly described in § 2. In § 3 we describe the main modifications of the model. Results of the modeling are presented in § 4. Discussion and conclusion follow in § 5. In Appendix, we present the treatment of the effect of the orbital eccentricity on the NS equilibrium period at the settling accretion stage.

2 Refinement of the settling accretion theory

The settling accretion theory is designed to describe the interaction of accreting plasma with magnetospheres of slowly rotating NSs. The key feature of the model is calculation of the steady plasma entry rate into the magnetosphere due to the Rayleigh-Taylor (the RT) instability regulated by plasma cooling (Compton or radiative). RT instability can be suppressed by fast rotation of the magnetosphere (Arons & Lea, 1980). In the present study we will use the criterion from the latter paper:

[TABLE]

At faster NS periods, the accretion regime at any X-ray luminosity will be supersonic because of efficient plasma penetration into the magnetosphere via Kelvin-Helmholtz instability (Burnard et al., 1983). Here and below, $\dot{M}=10^{16}\dot{M}_{16}\mathrm{g\,s}^{-1}$ is the accretion rate, $\mu=10^{30}\mu_{30}\mathrm{G\,cm}^{3}$ is the NS magnetic moment.

It was found (Elsner & Lamb, 1984; Shakura et al., 2012) that the Compton cooling of accreting plasma by X-ray photons generated near the NS surface enables a free-fall (Bondi-type) supersonic accretion onto the NS magnetosphere provided that the X-ray luminosity is above some critical value $L_{x}^{\dagger}\simeq 4\times 10^{36}$ erg s*-1*. X-ray luminosity $L_{\mathrm{x}}$ of an accreting NS is related to the mass accretion rate $\dot{M}_{\mathrm{x}}$ as $L_{\mathrm{x}}\approx 0.1\dot{M}_{\mathrm{x}}c^{2}$ ($c$ is the speed of light), and this critical $L_{\mathrm{x}}$ corresponds to a mass accretion rate onto NS $\dot{M}_{\mathrm{x}}\simeq 4\times 10^{16}$ g s*-1*.

At lower X-ray luminosities, accreting plasma above the magnetosphere remains hot enough to enable effective inflow into the magnetosphere due to RT instability, and the plasma entry rate is controlled by the plasma cooling rate. This results in the formation of a hot, convective shell above the magnetosphere. In this shell, the plasma gravitationally captured by the NS from the stellar wind of the companion (basically, at the Bondi-Hoyle-Lyttleton rate, $\dot{M}_{\mathrm{B}}$) subsonically settles down towards the NS magnetosphere at a rate $\dot{M}_{\mathrm{x}}\approx f(u)\dot{M}_{\mathrm{B}}$, with the dimensionless factor $f(u)\lesssim 0.5$111In a radial accretion flow with an account of the Compton cooling, the value $f(u)\gtrsim 0.5$ also corresponds to the location of the sonic point in the flow below the magnetospheric boundary $R_{\mathrm{A}}$ (Shakura et al., 2012) enabling a free-fall plasma flow down to $R_{\mathrm{A}}$ and the formation of a shock above the magnetosphere (Arons & Lea, 1976)., whose precise value depends on the plasma cooling regime (Compton or radiative). With good accuracy, this factor can be written as

[TABLE]

where $t_{\mathrm{ff}}$ is the free-fall time at the magnetospheric boundary (the Alfvén radius, $R_{\mathrm{A}}$), $t_{\mathrm{cool}}$ is the characteristic plasma cooling time. Clearly, in the case $t_{\mathrm{cool}}\gg t_{\mathrm{ff}}$, this factor can be very small, leading to an effective decrease in the mass accretion rate onto the NS surface compared to the maximum available Bondi rate, $\dot{M}_{\mathrm{B}}$.

In a circular binary system with component separation $a$, the Bondi-Hoyle-Lyttleton accretion rate can be estimated as

[TABLE]

where $\dot{M}_{\mathrm{o}}$ is the mass-loss rate from the optical star, the gravitational Bondi radius is

[TABLE]

Here $M_{\mathrm{x}}$ is the NS mass, $v_{\mathrm{rel}}^{2}=v_{\mathrm{orb,x}}^{2}+v_{\mathrm{w}}^{2}$ is the relative stellar wind velocity, $v_{\mathrm{orb,x}}$ is the NS orbital velocity, $c_{\mathrm{s}}$ is the sound velocity in the matter. The latter can be ignored in cold stellar winds of late-type stars. The estimate  (3) assumes a spherically-symmetric wind outflow from the optical star and is applicable for $R_{\mathrm{B}}\ll a$. The numerical factor $\delta\sim 1$ in Eq. (4) takes into account the actual location of the bow shock in the stellar wind (Hunt, 1971). Some modern numerical calculations (see, e.g. Liu et al., 2017) suggest a smaller (up to an order of magnitude) mass accretion rates onto a gravitating mass from the stellar wind than given by the standard Bondi-Hoyle-Lyttleton formula, which can be reformulated in terms of smaller values of the parameter $\delta\simeq 0.3-0.5$. Other studies (e.g., de Val-Borro et al., 2017) claim that accretion rate can exceed the Bondi one due to gravitational focusing of the wind. Therefore, it should be born in mind that the numerical factor $\delta$ in Eq. (4) can differ from unity by a factor of a few.

The Alfvén radius is defined from the pressure balance between the accreting plasma and the magnetic field at the magnetospheric boundary and depends on the actual mass accretion rate, $\dot{M}_{\mathrm{x}}$ (i.e., the observable X-ray luminosity $L_{\mathrm{x}}$), and the NS magnetic momentum, $\mu$:

[TABLE]

The factor $f(u)$ equals

[TABLE]

or

[TABLE]

for the Compton and radiative cooling, respectively (Shakura & Postnov, 2017; Shakura et al., 2018). Here $\zeta\lesssim 1$ is the numerical factor determining the characteristic scale of the growing RT mode (in units of the Alfvén radius $R_{\mathrm{A}}$).

Thus, specifying the plasma cooling mechanism, we are able to estimate the expected reduction in the mass accretion rate at the settling accretion stage, $f(u)=F(\dot{M}_{\mathrm{x}},...)=F(f(u)\dot{M}_{\mathrm{B}},...)$, and, by solving for $f(u)$, to find the explicit expression for $\dot{M}_{\mathrm{x}}$ as a function of $\dot{M}_{\mathrm{B}}$ and other parameters.

Putting all things together, we are able to express the expected accretion rate onto the NS (which can be directly probed by the observed X-ray luminosity $L_{\mathrm{x}}$) using the Bondi capture rate $\dot{M}_{\mathrm{B}}$, which can be calculated from the known mass-loss rate of the optical companion $\dot{M}_{\mathrm{o}}$, stellar wind velocity $v_{\mathrm{w}}$ and binary system parameters:

[TABLE]

for the Compton cooling and

[TABLE]

for the radiative cooling. The actual accretion rate onto the NS is taken to be $\dot{M}_{\mathrm{x}}=\max\{\dot{M}_{x}^{\mathrm{Comp}},\dot{M}_{x}^{\mathrm{rad}}$} and determines the actual plasma cooling regime. Matching of Eq. (8) and Eq. (9) shows that with a given NS magnetic momentum, the Compton cooling of plasma dominates for the gravitational capture mass rate

[TABLE]

As the parameter $\zeta$ is always less than 1, this explains why in our previous studies (Kuranov & Postnov, 2015) the Compton cooling was ignored for low mass-accretion rates. However, it becomes important during outbursts caused by the stellar wind parameters variations (see below).

During the settling accretion, a hot convective shell formed above the NS magnetosphere mediates the angular momentum transfer to/from the rotating NS, enabling long-term spin-down episodes with spin-down torques correlated with the X-ray luminosity, as observed, for example, in GX 1+4 (Chakrabarty et al., 1997; González-Galán et al., 2012). Turbulent stresses acting in the shell lead to an almost iso-momentum angular velocity radial distribution, $\omega(r)\sim 1/r^{2}$, suggesting the conservation of the specific angular momentum of gas captured near the Bondi radius $R_{\mathrm{B}}$, $j_{\mathrm{w}}=\eta\omega_{\mathrm{B}}R_{\mathrm{B}}^{2}$, with $\eta\approx 1/4$ (Illarionov & Sunyaev, 1975), where $\omega_{\mathrm{B}}=2\piup/P_{\mathrm{B}}$ is the orbital angular frequency. The numerical coefficient $\eta$ here can vary in a wide range due to inhomogeneities in the stellar wind and can be even negative (Ho et al., 1989); in our simulations, we varied this parameter in the range $\eta=[0,0.25]$.

The condition for quasi-spherical accretion to occur is that the specific angular momentum of a gas parcel is smaller than the Keplerian angular momentum at the Alfvén radius: $j_{\mathrm{w}}\leq j_{K}(R_{\mathrm{A}})=\sqrt{GMR_{\mathrm{A}}}$. In the opposite case, $j_{\mathrm{w}}>j_{K}(R_{\mathrm{A}})$, the formation of an accretion disc around the magnetosphere is possible.

The NS spin evolution at the settling accretion stage can be written in the form of angular momentum conservation equation (Shakura et al., 2012, 2018)

[TABLE]

where $I$ is the NS momentum of inertia, $Z=1/f(u)(u_{c}/u_{\mathrm{ff}})$ is the coupling coefficient of the plasma-magnetosphere interaction, $u_{c}/u_{\mathrm{ff}}\lesssim 1$ is the ratio of the convective velocity of gas in the shell and the free-fall velocity, $z<1$ is a numerical coefficient which takes into account the specific angular momentum of the matter. The variable specific angular momentum of the gravitationally captured stellar wind matter leads to the correction of the equilibrium NS period derived in Shakura et al. (2012, 2018) by the factor $0.25/\eta$. In the possible case of a negative value of $\eta$, the NS could rapidly spin down and even start temporarily rotate in the retrograde direction. However, it is difficult to imagine that such a situation could hold much longer than the orbital binary period. Therefore, possible episodes with negative $\eta$ would somewhat increase the NS equilibrium period $P_{\mathrm{eq}}$ which we recalculate at each time step of our population synthesis calculation. Below we will consider only positive values of $\eta$ and restrict the NS spin rotation by the orbital binary period, $P^{*}\leq P_{\mathrm{orb}}$.

It should be stressed that the hot magnetospheric shell can become dynamically unstable. For example, if large loops of magnetic field are present in the stellar wind, the magnetosphere can become unstable due to magnetic reconnection. This can give rise to short strong outbursts occurring in the dynamical (free-fall) time scale, such as, for example, giant flares observed in supergiant fast X-ray transients (SFXTs) (Shakura et al., 2014; Hubrig et al., 2018). Stellar wind inhomogeneities (especially, variations in the wind velocity $v_{\mathrm{w}}$) also can disturb the settling accretion regime and even lead to free-fall Bondi accretion episodes. Indeed, for long-period binaries with low orbital velocities of the NS, relative variations in the mass accretion rate are related to the wind velocity and density variations by the continuity equation:

[TABLE]

Wind velocity variations are most pronounced for low velocity winds with $\delta v_{\mathrm{w}}/v_{\mathrm{w}}\sim 1$, as in the case of SyXBs. Variations in the mass capture rate can give rise to significant relative variations in the X-ray luminosity, especially in the Compton cooling regime: $\delta L_{\mathrm{x}}/L_{\mathrm{x}}=(11/7)\delta\dot{M}_{\mathrm{B}}/\dot{M}_{\mathrm{B}}$ (see Eq. (8)). Therefore, flaring episodes of Bondi accretion (“high” states) on top of low-luminosity stable settling accretion intervals can be expected. This possibility should be borne in mind when comparing observational data with model calculations (see Section 4 below).

To conclude this Section, we summarise the main changes/additions to the description of the settling accretion regime compared to our previous studies (Lü et al., 2012; Kuranov & Postnov, 2015).

Criterion of slow NS rotation for plasma entry into the NS magnetosphere via RT instability (Arons & Lea, 1980) (Eq. (1)) is applied. For faster NS rotation, accretion onto the NS occurs in the Bondi (supersonic) regime. 2. 2.

Conditions for settling accretion to occur are supplemented by criteria for slow stellar winds. 3. 3.

The equilibrium NS period at the settling accretion stage takes into account variations in the specific angular momentum of the captured stellar wind (the parameter $\eta$). 4. 4.

Reduction of the NS equilibrium spin period due to orbital eccentricity (see Appendix A for more detail) is included.

3 The model

We use a modified version of the openly available BSE code (Hurley et al., 2000, 2002) appended by the block for calculation of spin evolution of magnetized NSs (Lipunov et al., 2009) with an account for settling accretion regime (see also Lü et al., 2012, for more detail).

3.1 Formation of SyXBs

The progenitor binary system of a SyXB with the NS accreting from stellar wind of an evolved low-mass optical companion should harbour a primary star more massive than $\sim 8M_{\odot}$. The primary evolves off the main-sequence, expands and can fill its Roche lobe. For a sufficiently high mass ratio, a common envelope (CE) forms during the first mass-transfer episode (see below). After the CE stage, the helium core of the primary evolves to a supernova that leaves behind a NS remnant (Wilson, 1971).222 Here and below we refer to the first numerical studies in which the formation of a bound remnant – a NS – was demonstrated. If the binary survives the supernova explosion, a young NS appears in a binary system accompanied by a low-mass main-sequence star. The latter, in turn, evolves off the main-sequence, and accretion onto the NS from the stellar wind of the red giant companion becomes possible, leading to the formation of a SyXB.

In a narrow range of primary masses close to the lower mass limit for NS production, after the first mass exchange the evolution of the helium core presumably ends up with the NS formation via an electron-capture supernova (ECSN) (Miyaji et al. (1980); see Siess & Lebreuilly (2018) for the recent study), which is not accompanied by a large NS kick. We assumed that NSs formed via ECSN get a kick of 30 km s*-1*. We found, however, that neither the mass range for ECSN nor the low kick amplitude noticeably affects the properties of the model SyXBs. For NSs formed from Fe-core collapse, the kick velocities of nascent NS were assumed to have a Maxwellian distribution with the 1D rms $\sigma=265$ km s*-1* (Hobbs et al., 2005). In any case, after the SN explosion, the young NS finds itself in an eccentric orbit that is gradually circularised due to the tidal interaction. The orbital circularisation is treated as implemented in the BSE code.

In addition, we consider the possibility of formation of NSs as a result of accretion-induced collapses (AICs) of oxygen-neon (ONe) white dwarfs which accumulated mass close to the Chandrasekhar limit (1.4$M_{\odot}$) via accretion of stellar wind (Canal & Schatzman, 1976). Since these events are expected (i) to be at least by an order of magnitude less energetic than typical iron core-collapse supernovae, (ii) to have a very low amount of mass lost, and (iii) to have a small collapse anisotropy (Dessart et al., 2006), we assign to the NSs produced via this channel kick velocities $30\,\>\mbox{km/s}$, like for descendants of ECSN.

3.2 Common envelope stage

Common envelope is the most prominent example of non-conservative binary evolution enabling an efficient angular momentum removal from the binary on a short time-scale (see e.g. Postnov & Yungelson (2014); Yungelson & Kuranov (2017) for more detail and references). Many details of the binary evolution during the CE stage remain unclear, but it is widely accepted to parameterize the CE evolution by the so-called $\alpha$-formalism (Webbink, 1984), in which the dimensionless parameter $\alpha_{\mathrm{CE}}$ means the fraction of the orbital energy difference between the onset and the end of the CE stage that is spent to unbind the envelope of the evolved (post main-sequence) companion $M_{1}$ from its core $M_{c}$:

[TABLE]

Here the second star $M_{2}$ is assumed to conserve mass during the CE stage, and the binding energy of the stellar envelope is parameterized by the dimensionless parameter $\lambda$ (see calculations for stars of different metallicities in Loveridge et al. (2011)). This model suggests that the product $\alpha_{\mathrm{CE}}\lambda$ determines the ratio of the final to initial binary separations $a_{f}/a_{i}$ after the CE stage: the smaller $\alpha_{\mathrm{CE}}\lambda$ (i.e., “more efficient CE”), the smaller $a_{f}/a_{i}$.

Alternatively, it was also suggested that the first CE-stage of evolution can be controlled by the system angular momentum $J$ conservation (the so-called $\gamma$-formalism that implicitly assumes energy conservation (Nelemans et al., 2000)):

[TABLE]

where $M$ is total mass of the system. Then, the CE evolution can be not accompanied by strong reduction of separation between the binary components. Nelemans & Tout (2005) have shown that the $\gamma$-formalism is able to explain the formation of a variety of binary systems, including low-mass X-ray binaries, if $\gamma$ is confined to a rather narrow range 1.5 – 1.75. If one assumes $\gamma$=1.75, the post-CE separation between the binary components is close to that obtained for $\alpha_{\mathrm{CE}}$=4 (for fixed $\lambda=0.5$), i.e. formally requires additional energy sources, apart from the orbital energy, for expulsion of CE. Clearly, the parameter $\alpha_{\mathrm{CE}}$ can (and probably should) vary from system to system and with evolutionary stage of the donor, but in the absence of a firm theory and any realistic numerical simulations, in our calculations we will treat $\alpha_{\mathrm{CE}}$ as a free parameter, one and the same for all binaries, that can take values both less and higher than unity (see also the discussion in Postnov & Yungelson (2014); Yungelson & Kuranov (2017)).

3.3 Stellar winds from

evolved low-mass stars

The wind mass-loss from red-giants is calculated using the formula from Kudritzki & Reimers (1978) with the parameter $\eta_{\mathrm{w}}=0.5$. Stellar wind velocities from red giants and AGB-stars are described as in Lü et al. (2006, (Eqs. (13)-(17), the terminal wind velocity 5 km s*-1*)). Stellar winds from stripped helium stars are treated according to Vink (2017).

SyXBs usually have moderately non-circular orbits. For low wind velocities, the most important accretion effects are due to the NS orbital velocity variations. In the population synthesis simulations, the orbital eccentricity effects in stellar winds can be taken into account, for example, by averaging the relevant quantities over the orbital period of the binary. For noticeable orbital eccentricities this results in a substantial decrease of the NS equilibrium periods at the settling accretion stage (Postnov et al., 2018). For low eccentricities ($e\lesssim 0.2$) the eccentricity effects are found to be not very significant (see Appendix A for more detail).

3.4 Conditions for settling accretion

The typical NS spin evolution proceeds as follows (see Lipunov et al. (2009) for more details and Lü et al. (2012) for the description of NS spin evolution in SyXBs). Initial NS spin periods are assumed to follow a Gaussian distribution centered at $P_{0}=100$ ms with comparable dispersion (Postnov et al., 2016). Magnetic fields of NS are assumed to be distributed log-normally, $f(\mu)d\mu\sim\exp[-(\log\mu-\log\mu_{0})^{2}/\sigma_{\mu}^{2}]$, with the mean value $\log\mu_{0}=30.35$ and dispersion $\sigma_{\mu}=0.55$. Accretion-induced magnetic field decay is taken into account as in Lü et al. (2012). Initially, the NS is in the “ejector” stage when the pressure of the relativistic wind generated near the NS surface dominates that of the surrounding plasma. Then the NS spins down to the “propeller” stage in which the centrifugal forces at the magnetospheric boundary prevent matter accretion. The propeller stage terminates when the corotation radius $R_{c}=(GM_{\mathrm{x}}P^{*2}/4\piup^{2})^{1/3}$ becomes larger than the magnetosphere radius $R_{\mathrm{A}}$, $R_{c}\geq R_{\mathrm{A}}$, enabling accretion.

In the binary systems without Roche-lobe overflow, we check which type of accretion onto the NS is possible – disc or quasi-spherical one. If the condition for the spherical accretion ($j_{\mathrm{w}}<j_{K}(R_{\mathrm{A}})=\sqrt{GM_{\mathrm{x}}R_{\mathrm{A}}}$) is met, we check the conditions for the settling accretion to occur:

(i) Sufficiently slow rotation of the NS enabling RT instability to provide the plasma entry into the NS magnetosphere (see Eq. (1) above). For faster spinning NS, plasma entry rate into the magnetosphere is regulated by the Kelvin-Helmholtz instability, and free-fall Bondi accretion is realised (Burnard et al., 1983).

(ii) Average X-ray luminosity is below $L^{\dagger}\simeq 4\times 10^{36}$ erg s*-1* to prevent rapid plasma cooling.

(iii) If the relative wind velocity $v_{\mathrm{rel}}$ is below $\simeq 80$ km s*-1*, the photoionization heating of gas above the Bondi radius $R_{\mathrm{B}}$ up to $T_{max}\simeq 5\times 10^{5}$ K is significant (Shakura et al., 2012), no strong bow shock is formed near $R_{\mathrm{B}}$, and the accretion proceeds from the effective Bondi radius $R_{\mathrm{B}}^{*}\simeq 3.5\times 10^{12}$ cm. Therefore, for low wind velocities and wide binary systems, when the Bondi radius calculated by formula  (4) can be formally very large, we set $R_{\mathrm{B}}=\min(R_{\mathrm{B}},50R_{\odot})$.

In binaries with Roche-lobe overflow, only disc accretion takes place.

3.5 Initial distributions

The masses of primary stars are distributed according to the Salpeter law, $dN/dM_{1}\propto M_{1}^{-2.35}$ ($0.1M_{\odot}\leq M_{1}\leq 100M_{\odot}$). Orbital separations of binaries are distributed following Sana et al. (2012). Binary mass ratios $q=M_{2}/M_{1}\leq 1$ and orbital eccentricities $e$ are assumed to have flat distributions: $dN/dq$=const., $dN/de$=const. in the range [0,1]. We assume that Galactic stellar binarity rate is 50%, i.e. 2/3 of stars enter binary systems.

3.6 Galactic star-formation history

To compare the population synthesis results with observations, we convolved the formation rate of SyXB with the Galactic star formation rate (SFR). We adopted a star-formation history in the Galactic bulge and thin disc in the form suggested by Yu & Jeffery (2010):

[TABLE]

with time $t$ in Gyr, $t_{0}$=4 Gyr, $\tau=9$ Gyr, Galactic age 14 Gyr. This model gives a total mass of the Galactic bulge and thin disk $M_{G}=7.2\times 10^{10}M_{\odot}$, which we use for the normalisation of our calculations.

4 Results

In our simulations, we have run $5\times 10^{6}$ binary systems for different values of the parameter $\alpha_{\mathrm{CE}}$ that mostly influences the number of SyXBs and their observed characteristics — NS spin period at the wind-accretion stage $P^{*}$, X-ray luminosity $L_{\mathrm{x}}$ and binary orbital period $P_{\mathrm{B}}$. The results for $\alpha_{CE}=4$, which provides the best agreement with observations, are presented in Figs. 1 – 3.

4.1 Galactic number of SyXBs

In our calculations, we selected all detached systems with the NS components that have companions that reached the first red giant branch (FGB) or are in the later stages of the evolution (core He-burning stage – CHeB, early-AGB stars – EAGB), but do not overflow their Roche lobes. By virtue of such a selection, our sample contains a number of binaries with orbital periods well below several hundred day (see Fig. 3). Strictly speaking, such objects do not comply with the “standard” definition of symbiotic stars that defines them as systems with donors – late giants (e.g., Kenyon & Webbink, 1984), because in binaries with such short orbital periods there is no room for late giants. Examination of the distribution over orbital periods of about 100 symbiotic stars with white dwarf components with known $P_{\mathrm{orb}}$ (Mikołajewska, 2012, Fig. 1) tells that it is justified to consider as SyXB systems only those with $P_{\mathrm{orb}}$  $\buildrel>\over{\sim}$  200 day. This suggests us to consider NS+FGB systems as a separate subclass of wind-accreting low-mass X-ray binaries333Note that these objects have, predominantly, hardly detectable low $L_{\mathrm{x}}$. and leaves as SyXBs predominantly systems with CHeB components. For this reason, we omit NS+FGB stars in Fig. 1 that shows the evolution of the number of Galactic SyXBs with time, calculated by convolution of the SyXB formation delay-time distribution ${\cal R}(t)$ with the assumed SFR rate (Eq. (14)): $N(t)=\int_{0}^{t}\mathrm{SFR}(t-\tau){\cal R}(\tau)d\tau$. Note that even not all CHeB-donor systems fit definition of SyXBs. Symbiotic X-ray binaries with CHeB components form, typically, in less than 1.5 Gyr after the formation. Therefore, their Galactic number evolution simply follows SFR set by Eq. (14).

Figure 1 shows that most of the NSs in SyXBs originate from ECSNs. The reason for this is clear: in their precursors, companions of the NS progenitors are (1-2) $M_{\odot}$ stars. In such a case, a natal kick of only $\sim 100\>\mbox{km/s}$ is sufficient to disrupt the binary in the SN event.

In Fig. 1, we also present the number of SyXBs with NS-components descending from white dwarfs via accretion-induced collapses, which is obtained assuming a 100% efficiency of accretion onto WD. Evidently, this number is an upper limit and can very strongly overestimate the real number of such systems. One of the main reasons is low accretion efficiency due to the well known mismatch between the ranges of accretion rates allowing stable hydrogen and helium burning, which prevents effective growth of mass of the accreting WD (see, e.g., Iben & Tutukov, 1996). As yet, it is impossible to quantify the effect of retention of matter accreted by ONe WDs (which, in principle, may be even negative) in the population synthesis code because of the absence of systematic computations for grids of models in the ($M_{\rm WD},\dot{M}_{\rm accr}$) space. Extrapolation of the extant data for the models of CO WDs may result in systematic errors, since for ONe WDs the masses of the critical shells of hydrogen needed for explosive burning are about twice as large as for CO WDs (Jose & Hernanz, 1998; Gil-Pons et al., 2003), hence the outbursts may be less frequent but more violent and to eject more mass. As follows from the estimates of Iben & Tutukov (1996) and calculations of Piro & Thompson (2014), efficient accumulation of matter by massive WDs in symbiotic systems is possible for accretion rates $\sim 10^{-7}-10^{-8}$  ${\rm M_{\odot}~{}yr^{-1}}$. However, such a relatively high $\dot{M}_{\rm accr}$ is not feasible for about a half of symbiotic stars, if accretion occurs from stellar wind (Lü et al., 2006). Even if the mass of an ONe WD attains $\mathrm{M_{Ch}}$, it is not clear as yet, whether the collapse will result in the formation of a NS or in the white dwarf destruction. According to Schwab et al. (2015); Jones et al. (2016), the collapse occurs if oxygen deflagration starts at the central density not lower than about $10^{10}\mathrm{g\,cm^{-3}}$. The collapse may be prevented by ignition of residual 12C which may be present in WD in non-negligible amount (Schwab & Akira Rocha, 2018). Finally, we note that most of the model systems with NSs resulting from AICs have orbital periods too short to be identified with symbiotic stars. This is in qualitative agreement with Hurley et al. (2010) who have shown, using treatment of retention implemented in BSE, that wide systems harbouring post-AIC NSs are exceptionally rare. Below we will omit post-AIC NSs from consideration because the modeling of their magnetic fields and initial spins requires the study of evolution of their WD progenitors, which is beyond the scope of the present paper, and their vanishingly small number.

To summarise, we can expect that the current number of SyXBs in the Galaxy is slightly less than 40 objects per $10^{10}$ $M_{\odot}$.

4.2 NS spin period-luminosity distribution

Figure 2, left panel, shows the $P^{*}-L_{\mathrm{x}}$ distribution of model wind-accreting NSs with evolved companions with masses below 3 $M_{\odot}$. Here $L_{\mathrm{x}}$ is the X-ray luminosity at the settling accretion stage. In the grey scale shown is the model number of systems normalized to $10^{10}$ $M_{\odot}$. Observed SyXBs from Table 1 with known spin periods and X-ray luminosities are shown by open circles with error bars. Positions of the systems with known luminosities or periods only are indicated by thin dashed lines. Wind-fed low-mass X-ray binaries with optical components of different types are shown separately, from left to right, top to bottom: NS with optical star at the first giant branch (FGB), NS with companion in the core He burning stage (CHeB) and NS with companion at the early AGB stage (EAGB). As discussed above, mostly, the population with CHeB-components may be identified with SyXBs. In the right bottom panel, three sub-populations are over-plotted. In Fig. 2, several regions are clearly seen. The most populated region is occupied by the model systems at the settling accretion stage with the mean spin periods 100–10000 s and $L_{\mathrm{x}}\sim 10^{32}-10^{34}$ erg s*-1*. The population of systems with $P^{*}\lesssim 30$ s best seen in the rightmost bottom panel in the plot at the left corresponds to NSs on the Bondi accretion stage. The strip at the bottom of the plot at $P^{*}\sim 10$ ms is formed by disc-accreting millisecond recycled NSs. The vertical strip at $L_{\mathrm{x}}\sim 10^{38}$ erg s*-1* is populated by super-Eddington accreting NSs. The last two types of accreting NSs should not be related to SyXBs, but they are descendants of SyXBs and are shown for completeness.

The right panel of Fig. 2 shows the populations of stars with FGB and CHeB components in the regimes of settling accretion and disc-accretion and the corresponding differential X-ray luminosity functions $dN/d\log L_{\mathrm{x}}$. For the latter systems, $L_{\mathrm{x}}$ corresponds to the outbursts, i.e., represents the maximum of $L_{\mathrm{x}}$.

Note also that for the $\alpha$-formalism of common envelopes, even if we would assume that the entire orbital energy is spent to expel the CE ($\alpha_{\mathrm{CE}}=1$), most of the model binaries would have parameters different from those of observed systems. However, the position of the model systems in the $P^{*}-L_{\mathrm{x}}$ diagram well overlaps with that of observed ones if the $\gamma$-formalism with parameter 1.75 is applied for the CE stage (see also Fig. 2, right panel). Formally, this is equivalent to $\alpha_{\mathrm{CE}}\approx 4$ and may mean that energy sources other than the binary orbital energy are necessary for expulsion of the CE (see also Yungelson & Kuranov (2017) and the discussion of low CE efficiency found in 3D-simulations by Ricker et al. (2018)).

Figure 2 suggests that systems with CHeB donors mostly contribute to the Galactic population of SyXBs; the number of SyXBs with EAGB donors is insignificant because of the short duration of this stage. Generally, the location of the observed SyXBs from Table 1 in this diagram is in agreement with calculations. We also stress here that the intrinsically unstable character of the wind subsonic accretion we discussed above enables the model systems to shift by an order of magnitude in the X-ray luminosity up to the maximum available value for the Bondi-Hoyle wind accretion $L_{\mathrm{x},B}$, which would result in better agreement with observations (see the upper rows of right panels in Fig. 2).

4.3 X-ray luminosity – orbital period distribution

The left panel of Fig. 3 shows the model $L_{\mathrm{x}}-P_{\mathrm{b}}$ distribution of wind-accreting systems. As in the case of $P^{*}-L_{\mathrm{x}}$ plot, the observed systems (shown by open circles and dashed lines) are better reproduced by the model with low-efficiency CE (large $\alpha_{\mathrm{CE}}=4$). A more efficient CE would give rise to model SyXBs with shorter orbital periods than observed.

There is a pronounced feature in the $L_{\mathrm{x}}-P_{\mathrm{B}}$ plot, a region at short orbital periods and low X-ray luminosities, quite densely populated by model systems, which persists for all model parameters. These model sources are produced by NS+FGB systems (see leftmost upper panel of Fig. 3). There can be also a small contribution of high-luminosity long-period NS+CHeB SyXBs (see Fig. 2). However, as discussed above, we do not identify these systems with SyXBs. Presently, no wind-fed LMXBs with such parameters are known. Detection of such sources by future more sensitive X-ray observations would provide an important test of the model.

4.4 NS spin period – orbital period distribution (the Corbet diagram)

The model Corbet diagrams ($P^{*}-P_{\mathrm{b}}$) are presented in the right panel of Fig. 3, for the same model parameters as in Fig. 2. Like Fig. 2, the model predicts appearance of a sub-population of rapidly rotating disc-accreting NSs (recycled millisecond X-ray pulsars), which are descendants of SyXBs. In addition, very slowly rotating low-luminosity NSs are also present in this diagram. Clearly, the observed SyXBs from Table 1 fall within model regions on this plot. However, small number statistics prevents us from obtaining firm conclusions.

5 Discussion and conclusions

In this paper, we present the results of population synthesis of Galactic symbiotic X-ray binaries harbouring a neutron star accreting from the wind of low-mass evolved companion. We have used a modified version of the open population synthesis code BSE supplemented by a detailed treatment of the evolution of rotating magnetized neutron stars. In comparison to our previous studies (Lü et al., 2012; Kuranov & Postnov, 2015), the code was upgraded for an accurate treatment of different regimes of wind accretion (subsonic and supersonic quasi-spherical and disc regimes). This enabled us to reach better consistency of the model with observed characteristics of known SyXBs.

We have found that under the usual assumptions about the binary evolution and NS formation parameters, the model NS spin period – X-ray luminosity (Fig. 2, left panel), X-ray luminosity – orbital period (Fig. 3, left panel) and NS spin period – orbital period (Fig. 3, right panel) distributions are simultaneously in agreement with observations if we assume a low efficiency of matter expulsion in the common envelope stage characterized by the CE efficiency parameter $\alpha_{\mathrm{CE}}=4$, suggesting that other energy sources than only orbital energy of the companions are required to expel the common envelope. This finding is in line with independent indications of low CE efficiency reported in Ohlmann et al. (2016); Ricker et al. (2018) and discussed in Yungelson & Kuranov (2017) for possible progenitors of SNe Ia. Other evolutionary parameters (mass range of progenitors of NSs formed via electron-capture supernovae and the value of the NS kick velocity, initial NS magnetic field distribution and parameters of the NS field decay, etc.) are found to have minor effect on the model properties of SyXBs.

We also stress that the observed X-ray luminosity of SyXBs is better reproduced when the model systems are at high (flaring) stage, with the mass accretion rate onto NS being equal to its maximum possible (Bondi) value $\dot{M}_{\mathrm{B}}$ from the stellar wind of the companion. This fact has a clear physical explanation: while the settling accretion stage (with reduced mass accretion rate $\dot{M}_{\mathrm{x}}=f(u)\dot{M}_{\mathrm{B}}$, see Eq. (9)) seems to be unavoidable to explain the long spin periods of accreting NSs and observed correlations, e.g. in GX 1+4 (González-Galán et al., 2012; Shakura et al., 2012), settling accretion becomes progressively unstable in the radiative cooling regime. That is, at $\dot{M}_{\mathrm{B}}\ll 4\times 10^{36}$ erg s*-1*, the source should demonstrate a flaring behaviour, with the mass accretion rate reaching $\dot{M}_{\mathrm{B}}$ in outbursts. An extreme example of such a behaviour could be supergiant soft X-ray transients (SFXTs) (Shakura et al., 2014), with short duty-cycle of outbursts. However, in the case of SyXBs, the properties of stellar winds are very different from those of early type supergiants in SFXTs, and the outburst duty cycle can be longer (e.g., as in GX 1+4). Therefore, it is more likely to observe a SyXB in the flaring state.

The Galactic number of model SyXBs was calculated using the convolution of the population calculated for an instantaneous star formation burst with the Galactic star-formation rate history (see Fig. 1). For the adopted Galactic SFR, this number is slightly below 40 per $10^{10}$ $M_{\odot}$ and was almost constant during the last 10 Gyr. Of all initial binary parameters, it is mostly affected by the shape of the initial binary mass ratio distribution and decreases by factor of a few when changing from a flat mass ratio distribution to that peaked at equal initial masses of binary components. In our model, about 80% of the NS components of SyXBs are formed by electron-capture supernovae and the rest – by the iron-core collapse SN. The number of SyXBs with the post-accretion-induced-collapse NSs is found to be vanishingly small. We expect that the majority of wind-emitting components in SyXBs should be core-helium burning stars.

The computed number of currently existing Galactic SyXBs (Fig. 1) is found to be substantially lower than predicted in our previous studies and is more consistent with the low number of observed objects, though, still exceeding the latter. We should note, however, after Masetti et al. (2012), that identification of new SyXBs requires a very precise (smaller than a few arc sec) localisation of X-ray sources in order to find their optical counterparts, which is a hard task in the crowded star fields, especially in the Galactic bulge.

The model also predicts the existence of very slowly rotating NSs in SyXBs, which spin down at the low-luminosity settling accretion stage. Detailed studies of the X-ray properties of Galactic SyXBs can be used as an independent instrument to test the evolution of binary stars in the Galaxy. We expect that the forthcoming all-sky X-ray survey like the SRG-eROSITA mission (Predehl et al., 2016, 2018) will be able to discover a dozen of Galactic SyXBs similar to the objects described in the present paper.

Our calculations predict the existence of a population of wind-fed LMXB also subject to settling accretion, which hosts FGB stars in the systems with orbital periods below about 200 – 300 days. This range of periods is shorter than the one typical for “classical” symbiotic binaries. Most of them are also X-ray dim sources with L_{\mathrm{x}}\ {\raise-2.15277pt\hbox{\buildrel<\over{\sim}}}\ 10^{32}\,\mathrm{erg~{}s^{-1}}, lower than $L_{\mathrm{x}}$ of typical known SyXBs. We expect that these binaries in the course of evolution will form common envelopes upon Roche lobe overflow by the optical components and merge or turn into NS+WD systems.

Acknowledgements

We thank the anonymous referee for useful remarks. This study was partially supported by RFBR grant No. 19-02-00790. The work of AGK is supported by RSF grant No. 14-12-00146 (population synthesis calculations). The work of KAP (Section 2 and Appendix A) is partially supported by RFBR grant 18-502-12025. KAP and AGK also acknowledge the support from the Program of development of M.V. Lomonosov Moscow State University (Leading Scientific School ’Physics of stars, relativistic objects and galaxies’).

Appendix A Effect of the orbital eccentricity on the NS equilibrium period at the settling accretion stage

Here we describe the effect of orbital eccentricity on the value of the equilibrium NS period at the settling accretion stage. Because of the orbital eccentricity, the spin-up and spin-down torques applied to the NS change during the orbital motion and should be averaged over the orbital period. In the Keplerian two-body problem, the separation between the barycentres of the binary components is

[TABLE]

where $p=a(1-e^{2})$ is the orbital semilatus rectum, $a$ is the semi-major axis of the orbit, $e$ is the orbital eccentricity, $\theta$ is the true anomaly. The stellar wind velocity is assumed to be spherically symmetric and centered on the optical star barycentre. If the radius of the optical star is $R_{\mathrm{o}}$, the stellar wind velocity as a function of distance from the star, normalized to the escape velocity at the optical star’s photosphere, $v_{\mathrm{esc}}=\sqrt{2GM_{\mathrm{o}}/R_{\mathrm{o}}}$, can be written as

[TABLE]

For example, for radiative-driven accelerating winds of early type stars

[TABLE]

where $\alpha=0.7...3.8$ and $\beta=0.7...2.0$ are obtained from recent numerical simulations (Vink, 2018).

For decelerating winds from evolved red giants – optical components of SyXBs, we follow the prescription given in Lü et al. (2006), which can be conveniently written as a continuous function

[TABLE]

In the case of an eccentric orbit, the specific angular momentum of captured wind matter is determined by the tangent component of the relative orbital velocity of the NS,

[TABLE]

while the gravitational capture Bondi radius $R_{\mathrm{B}}$ depends on the module of the sum of the orbital velocity vector $\bf{v}=\bf{v}_{n}+\bf{v}_{r}$ and the wind velocity vector $\bf{v}_{\mathrm{w}}(r)$ that has the radial component only:

[TABLE]

Here $v_{r}$ is the radial orbital velocity component:

[TABLE]

In the case of an eccentric orbit, the spin-up torque in Eq. (11) takes the form:

[TABLE]

and therefore (neglecting the term $z/Z\ll 1$) the equilibrium NS period averaged over the orbit can be determined from the condition $\langle dI\omega^{*}/dt\rangle=0$:

[TABLE]

Here $\langle\ldots\rangle=2\int_{0}^{\piup}\ldots d\theta$.

If we take into account mass conservation of spherically symmetric stellar wind of an optical star losing mass at the rate $\dot{M}_{\mathrm{o}}$, i.e.,

[TABLE]

and keep in mind that the Alfvén radius $R_{\mathrm{A}}\sim\dot{M}_{\mathrm{B}}^{-2/7}$, we arrive at the ratio of the NS equilibrium period in an eccentric orbit to that in the circular orbit ($e=0$, $P_{\mathrm{eq}}^{0}=P_{\mathrm{B}}(R_{\mathrm{A}}/R_{\mathrm{B}})^{2})$:

[TABLE]

Here the relative NS orbital velocity in the circular orbit $r=a$ is

[TABLE]

and $v_{w,0}=\sqrt{\frac{2GM_{\mathrm{o}}}{R_{\mathrm{o}}}}f(a)$ is the wind velocity for the circular orbit.

In Fig. 4 we plot the reduction factor $F(e)$ as a function of the orbital eccentricity $e$ for several values of $a/R_{\mathrm{o}}$, $M_{\mathrm{o}}/M$ and different stellar wind velocity laws [Eq. (17) and Eq. (18) for accelerating and decelerating stellar winds, respectively]. It is seen that in eccentric binaries the NS equilibrium period during the quasi-spherical settling accretion stage can be reduced by a factor of 10-100. This reduction is generally stronger for accelerating radiation-driven winds (left panel of Fig. 4) and is important for high-mass X-ray binaries (Postnov et al., 2018). For decelerating winds from late-type giants (this paper) this factor is less than 10 for orbital eccentricities below about 0.6.