Retrograde discs around one component of a binary are unstable to tilting

With hydrodynamic simulations we show that a coplanar disc around one component of a binary can be unstable to global tilting when the disc orbits in a retrograde direction relative to the binary. The disc experiences the largest inclination growth relative to the binary orbit in the outermost radii of the disc, closest to the companion. This tilt instability also occurs for test particles. A retrograde disc is much larger than a prograde disc since it is not tidally truncated and instead spreads outwards to the orbit of the companion. The coplanar retrograde disc remains circular while a coplanar prograde disc can become eccentric. We suggest that the inclination instability is due to a disc resonance caused by the interaction of the tilt with the tidal field of the binary. This model is applicable to Be/X-ray binaries in which the Be star disc may be retrograde relative to the binary orbit if there was a sufficiently strong kick from the supernova that formed the neutron star companion. The accretion on to the neutron star and the resulting X-ray outbursts are weaker in the retrograde case compared to the prograde case.


INTRODUCTION
Be/X-ray binary systems consist of a high mass Be-type main sequence star with a compact object companion.The companion is often a neutron star (Reig 2011), but may also be a black hole (Munar-Adrover et al. 2014).The Be star rapidly rotates (Slettebak 1982;Porter 1996;Porter & Rivinius 2003) and possesses a viscous Keplarian decretion disc of material flowing out from the star (Pringle 1991;Lee et al. 1991;Hanuschik 1996;Carciofi 2011).Observations show that decretion discs that are tilted relative to the orbital plane of the binary are common (e.g.Hummel 1998;Hughes & Bailes 1999;Hirata 2007;Martin et al. 2011).The tilt is likely due to a kick from the supernova from which the companion formed (Brandt & Podsiadlowski 1995;Martin et al. 2009;Salvesen & Pokawanvit 2020).
Material in the Be star disc may be captured by the companion neutron star and accreted onto its surface.Accretion results in two types of X-ray outbursts (Stella et al. 1986;Negueruela et al. 1998).Type I outbursts occur each binary orbital period.For an eccentric orbit binary, this occurs when the companion approaches periastron (e.g.Negueruela et al. 2001;Negueruela & Okazaki 2001;Okazaki et al. 2002Okazaki et al. , 2013)), while for a circular binary, this can occur when the companion approaches the apastron of an eccentric disc (Franchini & Martin 2019) at which point the companion is closest to the disc material.Type II outbursts are a few orders of magnitude brighter than type I outbursts and they last for several orbital periods (Kretschmar et al. 2013).Type II outbursts can be the result of large eccentricity growth in the circumprimary disc (Martin et al. 2014a;Martin & Franchini 2019;Suffak et al. 2022).This can occur through ★ E-mail: madeline.overton@unlv.eduvon Zeipel-Kozai-Lidov (ZKL, von Zeipel 1910;Kozai 1962;Lidov 1962) oscillations in a highly misaligned disc around one component of a binary (Martin et al. 2014b;Fu et al. 2015;Doğan et al. 2015;Lubow & Ogilvie 2017;Zanazzi & Lai 2017;Franchini et al. 2019).
Observations of the spin rate of neutron stars in Be/X-ray binaries have shown roughly equal numbers of spin up and spin down cases (Klus et al. 2014) and it has been suggested that retrograde accretion may be common (Christodoulou et al. 2017).Retrograde circumprimary discs may also occur in merging supermassive binary black hole systems (e.g.King & Pringle 2006;Nixon et al. 2011;Pereira et al. 2019).However, we are concerned with a time far before the merger phase where the companion drifts very slowly.There is also likely to be a circumbinary disc.Most simulations of mergers involving a retrograde disc have concentrated on the circumbinary disc (e.g.Nixon & Lubow 2015;Bankert et al. 2015).Retrograde particle orbits around one component of a binary have been considered previously and have been found to be stable closer to the companion than prograde orbits (Morais & Giuppone 2012) however, the orbits have been restricted to be in the binary orbital plane (e.g.Morais & Namouni 2022).
Motivated by the observations of Be/X-ray binaries, in this Letter, for the first time we explore the dynamics of a retrograde disc around one component of a binary.In section 2 we describe and compare the numerical experiments of Be/X-ray binaries with prograde and retrograde circumprimary discs.We find an inclination instability in retrograde discs that has not been seen before.To aid our understanding, in Section 3 we also investigate the behavior of test particles at key radii in the disc and discuss resonances in retrograde discs.We suggest that the inclination instability can be explained by a retrograde Lindblad resonance acting over the outermost regions of the disc, which arises from interaction of a small tilt of the disc away from the binary orbital plane with the tidal field of the binary.We draw our conclusions in Section 4.

HYDRODYNAMICAL SIMULATIONS
In this section we first discuss the setup of the hydrodynamical disc simulations and then the results of prograde and retrograde disc simulations.

Simulation setup
We use the phantom (Price & Federrath 2010;Lodato & Price 2010;Price et al. 2018) smoothed particle hydrodynamics (SPH) code (Monaghan 1992;Bate et al. 1995;Price 2012) to model prograde and retrograde discs in Be/X-ray binary star systems.The binary system consists of a  ★ = 18 M ⊙ Be star and a  NS = 1.4 M ⊙ neutron star companion, each represented by a sink particle.Any particles that move inside the sink radius are accreted and their mass and angular momentum are added to the sink particle (Bate et al. 1995).The Be star radius is  * = 8 R ⊙ and the neutron star radius is  NS = 1 R ⊙ , therefore, we do not model the inner region of the neutron star's accretion disc.The binary orbit is initially circular with a separation of  = 95 R ⊙ , which gives an orbital period of  orb = 24.3days.These parameters are motivated by Be/X-ray binaries with nearly circular orbits such as 2S 1553-542, which has an orbital eccentricity  < 0.09 and an orbital period of 30.6 days (Reig 2007) as well as the Be/X-ray binary transient 4U 0115+63 which has an orbital period of 24.3 days (Giacconi et al. 1972;Rappaport et al. 1978).
We initialize a circular disc of total mass  di = 10 −8 M ⊙ around the Be star with 5 × 10 5 particles.The gas particles are initially in Keplerian orbits with a small correction due to pressure effects.The disc initially extends from  = 8 to 30 R ⊙ .The surface density of the disc initially follows a profile Σ = Σ 0 ( ★ /) 3/2 [1 − ( ★ /) 1/2 ], where Σ 0 is a scaling constant.In the vertical  direction, the particles are initially distributed with a particle density  ∝ exp(− 2 /(2 2 )), where  is the disc scale height.The later disc evolution is not very sensitive to the initial surface density profile because the surface density evolves significantly during the simulation.The simulations are independent of the total disc mass and, by extension, Σ 0 since the disc mass is much smaller than the binary mass, and we do not include self-gravity of the disc.For example, if the initial total disc mass was increased by some factor while keeping all other parameters fixed, then the surface density profiles would increase by the same factor at all times.This property holds because we apply a locally isothermal equation of state so that the disc temperature is independent of density.However, the disc density would have an effect on the predicted observables, including the accretion rates and the resulting X-ray flux from outbursts.
We adopt a locally isothermal equation of state of the form

𝑞
(1) (Farris et al. 2014) with  = 0.5, where  and  NS are the distances from the primary and secondary objects, respectively.The disc temperature scales with  ∝  2 s .The constant of proportionality  s0 is chosen to give the disc aspect ratio, /, at the inner disc radius ( =  * = 8 R ⊙ ) to be 0.01.This sound speed gives  s ∝  − close to the Be star and  s ∝  − NS close to the neutron star.The sound speed is  s = Ω, where the Keplerian frequency is Ω = √︁   * / 3 and therefore the circumprimary disc has an nearly constant disc aspect ratio with radius of / ≈ 0.01.These choices give an approximate initial mean smoothing length over the disc scale height of ⟨ℎ⟩ / = 0.9.The Shakura & Sunyaev (1973) viscosity parameter  ≈ 0.3 in Be star discs (e.g.Rímulo et al. 2018;Ghoreyshi et al. 2018;Martin et al. 2019).We adapt the SPH artificial viscosity according to the recommended procedure given by Lodato & Price (2010).We use an artificial viscosity parameter of  AV = 3.32.We can calculate an approximate Shakura & Sunyaev (1973) viscosity parameter with (e.g.Okazaki et al. 2002).In addition, we take  AV = 2, which correspond to the non-linear component of the artificial viscosity that prevents particle penetration in high Mach number flows (Monaghan 1989).We consider two initial inclinations for the circumprimary disc with respect to the binary orbital plane,  = 0 • (prograde) and 180 • (retrograde).This is the only difference between the setup of the prograde and retrograde cases.We run the simulations for a time of 50 P orb .The binary and circumprimary disc parameters are summarized in Table 1.
There is no addition of mass to the disc during the simulation and therefore we model an accretion disc rather than a decretion disc.However, the outer parts of the disc behave like a decretion disc as they spread outwards (e.g.Pringle 1981;Martin & Lubow 2011).This difference does not play a role in the tilt instability that we discuss.Addition of material would occur close to the Be star in the very inner parts of the disc (e.g.Okazaki et al. 2002;Suffak et al. 2022) that are not affected by the tilt instability that we later describe.The disc mass decreases during the simulation as particles can be accreted on to either star and ejected to form circumbinary material.

Prograde disc
The orbital evolution of a prograde disc has been previously extensively studied (e.g.Lubow 1991Lubow , 1992b;;Murray 1998;Goodchild & Ogilvie 2006), we consider it here just for comparison to the retrograde case.We bin the particles into 100 radial bins evenly distributed between  = 8  ⊙ and  = 100  ⊙ and average the particle properties within each bin.The left panels of Fig. 1 show the surface density and inclination to the binary orbital plane of the prograde disc as a function of radius at different times.The prograde disc viscously expands outwards to a radius that is close to the 3:1 resonance given by (Goodchild & Ogilvie 2006).For our parameters ( b =  NS / * = 0.078), the tidal truncation radius is  res = 44.5  ⊙ .The disc surface density drops off outside of this radius.There is no inclination growth relative to the binary orbital plane in the prograde disc.Fig. 2 shows the column density of the circumprimary disc from different viewpoints at a time  = 41.6  orb .The top left panel shows the prograde disc in the  −  plane, the plane of the binary orbit.The -axis is along the direction of the binary angular momentum.The prograde disc becomes eccentric, reaching a maximum eccentricity of  ≈ 0.2.The disc eccentricity growth is explained by the effects of the 3:1 resonance (Lubow 1991).As the eccentricity in the disc increases, the neutron star is able to capture material during each orbital period when it passes close to the disc apastron (Franchini & Martin 2019).We see evidence of type I outbursts for time  ≳ 20  orb .The top right panel shows the prograde disc in the  −  plane, demonstrating the coplanarity of the disc.

Retrograde disc
The right panels of Fig. 1 show the surface density, inclination relative to the binary orbital plane and nodal phase angle of the retrograde disc.Comparing the top panels of Fig. 1, we see that the retrograde disc expands to larger radii than the prograde disc.This is explained by the lack of ordinary Lindblad resonances in a circular, coplanar retrograde disc (Nixon et al. 2011;Lubow et al. 2015).The interaction times between the orbiting companion and nearby gas in the disc are very different depending on the direction of rotation of the disc.The interaction time combined with the viscous expansion of the disc controls the effect of the binary on the disc.Therefore, there is no tidal truncation acting on the retrograde disc due to the very short companion interaction time.Instead, the retrograde disc expands until it reaches the orbit of the neutron star.At later times in the simulation the surface density drops off with a steeper gradient around  ≈ 70  ⊙ .This indicates there is another truncation mechanism acting on the retrograde disc related to the close interaction of the material with the neutron star.The gas particles can be either accreted on to the neutron star or ejected to circumbinary material due to this interaction (Franchini & Martin 2019).The tidal truncation radius of the prograde disc is smaller than the truncation radius of the retrograde disc.
Over time, the outer radii ( ≳ 50  ⊙ ) of the disc become significantly misaligned to the binary orbital plane.By time  = 25.05  orb the inclination of the retrograde disc reaches a quasi steady state in inclination between approximately 50  ⊙ and 70  ⊙ , where the inclination increases nearly linearly with radius.Beyond 70  ⊙ the inclination continues to increase with radius for all times.At 80  ⊙ , the inclination reaches a maximum of approximately  = 165 • .The bottom right panel of Fig. 1 shows the nodal phase angle (longitude of ascending node) as a function of radius for portions of the disc that are inclined less than 179 • .Once the disc becomes tilted away from the binary orbital plane, it undergoes nodal precession driven by the binary companion (e.g.Larwood et al. 1996;Papaloizou & Terquem 1995a;Lubow & Ogilvie 2000;Bate et al. 2000).Because the disc orbits in a retrograde direction, the nodal precession is prograde with respect to the binary orbit.The nodal phase angle increases with radius before reaching a plateau.This indicates that the disc becomes twisted due to differential precession and the outer parts of the disc precess together.
The bottom two panels of Fig. 2 show the column density of the retrograde disc from a top down (left panel) and edge on (right panel) view of the binary.Unlike the prograde disc case, the retrograde disc remains circular throughout the simulation.Additionally, as visible in the edge-on view, the outer radii begin experiencing inclination growth out of the binary orbital plane and the disc becomes warped.
The disc mass decreases throughout the simulations and particles from the disc may be accreted by the neutron star or ejected into circumbinary material (e.g.Franchini et al. 2019).There is less accretion onto the neutron star in the retrograde case than in the prograde case.In the retrograde case, the average accretion rate onto the neutron star over the entire duration of the simulation is approximately two orders of magnitude lower than in the prograde case.For the retrograde disc, more mass is ejected into the circumbinary material than accreted onto the neutron star by a factor of about 12.In contrast, this factor is only about 0.1 in the prograde case meaning that more material is accreted onto the neutron star than forming circumbinary material.In the prograde case there are accretion streams from the disc through which material is accreted on to the neutron star.These streams do not occur in the retrograde case.
The simulations presented here are independent of the total disc mass.However, the quasi steady state surface density profile seen at later times is affected by the temperature profile.Therefore we have also run a simulation with an alternative value for the  value given in equation (1).Increasing the value of  creates a steeper drop off in the steady-state surface density.With  = 1 we find a similar magnitude of tilting to the main results of this work (where  = 0.5).

RETROGRADE DISC TILT INSTABILITY
In order to better understand the cause of the warping in the disc simulation, we use the -body integrator, Rebound (Rein & Liu 2012;Rein & Spiegel 2015) to examine retrograde test particle orbits around the Be star.We simultaneously simulate test particles at varying initial separation.These are massless test particles which do not interact with each other.The particles are placed in a retrograde Keplarian orbit which is coplanar with respect to the binary orbit.Fig. 3 shows the inclination relative to the binary orbital plane and eccentricity of the test particles at an initial separation of  = 70, 80, and 82  ⊙ respectively.The figure shows that particles interior to 82  ⊙ experience large but chaotic oscillations in inclination and eccentricity.At  ≳ 82  ⊙ , the inclination instability becomes strong enough to eject the test particles.This indicates that the inclination instability is stronger at radii close to the neutron star orbit.The retrograde disc inclination at  = 70 and 80 R ⊙ is also shown in Fig. 3 for comparison.In the case of a disc, material can exist in the exterior region to 82  ⊙ (where test particles are ejected) because there is radial communication through the disc that stabilises the outer parts.The radial communication is likely dominated by the effects of viscosity rather than pressure because we choose a large value ( ≈ 0.3) for the Shakura & Sunyaev (1973) viscosity parameter (e.g.Nixon & King 2016).The vertical instability of particle orbits close to a perturber is indicative of a possible inclination instability of a fluid disc.However, the particle instability is dominated by the effects of incoherent motions that differ between nearby particles.To model a disc tilt, it is important include collective effects such as viscosity and pressure that cause the disc to respond coherently and limit the tilt near the resonance (e.g.Meyer-Vernet & Sicardy 1987).
A disc around one component of a circular orbit binary that undergoes nodal precession in the presence of dissipation moves towards retrograde alignment (e.g.Papaloizou & Terquem 1995b;Lubow & Ogilvie 2000;Bate et al. 2003).Therefore there is a balance between the viscosity that acts to align the disc to retrograde coplanar, and the tilting mechanism.A larger tilt would therefore be expected with a smaller viscosity.
Ordinary Lindblad resonances due to the companion are not present in a retrograde disc.At such resonances we have for azimuthal wavenumber  which is a positive integer and Ω = √︁   * / 3 is always positive, which corresponds to prograde orbits.The angular velocity of the binary companion is Ω b = 2/ orb .Lindblad resonances arise for tilted discs due to their interaction with the tidal potential of the binary.These resonances satisfy the condition that (5) (Lubow 1992a).For  = 1, a retrograde resonance is possible for Ω = −Ω b .At such a resonance a fluid disc can be tilt unstable.Although we simulate an initially coplanar disc, small numerical departures from exact coplanarity of the disc midplane inevitably occur that can be amplified to large values by an exponentially growing instability.
Unfortunately there is a singularity in the companion potential in a model involving an extreme binary mass ratio that occurs at a radius Ω = −Ω b and so the stability analysis in Lubow (1992a) cannot be directly applied.The disc does not extend to the resonance, but the resonance has some width over which it could act.Near resonant effects may then play a role.This provides a possible explanation for the tilting of a retrograde disc as well as the inclination of the retrograde test particle orbits.This explanation should be explored further.An estimate for the width of the resonance can be found from the outer radius of the disc at the time it first experiences inclination growth.In the top right panel of Fig. 1 the line corresponding to 5.01  orb has a small rise in the inclination at the very outer edge (at approximately 65  ⊙ ).At this point, there is no influence from exterior material which may cause the inclination growth, indicating this radius is where the resonance becomes strong enough to influence the disc.

CONCLUSIONS
We have found that a retrograde disc orbiting around one component of a binary is unstable to tilting.With hydrodynamical simulations we have modelled coplanar and circular binaries.For a retrograde disc we find that the disc becomes tilted away from retrograde coplanar with the largest inclination growth relative to the binary orbital plane occurring in the outermost radii.We have examined test particle simulations and seen significant inclination growth in their orbits that increase in magnitude with radius from the primary star.The inclination growth may be due to a retrograde Lindblad resonance occurring at Ω = −Ω b which acts over some width.In a disc, the large tilt in the outermost radii is spread over the radial width of the resonance and propagates through the disc resulting in inclination growth away from the location of the resonance.
This tilt instability has implications for accretion in Be/X-ray binaries.The strong warping in the outer parts of the Be star disc means that there is less material in the binary orbital plane close to the companion.There is therefore less accretion on to the neutron star in a retrograde system than would otherwise be expected.Since X-ray outbursts are driven by accretion on to the neutron star, these outbursts are expected to be weaker in a retrograde system.

Figure 1 .
Figure1.Surface density (upper panels), inclination to the binary orbital plane (middle panels) and nodal phase angle (lower right panel) of the prograde (left panels) and retrograde (right panels) discs as a function of radius at different times in the simulations.The surface densities are shown where they are greater than 5 × 10 −4 g cm −2 .The nodal phase angle is not well defined for a coplanar disc, so we only plot the values for the portions of the disc with inclination  < 179 • for a retrograde disc.

Figure 2 .
Figure 2. Column density of the Be star disc from the SPH simulations at a time of  = 41.6  orb .The top panels show the prograde disc and the bottom panels show the retrograde disc.The Be star is represented by the large white circle while the small white dot surrounded by a red circle represents the companion neutron star.The binary separation is 95  ⊙ .The size of the white circles are scaled to the size of the sink radius.The left panels show the  −  plane in which both the Be star and the neutron star orbit.The right panels show the  −  plane.

Figure 3 .
Figure 3. Test particle inclination relative to the binary orbital plane (upper panel) and eccentricity (lower panel) as a function of time.The initial orbital radius is  = 70 (gray), 80 (black) and 82  ⊙ (blue).The retrograde disc properties (from Section 2.3) are shown at  = 70  ⊙ (orange) and 80  ⊙ (green).The blue line is truncated at the time when the particle is ejected (eccentricity > 1).

Table 1 .
Parameters for the binary and the circumprimary disc.