The frequency of Kozai-Lidov disc oscillation driven giant outbursts in Be/X-ray binaries

TL;DR

This paper investigates how the Kozai-Lidov mechanism causes eccentricity growth in Be-star discs, leading to giant outbursts in Be/X-ray binaries, with models matching observed outburst timescales.

Contribution

It demonstrates the conditions under which Kozai-Lidov oscillations induce giant outbursts, linking disc properties and binary parameters with observed phenomena.

Findings

KL instability depends on disc inclination and aspect ratio.

Model predictions align with observed outburst timescales.

Flaring discs are more unstable in shorter orbital period binaries.

Abstract

Giant outbursts of Be/X-ray binaries may occur when a Be-star disc undergoes strong eccentricity growth due to the Kozai-Lidov (KL) mechanism. The KL effect acts on a disc that is highly inclined to the binary orbital plane provided that the disc aspect ratio is sufficiently small. The eccentric disc overflows its Roche lobe and material flows from the Be star disc over to the companion neutron star causing X-ray activity. With N-body simulations and steady state decretion disc models we explore system parameters for which a disc in the Be/X-ray binary 4U 0115+634 is KL unstable and the resulting timescale for the oscillations. We find good agreement between predictions of the model and the observed giant outburst timescale provided that the disc is not completely destroyed by the outburst. This allows the outer disc to be replenished between outbursts and a sufficiently short KL oscillation timescale. An initially eccentric disc has a shorter KL oscillation timescale compared to an initially circular orbit disc. We suggest that the chaotic nature of the outbursts is caused by the sensitivity of the mechanism to the distribution of material within the disc. The outbursts continue provided that the Be star supplies material that is sufficiently misaligned to the binary orbital plane. We generalise our results to Be/X-ray binaries with varying orbital period and find that if the Be star disc is flared, it is more likely to be unstable to KL oscillations in a smaller orbital period binary, in agreement with observations.

Equations

$$\sqrt{1-e^2}\cos i\approx const.$$

$$\frac{{\tau }_{\mathrm{KL}}}{P_{\mathrm{orb}}}\approx \frac{M_1+M_2}{M_2}\frac{P_{\mathrm{orb}}}{P_{\mathrm{p}}}(1-e_{\mathrm{b}}^2{)}^{\frac{3}{2}},$$

$$\varphi ={\tan }^{-1}\left(\frac{l_y}{l_x}\right).$$

$$\nu =\alpha {\left(\frac{H}{R}\right)}^2{R}^2\Omega ,$$

$$\frac{H}{R}={\left(\frac{H}{R}\right)}_0{\left(\frac{R}{R_0}\right)}^s,$$

$${\left(\frac{H}{R}\right)}_{\mathrm{crit}}=\frac{\sqrt{M_{2}M}}{M_1}\sqrt{\frac{R^3}{a^3}}$$

$${\tau }_{\nu }=2.2{\left(\frac{\alpha }{0.3}\right)}^{-1}{\left(\frac{(H/R{)}_{\mathrm{t}}}{0.08}\right)}^{-2}{\left(\frac{M_1}{18{\mathrm{M}}_{\odot }}\right)}^{-\frac{1}{2}}{\left(\frac{R_{\mathrm{t}}}{50{\mathrm{R}}_{\odot }}\right)}^{\frac{3}{2}-2s}yr.$$

$$\Sigma ={\Sigma }_0{\left(\frac{R}{R_{\mathrm{in}}}\right)}^{-n+s+1}$$

$$n=3s+2$$

$$⟨{\tau }_{\mathrm{KL}}⟩=\frac{{\int }_{R_{\mathrm{in}}}^{R_{\mathrm{t}}}\Sigma R^3\sqrt{\frac{G{M}_1}{R^3}}dR}{{\int }_{R_{\mathrm{in}}}^{R_{\mathrm{t}}}{\tau }_{\mathrm{KL}}^{-1}\Sigma R^3\sqrt{\frac{G{M}_1}{R^3}}dR}$$

$$P\lesssim P_{\mathrm{crit}}\approx 150{\left(\frac{\eta }{50/95}\right)}^{-3/2}{\left(\frac{(H/R{)}_0}{0.06}\right)}^{-3}\mathrm{day},$$

Full text

††pagerange: The frequency of Kozai–Lidov disc oscillation driven giant outbursts in Be/X-ray binaries–The frequency of Kozai–Lidov disc oscillation driven giant outbursts in Be/X-ray binaries††pubyear: 2019

1 Introduction

Be stars are rapidly rotating early type stars with a viscous thin, Keplerian disc (Lee, Osaki & Saio, 1991; Porter, 1996; Hanuschik, 1996; Hummel, 1998; Quirrenbach et al., 1997; Porter & Rivinius, 2003; Rivinius, Carciofi & Martayan, 2013; Okazaki, 2016). Material is added to the inner parts of the disc from the Be star (Cassinelli et al., 2002) and thus the disc is a decretion disc, rather than an accretion disc (Pringle, 1991). The stellar material spirals outwards through the disc due to the action of viscosity that is thought to be driven by the magnetorotational instability (Balbus & Hawley, 1991).

Be/X-ray binaries typically consist of a Be star and a neutron star in an eccentric orbit binary. The orbit became eccentric, and most likely inclined with respect to the spin of the Be star, when the neutron star formed in a supernova explosion. A slight asymmetry to the explosion caused a kick on the newly formed neutron star (Brandt & Podsiadlowski, 1995; Martin, Tout & Pringle, 2009). The tidal torque exerted by the neutron star companion truncates the Be star disc (Okazaki et al., 2002; Hayasaki & Okazaki, 2004; Reig, 2007). This is observationally confirmed since the disc is more dense when the Be star has a binary companion (Zamanov et al., 2001; Cyr et al., 2017).

Some Be/X-ray binaries are transient systems, meaning that they alternate between long quiescent periods and short outbursts, while others are persistent X-ray sources (Okazaki et al., 2002). There are two types of X-ray outbursts that occur through accretion on to the neutron star (e.g. Kuehnel et al., 2015). Type I X-ray outbursts in an eccentric orbit binary occur each orbital period when the neutron star is at periastron and able to capture material from the Be star disc (Okazaki & Negueruela, 2001; Negueruela et al., 2001; Okazaki, Hayasaki & Moritani, 2013). In a low eccentricity binary, type I outbursts may be driven by a disc that becomes eccentric through mean motion resonances with the binary (Franchini & Martin, 2019). Type II (or giant) outbursts are much brighter, last for longer and occur less frequently (Stella, White & Rosner, 1986; Negueruela et al., 1998; Kretschmar et al., 2013; Cheng, Shao & Li, 2014; Monageng et al., 2017).

In this work we focus on the giant outbursts and in particular the Be/X-ray binary 4U 0115+634. This object is one of the first discovered (Giacconi et al., 1972; Whitlock, Roussel-Dupre & Priedhorsky, 1989) and best-studied Be X-ray binaries (e.g. Campana, 1996; Negueruela et al., 1997). This a classical system with well constrained orbital parameters on which many models for the behaviour of Be X-ray binaries are based (e.g. Negueruela & Okazaki, 2001; Okazaki & Negueruela, 2001).

Giant outbursts in this system have been associated with the presence of a tilted disc around the Be star (Reig et al., 2007; Moritani et al., 2011; Martin et al., 2011; Moritani et al., 2013; Kato, 2014). The observed timescale between giant outbursts is around $3\,\rm yr$ (Whitlock, Roussel-Dupre & Priedhorsky, 1989; Negueruela et al., 2001; Reig et al., 2007; Reig & Blinov, 2018). However, sometimes these outbursts are separated by only 1-1.5$\,\rm yr$, and in this case the second outburst is of lower luminosity than the first (e.g. Rouco Escorial et al., 2017, 2019). After major outbursts a decrease in the size of the Be star disc is observed (Reig et al., 2007, 2016).

A sufficiently misaligned disc around a Be star can become highly eccentric (Martin et al., 2014a). This is a result of Kozai–Lidov oscillations that exchange inclination and eccentricity of a misaligned orbit around one component of a binary (Kozai, 1962; Lidov, 1962). In a fluid disc, if the radial sound crossing timescale is shorter than the period of the KL oscillations then the disc can undergo global KL oscillations (Martin et al., 2014b; Fu, Lubow & Martin, 2015a, b, 2017). The eccentricity growth causes the disc to overflow its Roche lobe and transfer material to the companion neutron star (Franchini, Martin & Lubow, 2019). This may be the mechanism behind type II X-ray outbursts observed in Be/X-ray binaries.

In this work we consider the frequency of type II outbursts as a result of KL oscillations of the Be star disc. In Section 2 we first consider the timescale for KL oscillations of a test particle orbit around the Be star in the binary 4U 0115+634. In Section 3 we explore the KL oscillation timescale for a steady state Be star decretion disc around the Be star in the binary 4U 0115+634. We examine parameters for which the disc is KL unstable and find that our model agrees well with the numerical simulation presented in Martin et al. (2014a). The advantage of our analytic approach is that we can examine a large parameter space. In Section 4 we generalise our results to a wider range of parameters and show that giant outbursts driven by KL oscillations are more likely to occur for smaller orbital period binaries. In Section 5 we discuss the implications of our results and we draw our conclusions in Section 6.

2 KL oscillation timescale for a particle orbit

We consider a binary with components of mass $M_{1}$ and $M_{2}$, orbital period $P_{\rm orb}=2\pi/\sqrt{G(M_{1}+M_{2})/a^{3}}$, (where $G$ is the gravitational constant), semi–major axis, $a$, and orbital eccentricity $e_{\rm b}$. The KL effect causes oscillations of the inclination, $i$, and eccentricity, $e$, of a misaligned test particle around one component of a binary system (Kozai, 1962; Lidov, 1962). For an initially circular test particle orbit around one component of a circular orbit binary, oscillations occur when the initial inclination satisfies $39^{\circ}\lesssim i\lesssim 141^{\circ}$. During the oscillations, the component of the angular momentum of the particle that is perpendicular to the binary orbital plane is conserved so that

[TABLE]

The KL oscillation timescale is approximately

[TABLE]

where $P_{\rm p}=2\pi/\sqrt{GM_{1}/R^{3}}$ is the orbital period of the particle that orbits around the Be star with semi-major axis $R$. We consider here first how properties of the binary orbit might change the KL oscillation timescale and then how properties of the particle orbit change the timescale.

2.1 Binary parameters

We consider a Be star binary model based on the parameters of 4U 0115+634. For our standard model we follow Negueruela et al. (2001) and take the Be star to have a mass $M_{1}=18\,\rm M_{\odot}$ and radius $R_{1}=8\,\rm R_{\odot}$ (Vacca, Garmany & Shull, 1996). The neutron star companion has a mass $M_{2}=1.4\,\rm M_{\odot}$ and is at an orbital semi–major axis of $a=95\,\rm R_{\odot}$ with orbital eccentricity $e=0.34$ (Rappaport et al., 1978). The orbital period is $P_{\rm orb}=24.3\,\rm day$. The errors in the orbital properties of 4U 0115+634 are very small. The orbital eccentricity is $0.3402\pm 0.0004$ and the orbital period $P_{\rm orb}=24.309\pm 0.021\,\rm day$ (Rappaport et al., 1978). The masses we take for each component in our model are on the low end of the possible ranges. The mass of the Be star is determined through its spectral type (B0.2Ve) and could be up to about $19.5\,M_{\odot}$ (Vacca, Garmany & Shull, 1996). The mass of a neutron star is in the range $1.1-2.1\,\rm M_{\odot}$ (e.g. Özel et al., 2012).

Fig. 1 shows the change in the KL oscillation timescale (calculated with equation (2), for fixed binary orbital period and for a particle of fixed semi–major axis) as a function of the neutron star mass. The two lines show different Be star masses, $M_{1}=18\,\rm M_{\odot}$ (solid line) and $M_{1}=19.5\,\rm M_{\odot}$ (dashed line). The possible error in the Be star mass does not significantly affect the KL oscillation timescale. However, the error in the unknown mass of the neutron star could lead to a decrease in the KL oscillation timescale of up to a maximum of about $30\%$. We discuss this further in Section 3.3.

2.2 Test particle orbital parameters

We consider here test particle orbits around the Be star with varying orbital properties. The left hand panel of Fig. 2 shows test particle orbits around the Be star initially inclined by $50^{\circ}$ and starting at different orbital separations. The binary orbits in the $x-y$ plane and begins at periastron along the $x$ axis. Since the particle is massless the binary orbit remains fixed. The direction of the angular momentum of the particle is given by the unit vector $\bm{l}=(l_{x},l_{y},l_{z})$ and its phase angle is defined as

[TABLE]

The initial phase angle is $\phi=90^{\circ}$. Particles closer to the Be star have a longer KL oscillation timescale, as predicted by equation (2). The farther out a particle is, the more chaotic its evolution. A particle that begins at a separation of $40\,\rm R_{\odot}$ is ejected from the system. The right hand panel of Fig. 2 shows the analytic prediction for the time of the first peak eccentricity given in equation (2). The blue points show the results of test particle simulations with an initially circular orbit inclined by $i=50^{\circ}$ with $\phi=90^{\circ}$. We note that the oscillation timescale is not uniform because of the eccentricity of the binary (e.g. Naoz et al., 2013). Also, the orbits we consider are not in the hierarchical triple body limit in which the separation of the particle from the Be star is much smaller than the separation of the binary.

An initially eccentric particle orbit has a shorter KL oscillation timescale (Franchini, Martin & Lubow, 2019). In the right hand panel of Fig. 2, the purple points show particles with initial eccentricity $e_{\rm p}=0.2$ and the red points with $e_{\rm p}=0.4$. The initial inclination is again $50^{\circ}$, the phase angle is $\phi=90^{\circ}$ and the argument of periapsis is $180^{\circ}$. For small separation, the KL oscillation timescale is a factor of 4.1 shorter for an initial eccentricity of 0.4 and a factor of 2.7 shorter for initial eccentricity of 0.2.

Fig. 3 shows test particle simulations that begin at different nodal phases, $\phi=0^{\circ}$ and $\phi=90^{\circ}$ for two different initial inclinations. These show that the time to the first peak is relatively insensitive to the initial inclination or the nodal phase angle. We also found that the KL oscillation timescale is rather insensitive to the initial argument of periapsis. The parameters that the timescale is most sensitive to are the semi–major axis of the particle and the initial eccentricity of the particle orbit.

For a circular orbit disc composed of test particles, each particle undergoes KL oscillations and nodal precession on a timescale that most strongly depends on its separation from the Be star. The disc of particles forms a thick structure or swarm on a short timescale. However, the pressure in a hydrodynamical gaseous disc connects different orbital radii of the disc together allowing it to undergo global KL oscillations. The decretion disc can extend to larger radii than those of stable test particle orbits. The timescale of the oscillations depends upon the distribution of material within the disc. In the next Section we explore the KL oscillation timescale for a gas disc.

3 Decretion disc model

In this Section we first consider the Be star decretion disc temperature structure and the surface density profile. We focus on the orbital parameters of the Be/X-ray binary 4U 0115+634. We estimate the global KL oscillation timescale for a Be star decretion disc as a means to calculate the frequency of type II X-ray outbursts. As the disc undergoes a KL oscillation, it becomes very eccentric. However, the inner parts of the disc continue to accrete circular orbit material from the star. The interaction of the circular orbit material being added to an eccentric disc is beyond the scope of this work but will be investigated in the future. Here, we consider a steady state decretion disc model for the Be star disc. This represents the quasi–steady state disc that is truncated by the tidal torque from the companion neutron star.

3.1 Disc temperature structure

The material in the disc orbits the Be star at Keplerian angular frequency $\Omega=\sqrt{GM_{1}/R^{3}}$. The disc extends from the stellar radius, $R_{\rm in}=8\,\rm R_{\odot}$, up to $R_{\rm t}$. The outer radius of the disc, $R_{\rm t}$, is determined by where the disc is tidally truncated by the companion neutron star. The outer edge of the disc for 4U 0115+634 in the smoothed particle hydrodynamical simulations of Martin et al. (2014a) is $R_{\rm t}\approx 50\,\rm R_{\odot}$. The viscosity of the disc is given by

[TABLE]

where $H$ is the disc scale height and $\alpha$ is the Shakura & Sunyaev (1973) viscosity parameter. The value of $\alpha$ in hot Be star discs is likely around $0.3$ (Jones, Sigut & Porter, 2008; Carciofi et al., 2012; Ghoreyshi & Carciofi, 2017; Rímulo et al., 2018). This is typical for a fully ionised accretion disc in which the viscosity is driven by the magneto–rotational instability (Martin et al., 2019). We assume that the aspect ratio is a power law in radius

[TABLE]

where $s$ is a constant and $(H/R)_{0}$ is the disc aspect ratio at $R=R_{0}$. We choose to scale the disc aspect ratio at an orbital radius of $R_{0}=50\,\rm R_{\odot}$. Fig. 8 shows the disc aspect ratio for three values for the power index, $s=0$ (a constant disc aspect ratio), $s=0.25$ and $s=0.5$ (an isothermal disc).

Observationally, the disc aspect ratios of Be stars discs are not well constrained. Wood, Bjorkman & Bjorkman (1997) found the aspect ratio at the inner edge of $\zeta$ Tauri to be 0.04 while other estimates suggest the disc aspect ratio may be much larger farther out in the disc (e.g Porter, 1996; Quirrenbach et al., 1997). Hanuschik (1996) found that the disc flares at large radii. Negueruela et al. (2001) modelled the inner disc edge to be $0.026$ for 4U 0115+634. Due to these uncertainties we have chosen to parameterise our models by the disc aspect ratio at the outer edge since this is more important for the dynamical evolution of the disc. Given the uncertainty in the disc aspect ratio, we consider two different values of $(H/R)_{0}=0.1$ and $(H/R)_{0}=0.06$.

3.1.1 Criteria for the disc to be KL unstable

The KL oscillation of a disc may be suppressed if the value for $H/R$ at the disc outer edge becomes too large (Lubow & Ogilvie, 2017; Zanazzi & Lai, 2017). For a small mass perturber, the critical value at the disc outer edge is given approximately by

[TABLE]

(Lubow & Ogilvie, 2017), where $M=M_{1}+M_{2}$ is the total mass of the binary. For a radius of $R=50\,\rm R_{\odot}$, we find $(H/R)_{\rm crit}=0.11$. For the parameters of 4U 0115+634, we plot the critical aspect ratio as a function of radius in the solid red line in Fig. 4. Since the disc models we consider all have smaller disc aspect ratio in the outer parts of the disc, they are unstable to KL disc oscillations.

The critical disc aspect ratio estimate for the KL mechanism to operate is in agreement with simulations of isothermal discs published in Martin et al. (2014a). We found that the KL mechanism did not operate for $H/R(R_{\rm in})=0.04$, but operated for smaller values considered. This corresponds to disc aspect ratio at the outer edge of about $H/R(R=50\,\rm R_{\odot})=0.1$. Thus the numerical simulations are in good agreement with the analytic estimate in equation (6). The two isothermal disc models shown in Fig. 4 (solid blue and black lines) both have $H/R(R_{\rm in})<0.04$.

3.1.2 Disc viscous timescale

The viscous timescale at the outer edge of the disc is defined as $\tau_{\nu}=R_{\rm t}^{2}/\nu$, where $R_{\rm t}$ is the outer truncation radius for the disc. For typical parameters this is given by

[TABLE]

Fig. 5 shows the viscous timescale as a function of radius for varying $s$. Since the outer parts of the disc are depleted during the outburst, material needs to spread outwards on a timescale shorter than the type II outburst timescale. This requires the disc aspect ratio at the outer disc edge to be $H/R\gtrsim 0.06$.

We have constrained the disc aspect ratio at the outer edge of the disc to be in the range $0.06\lesssim H/R\lesssim 0.11$. Even the shortest possible viscous timescale (for small outer truncation radius and $(H/R)_{0}=0.1$) is only about a factor of $1.5-3$ times smaller than the outburst frequency. Thus, we suggest that the disc cannot be completely depleted during an outburst in order to have such frequent outbursts.

3.2 Surface density of the disc

For a circular orbit binary, test particle orbits cross at a distance of about $0.48\,a=45.6\,\rm R_{\odot}$ (Paczynski, 1977) and this corresponds to the maximum radius at which a cold disc can exist. However a Be star disc may extend to larger radii than this as it is stabilised by pressure effects within it. The eccentricity of the binary orbit works in the opposite direction and leads to a smaller disc (Artymowicz & Lubow, 1994). However, an inclined disc feels a weaker binary torque than a coplanar disc and thus may be larger (Lubow, Martin & Nixon, 2015). Hydrodynamical simulations of such a disc find that the tidal truncation radius varies in time, but for a disc that is initially misaligned by $70^{\circ}$, it is typically around $50\,\rm R_{\odot}$ (Martin et al., 2014a). The tidal torque increases very strongly with radius (e.g. Papaloizou & Pringle, 1977; Martin & Lubow, 2011) and the surface density profile is sharply truncated at the outer disc edge.

In this work, we do not explicitly include the tidal torque from the neutron star that truncates the outer edge of the disc but the surface density profile is truncated at radius $R_{\rm t}$. For the Be star decretion disc the steady state surface density is

[TABLE]

for $R_{\rm in}<R<R_{\rm t}$, where $\Sigma_{0}$ is a constant that is determined by the total mass of the disc and $n$ is a constant. We chose this power law so that the density scales as $\rho\propto\Sigma/H\propto R^{-n}$. The value for $\Sigma_{0}$ changes in time depending on the total mass of the disc (e.g. Bjorkman & Carciofi, 2005; Carciofi et al., 2006). The amount of mass contained in the disc depends upon the accretion rate of material from the Be star to the disc and the mass loss rate from the disc. Since the disc KL oscillation timescale does not depend upon the total mass of the disc, but rather the distribution of mass, we do not make any assumption on the absolute value of the surface density of the disc.

The surface density of a steady state decretion disc satisfies $\nu\Sigma\propto R^{-1/2}$ (e.g. Pringle, 1981; Martin et al., 2011). For the disc to be steady we require the power law indices to be related by

[TABLE]

The more strongly flared the disc vertical structure, the sharper the drop off with radius in the surface density profile of the steady state disc.

Observationally, in the inner parts of Be star discs, the density at the midplane scales as $\rho\propto R^{-n}$ where $n$ is in the range $2-3.5$ for an isothermal disc (e.g. Cote & Waters, 1987; Porter, 1999). The sound speed is $c_{\rm s}=H\Omega$ and thus an isothermal disc has $H/R\propto R^{1/2}$, or $s=0.5$. With equation (9), the steady state isothermal disc has $n=3.5$, in agreement with the observations.

3.3 KL oscillation timescale

If the KL disc mechanism drives giant outbursts, then the KL oscillation timescale for the disc corresponds to the timescale between subsequent type II outbursts. For an initially circular test particle orbit, the eccentricity may go to zero between oscillations. However, because of dissipation within the disc, the eccentricity oscillation is damped and does not go all the way back to zero. Thus, if there is significant remaining material in the disc after an outburst, it may have high eccentricity without filling the Roche lobe and transferring to the neutron star. The KL oscillation timescale is shorter for test particles that start on an already eccentric orbit (see Figure 2). Thus, the KL oscillation timescale estimates in this Section are upper limits to the predicted timescale between type II outbursts.

The timescale for global KL oscillations of a disc is approximated by

[TABLE]

(Martin et al., 2014b). The total mass of the disc does not affect the timescale for KL oscillations. The KL oscillation timescale is sensitive to the size of the disc and the distribution of material within it. We consider the steady state Be star disc surface density solution given in equations (8) and (9).

Fig. 6 shows the KL oscillation timescale for varying disc outer radius for different values of the disc aspect ratio power law, $s$. Note that the timescale does not depend on the scaling of the disc aspect ratio, $(H/R)_{0}$, only the power law, $s$. The larger the disc outer radius, the shorter the timescale. This is because closer to the neutron star perturber, the KL oscillation timescale is shorter (see equation (2)). The timescale is relatively insensitive to the value of $s$, but the smaller its value, the more mass that is located farther out in the disc and thus the shorter the oscillation timescale.

The observed timescale between two subsequent type II outbursts is around $3\,\rm yr$. The KL oscillation timescales predicted by the initially circular orbit disc model are longer than the observed outburst timescales except for very large truncation radius and for a constant disc aspect ratio ($s=0$). We note that the error in the masses of the binary components may lead to a decrease of a few tens of percent in the timescale (see Section 2.1). This is not sufficient to explain this discrepancy. However, we showed in Section 2.2 that the KL oscillation timescale may be up to a factor of a few shorter for an initially eccentric disc. Thus, in order to reproduce the timescale between two subsequent type II outbursts, there must be eccentric orbit material left in the disc after the outburst. This is in agreement with our finding in Section 3.1 that the disc cannot be completely destroyed during the outburst as the viscous timescale in the disc is long.

The isothermal ($s=0.5$) 3D hydrodynamical simulation described in Martin et al. (2014a) had a peak eccentricity at a time of about $22\,\rm P_{\rm orb}\approx 1.5\,\rm yr$. Since this is half of the KL oscillation timescale, $\left<\tau_{\rm KL}\right>/2$, this appears to be consistent with the recurring type II outburst timescale observed, despite the fact that the disc was initially circular. However, the disc was able to undergo the first KL oscillation on an artificially short timescale due to the initial surface density profile which had $\Sigma\propto R^{-1}$. Since there is no source of mass in the disc simulation, the surface density evolves towards a steady state accretion disc with $\nu\Sigma$=const which gives $\Sigma\propto R^{-3/2}$ for an isothermal accretion disc. The isothermal decretion discs described in this work have $\Sigma\propto R^{-2}$. Thus, the initial disc set up had artificially too much mass in the outer parts of the disc and this led to a shorter KL oscillation timescale. A second outburst was not seen in the simulation because there was no material added at the disc inner edge and the resolution of the disc became too poor after the first outburst. In the future we will model the evolution of a disc that includes the accretion of material at the inner edge from the star. Three dimensional simulations would take into account not only the surface density evolution, but non–axisymmetric spiral arms (e.g. Panoglou et al., 2018, 2019) that cannot be included into our one–dimensional approach. However, the advantage of our analytic model is that we can explore a wide range of parameters.

Simulations of accretion discs undergoing KL oscillations previously have shown damped oscillations that decay in time. The inclination of the disc gradually decreases until the inclination is below the critical KL angle (e.g. Martin et al., 2014b; Fu, Lubow & Martin, 2015a, b). However these simulations did not have any source of material to the disc that would keep the inclination of the disc high. We suggest that the Be star disc can continue to undergo KL oscillations provided that the spin of the Be star remains misaligned to the orbital plane and thus provides a source of highly inclined material to the disc. Damping of the spin–orbit misalignment is expected to occur on the same timescale as rotation synchronisation (Hut, 1981; Eggleton & Kiseleva-Eggleton, 2001). The synchronisation timescale (Hurley, Tout & Pols, 2002) is much longer than the lifetime of the Be star except for the smallest orbital period Be/X-ray binaries (Stoyanov & Zamanov, 2009). Thus we expect that a spin–orbit misalignment imparted by the supernova that formed the neutron star to last for the lifetime of the Be/X-ray binary.

4 Disc KL oscillation timescale dependence on binary orbital period

The models presented so far have focused on the orbital parameters of the Be/X-ray binary system 4U 0115+634. We now examine how the mechanism operates for parameters relevant for Be/X-ray binary systems with different orbital period. Fig. 7 shows the KL oscillation timescale for a disc in a system with varying orbital period. We assume that the disc is truncated at $R_{\rm t}=(50/95)\,a$, which is relevant for binary eccentricity $e_{\rm b}=0.34$ and disc inclination of $i=70^{\circ}$. We vary the semi–major axis $a$ such that the orbital period changes in the range $10{\,\rm day}<P_{\rm orb}<500\,\rm day$. Note that smaller binary eccentricities and larger disc inclinations lead to wider discs (Artymowicz & Lubow, 1994; Lubow, Martin & Nixon, 2015; Miranda & Lai, 2015; Brown et al., 2019) and hence shorter KL oscillation timescales. Generally, the KL oscillation timescale is quite insensitive to the orbital period of the binary, within a factor of 2, assuming that the disc has expanded out to the tidal truncation radius. The viscous timescale at the outer edge of the disc increases with binary separation. This viscous timescale exceeds the KL oscillation timescale for longer orbital period binaries.

The flared disc model lines in Fig. 7 are truncated where $H/R(R=R_{\rm t})=(H/R)_{\rm crit}$. For the isothermal disc model, we estimate the critical binary orbital period below which the Be star disc is unstable to the KL mechanism. We require that the disc aspect ratio (Eq. 5 with $s=0.5$) is less than or equal to $(H/R)_{\rm crit}$ at the tidal truncation radius and find

[TABLE]

where $(H/R)_{0}$ is the disc aspect ratio at $R_{0}=50\,\rm R_{\odot}$ and we defined the ratio $\eta=R_{\rm t}/a$. For longer binary orbital period, $P>P_{\rm crit}$, (i.e. higher aspect ratios) the disc is stable against KL oscillations. We have assumed that the disc is wide and viscous enough to reach its tidal truncation radius. If the disc does not expand out to the truncation radius, type II outbursts are not likely to occur. Therefore equation (11) represents a necessary but not sufficient condition for the disc to be KL unstable and thus for type II outbursts to occur.

Note that this critical period should also depend on both the eccentricity and the inclination of the binary orbit but we have only taken that into account through the parameter $\eta$. If the binary is less eccentric, or the disc is more inclined, the size of the disc relative to the binary separation increases, i.e. $\eta$ increases (Artymowicz & Lubow, 1994; Lubow, Martin & Nixon, 2015). This leads to a smaller critical orbital period since the more extended disc has larger H/R in the disc outer parts and is more stable to KL oscillations.

In summary, if the disc has a flared vertical structure, the disc aspect ratio at the outer edge may exceed the critical value required for it to be unstable to KL oscillations in systems with long binary orbital periods. Thus, KL oscillations are more likely in smaller orbital period binaries. This is in agreement with observations that show that type II outbursts are more likely to occur in short period binaries (Cheng, Shao & Li, 2014). The critical orbital period depends on the disc aspect ratio and the disc radial extent, which depends on the eccentricity of the binary orbit, the inclination of the disc and the mass ratio of the binary.

5 Discussion

If KL disc oscillations drive type II outbursts, then any Be star X-ray binary is expected to have type II outbursts providing that two conditions are satisfied. First, the Be star decretion disc must be misaligned to the binary orbital plane by more than the critical angle required for KL oscillations. For a test particle in an initially circular orbit around one component of a circular orbit binary, this is $39^{\circ}$. However this angle varies for thick decretion discs depending upon the properties of the disc (Lubow & Ogilvie, 2017). Secondly, the orbital period of the binary must be short enough to allow the Be star disc to expand sufficiently far to overflow its Roche lobe during a KL oscillation while still being unstable to KL oscillations (see equation (6)). This is in agreement with the observational data provided in Cheng, Shao & Li (2014) that show that Be/X-ray binaries with short orbital period are more likely to undergo type II outbursts. However, there are some observed systems with long orbital periods that exhibit type II outbursts. GRO J1008-57 has an orbital period of $250\,\rm days$ and shows giant outbursts (e.g. Kühnel et al., 2017). For this system to undergo KL oscillations, the disc aspect ratio must be close to constant with radius for the mechanism to operate. Swift J1626.6-5156, with an orbital period of $132\,\rm days$ (e.g. Reig et al., 2011), and 1A 0535+262, with an orbital period of $111\,\rm days$ (Acciari et al., 2011) also show giant outbursts. In order for the KL disc mechanism to drive giant outbursts in these systems the disc must not be very flared. The disc aspect ratio must scale $H/R\propto R^{s}$ with $s\lesssim 0.25$. The KL disc driven giant outburst mechanism could be ruled out through observations of a disc with a low misalignment or a high disc aspect ratio.

A large portion of the Be star disc may be destroyed during the KL oscillation. The material can be re-accreted on to the Be-star (because the disc is highly eccentric), transferred to the neutron star or to a circumbinary orbit (Franchini, Martin & Lubow, 2019), or be ejected from the binary star system. However, the disc continues to be fed at the inner edge from the Be star. The eccentricity of the disc evolves due to the KL mechanism and due to the addition of circular orbit material at the inner edge. The KL oscillations occur on a shorter timescale with smaller amplitude for a larger initial particle eccentricity. Thus in the disc, this leads to a shorter timescale before the next outburst. Since this leaves less time for material to accumulate in the outer parts of the disc, we expect a smaller outburst, as is observed for closely spaced outbursts in 4U 0115+634 (e.g. Reig & Blinov, 2018). Thus, in a disc the KL oscillations depend sensitively on how much material is left over after the outburst, the eccentricity of the material and the decretion rate of material into the disc from the star.

The amount of material that is transferred from the Be star disc to the neutron star depends on several parameters. In particular Franchini, Martin & Lubow (2019) found that for a globally isothermal accretion disc ($s=0.5$) the amount of material around the secondary star was smaller compared to the case with a constant disc aspect ratio ($s=0$). Also, a higher disc aspect ratio as well as higher initial circumprimary disc inclination leads to larger amount of material being transferred. Furthermore the lower the binary mass ratio, the higher the mass transfer to the secondary star. This is relevant to this work because of the extreme mass ratios of these systems. We expect the mass transfer to be very high in Be star X-ray binaries undergoing KL oscillations and with Be star discs filling the Roche lobe.

The inclination of the disc also has a significant affect on the amount of mass transfer. The higher the initial disc inclination, the stronger the eccentricity growth and the more mass that is transferred to the companion. For an inclination of $70^{\circ}$ (but for an equal mass binary), most of the material in the Be star disc was transferred to the companion (see Figure 10 in Franchini, Martin & Lubow, 2019).

6 Conclusions

We have found that KL oscillations of the Be star disc in the Be/X-ray binary 4U 0115+364 operate on a timescale that matches the observed frequency of type II (giant) X-ray outbursts provided that the disc is not completely destroyed during the outburst. A highly inclined disc around a Be star becomes eccentric and overflows its Roche lobe so that material flows on to the companion neutron star. The KL oscillation timescale for a disc depends sensitively on the distribution of material within the disc and the eccentricity of the material. Since the disc is depleted during a giant outburst, material must be replenished between outbursts. The viscous timescale must be short enough for the disc to spread outwards between outbursts. For 4U 0115+634 this requires the disc aspect ratio at the outer disc edge to be $H/R\gtrsim 0.06$. However, the disc aspect ratio at the outer disc edge must be small enough for the KL mechanism to operate, $H/R\lesssim 0.11$. KL mechanism driven type II X-ray outbursts are more likely to occur in shorter period binaries, if the disc is flared.

Acknowledgements

We thank an anonymous referee for a thorough review. We acknowledge support from NASA through grant NNX17AB96G.