Mass Transfer and Disc Formation in AGB Binary Systems

TL;DR

This paper uses numerical simulations to study how mass transfer and disc formation occur in binary systems with AGB stars, revealing complex outflow morphologies and accretion behaviors influenced by system parameters.

Contribution

It introduces a comprehensive simulation approach combining theory and adaptive-mesh refinement to analyze mass transfer mechanisms in AGB binary systems with varied configurations.

Findings

Accretion rates vary non-linearly with binary separation.

Different outflow morphologies observed, including wind Roche lobe overflow and Bondi-Hoyle accretion.

Potential implications for planetary nebulae shaping.

Abstract

We investigate mass transfer and the formation of disc in binary systems using a combination of numerical simulations and theory. We consider six models distinguished by binary separation, secondary mass and outflow mechanisms. Each system consists of an asymptotic-giant-branch (AGB) star and an accreting secondary. The AGB star loses its mass via a wind. In one of our six models, the AGB star incurs a short period of outburst. In all cases, the secondary accretes part of the ejected mass and also influences the mass-loss rate of the AGB star. The ejected mass may remain gravitationally bound to the binary system and form a circumbinary disk, or contribute to an accretion disk around the secondary. In other cases, the ejecta will escape the binary system. The accretion rate on to the secondary changes non-linearly with binary separation. In our closest binary simulations, our models exemplify the wind Roche lobe overflow while in our wide binary cases, the mass transfer exhibits Bondi-Hoyle accretion. The morphologies of the outflows in the binary systems are varied. The variety may provide clues to how the late AGB phase influences planetary nebulae shaping. We employ the adaptive-mesh-refinement code ASTROBEAR for our simulations and include ray-tracing, radiation transfer, cooling and dust formation. To attain the highest computational efficiency and the most stable results, all simulations are run in the corotating frame.

Tables

Table 1. Table 1: Planck mean opacity of olivine and pyroxene at different temperature of blackbody radiation. The energy distribution of the radiation is assumed to be in Planck distribution of that temperature. The dust grains are spherical and have MRN size distribution.

$T_{\text{bb}}$(K)	${\kappa }_{\text{olivine}}$( cm2 g-1)	${\kappa }_{\text{pyroxene}}$( cm2 g-1)
2800	1787	461
2600	1540	356
2400	1303	269
2200	1076	197
2000	858	140

Table 2. Table 2: AGB star’s physical properties.

Mass	$1{\text{M}}_{☉}$
Luminosity	$2342{\text{L}}_{☉}$
Pulsating boundary density	$1.2\times {10}^{-8}$gcm-3
Outflow temperature	$2900$K

Table 3. Table 3: The first column lists the model number used throughout this paper; d ​ ( au ) 𝑑 au d(\text{au}) is the binary separation; d x , d y ​ ( au ) subscript 𝑑 𝑥 subscript 𝑑 𝑦 au d_{x},d_{y}(\text{au}) and d z ​ ( au ) subscript 𝑑 𝑧 au d_{z}(\text{au}) are the physical dimensions of x,y and z of the simulation box and we take d x , d y ​ ( au ) subscript 𝑑 𝑥 subscript 𝑑 𝑦 au d_{x},d_{y}(\text{au}) for all simulations; N x , N y subscript 𝑁 𝑥 subscript 𝑁 𝑦 N_{x},N_{y} and N z subscript 𝑁 𝑧 N_{z} are the number of base cells in x,y and z dimensions. The AMR column lists the highest level of AMR for each model. Burst refers to the bursting state of the AGB star in section 2.2.2 . Only the AGB star in model 6 has a burst. r acc subscript 𝑟 acc r_{\text{acc}} is the accretion radius of the sink particle. r flux subscript 𝑟 flux r_{\text{flux}} is the radius of the shell through which the escaping flux is measured.

model	$M_{\text{AGB}}$	$M_{\text{sec}}$	$d$	$d_x,d_y$	$d_z$	$N_x,N_y$	$N_z$	AMR	Burst	$r_{\text{acc}}$	$r_{\text{flux}}$
	M☉	M☉	au	au	au					au	au
1	1.0	0.1	3	40	40	100	100	4	No	0.1	3
2	1.0	0.5	4	96	48	60	30	6	No	0.1	4
3	1.0	0.5	6	80	40	100	50	5	No	0.1	7
4	1.0	0.5	8	128	48	80	30	6	No	0.1	7
5	1.0	0.5	10	192	64	60	20	7	No	0.1	9
6	1.0	0.5	4	96	48	60	30	5	Yes	0.2	N.A.

Table 4. Table 4: Model numbers and the average mass loss rate of the AGB star and the mass accretion rate of the secondary. We denote the percentage of the accretion rate to the mass loss rate as ζ 𝜁 \zeta . We also show the percentage ϵ italic-ϵ \epsilon of the changed mass loss rate to the mass loss rate in isolated AGB model.

model#	${\dot{M}}_{\text{AGB}}$	${\dot{M}}_{\text{acc}}$	$\zeta$	$ϵ$
	M☉ yr-1	M☉ yr-1	per cent	per cent
1	$2.960\times {10}^{-7}$	$1.159\times {10}^{-7}$	$39.2$	$28.1$
2	$3.325\times {10}^{-7}$	$1.224\times {10}^{-7}$	$36.8$	$43.9$
3	$2.208\times {10}^{-7}$	$4.513\times {10}^{-8}$	$20.4$	$-4.42$
4	$2.658\times {10}^{-7}$	$8.097\times {10}^{-9}$	$3.05$	$15.1$
5	$2.606\times {10}^{-7}$	$5.848\times {10}^{-9}$	$2.24$	$12.8$

Equations

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \mathbf{u}\right)$$

$$\frac{\partial \rho \mathbf{u}}{\partial t}+\nabla \cdot \left(\rho \mathbf{uu}\right)$$

$$\frac{\partial E}{\partial t}+\nabla \cdot \left[\left(E+p\right)\mathbf{u}\right]$$

$$p$$

$$n$$

$${\mathbf{f}}_{\text{g}}$$

$${\mathbf{f}}_{\text{rad}}$$

$${\mathbf{f}}_{\text{ni}}$$

$$\mathbf{\Omega }$$

$$E$$

$$e$$

$${\hat{\mathbf{r}}}_{\text{AGB}}$$

$${\hat{\mathbf{r}}}_{\text{sec}}$$

$$L_{\text{acc}}=214\frac{M_{\text{WD}}}{{\text{M}}_{\text{sun}}}\frac{10^9\text{cm}}{r_{\text{WD}}}\frac{{\dot{M}}_{\text{acc}}}{10^{-7}({\text{M}}_{\text{sun}}{\text{yr}}^{-1})}{\text{L}}_{\text{sun}}$$

$$\tau (r_0,\theta )={\int }_{r_{\text{ps}}}^{r_0}\kappa (r,\theta )\rho (r,\theta )\text{d}r$$

$$\rho \left(r,\theta \right)=\frac{{\int }_0^{2\text{uppi}}\rho \left(r,\theta ,\varphi \right)\text{d}\varphi }{2\text{uppi}}$$

$$\begin{gathered} {\kappa }_{\text{dust}}\left(T_{\text{bb}}\right)= \\ \frac{{\int }_0^{\infty }f\left(\lambda ,T_{\text{bb}}\right)\left({\sigma }_{\text{abs}}\left(\lambda \right)+{\int }_0^{\text{uppi}}{\sigma }_{\text{sca}}(\lambda ,{\theta }_{\text{s}})(1-\cos {\theta }_{\text{s}})\text{d}{\theta }_{\text{s}}\right)\text{d}\lambda }{m_{\text{p}}{\int }_0^{\infty }f\left(\lambda ,T_{\text{bb}}\right)\text{d}\lambda } \end{gathered}$$

$${\kappa }_{\text{total}}=\frac{{\rho }_{\text{gas}}{\kappa }_{\text{gas}}+{\rho }_{\text{dust}}{\kappa }_{\text{dust}}}{\rho },$$

$$T_{\text{eq}}=\sqrt[4]{\frac{\beta L\exp \left(-\tau \right)}{4\text{uppi}{\sigma }_{\text{sb}}{r}^2}}$$

$${\mathbf{f}}_{\text{gAGB}}=\alpha \frac{G{M}_{\text{AGB}}}{r_{\text{AGB}}^2}{\hat{\mathbf{r}}}_{\text{AGB}},$$

$$\frac{v_{\text{burst}}}{\text{km }{\text{s}}^{-1}}=5+35\exp \left(-\frac{((t/\text{yr}-32.5)/2{)}^{20}}{(5/2.2{)}^{20}}\right),$$

$$u_{\text{p}}=\mathbf{u}\cdot {\hat{\mathbf{n}}}_{\text{b}},$$

$${\dot{M}}_{\text{AGB}}={\dot{M}}_{\text{esc}}+{\dot{M}}_{\text{cd}}+{\dot{M}}_{\text{ad}}+{\dot{M}}_{\text{acc}}$$

$${\dot{M}}_{\text{AGB}}\approx {\dot{M}}_{\text{esc}}+{\dot{M}}_{\text{acc}}$$

$${\rho }_x\left(y_0,z_0\right)=\frac{{\sum }_{x=x_{\text{min}}}^{x=x_{\text{max}}}\rho \left(x,y_0,z_0\right)V\left(x,y_0,z_0\right)}{\left(x_{\text{max}}-x_{\text{min}}\right)S_{y_0,z_0}}$$

$${\rho }_z\left(x_0,y_0\right)=\frac{{\sum }_{z=z_{\text{min}}}^{z=z_{\text{max}}}\rho \left(x_0,y_0,z\right)V\left(x_0,y_0,z\right)}{\left(z_{\text{max}}-z_{\text{min}}\right)S_{x_0,y_0}}$$

Full text

††pubyear: 2017††pagerange: Mass Transfer and Disc Formation in AGB Binary Systems–A

1 Introduction

More than 50 per cent of the total population of stars are found in binary systems. During most of their main sequence lifetime, they do not exchange mass but when one of the star progresses on to the giant branch [the red giant branch (RGB) or later asymptotic giant branch (AGB) phases], its radius grows large enough that it may interact with the secondary. Mass transfer can ensue.

Single AGB stars will lose mass by driving AGB winds modulated with pulsations (Habing & Olofsson, 2013). Typical mass loss rates and wind speeds are $10^{-7}-10^{-5}$ M*☉* yr*-1and $10-20$ km s-1*, respectively (Herwig, 2005). The standard theory of AGB wind formation can be described as follows: the pulsation in the AGB star’s atmosphere drives shocks to large radii. The post-shock, compressed gas will form dust grains when cooling drops the equilibrium temperature below the dust condensation temperature. In turn, the dust grains will absorb radiation from the AGB star and accelerate. As the accelerated dust grains carry the gas with them, the AGB wind is formed. This theory has found considerable phenomenological success (Bowen, 1988; Liljegren et al., 2016). A lingering question is where does the wind pulsation come from? The answer likely involves convection (Freytag & Höfner, 2008; Arroyo-Torres et al., 2015) and radiation transfer. Convection is a complex problem and in the present absence of a detailed convective AGB stellar atmosphere model we adopt a phenomenological AGB pulsation wind model in our simulations.

Gas ejected by the AGB wind in binary can be accreted by the secondary. Gas that is not accreted can remain gravitationally bound to the binary system or eventually escape. There are three important length scales in this context: the radius of the giant star, the dust formation radius, and the binary separation. Different combinations of these length scales lead to fundamentally different mass transfer mechanisms. If the binary separation is much larger than the dust formation radius, the mechanisms of accretion on to the secondary will be BH accretion. If the binary separation is slightly larger than the dust formation radius, the accretion process will be wind Roche lobe overflow (WRLOF). And if the binary separation is less than the dust formation radius, Roche lobe overflow (RLOF) may incur. When both stars’ Roche-lobes are filled with gas, a common envelope (CE) phase will commence.

BH accretion is perhaps the most thoroughly studied of the aforementioned accretion mechanisms (Edgar, 2004). The classic BH problem (Bondi & Hoyle, 1944) has an accreting object in parallel flow and predicts that the gas within a cylinder will be focused towards the back of the object. The flow rate within that cylinder is the upper limit of the accretion rate. This model is suitable for an isolated star travelling through interstellar medium or X-ray binaries (Pfahl et al., 2002). In intermediate separation binary systems, gas flow (wind) become non-parallel and the full picture of binary systems is needed. Huarte-Espinosa et al. (2013) considered the wind trajectory and simulated the 3D wind-capture disc in binary systems. However, they solved the problem in a corotating box just around the secondary. In this paper, we focus on the holistic processes involving both the mass-losing star and the accreting star. The dynamics of the wind will be self-consistently solved.

Podsiadlowski & Mohamed (2007) put forward the concept of WRLOF. They argued that WRLOF could be a very efficient way of mass transfer mechanism in binaries. In a more general picture, Dermine et al. (2009) concluded that the Roche lobe can be substantially modified by radiation force. They also realized that the Roche lobe model may not be useful at all in very luminous systems (e.g. AGB stars). Despite the possible intricacy, WRLOF model was applied to study the population of carbon-enhanced metal-poor stars and got positive results (Abate et al., 2013). However, a more precise and global picture of WRLOF is still needed. Specifically, we treat the binary stars and the AGB wind as an integrated system and study their mutual influence. That is the interaction we want to address.

RLOF only happens when the Roche lobe surface is connected and WRLOF happens when the Roche lobe surface is detached. Schuerman (1972) illustrated the boundary of the two scenarios with a critical value $\delta_{\text{c}}$. When the ratio of radiation force to gravity is larger than $\delta_{\text{c}}$, the Roche lobe surface is detached (Dermine et al., 2009). To estimate the size of the Roche lobe, Eggleton (1983) approximated the radii of the Roche lobe with a smooth function that became useful in many subsequent studies (Staff et al., 2016a). RLOF model was also applied to X-ray binaries by Savonije (1978) and was able to match the observation. A pioneering hydrodynamical study on RLOF was performed by Sawada et al. (1986). Their 2D hydrodynamical simulations revealed the complex flow pattern in a close binary system.

The concept of CE evolution was formally proposed by Paczynski (1976). It is a short lived and rapidly evolving (e.g. rapidly decreasing orbital separation, ejection of the envelope or tidal destruction) phase (Ivanova et al., 2013) in binary evolution.The progenitors of CE and planetary nebulae (PNe) are possibly very similar and a link has been confirmed observationally (Hillwig et al., 2016). However, the CE phase is still poorly understood and attracting increasing attention (Passy et al., 2012; De Marco et al., 2016). Nordhaus & Blackman (2006), Ivanova & Nandez (2016) and Nandez & Ivanova (2016) investigated the dynamics of in-spiral of the companion and the ejection of the envelop during the in-fall. Staff et al. (2016b) performed a 3D hydrodynamic simulation and validated the in-fall process of a $10M_{\text{J}}$ planet into the envelope of an RGB star and an AGB star.

As a consequence of the aforementioned four mass transfer mechanisms, binary systems are the progenitors of many other objects.

PNe are one class of successor of binary systems. Various PNe shapes and shaping mechanism were discussed(Balick & Frank, 2002). It is believed that disc and tori are common in the progenitors of PNe (van Winckel, 2003; de Ruyter et al., 2006; van Winckel et al., 2009; Bujarrabal et al., 2013) and they are important in the evolution of PNe. A pioneering work in validating this assumption was done by Mastrodemos & Morris (1998) and continuous efforts have been devoted to study the formation of circumbinary discs in AGB binary systems (Akashi & Soker, 2008; Akashi et al., 2015; Chen et al., 2016b). Besides PNe, Type Ia supernovae are also realized to be important successors of binary systems (Whelan & Iben, 1973; Iben & Tutukov, 1984). Researches on Type Ia supernovae focus on the explosion more than on the initial condition and a satisfactory answer can only be found with realistic binary models.

For our binary calculations, we use ASTROBEAR 111https://astrobear.pas.rochester.edu/ (Cunningham et al., 2009; Carroll-Nellenback et al., 2013). ASTROBEAR is a magneto-hydrodynamic code that has adaptive-mesh-refinement (AMR) capability. AMR allows us to only refine the region that has the highest variable (density, pressure etc.) gradients. Every additional level of AMR means the doubling of resolution of every dimension. Thus the effective resolution of a 3D simulation with n levels of AMR is $8^{n}$ times higher than the same one without AMR but the with the same base grid. This technique can speed up the binary calculation by more than $10^{5}$ times compared to a uniform grid code. We perform our calculations on University of Rochester’s BlueHive supercomputer cluster 222https://info.circ.rochester.edu/. Each binary model (without a mass loss ’burst’) takes roughly 240 h on 240 cores and the binary model with a burst takes 300 h on 240 cores. The equivalent total computational time for the six simulations is $4.18\times 10^{5}$ CPU hours.

In this paper, we present the first global 3D radiation hydrodynamic calculation of mass transfer in AGB binary systems. Section 2 contains the equations and physics in our calculation. A description of the isolated AGB star is also included in section 2. Section 3 will present the binary simulation results: the affected mass loss rate , the accretion rate and the morphology of the outflow. We also qualitatively analyse and classify the simulations in Section 3. Section 4 has the conclusions and discussion.

2 model description

We now present the equations used in our calculation and explain the underlying physics and approximations. In what follows, we describe both our isolated AGB wind model and the implementation of the ’burst’.

2.1 Equations and physics

2.1.1 Equations

We use radiation hydrodynamic equations of an ideal gas:

[TABLE]

In above equations, $\rho,n,\bm{u},p,T,E,\dot{Q}$ (see section 2.1.5) are the density, particle number density, velocity, pressure, temperature, total energy and cooling rate in the corotating frame respectively. $k_{\text{b}}$ is the Boltzmann constant. The subscript ni stands for non-inertial terms due to carrying out the calculation in the rotating frame.

[TABLE]

$m_{H}=1.673\times 10^{-24}$g is the mass of hydrogen atom and $\mu=1.3$ is the mean atomic weight. $d$ is the initial binary separation. The origin is set at binary system’s centre of mass. The unit vectors in their explicit forms are:

[TABLE]

The quantity $\bm{f}_{\text{g}}$ is the binary stars’ force on the gas; $\bm{f}_{\text{rad}}$ is the radiation force on the gas; $M_{\text{AGB}},M_{\text{sec}}$ are the masses of the AGB star and the secondary star, respectively; $\alpha$ is a coefficient which will be explained in section 2.2; $\kappa$ (see section 2.1.4) is the Planck mean opacity that changes when dust forms; $L$ is the AGB star’s luminosity; $c$ is the speed of light; $\tau$ is the optical depth (see section 2.1.3); $e$ is the internal energy of the gas; and $\gamma=\frac{5}{3}$ is the ratio of specific heat.

2.1.2 Gravity

We solve the gravity interaction using Newton’s laws. We include the gravitational influence of the stars on the gas and the stars on each other but the stars are not influenced by gravity of the gas in our model. The calculation would be much less efficient and the model would be more complicated were we to take the gas self-gravity and its influence on the stars into account. The gravitational feedback from gas to stars is, however, unimportant on time-scales of our simulation ($t\sim 100$yr). It may be important over longer times such as that of the stellar lifetime. Gravitational interaction between gas and stars can, for example, contribute to the binaries’ orbital decay. We will discuss the binaries’ orbital change problem in a subsequent paper.

2.1.3 Radiation transfer

The typical luminosity of an AGB star is $\sim 10^{37}$erg s*-1*. In our model, the AGB star’s luminosity is $9.13\times 10^{36}$erg s*-1or $L=2342$L☉*. This value is determined by calculating the radiation power at the hard pulsating sphere that underlies the photosphere i.e. $L=4\uppi\sigma_{\text{sb}}(0.9\text{au})^{2}\times(2900\text{K})^{4}$. If the secondary is a low-mass main-sequence star, its luminosity is $L_{\sun}$, which is small compared to the AGB star’s luminosity. If the secondary is a white dwarf (WD), its luminosity will be close to the accretion luminosity. The upper limit of the accretion luminosity will be

[TABLE]

The maximum stable accretion rate measured in our calculation is $\sim 1.2\times 10^{-7}$ M*☉* yr*-1*(table 4) and the maximum secondary’s mass is $0.5$ M*☉. Thus the accretion luminosity will be no greater than $200$L☉*. We note that any jets that form from the secondary will further decrease the available energy in luminosity (Soker & Livio, 1994). With all these estimations, it is reasonable to consider the AGB star as the only radiation source in the binary system. Models that incorporates the feedback from the secondary can be studied in the future.

We apply 2D ray-tracing algorithm to calculate the radiation transfer. We notice that the binary system has periodic or symmetric outflows due to the circular orbital motion, so we average the radiation over the azimuthal angle where the axis of symmetry is crossing the centre of the AGB star and parallel to the z-axis. We assume single interacting photons in our calculations, i.e. photons will only interact with the dust or the gas once. In reality, there could be re-emitted photons and they could still interact with the dust or the gas but neglect that complication at present.

In equation (7), the optical depth $\tau$ is calculated by:

[TABLE]

[TABLE]

where $r_{\text{ps}}$ is the photosphere’s radius, $r_{0}=\left|\bm{r}-\bm{r_{\text{AGB}}}\right|$, $\theta\in\left[0,\uppi\right]$ is the polar angle, $\phi\in\left[0,2\uppi\right]$ is the azimuthal angle and $\kappa=8$ cm2 g*-1if dust is present or $\kappa=2\times 10^{-4}$ cm2* g*-1*if dust is not present (see section 2.1.4).

2.1.4 Dust formation and opacity

Dust can absorb and scatter photons effectively. When the photons leave the photosphere, they enter the optically thin regime and we use the Planck mean opacity calculated as (Chandrasekhar, 1950):

[TABLE]

where $\kappa_{\text{dust}}\left(T_{\text{bb}}\right),f(\lambda,T_{\text{bb}}),\lambda,T_{\text{bb}},\theta_{\text{s}},\sigma_{\text{abs}},\sigma_{\text{sca}},m_{\text{p}}$ are the temperature-dependent Planck mean opacity of the dust, energy distribution of blackbody radiation, wavelength of the photon, blackbody temperature, scattering angle, absorption cross-section scattering cross-section and the mass of the scattering particle, respectively. The scattering for dust is forward scattering, that is, $\sigma_{\text{sca}}(\lambda,\theta_{\text{s}})$ approaches [math] rapidly with the increase of $\theta_{\text{s}}$. Thus the scattering only contribute a small fraction of Planck mean opacity in dust and we can ignore it for simplicity.

The cross-section $\sigma_{\text{abs}}$ can be solved with Mie theory (van de Hulst, 1957; Mätzler, 2002) but it is dust species and dust size dependent. To carry out a quantitative study, we assume the dust is spherical and the size has the distribution of $\text{d}n\sim a^{-3.5}\text{d}a$ between a dust radius of $0.1-0.3\mu$m where $a$ is the dust radius and $n$ is the number density (Mathis et al. (1977), MRN distribution hereafter). We take the most common species of dust around oxygen-rich AGB stars (olivine and pyroxene Jaeger et al. (1994); Dorschner et al. (1995)). In table 1 we present the Planck mean opacity for different temperature black bodys.

The opacity is decreasing with the decreasing radiation temperature. Note that the black body radiation temperature is not necessarily equal to the equilibrium temperature of the dust. The latter is determined by the total flux of radiation while the former is not. In our simulation, the temperature at the photosphere is $2684$K ($T=\left(\frac{L}{4\uppi\sigma_{\text{sb}}(1.05\text{au})^{2}}\right)^{1/4}$), if we take the atmospheric reddening into consideration, the ’effective blackbody temperature’ could be lower but that is beyond this paper’s discussion. The total opacity of the gas-dust mixture is calculated by:

[TABLE]

where $\rho=\rho_{\text{gas}}+\rho_{\text{dust}}$ is the density used in the radiation transfer equations. $\kappa_{\text{gas}}=2\times 10^{-4}$ cm2 g*-1*, and $\rho_{\text{dust}}/\rho=0.01$ is used in our calculation (Bowen, 1988). Given our assumptions, we choose $\kappa_{total}=8$ cm2 g*-1*. This value may be slightly larger than the real opacity of gas-dust mixture around an AGB star but it lies within a reasonable range for the calculations.

The dynamic dust formation is also included in our simulation. Dust forms when the equilibrium temperature is below the dust’s condensation temperature; we calculate the equilibrium temperature by

[TABLE]

Where $\beta=0.3,\sigma_{\text{sb}},r$ are the coefficient of absorption, Stefan-Boltzmann constant and the distance to the centre of the AGB star respectively. $\tau$ will be calculated by equation 15. We add $\beta$ in our equation because the dust can only absorb short and specific wavelength of radiation. We choose this value because it gives us reasonable AGB winds (section 2.2). Finally, the condensation temperature of dust is assumed to be $1500$K.

2.1.5 Cooling

A very important physical ingredient in binary simulation is cooling. The accretion on to the secondary star and the strong shocks around the AGB star will compress and heat the gas. Disc formation and accretion rates are dependent on cooling (Lee & Chiang, 2015).

We use the same cooling mechanisms as Mastrodemos & Morris (1998): electron excitation of H i and molecular cooling between molecular hydrogen and water. We updated a few input parameters to make the cooling more realistic.

When outside the photosphere and when $T>4000$K, the electron excitation of H i is considered. The cooling rate is calculated by $\dot{Q}=7.3\times 10^{19}n_{\text{e}}n_{\text{H}}\exp\left(-118400/T\right)$ (Spitzer, 1978). The number density of electrons is characterized by a step function to further distinguish high temperature cooling and low temperature cooling. For $4000<T<8000\text{K}$, $n_{\text{e}}=6.0\times 10^{-4}n_{\text{H}}$ and for $T>8000$K, $n_{\text{e}}=1.0\times 10^{-3}n_{\text{H}}$.

Molecular cooling is more complicated. In reality, the formation of molecules requires both low temperature and consideration of relaxation times. H2 is scarce in the constantly shocked AGB atmosphere since the timescale of H2 formation is longer than the time between successive shocks (Bowen, 1988; Cherchneff, 2006). We empirically choose the H2 formation region to be close to the dust formation region, thus molecular cooling will not be considered within a $1.8$au radius of the AGB star. Outside the $1.8$au sphere, when $T<4000$K, the molecular cooling is calculated by following Neufeld & Kaufman (1993) and Mastrodemos & Morris (1998).

2.2 Isolated AGB star

Before discussing the binary model we first describe our representation of the isolated AGB star. The test results of normal pulsating AGB star and the AGB star with burst will be presented in separate subsections. Our isolated AGB star model is phenomenological and mainly inspired by Bowen (1988) and Liljegren et al. (2016). A sketch of the AGB star’s structure is shown in Fig. 1 and its physical properties are listed in Table 2.

There are three important shells in the isolated AGB star model: the ’hard’ pulsating inner boundary, the photosphere and the dust formation zone. The radius of hard pulsating inner boundary is the hydrodynamic boundary fixed in radius, density and temperature but varying in velocity. It mimics the convective instability in the atmosphere (Freytag & Höfner, 2008) and drives shocks (or potentially a burst) outwards. This boundary is ’hard’ because it will not be affected by the condition on the outside and therefore the dynamics of this shell are predetermined. The photosphere is presumed to be shell with inner edge $r=1.05$au. In reality, the photosphere boundaries should be dynamic but that requires a more detailed radiation transfer calculation (and computational resources) than we attempt here. Between the pulsating boundary and the photosphere we use a reduced gravity described as

[TABLE]

where $\alpha=0.8$. This reduced gravity could come from a combination of radiation pressure and instabilities. In our model testing experience, different $\alpha$ value will give us different wind properties and there is no simple linear relation of $\alpha$ to wind speed and mass-loss rate. For example, in the small $\alpha$ situation, the pulsation at the hard pulsating sphere may push more material to the photosphere and dust formation zone with higher speed but the optical depth will increase dramatically at the dust formation zone. Then the radiation pressure may drop fast with the increase of radii and the wind may decelerate. How useful is $\alpha$ and how to justify $\alpha$ more rigorously remains a question in phenomenological AGB wind models. We also allow for the atmospheric reddening in our binary model using $\kappa_{gas}=2\times 10^{-4}$ cm2 g*-1*between the pulsating boundary and the outer edge of photosphere. Outside the photosphere, the gas opacity is kept the same with $\alpha=1.0$ and the cooling from electron excitation turned on. Finally, we have the dust formation zone whose inner edge is dynamic due to the atmospheric reddening. In the 3D isolated AGB star simulation, the dust formation radius is between $1.8\text{au}-2.0\text{au}$.

We now present the two AGB sub-models used in the simulations.

2.2.1 Pulsation

We will always assume pulsation in our binary simulations and we use it to solve the isolated spherically symmetric pulsating AGB wind problem in 3D with ASTROBEAR. The physical dimension of the simulation box is $(24\text{au}^{3}$. The base resolution is $240^{3}$ and we add one level of AMR; therefore, the physical dimension of smallest cell is $(0.05\text{au})^{3}$. The numerical solution for the outflow from one of our pulsating isolated AGB stars is shown in Fig. 2.

The upper panel of Fig. 2 shows that the wind velocity approaches $15$km s*-1*asymptotically. The flatness of the velocity profile implies that the radiation pressure on the wind balances the gravitational force. This balancing is dynamic because the optical depth (Eq. 15) changes within one period of pulsation. Thus the radiation force (Eq. 7) also changes. Between the surface of the pulsating boundary and the dust forming region, the velocity can go negative, so fall back can occur. In this region, gas is constantly shocked as material rises and falls. The escape velocity to the AGB star’s gravity is shown with dashed line.

In the lower panel of Fig. 2 we see the radial velocity at the pulsating boundary which is $v_{p}=5\cos{\left(2\uppi t/P\right)}$km s*-1where $P=365$d. The mass loss rate is measured at $10$au, that is, we summed the mass flux through the $10$au shell. Thus we expect to see the mass loss rate out of phase with the pulsation. We can however see from Fig. 2 that the mass loss rate has the same period as the pulsation. The average mass loss rate during normal pulsation is $\dot{M}_{\text{avg}}=2.31\times 10^{-7}$ M☉* yr*-1*. This is smaller than the mass loss rate predicted by the empirical formula given by van Loon et al. (2005). Their formula predicts a mass loss rate of $2.57\times 10^{-6}$ M*☉* yr*-1*. Our lower mass loss rate is mainly a consequence of not constructing a fully self-consistent atmosphere which would at least require a more accurate equation of state and radiation transfer. Enough resolution should also be contained within each scale height. Since the intention of our present paper is to study wind binary interactions we only need to get the general characteristics of the wind correct (dust driving etc).

2.2.2 Burst

In our isolated AGB star model, a ’burst’ is an increase in mass outflow over a short period of time of the AGB star. In reality, the outflow could be induced by the in plunge of a planet with $M_{\text{p}}\in\left[M_{\text{J}},50M_{\text{J}}\right]$ (Nordhaus & Blackman, 2006; Nordhaus et al., 2007; Nordhaus & Spiegel, 2013). 3D hydrodynamic calculation (Staff et al., 2016b) shows that a $10M_{J}$ planet could unbind $10^{-4}-10^{-3}$ M*☉of the giant star’s envelope and increase the giant’s luminosity. Ivanova & Nandez (2016) use 3D hydrodynamic model and shows that a $0.05$ M☉planet could unbind $0.02$ M☉*of the giant’s envelop. They both see ’fallback’ phenomenon (Chen et al., 2016a) in the simulations. In the hope of testing the dynamics of such burst in binary systems, we build our phenomenological model to mimic their solutions. The details of plunge in process of a planet is not included in our phenomenological model, thus the ’burst’ exists numerically. Next we will explain how we induce the ’burst’ in our isolated AGB star model. It remains the same in the binary simulation.

The total span of the burst is $5$yr. During the burst, the outflow density at the hard pulsating sphere will change from $1.2\times 10^{-8}$ to $1.2\times 10^{-12}$ g cm*-3*, the outflow velocity which also lies at the hard pulsating sphere will vary with time:

[TABLE]

where $30<t<35\text{yr}$. This expression is chosen to be generic and any smooth function with similar shape would work. The outflow temperature will increase from $2900$ to $3450$K and the luminosity will be doubled to $4684$L*☉*. We also stretch the reduced gravity region to $2$au from $1.05$au. Our goal in driving the burst is to have an outflow that captures some of the key features of what would be found in more from more detailed models. The hydrodynamics deep inside the burst is not our primary interest and we chose parameters that heuristically allow the wind to increase its mass lose rate. Fig. 3 presents the numerical solution of the burst.

Like the pulsating AGB model, the burst model is also solved with ASTROBEAR in 3D. The burst lasts from $30$ to $35$yr. After the burst, the AGB star returns to the normal pulsating state. We still measure the mass loss rate at $10$AU, so there is lag between the burst and the peak of the mass loss. The mass-loss rate ramps up to the peak ($2.3\times 10^{-4}$ M*☉* yr*-1*) rapidly then decays and finally some leftover material starts to fall back. The total mass loss driven by the burst is roughly $2.9\times 10^{-4}$ M*☉*.

2.3 Accreting star

The accreting secondary is the other important object in binary simulations. We first consider the range of physics which might be included in including the secondary. As mentioned before, a low-mass-main-sequence stellar accretion disk is not luminous, but it might have an extended magnetic field. In that case, the magnetic field and the disk could interact with each other. If the secondary is a WD, the radiation from the disk could be strong enough to change the dynamics of the accretion process, especially for possibly super-Eddington accretion in the RLOF. The temperature of the innermost region of the disk can be up to $\sim 10^{6}$K and the lowest temperature in an AGB binary system can be $\sim 10$K. Thus one needs a continuous tablular equation of state from $10-10^{6}$K (Pejcha et al., 2016; Tomida et al., 2015). The density is also spans many orders of magnitude ($10^{-20}-10^{0}$ g cm*-3*). The possible emerging outflows could change the morphology of the binary system and the accretion disk (Blackman et al., 2001; Frank & Blackman, 2004; Blackman & Pessah, 2009). Diskoseismology (Nowak & Wagoner, 1991) is relevant if self-gravity is to be considered in the binary simulation. All these questions converge to the question of how to treat the inner boundary condition of the accreting star.

Given resolution constrains (even with AMR) in our calculation the accreting star is a gravitating sink particle and we adopt the Krumholz et al. (2004) prescription. The sink particle has been tested via the BH accretion problem where the accretion rate, numerical viscosity and the accretion flow under various conditions were examined and found to produce excellent agreement with analytic results. One must be careful in that this prescription for sink particles can produce errors in the accretion rate when the flow is transitioning from subsonic to supersonic at the accretion radius.

2.4 Outer boundary condition

For our outer boundary we carve out a cylinder with its axis of symmetry parallel to the z-axis inside the cuboid. The radius of the cylinder is one cell smaller than the half length of the x or y dimension. In the lab-frame, we let the region outside the cylinder have high Mach number supersonic outward flows. The outward flow condition is reset before every step of the simulation. Then we transform the Galilean variant quantities ($\bm{u},E$) of the outside cylinder region to the corotating frame. The advantage of this boundary condition is that it prevents unphysical inflow from entering the grid due to the frame rotation. For the top and bottom surfaces of the cylinder, the physical quantities in the corotating frame ($\rho,\bm{u},p,E,e,T$) outside the boundaries are first extrapolated. The perpendicular component of the velocity is then calculated as

[TABLE]

where $\hat{\bm{n}}_{b}$ is the unit outward normal vector of that boundary. If $u_{p}\geq 0$, we do nothing; if $u_{p}<0$, we set it to [math]. This prescription eliminates unphysical fallback of material into the simulation.

3 binary simulations

3.1 Numerical setup of binary simulations

The setups of the six runs are listed in table 3. The column named AMR represents the highest level of AMR in the simulation. This highest level is only applied to the region around the secondary in the ’no burst’ models. The AGB star is kept at the same resolution as in the isolated model. We do this to focus the computational power on resolving the accretion disk around the secondary to the finest level. For the binary model with the burst, we do not use the same level of AMR as model 2 so that we can compare the accretion rate before burst to model 2. We find that their accretion rate is nearly the same, this feature ensures that our resolution around the secondary does not influence the accretion rate, i.e. we have enough resolution.

We now explain the models in the table by taking model 1 as an example. For model 1, we simulate a binary system with $1M_{\sun}$ AGB star and a $0.1\text{M}_{\sun}$ companion. The separation of the binary stars is $3$au. The physical dimension of the simulation box is $(40\text{au})^{3}$. We divide the space into $100^{3}$ base cells so that each base cell has a physical volume of $(0.4\text{au})^{3}$. With 4 levels of AMR, the physical volume of the finest cell is $(0.025\text{au})^{3}$. The region with highest level of AMR is restricted to the Roche lobe of the secondary. The resolution of the AGB star is the same as in section 2.2.1, i.e. the physical volume of the finest cell is $(0.05\text{au})^{3}$. It is assumed that the secondary will accrete the gas within $r_{acc}=0.1$au. The AGB star is in its normal pulsating state throughout the simulation.

We do not deduct the mass lost from the AGB wind from the mass of the star in our calculation because the total mass loss in the entire simulation is small compared to the mass of the AGB star. The total mass loss $\text{d}M$ in a $100$yr simulation is $<3\times 10^{-5}\text{M}_{\sun}$. This value will not exceed $3\times 10^{-4}\text{M}_{\sun}$ even in the simulation with burst. We expect that such a small mass change will not signficantly alter the dynamics of the AGB wind, the accretion of the secondary, or the outflow morphology or orbital dynamics. For the same reason, the angular frequency in the corotating frame simulations and the center of mass are fixed.

For longer simulations (e.g. $>>100$yr) than considered here however, the mass loss will be important over the lifetime of the AGB star ($M_{\text{AGB}}/\dot{M}_{\text{AGB}}$) and is a topic of future work.

3.2 Dynamics of the mass loss rate and the accretion rate

In our calculations, we find that the secondary influences the mass loss rate from the primary and that the accretion rate on to the secondary varies dramatically with different binary separations, $d$. The results are listed in table 4 where $\zeta=\frac{\dot{M}_{\text{acc}}}{\dot{M}_{\text{AGB}}}$ and $\epsilon=\frac{\dot{M}_{\text{AGB}}-\dot{M}_{\text{avg}}}{\dot{M}_{\text{avg}}}$.

How do we calculate $\dot{M}_{\text{AGB}}$? The mass lost from the AGB star can go to each of four places:

It can escape the binary system at rate of $\dot{M}_{\text{esc}}$. 2. 2.

It can fall back into a circumbinary disk $(\dot{M}_{\text{cd}})$ 3. 3.

It can fall into an accretion disk around the secondary ($\dot{M}_{\text{ad}}$). 4. 4.

It can fall directly on to the secondary and be accreted at rate of $\dot{M}_{\text{acc}}$.

The mass loss rate of the AGB star is thus the sum of these four contributions:

[TABLE]

When the simulation has come to stable state, the mass in the circumbinary disk and in accretion disk around the secondary are nearly constant so that ($\dot{M}_{\text{cd}}\approx 0$ and $\dot{M}_{\text{ad}}\approx 0$). Therefore:

[TABLE]

The accretion rate $\dot{M}_{\text{acc}}$ is calculated by ASTROBEAR using the prescription introduced by Krumholz et al. (2004) when the simulation is in stable state. The detail algorithm can be found in section 2.4 in Krumholz et al. (2004). When the simulation has reached its stable state,$\dot{M}_{\text{acc}}$ does not vary much (see the first 30 years in Fig. 7) during pulsation but we still take the average value over 10 pulsation period. $\dot{M}_{\text{esc}}$ is calculated by summing the flux through a shell of radius $r_{\text{flux}}$ (see table 3) centered at the origin. Since the flux and the accretion rate vary with the pulsation, we take the average value over 10 pulsation periods when the binary system has come to stable state.

As defined above, $\epsilon$ measures the change in mass loss rate compared to the isolated AGB model. It is evident from the table that the AGB mass loss rate is different when there is a companion. There are at least three ways that the secondary can affect $\dot{M}_{AGB}$:

The secondary can attract material around it to make the optical depth rise. The reduces the radiation force on the radiation-shielded gas and some of it will fall back toward the AGB star. 2. 2.

The secondary may not be able to accrete material at its inner boundary as fast as it accretes material at the outer boundary of its disk. Material will then pile up in the accretion disk. Some of it will flow through the L2 point (Pejcha et al., 2016) and leave the system. 3. 3.

The secondary can pull some of the material inside the dust fomrmation radius (i.e. before it is radiatively accelerated) that would otherwise fall back on to the AGB star in the isolated model.

The first mechanism above (optical depth mechanism) reduces $\dot{M}_{\text{AGB}}$ while the second (L2 mass loss mechanism) and the third mechanisms increase $\dot{M}_{\text{AGB}}$. In model 1, 2 and 3, all three mechanisms operate; L2 and L3 mass loss are found in model 1 and 2 and L2 mass loss is found in model 3 (Fig. 4). For model 4 and 5 the third mechanism above dominates via BH accretion on to the secondary.

The accretion rate varies by a factor of $20$ between different models. In fact the accretion rate can be up to $\zeta=39.2$ per cent (model 1) of $\dot{M}_{\text{AGB}}$ in our simulations. In that model, the dust formation radius is $r_{\text{dust}}\approx 1.8$au. The pulsation makes the optical depth oscillate and thus the radiation force from the AGB star also oscillates. The AGB star’s Roche lobe deforms and can detach from that of the secondary under this oscillating radiation force and the L1 point may vanish in some circumstances (Dermine et al., 2009). Given that the typical BH accretion gives $\zeta\ll 10$ per cent whereas RLOF can produce $\zeta\approx 100$ per cent, we interpret this high but not $100$ per cent mass transfer efficiency to exemplify WRLOF. For models 4 and 5 respectively, $\zeta=3.25$ per cent and $\zeta=2.24$ per cent. These values and the spiral outflow wind from the binary (Fig. 4) and suggest that models 4 and 5 exhibit BH accretion.

3.3 Results from the ’no burst’ binary simulations

We present the figures from model 1 - model 5 in Fig. 4. The figures show results from simulations that which achieve a stable state for at least $10$yr.

3.3.1 Bipolar outflow

The low mass secondary in orbit around AGB star can shape the AGB wind into a bipolar outflow structure (Nordhaus & Blackman, 2006). An observational example is L2 Puppis (Kervella et al., 2015; Kervella et al., 2016). Chen et al. (2016b) found that 3D hydrodynamic simulation of L2 Puppis and found very good agreement with the morphology and spectral energy distribution. That work invoked an outburst from the AGB star. Here we see that model 1 and 2 conspicuously show a bipolar outflow even without the outburst. Repeating extended lobes are found both in the north and south polar directions. The lobes can extend to $10-20$au with ripple-like structures caused by the pulsation. These structures can also be affected by the binaries’ orbital motion when we compare model 1 and model 2. In model 2, the secondary is more massive so its gravity has a more conspicuous imprint on polar outflows. In model 1 the ripples are more symmetric because the secondary mass is only $0.1\text{M}_{\sun}$ and the dynamics of the polar outflow is dominated by the radiation force from the AGB star.

3.3.2 Circumbinary disc and accretion disc

Circumbinary disc are common around AGB and RGB stars (Bujarrabal et al., 2013; Kamath et al., 2016). In our simulations, conspicuous circumbinary disc or tori have been found in models 1, 2 and 3. The circumbinary disk in model 2 resembles the simulation in Mastrodemos & Morris (1998). The circumbinary disk in model 1 is optically thick in the equatorial plane and the material can fall back to the AGB star. The gap of the circumbinary disk in model 1 resembles a pioneer SPH simulation done by Artymowicz & Lubow (1996). The circumbinary disc in all the simulations are dynamically forming. They expand outward and are replenished at their inner radii. Models 4 and 5 do not show circumbinary disc but instead exhibit spiral structure.

Accretion disc are also found in models 1, 2 and 3. The x-plane cuts of the three models show that the accretion disc are stable in the equatorial plane and are thin disc. The accretion disc do not exhibit ’flip-flop’ instability in these three models. Possible reasons may be that the mass transfer mechanism in these models is not BH accretion (Matsuda et al., 1991; Ruffert, 1996; Bate, 1997; Bate & Bonnell, 1997; Foglizzo & Ruffert, 1999), 3D simulation can suppress the instability (Soker, 1990, 1991) and some physical process (i.e cooling, it somewhat comparable to making $\gamma$ closer to 1) is removing the excess energy thus may also suppress the instability. In all three models, the accretion disc are filling their Roche lobes and losing mass via the L2 point. The gas leaking from the L2 point merges into the circumbinary disc and can be pushed to higher orbit or escape. However, in models 4 and 5, no obvious accretion disc are found.

Huarte-Espinosa et al. (2013) carried out local simulations around the accreting star. finding $\zeta\approx 3$ per cent and accretion disc. Comparing our models 4 and 5 to their simulations, we identify three reasons that conspire to prevent accretion disk formation:

Although we have similar $\zeta$ values, our AGB wind mass loss rate ($\dot{M}_{\text{avg}}=2.31\times 10^{-7}$ M*☉* yr*-1*) is much smaller than theirs ($10^{-5}$ M*☉* yr*-1*). 2. 2.

Our AGB wind velocity is higher and the mass of the secondary is smaller in our simulations than theirs. Both result in less captured gas. 3. 3.

The azimuthally averaged radiation force from the AGB star opposes accretion disk formation.

3.3.3 Spiral structure

Spiral structure around AGB stars has been discovered observationally (Maercker et al., 2012; Mauron & Huggins, 2006; Kim et al., 2015) and validated by simulations (Mastrodemos & Morris, 1998; Edgar et al., 2008). Figures in models 4 and 5 also successfully reproduce the spiral structure which is similar to CIT 6 and R Sculptoris. The spiral structure in our simulation is caused by BH accretion. In BH accretion, gas is focused behind the secondary due to under the gravitational force and cooled quickly due to the increased density. The gas conserves its outward momentum and angular momentum but loses its internal energy through cooling. The cold gas can keep its sharp density gradient for very long time. The gas is focused into the equatorial plane where it sweeps up more gas. Under some circumtances, L2 mass loss can also produce spiral structure (Pejcha et al., 2016).

The $x=0$ mid-plane cut of the spiral structure shows that gas preferentially accumulates in the equatorial plane. The morphology of the $x=0$ plane cut is different from a torus or circumbinary disk. The circumbinary disk and spiral likely represent two cases of a ’continuous’ morphological sequence that depends on radiation pressure. To see this, first compare the $z=0$ plane cut of models 4 and 5, we see that the curvature of the spiral structure in model 4 is greater than that of model 5. This spiral curvature is influenced by the radiation pressure on the spirals; the larger the radiation pressure, the smaller the curvature. The radiation pressure is in turn, determined by the optical depth and density in the equatorial plane. Thus the equatorial plane density is directly linked to the curvature of the spirals. Since the secondary in model 4, with its $8$au separation, focuses more material into the equatorial plane than model 5 (with $10$au separation), the curvature of the model 4 spiral should be larger than that of model 5.

Furthermore, if the secondary focuses so much gas in its tail that the radiation force decreases more steeply with increasing radius than the gravitational force, the radius of the curvature of the spiral may become small enough (comparable to the binary separation), that the head and the tail of the spiral connect. The spiral structure will then become a circumbinary disk. Although we describe this transition between spiral and circumbinary disk as a ’continuous’ sequence determined by the radiation pressure, the interaction between gas focused by the secondary, radiation force, and curvature of the spiral may make this transition range small and highly nonlinear. There remains opportunity to study this question further.

3.4 The burst

The numerical setup of the burst simulation is represented by model 6 in table 3. For this case, the AGB star incurs a burst from $t=30yr$ to $t=35yr$ (section 2.2.2) and maintains its normal pulsating state before and after the burst. To demonstrate the distribution of the gas before and after the burst, we calculate average volume densities along $x$ and $z$ axes respectively. By computing the respective local surface densities and dividing by the domain range of the projected axis over the region of averaging.

Mathematically, we have

[TABLE]

[TABLE]

where $V\left(x,y_{0},z_{0}\right)$ and $V\left(x_{0},y_{0},z\right)$ are the volume of the cubes, $S_{y_{0},z_{0}}$ and $S_{x_{0},y_{0}}$ are the projected area of the columns. $\rho_{x}\left(y_{0},z_{0}\right)$ and $\rho_{z}\left(x_{0},y_{0}\right)$ ($\rho_{x}$ and $\rho_{z}$ hereafter) are the average density of x axis and z axis correspondingly.

Fig. 5 shows the the average density of $x$ and $z$ projections at different times. At $t=30.0$yr, the burst is just about to start and the binary in the same state as model 2. We can therefore compare this state with the mid-plane-cut figures of model 2. The figure in the left column of the first row shows $\rho_{z}$, and reveals the accretion disk around the secondary and the circumbinary disk. There are two arms stretching out from the L2 and L3 points. The shape of the outer boundary is also seen. The figure in the left column of the second row shows $\rho_{x}$ in which the accretion disk and the circumbinary disk are now seen edge-on. The figures also are also spiral arms in the polar direction which look very similar to those in L2 Puppis (Kervella et al., 2015; Kervella et al., 2016).

The first figure in the right column of Fig. 5 shows $\rho_{z}$ and $\rho_{x}$ at the later time $t=35.3$yr. At this time, the burst is over but we can see that the ejected material during the burst was asymmetrically distributed. In the equatorial plane, there is a ring-like structure and the accretion disk is temporarily absent. The polar outflow has a cigar-shaped structure. There is a low density ’hollow’ region around the center in the $x=0$ plane cut residing inside the high density ejecta (see Appendix). The ejected material will continue to expand while gas farther in is falling toward the binary stars, resulting in the cavity. This gas can be accreted by the secondary, re-enter the AGB star’s boundary and disappear from the simulation, or remain in the equatorial plane and form a larger circumbinary disk.

By $t=58.1$yr (third row of Fig. 5) the gas in the polar direction has been falling back for many years. The flow become turbulent and shocks created by the colliding flows are evident throughout the region. The thin accretion disk around the secondary has fully formed. The region very close to the binary star no longer resembles that of model 2 as there are no arms stretching arms from the L2 and L3 points. This region will continue to evolve and eventually return to a stable state.

At $t=67.9$yr, the accretion disk can be seen in both $x$ and $z$ axis projections. There is a circumbinary disk forming and we can see the two arms again stretching from the L2 and L3 points. The center region is now returning to the stable state that resembles this region model 2. However, much of the fall back material is orbiting the binary at larger radii, which is very different from the result of model 2. The polar gas is still falling back on to the equatorial plane and colliding with the AGB wind pulses. This create shocks which cool and lead to filament structure in the density projections. We believe that the simulation is still evolving but following the simulation further was too computational demanding for present purposes: at current resolution, to simulate $1$yr of this model cost 1400 CPU hour.

3.5 Accretion rate versus time

We present the time dependent accretion rate of Model 6 in Fig. 7. Before the burst ($t=30.0$yr), the accretion rate is roughly $\dot{M}_{\text{acc0}}=1.1\times 10^{-7}$ M*☉* yr*-1*and oscillates slightly with the AGB pulsation. This value is the comparable to the accretion rate for model 2 in table 4. Since we use one more level of AMR in model 2 around the accreting star, the similar accretion rate implies that the resolution around the secondary is not influencing this quantity.

During the burst, indicated by the yellow strip in 7, the accretion rate reaches $5.0\times 10^{-7}$ M*☉* yr*-1*, or 5 times of the $\dot{M}_{\text{acc}}$. The ejected material during the burst leaves a small ’hollow’ region within. This causes the accretion rate to drop, reaching as low as approximately half of $\dot{M}_{\text{acc0}}$.

The final fall-back phase is perhaps the most interesting. Here the accretion rate can surge to more than $8.0\times 10^{-7}$ M*☉* yr*-1*or 8 times $\dot{M}_{\text{acc0}}$. As the material continues to fall back, the accretion rate decreases slowly, sustaining for more than $20$yr.

The two peaks and the prolonged decay of the accretion rate may lead to interesting time a of the binary light curves or temporal variability and spatial variability (i.e. knots) in jets. (Lopez et al., 1993; Soker & Livio, 1994; Corradi et al., 1999; Blackman et al., 2001; Balick & Frank, 2002).

4 conclusions and discussion

To study mass transfer in binary systems with an AGB primary, we perform 3D radiation hydrodynamic simulations with ASTROBEAR. We carry out analysis from first principles and examine the physical properties of isolated AGB stars (section 2.2). Both normal pulsating states and bursting states are explained in detail and discussed (van Loon et al., 2005; Ivanova & Nandez, 2016; Staff et al., 2016b). We perform six binary simulations with ASTROBEAR. Five of the binary models have a normal pulsating AGB star and one of the binary models has a bursting AGB star (Sec. 3).

Our isolated AGBs star produce AGB winds with reasonable terminal velocities and velocity profiles (Fig. 2) but low mass-loss rates. We think it is mainly the incorrect EOS and lack of radiative-convective envelop that cause this small mass loss rate. In radiation hydrodynamic binary simulations, EOS and radiation pressure can induce convection and determine density profiles, which in turn sets the optical depth which effects the dynamics and mass-loss rate. Therefore, in principle, different EOS and opacity profile would yield different dynamics and mass loss rate. On the other hand, there is a similarity between high mass-loss rates in wide binaries and low-mass loss rates in close binaries as the density drops with distance. Therefore it is still worthwhile to carry out simulations in binary systems with low mass-loss rates. In particular, self-consistent AGB wind models are important and are of interest.

All binary models with the normal pulsating AGB star exemplify BH accretion and WRLOF. The ratio of accretion rate to mass-loss rate, $\zeta$, and the mass-loss rate enhancement, $\chi$, are calculated for each model. At least three mechanisms are responsible for the non-zero $\zeta$ and $\chi$ values in close binary stars (Sec 3.2). The results from model 1 to model 5 are displayed in Fig. 4. Models 1 and 2 exhibit clear evidence of bipolar outflows. Circumbinary disc and accretion disc can be found in models 1, 2 and 3. Spiral structure is observed in model 4 and model 5. We do not find any obvious accretion disc in model 4 and model 5. By comparing them to Huarte-Espinosa et al. (2013), we list the possible reasons why we do not see the accretion disc in our models. We successfully reproduce all three types of structures (bipolar, circumbinary disc and accretion disc and spiral structure) with the same mass-losing AGB star. We believe the versatility proves the broad correctness of our binary model.

The binary simulation with the bursting AGB star (Sec. 3.4) has the most intricate structure. We present both plots of the middle-plane cuts (Appendix) and projections. At first, it is similar to model 2. We do this with intent, as this demonstrates that the resolution is not adversely effecting the accretion rate and the dynamics. After the burst, the ejecta pushes the ambient material significantly outward, thereby creating a hollow inner region. The accretion rate drops during this time (Fig. 7), afterwhich, we observe significant amounts of gas falling back. During this time, the infalling gas collides and shocks. The high density regions cool faster and conserve their filament structures. At the end of our simulation, we can see the actively forming circumbinary disk and accretion disk at the inner region and the evidence of the burst at large radii in the equatorial plane. The accretion rate shows two peaks, the first one corresponds to the burst enhancement and the second one corresponds to the fall back process. We can see prolonged high accretion rates due to the fall back process. This accretion rate change may provide clues on the formation of jets and knots in binary stars.

Lastly, we address two numerical difficulties in binary simulations.

Two important structures in AGB binary systems are accretion disc and the AGB star itself. Radiation transfer is equally important in both objects. For now, we only resolve the radiation transfer around the AGB star with 2D ray-tracing method because of efficiency. Global 3D radiation transfer will be the best way to resolve the radiation transfer process in close binary systems but is currently limited due to computational complexity. 2. 2.

In corotating simulations, we find that numerical viscous heating is a spurious, adverse effect. Numerical heating is rooted in Eulerian codes as they suffer from the so called ’Galilean non-invariance’ problem (Springel, 2010). Robertson et al. (2010) did a critically study on this problem and claimed that simulations done by Eulerian codes can converge to the reasonable solution if the resolution is high enough. In our simulations, we enlarge the resolved area and increased the level of AMR as the ejecta expand to keep the numerical viscosity low. We do this because the Mach number in the lab-frame at large radii is usually around $10$, but it may increase to $20$ in the corotating frame as the bulk velocity is Galilean variant. If the numerical viscous heating converts $1$ per cent of the kinetic energy to internal energy in $1$yr of simulation, the temperature will rise $4$ times in $1$yr. However, the increasing of temperature will be slowed down when the Mach number is lowered. One obvious effect is that the circumbinary disk may be puffed up and even disappear resulting in no Keplerian disc (de Ruyter et al., 2006; Bujarrabal et al., 2013) which should be observed at large radii. This effect is stronger in close binary simulations as they have higher angular frequency. Cartesian meshes may also worsen the situation. To better understand dynamics of circumbinary disc, we think both Lagrangian fluid codes or Eulerian codes with cylindrical or spherical meshes would be beneficial.

Acknowledgements

We gratefully acknowledege all the staffs in Center for Integrated Research Computin (CIRC) for supporting us with the amazing computer resources, we can not finish the simulations so fast without their help. Special thanks should be given to Dr. Baowei Liu for helping us setting up the code. We also thank Mr. Richard E. Sarkis and Mr. David Munson for technical support. We appreciate the discussion with Prof. Bruce Balick and Prof. Hans van Winckel on AGB binary problems and Prof. James Stone on the ’flip-flop’ instability.

Appendix A the ’Hollow’ in model 6

Figure 8 shows the middle plane cut plots from $t=35.3$yr to $t=58.1$yr. During this time, the gas is falling back and it is more clear to see the structure in middle plane cut plots than the projections.

We can clearly see the ’eye shape’ low density region inside the ejecta at $t=35.3$yr. It is ’closing’ with the fall back of gas. The fall back gas will collide into each other and create shocks. It can also interact with the successive AGB wind as can be seen at $t=58.1$yr.