© The Authors 2025

1. Introduction

Approximately 50 percent of low- and intermediate-mass stars (1.5 M_⊙≲M≲5 M_⊙) in the Galaxy are part of a binary system (Duchêne & Kraus 2013), where interactions between the two components can impact their evolution and orbital parameters. An example of such systems is barium (Ba) stars, which are main-sequence or giant stars with surface enrichment in heavy elements produced during the so-called s-process (Jorissen et al. 1998; Van der Swaelmen et al. 2017). Nearly all Ba stars are found in binary systems with a white dwarf (WD) companion, leading to the widely accepted paradigm that the observed pollution in heavy elements occurred when the WD was on the asymptotic giant branch (AGB) and transferred the products of its nucleosynthesis to its binary companion (for example, McClure et al. 1980; McClure 1984; Jorissen et al. 2019). As the AGB companion evolves into a faint WD, only the Ba star remains detectable. The mass transfer responsible for the chemical pollution is also likely to affect the orbital properties of these systems. Observations of Ba stars show that systems with orbital periods less than 2000–3000 days are able to maintain a high eccentricity, at odds with current predictions from tidal models (Zahn 1977), which predict efficient circularization. These systems inherited their orbital properties from past binary interactions that most likely occurred during the AGB phase, when the primary reached its largest extent and was ejecting its envelope via dust-driven winds.

Hydrodynamical simulations of binary systems involving a primary AGB star (for example, Theuns & Jorissen 1993; Mastrodemos & Morris 1998; Kim & Taam 2012; Maes et al. 2021) have shown that if the AGB wind velocity is slower than the orbital speed, the gravitational pull of the companion star alters the shape of the outflow. The wind outflow shifts from an isotropic morphology to one that is directed toward the companion star, and for this resemblance to RLOF it earned the name wind Roche lobe overflow (WRLOF; Mohamed & Podsiadlowski 2007). The modified outflow pattern results in an increased accretion efficiency onto the companion star. Furthermore, winds may escape through the second Lagrangian point (L₂), carrying an increased amount of angular momentum, leading to a shortening of the orbital period, contrary to the case of isotropic, fast winds (Saladino et al. 2019; Chen et al. 2018). Abate et al. (2018) argued that increased angular momentum extraction by winds could be an important factor explaining the orbital properties of CEMP-s stars, the low-metallicity counterparts of Ba stars.

Winds escaping through L₂ also give rise to the formation of a spiral shock behind the companion star, potentially evolving into a circumbinary disk (CB disk) (Mohamed & Podsiadlowski 2007; Chen et al. 2018). Moreover, the spectral energy distributions of all post-AGB binaries indicate the presence of a Keplerian CB disk (Bujarrabal et al. 2013; Oomen et al. 2018). The interaction between the binary system and the CB disk is commonly modeled using linear perturbation theory involving Lindblad resonances, which are known to induce eccentricity in the orbits (Goldreich & Tremaine 1979; Artymowicz & Lubow 1994). However, many studies have struggled to explain the eccentricities of post-AGB binaries with periods below 1000 days without assuming disk masses that exceed observed values (Rafikov 2016; Oomen et al. 2020). Although no better model for long-term binary–CB disk interaction currently exists, several recent smoothed-particle hydrodynamics (SPH) simulations suggest that the process may be more complex, involving accretion streams from the CB disk that enter the binary cavity, leading to mass deposition on both binary components (Muñoz et al. 2019; Tiede et al. 2020; D’Orazio & Duffell 2021; Zrake et al. 2021).

When initiated on an eccentric orbit, RLOF mass transfer is phase dependent, increasing near periastron and decreasing near apastron. This phase dependence can generate eccentricity that effectively counteracts the circularization induced by tidal forces. Vos et al. (2015) suggested that this principle could explain the properties of long-period sdB binaries. A similar mechanism, known as tidally enhanced wind mass loss, proposed by Tout & Eggleton (1988), introduces a phase dependency in the wind mass loss, resulting in an increased ejection rate at periastron. Akin to the phase-dependent RLOF, this modified wind mass loss can generate eccentricity in binary systems. Furthermore, the enhanced mass loss accelerates the evolutionary pace of the AGB star, reducing the time available for tidal forces to circularize the orbit. Consequently, the star's evolution may terminate earlier with a smaller radius, enabling it to reach shorter orbital periods before initiating RLOF. This mechanism was invoked by Siess et al. (2014) to explain the properties of the He white dwarf system IP Eri. Additionally, Bonačić Marinović et al. (2008) utilized this mechanism, coupled with phase-dependent RLOF, to elucidate the origin of eccentric post-AGB binaries.

This study provides one of the first investigations into the long-term impact of WRLOF, CB disk, and tidally enhanced wind mass loss on the orbital evolution of potential progenitor systems of Ba stars. To achieve this, we use the binary evolution code BINSTAR (Siess et al. 2013), which incorporates these key mechanisms. The paper is organized as follows: In Sect. 2, we present our modeling of WRLOF, phase-dependent RLOF, and our model of the CB disk and its interaction with the binary. In Sect. 3, we estimate the range of initial orbital parameters for the systems to experience WRLOF. In Sect. 4, we present our results, beginning with an overview of the observational data sample, followed by an analysis of a representative system. We then present the evolutionary tracks in the e−log P diagram for a particular set of parameters, followed by an exploration of the effects of varying these parameters. We also analyze the effect of tidally enhanced wind mass loss and of changing the initial mass of the primary. Finally, we estimate the effects of the CB disk after the end of the AGB phase until its dissipation, and the effect of non-conservative RLOF. In Sect. 5, we discuss our findings, and in Sect. 6, we provide a conclusion.

2. The binary paradigm

2.1. The Wind Roche Lobe Overflow (WRLOF)

Hydrodynamical simulations of the evolution of an AGB primary with a lower-mass companion (for example, Mohamed & Podsiadlowski 2007; de Val-Borro et al. 2009; Chen et al. 2017; Saladino et al. 2019) indicate that if dust forms near or beyond the AGB's Roche lobe (i.e., $R_{\text{dust}}/R_{L_1}\gtrsim 1$), then the radiatively driven wind is not accelerated to the AGB's escape velocity. Consequently, the wind is not ejected isotropically, but channels through the inner Lagrangian (L₁) point and focuses onto the companion star. The flow geometry resembles the classical RLOF scenario and has since been coined wind Roche lobe overflow. A fraction of the wind is accreted onto the companion, leading to a significant increase in the companion's mass. The remaining material escapes the system through the L₂ point behind the AGB's companion, carrying more angular momentum compared to the isotropic wind mass loss case, producing a shrinkage of the orbit (Jahanara et al. 2005; Chen et al. 2018; Saladino et al. 2019). Mohamed & Podsiadlowski (2012) also found that the gas escaping through L₂ leads to the formation of a single spiral arm around the binary. The spiral is compressed, and the enhanced density can trigger dust formation. Whether this structure subsequently evolves into a stable CB disk has been claimed (Chen et al. 2017, 2020).

In this study, our paradigm is that WRLOF leads to the formation of a CB disk as a result of material channeling through the L₂ point. In our models, we assume that WRLOF and disk formation occur once $R_{\text{dust}}/R_{L_1}>1$. We estimate the dust formation radius R_dust using the formulation of Höfner (2007):

$$R_{\text{dust}}=\frac{1}{2}{R}_1{\left(\frac{T_{\mathrm{eff},1}}{T_{\text{dust}}}\right)}^{2.5},$$(1)

where T_{eff, 1} and R₁ are the AGB's effective temperature and radius, and T_dust the dust condensation temperature. For carbon-rich (C/O>1) stars, T_dust≈1500 K, and for oxygen-rich stars (C/O <1), T_dust≈1000 K.

2.2. Evolution of the orbital parameters

The change in the orbital separation is determined from the rates of change of the donor and gainer masses, ${\dot{M}}_1$ and ${\dot{M}}_2$ respectively, the orbital angular momentum loss rate, ${\dot{J}}_{\text{orb}}$, and the rate of change of the eccentricity ė, i.e.

$$\frac{\dot{a}}{a}=2\frac{{\dot{J}}_{\text{orb}}}{J_{\text{orb}}}-2\frac{{\dot{M}}_1}{M_1}-2\frac{{\dot{M}}_2}{M_2}+\frac{{\dot{M}}_1+{\dot{M}}_2}{M_1+M_2}+\frac{2e\dot{e}}{1-e^2},$$(2)

where

$$\dot{e}={\dot{e}}_{\text{tides}}+{\dot{e}}_{\text{res}}+{\dot{e}}_{\text{winds}}+{\dot{e}}_{\text{RLOF}}$$(3)

is the sum of the contributions from tides ${\dot{e}}_{\text{tides}}$, calculated using the prescriptions of Zahn (1977, 1989), the resonant interactions with the CB disk ė_res (Eq. (21)) and the exchanges of mass between the stars by the winds ${\dot{e}}_{\text{winds}}$ and RLOF ${\dot{e}}_{\text{RLOF}}$ (Eq. (16)).

The total angular momentum of the binary, J_Σ, is the sum of the orbital angular momentum, J_orb and the spin angular momenta of each star, J_{1, 2}

$$J_{\mathrm{\Sigma }}=J_{\text{orb}}+J_1+J_2.$$(4)

Taking the time derivative of Eq. (4) and solving for ${\dot{J}}_{\text{orb}}$ gives

$${\dot{J}}_{\text{orb}}=-{\dot{J}}_1-{\dot{J}}_2+{\dot{J}}_{\mathrm{\Sigma }}.$$(5)

The torques on the stellar spins ${\dot{J}}_{1,2}$ come from tidal interaction and mass transfer and are given in Siess et al. (2013).

2.3. Wind mass transfer

In the WRLOF regime, the wind accretion efficiency onto the companion, ${\beta }_{\text{wind}}=|{\dot{M}}_{\text{acc},2}/{\dot{M}}_1|$, and the angular momentum carried away from the system by the wind, ${\dot{J}}_{\text{wind}}$, are increased compared to the standard (non-dusty) prescription. To parameterize ${\dot{J}}_{\mathrm{\Sigma }}$, we use the prescription derived by Saladino et al. (2019) based on SPH simulations. In their formulation, the systemic angular momentum loss is given by

$${\dot{J}}_{\mathrm{\Sigma }}=\eta \frac{J_{\text{orb}}}{\mu }\dot{M},$$(6)

where μ=M₁M₂/(M₁+M₂) is the reduced mass and η a function of the ratio of the wind to orbital velocities (v_wind/v_orb) and of the mass ratio q=M₁/M₂ given by

$$\eta \left(\frac{{\mathit{v}}_{\text{wind}}}{v_{\text{orb}}},q\right)=\text{min}\left(\frac{1}{c_1+(c_2{\mathit{v}}_{\text{wind}}/v_{\text{orb}}{)}^3}+{\eta }_{\text{iso}},0.6\right),$$(7)

where c₁=max(q, 0.6q^1.7), c₂ = 1.5+0.3q and η_iso = 1/(1+q)² is the fraction of the specific orbital angular momentum lost by the wind in the isotropic (fast) wind case.

We can assess the effect of fast winds on the orbital period by substituting $\dot{J}\text{orb}$ in Eq. (2) with ${\dot{J}}_{\mathrm{\Sigma }}$ in Eq. (6), and neglecting accretion onto the companion, ${\dot{M}}_2=0$. For a circular orbit, we find

$$\frac{\dot{a}}{a}=\frac{2\hspace{0.5em}\dot{M}}{q\hspace{0.5em}M}\left(\eta \hspace{0.5em}(1+q{)}^2-(1+\frac{q}{2})\right).$$(8)

In the limit v_wind≫v_orb, η approaches the fast (isotropic) wind case η_iso, resulting in

$$\frac{\dot{a}}{a}=-\frac{\dot{M_1}}{M},$$(9)

showing that the orbit widens as the primary loses mass, independently of the mass ratio. Following Siess et al. (2014), this equation can be integrated, and using Kepler's third law, the evolution of the orbital period can be expressed in terms of the mass ratio

$$P=P_0{\left(\frac{1+q_0}{1+q}\right)}^2,$$(10)

where P, P₀, q and q₀ are the current and initial period and mass ratio, respectively. Equation (8) also tells us that if $\eta >{\eta }_{\text{crit}}=(1+\frac{q}{2}){\eta }_{\text{iso}}$, systemic wind mass loss produces a shrinkage of the orbit.

The top two panels of Fig. 1 show the dependence of η on the orbital period for two systems characterized by (M₁, M₂) = (2,1) and (1,1), respectively (with masses expressed in  M_⊙). The wind velocity is a free parameter, and for AGB stars, we use representative values of 6, 10 and 15 km s⁻¹. For both values of q and for v_wind = 6 and 10 km s⁻¹, systems with P_orb<20 000 d have η>η_crit, meaning that systemic wind mass loss will tighten the orbit. For v_wind = 15 km s⁻¹, η<η_crit only for systems with periods ≳6000 d. In wider binaries, v_wind≫v_orb so the fast wind regime prevails and mass loss pushes the stars apart.

Fig. 1. Top two panels: Fraction of the specific angular momentum η carried away by winds during WRLOF as a function of orbital period for a 1+1 M⊙ (top panel) and 2+1 M⊙ (top-middle panel) binary system, and for different wind velocities (vwind, as labeled by the different colors). The horizontal long-dashed line represents the value of η in the case of isotropic wind emission (Jeans mode) and the short-dashed line the value of ηcrit above which the orbit shrinks as a result of wind mass loss. Bottom two panels: Accretion efficiency onto the companion, $\beta =|{\dot{M}}_{\text{acc},2}/{\dot{M}}_1|$ for the same two systems as in the upper panels, either through wind (βBHL, dashed lines) or through WRLOF (βWRLOF, solid lines).

The wind accretion efficiency onto the companion is generally estimated from the Bondi-Hoyle-Littleton prescription (e.g. Boffin & Jorissen 1988)

$${\beta }_{\text{BHL}}=\frac{{\alpha }_{\text{BHL}}}{(1-e^2{)}^{\frac{1}{2}}(1+q{)}^2}{\left(\frac{{\mathit{v}}_{\text{orb}}}{{\mathit{v}}_{\text{wind}}}\right)}^4{\left[1+{\left(\frac{{\mathit{v}}_{\text{orb}}}{{\mathit{v}}_{\text{wind}}}\right)}^2\right]}^{-\frac{3}{2}},$$(11)

where α_BHL = 0.75 is a parameter, ${\mathit{v}}_{\text{orb}}=(\mathit{GM}/a{)}^{\frac{1}{2}}$ the average orbital velocity, and v_wind is the wind velocity given by Eq. (3) of Vassiliadis & Wood (1993). This prescription should be restricted to fast winds, where v_wind≫v_orb. As shown in hydrodynamical simulations, when v_wind≲v_orb, the wind is focused towards the orbital plane, leading to an increased accretion efficiency onto the companion. Based on their SPH simulations, Saladino et al. (2019) derived a prescription for the accretion efficiency β_WRLOF onto the companion that is given by

$${\beta }_{\text{WRLOF}}=\text{min}(\alpha {\beta }_{\text{BHL}},{\beta }_{\text{max}}),$$(12)

where β_BHL is computed using Eq. (11) with α_BHL = 1, and

$$\alpha =0.75+\frac{1}{k_1+(k_2{\mathit{v}}_{\text{wind}}/{\mathit{v}}_{\text{orb}}{)}^5},$$(13)

where k₁ = 1.7+0.3 q, k₂ = 0.5+0.2 q, and β_max=min(0.3, 1.4 q⁻²), the maximal accretion efficiency found in their simulations. The bottom two panels in Fig. 1 show β_WRLOF and β_BHL as a function of the orbital period for two different mass ratios and different wind velocities. As expected, the accretion efficiency decreases for higher values of v_wind and in systems with longer orbital periods, where v_orb is slower.

Variations of the mass loss rate along an elliptical orbit induce a change in eccentricity that can be expressed as (Eggleton 2011)

$$\begin{array}{cc}\hfill {\dot{e}}_{\text{wind}}& =[\left(1-{\beta }_{\text{wind}}\right)\frac{|{\dot{M}}_{\text{wind}}|}{M}+2\hspace{0.5em}{\beta }_{\text{wind}}\hspace{0.5em}|{\dot{M}}_{\text{wind}}|\left(\frac{1}{M_1}-\frac{1}{M_2}\right)]\hfill \\ \hfill & \times (e+cos\nu ),\hfill \end{array}$$(14)

where ν is the true anomaly and β_wind the efficiency of wind accretion onto the companion. When the evolutionary timestep, Δt>P_orb, Eq. (14) is averaged over the orbital period using the procedure outlined in Siess et al. (2014). The net change in eccentricity is nonzero only if the bracketed term depends on the orbital phase. This happens for example for tidally enhanced winds (Sect. 4.5) and, in this case, the systemic mass loss (first term in brackets) always leads to eccentricity pumping, while wind accretion onto the companion (second term) has a positive contribution only when M₁<M₂ and negative in the opposite case. The total change in eccentricity due to wind mass loss should include the contributions from both stellar components.

2.4. Change in a and e due to Roche lobe overflow

As the primary ascends its AGB, it may eventually fill its Roche lobe and initiate RLOF. The expansion of the convective envelope of the AGB upon mass loss may lead to dynamically unstable mass transfer and common envelope evolution, unless the Roche lobe expands at a similar rate as the star's radius.

Assuming conservative RLOF (${\dot{M}}_1={\dot{M}}_{\text{RLOF}}=-{\dot{M}}_2$ and ${\dot{J}}_{\text{orb}}=0$), Eq. (2) writes

$$\frac{\dot{a}}{a}=2\hspace{0.5em}{\dot{M}}_{\text{RLOF}}\left(\frac{1}{M_1}-\frac{1}{M_2}\right)+\frac{2\hspace{0.5em}e\hspace{0.5em}\dot{e}}{1-e^2}.$$(15)

If the orbit is eccentric, the change in e due to mass loss is also given by Eq. (14), replacing ${\dot{M}}_{\text{wind}}$ by ${\dot{M}}_{\text{RLOF}}$ and β_wind by ${\beta }_{\text{RLOF}}=-{\dot{M}}_1/{\dot{M}}_2$ ( = 1 in our conservative approach), to give (see also Soker 2000)

$${\dot{e}}_{\text{RLOF}}=[2{\dot{M}}_{\text{RLOF}}\left(\frac{1}{M_1}-\frac{1}{M_2}\right)]\times (e+cos\nu ).$$(16)

It follows that if RLOF occurs when M₁>M₂, both the separation and the eccentricity decrease (the second term in Eq. (15) is generally negligible for this conclusion to hold). However, when the mass ratio is reverted, RLOF widens the orbit and the eccentricity increases.

2.5. Resonant interactions with the circumbinary disk

The resonant interactions between the CB disk and the binary have been studied by Goldreich & Tremaine (1979) and Artymowicz & Lubow (1994) using linear perturbation theory. They found that the binary can exert a positive or negative torque on the CB disk at resonant locations, which, in turn, affects the orbital parameters. The potential of the binary Φ(r, θ, t) (where r is the distance from the CB disk center located at the binary's mass center, and θ is the azimuthal angle) is developed in series as

$$\mathrm{\Phi }(r,\theta ,t)={\displaystyle \sum _{m\ell }}{\varphi }_{m\ell }(r)exp[i(m\theta -\ell {\mathrm{\Omega }}_{\text{B}}t)],$$(17)

where l is the time-harmonic number, m≥0 is the azimuthal number, ϕ_mℓ(r) is the radially varying potential component for a given (m, ℓ) and Ω_B = 2π/P_orb is the binary's mean orbital angular speed, with P_orb the orbital period. Each harmonic ϕ_mℓ(r) is associated with an inner Lindblad resonance (ILR) that produces a positive torque on the binary and an outer Lindblad resonance (OLR) that exerts a negative torque. The resonances are located at radii

$$r_{\mathit{LR}}={\left(\frac{m\pm 1}{l}\right)}^{\frac{2}{3}}a,$$(18)

where the positive and negative signs refer to the OLR and ILR, respectively, and a is the binary semi-major axis.

For small eccentricities (e≲0.2), SPH simulations (Artymowicz & Lubow 1994) show that the (2, 1) OLR is the dominant component for the angular momentum exchange between the CB disk and the binary. From now on, we write results relevant only to this resonance. A more general form of the equations, valid for any other resonance is presented in Appendix A. Goldreich & Tremaine (1979) computed the torque ${\dot{J}}_{\text{res}}$ exerted by the CB on the binary assuming that (i) the disk is thin (H/r≲0.1, where H is the disk's pressure scale height), (ii) the disk is coplanar with the binary orbit, and (iii) the potential of the non-axisymmetric perturbations of the disk's structure is small compared to that of the binary. The torque, evaluated at the (2, 1) outer resonance, and in the limit of small e reads

$${\dot{J}}_{\text{res}}=-\frac{49{\pi }^2}{8}G\frac{{\mu }^2}{M}a\mathrm{\Sigma }(r_{\text{res}})e^2,$$(19)

where Σ(r_res) is the disk's surface density at the (2, 1) OLR located at r_res = 3^2/3a≈2.08 a. The effect of the torque on the binary separation writes

$${\left(\frac{\dot{a}}{a}\right)}_{\text{res}}=-\frac{49{\pi }^2}{8}\mathrm{\Sigma }(r_{\text{res}})\frac{\mu }{M}\sqrt{\mathit{Ga}}\hspace{0.5em}{e}^2\hspace{0.5em},$$(20)

where μ=M₁M₂/(M₁+M₂) is the reduced mass. Finally, the associated rate of change of the eccentricity is

$${\dot{e}}_{\text{res}}=-\left(\frac{1-e^2}{e}\right)\left[\frac{1}{2}-\frac{1}{(1-e^2{)}^{\frac{1}{2}}}\right]{\left(\frac{\dot{a}}{a}\right)}_{\text{res}}.$$(21)

The (2, 1) OLR acts to shrink the orbit and increase the eccentricity of the binary. Substituting Eq. (20) into Eq. (21), it becomes apparent that the eccentricity pumping is proportional to e and vanishes for circular orbits.

In more eccentric systems, higher order OLRs and ILRs are activated. The inner Lindblad resonances have the opposite effect on the orbit than the outer ones, and taking them into account, the eccentricity pumping is damped. Roedig et al. (2011) found that for black hole binaries with e≈0.6−0.8 the eccentricity pumping roughly cancels and e saturates. A similar behavior is also found in the 2D SPH simulations performed by D’Orazio & Duffell (2021), where they explain this saturation by the deformation of the disk when e>0.4 (see also Valli et al. 2024). We account for the damping effects of the LRs at high eccentricities by multiplying Eqs. (19) and (20) by a Fermi-like function, ℱ(e), given by

$$\mathcal{F}(e)=\frac{1}{exp\left(\frac{e-0.3}{0.04}\right)+1}.$$(22)

This allows us to reproduce the saturation effect found by Valli et al. (2024), i.e. $\dot{e}\to 0$ when e>0.4.

The disk's inner edge r_in is localized where the outwards torque of the (2, 1) OLR compensates the disk's inwards viscous torque. To calculate r_in, we use Eq. 10 of Dermine et al. (2013), which reads¹

$$\frac{r_{\text{in}}}{a}=1.7+\frac{3}{8}\sqrt{e}log(\mathcal{R}),$$(23)

where ℛ=(αγ²)⁻¹ is the Reynold's number, α is the viscosity parameter of Shakura & Sunyaev (1973), and γ=H/r_in with H the disk's pressure scale height.

Our expression for ė_res contrasts with those of Dermine et al. (2013), Vos et al. (2015) and Oomen et al. (2020), where in their formulation ė_res depends on the total angular momentum of the disk and on the viscous timescale. They used an approximation for ė_res proposed by Lubow & Artymowicz (1996) that is valid when the disk structure has relaxed and is close to a steady state, which may be appropriate for a CB disk in post-AGB systems. The fact that in our model, matter is continuously fed into the disk, often faster than the viscous timescale, implies that the disk may not be in a steady state.

2.6. The surface density at the OLR

The value of the surface density at the resonance Σ(r_res) can be calculated by solving the time-dependent disk evolution equations, with appropriate boundary conditions (see, e.g. Pringle 1991). The modeling of the disk's structure in such detail is beyond the scope of the current investigation. Instead, we use a simple model to simulate the evolution of Σ(r_res):

$$\frac{\partial \mathrm{\Sigma }(r_{\text{res}},t)}{\partial t}=\xi \hspace{0.5em}\frac{\left|\dot{M}(t)\right|}{4\pi a^2}-\frac{\mathrm{\Sigma }(r_{\text{res}},t)}{{\tau }_{\text{visc}}},$$(24)

where $\dot{M}(t)<0$ is the systemic mass loss rate. The first term on the right-hand side describes the injection of matter at the inner edge of the disk in a shell of radius Δr∝a. The parameter ξ describes both the fraction of mass lost by the binary that is used to feed the disk, and the unknown distribution of matter near the disk edge. As will be seen in Sect. 4.4.2, typical values for ξ range between 10⁻⁴ and 5×10⁻⁴. For smaller values of ξ the effect of the CB disk is negligible, while for larger ones, the density at the inner rim (Σ_in) is higher and leads to disk masses above 10⁻²  M_⊙, in contradiction with current estimates in post-AGB stars (Bujarrabal et al. 2013; Oomen et al. 2020; Bollen et al. 2022).. The second term expresses the spreading of the disk and the associated decrease in Σ(r_res) by viscous torques that transport mass and angular momentum outward. The viscous timescale (Spruit & Taam 2001)

$${\tau }_{\text{visc}}=\frac{4r_{\text{in}}^2}{3{\nu }_{\text{in}}}=\frac{2}{3\pi \alpha {\gamma }^2}{\left(\frac{r_{\text{in}}}{a}\right)}^{\frac{3}{2}}{P}_{\text{orb}}$$(25)

is evaluated at the disk's inner edge r_in. Here, α is the Shakura & Sunyaev (1973) disk parameter, ν_in the viscosity at r_in given by

$${\nu }_{\text{in}}=\alpha c_{s}H=\alpha {\gamma }^2\sqrt{{\mathit{GMr}}_{\text{in}}},$$(26)

and c_s is the sound speed. Integration of Eq. (24) yields

$$\mathrm{\Sigma }(r_{\text{res}},t+\mathrm{\Delta }t)=\mathrm{\Sigma }(r_{\text{res}},t)e^{-\mathrm{\Delta }t/{\tau }_{\text{visc}}}+\xi \frac{\left|\dot{M}\right|}{4\pi a^2}{\tau }_{\text{visc}}\left(1-e^{-\mathrm{\Delta }t/{\tau }_{\text{visc}}}\right)$$(27)

where Δt is the time-step. The disk temperature was estimated by Hillen et al. (2014) using radiative-transfer calculations on a parametric CB disk structure. Their calculations indicate a temperature at the inner edge of ≈1600 K, declining with distance. Since this is below the ionization temperature of hydrogen, we assume α = 0.01 in Eq. (26), appropriate for cool disks where the magneto-rotational instability, thought to drive the viscosity mechanism, is suppressed (Gammie 1996). Additionally, studies of CB disk formation around cataclysmic variables (Dubus et al. 2002) suggest that γ≲0.05. We therefore take γ = 0.05 and find t_visc≈10⁵ yr for P_orb = 1000 d, which gives a sensible disk lifetime.

2.7. Stellar models

We examine the evolution of two representative binaries consisting of solar metallicity stars with initial masses of 1.2+0.8 M_⊙ and 2.0 + 1.0 M_⊙. We did not consider more massive primaries (M_i≥3  M_⊙) because, as recently analyzed by Karinkuzhi et al. (2018), the s-process abundances in Ba stars do not comply with the expected nucleosynthesis arising from the ²²Ne(α,n) neutron source operating in these more massive AGB stars. The first system (1.2+0.8) considers the least massive primary for the star to evolve off the AGB within the Galactic lifetime. These lowest-mass primaries are also those that reach the smallest radius at the AGB tip and thus may help populate the short-period range of Ba stars.

The simulations are carried out using the stellar evolution code BINSTAR (Siess et al. 2013). In the simulations, we adopt a solar composition (Z = 0.02), according to Grevesse et al. (1996). Convection is treated using the standard mixing length theory, with the mixing length parameter α_MLT = 1.75, which results from fits of solar models. The input stellar physics is the same as described in Siess & Arnould (2008). The mass loss is prescribed according to Schröder & Cuntz (2007) up to the end of core helium burning and then during the AGB phase, we switch to the Vassiliadis & Wood (1993) prescription. Radiative opacities are taken from Iglesias & Rogers (1996) above 8000 K and from Ferguson et al. (2005) at lower temperatures. The nuclear network includes 182 reactions coupling 55 species from H to Cl as described in Siess & Arnould (2008). We also treat overshooting at the base of the convective envelope as a diffusive process, following Herwig (2000) with f_over = 0.01. The evolution starts when the stars are on the zero-age main sequence (ZAMS) with the stellar spins synchronized with the orbital period.

The simulated 2.0+1.0 M_⊙ binaries have initial periods P_i={4000, 6000, 8000, 12 000, 16 000} days and initial eccentricities e_i={0.01, 0.1, 0.3}. For 1.2+0.8 M_⊙ binaries, we compute additional models for P_i = 2000 and 3000 d and the same eccentricities. As will be seen in Sect. 4, binaries with P_orb<2000 d completely circularize during their evolution, and above P_orb>16 000 d WRLOF is not initiated (Sect. 3). These systems are of limited interest and not considered in this study. The effects of the wind speed and CB disk are evaluated by varying v_wind and ξ, whose values are summarized in Table 1. The evolution of the systems was followed until the end of the AGB phase, when convergence problems arose. At this point, less than 0.05 M_⊙ remained in the envelope.

3. Conditions for WRLOF

Using Kepler's third law and the Eggleton (1983) formula for the Roche radius, the condition for WRLOF $R_{\text{dust}}/R_{L_1}=1$ evaluated at periastron when it is most favorable, can be expressed in terms of the orbital period as

$$P=(1-e{)}^{-3/2}{\left(\frac{R_{\text{dust}}}{f(q)}\right)}^{\frac{3}{2}}\frac{2\pi }{\sqrt{G(M_1+M_2)}},$$(28)

where $R_{L_1}\equiv a\hspace{0.5em}(1-e)\hspace{0.5em}f(q)$, and f(q) = 0.49q^2/3/(0.6q^2/3+ln(1+q^1/3)). The orbital period, eccentricity and masses appearing in the previous equation do not correspond to the initial values because mass and angular momentum have been removed from the system by stellar winds. In the absence of WRLOF, wind mass loss likely proceeds in the Jeans mode and the difference between the initial (P₀) and current period (P) at the onset of the WRLOF can be estimated using Eq. (10). Given the long periods of our systems, tides have a negligible effect prior to the onset of WRLOF, so the eccentricity can be considered constant, i.e. e=e₀. Together, Eqs. (10) and (28) give the initial period for the onset of WRLOF

$$P_0^{\text{WRLOF}}=(1-e_0{)}^{-3/2}{\left(\frac{R_{\text{dust}}}{f(q)}\right)}^{\frac{3}{2}}\frac{2\pi }{\sqrt{G(M_1+M_2)}}{\left(\frac{1+q}{1+q_0}\right)}^2.$$(29)

The dust condensation radius R_dust is given by Eq. (1) where in this equation the radius, effective temperature and mass of the primary star are taken from single stellar evolution models computed with the BINSTAR code. We further consider an O-rich dust with T_dust = 1000 K and a circular orbit (e₀ = 0).

Figure 2 shows the evolution of $P_0^{\text{WRLOF}}$ as a function of the primary mass for the two systems considered in this study. The critical initial periods P_crit below which a circular 2.0 + 1.0 M_⊙ binary initiates WRLOF during the red giant branch (RGB) and AGB phases are respectively 9700 d and 20 000 d (indicated by the horizontal dotted lines labeled RGB and AGB P_crit in Fig. 2). For a circular 1.2+0.8 M_⊙ binary, the critical periods are shorter because of the smaller radius reached by lower-mass stars and ranges around 11 000 d and 15 400 d for the RGB and AGB phases respectively. During the AGB phase, WRLOF may be activated episodically during the short-lived thermal pulses (for the 2.0+1.0 M_⊙ system, this situation would appear if 16 500 d ≲ $P_0^{\text{WRLOF}}\lesssim 20000$ d). However, due to the small amount of mass lost during these episodes, the orbital evolution is not significantly impacted, as shown by our tests.

Fig. 2. Evolution of the initial period $P_0^{\text{WRLOF}}$ for WRLOF as a function of the mass of the primary (Eq. (29)) for a 2.0 + 1.0 M⊙ (solid blue line) and a 1.2+0.8 M⊙ (dashed orange line) binary. The spikes in the curves are associated with the development of thermal pulses. The black dotted lines indicate the critical initial period of circular binaries Pcrit below which WRLOF first starts during the RGB or the AGB phase.

The dependence of the critical period for the onset of WRLOF on the RGB or AGB on the eccentricity and T_dust can be estimated from Eq. (29) and is illustrated in Fig. 3. With increasing eccentricity, binaries are more likely to experience WRLOF at any given period because of the reduction of the periastron distance. As a consequence the critical period to start WRLOF at a specific phase increases. For example, in the 2.0+1.0 M_⊙ system, P_crit is 70% longer for e = 0.3, compared to the circular case.

Fig. 3. Regions in the e−log P diagram where WRLOF is active. If the initial period Pi of a system lies to the left of the curve, a binary undergoes WRLOF during the RGB (orange dashed line), the AGB (solid blue line), or the AGB with C-rich dust (i.e. Tdust = 1500 K, green short-dashed line). Top and bottom panels refer to 1.2 + 0.8  M⊙ and 2.0 + 1.0  M⊙ binaries, respectively. The orange dashed curve and solid blue curve are computed using Tdust = 1000 K (i.e. O-rich dust).

Combining Eqs. (1) and (29) implies $P_i\propto T_{\text{dust}}^{-15/4}$. This strong dependence on the dust condensation temperature implies that the critical period for C-rich stars will be a factor of five smaller than for O-rich stars. The consequences of this composition change are quite significant. If the conditions for the occurrence of WRLOF in a system involving an O-rich AGB star are initially met, the fact that the AGB star becomes C-rich, brings R_dust within its AGB Roche lobe and WRLOF stops. To have WRLOF during the C-rich AGB phase, the systems require shorter initial periods that may not be compatible with the absence of RLOF.

4. Results

Before entering the detailed analysis of our models, we first describe our reference data sample of Ba stars that we use for comparison.

4.1. The observational data

The distribution of post-AGB stars, Ba dwarfs and Ba giants in the e−log P diagram is shown in Fig. 4. The data come from Jorissen et al. (2016, 2019), Van der Swaelmen et al. (2017), Escorza et al. (2019) and references therein. The post-AGB sample is complete up to a period of ≈2000 d while, for dwarf and giant Ba stars this limit is pushed to P≈10⁴ d.

Fig. 4. Distribution in the e−log P diagram of the observed post-AGB binaries (top), Ba dwarfs (middle), and Ba giants (lower panel). The different lines delineate the orbital parameters for which a primary of radius 250 R⊙ (solid line), 350 R⊙ (dashed line) and 450 R⊙ (dotted line) fills its Roche lobe. These radii correspond to the maximum radii typically attained by solar-metallicity AGB stars of initial masses 1.2, 2 and 3 M⊙, respectively. Both components are assumed to have 1 M⊙ at the time of the RLOF, although the exact masses do not meaningfully affect the lines. As expected, lower primary-mass systems can reach shorter periods without RLOF.

The distributions (Fig. 4) show a deficit of post-AGB and Ba dwarfs at large periods (above P≳2000 d) compared to the distribution of Ba stars which extends up to 30 000 days. The reason for this difference can be attributed to the incompleteness of the sample for large periods but also to biases. In particular, Ba dwarfs, being fainter, are more difficult to observe, and on the main sequence, intermediate-mass Ba-dwarfs appear as A-type stars which are faster rotators, often variable (e.g., δ Scuti, γ doradus and roAp stars), and known to show chemical anomalies attributed to mixing processes (e.g., Am and Ap stars). Under these conditions, it may be difficult to classify an A-type star in the Ba-dwarf family.

It is important to stress that the post-AGB sample is likely contaminated by post-RGB systems. Kamath et al. (2014) found that approximately two thirds of their Small Magellanic Cloud post-AGB binaries were in fact post-RGB systems. Therefore comparisons between the stellar models and the post-AGB or Ba dwarfs will only be meaningful for relatively long-period systems.

As far as eccentricity is concerned, it seems that between the post-AGB and Ba dwarf stages, the overall eccentricity of the systems has decreased. This could be attributed to tidal effects circularizing the orbit on the long lifetime of the Ba dwarf provided the two classes are evolutionary connected. The evolution between the Ba dwarfs and the Ba giants has been investigated by Escorza et al. (2020) using the BINSTAR code. Adopting as initial conditions the orbital and eccentricity distributions of Ba and CH dwarfs, these authors evolved those systems and were able to reproduce satisfactorily the e−log P distribution of Ba giants.

Another limitation of the modeling comes from the fact that our binary code cannot handle unstable mass transfer, when the convective envelope of the AGB star overfills its Roche lobe. The outcome of such an evolution is uncertain and depends on various parameters, including the mass ratio at the beginning of mass transfer and on the dynamical response of the AGB star. It is generally believed that the expansion of the AGB star in response to mass loss leads to a common envelope phase during which a fraction of the orbital energy is transferred to the envelope, powering its ejection. During this process, the orbital parameters are likely to be affected but the magnitude of the changes are rather hard to estimate.

From single-star models computed with the BINSTAR code (Siess et al. 2013), we find that in the mass range 1.2≤M₁/M_⊙≤3, the maximum radius at the AGB tip, R_AGB−tip, varies between 250 and 450 R_⊙. Therefore, in Fig. 4, systems that are located in the e−log P diagram to the left of the critical line defined by $R_{L_1}=min(R_{\mathrm{AGB},\mathrm{tip}})\approx 250R_{\odot }$ will overfill their Roche lobe either on the RGB or AGB. The location of this critical line slightly depends on the mass ratio and, for circular systems, corresponds to P≈2000 d. Therefore, our comparisons will be limited to those systems with P≳2000 d.

4.2. Analysis of a representative model

In this section, we study the orbital evolution of a 2.0+1.0 M_⊙ binary with P_i = 8000 d, e_i = 0.1, computed using the parameter set AV6 in Table 1. The analysis of this representative model illustrates the various processes at play, providing a foundation for understanding the evolutionary tracks discussed later. The evolution of this system is presented in Fig. 5 and can be divided into five stages, labeled S1–S5. Note that not all models proceed through all these stages. For example, if a binary never undergoes RLOF, it may skip the S4 stage.

Fig. 5. Evolution of a 2.0+1.0 M⊙ binary with initial period Pi = 8000 d and eccentricity ei = 0.1 using the parameter set AV6. We present, as a function of the primary mass M1, (a) the orbital period, Porb; (b) the eccentricity e; (c) the Roche-filling factor $R_1/R_{L_1}$; (d) the primary mass loss rates from the wind ${\dot{M}}_{\text{wind}}$ (dashed curve) and RLOF ${\dot{M}}_{\text{RLOF}}$ (solid curve); (e) the surface density at the CB disk's inner rim, Σin; (f) the change in the orbital eccentricity ė (black curve), which is decomposed into contributions from the CB disk (orange), tides (red), RLOF (green). The evolution of the system is divided into different stages, labeled S1–S5, and separated by the vertical black dotted lines.

S1: Main sequence, RGB and core He-burning phase (epoch mass M₁≳1.9  M_⊙).

The changes in orbital parameters during the main sequence and early RGB (M₁≳1.95  M_⊙) are small (panels a and b). At this stage, the dust condensation radius is within the primary's Roche lobe ($R_{\text{dust}}/R_{{\text{L}}_1}<1$), preventing WRLOF and the formation of a CB disk. The isotropic winds lead to a slight increase in the period, as expected from Eq. (10). On the RGB, both R and R_dust increase significantly and eventually, $R_{\text{dust}}/R_{{\text{L}}_1}>1$, briefly triggering WRLOF at M₁ = 1.95  M_⊙. In this regime, the winds extract a larger amount of angular momentum (Eq. (6)), the orbit tightens and P_orb returns close to its initial value. By the end of the RGB phase, the CB disk surface density Σ_in reaches 0.03 g cm⁻² (panel e). When helium ignites at M₁ = 1.9 M_⊙, the star contracts rapidly within its Roche potential, considerably reducing the wind mass loss rate and the effects of tides which become inefficient (panel f). When resonant interaction with the CB operates, it produces a small increase in eccentricity from 0.1 to 0.13 until the disk dissipates, and Σ_in vanishes (panel e) after a few viscous timescales (with τ_visc of the order of ∼10⁵ yr; Eq. (25)). The orbital evolution of the system up to the end of the core He-burning phase is thus characterized by a small increase in eccentricity by the CB disk and a small change in the orbital period given the small amount of mass lost during S1 (∼0.1  M_⊙).

S2: Disk-dominated phase (epoch mass 1.9>M₁/M_⊙>1.55).

During phase S2, the dust condensation radius moves again past the Roche lobe radius and $R_{\text{dust}}/R_{{\text{L}}_1}$ will remain larger than unity until the end of the simulation. The activation of WRLOF will allow the CB disk to reform (panel e) and produce a strong orbital contraction due to the systemic wind mass loss. Interestingly, the CB disk has a very weak effect on the separation. Mass feeding into the disk increases Σ_in (panel e) and leads to efficient eccentricity pumping (${\dot{e}}_{\text{CB}}$, orange curve in panel f). But, as the star expands on the thermally pulsing asymptotic giant branch (TP-AGB) and the separation keeps shrinking, tides become stronger (${\dot{e}}_{\text{tides}}$, red curve).

Thermal pulses can be identified as spikes separating the interpulse phases in panels c and d displaying to the overfilling factor and the wind mass loss rate, respectively. During the relaxation of the structure following the energy release by the pulse, the star contracts, reducing the efficiency of the tides and the mass loss rate. This allows the CB disk to act almost unconstrained, producing small jumps in e (panel b).

The TP-AGB phase is thus characterized by a decrease in the orbital period (panel a), caused primarily by WRLOF systemic mass loss. The CB disk acts efficiently during this stage, resulting in a growth of the eccentricity (panel b). It is important to remember that in binaries with a shorter initial period, tides will be stronger and circularization may already begin in this phase.

S3: Tidally dominated phase (epoch mass 1.55>M₁/M_⊙>1.15).

Throughout this phase, tides continue to strengthen as the star expands and the orbital period decreases, resulting in a near circularization of the system (e≈0.01 at M₁ = 1.15 M_⊙, panel b). The wind mass loss rate reaches its superwind value of ≈10⁻⁵ M_⊙ yr⁻¹, leading to the accretion of 0.25 M_⊙ onto the companion while in the meantime the primary loses 0.85 M_⊙. These numbers are in line with the ∼30% accretion efficiency predicted by our model (Fig. 1, third and fourth panels). The shrinkage of the separation is mainly due to systemic wind mass loss, the mass transfer onto the secondary providing a small positive contribution to ȧ.

S4: RLOF mass-transfer phase (epoch mass 1.15>M₁/M_⊙>0.8).

The expanding AGB star eventually fills its Roche lobe at M₁≈1.15 M_⊙(panel c) and mass transfer starts. It is important to stress that the companion is more massive than the AGB star (M₂ = 1.25 M_⊙while M₁ = 1.1 M_⊙) and that RLOF begins in an eccentric orbit (e = 0.015, panel b). This unusual situation leads to an overall increase in the separation, despite the tendency of tides to circularize the orbit. The phase-dependent RLOF in the eccentric orbit also produces a modest increase in eccentricity (Eq. (16) and panel f). The abrupt variation in $R_1/R_{L_1}$ and ${\dot{M}}_{\text{RLOF}}$ (panels c and d) at M₁ = 0.85 M_⊙ is due to the development of a thermal pulse and occurs just before the termination of RLOF at M₁≈0.8 M_⊙.

S5: End of the AGB phase (epoch mass 0.8>M₁>0.6 M_⊙).

During this final stage, the AGB star expels its remaining envelope through winds. As the convective envelope becomes thinner, tides become inefficient. The separation now increases due to wind accretion onto the companion, which has become stronger than the orbital shrinkage due to WRLOF. Furthermore, the eccentric RLOF prevented the complete circularization of the binary, enabling efficient resonant interactions with the CB disk in this final stage (panel f). This simulation ends at M₁ = 0.6 M_⊙ with an eccentricity of 0.1 and a period of 4300 d.

As for most stellar evolution codes, convergence difficulties arise at the tip of the AGB (e.g. Wagenhuber & Weiss 1994; Lau et al. 2012), preventing the primary from evolving on the post-AGB and WD cooling tracks. Given the significant decline in the mass loss rate and the short duration of the post-AGB phase, the impact of winds, tides, and RLOF on the orbital parameters is probably limited. However, the CB disk can still act on the eccentricity before dissipating, as will be discussed in Sect. 4.8.

4.3. Evolution in the e−log P diagram

The orbital evolution in the e−log P diagram of Fig. 6 can be understood from our previous analysis of the 8000 d+0.1 system. For systems with v_wind = 6 km s⁻¹and orbital periods shorter than 20 000 d, the activation of WRLOF during the RGB produces a decrease of the orbital period (because η>η_crit). The formation of the CB disk then allows the eccentricity to grow up to core He ignition. During the AGB, WRLOF resumes, the separation decreases and the eccentricity starts rising again until tides take over and start circularizing the system.

Fig. 6. Evolution in the e−log P diagram of the 2.0+1.0 binaries using parameter set AV6. The color of the points refers to the outcome for the system at the end of the simulation. Red represents eccentric final orbits (ef>0), green represents the same, but with the added condition that the binary undergoes RLOF. Orange corresponds to binaries that experienced RLOF, but ended up on a circular orbit. Star symbols show the final orbital parameters. The blue track identifies systems where the AGB radius exceeds the outer Lagrangian point (OLP), at which point, marked by the blue star, the simulation is stopped. The orbital evolution of binaries is illustrated with an evolutionary track, where the thickness of the line increases once the RLOF starts. The solid black line corresponds to orbital parameters such that a 1 M⊙ primary at the tip of the AGB (having a radius of 350 R⊙) fills its Roche lobe, assuming a 1 M⊙ companion.

Eventually, the binaries that come sufficiently close can undergo RLOF, as depicted in Fig. 6 by a thicker curve. As long as the mass ratio q stays above unity, the period and eccentricity decrease, but when q<1, both P_orb and e begin to increase again, producing a loop at the end of the evolutionary track. The final orbital evolution follows three possible paths, depending on whether the final orbit is circular or eccentric and RLOF is stable or not. We consider that RLOF is unstable if the radius of the AGB exceeds its outer Lagrange point (OLP) determined using Eq. (4) of Temmink et al. (2023). In that case, the simulation is stopped and the result discarded.

The first scenario involves binaries with P_i = 16 000 d, the 12 000d+0.01, and 12 000d+0.1 systems (red circles in Fig. 6). These systems end their evolution on eccentric orbits and avoid RLOF. The final period ranges between 6000 and 9000 d and owing to the CB disk, the final eccentricities hit values between 0.2 and 0.3.

The second scenario includes the 12 000d+0.3 binary, systems with P_i = 8000 d, and two with P_i = 6000 d (green triangles in Fig. 6). These binaries undergo stable RLOF in eccentric orbits and because of the phase-dependent mass transfer, the eccentricity is also pumped during RLOF. This process is effective in nearly circularized binaries such as the 8000 d+0.1 system. However, with the reversal of the mass ratio (q<1), the separation will increase again and prevent the system from reaching periods below the black line in Fig. 6, which roughly corresponds to the critical limit for the onset of RLOF.

The final scenario (orange squares) accounts for binaries with initial periods shorter than 6000 d, the exact value depending on the eccentricity. Owing to the greater proximity of the stars, tidal forces quickly circularize the orbit. Consequently, when RLOF begins, the binary system cannot regain some eccentricity, as in the previous scenario. The final period of these circular systems ranges between 2000 and 5000 d.

We stress that at the end of the simulation, the CB disk may still be present and may further impact the orbital evolution. This effect will be assessed in Sect. 4.8.

4.4. Effect of parameter variations

To investigate the effects of the wind velocity, v_wind, and disk parameter, ξ (with values as listed in Table 1), we computed additional models for the 2.0+1.0  M_⊙ system but keeping the same grid in orbital parameters, i.e. P_i={4000, 6000, 8000, 12 000, 16 000} d, and e_i={0.01, 0.1, 0.3}.

4.4.1. Effect of changing the wind speed v_wind

We first focus on models AV6, AV10, and AV15, which use the same value of ξ = 5×10⁻⁴ but differ in the selected wind speed, with v_wind = 6, 10 and 15 km s⁻¹ respectively. The results of these models are displayed in the left column of Fig. 7.

Fig. 7. Same as Fig. 6 but comparing the impact of different wind speed vwind = 6, 10, and 15 km/s and with ξ = 5×10−4 (left column, sets AV6, AV10, AV15) and ξ = 10−4 (right column, sets BV6, BV10 and BV15).

Increasing the wind speed reduces the specific angular momentum carried by the ejected material (η, see Fig. 1) as well as the accretion efficiency β onto the companion. Models AV6 and AV10 present essentially the same evolutionary tracks in the e−log P plane. In the AV10 models, the wind extracts less specific angular momentum but more material escapes from the system (due to the lower accretion efficiency onto the companion). The combination of these two antagonistic effects results in similar angular momentum losses for models AV6 and AV10, and in the end, they lead to comparable evolution.

On the contrary, the AV15 models behave quite differently: only the 6000 d+0.3 and P_i = 4000 d systems experience RLOF and in the remaining models (red track), only the 6000d+0.1 and 8000d+0.3 binaries stay close to their initial periods whereas wider systems expand during their evolution. This can be explained from inspection of Fig. 1 which shows that for v_wind = 15 km s⁻¹ and P_orb≳6000 d, η<η_crit so the systemic wind mass loss will cause the separation to increase. Because AV15 systems remain in long orbital periods, tides have a negligible effect and the resonant interactions with the CB disk can efficiently pump the eccentricity up to 0.3 in systems with initial periods between 6000 d and 8000 d. In systems with larger separations (P_i≥12 000 d), the increase in Σ_in is lower (Eq. (24)) and the effect of the disk weaker, resulting in a wider range of eccentricities, distributed between 0.1 and 0.35.

4.4.2. Effect of surface disk density parameter ξ

The models in the right column of Fig. 7 use the same values of v_wind as their left counterparts, but with ξ = 10⁻⁴ (ξ represents the fraction of the mass lost by the primary that feeds the disk). With decreasing values for ξ, the density at the disk inner rim Σ_in (Eq. (24)) is smaller and the impact of the CB disk (${\dot{e}}_{\text{res}}$) is reduced.

In models BV6, BV10, and BV15, the effects of the disk are negligible and the results are essentially the same as for the same models computed with ξ = 0. In the e−log P diagram, models BV6 and BV10 evolve toward shorter orbital periods at nearly constant e, until tides become efficient and circularization proceeds. In the BV15 models, binaries with P_i≥6000 d expand under wind mass loss and in this case the eccentricity is not reduced. It can even slightly increase due to phase-dependent wind accretion. The surface densities in these models are of the order of 0.1–1 g cm⁻² (compared to 1–10 g cm⁻² in models with ξ = 5×10⁻⁴). This gives an estimate of the lower limit on Σ_in that is needed for the effects of the CB disk to manifest during the ascent of the AGB.

The comparison between models A and B (left column of Fig. 7) shows that the occurrence of RLOF does not greatly depend on the value of ξ. However, due to the higher value of ξ, models A are more prone to start RLOF on an eccentric orbit, leading in the end to a more efficient eccentricity pumping.

4.4.3. Effect of changing the tidal strength

Nie et al. (2017) suggested that a reduction of the tidal circularization rate by up to two orders of magnitude is needed to explain the existence of ellipsoidal red giants in eccentric binary systems. Following Escorza et al. (2020), we consider this effect in models BT01, where the ${\dot{e}}_{\text{tides}}$ terms are reduced by a factor 10 (F_tide = 0.1).

Our results indicate that binaries with P_orb = 16 000 d are too wide for tidal effects to play any role and their tracks are identical to the corresponding systems in model B. For shorter periods, the evolution initially proceeds as in models B, with a contraction of the orbit due to systemic wind mass loss and an increase in the eccentricity resulting from interactions with the CB disk. When mass transfer starts (around 3000 d), the eccentricity is higher than in model AV6 and the e-pumping by RLOF more efficient. When the mass ratio inverts, the orbit widens and the eccentricity can reach 0.6. In conclusion, F_tide has a significant impact on the orbital parameters as it allows to reach higher eccentricities but it does not help in producing shorter period systems.

4.5. The tidally enhanced wind scenario

To account for the observed systems with P_orb≲2000 d, the evolution of the primary must be truncated before its reaches its maximum radius, thereby allowing the binary to shrink to shorter orbital periods before RLOF occurs. We investigate in this section if tidally enhanced wind loss could produce this effect.

Tout & Eggleton (1988) suggested that the presence of the companion star could enhance the primary's wind mass loss rate. They propose a formulation for the wind loss rate, ${\dot{M}}_{\text{wind}}$, given by

$${\dot{M}}_{\text{wind}}={\dot{M}}_{\text{wind}}^{\text{VW}}\times \left\{1+B_{\text{wind}}\times \text{min}\left[{\left(\frac{R_1}{R_{L_1}}\right)}^6,\frac{1}{2^6}\right]\right\},$$(30)

where ${\dot{M}}_{\text{wind}}^{\text{VW}}$ is the AGB wind-loss rate (Vassiliadis & Wood 1993, in our case) and B_wind is a free parameter. The last term in square brackets accounts for a saturation when the star has filled more than one half of its Roche lobe. Tout & Eggleton (1988) found that B_wind≈10⁴ could match the orbital parameters of the RS CVn binary Z Her, a value not dissimilar from that of Siess et al. (2014) (B_wind≈3.6×10⁴) to account for the properties of IP Eri, a long-period eccentric binary with a helium WD.

Figure 8 shows the evolution of the orbital period, eccentricity and primary-star radius for a 2.0 + 1.0  M_⊙, 6000d+0.1 system using the parameters from model AV6, but with B_wind = 10³ (orange) and B_wind = 10⁴ (green). The standard model with B_wind = 0 is shown in blue.

Fig. 8. Orbital period Porb (a), eccentricity (b), primary mass loss rate (c), radius R1 (d), and Roche filling factor $R_1/R_{L_1}$ (e) for the 2.0+1.0 M⊙ system with Pi = 6000 d, ei = 0.1 and different values of Bwind: 0 (dotted blue line), 103 (dashed orange line), and 104 (solid green line).

The tidally enhanced mass loss rate has two effects. First, it accelerates the AGB evolution of the primary and thus reduces the timescale during which tidal circularization and the interaction with the CB disk can operate. As a consequence, the eccentricity remains relatively constant (green and orange curves, panel b) contrary to the system computed without the wind enhancement which nearly circularizes (blue curve). Second, by peeling-off the surface layers more efficiently, the enhanced mass loss limits the expansion of the star which reaches a maximum radius of 260 R_⊙ for B_wind = 10⁴, compared to 350 R_⊙ in the standard (B_wind = 0) case (panel d). Consequently, in the B_wind = 10⁴ system, RLOF begins at a shorter period (∼2200 d) than in the standard case (at P_orb≈2400 d, panel a). During the mass-transfer episode, the B_wind = 10⁴ model experiences a modest increase in P_orb, to 2600 d, contrarily to the other two cases, where P_orb increases to 4000 d. This stems from the fact that in this model the rates of wind mass loss and RLOF mass loss are comparable, of the order of ∼10⁻³ M_⊙ yr⁻¹, so their effects on the orbit compensate. This situation, occurring for the largest value B_wind, is unusual because typically one of the two processes (wind mass loss or RLOF mass transfer) predominates. Therefore, the results for this model should be approached with caution, as the simple modeling applied here may not be physically sound.

The model with B_wind = 10³ (orange curves) has a special evolution in the sense that the expansion of the primary (up to 400 R_⊙, panel d) following the thermal pulse at M₁ = 1.43 M_⊙ leads to a quite large Roche-filling factor ($R_1/R_{L_1}=1.25$, panel e). In this case, the RLOF mass-transfer rate remains larger (5×10⁻³ M_⊙ yr⁻¹) than the wind mass loss rate (10⁻⁴ M_⊙ yr⁻¹, panel c). Therefore, contrarily to the B_wind = 10⁴ case, the orbital period increases to 4000 d after RLOF, similarly to the standard case.

We simulated a grid of binaries using the same parameters as for AV6, but with B_wind = 10³ (model AV6C01) and B_wind = 10⁴ (model AV6C1). The evolutionary tracks and final orbital parameters are shown in Fig. 9. When the tidally enhanced wind is present (panels AV6C01 and AV6C1), binaries with P_orb≲8000 d are more likely to initiate RLOF on eccentric (green triangles) rather than circular (orange squares) orbits. This results from the shortening of the AGB phase which reduces the time during which tidal circularization can act. These binaries initiate RLOF in more eccentric orbits and end up with higher final eccentricities, up to 0.4 and 0.5 for B_wind = 10³ and 10⁴ models respectively, compared to a maximum of ∼0.3 in the AV6 case. Additionally, because the primaries do not grow as big and start RLOF in tighter orbits, they can attain shorter final orbital periods, as low as 3000 d for B_wind = 10³ and 2000 d for B_wind = 10⁴ compared to 4000 d in the non-enhanced case. Remarkably, the region in the e−log P where the AV6C01 and AV6C1 models end their evolution coincides with the lines where a primary with a radius of 350 R_⊙ and 250 R_⊙ fills its Roche lobe, assuming at that time a mass ratio q∼1 (solid black line in Fig. 9). These radii agree well with the maximal radius attained by AV6C01 and AV6C1 models.

Fig. 9. Same as Fig. 6 but for the tidally enhanced wind simulations (parameter sets AV6, AV6C01, and AV6C1) and for the lower-mass system AV6, 1.2 + 0.8 M⊙, bottom right panel). For the massive systems, the three solid lines correspond to orbital parameters where a 250, 350, 450 R⊙ primary fills its Roche lobe, assuming a mass ratio of 1.

4.6. Non-conservative RLOF

During conservative RLOF, mass transfer from a lower-mass donor widens the orbit, preventing the binary from reaching orbital periods significantly shorter than the initial period at which RLOF began. In the case of non-conservative RLOF, if sufficient mass is lost with high specific angular momentum, the separation can keep shrinking even for a mass ratio less than one.

Non-conservative RLOF is modeled using a constant accretion efficiency β_RLOF≤1, such that the gainer accretes at a rate ${\dot{M}}_2=-{\beta }_{\text{RLOF}}\hspace{0.5em}{\dot{M}}_1$. When β_RLOF<1, mass is lost from the system at a rate $\dot{M}=(1-{\beta }_{\text{RLOF}})\hspace{0.5em}{\dot{M}}_1$ and carries away angular momentum at a rate given by Eq. 6. In this formulation we impose η = 0.6, corresponding to the maximum value of η found in the WRLOF simulations of Saladino et al. (2019) for very low wind velocities. We performed simulations for the 2+1 M_⊙ binary using the parameter set AV6, and considered three accretion efficiencies β_RLOF={0.25, 0.5, 0.75}.

Non conservatism only affects the final evolution of the models that cross the solid black line in Fig. 10 which represents the limit below which RLOF is expected to occur during the AGB phase. Models with β_RLOF = 0.75 are quite similar to the conservative ones. The binaries that go through the RLOF phase end up with similar eccentricities but they reach slightly shorter period (between 3500 d and 4200 d compared to 4000 d and 5000 d when β_RLOF = 1). The models with β_RLOF values of 0.5 and 0.25 reach periods as low as 3000 d and 2100 d, respectively. These results show that the degree of conservativeness during RLOF can play a significant role in determining the final orbital parameters. Nevertheless, the observed eccentric barium stars and post-AGB binaries with periods shorter than 2000 d remain unaccounted for.

Fig. 10. Same as Fig. 6 but for non-conservative RLOF. The accretion efficiencies βRLOF are indicated in the upper-left corner of each panel. Conservative mass transfer corresponds to βRLOF = 1 (i.e. AV6). The part of the tracks displayed with thick solid lines refers to ongoing RLOF.

4.7. Influence of the initial binary masses

To assess the impact of a lower-mass binary, we computed a new grid of models using the parameter set AV10 but with M₁ = 1.2 M_⊙ and M₂ = 0.8 M_⊙. Because of its lower mass, the primary reaches a smaller radius at the AGB tip (≈250 R_⊙ compared to 350 R_⊙ for M₁ = 2 M_⊙), implying that RLOF could occur at periods between 2000 d and 3000 d, depending on eccentricity. To consider these cases, we extended the range of initial periods by computing models with P_i= 2000 d and 3000 d.

In the e−log P diagram, the tracks (Fig. 9) follow an evolution similar to the 2.0+1.0 M_⊙ binaries (top panel) and can still be explained in terms of the detailed analysis performed in Sect. 4.2. However, because the evolution of the lower-mass primary is slower and the mass loss rate is smaller, the orbital shrinkage is happening more gradually, leaving more time for the CB disk to act before the tides take over. This translates into a faster increase in the eccentricity with decreasing period (Fig. 9). This also explains why most binaries that fill a significant fraction of their Roche lobe (i.e. P_i≤4000 d, orange squares) circularize before the start of RLOF.

Binaries with P_i≥6000 d do not initiate RLOF and remain on eccentric orbits (red circles). The eccentricity increases due to the CB disk and climbs up to 0.1–0.4, before ${\dot{e}}_{\text{CB}}$ saturates. Because of the lower amount of mass available to feed the CB disk, the surface density Σ_in reaches approximately half the value (1−5 g cm⁻²) of the M₁ = 2.0 M_⊙ model.

4.8. Orbital evolution during the post-AGB phase

Our simulations stop at the end of the AGB phase when less than ≈0.05 M_⊙ is left in the envelope. At that time, mass transfer by RLOF has ceased and tidal interactions have become very weak. However, the surface disk density Σ_in is still high and until the disk dissipates [after ∼120 000 years for P_orb = 2000 d (Eq. (25))], resonant interactions can still operate. To assess the impact of the CB disk during the post-AGB evolution, we resumed the simulations from the last computed model, freezing the structure (and mass) of the primary, letting the companion star evolve normally. The simulations are run for an additional 10⁶ yr until the disk has dissipated and the star stabilized in the e−log P plane. We consider two cases characterized by different CB disk efficiencies (ξ) and in a first step, neglect the effects of tides on the orbital evolution.

The orbital parameters of systems that ended their AGB phase on circular orbits remain almost unaffected during the disk dissipation phase (i.e., superimposed blue dots and stars with e≈0 in Fig. 11). For the AV6 systems finishing their AGB evolution in eccentric orbits (blue dots, top panel), an efficient e-pumping operates during the post-AGB phase leading to eccentricities e≈0.4 at which point the resonant interactions saturate (Eq. (22)). The orbital periods slightly decrease but remain too long to account for the distribution of eccentric barium stars with P_orb<3000 d. In models BV6 (having a reduced ξ), the CB disk has a weaker effect and the final eccentricities distribute more evenly between 0 and 0.35 with little change in the period (bottom panel).

Fig. 11. Orbital parameters of the 2.0+1.2 M⊙ binaries at the end of the simulations, at the AGB tip (circles), for the parameter sets AV6 (top) and BV6 (bottom), and after the dissipation of the CB disk (blue stars).

When the simulations are resumed, the AGB star has a radius of approximately 200 R_⊙, suggesting that tides, at least at the beginning of the post-AGB phase, could still be effective. While our approach does not account for the peel-off and the reduction of the primary radius during the post-AGB phase, we can estimate the maximum effect of tides by re-running the simulations keeping the tidal torques of the primary to its value at the start of the post-AGB phase. These calculations show almost no changes in the final orbital parameters compared to models without tidal torques. So we conclude that during this fast-evolving phase, within our formalism, tides are negligible whereas the CB disk can have a strong impact.

5. Discussion

5.1. Model limitations

5.1.1. WRLOF

Our paradigm rests on the assumption that WRLOF leads to additional extraction of angular momentum from the system and to the formation of a CB disk. The activation of WRLOF depends on the location of the dust-forming radius with respect to the Roche radius but also on the type of dust to be considered. In our simulations, the primary star remains oxygen-rich. This is a consequence of our use of a relatively low-mass primary, of the absence of extra-mixing processes (for example, overshooting at the base of the envelope), and/or because of a potential increase in the mass loss rate due to tidal interactions or RLOF. If the star turns C-rich, the change in condensation temperature will produce a drop in the dust radius by a factor ∼0.36 (Eq. (1)), potentially bringing R_dust within the AGB-star Roche lobe. With the deactivation of WRLOF, the Jeans mass loss would be restored and the orbit would widen upon mass loss. The C enrichment of the envelope by third dredge-up episodes will thus hinder the formation of short-period systems. However, the effect may be limited considering that not all stars become C-rich, and when it happens, it can be late in the evolution, with little time left for orbital changes. Marini et al. (2021) find that a 2 M_⊙, Z = 0.008 star spends 25% of its AGB lifetime as C-rich, this percentage dropping to 8% for a 1.25 M_⊙ star.

The efficiency of WRLOF also depends on the specific angular momentum carried away by the escaping material, as described by our parameter η, which remains uncertain. In our study, we use the prescription of Saladino et al. (2019), derived from SPH simulations of circular systems. As shown later by Saladino & Pols (2019), introducing eccentricity in the simulation reduces η by up to 20% compared to the circular case. These simulations however do not include dust which presence can have a strong impact on the flow, as it can provide additional acceleration, via radiation pressure on the grains. Dust may as well induce heat transfer as gas and dust may not be at the same temperature, which can potentially impact dust formation. Other physical ingredients, such as the equation of state (e.g., Theuns & Jorissen 1993; Malfait et al. 2024) or wind acceleration mechanisms (e.g., Esseldeurs et al. 2023) may potentially impact the parameter η as well.

5.1.2. Disk formation and evolution

This theory of excitation of resonances by the CB disk was originally developed within the context of satellites embedded in a circumplanetary disk. It assumes extreme mass ratios (q=M_S/M_P≪1, where M_P and M_S are the masses of the planet and satellite, respectively), along with the fact that particles within the disk follow circular near-Keplerian orbits, that the disk is axisymmetric and the eccentricity e≪1 (see Goldreich & Tremaine 1982, for a review).

SPH simulations performed by Shi et al. (2012) of a CB disk around an equal-mass binary showed that these assumptions break down. Indeed, they find that the CB disk becomes eccentric, in turn generating an asymmetric feature consisting of a high-density “lump”. Additionally, from the simulations they calculated the total torque acting on the disk via the binary, and compared this to the total torque associated with the OLR at the disk's inner edge, as predicted by the linear perturbation theory (i.e. our ${\dot{J}}_{m\ell }^{+\text{LR}}$). They found that the latter is about a factor of 4 larger than the simulated values. Since our progenitor binaries have q≈1, we may therefore overestimate the rate of eccentricity pumping by the same factor (since ${\dot{e}}_{m\ell }\propto {\dot{J}}_{m\ell }$). More theoretical work is therefore needed to extend the linear perturbation framework to the binary regime.

Simulations (Roedig et al. 2012; Miranda et al. 2017) also show the appearance of accretion streams emanating from the disk's inner edge, which penetrate the binary system and deposit mass and angular momentum onto the stars (see also Artymowicz & Lubow 1994). Shi et al. (2012) found that mass transfer could compensate for the orbital shrinkage produced by the binary's interaction with the CB disk, although the orbit contracted overall. Similar results were found by Roedig et al. (2012) and Valli et al. (2024). Whether the orbit expands or contracts, however, depends on the competition between the second, third and fourth terms in Eq. (2). In neither of these studies did the accretion streams alter the eccentricity, as also verified by Shi et al. (2012) using analytical arguments.

It is hard to quantify how much mass is exchanged between the binary and the CB disk and how it impacts our results. For post-AGB systems, with CB disks of masses ≲10⁻² M_⊙, Oomen et al. (2020) found that the accretion rates from the CB disk to the binary that reproduce observed depletion patterns, are not strong enough to significantly alter the orbit. Although it is beyond the scope of this study to investigate the impact of these streams, it is clear that they can substantially modify the evolution of the orbital separation.

5.1.3. Roche lobe overflow

The occurrence of RLOF mass transfer in eccentric orbits can have significant consequences, because the assumptions of the Roche models break down. Sepinsky et al. (2007) showed that the Roche radius is modified if the star is not synchronized and/or the orbit is eccentric. If the star is rotating sub-synchronously (f_rot≡Ω_*/ω_peri<1, where Ω_* is the star's spin angular speed and ω_peri is the orbital angular speed at periastron), the Roche radius $R_{L_1}$ is larger compared to the classical value given by Eggleton (1983). The opposite is true for a super-synchronously rotating star (f_rot>1). The deformation of the equipotentials opens the Roche lobe around the companion, potentially leading to nonconservative mass transfer through L₂. SPH simulations by Lajoie & Sills (2011) estimated a constant escape rate of 5% from L₂. The escaping matter carries angular momentum, which leads to changes in the orbital response of the binary to the RLOF, and possibly affects the stability of mass transfer.

With our η prescription for the angular-momentum loss from the system, the binary systems can undergo RLOF with initial periods ≲8000 d (for v_wind≈10 km s⁻¹). In their study of blue lurkers and blue stragglers, Sun et al. (2024) implemented wind accretion by WRLOF following the prescription of Abate et al. (2013), but assumed that the matter escapes in the Jeans mode, η=η_iso. As a result of this different treatment of angular momentum loss, they find that RLOF can only be triggered in 2.0 + 1.1 M_⊙ systems with initial periods shorter than about 900 d. The condition for RLOF thus depends critically on the treatment of η.

The RLOF episodes are stable in our simulations. Usually, it is assumed that upon mass loss, convective AGB envelopes expand adiabatically faster than the Roche radius, leading to unstable mass transfer if q>0.8 (Ivanova et al. 2013). Recent studies (e.g. Temmink et al. 2023) showed that the adiabatic response of the donor occurs only above a critical mass-transfer rate. Below this threshold, the superadiabatic surface layers can adjust thermally, resulting in self-regulating mass transfer. Performing the same analysis as in Temmink et al. (2023) for our 2 M_⊙ star at the tip of the AGB, we find that this critical rate is as high as 10⁻¹ M_⊙ yr⁻¹, similar to what was found by Ivanova et al. (2013) and much larger than the mass transfer rate in our simulations. We also checked that the AGB star does not overflow through the outer Lagrangian point and that upon accretion, the companion does not fill its own Roche lobe, thus avoiding evolution to a contact system.

5.2. Comparison to observations

Figure 12 compares the final orbital parameters of our simulated binaries with the observed values for post-AGB stars. There is an almost total incompatibility between the regions of the e−log P diagram reached by our simulations and that populated by post-AGB systems. The situation would be somewhat more favorable if the simulations were compared to the barium dwarf sample (see Fig. 4), because the post-AGB sample lacks long-period systems (be it a real lack or the result of insufficiently long radial-velocity monitoring). However, this comparison was not displayed on Fig. 12 because the observed properties of the two samples look different in Fig. 4, and it is currently unknown whether this difference is the result of an observational selection bias, of the signature of a further (unidentified) physical process operating between the post-AGB and barium-dwarf phases, or a genuine difference hinting at different progenitor populations. Indeed, binary post-AGB stars are often not s-process rich (van Winckel 2003). In any case, the observed systems with P_orb<2000 d cannot be explained by our models, because after mass reversal, RLOF leads to an increase in the orbital period.

Fig. 12. Comparison of the final orbital parameters obtained in our simulations with observed post-AGB orbital elements (red diamonds; see Sect. 4.1 for references of the data). Each panel corresponds to a different model. Green squares correspond to the orbits of the simulated systems at the end of the AGB phase. During the subsequent post-AGB phase, the CB disk continues to act until its dissipation, producing final orbits represented by yellow star symbols.

Our models with tidally enhanced winds (labeled using C01 or C1 suffix in Fig. 12; also Fig. 9) produce the shortest periods because they do not expand as much. Non-conservative evolution (Fig. 10) is able as well to produce tight orbits if a lot of mass is ejected during RLOF (β_RLOF≲0.5), but it is still insufficient to account for the observed distributions. In conclusion, another mechanism must be invoked to account for systems with P_orb<2000 d and the best candidate remains common envelope or grazing-envelope evolution (Kashi & Soker 2018). As revealed by Fig. 12, the sought mechanism should be able to produce systems with periods in the range 100–2000 d. A CB disk, either already present or formed as a result of mass loss during the CE evolution, could subsequently continue to drive eccentricity growth. However, if a system is circularized at the end of the CE phase, reproducing the observed eccentricities of post-AGB binaries would require CB disk masses exceeding those inferred from observations, as noted by Rafikov (2016) and Oomen et al. (2020). This suggests that at least a fraction of post-CE systems must retain some degree of eccentricity.

Observations have revealed the existence of a correlation between the mass of the helium white dwarf and the orbital period for helium white dwarfs and wide sdB binaries formed during the red giant branch phase which has been interpreted as a result of the strong dependence of a red giant's radius on its core mass (Rappaport et al. 1995). To investigate if such a relation exists among Ba stars, as some of these systems may have formed by RLOF, we used the white dwarf masses of barium star binaries with [Fe/H] > –0.5 determined by Escorza & De Rosa (2023). As shown in Fig. 13, there is no clear evidence of a relation between the WD mass and the orbital period. Two factors could explain this observation: (1) the occurrence of thermal pulses during the AGB phase produces a strong increase in the stellar radius (Fig. 2), triggering episodic RLOF potentially followed by a CE evolution that can impact the stellar masses and the orbital period and (2) the CO core mass is not a simple function of the stellar radius as it depends on the efficiency of the third dredge-up (Herwig et al. 1998), which in turn depends on internal mixing processes (potentially altered by the presence of the companion) and on the mass loss rate which is impacted by the envelope enrichment following the third dredge-up episodes. The situation is clearly more complicated when dealing with an AGB star and this may reflect in the (absence of a) white dwarf mass–orbital period relation.

Fig. 13. Orbital period-white dwarf mass relation from a sample of 60 Ba stars analyzed by Escorza & De Rosa (2023). The factor (1−e)3/2 has little impact and is introduced to determine the orbital period if the system had a separation given by the periastron distance where RLOF is most likely to occur. No clear correlation is present unlike in sdB and He white dwarf binary systems.

6. Summary and conclusion

In this study, we explored the effects of WRLOF, CB disk, and tidally enhanced winds primarily on 2.0+1.0 M_⊙ binaries, with initial periods between 4000 and 16 000 d. We found that WRLOF is initiated by the binary systems with initial periods up to 20 000 d as the primary ascends the AGB, however for 16 000 d < P_i<20 000 d, WRLOF occurs only during thermal pulses, which does not alter the orbital evolution. In this regime, winds escaping the systems efficiently extract angular momentum from the binary system and are the main driver for orbital contraction. This effect depends on the AGB wind speed, with lower values leading to faster orbit contraction. For typical AGB velocities (6–15 km s⁻¹), stars initiate RLOF for periods as long as 8000 d. WRLOF may also be initiated during the RGB phase of the primary. However, because of the lower mass-loss rate, its effect on the orbit is limited. By the end of the simulations, the companion has accreted between 0.2 and 0.5 M_⊙ by wind and RLOF.

In our model, the CB disk is fed when WRLOF starts. Its effects are strongest during the TP-AGB phase, where the binary gains most of its eccentricity from the CB disk. During the superwind phase, a higher mass loss rate leaves less time for the CB disk to act.

Binaries may initiate RLOF when they are in an eccentric orbit. The phase-dependent mass transfer produces a significant pumping of eccentricity up to 0.4, even if the orbit is nearly circularized. RLOF is stable for most of our binaries because the mass transfer remains below the threshold value (∼10⁻¹ M_⊙ yr⁻¹) for an effective adiabatic response of the AGB envelope. The few exceptions are for shorter periods and higher eccentricities (typically for the 2.0 + 1.0 M_⊙, P_i = 4000 d, e = 0.3 models). On the other hand, RLOF can also lead to an increase in the orbital period once the mass ratio is inverted. Therefore, systems that initiate (stable) RLOF cluster at the critical period of ∼4000 d, for a wide range of eccentricities.

We explored the effect of tidally enhanced winds on the orbit. The increased mass-transfer rate accelerates the evolution of the AGB star, leaving less time for tidal forces to circularize the orbit. For B_wind≳10⁻³, RLOF occurs in an eccentric orbit and is further pumped during mass transfer with eccentricities reaching e≈0.5. The accelerated evolution results in fewer thermal pulses which in turn limits the pollution of the AGB envelope by heavy elements and potentially causes difficulties in accounting for their large overabundances observed in the companion barium star. We also simulated several systems where tidal forces were weakened by a factor of 10. The majority of these systems did not circularize before RLOF, leading to similarly high final eccentricities. Introducing non-conservative RLOF allowed to reduce the period while maintaining some non-null eccentricity but this requires a low accretion efficiency which may not be realistic.

Finally, we also studied the effect of the disk during the post-AGB phase by freezing the structure of the primary and following the subsequent evolution of the system. We find that the action of the CB disk on the eccentricity is decisive during the post-AGB phase and partially shadows the previous evolution.

Our models successfully reproduce the observed orbital parameters of dwarf barium stars in the e−log P diagram only for systems with P_orb≳4000 d. Specifically, models incorporating tidally enhanced winds or non-conservative RLOF generate eccentric barium stars with periods down to 2000 d, somewhat extending the range of accountable eccentric systems. However, our current models cannot reproduce the eccentric dwarf Ba stars with P_orb≲2000 d, suggesting that these systems underwent unstable mass transfer followed by common-envelope (CE) evolution. In a recent study, Glanz & Perets (2021) showed that common envelope evolution of initially eccentric binaries can indeed produce post-CE binaries with relatively high eccentricities, up to e = 0.4. We could also speculate that some fraction of the ejected CE could form a circumbinary disk, that could then further pump the eccentricity. This remains however to be explored with dedicated simulations.

Acknowledgments

LS is senior FNRS research associate.

References

Appendix A: Variation of ė and ȧ via a circumbinary disk

In this section, we give a more general overview of the resonant interactions than is found in Goldreich & Tremaine (1979) and Artymowicz & Lubow (1994). We present the equations for the torque between the binary and the CB disk and the ensuing variations of ė and ȧ due to Lindblad resonances in the form usually found in the literature. The evaluation of these equations for the resonance (2, 1) gives Eqs. (19), (20) and (21), which correspond to ${\dot{J}}_{\text{res}}$, ȧ_res and ė_res.

In the orbital plane of the binary, centered at the mass center, we expand the time varying potential of the binary Φ(r, θ, t) in double series as

$$\mathrm{\Phi }(r,\theta ,t)={\displaystyle \sum _{m\ell }}{\varphi }_{m\ell }(r)exp[i(m\theta -\ell {\mathrm{\Omega }}_{\text{B}}t)],$$(A.1)

where (r, θ) are the polar coordinates, l is the time-harmonic number, m≥0 is the azimuthal number, ϕ_mℓ(r) is the radially varying potential component for a given (m, ℓ) and Ω_B = 2π/P_orb is the binary's mean orbital angular speed, with P_orb the orbital period. Each component ϕ_mℓ(r) perturbs the disk's structure, giving rise to an m-armed spiral that rotates with a pattern speed

$${\mathrm{\Omega }}_{\text{P}}=\left(\frac{\ell }{m}\right){\mathrm{\Omega }}_{\text{B}},$$(A.2)

i.e. the spiral is stationary with respect to a frame that rotates with angular speed Ω_P. Lindblad resonances (LRs) are localized where the particle of orbital velocity Ω(r) completes one epicyclic orbit of frequency κ(r) between each encounter of the m-armed spiral,

$$m[\mathrm{\Omega }(r)-{\mathrm{\Omega }}_{\text{P}}]=\pm \kappa (r),$$(A.3)

where henceforth, the upper (plus) sign corresponds to the outer LR (OLR), and the lower (minus) sign to the inner LR (ILR). Assuming that the particle is following a near-Keplerian orbit (κ(r)≈Ω(r)), then From Eq. (A.3), we find that the OLR and ILR are located at

$$r_{\pm \text{LR}}={\left(\frac{m\pm 1}{\ell }\right)}^{\frac{2}{3}}a.$$(A.4)

For the (2, 1) OLR, the resonance is at a radial distance r_res = 3^2/3a≈2.08a.

The angular momentum flux at the LRs, ${\dot{J}}_{m\ell }^{\pm \text{LR}}$ was determined by Goldreich & Tremaine (1978, 1979) using linear perturbation theory. The model assumes that (i) the disk is thin (0.01≲ H/r≲0.1, where H is the disk's pressure scale height), (ii) the disk is co-planar with the binary orbit and (iii) the potential of the non-axisymmetric perturbations of the disk's structure is small compared to that of the binary. They find that

$${\dot{J}}_{m\ell }^{\pm \text{LR}}=-m{\pi }^2\mathrm{\Sigma }(r_{\pm \text{LR}})|\mathrm{\Psi }{|}^2{\left(r\frac{\text{d}\mathcal{D}}{\text{d}r}\right)}^{-1},$$(A.5)

where Σ(r_±LR) is the disk's surface density at the LR, $\mathcal{D}={\kappa }^2(r)-m^2[\mathrm{\Omega }(r)-{\mathrm{\Omega }}_{\text{P}}{]}^2$ measures the “distance” from the LR, and Ψ=(Λ_mℓ∓2m)ϕ_mℓ is the forcing function, with Λ_mℓ≡(dlnϕ_ml/dlnr)_LR, also evaluated at the LR (Artymowicz & Lubow 1994). For a Keplerian disk, we have

$$r\frac{\text{d}\mathcal{D}}{\text{d}r}=\mp 3{\mathrm{\Omega }}_{\text{B}}^2\left(\frac{{\ell }^2}{m\pm 1}\right).$$(A.6)

In order to compute the torque for the (2, 1) resonance, Artymowicz & Lubow (1994) use a Taylor expansion for small eccentricities of the potential components

$${\varphi }_{21}(r)=\frac{9{\mathit{GMa}}^2}{4r^3}e\mathcal{M}(1-\mathcal{M}),$$(A.7)

where ℳ=M₁/M. The torque ${\dot{J}}_{\text{res}}$ (Eq. 19) is obtained by evaluating Eq. A.5 for (m, l) = (2, 1).

We derive now the rates of change of a and e due to the resonant interaction with a circumbinary disk. The total orbital energy of the binary system, ℰ, is

$$\mathcal{E}=-\frac{\mathit{GM}\mu }{2a},$$(A.8)

where M=M₁+M₂ is the total mass and μ=M₁M₂/M is the reduced mass. In case of a resonant interaction and in the absence of mass transfer, the rates of change of the semi-major axis, ȧ, and of the orbital energy $\dot{\mathcal{E}}$ are related by

$$\frac{\dot{a}}{a}=-\frac{\dot{\mathcal{E}}}{\mathcal{E}}.$$(A.9)

We can express $\dot{\mathcal{E}}$ in terms of the rate of change of the orbital angular momentum, ${\dot{J}}_{\text{orb}}$ by using the Jacobi constant $\mathcal{C}=\mathcal{E}-{\mathrm{\Omega }}_{\text{P}}{J}_{\text{orb}}$ where Ω_P=(ℓ/m) Ω_B is the pattern speed. Hence,

$$\dot{\mathcal{E}}={\mathrm{\Omega }}_{\text{P}}{\dot{J}}_{\text{orb}}=-{\mathrm{\Omega }}_{\text{P}}{\dot{J}}_{\mathrm{m}\ell }^{\pm \text{LR}},$$(A.10)

where we have made use of the fact that ${\dot{J}}_{\text{orb}}=-{\dot{J}}_{m\ell }^{\pm \text{LR}}$. Therefore, Eq. A.9 writes as

$$\frac{\dot{a}}{a}={\mathrm{\Omega }}_{\text{P}}\frac{{\dot{J}}_{m\ell }^{\pm \text{LR}}}{\mathcal{E}}.$$(A.11)

Substituting Eqs. (A.5), (A.6) and (A.8) into Eq. A.11 yields

$${\left(\frac{\dot{a}}{a}\right)}_{m\ell }=\mp \frac{2{\pi }^2}{3}\frac{\mathrm{\Sigma }(r_{\pm \text{LR}})}{\mu }\left(\frac{m\pm 1}{\ell }\right)\frac{|({\mathrm{\Lambda }}_{m\ell }\mp 2m){\varphi }_{m\ell }{|}^2}{{\mathrm{\Omega }}_{\text{B}}^3{a}^2},$$(A.12)

which, evaluated at the (m, l) = (2, 1) resonance, gives Eq. 20. The eccentricity can be expressed in terms of J_orb and ℰ using

$$e^2-1=\frac{2\mathcal{E}{J}_{\text{orb}}^2}{(\mathit{GM}{)}^2{\mu }^3},$$(A.13)

Taking the time derivative of Eq. A.13 yields

$$\frac{e\dot{e}}{1-e^2}=-\frac{1}{2}\frac{\dot{\mathcal{E}}}{\mathcal{E}}-\frac{{\dot{J}}_{\text{orb}}}{J_{\text{orb}}}=-\frac{1}{2}\frac{\dot{\mathcal{E}}}{\mathcal{E}}+\frac{{\dot{J}}_{m\ell }^{\pm \text{LR}}}{J_{\text{orb}}},$$(A.14)

where J_orb is

$$J_{\text{orb}}=\mu [\mathit{GMa}(1-e^2){]}^{\frac{1}{2}}.$$(A.15)

Inserting Eqs. (A.8), (A.10) and (A.15) into Eq. A.14 gives

$$\frac{e\dot{e}}{1-e^2}=-\frac{{\dot{J}}_{m\ell }^{\pm \text{LR}}}{\mu {\mathrm{\Omega }}_{\text{B}}^{2}a}\left[\frac{\ell }{m}-\frac{1}{(1-e^2{)}^{\frac{1}{2}}}\right].$$(A.16)

Finally, inserting Eq. A.5 into Eq. A.16 and using Eq. A.12 yields

$${\dot{e}}_{m\ell }=\pm \left(\frac{1-e^2}{e}\right)\left(\frac{m}{2\ell }\right)\left[\frac{l}{m}-\frac{1}{(1-e^2{)}^{\frac{1}{2}}}\right]{\left(\frac{\dot{a}}{a}\right)}_{m\ell }.$$(A.17)

Equation (21) is then obtained by evaluating for the (2, 1) resonance.

All Tables

All Figures

Fig. 1. Top two panels: Fraction of the specific angular momentum η carried away by winds during WRLOF as a function of orbital period for a 1+1 M⊙ (top panel) and 2+1 M⊙ (top-middle panel) binary system, and for different wind velocities (vwind, as labeled by the different colors). The horizontal long-dashed line represents the value of η in the case of isotropic wind emission (Jeans mode) and the short-dashed line the value of ηcrit above which the orbit shrinks as a result of wind mass loss. Bottom two panels: Accretion efficiency onto the companion, $\beta =|{\dot{M}}_{\text{acc},2}/{\dot{M}}_1|$ for the same two systems as in the upper panels, either through wind (βBHL, dashed lines) or through WRLOF (βWRLOF, solid lines).

In the text

Fig. 2. Evolution of the initial period $P_0^{\text{WRLOF}}$ for WRLOF as a function of the mass of the primary (Eq. (29)) for a 2.0 + 1.0 M⊙ (solid blue line) and a 1.2+0.8 M⊙ (dashed orange line) binary. The spikes in the curves are associated with the development of thermal pulses. The black dotted lines indicate the critical initial period of circular binaries Pcrit below which WRLOF first starts during the RGB or the AGB phase.

In the text

Fig. 3. Regions in the e−log P diagram where WRLOF is active. If the initial period Pi of a system lies to the left of the curve, a binary undergoes WRLOF during the RGB (orange dashed line), the AGB (solid blue line), or the AGB with C-rich dust (i.e. Tdust = 1500 K, green short-dashed line). Top and bottom panels refer to 1.2 + 0.8  M⊙ and 2.0 + 1.0  M⊙ binaries, respectively. The orange dashed curve and solid blue curve are computed using Tdust = 1000 K (i.e. O-rich dust).

In the text

Fig. 4. Distribution in the e−log P diagram of the observed post-AGB binaries (top), Ba dwarfs (middle), and Ba giants (lower panel). The different lines delineate the orbital parameters for which a primary of radius 250 R⊙ (solid line), 350 R⊙ (dashed line) and 450 R⊙ (dotted line) fills its Roche lobe. These radii correspond to the maximum radii typically attained by solar-metallicity AGB stars of initial masses 1.2, 2 and 3 M⊙, respectively. Both components are assumed to have 1 M⊙ at the time of the RLOF, although the exact masses do not meaningfully affect the lines. As expected, lower primary-mass systems can reach shorter periods without RLOF.

In the text

Fig. 5. Evolution of a 2.0+1.0 M⊙ binary with initial period Pi = 8000 d and eccentricity ei = 0.1 using the parameter set AV6. We present, as a function of the primary mass M1, (a) the orbital period, Porb; (b) the eccentricity e; (c) the Roche-filling factor $R_1/R_{L_1}$; (d) the primary mass loss rates from the wind ${\dot{M}}_{\text{wind}}$ (dashed curve) and RLOF ${\dot{M}}_{\text{RLOF}}$ (solid curve); (e) the surface density at the CB disk's inner rim, Σin; (f) the change in the orbital eccentricity ė (black curve), which is decomposed into contributions from the CB disk (orange), tides (red), RLOF (green). The evolution of the system is divided into different stages, labeled S1–S5, and separated by the vertical black dotted lines.

In the text

Fig. 6. Evolution in the e−log P diagram of the 2.0+1.0 binaries using parameter set AV6. The color of the points refers to the outcome for the system at the end of the simulation. Red represents eccentric final orbits (ef>0), green represents the same, but with the added condition that the binary undergoes RLOF. Orange corresponds to binaries that experienced RLOF, but ended up on a circular orbit. Star symbols show the final orbital parameters. The blue track identifies systems where the AGB radius exceeds the outer Lagrangian point (OLP), at which point, marked by the blue star, the simulation is stopped. The orbital evolution of binaries is illustrated with an evolutionary track, where the thickness of the line increases once the RLOF starts. The solid black line corresponds to orbital parameters such that a 1 M⊙ primary at the tip of the AGB (having a radius of 350 R⊙) fills its Roche lobe, assuming a 1 M⊙ companion.

In the text

	Fig. 7. Same as Fig. 6 but comparing the impact of different wind speed vwind = 6, 10, and 15 km/s and with ξ = 5×10−4 (left column, sets AV6, AV10, AV15) and ξ = 10−4 (right column, sets BV6, BV10 and BV15).
In the text

	Fig. 8. Orbital period Porb (a), eccentricity (b), primary mass loss rate (c), radius R1 (d), and Roche filling factor $R_1/R_{L_1}$ (e) for the 2.0+1.0 M⊙ system with Pi = 6000 d, ei = 0.1 and different values of Bwind: 0 (dotted blue line), 103 (dashed orange line), and 104 (solid green line).
In the text

	Fig. 9. Same as Fig. 6 but for the tidally enhanced wind simulations (parameter sets AV6, AV6C01, and AV6C1) and for the lower-mass system AV6, 1.2 + 0.8 M⊙, bottom right panel). For the massive systems, the three solid lines correspond to orbital parameters where a 250, 350, 450 R⊙ primary fills its Roche lobe, assuming a mass ratio of 1.
In the text

	Fig. 10. Same as Fig. 6 but for non-conservative RLOF. The accretion efficiencies βRLOF are indicated in the upper-left corner of each panel. Conservative mass transfer corresponds to βRLOF = 1 (i.e. AV6). The part of the tracks displayed with thick solid lines refers to ongoing RLOF.
In the text

	Fig. 11. Orbital parameters of the 2.0+1.2 M⊙ binaries at the end of the simulations, at the AGB tip (circles), for the parameter sets AV6 (top) and BV6 (bottom), and after the dissipation of the CB disk (blue stars).
In the text

Fig. 12. Comparison of the final orbital parameters obtained in our simulations with observed post-AGB orbital elements (red diamonds; see Sect. 4.1 for references of the data). Each panel corresponds to a different model. Green squares correspond to the orbits of the simulated systems at the end of the AGB phase. During the subsequent post-AGB phase, the CB disk continues to act until its dissipation, producing final orbits represented by yellow star symbols.

In the text

	Fig. 13. Orbital period-white dwarf mass relation from a sample of 60 Ba stars analyzed by Escorza & De Rosa (2023). The factor (1−e)3/2 has little impact and is introduced to determine the orbital period if the system had a separation given by the periastron distance where RLOF is most likely to occur. No clear correlation is present unlike in sdB and He white dwarf binary systems.
In the text