© ESO 2018

1. Introduction

Depending on their evolutionary state and orbital separation, stars in binary systems can interact in various ways. One of the main interaction processes is the exchange of mass between the stars, which strongly affects the stellar masses, spins, and chemical compositions. Furthermore, as a consequence of mass transfer and the loss of mass and angular momentum from the system, significant changes in the orbital parameters can occur (e.g. van den Heuvel 1994; Postnov & Yungelson 2014). In relatively close binaries mass transfer usually occurs via Roche-lobe overflow (RLOF), while in wide binaries mass transfer by stellar winds can take place. Low- and intermediate-mass binary stars usually interact when the most evolved star reaches the red-giant phase or the asymptotic giant branch (AGB). This is a consequence of (1) the large stellar radius expansion during these evolution phases, resulting in a wide range of orbital periods for RLOF to occur, and (2) the strong stellar-wind mass loss of luminous red giants. During the final AGB phase the mass-loss rate is of the order of 10⁻⁷–10⁻⁴ M _⊙ yr⁻¹ (Höfner 2015) and the velocities of the winds are very low (5–30 km s⁻¹), which allows efficient wind mass transfer even in very wide orbits.

Several classes of relatively unevolved stars show signs of having undergone such mass transfer in the form of enhanced abundances of carbon and s-process elements, which are the nucleosynthesis products of AGB stars. These include barium (Ba) and CH stars (Keenan 1942; Bidelman & Keenan 1951), extrinsic S stars (those in which the radioactive element technetium is absent, Smith & Lambert 1988), and carbon-enhanced metal-poor stars enriched in s-process elements (CEMP-s stars, Beers & Christlieb 2005). These stars are thought to be the result of the evolution of low- and intermediate-mass binaries that have undergone mass transfer during the AGB phase of the more massive star onto the low-mass companion, which is the star we are currently observing. The erstwhile AGB star is now a cooling white dwarf (WD), which is in most cases invisible, but which reveals itself by inducing radial-velocity variations of the companion (McClure et al. 1980; McClure 1984; Lucatello et al. 2005). These systems have orbital periods that lie mostly in the range 10²–10⁴ days and often show modest eccentricities (Jorissen et al. 1998, 2016; Hansen et al. 2016). They have these orbital properties in common with many other types of binary systems that have interacted during a red-giant phase, such as post-AGB stars in binaries (van Winckel 2003), symbiotic binaries (Mikołajewska 2012) and blue stragglers in old open clusters (Mathieu & Geller 2015).

These orbital properties are puzzling when we consider the expected consequences of binary interaction during the AGB phase. Interaction via RLOF from a red giant or AGB star in many cases leads to unstable mass transfer (Hjellming & Webbink 1987; Chen & Han 2008), resulting in a common envelope (CE) phase that strongly decreases the size of the orbit (Paczynski 1976). In addition, the orbit is expected to circularise due to tidal effects even before RLOF occurs (Zahn 1977; Verbunt & Phinney 1995). In the case of wind interaction, if the gas escapes isotropically from the donor as is commonly assumed, a fraction of the material will be accreted via the Bondi-Hoyle-Lyttleton (Hoyle & Lyttleton 1939; Bondi & Hoyle 1944; hereinafter BHL) process and the rest of the material escapes, removing some angular momentum and widening the orbit. Taking into account these two main mass-transfer mechanisms, binary population synthesis models predict a gap in the orbital period distribution of post-mass-transfer binaries in the range 1–10 yr, and circular orbits below this gap (Pols et al. 2003; Izzard et al. 2010; Nie et al. 2012). However, the observed period distributions of the various types of post-AGB binaries discussed above show no sign of such a gap; in fact their periods appear to concentrate in the predicted gap.

This discrepancy points to shortcomings in our understanding of the evolution of binary systems containing red giants and AGB stars. Some of these shortcomings are related to the inadequacy of the BHL accretion scenario to describe the wind interaction process. The BHL description is a good approximation for wind accretion only when the outflow is fast compared to the orbital velocity. These conditions are not met for AGB winds, which are slow compared to the typical orbital velocities and have substantial density gradients. Furthermore, the slow wind velocities allow the escaping gas to interact with the binary and can produce non-isotropic outflows that carry away more angular momentum than in the isotropic case. The complicated interaction between the slow wind and the binary cannot easily be described by an analytical model, and hydrodynamical simulations are needed to study this process.

Several such hydrodynamical studies of the interaction between a red giant undergoing mass transfer via its wind with a companion star have been performed over the last two decades, covering a wide range of phenomena (Theuns & Jorissen 1993; Theuns et al. 1996; Mastrodemos & Morris 1998; Nagae et al. 2004; Mohamed & Podsiadlowski 2007; de Val-Borro et al. 2009, 2017; Huarte-Espinosa et al. 2013; Liu et al. 2017; Chen et al. 2017). Theuns & Jorissen (1993) and Theuns et al. (1996) were the first to carry out 3D smoothed-particle hydrodynamics (SPH) simulations of wind mass transfer, in which they showed that the morphology of the flow differs from BHL expectations. They also showed that the equation of state (EoS) used is important to determine the formation of an accretion disk around the companion. Mohamed & Podsiadlowski (2007) and Mohamed (2010) studied the problem of AGB mass transfer in very wide binaries using SPH simulations, in which they took into account the acceleration of the wind by radiation pressure on dust particles. They proposed a new mode of mass transfer, intermediate between RLOF and wind accretion, which they called wind Roche-lobe overflow (WRLOF). They showed that if the dust is formed outside the Roche lobe of the AGB star, wind material can be transferred efficiently through the inner Langrangian point, leading to much higher accretion rates than predicted by the BHL prescription.

The studies described above have helped to understand the physical processes involved in AGB wind interaction. However, only a few studies have focused on the implications for the orbital evolution of binary systems undergoing wind interaction. Given the discrepancy between the timescales that are accessible in simulations and binary evolution timescales, the computation of angular-momentum loss is an important way to predict the change in the orbit. In the context of X-ray binaries, Brookshaw & Tavani (1993) developed the first numerical simulations to study the loss of angular momentum and its influence on the binary orbit, by means of ballistic calculations to model wind particles ejected by the mass-losing star. Similar calculations are presented in Hachisu et al. (1999) in the context of wide symbiotic binaries with a mass-losing red giant as potential progenitors of Type Ia supernovae. Both studies find that the specific angular-momentum loss depends on the orbital parameters and the wind ejection velocity. Although the ballistic treatment and its neglect of any wind acceleration mechanism appears to be adequate for fast winds, for low-velocity mass outflows a hydrodynamical treatment is needed. Jahanara et al. (2005) used a three-dimensional grid-based code to model binary systems with one star undergoing wind mass loss. They study the amount of angular momentum removed as a function of the wind speed at the Roche-lobe surface for different assumed mass-loss mechanisms. Their calculations show that for low wind velocities the specific angular-momentum loss is large (although less so than implied by the ballistic calculations mentioned above), whereas for high velocities the specific angular-momentum loss decreases to the value expected for a non-interacting, isotropic wind. Recently, Chen et al. (2018) performed grid-based hydrodynamics simulations which also include radiative transfer in order to determine the orbital evolution of binary systems interacting via AGB wind mass transfer. For relatively short orbital periods they find that mass transfer will occur via WRLOF, leading to a shrinking of the orbit and to a possible merger of the components, whereas for wider separations the mass transfer process resembles the BHL scenario and the orbit tends to expand.

In this paper we try to bridge the previously discussed gap in the orbital period diagram by computing the angular-momentum loss and its effect on the orbit of a binary containing a masslosing AGB star. To do so, we perform smoothed-particle hydrodynamical simulations including cooling of the gas to model the AGB wind. We present our results as a function of the ratio of the terminal wind velocity to the orbital velocity (v_∞/v_orb). In Sect. 2, we briefly review the equations governing the orbital evolution of a binary system. In Sect. 3 we describe the model we used and the numerical set-up for the cases we studied, and in Sect. 4 we show the results. In Sect. 5 we discuss the implications of our findings for the orbital evolution of binary systems. Finally in Sect. 6 we conclude.

2. Angular momentum loss and mass accretion rate

The total orbital angular momentum of a binary system in a circular orbit is given by:$$J={\mu a}^2\mathrm{\Omega },$$(1)

where μ = M_a M_d/(M_a + M_d) is the reduced mass of the system, M_a and M_d are the masses of the accretor and the donor star respectively, a the orbital separation of the system, and Ω the angular velocity of the binary.

If the donor star is losing mass at a rate Ṁ_d < 0 and mass transfer is non-conservative, the companion star will accrete a fraction β of the material and the rest will be lost, carrying away angular momentum from the system. The change in orbital angular momentum can be parametrised as:$$\dot{J}={\eta a}^2\mathrm{\Omega }{\dot{M}}_{\text{bin}},$$(2)

where Ṁ_bin = (1−β) Ṁ_d is the mass-loss rate from the system and η is the specific angular momentum lost in units of J/μ. Hence, the change in orbital separation for non-conservative mass transfer will be:$$\frac{\dot{a}}{a}=-2\frac{\dot{M_{\mathrm{d}}}}{M_{\mathrm{d}}}[1-\mathrm{\beta }q-\eta (1-\mathrm{\beta }\left)\right(1+q)-(1-\mathrm{\beta }\left)\frac{q}{2(1+q)}\right],$$(3)

where q = M_d/M_a is the mass ratio of the stars. This equation can be solved analytically only for a few limiting cases. In the Jeans or fast wind mode, mass is assumed to leave the donor star in the form of fast and spherically symmetric wind. Since the speed of the wind is much larger than the orbital velocity of the system, the wind does not interact with the companion, escaping and taking away the specific orbital angular momentum of the donor $({\eta }_{\text{iso}}=M_{\mathrm{a}}^2/M_{\text{bin}}^2)$. In the case that β = 0, the change in the orbit is ȧ/a = − Ṁ_d/M_bin. This mass transfer mode leads to a widening of the orbit. The mass transfer efficiency β is often described in terms of the BHL analytical model, which can be used as a reference with which to compare the simulation results. In the framework of a binary system, the BHL accretion rate is given by:${\dot{M}}_{\text{BHL}}={\alpha }_{\text{BHL}}\pi R_{\text{BHL}}^2{v}_{\text{rel}}\mathrm{\rho },$(4)

for a high-velocity wind and assuming a supersonic flow. Here, $R_{\text{BHL}}=2G{M}_a/v_{\text{rel}}^2$ is the BHL accretion radius, $v_{\text{rel}}^2=v_{\mathrm{w}}^2+v_{\text{orb}}^2$ is the relative wind velocity seen by the accretor (Theuns et al. 1996), ρ is the density at the position of the companion and α_BHL is the efficiency parameter, of order unity, that physically represents the location of the stagnation point in units of the accretion radius (Boffin et al. 2015). Using that ρ = Ṁ_d/4πa²v_w for a steady-state spherical wind and $v_{\text{orb}}^2=G(M_{\mathrm{a}}+M_{\mathrm{d}})/a,$, we can write the accretion efficiency in the BHL approximation as¹:$${\mathrm{\beta }}_{\text{BHL}}=\frac{{\mathrm{\alpha }}_{\text{BHL}}}{{(1+q)}^2}\frac{v_{\text{orb}}^4}{v_{\mathrm{w}}{(v_{\mathrm{w}}^2+v_{\text{orb}}^2)}^{3/2}}\mathrm{.}$$(5)

Theoretical considerations and numerical simulations of BHL accretion (e.g. see Edgar 2004; Matsuda et al. 2015) indicate that for a uniform flow, α_BHL has a value of about 0.8. In applications to wind accretion Eq. (5) is sometimes divided by a factor of two (e.g. Boffin & Jorissen 1988; Hurley et al. 2002; Abate et al. 2013). The corresponding value of α_BHL is then higher by factor of two (e.g. Abate et al. 2015b, use a standard value of α_BHL = 1.5 in their models).

In realistic situations, and especially in the case of AGB wind mass transfer in binaries, the values of β and η cannot be expressed in analytical form and have to be derived from, for example, hydrodynamical simulations. This is the purpose of this paper. We use the fast wind mode and the BHL accretion efficiency as a reference with which to compare our simulation results.

3. Method

3.1. Stellar wind

The mechanism driving the slow winds of AGB stars is not well understood. The pulsation-enhanced dust-driven outflow scenario describes it in two stages (Höfner 2015): in the first stage, pulsations and/or convection cause shock waves which send stellar matter on near-ballistic trajectories, pushing dust-free gas up to a maximum height of a few stellar radii. At this distance, the temperature has dropped enough (~1500 K) to allow condensation of gas into dust. In the second stage, the dust-gas mixture is accelerated beyond the escape velocity due to radiation pressure onto the dust. Simulating this process involves several physical mechanisms for which hydrodynamical and radiative transfer codes have to be invoked. This goes beyond the scope of this work. Instead we use a simple approximation to model the velocity profile of the AGB wind. In order to do so, we use the STELLAR_WIND.PY (van der Helm et al. 2018) routine of the AMUSE² framework (Portegies Zwart et al. 2009, 2013; Pelupessy et al. 2013). Below we briefly explain how this routine works, for a detailed description we refer the reader to van der Helm et al. (2018).

 STELLAR_WIND.PY has three different modes to simulate stellar wind. The mode which best approximates wind coming from an AGB star is the accelerating wind mode and it works in the following way. The user provides the stellar parameters, such as mass, effective temperature, radius and mass loss rate, which can be derived from one of the stellar evolution codes currently in AMUSE. The initial and terminal velocities of the wind, v_init and v_∞, are user input parameters. SPH particles are injected into a shell that extends from the radius of the star to a maximum radius defined by the position that previously released particles have reached after a timestep Δt. The positions at which these particles are injected correspond to those derived from the density profile of the wind ρ(r). Thus, it is assumed that the velocity and density profiles of the wind, which are related by the mass conservation equation, are known. The wind acceleration function a_w(r) is obtained from the equation of motion:$$v\frac{\mathrm{d}v}{\mathrm{d}r}=a_{\mathrm{w}}\left(r\right)-\frac{G{M}_{\star }}{r^2}+c^2(\frac{2}{r}+\frac{1}{v}\frac{\mathrm{d}v}{\mathrm{d}r})$$(6)

where v = v(r) is the predefined velocity profile for the wind, and c is the adiabatic sound speed in the wind. The second term on the right-hand side describes the deceleration by the gravity of the star and the last term represents the gas-pressure acceleration. By computing a_w in this way and applying this acceleration to the gas particles, we guarantee that the gas follows the predefined density profile. If v_init = v_∞, then the term on the left hand side will be zero and the acceleration of the wind balances the gravity of the star and the gas pressure gradient. We note that Eq. (6) assumes an adiabatic equation of state (EoS) and ignores cooling of the gas. We have verified that ignoring the cooling term in Eq. (6) does not change the physical properties of the wind such as density and velocity. The module STELLAR_WIND.PY offers a variety of functions that resemble different types of accelerating winds. Although the most appropriate prescription would be to choose a function with an accelerating region, for the present work we have chosen a constant velocity profile, i.e. v_init = v_∞, for two reasons: to reduce the computational time and to be able to compare our results to other work where similar assumptions have been made.

3.2. Cooling of the gas

The internal energy change of the gas plays an important role when modelling the interaction of a star undergoing mass loss with a companion star. Theuns & Jorissen (1993) found that the formation of an accretion disk around the companion star depends on the EoS. In their study they modelled two cases, one in which the EoS was adiabatic (γ = 5/3) and one in which the EoS was isothermal (γ = 1). In the first case, no accretion disk was formed: when material gets close to the companion star, it is compressed by the gravity of the accretor, increasing the temperature and enhancing the pressure, expanding it in the vertical direction. In the isothermal case the pressure does not increase too much, permitting material to stay confined in an accretion disk. Although illustrative, these cases are not realistic and a more physical prescription for cooling of the gas is needed.

Proper modelling of cooling of the gas requires the implementation of radiative transfer codes. For simplicity, in this study, we use a modified approximation for modelling the change in temperature in the atmospheres of Mira-like stars based on Bowen (1988). In his work, the cooling rate Q̇ is given by:$$\dot{Q}=\frac{3k}{2\mu {\mathrm{m}}_u}\frac{(T-T_{\text{eq}})\mathrm{\rho }}{C}+{\dot{Q}}_{\text{rad}}\mathrm{.}$$(7)

The first term in Eq. (7) assumes that cooling comes from gas radiating away its thermal energy trying to reach the equilibrium temperature given by the Eddington approximation for a grey spherical atmosphere (Chandrasekhar 1934):$$T_{\text{eq}}^4=\frac{1}{2}{T}_{\text{eff}}^4\left\{\right[1-{(1-\frac{R_{\star }^2}{r^2})}^{1/2}]+\frac{3}{2}{\displaystyle {\int }_r^{\infty }}\frac{R_{\star }^2}{r^2}({\kappa }_{\mathrm{g}}+{\kappa }_{\mathrm{d}})\mathrm{\rho d}r\},$$(8)

where T_eff is the effective temperature of the star, κ_g and κ_d the opacities of the gas and dust respectively, R_⋆ is the radius of the star, r the distance from the star and ρ the density of the gas. The first term, $W=\frac{1}{2}[1-{(1-R_{\star }^2/r^2)}^{1/2}],$, corresponds to the geometrical dilution factor, and the term in the integral plays the role of the optical depth in plane-parallel geometry. The constant C in Eq. (7) is a parameter reflecting the radiative equilibrium timescale. Following Bowen (1988), we adopt a value of C = 10⁻⁵ g s cm⁻³.

The second term in Eq. (7) corresponds to radiation losses for temperatures above ~7000 K, where the excitation of the n = 2 level of neutral hydrogen is mainly responsible for the energy loss (Spitzer 1978). In this work however, we use an updated cooling rate prescription for high temperatures (log T ≥ 3.8 K) by Schure et al. (2009). The advantage of this prescription is that contributions to the cooling rate come not only from neutral hydrogen but individual elements are taken into account according to the abundance requirements. We use the abundances for solar metallicity given by Anders & Grevesse (1989). Then, the second term in Eq. (7) becomes:$${\dot{Q}}_{\text{rad}}=\frac{{\mathrm{\Lambda }}_{\text{hd}}{n}_{\mathrm{H}}^2}{\mathrm{\rho }}$$(9)

with Λ_hd interpolated between the values given in Table 2 of Schure et al. (2009) and n_H = Xρ/m_H.

For low temperatures the first term of Eq. (7) is the dominant term. If this term is not taken into account the internal energy of the particles reaches very low values leading to unphysical temperatures. We note that in Eq. (8), the optical depth at large distances from the star will be very small regardless of the opacity values, thus for distant regions the equilibrium temperature will be mainly determined by the geometrical dilution factor W. For this reason, and since the opacities are only used to calculate the equilibrium temperature at given radius, the opacities in Eq. (8) are taken constant during the simulations with values of κ_g = 2 × 10⁻⁴ cm² g⁻¹ and κ_d = 5 cm² g⁻¹ (Bowen 1988).

3.3. Computational method

To calculate self-consistently the 3D gas dynamics of the wind in the evolving potential of the binary, we used the SPH code FI (Hernquist & Katz 1989; Gerritsen & Icke 1997; Pelupessy et al. 2004) with artificial viscosity parameters as shown in Table 1. The AGB star was modelled as a point mass, whereas the companion star was modelled as a sink particle with radius corresponding to a fraction of its Roche Lobe radius R_L,a. In order to model the dynamics of the stars, we coupled the SPH code with the N-body code HUAYNO (Pelupessy et al. 2012) using the BRIDGE module (Fujii et al. 2007) in AMUSE. The N-body code is used to evolve the stellar orbits, and to determine the time-dependent gravitational potential in which the gas dynamics is calculated. Since the interest in this work is to study the evolution of the gas under the influence of the binary system, the BRIDGE module allows the gas to feel the gravitational field of the stars, however the stars do not feel the gravitational field of the gas; also selfgravity of the gas is neglected. This is a fair assumption given that the total mass of the gas particles in the simulation (4 × 10⁻⁵M_⊙) is very small compared to the binary mass (4.5 M_⊙). This assumption does not have a significant effect on the evolution of the gas, as was verified in a test simulation.

The stellar parameters of our simulated systems are chosen to match those presented in Theuns & Jorissen (1993). This will allow a direct comparison between our results and their work. The mass of the AGB star is M_d = 3 M_⊙ and that of the low-mass companion is M_a = 1.5 M_⊙. The radius of the primary star is R_d = 200 R_⊙ and for the secondary we assume a sink radius equal to R_a = 0.1R_L,a. The orbit is circular and the separation of the stars is a = 3 AU in our test simulations. The mass-loss rate Ṁ_d = 10⁻⁶M_⊙ yr⁻¹ is constant during the simulation; the terminal velocity of the wind is v_∞ = 15 km s⁻¹ and the gas is assumed to be monoatomic γ = 5/3 with constant mean molecular weight μ = 1.29 corresponding to an atomic gas with solar chemical composition (X = 0.707, Y = 0.274, Z = 0.019; Anders & Grevesse 1989). The effective temperature of the AGB star, which is also the initial temperature given to the gas particles, is T_eff = 2430 K for all models except the isothermal test case, for which T_eff = 4050 K in order to compare to Theuns et al. (1996). We assume the AGB star to be non-rotating as a result of its history of expansion after the main sequence, combined with angular momentum loss resulting from earlier mass loss on the RGB and AGB. Thus, we ignore the possibility of subsequent spin-up by tidal interaction (see Sect. 5.1.4 for a discussion). The stars are placed in such a way that the centre of mass is located at the origin of the system of reference. The setup for the other simulations is the same except for the values of the separations and the terminal velocities of the wind. We also ran a few simulations with a smaller sink radius (see Table 2). Each simulation was run for eight orbital periods.

3.4. Testing the method

In this subsection, we compare our simulations (models T1, T2, T3 in Table 2) with those of Theuns & Jorissen (1993). The top panel of Fig. 1 shows the velocity field of the gas in the y = 0 plane, perpendicular to the orbital plane, for models T1, T2 and T3 after 7.5 orbits. The colourmap shows the temperature of the gas in the y = 0 plane. The donor star is located at x = 1 AU, z = 0 AU and the companion star at x = −2 AU, z = 0 AU. Gas expands radially away from the primary star. Similar to Theuns & Jorissen (1993), in the adiabatic case (model T1) we observe a bow shock close to the companion star where high temperatures are reached. The high-temperature regions are located symmetrically above and below the star in this plane. The gravitational compression of the gas by the companion also enhances the temperature in this region, expanding the gas in the vertical direction and preventing it from creating an accretion disk. On the other hand, when an isothermal equation of state is used (model T2), an accretion disk is formed around the companion. In addition, two spiral arms are observed around both systems similar to those observed by Theuns & Jorissen (1993). These features can be seen in the bottom panel of Fig. 1, where density in the orbital plane (z = 0) is shown. Model T3 shows that using an EoS that includes cooling of the gas (Eq. (7)) gives a similar outflow structure as in the isothermal case, including the formation of an accretion disk. However, in this model there is an increase in temperature where the spiral shocks are formed. These effects are also found in our science simulations, which we describe in detail in Sect. 4.2.

Fig. 1. Top: flow structure in the y = 0 plane for different EoS after 7.5 orbits. The colourmap shows the temperature at y = 0. The white arrow on top of the velocity field plots corresponds to the magnitude of a 50 km s−1 velocity vector. In the three images we can observe the wind leaving the donor star radially at x = 1 AU. The accretor star is located at x = −2 AU. For model T1, the temperatures reached in the region near the companion star are very high preventing material to settle into an accretion disk. Bottom: gas density in the orbital plane (z = 0) is shown. The accretor is located at x = −2, y = 0, and the AGB star at x = 1, y = 0. In the middle panel we see that for the isothermal EoS (T2) gas in the vicinity of the accretor settles into an accretion disk. When cooling is included (T3), an accretion disk around the companion is also formed.

3.5. Mass accretion rate and angular momentum loss

We measure the mass accretion rate by adding the masses m_g of the gas particles that entered the sink during each timestep δt in the simulation. Due to the discreteness of the SPH particles, the resulting accretion rates will be subject to shot noise. This can be seen as fluctuations in the mass accreted as a function of time. In order to suppress statistical errors in the results, we average the mass accretion rate over longer intervals of time.

Measuring the angular momentum loss rate from the system is more complicated. An advantage of SPH codes is that they conserve angular momentum extremely well compared to Eulerian schemes (Price 2012). One question that arises is at which distance an SPH particle is no longer influenced by the gravitational potential of the binary system. In order to determine this, we construct spherical shells at various radii r_S from the centre of mass of the binary. In the same way as for the accretion rate, we compute the net mass-loss rate through such a surface. Each particle i carries angular momentum J_i , of which we add up the components J_{z ,i} along the z-axis, perpendicular to the orbital plane. We compute the specific angular momentum of the outflowing mass as follows,$${\dot{J}}_{\text{loss}}=\frac{{\mathrm{\Sigma }}_i^N{J}_{z,i}}{\mathrm{\Delta }N{m}_{\mathrm{g}}},$$(10)

with ΔN the number of particles that crossed the shell during a timestep Δt. Similar to the accretion rate, we average the mass and angular-momentum loss over longer intervals of time to suppress statistical fluctuations. We also verified that, once the simulation has reached a quasi-steady state, the orthogonal (x and y) components of the angular momentum-loss are negligible.

In Fig. 2 we show the resulting average mass and angular momentum-loss per orbit in our standard model V15a5 (see Table 2) as a function of time for different shell radii r_S. After several orbits, both quantities become approximately constant in time and almost independent of the chosen shell radius. This steady state is reached later for larger radii, corresponding to the longer travel time of the wind particles to reach this distance from the donor star. The near-constancy of the angularmomentum loss rate as a function of r_S indicates that the torque transferring angular momentum from the orbit to the escaping gas operates inside a few times the orbital separation a. Beyond r_S ≈ 10 AU (≈2a) the transfer of angular momentum to the gas is essentially complete, and angular momentum is simply advected outwards with the flow. We therefore adopt the mass and angular momentum loss values measured at r_S = 3a, which reach a steady state after about 4 orbits in model V15a5. The time taken to reach the quasisteady state depends on the terminal velocity of the wind, being shorter for higher velocity.

Fig. 2. Top: average mass loss per orbit for the standard model V15a5 measured at different radii from the centre of mass of the binary as indicated by different colours. Bottom: corresponding angular-momentum loss. The black thick curve for rS = 15 AU corresponds to the radius at which we measure J̇bin.

Once a simulation has reached a quasi-steady state, we compute the mass accretion efficiency β = Ṁ_a/ Ṁ_d from the average mass accretion rate over the remaining N orbits in the simulation. Similarly, we average the angular momentum loss (Eq. (10)) over the last N orbits to compute the value of η, i.e. the specific angular momentum of the lost mass in units of a²Ω.

4. Results

We performed 13 different simulations, listed in Table 2, for which we discuss the results in the following subsections. The first three simulations (T1–T3) are performed to test the EoS (see Sect. 3.4). The following three simulations (R1–R3) correspond to the test of different SPH resolutions, compared to simulation V15a5 which we consider the standard model. For this model we chose an orbital separation of 5 AU and a terminal velocity of the wind of 15 km s⁻¹. V10a5, V30a5 and V150a5 correspond to simulations with the same orbital separation but different velocities of the wind, either higher or lower than for the standard case. V10a5s2, V15a5s2 and V30a5s2 correspond to test simulations in which we study the effect of assuming a smaller sink radius (R_a = 0.06R_L,a) on the mass-accretion rate and radius of the accretion disk. It should be noted that simulation V10a5s2 was performed at a lower resolution than the others. Finally in V19a3 we change the orbital separation to 3 AU, keeping the same ratio of the velocity of the wind to the orbital velocity of the system as in V15a5.

4.1. Convergence test

Choosing the best resolution for SPH simulations is not simple and depends on the physical process under investigation. To determine the optimal resolution within a reasonable computational time, we performed four simulations of our standard model in which we vary the mass of the SPH particles m_g (R1, R2, R3 and V15a5 in Table 2). We note that a different m_g also implies a change in the average smoothing length h̄ of the SPH particles. Table 2 shows the value of h̄ within a radius a of the centre of mass of the binary, and the 95% confidence interval of smoothing-length values within this radius. The total number of particles n generated during one simulation is given by n = Ṁ_dt_end/m_g. Thus, the number of particles generated for the lowest resolution (largest m_g) is n = 5.6 × 10⁵ and n = 3.3 × 10⁶ for the highest resolution (smallest m_g). This last simulation ran only for six orbital periods.

Since the interest of this study is to obtain numerical estimates for the average mass-accretion efficiency and average angular-momentum loss during the simulation, these were the quantities we used to check for convergence. Figure 3 shows the values obtained for these quantities (as explained in Sect. 3.5) as a function of the mass of the SPH particles m_g. The error bars in this figure correspond to five times the standard error of the mean σ_m, which was estimated by means of ${\sigma }_{\mathrm{m}}^2=1/N{\mathrm{\Sigma }}_{i=1}^N{\sigma }_{\mathrm{m},{\text{orb}}_i}^2,$, where σ_m,orbi is the standard error of the mean per orbit³ and N the number of orbits during which the quasi-steady state has been reached. Since the number of timesteps over which the average is taken is large, the standard error turns out to be very small. This is the reason why we chose to plot five times the value of σ_m in Fig. 3. Because the simulation with the highest resolution (smallest m_g) only ran for six orbital periods, and the system reaches the quasi-steady state after four orbits, the error in the mean is larger than for the other simulations.

Fig. 3. Mean values of the β and η parameters, relating to mass and angular momentum loss, as a function of the mass of the SPH particles used to test different resolutions for the simulations. The error bars correspond to the overall standard error of the mean σm multiplied by a factor of 5 for better appreciation.

One thing that can be seen from Fig. 3 is that, even though by only a small amount, the accretion efficiency increases with increasing resolution, contrary to what Theuns et al. (1996) observed. However, we note that in their work the accretion rate was computed by setting an accretion radius proportional to the resolution of the SPH particles (in terms of the smoothing length h), whereas in our case the radius of the sink has the same constant value for all our models. On the other hand, the angular momentum per unit mass lost by the system is approximately independent of the resolution.

We also find that the time variability in the mass-accretion rate (described in Sect. 4.2.3) increases with increasing resolution. It is not clear why this takes place, but we should bear in mind that this variability could also be an artefact of the numerical method. As will be discussed in Sect. 4.2.3 there appears to be a correlation between the accretion disk mass and the mass-accretion rate. Given that the results of SPH simulations for accretion disks depend on the formulation of artificial viscosity, we performed a test simulation with the lowest resolution (m_g = 8 × 10⁻¹¹) adopting artificial viscosity parameters in FI equal to those used by Wijnen et al. (2016) for their protoplanetary disk models (α_SPH = 0.1 and β_SPH = 1). We find almost the same value for the average mass-accretion rate (β = 0.223 ± 0.002, versus β = 0.221 ± 0.002 for simulation R3). With this alternative prescription of artificial viscosity we find similar variability in the mass-accretion rate.

In order to guarantee reliable results within a reasonable computational time, the chosen resolution for the other simulations is the one in which the mass of the individual SPH particles is m_g = 2 × 10⁻¹¹M_⊙. For this simulation, the value of β differs by 9% compared to the value obtained for the simulation with highest resolution (R1). This is indicative of the numerical accuracy of the accretion rates obtained in our simulation. By contrast, the angular-momentum loss rates (η values) we measure differ by <0.4% between simulations at different resolution.

4.2. Summary of model results

4.2.1. Morphology of the outflow

Since the velocity of the wind is of the same order as the orbital velocity of the system, Coriolis effects play an important role in shaping the outflow. Figure 4 shows the density of the gas in the orbital plane for different simulations after 7.5 orbits. The orbital motion is anticlockwise in this view. We observe that in the simulations with relatively low wind velocities, v_∞/v_orb < 1 (V10a5, V15a5, and V19a3), two spiral arms are formed around the system. These spiral arms delimit the accretion wake of the companion star. Similar to Theuns & Jorissen (1993), we find that the inner spiral shock, wrapped around the donor star, is formed by stellar-wind material leaving the AGB star colliding with gas moving away from the companion star with v_x > 0. Near the companion star a bow shock is created by the wind coming from the AGB star as it approaches the accretor. The temperature of the gas in this region increases and material is deflected behind the star forming the outer spiral arm. The Mach number of the flow near the accretor has values between 2 and 6. Because of its proximity to the donor star, where the gas density is high, the inner spiral arm has higher densities than the outer spiral arm. In these systems we also observe an accretion disk around the companion (Sect. 4.2.2). The disk is asymmetric and is rotating anticlockwise.

Fig. 4. Gas density in the orbital plane (z = 0) for V10a5, V10a5s2, V15a5, V19a3, V30a5 and V150a5 after 7.5 orbits.

Figure 4 shows that for model V30a5 the density of the gas in the accretion wake is lower than in models with v_∞ < v_orb (we note that a different colour scale is used in this panel). The density of the accetion wake is high close to the companion star, but it decreases along the stream. Two spiral arms delineate the edge of the accretion wake, which converge into what appears to be a single arm close to the companion star. No accretion disk is observed in this model. For model V150a5 a single spiral arm is distinguished. The density along this stream is very low, because for simulation V150a5 the gas density is much lower than for the other simulations and only a small fraction of the stellar wind is focused into the accretion wake. Also notice that given the low density of the gas, the effective resolution of this model is lower. Given the large velocity of the wind compared to the orbital velocity of the system, the wind escapes the system without being affected by the presence of the companion star and remains spherically symmetric.

If we compare simulations V15a5 and V30a5 we see that the angle formed by the inner high-density spiral arm and the axis of the binary changes as a function of the velocity of the wind. For the low-velocity model, this angle is very sharp and as the wind velocity increases the angle of the stream becomes more oblique. This is expected given that the angle of the accretion wake depends on the apparent wind direction seen by the accretor in its orbit, such that tan θ = v_orb/v_wind.

4.2.2. Accretion disk

An accretion disk is formed around the accretor in the simulations with v_∞/v_orb < 1 (V10a5, V10a5s2 V15a5, V15as2 and V19a3). Figure 5 shows the density of the gas in the orbital plane for the accretion disks formed around the secondaries for simulations V10a5, V15a5 and V19a3. One common feature is that the accretion disks are not symmetric. This asymmetry is associated with wind coming from the donor star colliding with material already present in the disk. Another feature in the flow pattern of the simulations is formed by two streams of gas feeding the disk: as seen from the velocity field in the corotating frame (not shown), one stream arises from material moving along the inner spiral arm with v_x , v_y < 0 and the second one comes from material moving along the outer spiral arm with v_x, v_y > 0.

Fig. 5. Left panel: density in units of g cm−3 in the orbital plane of the region centred on the accreting star for three models with v∞/vorb < 1, zooming in on the accretion disks. The size of the box corresponds to the Roche-lobe diameter of the companion star. For models V15a5 and V19a3, the situation after 7.5 orbits is shown, whereas for model V10a5 we show the accretion disk before it disappears, but in the same orbital phase. Right panel: mass of the accretion disk, i.e. the mass contained in a sphere of radius Rdisk (see Table 3) as a function of time.

Even though the accretion disk is not symmetric, we approximate it as such in order to estimate its size and mass. We construct spheres of different radii centred at the position of the companion star and measure the mass within each sphere as a function of time. The size of the disk is taken as the radius of that sphere outside which the mass is approximately constant as a function of radius. The disk mass is taken as the total gas mass within that sphere. This approach is reasonable because the gas mass inside the spheres is dominated by the mass of the disk. Table 3 shows the resulting values of the accretion disk radius and mass for the simulations in which a disk was present after eight orbital periods, except for model V10a5 where the given radius and mass correspond to a time of two orbits.

In the right-hand panels of Fig. 5, we show the mass within the disk as a function of time for the same models. For model V15a5 (middle right), it is clear that the disk mass is variable over time. It is not obvious where this variability comes from, but we note that the outer disk radius is not much larger than the sink radius (34.7 R_⊙) which determines the inner edge. To investigate this behaviour we performed a simulation with the same parameters but a smaller sink radius (20.8 R_⊙, V15a5s2). In this simulation the disk reaches a larger radius and gradually grows in mass, up to a value at the end of the simulation about 9 times larger than in V15a5 (see Table 3). This can be understood as the smaller sink allows gas in the disk to spiral in to closer orbits around the secondary by viscous effects, transferring angular momentum to the outer regions of the disk which thus grows in size. Such behaviour is physically expected, but simulations with even smaller sink radii and extending over a longer time are required to investigate the mass and radius of the disk and their possible variability.

For model V10a5 (top right panel of Fig. 5), we observe a highly dynamic accretion disk, which disappears after two orbital periods. The moment when the disk vanishes can be seen as a steep decrease in the mass of the accretion disk. The reason for this disappearance is that after two orbital periods the size of the outer radius of the disk becomes smaller than the radius of the sink particle. When performing a similar simulation but with a smaller sink radius, V10a5s2, the accretion disk remains for the rest of the simulation and gradually increases in mass, similar to what was seen in simulation V15a5a2 discussed above. The size of the disk remains smaller than in the models with a wind velocity of 15 km s⁻¹. A similar pattern is found in model V19a3 (bottom panel of Fig. 5), in which the disk mass increases with time. A longer simulation will be needed to see if it converges towards a constant value.

For the models in which v_∞/v_orb > 1, no accretion disk is formed around the companion. Also in a test simulation with a smaller sink particle (V30a5s2) the accreted matter is swallowed by the sink particle without forming an accretion disk. We conclude that if an accretion disk were to form in this system it will be smaller than the assumed sink radius.

4.2.3. Mass-accretion rate

Figure 6 shows the average mass-accretion rate over intervals of one year for models V10a5, V15a5, V30a5 (solid lines), V10a5s2, V15a5s2 and V30a5s2 (dashed lines).

Fig. 6. Average accretion rate over one-year intervals for models V15a5 (top), V10a5 (middle) and V30a5 (bottom) is shown in blue. The dashed red lines correspond to the year-averaged mass accretion rate for the corresponding models with a smaller sink radius.

In models V15a5 (middle panel) and V19a3 (not shown), we find a lot of variation in the accretion rate, which clearly exceeds the noise associated with the finite number of particles accreted. This suggests a correlation between the mass of the accretion disk and the accretion rate. This can be observed if we compare the middle right panel of Fig. 5 and the middle panel of Fig. 6, where the peaks and troughs in the mass of the disk coincide in time with those observed in the mass-accretion rate. In the case of model V15a5s2, where the mass in the disk gradually increases with time, the accretion rate as a function of time also shows an increase (dashed line in the middle panel of Fig. 6). The average accretion rate for model V15a5s2 is a factor of 1.5 lower than for model V15a5 where the sink radius is larger. For model V10a5 (top panel of Fig. 6), we find that after the accretion disk becomes smaller than the sink and is swallowed, the variation in the accretion rate decreases. For model V10a5s2 with a smaller sink, although the accretion disk is present for the entire simulation, the variability in the accretion rate is lower than for model V15a5s2 (where the accretion disk is also present over the entire simulation). However, we should note that the resolution for model V10a5s2 is much lower. Similar to models V15a5 and V15a5s2, the massaccretion efficiency is lower for model V10a5s2 than for model V10a5. In the bottom panel of the same figure, we show the accretion rate for model V30a5. Apart from the Poisson noise associated with the finite resolution of the SPH simulations, the accretion rate is constant in time. For the same simulation but with a smaller sink we find a slightly lower accretion efficiency.

We find that the mass that is not accreted in models V10a5s2 and V15a5s2 is stored in the accretion disk. Figure 7 shows the net mass inflow rate into a shell of radius 0.4R_L,2 around the secondary for models V10a5, V15a5 and V30a5 (solid lines) and their corresponding models with a smaller sink radius (dashed lines). The chosen shell radius is larger than the size of the accretion disk in all the simulations, so that the inflow rate corresponds to the total mass-accumulation rate of the sink and the disk, Ṁ_0.4RL = Ṁ_sink + Ṁ_disk, since the amount of gas within this volume that is not in the disk is very small. In the steady state, the inflow rate we find in this way is insensitive to the chosen sink radius. For simulations V15a5 and V15a5s2, this combined accretion rate is seen to be very similar and also shows similar variability over time. A similar feature holds for simulation V10a5 and V10a5s2. Since it is likely that the mass in the disk will eventually reach the stellar surface, this may be a better measure of the steady-state accretion rate that is insensitive to the assumed sink radius. In the case of models V30a5, V30a5s2, we observe the same small discrepancy in the inflow rate at R = 0.4R_L,2 that we find in Fig. 6. The net influx of mass within the shell is the same in both simulations, reflecting the steadystate condition reached. However, when the sink radius is larger it captures material that may otherwise escape. Therefore, for simulations such as V30a5 in which no accretion disk is formed, our results provide only an upper limit to the amount of mass accreted.

Fig. 7. Net inflow rate into a sphere of radius 0.4RL,2 around the accretor star averaged over intervals of one year for the same models as in Fig. 6.

In order to determine the long term evolution of the orbit, one of the quantities of interest from these simulations is the average mass-transfer efficiency β. Given the results discussed above and shown in Figs. 6 and 7, we take the mass-transfer efficiency as β_0.4RL = Ṁ_0.4RL/Ṁ_d. In Table 3, we provide the mean values of this quantity for the science simulations, as well as the mean values for the material captured by the sink only, i.e. β_sink = Ṁ_sink/Ṁ_d. In Fig. 8, we show the corresponding values for the accretion efficiency β_0.4RL as a function of the ratio v_∞/v_orb for the science simulations. The dotted line corresponds to β_BHL in Eq. (5) with α_BHL = 1. Solid dots in the figure correspond to models V15a5–V150a5 and stars to models T3 and V19a3, both of which have a = 3 AU.

Fig. 8. Fraction of mass accreted by the companion as a function of v∞/vorb. The dotted line corresponds to BHL accretion rate with BHL = 1, and the dashed line to the accretion rate for a geometrical cross section of size Rsink. Dots correspond to models in which the orbital separation is 5 AU and stars to models in which a = 3 AU. The values shown are measured at a radius of 0.4 RL,2 (β0.4RL,2).

Model V150a5, in which we use a wind velocity ten times higher than the typical velocities of AGB winds, was set up to approximate the isotropic wind mode. Figure 8 shows that for this model, our numerical result for the accretion efficiency exceeds the expected value for BHL accretion by a factor of 2.5 (see also Table 3). However, we note that the BHL accretion radius for this simulation is smaller than the radius of the sink, which will result in a discrepancy with the BHL prediction for the accretion efficiency. The dashed line in the same figure corresponds to the accretion efficiency β_sink assuming the geometrical cross section of the sink instead of the BHL cross section. The difference between the simulation result and β_sink is reduced to a factor of 1.3.

Figure 8 shows several interesting features. As the ratio v_∞/v_orb increases, the value of β decreases, following a similar trend as expected from BHL accretion. However, for v_∞/v_orb < 1 we find values that exceed the β_BHL by up to a factor of 1.5–2.3. For the lowest wind velocity, the mass transfer efficiency is approximately 40%. Finally, for models V15a5 and V19a3, which have the same v_∞/v_orb ratio, the accretion efficiencies shown in Fig. 7 differ by a factor of 1.4.

4.2.4. Angular-momentum loss

The second quantity of interest resulting from our simulations is the value of the specific angular momentum of the material lost by the system. Figure 9 shows the values of this quantity in units of J/μ as a function of v_∞/v_orb for the different models. We see that the lower the velocity of the wind, the larger the specific angular momentum of the material lost. This is consistent with the expectation that when the velocity of the wind is small compared to the relative orbital velocity of the system, more in teraction between the gas and the stars occurs and more angular momentum can be transferred to the gas. This strong interaction is also seen in the high mass-accretion efficiency, as well as in the existence of a large accretion disk and in the structure of the spiral arms of the systems. The simulations with v_∞ < v_orb show spiral arms that are more tightly wound around the system, and with higher gas densities, than the simulations for higher wind velocities (see Sect. 4.2.1 and Fig. 4). These spiral arms correspond to the accretion wake of gas interacting with the secondary star. We interpret the loss of angular momentum as the result of a torque between the gas in this accretion wake and the binary system. The magnitude of this torque depends on the orientation of the wake and the density of the gas behind the companion star. When the accretion wake is approximately aligned with the binary axis and the density in the wake is relatively low, as in simulations with v_∞/v_orb > 1, the torque between the gas and the binary will be small. The torque exerted on the gas will increase as the wake misaligns with the binary axis, and as the gas density in the accretion wake becomes higher. Both these effects occur for low values of v_∞/v_orb, in particular in the inner spiral arm that has a high density and is oriented almost perpendicular to the binary axis (Fig. 4).Therefore, the transfer of angular momentum from the orbit to the outflowing gas increases strongly with decreasing v_∞/v_orb.Furthermore, it implies that most of the angularmomentum transfer occurs at short distances from the binary system, as is confirmed by the test discussed in Sect. 3.5 (see Fig. 2).

Fig. 9. Specific angular momentum of the mass lost from the system as a function of v∞/vorb. The lower the velocity of the wind compared to the orbital velocity, the higher the specific angular momentum loss. For equal velocity ratios, the specific angular momentum is very similar. The dashed line shows the value ηiso expected for an isotropic wind for a mass ratio q = 2.

Model V150a5 approximates the isotropic or fast wind mode fairly well. In this case, the ejected matter escapes from the system with very little interaction, only removing the specific angular momentum of the orbit of the donor star. For this simulation, the value we obtained for η is η_num = 0.112 ± 0.001, very close to the expected value for the mass ratio of the stars in our models, η_iso = 0.111 (dotted line in Fig. 9).

An interesting result of our models is that for the same value of v_∞/v_orb, but different orbital separations, the specific angular momentum of the mass lost is very similar. For models V15a5 η ≈ 0.52 and for V19a3, η ≈ 0.54. Model V10a5 has the largest loss of angular momentum among the high-resolution simulations with η = 0.64, although the low-resolution simulation T3 has an even large amount of angular-momentum loss, η = 0.68.

5. Discussion

5.1. Assumptions and simplifications

A complete study of binary systems interacting via AGB winds requires the simultaneous modelling of hydrodynamics, radiative transfer, dust formation and gravitational dynamics. In this study, several simplifications have been made.

5.1.1. Wind acceleration

Dust formation and radiative transfer are not included in our models (except for the simplified effective cooling model) although these processes play a major role in accelerating the gas away from the star. Instead of explicitly computing the gas acceleration, a constant velocity of the gas is assumed with values typical of the terminal velocities of AGB winds. This method is chosen to simplify the calculations and to guarantee that the velocity of the wind v_∞ at the location of the companion star has a predefined value as we want to study the effect of the velocity of the wind on the interaction of the gas with the stars. On the other hand, observations and detailed wind models both indicate that most of the acceleration of AGB winds occurs in the dustformation zone located at R_dust ~ 2 − 3R_⋆ (Höfner & Olofsson 2018). For the orbital parameters assumed in our models, R_dust is similar to or larger than the Roche-lobe radius of the donor star, R_L,1. As described in Mohamed (2010) and Abate et al. (2013), in this situation mass transfer can occur by WRLOF, i.e. by a dense flow through the inner Lagrangian point into the Roche lobe of the companion. However, given that the wind particles in our simulations are injected with predefined terminal velocities from the radius of the donor star, no WRLOF is observed.

If the gradual acceleration of the wind is taken into account, the density and velocity structure of the gas inside the orbit of the companion will be different which may affect the outflow morphology in our low wind-velocity simulations. To some extent the effects of wind acceleration can be accounted for by assuming an increasing velocity profile in Eq. (6), mimicking the velocity profile of an AGB wind. The resulting acceleration in such models would still be spherically symmetric. In addition, the modified outflow in a binary can change the optical depth of the dust and thus affect the acceleration of the wind itself, making it non-isotropic. This can only be studied by explicitly modelling dust formation and radiative transfer in wind masstransfer simulations (e.g. see Mohamed & Podsiadlowski 2012; Chen et al. 2017) and is beyond the scope of the present study. For this reason, and although the simulations presented here give insight into the dependence of the angular-momentum loss and the accretion efficiency on the v_∞/v_orb ratio, our results should be interpreted with some caution.

5.1.2. Accretion

Given the small radius of a main-sequence or white-dwarf companion star, resolving it spatially would be very computationally demanding. For this reason the accretor was modelled as a sink particle with radius equal to a fraction of the size of its Roche lobe. In the simulations where an accretion disk is formed, the inner disk radius is limited by the sink radius. However, physically there is no reason to expect that the disk will not extend inwards to the stellar radius by viscous spreading, and any gas captured in the disk will be transferred towards the stellar surface. As we have shown in Sect. 4.2.3, the sum of the accretion rate onto the sink and the mass accumulation rate in the disk appears to be independent of the assumed sink radius, and should give a reliable measure of the mass transfer rate to the companion. On the other hand, we cannot be sure that all matter that approaches the surface of the star through the disk will be accreted, as some fraction of it may be expelled by a wind or jet formed in the inner part of the accretion disk. In simulations with higher wind velocities where we do not find an accretion disk, it is possible that a disk is still formed at a smaller radius within the sink. However, some of the mass captured within the sink may not end up in the disk and reach the surface of the companion. For these reasons, and although our results for the mass accretion efficiency with different prescriptions for the artificial viscosity and sizes of the sink particle seem to converge, the values provided should be taken as upper limits.

5.1.3. Angular-momentum transfer

In our simulations the gravitational field of the gas is neglected, i.e. the gas feels the gravitational potential of the stars but the stars do not feel the gravity of the gas. This does not prevent us from inferring the gravitational influence of the gas on the binary system, because conservation of angular momentum dictates that the angular momentum transferred to the gas, which we measure from the simulation, is taken out of the orbit. This is expressed by Eq. (3) in Sect. 2. We have performed test simulations in which the gravitational forces are symmetric, i.e. the stars do feel the gravity of the stars, and find that the results for the gas dynamics and angular-momentum loss are not affected. We also find (see Sect. 4.1) that the specific angular momentum of the outflowing gas is very well conserved once it travels beyond a few orbital radii. This is an advantage of using an SPH code for this work, because SPH codes are much better at conserving angular momentum than grid-based hydrodynamical codes. Therefore, we are confident that our results for the angular-momentum loss are robust and insensitive to the numerical assumptions in the simulations.

5.1.4. Rotation

The stars in our simulations are assumed to be non-rotating. As mentioned in Sect. 3, given the evolutionary history of expansion and mass loss of the AGB star, we expect the rotation rate of the donor to be negligible compared to the orbital frequency, as long as tidal interaction can be ignored. However, the tidal synchronisation timescale of the donor can be fairly short. Applying the equilibrium tide model as described in Hurley et al. (2002) and using the orbital and stellar parameters of the simulated systems (in particular, a stellar radius of 200 R_⊙) we find a synchronisation timescale of several times 10⁴ yr for an orbital separation of 5 AU, and about 10³ yr for 3 AU. This is much shorter than the expansion timescale of the AGB star, which is of the order of 10⁶ yr, indicating that the donor is likely to rotate synchronously with the orbit. As a consequence, an additional transfer of angular momentum to the gas will take place at the expense of the rotational angular momentum of the donor. This occurs at a rate:$${\dot{J}}_{\text{rot}}=\frac{2}{3}{R}_{\mathrm{d}}^2\mathrm{\Omega }{\dot{M}}_{\mathrm{d}},$$(11)

if we assume the wind decouples from the star at a spherical shell of radius R_d. In the absence of magnetic or hydrodynamic coupling between the outflowing wind and the star, R_d can be identified with the stellar radius. If the donor remains tidally locked to the orbit, this angular momentum is continually supplied to the gas from the orbit, leading to additional orbital angular-momentum loss compared to the non-rotating case. We choose not to include this in our simulations, because the effects of spin angular-momentum loss and tidal interaction can be taken into account separately, as is often done in binary population synthesis modelling (e.g. Hurley et al. 2002; Abate et al. 2015c). For the systems we have simulated, the rotational angular-momentum loss implied by Eq. (11) is only a small fraction of the orbital angular-momentum loss, so the effect on the orbit will be small. By ignoring rotation, however, we also neglect the possibility that the morphology of the outflow itself is modified, and thereby the way it interacts with the companion. However, if the wind decouples from the donor star at a radius much smaller than the orbital separation, the flow in the vicinity of the companion will not be strongly affected.

On the other hand, the companion star will gain not only mass but also angular momentum from the material accreted which will lead to spin-up of the star. Once the star is spun up to its critical rotation, no more accretion can take place (Packet 1981; Matrozis et al. 2017), imposing a constraint on the amount of accreted material. In our simulations we keep track of the angular momentum accreted by the sink particle. However, this is not representative of the true angular momentum added to the companion because we use a sink radius that is much larger than the expected stellar radius. In cases where an accretion disk is formed, the angular momentum of the accreted gas corresponds to a Kepler orbit at the sink radius, whereas in reality by the time the gas reaches the surface of the star it would have transferred angular momentum to the outer regions of the disk. Therefore our simulations do not allow us to study the effects of angularmomentum accretion on the secondary (Liu et al. 2017).

5.2. Comparison to other work

5.2.1. Angular-momentum loss

Our results for the angular-momentum loss are in approximate agreement with other work which uses different methods. Jahanara et al. (2005) performed 3D hydrodynamical grid simulations of a star undergoing mass loss and interacting with a companion star. They study the amount of angular-momentum loss as a function of the wind speed at the surface of the donor for various mass ratios and different assumed mass-loss mechanisms. The mechanism that best approximates the mass loss from an AGB star is the radiatively driven (RD) wind mechanism. This is roughly equivalent to our assumptions, although Jahanara et al. (2005) use an adiabatic EoS without cooling and they assume a much higher sound speed than in our simulations, which leads to substantial gas-pressure acceleration in their models at low wind velocities. Despite these differences, and even though the mass ratio in our simulations, q = 2, is different from the mass ratios assumed in Jahanara et al. (2005), we can make a rough comparison of our results with their Fig. 7, in which q = 1. In that figure v_R = v_w,RL/v_orb corresponds to the average wind velocity at the Roche-lobe surface of the donor, which can be compared to our v_∞/v_orb ratio. The parameter ℓ_w corresponds to the specific angular momentum of the material lost in units of J/μ and is equivalent to our , although in Jahanara et al. (2005) the donor is kept in corotation with the orbit, so that ℓ_w includes the spin angular momentum loss from the donor star. We compare these results to our simulations in Fig. 10, where we have subtracted the small amount of spin angular momentum-loss from the ℓ_w values using Eq. (11), so that they only represent the orbital angular momentum loss. Jahanara et al. (2005) find that the strongest angular-momentum loss, with ℓ_w ≈ 0.6, occurs for the lowest wind velocities, corresponding to v_R ≈ 0.7, which is comparable but slightly larger than our results (see Fig. 10). For increasing wind velocity the specific angular momentum decreases, converging to values equal to the isotropic case for q = 1, i.e. η = 0.25, which is consistent with our results for model V150a5.

Fig. 10. Comparison of our results (Fig. 9) to other work. The data points correspond to the specific orbital angular momentum lost from the system via winds obtained using different methods. The open squares correspond to the results of Chen et al. (2018), triangles show the results of Jahanara et al. (2005) for the radiatively driven wind mode. In the latter case we have interpreted the wind velocity at the Roche-lobe surface as corresponding to v∞. The results from both papers have been corrected for the spin angular-momentum loss (see text) to compare with our results. The colours correspond to different mass ratios q and the lines to the fit formula described by Eq. (13). The blue dotted line is the fit for q = 2, and the orange dashed line for q = 1.

Chen et al. (2017, 2018) also performed grid code simulations of binary systems interacting via AGB wind mass transfer in order to study the orbital evolution, but including pulsations of the AGB star as well as cooling, dust formation and radiative transfer. They modelled systems consisting of a primary star of 1 M_⊙ with a terminal wind velocity v_∞ ≈ 15 km s⁻¹ and orbital separations between 3 and 10 AU. Most of their models have a secondary mass of 0.5 M_⊙, i.e. the same mass ratio as in our simulations. Their models with a > 6 AU display a similar morphology to our simulations with v_∞/v_orb ≳ 1, showing a spiral arm structure corresponding to the BHL accretion wake of the secondary. On the other hand, their models with smaller separation (and v_∞/v_orb ≲ 1) show a flow morphology resembling WRLOF and appear to be forming a circumbinary disk, which we do not find in our low-velocity simulations. These differences are likely due to the differences in modelling the wind acceleration process. Chen et al. (2018) express the angular momentum loss from the system in terms of a parameter γ, which is the specific angular momentum of the matter lost in units of J/M_bin. This is equivalent to our description in terms of η, using the transformation$$\eta =\frac{\mu }{M_{\text{bin}}}\gamma =\frac{q}{{(1+q)}^2}\gamma \mathrm{.}$$(12)

They find larger angular momentum loss from systems with smaller a, i.e. with smaller v_∞/v_orb, similar to our results. For their models with q = 2 and a ≥ 6 AU (which have 1.0 ≲ v_∞/v_orb ≲ 1.3) they find similar values of η to our model V30a5 with v_∞/v_orb = 1.06. However, for the model with the same mass ratio and a = 4 AU (v_∞/v_orb ≈ 0.8, intermediate between our models V15a5 and V30a5), they obtain a higher ≈ 0.6 than we find for V15a5 which has lower v_∞/v_orb.⁴ In making the comparison with our models we have to take into account that in the Chen et al. (2018) simulations the AGB star spin is synchronised with the orbit. The spin angular momentum transferred to the gas in their simulations is substantially larger than given by Eq. (11), which they ascribe to subsonic turbulence between the photosphere and the dust formation radius. When we subtract the contribution of the stellar spin (as given in their Table 4) from the total angular-momentum loss, the resulting η values (in particular for a = 4 AU) are in better quantitative agreement with our results (see open squares in Fig. 10).

5.2.2. Accretion efficiency

Several previous studies have investigated the mass-accretion efficiency during wind mass transfer. Theuns et al. (1996) and Liu et al. (2017) performed SPH simulations of a binary with exactly the same parameters as our test models (T1–T3), using an EoS without cooling. They find accretion efficiencies between 1% and 2.3% when using an adiabatic EoS, which is quite comparable to the low β value of our model T1 (and substantially less than expected from BHL accretion, see also Nagae et al. 2004). On the other hand, when they apply an isothermal EoS the accretion efficiency increases to 8%–11%, respectively. This is about a factor of three smaller than we find in our isothermal test model T2 and in model T3 that includes gas cooling explicitly. The reason for this difference is unclear, but we note that the algorithm for computing the accretion rate in these papers is different than the method we use in this work. Besides the discrepancies in the accretion efficiencies, these studies and ours both indicate that allowing for gas cooling results in substantially higher accretion efficiencies.

As discussed in Sect. 5.1.1, our simulations lie in a regime where WRLOF might be expected if wind acceleration were taken into account. In the WRLOF regime accretion onto the companion star occurs by a combination of material flowing through the inner Lagrangian point and gravitational focusing towards the companion star. Our simulations show the latter effect, but not the former. In the WRLOF simulations by Mohamed (2010; see also Abate et al. 2013) and Chen et al. (2017), which also include gas cooling, high accretion efficiencies of up to 40–50% are found. It is interesting that we find quite similar β values in our lowest wind-velocity models, even though there is no significant flow though the inner Lagrangian point in our models. This suggests that the high accretion efficiencies found in WRLOF simulations may be caused predominantly by gravitational focusing onto the companion star. However, given the uncertainties in reliably determining accretion rates as discussed in Sects. 4.2.3 and 5.1.2, one should be cautious in making such comparisons.

5.3. Implications for binary evolution

A very interesting result of our simulations is the fact that for models with different orbital separations but the same v_∞/v_orb, we find that a similar specific angular momentum is removed from the system. This suggests that the angular momentum loss, as expressed in the parameter η, may depend primarily on v_∞/v_orb and relatively little on other parameters of the system. The comparison made above with the results of Jahanara et al. (2005) and Chen et al. (2018) strengthens this tentative conclusion, and also suggests that η may be relatively insensitive to the binary mass ratio.

We find that the results of all three sets of simulations, after correcting η for the spin contribution as described in Sect. 5.2, are fairly well described by the following simple function:$$\eta ={\eta }_{\text{iso}}+\frac{1.2-{\eta }_{\text{iso}}}{1+{\left(2.2\frac{v_{\infty }}{v_{\text{orb}}}\right)}^3}\mathrm{.}$$(13)

This relation is shown in Fig. 10 for two values of the mass ratio, q = 2 (corresponding to our simulations and those of Chen et al. 2018) and q = 1 (corresponding to the results of Jahanara et al. 2005). The function is constructed to converge for very small v_∞/v_orb at η = 1.2, which is the maximum value found by Jahanara et al. (2005) in their low-velocity, “mechanically driven” wind models and appears to be independent of q. Most of the points are reasonably well fitted by this relation, the main exception being the Chen et al. (2018) result for v_∞/v_orb ≈ 0.8 which is well above the line. However, the sparse data from hydrodynamical simulations available so far do not warrant a fitting function with a larger number of adjustable parameters. We stress that we consider Eq. (13) to be very preliminary. A larger grid of simulations with different mass ratios, separations and mass-loss rates will be needed in order to investigate if in general the specific angular-momentum loss of the material can be written simply in terms of the velocity ratio.

Knowing the amount of angular momentum-loss and the fraction of mass accreted by the companion, i.e. the values of η and β, we can predict the evolution of the orbital separation by means of Eq. (3). Figure 11 shows the theoretical prediction for ȧ/a based on the values obtained for β and η for a mass ratio q = 2. For velocity ratios v_∞/v_orb > 1, we conclude that the orbit should widen, and for high terminal wind velocities relative to the orbital velocity the isotropic-wind regime is approached, with an accretion rate similar to the BHL approximation. On the other hand, when v_∞/v_orb < 1 the orbit will shrink, on a timescale similar to the mass-loss timescale of the AGB star. As an example, by integrating Eq. (3) over time and keeping constant the value of the mass-loss rate and the values of β and η obtained for model V15a5, we find that by the time the AGB star reaches a typical WD mass of 0.6 M_⊙, the orbital separation will shrink by a factor of ≈0.6, from 5 AU to ≈3 AU.

Fig. 11. Evolution of the orbital separation as a function of v∞/vorb for the systems simulated in this paper. The systems with ȧ/a < 0 will shrink whereas the ones with ȧ/a > 0 will expand.

This has important consequences for the evolution of AGB binaries and may help to explain the puzzling orbital periods of Ba stars, CEMP-s stars and other post-AGB binary systems as discussed in Sect. 1. Abate et al. (2015a,b) showed that in order to simultaneously explain the observed abundances and the short orbital periods of individual CEMP-s stars, efficient mass accretion and enhanced angular-momentum loss is needed compared to the predictions of BHL accretion and an isotropic stellar wind. In their models they used a value for the specific angular momentum of the escaping gas equal to two times the average specific angular momentum of the binary, i.e. γ = 2. Using Eq. (12) and a mass ratio q ≈ 2, this translates into values compatible with our simulations for v_∞/v_orb < 1. In their population synthesis models of Ba stars Izzard et al. (2010) found that including angular-momentum loss with the same value γ = 2 helps to reproduce the observed period distribution, while the isotropic-wind model produces periods that are too long. However, the impact of our results, in particular the relation between η and the wind velocity ratio, on the period distributions of post-AGB binaries has to be verified by population synthesis modelling.

Our results also suggest that the number of systems entering a CE phase will increase as a result of the shrinking orbits and high accretion efficiency during low-velocity wind interaction. This is important given that many classes of evolved binaries are thought to be the product of a CE phase, such as cataclysmic binaries consisting of an accreting WD and a Roche-lobe filling low-mass main-sequence star (Knigge et al. 2011), binary central stars of planetary nebulae (Miszalski et al. 2009) and close double white-dwarf binaries (Iben & Tutukov 1984a). All these systems currently have orbital periods ranging from hours to a few days, but must initially have been wide enough to accommodate the red-giant progenitor of a WD, so that drastic orbital shrinkage must have occurred during a CE phase (Paczynski 1976). Furthermore, the double-degenerate formation scenario for type Ia supernovae (SN Ia) invokes the merger of two sufficiently massive carbon-oxygen WDs, brought together into a close orbit during a CE phase and eventually merging due to gravitational-wave emission (Iben & Tutukov 1984b; Webbink 1984). The formation rates of all these systems and events may increase as a result of angular momentum-loss during wind interaction of their wide-orbit progenitors. In addition, enhanced angular-momentum loss and high accretion efficiencies during wind mass transfer may increase the number of potential SN Ia progenitors via the wide symbiotic channel in the single-degenerate scenario (Hachisu et al. 1999).

Finally, we note that although our simulations strictly apply to low-mass binaries with AGB donor stars, the same physical processes are likely to occur in other binaries in which a star loses mass via stellar winds. Whenever the wind velocity is similar to or lower than the orbital velocity, enhanced angular-momentum loss from the orbit may occur. In particular, this may apply to massive binaries under two circumstances: when the mass-losing star is a red supergiant in a very wide binary, or when it is a compact Wolf-Rayet star in a very close binary orbiting a compact object or another Wolf-Rayet star. These binary configurations occur as intermediate stages in several of the progenitor scenarios proposed for the mergers of binary neutron stars and black holes (e.g. Mandel & de Mink 2016; Tauris et al. 2017; van den Heuvel et al. 2017). Dedicated simulations of stellar-wind interaction in such massive binaries are needed to quantify the effect on the evolution of their orbits and the possible consequences for the detection rates of gravitational waves by compact binary inspirals (Belczynski et al. 2002; Chruslinska et al. 2018).

6. Summary and conclusions

We have performed hydrodynamical simulations of low-mass binaries in which one of the components is an AGB star losing mass by a stellar wind at a constant rate and with constant velocity. The companion star is represented by a sink mass with radius equal to a fraction (0.06–0.1) of its Roche lobe and is located at separations of 3–5 AU from the AGB star. We have performed simulations for different ratios of the terminal wind velocity to the relative orbital velocity v_∞/v_orb in order to study the effect of wind mass loss on the orbits of the binary by determining the specific angular momentum of the mass that is lost and the mass-accretion rate onto the companion.

We find two regimes of interaction in terms of v_∞/v_orb. For cases in which v_∞ < v_orb an accretion disk is formed around the companion star, as well as two spiral arms around the system that merge at larger distances. The inner spiral arm wraps closely around the mass-losing star and a consists of high-density gas. In these systems, we also find a high value for the accretion efficiency β, and a high angular-momentum loss per unit mass of material lost from the system. The values of both quantities increase with decreasing velocity of the wind relative to the orbital velocity. For v_∞ > v_orb the BHL accretion regime is approached and the angular-momentum loss is smaller than for cases with v_∞ < v_orb. In these models, only one or two spiral arms are observed with relatively low gas density. No accretion disk is found in any of these models; if a disk forms it must be smaller than the assumed sink radius. Our models indicate that the exchange of angular momentum occurs at close distances from the centre of mass of the binary and is mainly caused by the torque between the binary and the gas in the accretion wake that forms behind the accretor. The strength of this torque increases with the gas density in the wake and its misalignment angle, both which are larger for lower values of v_∞ < v_orb.

Even though the orbital separations chosen in our simulations are in a regime where WRLOF is expected, we do not find the characteristic flow through the inner Lagrangian point encountered in WRLOF simulations because we impose a constant-velocity outflow from the AGB star. Nevertheless, we find similarly high accretion efficiencies in our low windvelocity models as in WRLOF simulations, in which the gradual acceleration of the wind is modelled explicitly. Our work suggests that gravitational focusing by the companion, in combination with efficient gas cooling, are the main processes that result in a high accretion efficiency during wind mass transfer. We note, however, that accretion rates are difficult to determine accurately from our simulations, in contrast to the angular momentum loss rates which are computationally very robust.

Based on the results we obtain for the mass-accretion efficiency and the specific angular momentum of the material lost, we predict the effect on the orbital separation. We find that for v_∞ < v_orb the orbits will shrink and when v_∞ > v_orb the orbits will widen. Furthermore, for the same ratio of wind velocity to orbital velocity we find approximately the same value for the specific angular momentum of the material lost. This suggests that the velocity ratio is the main factor determining the orbital evolution in systems undergoing wind mass transfer. Our results can help explain the puzzling orbits of post-AGB binaries, such as Ba stars and CEMP-s stars, but this has to be verified by binary population synthesis models. Our results also suggest that systems that are initially too wide to undergo Roche-lobe overflow can enter a CE phase, which will have consequences for the expected formation rates of systems such as cataclysmic variables, Type Ia supernovae and white-dwarf mergers producing gravitational waves.

Acknowledgments

We would like to thank the anonymous referee for the very helpful comments that improved this paper. We also would like to thank Carlo Abate, Zhengwei Liu, and Richard Stancliffe for the interesting discussions about the method and the results presented in this paper. Additionally, we would like to thank Sebastian Ohlmann for discussions on the hydrodynamics method. MIS would like to thank Frank Verbunt, Thomas Wijnen, and Andrei Igoshev for the interesting discussions on the topic of this paper.

References

All Tables

All Figures

Fig. 1. Top: flow structure in the y = 0 plane for different EoS after 7.5 orbits. The colourmap shows the temperature at y = 0. The white arrow on top of the velocity field plots corresponds to the magnitude of a 50 km s−1 velocity vector. In the three images we can observe the wind leaving the donor star radially at x = 1 AU. The accretor star is located at x = −2 AU. For model T1, the temperatures reached in the region near the companion star are very high preventing material to settle into an accretion disk. Bottom: gas density in the orbital plane (z = 0) is shown. The accretor is located at x = −2, y = 0, and the AGB star at x = 1, y = 0. In the middle panel we see that for the isothermal EoS (T2) gas in the vicinity of the accretor settles into an accretion disk. When cooling is included (T3), an accretion disk around the companion is also formed.

In the text

	Fig. 2. Top: average mass loss per orbit for the standard model V15a5 measured at different radii from the centre of mass of the binary as indicated by different colours. Bottom: corresponding angular-momentum loss. The black thick curve for rS = 15 AU corresponds to the radius at which we measure J̇bin.
In the text

	Fig. 3. Mean values of the β and η parameters, relating to mass and angular momentum loss, as a function of the mass of the SPH particles used to test different resolutions for the simulations. The error bars correspond to the overall standard error of the mean σm multiplied by a factor of 5 for better appreciation.
In the text

	Fig. 4. Gas density in the orbital plane (z = 0) for V10a5, V10a5s2, V15a5, V19a3, V30a5 and V150a5 after 7.5 orbits.
In the text

Fig. 5. Left panel: density in units of g cm−3 in the orbital plane of the region centred on the accreting star for three models with v∞/vorb < 1, zooming in on the accretion disks. The size of the box corresponds to the Roche-lobe diameter of the companion star. For models V15a5 and V19a3, the situation after 7.5 orbits is shown, whereas for model V10a5 we show the accretion disk before it disappears, but in the same orbital phase. Right panel: mass of the accretion disk, i.e. the mass contained in a sphere of radius Rdisk (see Table 3) as a function of time.

In the text

	Fig. 6. Average accretion rate over one-year intervals for models V15a5 (top), V10a5 (middle) and V30a5 (bottom) is shown in blue. The dashed red lines correspond to the year-averaged mass accretion rate for the corresponding models with a smaller sink radius.
In the text

	Fig. 7. Net inflow rate into a sphere of radius 0.4RL,2 around the accretor star averaged over intervals of one year for the same models as in Fig. 6.
In the text

Fig. 8. Fraction of mass accreted by the companion as a function of v∞/vorb. The dotted line corresponds to BHL accretion rate with BHL = 1, and the dashed line to the accretion rate for a geometrical cross section of size Rsink. Dots correspond to models in which the orbital separation is 5 AU and stars to models in which a = 3 AU. The values shown are measured at a radius of 0.4 RL,2 (β0.4RL,2).

In the text

	Fig. 9. Specific angular momentum of the mass lost from the system as a function of v∞/vorb. The lower the velocity of the wind compared to the orbital velocity, the higher the specific angular momentum loss. For equal velocity ratios, the specific angular momentum is very similar. The dashed line shows the value ηiso expected for an isotropic wind for a mass ratio q = 2.
In the text

Fig. 10. Comparison of our results (Fig. 9) to other work. The data points correspond to the specific orbital angular momentum lost from the system via winds obtained using different methods. The open squares correspond to the results of Chen et al. (2018), triangles show the results of Jahanara et al. (2005) for the radiatively driven wind mode. In the latter case we have interpreted the wind velocity at the Roche-lobe surface as corresponding to v∞. The results from both papers have been corrected for the spin angular-momentum loss (see text) to compare with our results. The colours correspond to different mass ratios q and the lines to the fit formula described by Eq. (13). The blue dotted line is the fit for q = 2, and the orange dashed line for q = 1.

In the text

	Fig. 11. Evolution of the orbital separation as a function of v∞/vorb for the systems simulated in this paper. The systems with ȧ/a < 0 will shrink whereas the ones with ȧ/a > 0 will expand.
In the text