Be discs in binary systems: II. Misaligned orbits

TL;DR

This study uses SPH simulations to explore how misaligned binary companions influence Be star disc structures, revealing effects on truncation, tilt, and density, with implications for understanding disc dynamics in binary systems.

Contribution

It provides a systematic analysis of the impact of misalignment, viscosity, and orbital period on Be star disc structure and tilt, extending previous models with detailed simulation results.

Findings

Misalignment affects outer disc density and truncation radius.

Disc tilt increases with radius, peaking around 5 stellar radii.

Analytical scale height formula matches well except in density-enhanced regions.

Abstract

We use a smoothed particle hydrodynamics (SPH) code to examine the effects of misaligned binary companions on Be star discs. We systematically vary the degree of misalignment between the disc and the binary orbit, as well as the disc viscosity and orbital period to study their effects on the density in the inner and outer parts of the disc. We find that varying the degree of misalignment, the viscosity, and the orbital period affects both the truncation radius and the density structure of the outer disc, while the inner disc remains mostly unaffected. We also investigate the tilting of the disc in the innermost part of the disc and find the tilt increases with radius until reaching a maximum around 5 stellar radii. The direction of the line of nodes, with respect to the equator of the central star, is found to be offset compared to the orbital line of nodes, and to vary periodically in time, with a period of half a orbital phase. We also compare the scale height of our discs with the analytical scale height of an isothermal disc, which increases with radius as $r^{1.5}$. We find that this formula reproduces the scale height well for both aligned and misaligned systems but underestimates the scale height in regions of the disc where density enhancements develop.

Tables

Table 1. Table 1: Adopted parameters of the binary system.

Parameters	Values
Mass of primary	2.9 M☉
Radius of primary	3.67 R☉
Effective temperature	12000 K
Mass of secondary	0.31 M☉
Radius of secondary	0.38 R☉
Eccentricity	0
Mass injection rate	10-8 M☉/yr
Particles mass	$4\times {10}^{-15}{M}_{☉}$
Viscosity parameters	0.1 / 0.5 / 1.0
Misalignment angles	0$°$ / 30$°$ / 45$°$ / 60$°$
Period (days)	30 / 60s

Table 2. Table 2: Approximate number of active particles after 50 P o ​ r ​ b subscript 𝑃 𝑜 𝑟 𝑏 P_{orb} for short period (30 day) systems.

Viscosity	Misalignment	Number of active
parameter	angle	particles
	0$°$	100,000
${\alpha }_{\mathrm{SS}}=0.1$	30$°$	79,000
	45$°$	75,000
	60$°$	18,000
	0$°$	37,000
${\alpha }_{\mathrm{SS}}=0.5$	30$°$	36,000
	45$°$	37,000
	60$°$	39,000
	0$°$	25,000
${\alpha }_{\mathrm{SS}}=1.0$	30$°$	25,000
	45$°$	25,000
	60$°$	26,000

Table 3. Table 3: Steady-state values for the short period (30 day) systems.

Viscosity	Misalignment	Inner slope	Outer slope	Truncation radius
parameter	angle	$m$	$n$	$R_T$ (stellar radii)
	0$°$	1.49	7.78	6.91
${\alpha }_{\mathrm{SS}}=0.1$	30$°$	1.58	6.60	6.97
	45$°$	1.67	5.42	7.00
	60$°$	1.65	4.97	7.72
	0$°$	1.92	5.48	8.03
${\alpha }_{\mathrm{SS}}=0.5$	30$°$	1.99	3.95	8.41
	45$°$	1.99	3.43	9.27
	60$°$	2	3.81	10.9
	0$°$	1.93	3.31	7.68
${\alpha }_{\mathrm{SS}}=1.0$	30$°$	2.06	2.47	8.62
	45$°$	2.06	2.22	10.1
	60$°$	2.06	3.17	12.3

Table 4. Table 4: Steady-state values for the long period (60 day) systems.

Viscosity	Misalignment	Inner slope	Outer slope	Truncation radius
parameter	angle	$m$	$n$	$R_T$ (stellar radii)
	0$°$	1.61	5.23	10.34
${\alpha }_{\mathrm{SS}}=0.1$	30$°$	1.8	4.41	10.81
	45$°$	1.81	3.96	11.34
	60$°$	1.78	3.55	10.69
	0$°$	1.95	3.32	12.25
${\alpha }_{\mathrm{SS}}=0.5$	30$°$	1.95	2.75	12.32
	45$°$	1.96	2.54	12.89
	60$°$	2.01	0.22	13.37
	0$°$	2.08	2.27	13.30
${\alpha }_{\mathrm{SS}}=1.0$	30$°$	2.09	2.01	13.34
	45$°$	2.12	1.81	13.84
	60$°$	2.11	1.55	14.23

Equations

$$R_{\mathrm{L}}=\frac{0.49q^{2/3}}{0.6q^{2/3}+\ln (1+q^{1/3})}a,$$

$${⟨\Sigma (r)⟩}_{\varphi }\propto \frac{(r/R_t{)}^{-m}}{1+(r/R_t{)}^{n-m}},$$

$$z_{\mathrm{CoM}}(r,\varphi )=\frac{{\sum }_k{\rho }_k(r,\varphi )z_k}{{\sum }_k{\rho }_k(r,\varphi )},$$

$$\gamma (r,\varphi )=\arctan \frac{z_{\mathrm{CoM}}(r,\varphi )}{r}.$$

$$\gamma (\varphi {)}_r={\tau }_r\sin (\varphi +{\beta }_r),$$

$$H_{\mathrm{th}}\equiv H(r)=c_s/\Omega (r),$$

$$\rho (z)\propto e^{-z^2/2H_{\mathrm{sim}}^2},$$

$$\rho (z)\propto e^{-(z-\lambda {)}^2/2H_{\mathrm{sim}}^2}.$$

Full text

††pagerange: Be discs in binary systems: II. Misaligned orbits–Be discs in binary systems: II. Misaligned orbits

1 INTRODUCTION

The B-emission or Be stars are predominantly characterized by emission in the Balmer lines and have been studied since they were first detected by Secchi (1867). In addition to the emission in the hydrogen lines, emission is often observed in singly or doubly ionized metals and is a direct result of radiative recombination in a disc-like distribution of gas. Other defining features include continuum excess in the visible up to the radio, linear polarization and variability over a significant range in period. These features originate in a geometrically thin disc of material ejected from the rapidly rotating central star. During the time since Secchi’s discovery, our understanding of these star/disc systems has steadily improved but there are a significant number of remaining puzzles. The current status of Be-star research has been recently reviewed by Rivinius et al. (2013).

In a pioneering paper, Shakura & Sunyaev (1973) investigated the transfer of angular momentum in discs surrounding black holes and in doing so established a framework for future studies of hydrodynamic discs. The review paper by Pringle (1981) summarises the major properties of accretion discs including detailed discussion about the role of viscosity and emission processes. Since that time there have been a multitude of investigations that adopt and apply accretion theory to model astrophysical discs that surround a variety of different types of objects.

The same prescription originally developed for viscous accretion discs can be applied to study Be star discs except, in this case, the direction of the radial motion can switch between outward flowing, when the disc is being actively fed by the star, to inward flowing, as a pre-existing disc passively dissipates. In addition to the differences in the radial velocity, accretion discs fuel jets to shed angular momentum while decretion discs take the angular momentum that has been contributed from the central star and put material, at least some of it, into Keplerian orbit (Baade et al., 2016). The viscous decretion disc (VDD) scenario for Be stars, originally adopted by Lee et al. (1991) assumes that angular momentum is transferred from the star to the disc via an unknown mechanism. His work was followed up by Okazaki (2001) who found that the outflow in the inner disc is subsonic and rotating in a Keplerian fashion, consistent with current findings about Be stars. The viscosity parameter is given by the usual $\alpha_{\mathrm{SS}}$ viscosity parameter and for Be stars is usually assumed to be in the range of 0.1 to 1 (see Lee et al., 1991). The VDD has been the focus of an increasing number of studies. For example, Martin et al. (2011) studied warping and precession of Be star discs and find that in non-coplanar binary systems the disc can become distorted and warped. Fu et al. (2015) studied oscillations, in particular the Kozai-Lidov mechanism, in hydrodynamic discs in binary systems by considering the effect of disc pressure and viscosity. They find that given sufficient time the disc becomes co-planar with the binary companion. Lubow et al. (2015) find that in misaligned systems the disc can be more extended in radial distance than a disc in a co-planar system. Despite the fact the focus of Fu et al. (2015) and Lubow et al. (2015) studies are about accretion discs, their work can help provide insight to the results presented in this investigation.

For the case of an isolated Be star, the central star supplies matter and angular momentum to the disc and is the primary energy input to the disc through its radiation field. The stellar wind potentially could ablate the disc material or help to constrain material to the equatorial regions especially in the inner parts of the disc (Okazaki, 2012). However, in a binary system the situation is more complex. In addition to tidally truncating the disc, the companion interacts with the disc through resonant torques which could cause the disc to warp and/or precess (Okazaki, 2012; Kee et al., 2016) and these effects may complicate the density structure for the case of both aligned and misaligned discs.

The innermost disc is dominated by viscous torques so that isolated and binary systems may be similar in the inner regions (Okazaki, 2012). However, as explained by Okazaki et al. (2002), with increasing radial distance the resonant torques due to the binary companion begin to play a bigger role. Okazaki et al. (2002) studied the gravitational effects of close companions on the structure of the discs in Be/X-ray binary systems and found that the disc was truncated, at a certain point, called the truncation radius, by the companion. Similar disc truncation have also been predicted in accretion systems with eccentric and misaligned binary companions by Artymowicz & Lubow (1994) and Larwood et al. (1996), respectively. This was also suggested by Reig et al. (1997) who studied Be/X-ray systems and found a relationship between the size of the H$\alpha$ emitting region and orbital period, demonstrating that the companion, in their case a neutron star, halted the flow of material and effectively truncated the disc. In a study of Be discs in coplanar binary systems, Panoglou et al. (2016) (hereafters Paper I) found that the position of the truncation radius depends on the parameters of both the disc (kinematic viscosity) and the binary (orbital period, mass ratio, eccentricity).

Just inside the truncation radius, it was found that the disc density decreases at a slightly slower rate with increasing distance from the star than what is observed in discs around isolated Be stars (Okazaki et al., 2002, Paper I). As a result, the radial density distribution is flatter in binary systems and consequently discs in binary systems may be denser, and therefore more massive, than those in isolated systems. This phenomenon was called the “accumulation effect”. In Paper I, the effect was found to be stronger the smaller the viscosity and the orbital period and larger the mass ratio. Outside the truncation radius (i.e. the outermost portions of the disc), however, the density falls off at a substantially increasing rate. Okazaki et al. (2002) as well as Paper I also observed that $\alpha_{\mathrm{SS}}$ had an effect on the rate of the density fall-off in the outer part of the disc, with high viscosity discs having slower drop-off rates than low viscosity discs (see figure 11 of Paper I).

In addition to these changes in density, Okazaki et al. (2002) found that the surface density of the disc revealed the development of a two-armed spiral density wave at periastron (figures 10 and 11 of Okazaki et al., 2002) with one arm preceding the companion and the other located on the opposite side of the disc. Similar spiral waves were also observed in Paper I.

Some progress has been made in recent years regarding the value of $\alpha_{\mathrm{SS}}$ in the disc. Carciofi et al. (2012) modelled the 2003 dissipation phase of 28 CMa and found that $\alpha_{\mathrm{SS}}=1.0\pm 0.2$. Subsequently, Ghoreyshi & Carciofi (2017) revisited this value to about 0.3, and extended the modelling to 4 complete cycles of disc formation and dissipations, spanning almost 4 decades of disc activity. They found that the values of $\alpha_{\mathrm{SS}}$ are systematically larger during disc build-up than at dissipation, and typically range between 0.1 and 1. In a recent work, Rímulo et al. (2017) modelled the lightcurves 54 Be stars in the SMC and found that the viscosity parameter is roughly two times larger at build-up ($\left\langle\alpha_{\rm SS}^{\rm bu}\right\rangle=0.63$) than at dissipation ($\left\langle\alpha_{\rm SS}^{\rm d}\right\rangle=0.29$)

In this study we use a numerical hydrodynamics code to expand on the works of Okazaki et al. (2002) and in Paper I by investigating how these system evolve when the plane of the disc and binary are misaligned. Our goal is to understand how this misalignment, as well as various disc and orbital parameters, affect the density structure and dynamics of the disc. While much of the recent work in the literature is focused on X-ray and $\gamma$-ray Be binary systems, in this study we investigate, in a systematic fashion, the effect of a low mass, main sequence binary companion on Be star discs for misaligned orbits over a range of $\alpha_{\mathrm{SS}}$ and orbital period. The organization is as follows: Section 2 describes our methodology, results are presented in Section 3, and a discussion and summary is provided in Section 4.

2 Methodology

2.1 Disc Modelling

In order to study the hydrodynamical effects of a misaligned companion on the structure and dynamics of the disc, a three-dimensional (3D) smoothed particle hydrodynamics (SPH) code was used. This type of code solves the fluid equations by dividing the fluid into discrete elements called particles. The properties of each particle are smoothed over a finite spacial distance (smoothing length) using a kernel function in order to simulate a continuous fluid. The SPH code used in this work is based on the code developed by Benz (1990) and Benz et al. (1990). It was later modified by Bate (1995) to include a second-order Runge-Kutta-Fehlberg integrator which uses individual time steps for each of the particles in order to decrease computational time. This code was later refined by Okazaki et al. (2002) for use in simulations of Be-binary systems.

The disc systems constructed and followed in this work are binary systems where the primary is a Be star and the secondary is a low mass non-degenerate star assumed to be in a circular orbit. Each simulation starts as a disc-less system at the beginning of its disc building phase. The injection of mass into the disc is accomplished by placing a prescribed number of particles around the star at a radius of $r_{\mathrm{inj}}=1.04R_{\star}$, referred to as the injection radius. The position of these particles along $r_{\mathrm{inj}}$ is chosen randomly while their velocity is set to be Keplerian (no radial momentum). This is done for each time step of the simulation. All particles were given the same mass, which remained constant throughout the simulations. These masses, $m_{\mathrm{part}}$, are calculated based on the mass injection rate ($\dot{M}_{\mathrm{inj}}$) of the system and the number of particles ($N_{\mathrm{part}}$) created during each mass injection event, both of which are constant throughout the simulations. For this work these values were set to $\dot{M}_{\mathrm{inj}}=10^{-8}$ M*☉/yr and $N_{\mathrm{part}}=6000$, resulting in $m_{\mathrm{part}}\approx 4\times 10^{-15}$ M☉*.

Once injected into the disc, particles interact with each other through viscosity, which transfers angular momentum throughout the disc. Most particles eventually fall back onto the star while the remainder move outward, carrying angular momentum with them away from the central star. As shown by Rímulo et al. (2017), the fraction that falls back is about $\dot{M}_{\rm inj}[2.-(R_{\mathrm{inj}}/R_{\mathrm{star}})^{1/2}]$, or about 98$\%$ in our case. The Shakura-Sunyaev viscosity $\alpha_{\mathrm{SS}}$-prescription (Shakura & Sunyaev, 1973) is used to define the viscosity of the system. Details about the implementation of the viscosity into the SPH code can be found in Okazaki et al. (2002). The disc is assumed to be isothermal, therefore the gas particles are all set to the same temperature. The disc temperature ($T_{\mathrm{disc}}$) was set to 60% of the effective temperature of the star, which was found to be a good approximation for the temperature structure of the disc (Carciofi & Bjorkman, 2006; Millar & Marlborough, 1998).

The stars are modelled using sink particles delimited by an accretion radius inside which any particle is assumed to have accreted onto the star and is therefore removed from the simulation. The accretion radius of the primary (Be star) is equal to the radius of the star, which in this work was set to 3.67 $R_{\star}$, typical of a main-sequence B5 star. The accretion radius for the secondary is determined using the Eggleton approximation of the Roche lobe $R_{\mathrm{L}}$.

[TABLE]

where q is the mass ratio of the stars ($q=M_{\mathrm{sec}}$/$M_{\mathrm{prim}}<1$) and a is the separation distance between them. Note that a changes over time for eccentric systems, as is discussed in Paper I, however since we focus on circular orbits in this work $R_{\mathrm{L}}$ remains fixed throughout each simulation. The masses, radii, and effective temperatures of the stars were based on the parameters of the Pleione binary system (Harmanec, 1988; Hirata, 2007; Nemravová et al., 2010). It’s important to note that this study does not focus on the modelling of this particular binary system, but instead aims to study the general behaviour of Be discs in close binary systems. For this purpose, we have modified a few parameters, such as the orbital period and eccentricity of the system in order to reduce the complexity of the simulations. All these parameters can be found in Table 1.

In order to investigate the effect of the companion on the disc, we varied three parameters: the misalignment angle $\theta$, defined as the angle between the stellar equator of the primary and the plane of the binary orbit, the viscosity parameter, $\alpha_{\mathrm{SS}}$, and the orbital period of the binary. Simulations were run for four misalignment angles ($\theta$ = 0$\degr$, 30$\degr$, 45$\degr$, 60$\degr$), three values of viscosity ($\alpha_{\mathrm{SS}}$ = 0.1, 0.5, 1.0), and for binary systems with short (30 days) and long (60 days) orbital periods, corresponding to a semi-major axis of 59.9 $R_{\sun}$ and 95.1 $R_{\sun}$, respectively.

Following Paper I, we define $p$ as the orbital phase of the secondary during a single orbit, with $0\leq p<1$. We examine four orbital phases, $p$ = 0.00, 0.25, 0.50, and 0.75. Figure 1 shows a graphical depiction of this notation scheme. Here the $x-y$ plane is defined as the plane of the stellar equator with the $z$ axis pointing out of the figure. The Be star is represented by the gray circle at the center of the figure while the solid black line shows the path of the secondary111Note that this is the orbital path projected onto the plane of the disc, explaining why the orbit appears elliptical.. For each position of the secondary (black circles) we have indicated the corresponding orbital phase $p$ and whether the secondary is above ($z>0$), below ($z<0$), or in ($z=0$) the plane of the disc. Note that in the aligned case ($\theta=0$) the path of the secondary is circular and the secondary stays in the plane of the disc at all times. A cylindrical coordinate system ($r$,$\phi$,$z$) is used when describing specific locations or cross-sections inside the disc, with the origin defined at the primary and the ($r$,$\phi$,0)-plane defined as the plane of the stellar equator. The angle $\phi$ is defined from the positive $x$-axis moving counter-clockwise. The dashed lines in Figure 1 show the directions of the azimuthal angles $\phi$ = 0$\degr$, 90$\degr$, 180$\degr$, and 270$\degr$.

3 Results

3.1 Average Surface Density

It is informative to calculate some of the average features of the disc prior to investigating localized changes or perturbations. We start with an investigation of the average disc density. This is accomplished by transferring the particles to a cylindrical grid centred at the primary. The density is then integrated vertically to obtain the azimuthally averaged surface density ($\left<\Sigma(r)\right>_{\phi}$). This procedure was adopted for Paper I.

One important aspect we first consider is whether our simulations are sufficiently developed so that the only major variations are those caused by the motion of the secondary. As Haubois et al. (2012) showed, discs around isolated Be stars will continue to grow as long as there is a continuous and constant injection of mass and angular momentum from the star, meaning they will never reach a state of perfect equilibrium. However we also know that the growth of the disc, although fast in the beginning of its building phase, slows dramatically as the disc expands outward (Haubois et al., 2012). The discs should therefore reach a quasi-static state (QSS), where the disc growth become insignificantly small within the span of one orbital period. It is also important to note that, since Be star discs are built from the inside out, the inner parts reaches QSS much faster than the outer parts. See section 2.1 of Paper I for a more in depth discussion.

Figures 2(a) through (f) show the temporal evolution of the azimuthally averaged surface density, $\left<\Sigma(r)\right>_{\phi}$, at various times during the disc building phase (see legend). The leftmost and rightmost panels show the results for short and long period systems, respectively, with each row, from top to bottom, for $\alpha_{\mathrm{SS}}$ of 0.1, 0.5, and 1.0, respectively. Each simulation was run up to 50 orbital periods ($P_{orb}$). Note that only the aligned systems ($\theta=0$) are shown in Figure 2, however we have found the QSS time is virtually the same for the corresponding misaligned cases.

We see that all simulations have reached QSS by 50 $P_{orb}$, although the higher viscosity discs reach their QSS much faster than lower viscosity discs. This result is similar to the temporal evolution of the surface density profile found in Paper I for aligned systems and by Haubois et al. (2012) for isolated Be stars; both found that discs with higher viscosity reach their steady-state faster than those with low viscosity. This happens because the $\alpha_{\mathrm{SS}}$ acts as a timescale factor in the hydrodynamical equations governing theses discs (Pringle, 1981). The number of active particles for the short period systems after 50 $P_{orb}$ is presented in Table 2.

Figures 3(a) through (f) show $\left<\Sigma(r)\right>_{\phi}$ after 50 full orbital periods. We see that the $\left<\Sigma(r)\right>_{\phi}$ profiles in each case are similar in shape as those found in Okazaki et al. (2002) in his study of Be/X-ray binaries, with shallower slopes in the inner part of the disc and steeper slopes in the outer parts. As Okazaki et al. (2002) explains this drop in surface density, or truncation, is due to the tidal forces generated by the secondary as it orbits the primary. Furthermore, similar to the findings of Okazaki et al. (2002), we notice that surface density profiles of the outer part of the disc become shallower with increasing $\alpha_{\mathrm{SS}}$. This can be attributed to the fact that high viscosity discs have a faster recovery time than low viscosity discs, since viscosity is the primary means with which angular momentum in the disc is transported away from the star.

Now we move onto a systematic investigation of the changes in the disc over the range in parameters described previously. The misalignment angle also affects the surface density profile of the disc as seen in Figure 3. In all six panels we see that the greater the misalignment angle, the denser the outermost part of the disc. Conversely, the density in the inner parts appears to be mostly unaffected by the misaligned companion. In order to better investigate this, three properties of these profiles were examined, the truncation radius and the fall-off rate of the surface density in the inner and outer parts of the disc. The inner and outer parts are defined with respect to the truncation radius, the inner disc being the region inside the truncation radius and the outer part exterior to it. For this analysis the azimuthally averaged surface density profiles were fitted using the same function adopted by Okazaki et al. (2002).

[TABLE]

where $R_{t}$ is the truncation radius, m is the power-law exponent of the inner disc, and n is the power-law exponent of the outer disc. Hereafter we will refer to m and n as the inner and outer exponent, respectively. These three parameters are obtained from this fit.

Figure 4 shows the temporal evolution of the surface density fall-off rate, $m$, in the inner parts of the disc. We see that the fall-off rate is initially high in most cases, but rapidly decreases until a steady-state value is reached. The estimated steady-state values for $m$, obtained by an exponential fit, can be found in Tables 3 and 4 for the short and long period systems, respectively.

However, where Haubois et al. (2012) found that the power law for isolated Be-stars tended toward $m=2$, we found, in almost all cases, that the power law tended toward values of $m\leq 2$, which indicates the inner part is denser than a typical disc in an isolated Be star system. The same accumulation effect was also observed in Paper I for aligned systems. As we see from Figure 4 as well as Tables 3 and 4, the accumulation effect is stronger, overall, in aligned systems compared to misaligned systems, with little variation in the latter. We notice that the difference in accumulation between the aligned and misaligned systems increases substantially with viscosity. We also see that both the viscosity coefficient and orbital period have effects on the steady-state value of $m$. These dependences are in agreement with the results of Paper I. The only system where this effect is absent is in the system with both high viscosity and long orbital period (panel f). Here we see that $m$ tends toward values equal or greater than 2. This indicates that the secondary is distant and viscosity is large enough that the torque is too weak to affect the density distribution in the inner disc, which then behaves like an isolated Be disc.

Figure 5 shows the temporal evolution of the surface density fall-off rate, $n$, in the outer parts of the disc. We see here that the slope $n$ reaches a steady state value more quickly than $m$. However it is important to note that although the outer density slope reaches QSS faster, the density scale still requires a longer time to reach QSS compared to the inner disc. The estimated steady-state values for $n$ are also presented in Tables 3 and 4. As seen in Figures 3 and 5 the fall-off rate in the outer disc is dependent on the misalignment angle, i.e. slower fall-off rates for highly misaligned systems.

The results for the truncation radius, obtained from Equation 2, are shown in Figure 6 with the steady-state values in Tables 3 and 4. The panel layout is the same as Figures 3 and 5. We see that the steady-state values of the truncation radius tends to be closer to the star for low misalignment angles and move farther out as the misalignment angle increases. We also notice that the truncation radius is farther away from the star in systems with higher disc viscosity and longer orbital periods. This is similar to the findings presented in Paper I for aligned systems.

The relationship between $R_{t}$, $m$, and $n$ with regards to viscosity, orbital period, and misalignment angle is relatively intuitive. As mentioned earlier, viscosity determines how fast the disc can grow and recover between each interaction with the secondary. We therefore expect high viscosity discs to have more time to relax back to their unperturbed state than less viscous ones, reducing the impact of the accumulation effect (higher $m$) and allowing more material to flow into the outer disc (lower $n$ and higher $R_{t}$). A longer orbital period also means that the outer disc has more time the rebuild itself between passages, and since the secondary is radially farther away from the disc, it has a smaller impact on the truncation of the disc. Similarly, increasing the misalignment angle also increases the vertical distance between the disc and the secondary when the latter is at its maximum elevation. This increased separation weakens the tidal torque on the disc resulting in a weaker truncation.

3.2 Disc Warping

In order to study the warping of the disc, we consider the vertical displacement of the disc above and below the equatorial plane along various radial directions. To accomplish this, we first determine the position of the centre of mass along the vertical axis ($z_{\mathrm{CoM}}$) along the $r$-axis for various values of $\phi$ (dashed lines in Figure 1), using the following equation:

[TABLE]

where $\rho_{k}(r,\phi)$ and $z_{k}$ are the density and corresponding vertical position of a column density located at $(r,\phi)$. Next we convert the vertical position $z_{\mathrm{CoM}}$ into an angular position $\gamma$ by trigonometry:

[TABLE]

We then calculated the tilt of each ring section of the disc by fitting the following equation using the results obtained by Equation 4 for multiple values of $r$:

[TABLE]

where $\tau_{r}$ the tilting angle of the disc at radius $r$, defined as the angle formed by the normal of the disc section and the $z$-axis. Similarly, $\beta$ represents the direction of the line of node of the ring section in relation to the equator of the primary.

An example of this procedure can be found in Figure 7, which shows the relation between the angular position, $\gamma$, of the center of mass, defined from the equatorial plane, and the azimuthal angle, $\phi$, at a radial distance of $r=2.79$ stellar radii. Data taken from the simulation are represented by the black dots while the solid gray line shows the result of the best fit using Equation 5. Note that, close to the star, the scatter on the data from the simulation was higher due to the constant injection of new mass, however we managed to obtain a good fit.

Figure 8 shows the tilting angle, $\tau$, as a function of radius at various orbital phases ($p$ = 0.0, 0.25, 0.5, 0.75) for our three misaligned systems.

As expected, we see clear evidence of vertical displacement in the disc as well as other interesting features. In all three cases, the tilting increases steadily in the first few stellar radii until reaching a peak around 4-5 stellar radii. Furthermore this increase does not appear to vary with orbital phase suggesting that the tilting in the very inner region is quite stable. We also notice that the tilt becomes more significant as the misalignment angle increases from 30$\degr$ to 45$\degr$, but then decreases slightly when increased to 60$\degr$.

Similarly to Figure 8, Figure 9 shows the direction of the line of nodes, $\beta$, as a function of radius at various orbital phases ($p$ = 0.0, 0.25, 0.5, 0.75) for our three misaligned systems. The line of node is defined with respect to the plane of the equator of the central star and its direction is defined in the same way as the azimuthal angle $\phi$, as depicted in Figure 1. Also note that due to we will only use one angular value to define lines of nodes to avoid confusion. We notice that the disc sections tend to tilt toward $\phi=$ 80$\degr$ to 100$\degr$. This is an offset of about 90$\degr$ compared to the line of nodes of the orbital plane, which is $\phi=0\degr$ for all misaligned systems in this study. Furthermore, we see that this offset is greater for lower misalignment angles and decreases as the misalignment angle increases. This offset also increases as we move farther from the star.

Looking at each panel we also notice small changes in the line of nodes over one orbital period. In all three cases the direction of the tilt seems to oscillate with a period of half a phase. The amplitude of this oscillation also appears to increase with misalignment angle.

Finally, we note that both the tilting angle and line of nodes appears to fluctuate at distances greater than 5 stellar radii. This is likely do to the fact that Equation 5 assumes that the tilted sections are perfectly circular, which is a good approximation close to star but, in the outer discs, the two-armed spiral arms develop as a result of the individual orbits of the particles becoming elliptical. Therefore, a circular fit fails to track the motions of the ensemble of particles around the star.

3.3 Scale Height

The scale height (or effective thickness) of the disc is an important quantity sometimes assumed or estimated for modelling discs. Variations in scale heights can have a significant effect on the density structure and the temperature distribution of the disc, which, in turn, will have an effect on the estimation of observables such as emission line profiles and IR continuum levels (Carciofi & Bjorkman, 2008). Therefore we investigate the disc scale height in this work. For an isolated Be star, the scale height is defined as

[TABLE]

where $\Omega$ is the Keplerian angular velocity at a distance $r$ from the star (Lightman, 1974). This relationship is obtained by assuming that the disc is in vertical hydrostatic equilibrium. A derivation of this equation can be found in Carciofi & Bjorkman (2008). Equation 6 assumes the presence of a single source of gravitational potential, the Be star. It is therefore important to determine whether or not Equation 6 is appropriate to estimate the scale height in Be binary systems that are aligned and misaligned.

In order to extract the scale height throughout the disc, we computed the vertical density structure of the gas at various points in the disc. The vertical density structure of the disc can be described using a Gaussian function

[TABLE]

where $H_{\mathrm{sim}}$ is the scale height of the disc. Equation 7, however, assumes a Gaussian centred at $z=0$, which is a good assumption for isolated Be star systems or aligned systems but not for misaligned systems where we have a warped disc. We therefore modified Equation 7 to include a second fitting parameter, $\lambda$, which represents how shifted the center of mass is at this particular location compared to $z=0$:

[TABLE]

We then use this relation to fit the vertical density structure of the disc at various locations in order to calculate the scale height for our simulated discs.

Figure 10 shows the ratio between the scale height obtained through simulations ($H_{\mathrm{sim}}$) and the theoretical scale height ($H_{\mathrm{th}}$) computed from Equation 6, for short period systems, $\alpha_{\mathrm{SS}}=0.5$, and with misalignment angles of 30$\degr$.

Overall we see that Equation 6 is relatively good at estimating the scale height. However, deviations between simulations and theory can be seen in regions where density waves typically appear in Be-binary systems (see figure 4 of Paper I and figures 10 and 11 of Okazaki et al., 2002). This indicates that Equation 6 might not adequately describe the scale height in the regions of the disc where these density enhancements appear. Note that these results are consistent for all misalignment angles.

4 Discussion and Conclusion

In this investigation we use numerical simulations to study the behaviour of Be star discs in aligned and misaligned binary systems. Simulations were obtained using a smoothed-particle hydrodynamics (SPH) code following previous methods developed and described in the literature (Benz, 1990; Benz et al., 1990; Okazaki et al., 2002). In our simulations, we varied the degree of misalignment between the plane of the orbit and the plane of the disc, the viscosity parameter, $\alpha_{\mathrm{SS}}$, and the orbital period of the system, $P_{orb}$.

First we studied the azimuthally averaged surface density, $\left<\Sigma(r)\right>_{\phi}$, by investigating the density fall-off rate in the inner and outer parts of the disc ($m$ and $n$ respectively), and the truncation radius ($R_{t}$). We find, as reported in Paper I for coplanar orbits, that material tends to build up in the innermost part of the disc, a phenomenon called the accumulation effect, for our misaligned systems with the density fall-off values of $m$ reaching values below 2 for most systems. Furthermore we find that the degree of misalignment does have an effect on the magnitude of the accumulation, with a greater accumulation (lower $m$) in aligned systems and with less accumulation in misaligned systems. This difference between aligned and misaligned systems is found to also increase with viscosity. Only for discs with the highest viscosity ($\alpha_{\mathrm{SS}}=1.0$) and the longest orbital period ($P_{orb}$ = 60 days) do we find little evidence of accumulation. As reported in Paper I, both $\alpha_{\mathrm{SS}}$ and $P_{orb}$ also have an influence on the accumulation found in the inner disc. Similarly we find that $R_{t}$ also has a dependence on the misalignment angle: discs in aligned systems are found to be truncated closer to the star compared to discs in misaligned systems, with the truncation farther from the primary as the misalignment increases. This is in agreement with the findings Lubow et al. (2015) that discs in misaligned systems have a greater radial extent than their counterparts in aligned systems. We also find that $\alpha_{\mathrm{SS}}$ and $P_{orb}$ have similar effects on the truncation, that is, the truncation is farther from the primary as both $\alpha_{\mathrm{SS}}$ and $P_{orb}$ increase. This is also in agreement with the results of Paper I. Finally, $n$ is found to vary with the misalignment of the system as well. Larger misalignment angles result in denser outer discs, i.e. smaller $n$.

We also studied the warping of the disc by computing the tilting angle and line of nodes of ring segments in the inner part of the disc. We find that the tilt of the disc increases steadily as we move farther from the star, reaching a peak around 4-5 stellar radii. The maximum tilt is found to be greater at a misalignment angle of 45$\degr$ compared to 30$\degr$ and 60$\degr$. The direction of the line of noes is found to vary with both distance from the primary star and with orbital phase. We also find that an offset of about 90$\degr$ between the line of nodes of the disc and of the orbit. This offset is found to decrease when the misalignment angle of the system increases.

Our findings about disc tilting are interesting as they may help to understand the variations in the polarimetric position angle observed in many Be stars. Linear polarization in Be stars is due primarily to the scattering of starlight by free electrons in the innermost region of the disc (Coyne & Kruszewski, 1969; Zellner & Serkowski, 1972) and is used to measure the position angle of Be stars in the sky. Changes in position angle are due to either variations in mass injection from central star or geometric changes in innermost part of disc. Our results suggests that these variation can also be explained by the presence of a misaligned binary companion. For example, Vince et al. (1995) found evidence of small cyclical variations, 30$\degr$ or smaller, in the polarimetric position angle of multiple Be stars, which is in agreement with our results.

Finally we compared the scale heights of our discs with the theoretical scale heights calculated using the standard equation for the case of an isolated Be star (Carciofi & Bjorkman, 2008). We find that overall this approximation is relatively good at describing the scale height for all misalignment angles. However, this standard prescription is too simplistic for describing the scale height in regions where density enhancements are present. This should be given careful consideration in future studies.

Martin et al. (2014) studied Kozai-Lidov oscillations in hydrodynamic discs and showed that for highly misaligned systems this mechanism could be set up. Figure 1 in Fu et al. (2015) shows that conspicuous disc oscillations can occur for orbital warpings that are large (i.e. 60 and 70*∘). During these oscillations the disc can become quite eccentric over time as disc tilting and eccentricity are interchanged. Martin et al. (2014) suggest that these oscillations may be important for transient discs in Be/Xray systems. Using equation 3 in Fu et al. (2015) we find that Kozai oscillations would not be expected for our simulations computed to 50 Porb*. According to our calculations the simulations would need to be followed for approximately 80 Porb in order to see evidence of Kozai oscillations in the inner disc. Nevertheless, in a future study we plan to study binaries in eccentric orbits and we plan to closely examine the disc shape and variations in inclination to search for the presence of Kozai oscillations.

Given we find that the disc density distribution and the geometry of the disc is altered in misaligned systems, we plan to adapt our codes so that we can compute observables directly from our simulations. Changes in disc density in the inner disc would be seen in observables such as polarization, UBV photometric colours, and spectroscopy, while interferometry, infrared and radio measurements could be used to predict the overall size and density in the disc. Changes in shape and strength of spectral lines could occur, depending on where they are formed in the disc, especially if portions of the disc or central star become obscured along the observers line of sight due to warping or tilting of the disc. As well, in future, we plan to extend this work to investigate the effects on the disc for non-circular binary systems.

Acknowledgements

We thank the anonymous referee for comments and suggestions which helped us improve the quality of this paper. This work was made possible by the facilities of the Shared Hierarchical Academic Research Computing Network (SHARCNET:www.sharcnet.ca), Compute/Calcul Canada, and Laboratory of Astroinformatics (IAG/USP, NAT/Unicsul, FAPESP grant No 2009/54006-4). CEJ acknowledges support from NSERC, the Natural Sciences and Engineering Research Council of Canada. DP acknowledges support from FAPESP grant No 2013/16801-2 and CNPq (300235/2017-8). ACC and CEJ thanks UofT-FAPESP-UWO joint researcher program. ACC acknowledges the support from CNPq (grant 307594/2015- 7) and FAPESP (grant 2015/17967-7).