Eccentric neutron star disk driven type II outburst pairs in Be/X-ray binaries

Be star X-ray binaries are transient systems that show two different types of outbursts. Type I outbursts occur each orbital period while type II outbursts have a period and duration that are not related to any periodicity of the binary system. Type II outbursts may be caused by mass transfer to the neutron star from a highly eccentric Be star disk. A sufficiently misaligned Be star decretion disk undergoes secular Von Zeipel-Lidov-Kozai (ZLK) oscillations of eccentricity and inclination. Observations show that in some systems the type II outbursts come in pairs with the second being of lower luminosity. We use numerical hydrodynamical simulations to explore the dynamics of the highly misaligned disk that forms around the neutron star as a consequence of mass transfer from the Be star disk. We show that the neutron star disk may also be ZLK unstable and that the eccentricity growth leads to an enhancement in the accretion rate onto the neutron star that lasts for several orbital periods, resembling a type II outburst. We suggest that in a type II outburst pair, the first outburst is caused by mass transfer from the eccentric Be star disk while the second and smaller outburst is caused by the eccentric neutron star disk. We find that the timescale between outbursts in a pair may be compatible with the observed estimates.


INTRODUCTION
Be star X-ray binaries are double systems composed of a neutron star that is accreting material from a companion Be star. The orbital periods are typically tens to hundreds of days. Be stars are hot massive early type main sequence stars with a B spectral type (mass in the range 3 − 20 M ) whose spectrum has shown one or more Balmer H-α emission lines (Porter & Rivinius 2003). They are also rapidly rotating close to their break up velocity (Slettebak 1982). The stars brightness and spectra have been found to vary in time. Their emission is associated with the presence of circumstellar gas in the form of a relatively thin decretion disk (Pringle 1991;Porter & Rivinius 2003). Rajoelimanana et al. (2011) studied the properties of the optical emission from Be stars and found superorbital quasi-periodic variations on timescales of hundreds to thousands of days that are believed to be related to the formation and the depletion of the circumstellar disk.
Some of these systems have a transient nature owing to the interaction between the material from the Be star disk flowing onto and being accreted by the neutron star. Be/X-ray transients typically show two types of outburst (Stella 1986;Negueruela 1998). Type I outbursts are less luminous and occur each orbital period. The occurrence of type I outbursts can be explained either by the eccentric nature of the binary orbit that allows the neutron star to rip off material from the Be star disk close to periastron passage Okazaki et al. 2013) or an eccentric disk in systems with small binary eccen-tricity Martin & Franchini 2021).
Giant or type II outbursts have luminosities close to the Eddington limit and no periodic modulation. In general, type II outbursts last several binary orbital periods (Kretschmar et al. 2012). Reig & Blinov (2018) found that type II outbursts often occur in pairs with the second one being less luminous than the first. Type II outbursts may be explained by the presence of a highly misaligned Be star disk that becomes eccentric due to the Von Zeipel-Lidov-Kozai (ZLK) mechanism (von Zeipel 1910;Kozai 1962;Lidov 1962) and transfers mass through Roche lobe overflow onto the neutron star (Martin et al. 2014b,a;Fu et al. 2015a;. The ZLK oscillations and mass transfer timescales match the observed intervals between type II outbursts provided the binary orbital period is 150 days and the Be star disk is flared . As material is transferred from the Be star disk to the neutron star, a compact disk forms around the neutron star (Hayasaki & Okazaki 2004, 2006. Because of the inclination of the Be star disk, the neutron star disk is likely misaligned and initially eccentric . This idea is supported by the detection of quasiperiodic oscillations (QPOs) in some of the outbursts that are believed to be associated with the formation of an accretion disk around the neutron star (Finger et al. 1996;Wilson et al. 2002).
In this paper we extend these previous works to explore the dynamics of the accretion disk that forms around the neutron star. In Section 2 we present the results of a simulation of a misaligned Be star disk. We show that the disk can form around the neutron star with a sufficiently high inclination to undergo ZLK oscillations itself. In Section 3 we show a second simulation of the neutron star disk evolution with a higher resolution. We find that the ZLK oscillation of the neutron star disk causes a long lasting accretion phase onto the compact object that resembles a type II outburst. We draw our conclusions in Section 4.

MISALIGNED BE STAR DISK SIMULATION
We use the smoothed particle hydrodynamics (SPH) code phantom (Price & Federrath 2010;Lodato & Price 2010) to model the system composed of a Be star surrounded by a highly misaligned accretion disk with initial inclination i = 70 • and a companion neutron star. The binary components are modelled as sink particles with masses of M * = 18 M and M ns = 1.4 M , respectively. We choose a binary separation of a = 95 R that corresponds to an orbital period of P b = 24 days and an eccentricity of e = 0.34.
These values correspond to the orbital parameters of the most active Be/X-ray binary transient 4U 0115+63. This is one of the first discovered (Giacconi et al. 1972;Whitlock et al. 1989) and extensively observed (Campana 1996;Negueruela et al. 1997) Be/X-ray binaries. Thanks to its very well constrained orbital parameters, many models for Be/X-ray binaries are based on this source . Giant outbursts in this system have been associated with the presence of a tilted disk around the Be star (Reig et al. 2007;Martin et al. 2011;Kato 2014) and the timescale between subsequent giant outbursts has been estimated to be 3 yrs, even though sometimes these outbursts can occur 1-1.5 yrs apart (Reig et al. 2007;Reig & Blinov 2018). Note that, for this system, the ZLK oscillations of the Be star disk might operate on a timescale that matches the observed 3 yr frequency of type II outbursts, provided that the disk is not completely destroyed during the outburst .
The accretion radii are chosen to be R acc, * = 8 R and R acc,ns = 1 R for the Be star and the neutron star respectively. Particles inside these radii are accreted onto the respective sink particle and their mass and angular momentum added to the sink (Bate et al. 1995). The accretion disk around the Be star extends initially from R in = 8 R to an outer radius R out = 50 R and has a total initial mass M d = 1 × 10 −8 M , distributed with a surface density profile Σ ∝ R −1 (Cote & Waters 1987;Porter 1999). The neutron star does not have an accretion disk at the beginning of the simulation.
The Be star accretion disk initially expands out to the tidal truncation radius. Note that given the moderate binary eccentricities of these systems, the gravitational effect of the neutron star depends on the orbital phase. Therefore the disk radius and the truncation radius are also phase dependent. As expected, the Be star disk is smaller when the neutron star is at periastron and the truncation radius reaches its minimum value, while it spreads outwards when the neutron star is at apastron (Goldreich & Tremaine 1979;Artymowicz & Lubow 1994;.
We use N = 1 × 10 6 particles to resolve the Be star gaseous disk in our simulation. We model the disk viscosity using a Shakura-Sunyaev α viscosity which is implemented via a direct Navier-Stokes viscosity (Flebbe et al. 1994). We therefore use the switches for shocks and to prevent particle interpenetration (Cullen & Dehnen 2010) with α AV between 0.1 and 1 and β AV = 2 to mini-mize the amount of numerical viscosity in the simulation (Meru & Bate 2012). We choose α = 0.2 to model the Be star accretion disk (King et al. 2007;Rímulo et al. 2018;. We use the locally isothermal equation of state for the gas adopted in (Farris et al. 2014) where R * , R ns are the radial distances of the gas particles from the Be star and the neutron star respectively and c s0 is a constant of proportionality. We reduce the temperature of the gas around the neutron star through the constant F by two orders of magnitude at radial distance less than 35 R from the neutron star, similar to Smallwood et al. (2021). This is necessary to better resolve the disk around the neutron star since the viscosity in SPH depends on the number of particles and is therefore prohibitively large for the particles to form a proper disk without being promptly accreted onto the neutron star. Reducing the temperature around the neutron star reduces the disk aspect ratio and therefore the viscosity. We take q = 0 to model a flared circumstellar disk. Martin & Franchini (2019) showed that flared circumstellar disks are unstable against ZLK oscillations for binaries with orbital period up to roughly 100 days. According to the analytical calculations outlined in , since the viscous timescale must be short enough for the disk to spread outwards between outbursts the aspect ratio at the disk outer edge should be ≥ 0.06. At the same time the aspect ratio at the disk outer edge cannot be too large (H/R ≤ 0.11), otherwise the disk would be stable against ZLK oscillations (Lubow & Ogilvie 2017). Therefore we take the disk aspect ratio at the inner edge to be H/R (R in ) = 0.024 which corresponds to H/R (R = 50 R ) = 0.06. Figure 1 shows the eccentricity and inclination evolution of the Be star accretion disk (left) and of the material that forms a disk around the neutron star (right) averaged over the radial extent of each respective disk. At a time of roughly 30 P b the Be star disk undergoes a large ZLK oscillation and its eccentricity increases up to about e = 0.6. The disk fills the Roche lobe and mass is transferred to the neutron star. This is reflected by the Be star disk mass evolution, shown in the bottom left panel of Fig. 1. Since the Be star disk is misaligned with respect to the binary plane, the mass transfer to the neutron star does not necessarily occur in the binary orbital plane, but instead it can take place on a highly misaligned orbit (see upper panels in Fig. 3). An accretion disk with mass ∼ 10 −10 M forms around the neutron star at a time of around 30 P b with initial eccen-  tricity e 0.4 and inclination ∼ 70 • . While the disk is poorly resolved, its evolution is shown in the right panels of Fig. 1. The neutron star disk initially circularises to an eccentricity of about 0.1 over the first few orbital periods. Because the neutron star disk forms above the critical angle required for ZLK oscillations, it is also unstable and its eccentricity increases up to a maximum of about e = 0.6 at a time of around 38 P b . However we can only resolve this newly formed disk for about 12 P b since the disk is accreted rapidly during the ZLK eccentricity peak.
The black lines in the left panels of Figure 2 show the accretion rate (upper panel) and the total accreted mass (lower panel) onto the neutron star as a function of time in units of the binary orbital period. Both quantities are scaled to the newly formed neutron star disk mass measured at t = 36 P b , i.e. M nsdi = 8 × 10 −11 M . At 30 P b the Be star disk undergoes the first large ZLK oscillation causing a high amplitude, longer accretion episode that resembles a Type II outburst. This is due to the mass transfer between the Be star disk and the neutron star, which occurs via Roche lobe overflow (Martin et al. 2014a;. Figure 3 shows the column density maps of the simulation at times of 29, 33, 38 P b from the upper to the lower panel in the x-y (left panels) and x−z plane (right panels). We can see the neutron star disk formation as a consequence of the first ZLK oscillation of the Be star disk in the upper panel. The neutron star disk becomes almost circular in the middle panel (consistent with the upper right panel of Fig. 1). The neutron star disk eccentricity increases again due to the ZLK mechanism at around 38 P b (lower panel).

ECCENTRIC NEUTRON STAR DISK SIMULATION
Previous hydrodynamical simulations performed by Hayasaki & Okazaki (2004, 2006 investigated the accretion flow around the neutron star induced by both Roche lobe overflow and from the pericentre passage. They found that a time persistent accretion disk forms around the compact object in both cases. However in their simulations the mass flow onto the neutron star is assumed to be coplanar with the binary orbital plane.
Even though we tuned the resolution around the neutron star in the simulation by changing the equation of state, this was not sufficient to resolve accurately the neutron star disk. Thus, in order to understand the neutron star disk dynamics, we now consider a second simulation in which we take the same binary parameters but have no disk around the Be star and an initially circular accretion disk inclined by 70 • around the neutron star modelled with N = 500, 000 particles. The initial disk mass is M d = 1 × 10 −8 M distributed with a surface density profile Σ ∝ R −3/2 between R in = 3 R and R out = 10 R . Note that, since the disk mass is very small compared to the mass of the binary the choice of its initial value does not influence the dynamics of the neutron star disk. In addition, we scale the accretion rate and the total accreted mass to the initial neutron star disk mass. The accretion radii of the two sinks are the same as in the previous simulation. The disk aspect ratio is H/R (R in ) = 0.01 and the viscosity parameter is α = 0.1. We modelled the neutron star disk using a global isothermal equation of state, i.e. we use Eq. 1 with q = 0 and F = 1 in the entire domain. We chose the initial inclination and the disk radial extent based on the characteristics of the neutron star disk that forms from the Be star disk ZLK oscillations in the previous simulation (see right panel of Fig. 1). The timescale for ZLK oscillations depends on the disk radial extent and on the distribution of mass within it (Martin et al. 2014a;. In particular, larger disk outer radii and higher mass contents in the outer regions lead to shorter timescales. The red lines in Figure 2 show the accretion rate (upper panel) and the accreted mass (lower panel) onto the neutron star in this second simulation. These quantities are both scaled to the initial mass of the neutron star disk in this second simulation, i.e. 10 −8 M . We start the red lines at t = 36 P b since by this time the first ZLK oscillation of the Be star disk ceases to cause mass transfer to the neutron star disk. The dynamics of the neutron star disk are dominated by the inflow of material until this time. Furthermore the material in the disk at this time is on almost circular orbits (see upper right panel in Fig. 1), which corresponds to the initial condition in our second simulation.
The neutron star disk is initially circular but rapidly, within the first few binary orbits, reaches an eccentricity of e = 0.8 while its inclination drops to ∼ 39 • due to the ZLK mechanism. We find this eccentricity increase to lead to an episode of enhanced accretion rate onto the neutron star, represented by the peak at 40 P b , shown by the red line in the upper right panel of Fig. 2. The characteristics of this second peak are compatible with type II outbursts since the accretion rate increases by more than an order of magnitude and the peak lasts for many orbital periods. This second outburst amplitude is about a factor of more than 2 smaller than the first outburst seen in the Be star disk simulation. The relative magnitude of the two outbursts essentially depends on how much mass is directly accreted with compared to the mass that is transferred and forms the disk. However, since we take the initial neutron star disk mass from the poorly resolved material in the Be star disk simulation, what we show is likely a lower limit. If we were able to resolve the neutron star disk better in the first simulation we would expect a higher disk mass and therefore a larger second outburst.
The timescale for the neutron star disk ZLK oscillation depends on the disk viscosity and on the distribution of the gas within the disk. If the disk is viscous enough to expand rapidly, the type II outburst occurs after a few binary orbits because the disk outer parts become ZLK unstable fast. The timescale for global ZLK oscillations can be written as (Martin et al. 2014a; (2) where p is the surface density power law index, M is the total mass of the binary and M 1 , M 2 are the Be star or neutron star disc masses depending on which disc is undergoing ZLK oscillations.
The neutron star disk in the previous simulation (see Section 2) forms with an outer radius of approximately 13 R . The corresponding global ZLK timescale is, assuming a surface density profile Σ ∝ R −3/2 , τ ZLK 11 P b (Martin et al. 2014a;. We note that the surface density profile of the neutron star disc might deviate from the initial power law as a result of the interaction with the companion and of its own viscous evolution. A steeper (shallower) profile leads to a longer (shorter) ZLK timescale (see Eq. 2). The eccentricity oscillation shown in the upper right panel of Figure 1 occurs at around ∼ 3 P b which is slightly shorter than the predicted τ ZLK /2. Since the disk is initially eccentric, the ZLK period is indeed expected to be shorter than the analytical prediction .
The neutron star disk in this second simulation reaches instead the peak eccentricity at 6 P b which is consistent with the theoretical estimate of the ZLK timescale for an initially circular disk based on its outer radius and surface density profiles, i.e. τ ZLK 11 P b . The recurrence time between subsequent type II outbursts has been estimated to be around 3 yrs from observations (Whitlock et al. 1989). This periodicity has also been confirmed numerically by . However Reig et al. (2007) showed that this three-years cycle might be interrupted by two close, 1-1.5 yrs apart, type II outbursts When this occurs, after the second outburst the disk might be completely dissipated (Reig & Blinov 2018). For the system 4U 0115+63 this translates to ∼ 15 − 23 P b between the two outbursts.
The timescale between the two outbursts corresponds to the time to form the neutron star disk (about half the ZLK timescale of the Be star disk) plus half the ZLK oscillation timescale of the neutron star disk in the second simulation. The analytically estimated value is consistent with the timescale we show in the upper panel of Figure 2, i.e. ∼ 12 P b . This is slightly shorter than the expected time for 4U 0115+63. However, since the ZLK period is very sensitive to the disk surface density profile, a slightly steeper profile, e.g. Σ ∝ R −2 , leads to a longer ZLK timescale. The estimated timescale for the ZLK oscillation of the neutron star disk for such profile is τ ZLK /2 7 P b . This added to the ZLK timescale of the Be star disk, results in a longer interval between outbursts, closer to the 1 yr interval estimated for this source.
Furthermore, the circularization radius of the material coming from the Be star disk, i.e. the neutron star disk outer radius, might also depend on the viscosity, equation of state and mass distribution in the simulation. In particular, smaller outer radii lead to longer ZLK oscillation timescales. Therefore in order for the type II outbursts to be further apart we expect a steeper surface density profile or a much smaller outer radius of the neutron star disk. The timescale for the first ZLK oscillation peak is also longer with a smaller initial disk inclination or a smaller disk aspect ratio (Fu et al. 2015b).
Thus, while our simulation shows a timescale that is a little short for the specific case of 4U 0115+63, we suggest that the neutron star disk might form with a lower inclination than 70 • .
Multiple type II outbursts have been also observed in the X-ray and optical light curves of the two very well studied sources 1A 0535+262 (Rosenberg et al. 1975) and SAX J2103.5+4545 (Hulleman et al. 1998).
1A 0535+262 has an inferred orbital period P b = 111 days (Finger et al. 1994) and underwent the first type II outburst after 10 yrs of quiescence in 2005 (Smith et al. 2005), followed by another one in 2009 (Caballero et al. 2010) and a third one 2 yrs later in 2011 (Monageng et al. 2017). The interval between the 2005 and 2009 outbursts corresponds to ∼ 13 P b for this particular source and is consistent with the timescale inferred from our simulations while the 2011 outburst occured after ∼ 7 P b .
The source SAX J2103.5+4545 has one of the shortest known Be/X-ray binaries orbital period at P b = 12.68 days (Baykal et al. 2000). In the optical band, this source showed two type II outbursts 2.7 yrs apart between 2007 and 2010 and then another pair 0.7 yrs apart (Camero et al. 2014). As these timescales translate to ∼ 75 and ∼ 19 P b , our model could in principle only explain the latest pair of outbursts.

CONCLUSIONS
In this paper we have investigated the dynamics of a highly misaligned disk surrounding the neutron star in a Be/X-ray binary system using 3D SPH simulations. . Column density plots of the Be star and neutron star disks. The binary components are shown in green and their size is determined by their accretion radii. The view is of the x-y plane (i.e. the binary orbital plane) (left panels) and of the x − z plane (right panels). The color scale spans about three orders of magnitude in density and is the same for all the plots. The first, second and third row correspond to 29, 33, 38 P b .
We first ran a simulation of a Be star disk and showed that during ZLK oscillations of eccentricity and inclination, mass transfer produces a highly inclined eccentric neutron star accretion disk. We were unable to obtain a good enough resolution of the neutron star disk by forming it in this way and so we therefore used the parameters of this newly formed neutron star disk as the initial conditions for a second simulation involving only the neutron star disk.
We have shown that the dynamics of the two disks can lead to type-II outburst pairs. The first type II outburst is triggered by the transfer of Be star disk material to the neutron star. The subsequent excitation of the neutron star disk eccentricity, as a result of the ZLK mechanism, causes another large increase of accretion rate onto the compact object. We interpret this as a possible second type II outburst. We finally find the timescale between the two giant outbursts to be ∼ 12 P b , compatible with the 1−1.5 yrs interval observed in 4U 0115+63 (Reig et al. 2007;Reig & Blinov 2018) and possibly also with outbursts pairs in the sources 1A 0535+262 and SAX J2103.5+4545. A more thorough exploration of the parameter space for this mechanism to occur in other Be/X-ray binaries is deferred to a future work.
It is worth noting that whether the Be star disc driven type II outburst occurs as a single event or produces a pair of outbursts does not influence the occurrence rate of type II outbursts in general, as the pair is driven by the same mass transfer mechanism and the neutron star disc driven outburst does not influence the Be star accretion disc. We can therefore say that the occurrence rates of type II outbursts (or outburst pairs) driven by ZLK oscillations are the rates estimated previously in  and depend on the discs viscosity and temperature, that essentially regulate the mass distribution within the disc, as well as the timescale for Be star disc refilling (e.g. Suffak et al. 2022).
We finally note that the resolution in the first simulation was not high enough to completely resolve the second outburst due to the eccentric neutron star disk. This might due to the very large viscosity caused by the small number of particles in the disk coupled with the fact that type I outbursts continue to occur because of the presence of the Be star disk. The improvement of the resolution of the disk around the neutron star using more advanced numerical techniques will be the subject of a future work.
We warmly thank the referee for very useful comments that improved the quality of the manuscript. AF acknowledges financial support provided under the European Union's H2020 ERC Consolidator Grant "Binary Massive Black Hole Astrophysics" (B Massive, Grant Agreement: 818691). RGM acknowledges support from NASA through grant 80NSSC21K0395. We acknowledge the use of SPLASH (Price 2007) for the rendering of the figures.

DATA AVAILABILITY
The SPH simulations results in this paper can be reproduced using the phantom code (Astrophysics Source Code Library identifier ascl.net/1709.002). The data underlying this article will be shared on reasonable request to the corresponding author.