© ESO 2020

1. Introduction

The evolution of evolved stars in binary systems is still a poorly understood topic in astronomy (for a recent review, see De Marco & Izzard 2017; Boffin & Jones 2019). Consequently, understanding the different possible interactions between binary stars is an active field of research. In particular the interactions occurring during late stages of stellar evolution, when one of the stars has expanded to giant dimensions, still pose many uncertainties. Post-asymptotic giant branch (post-AGB) stars in binary systems are optimal candidates to study these interactions since they left their phase of strong interaction very recently. Previously, we analysed the orbital properties of a sample of 33 Galactic post-AGB binaries (Oomen et al. 2018). The orbital periods are in the range of 100−2500 days, with the majority having non-zero eccentricities even though tides are expected to be strong enough to circularise the orbit in the AGB phase.

Remarkably, the orbital properties of post-AGB binaries show many similarities with other post-interaction systems, such as barium stars (Jorissen et al. 2019; Escorza et al. 2019), carbon-enhanced metal-poor stars (Jorissen et al. 2016; Hansen et al. 2016), S-type symbiotics (Fekel et al. 2000), blue stragglers (Mathieu & Geller 2015), and wide subdwarf B-type (sdB) binaries (Vos et al. 2017). To explain the eccentric orbits in these systems, several eccentricity-pumping mechanisms have been proposed, such as Lindblad resonances in circumbinary discs (Dermine et al. 2013; Vos et al. 2015; Izzard & Jermyn 2018), phase-dependent mass loss (Bonačić Marinović et al. 2008; Kashi & Soker 2018), and white-dwarf kicks (Izzard et al. 2010). However, binary population synthesis models are still unable to reproduce the general orbital properties of these various classes of post-interaction binaries.

The binary interaction between an AGB star and its companion leads to an epoch of strong mass loss. When most of the envelope mass of the donor is lost, the star contracts as a result of mass loss, which sparks the end of the strong interaction with its companion. At this point, the star enters the post-AGB phase, in which the star makes a transition from the AGB to a central star of a planetary nebula (CSPN) on a very short timescale (Van Winckel 2003; Miller Bertolami 2016). The stellar structure consists of a degenerate CO core with a mass larger than 0.51 M_⊙, a He shell, and a low-mass (but extended) H-rich envelope of mass less than 0.02 M_⊙. During the post-AGB phase, this remaining envelope mass is lost and the stellar radius decreases. Because this evolution happens at an almost constant luminosity, the effective temperature of the star increases from about 3000 K to 30 000 K, at which point the central star starts to ionise its surroundings and becomes a CSPN. This evolution continues until the star becomes a hot white dwarf with T_eff ≳ 100 000 K.

A significant fraction of binaries interact on the red giant branch (RGB) instead of the AGB. The evolution of those systems can lead to the formation of the so-called post-RGB binaries (Kamath et al. 2014, 2015, 2016). The properties of those systems are very similar to those of post-AGB binaries, except for a lower luminosity (≲2500 L_⊙), and in some cases different chemical abundances due to their different nucleosynthetic history. Because post-RGB stars avoided core He burning, these stars have degenerate He cores (with M_core <  0.47 M_⊙) and will evolve to eventually become He white dwarfs. Throughout this work, we will refer to the combined population of post-RGB and post-AGB stars as the post-AGB stars, unless specifically stated otherwise.

In post-AGB stars, binarity is related to the presence of a circumbinary disc, since all confirmed post-AGB binaries have a disc-type spectral energy distribution (SED; Waelkens et al. 1991; Oomen et al. 2018). Disc-type post-AGB stars are relatively common and some 80 sources are currently known in our Galaxy (De Ruyter et al. 2006; Gezer et al. 2015; Van Winckel 2019). Discs manifest themselves as a near- or mid-infrared excess, because hot dust is kept close to the central star. Moreover, the formation of such discs can be related to binarity, since mass loss from the AGB star will mainly be focused in the orbital plane in case of a binary. Mass flowing towards the companion leaves the system via the L₂ Lagrangian point, thereby creating a spiral arm (Frankowski & Jorissen 2007; Chen et al. 2017). However, if the companion is massive enough, the spiral arm might remain bound and wrap around the binary, creating a ring as it intersects itself (Shu et al. 1979; Pejcha et al. 2016; MacLeod et al. 2018; Bermúdez-Bustamante et al. 2020).

The presence of discs in post-AGB binaries has also been related to a chemical peculiarity called depletion (Van Winckel et al. 1995; Van Winckel 1997; Maas et al. 2002; Kamath & Van Winckel 2019). Depletion is the systematic underabundance of refractory elements in the photospheres of post-AGB stars, in which the underabundance often scales with the condensation temperature of an element (Maas et al. 2005; Giridhar et al. 2005). The current consensus is that this depletion is caused by the accretion of metal-poor gas from the circumbinary disc (Waters et al. 1992). In this process of dust-gas winnowing, dust is subject to a much stronger radiation pressure than gas, hence does not get accreted. The photosphere of the post-AGB star becomes diluted with gas poor in refractory elements, which leads to the observed abundance patterns (Oomen et al. 2019).

Hydrodynamical models of circumbinary discs show that accretion of gas into the binary cavity is a natural process (MacFadyen & Milosavljević 2008; Shi et al. 2012; D’Orazio et al. 2013). The viscous evolution of the disc leads to accretion streams, taking gas from the inner rim of the circumbinary disc to the circumstellar accretion discs (Mösta et al. 2019; Muñoz et al. 2019). Most studies show that the accretion process in binaries is quite efficient, with a mass flux that is equal to or larger than in the single-star case (Farris et al. 2014; Shi & Krolik 2015). In a previous work (Oomen et al. 2019, which we will refer to as Paper I), we have evolved a grid of post-AGB models while accreting depleted gas at different rates. By comparing this to the observed post-AGB sample, we have shown that accretion rates of the order of 10⁻⁶ M_⊙ yr⁻¹ are required to explain the depletion patterns of the higher-mass post-AGB population.

In this work, we build on the results of Paper I in order to investigate the peculiar eccentric orbits of post-AGB binaries. Instead of modelling the depletion in the general population of disc-type post-AGB stars, we now focus on a more detailed approach for three specific systems. To this end, we make use of the constraints on the disc model from the depletion analysis, which provide us with important information on disc masses, accretion rates, and evolution timescales. We use this to calculate the effect of the circumbinary disc on the binary orbit. The goal is to evaluate whether current theories of disc-binary interactions are capable of reproducing both the depletion and the observed orbital eccentricities, assuming that the progenitor system leaves its strong mass-loss phase in a circular orbit at the start of the post-AGB phase.

We discuss the different observational constraints on our stars in Sect. 2 and we present the data on those stars in Sect. 3. In Sect. 4, we provide the methodology of our stellar evolution and binary evolution calculations, and we introduce the disc-binary interactions used in this work. We present the results in two parts. Firstly, the depletion analysis is given in Sect. 5. Secondly, the results of the orbital modelling are presented in Sect. 6. Finally, we discuss the results in Sects. 7 and 8.

2. Observational constraints

2.1. Luminosities

Luminosity is one of the most important parameters for post-AGB evolution. It determines the evolutionary timescale of the post-AGB phase, since it both influences the wind mass-loss rate and the nuclear burning rate, two processes that directly change the envelope mass of the star. The poorly determined luminosities of post-AGB stars are still one of the main bottlenecks in modelling their evolutionary properties. Below, we discuss two methods by which we can estimate the luminosities of the stars in our sample.

Many post-AGB stars exhibit strong radial pulsations because they are located in the Cepheid instability strip. Post-AGB stars populate the long-period end of the type II Cepheids also known as RV Tauri pulsators (Wallerstein 2002). The radial modes of these stars can be directly related to their physical sizes, which together with a colour or temperature term, allows one to relate the luminosity of the star to its pulsation period. This is known as the period-luminosity-colour (PLC) relation (Alcock et al. 1998; Ripepi et al. 2015).

We use the PLC relation calibrated for RV Tauri stars in the Large Magellanic Cloud (LMC) of Manick et al. (2018) to determine the luminosities of post-AGB stars in our sample. This relation uses the colour-corrected V-band magnitude known as the Wesenheit index (Ngeow & Kanbur 2005) and is given by

$$\begin{array}{c}\hfill M_{\mathrm{bol}}=-3.75log(P_0)+0.55+BC+2.55{(V-I)}_0,\end{array}$$(1)

where P₀ is the observed fundamental pulsation period in days, BC is the bolometric correction in the V band taken from Flower (1996) based on the stellar effective temperature, and (V − I)₀ is the intrinsic colour of the star. The latter is determined from the (dereddened) photospheric model of the SED fit, which is constrained using spectroscopic data from the literature (see Oomen et al. 2018, for more information).

In addition to the pulsation luminosities, we can derive luminosities based on their Gaia parallaxes (Gaia Collaboration 2018) and the integrated photospheric fluxes of the SEDs. To do this, we use distances derived by Bailer-Jones et al. (2018) for the stars in our sample. The luminosities based on this analysis are published in Paper I for all Galactic disc-type post-AGB stars with chemical abundance data available in literature.

However, we consider the luminosity from the pulsation analysis to be more robust, since the second data release of Gaia does not take into account the orbital motion of the binary. This leads to systematic errors on the parallax determination, hence on the distances to our stars. Binaries with orbital periods in the range of 100−700 days are particularly prone to this effect, and often lead to an underestimation of the parallax (Pourbaix 2019).

2.2. Binary properties

Using the luminosity of the post-AGB star, it is possible to estimate its mass from the core mass-luminosity relation (Vassiliadis & Wood 1994). For post-RGB stars, this relation is quite strict since the core mass directly determines the pressure gradient in the H-burning shell. For post-AGB stars, there is more scatter in this relation because the mass of the He shell influences the luminosity. However, it is still possible to get a good estimate for the core mass¹. Once the core mass is known, the total mass is known as well since post-AGB stars have very little mass left in their envelopes (< 0.01 M_⊙). For the stars in our sample, we use a grid of MESA models (see Sect. 4.2) to estimate the core masses based on the observed luminosities.

We have orbital elements available for many post-AGB stars based on a long-term, systematic, radial-velocity monitoring programme with the HERMES spectrograph (Raskin et al. 2011) on the 1.2 m Mercator telescope. This provides us with orbital periods, eccentricities, projected semi-major axes, and mass functions (Oomen et al. 2018; Manick et al. 2019). The inclination angles are not known in general. However, based on the observed properties of the systems, it is possible to define a range of possible inclination angles in order to derive statistical distributions of the orbital properties of the binaries.

We derive several useful system properties, such as the mass of the companion (M₂) and the semi-major axis (a_b). We use a statistical distribution for the inclination angle that represents random orientations in three-dimensional space given by i = arccos(z), where z is the projection of the unit normal vector of the orbital plane towards the observer, which is uniformly distributed in the range of [cos(i_min),cos(i_max)]. We derive the distribution of M₂ by taking 100 000 random values from the distributions of sin i, M₁, and f(m). In a similar way, we used the resulting distribution of M₂ along with the distributions of a₁ sin i, sin i, and M₁ to find the distribution of a_b.

Finally, the Roche-lobe radius is derived using the semi-major axis and mass ratio of the system. The Roche lobe is strictly defined only for circular orbits and synchronised rotation rates. Because we are dealing with eccentric orbits, we evaluate the Roche-lobe radius at periastron using the assumption of circularity and synchronicity. Consequently, we write

$$\begin{array}{c}\hfill R_{\mathrm{L},\mathrm{peri}}=\frac{0.49q^{2/3}}{0.6q^{2/3}+ln(1+q^{1/3})}{a}_{\mathrm{b}}(1-e),\end{array}$$(2)

where R_L, peri is the radius of a sphere with the same volume as the Roche lobe of a binary with separation equal to the periastron distance and q is the mass ratio defined as M₁/M₂ (Eggleton 1983, modified). This is a reasonable approximation that overestimates the Roche-lobe radius by up to 5% (Sepinsky et al. 2007). If the stellar radius remains smaller than the Roche-lobe radius at periastron, the star will not fill its Roche lobe at any orbital phase, even if it is eccentric. We use the median values and 1σ confidence intervals of the distributions of these parameters in our analysis.

2.3. Depletion

In this work, we want to quantify the level of chemical depletion in post-AGB binaries by estimating the dilution of the photosphere by accreted gas. In general, the dilution factor in the photosphere of the post-AGB star can be derived by the ratio of a volatile element, such as Zn or S, to a refractory element, such as Ti or Sc. In plateau depletion patterns, the abundances of moderate to highly refractory elements have similar abundances, such that other elements like Fe can be used together with Ti. In some cases, however, the abundances of volatile elements like Zn or S are significantly lower than those of C, N, and O. This could imply that even S and Zn are attenuated by the depletion process, hence are not suitable to measure the magnitude of depletion (Rao et al. 2012; Gezer et al. 2019). Therefore, we will use the O abundance in those cases to constrain the initial metallicity of the star, because oxygen will not be affected by the depletion process since most of the atoms are locked in CO molecules in gas phase and get accreted by the star.

Furthermore, O will not be affected much by nuclear burning processes and subsequent dredge-ups, as opposed to the N and C abundances. Indeed, it is expected that the envelope composition has been modified by the nuclear evolution of these evolved stars (Karakas & Lattanzio 2014, and references therein). However, none of the confirmed post-AGB binaries show enrichment of s-process elements (with the exception of HD 158616, De Smedt et al. 2016). Additionally, the dust around post-AGB binaries consists of silicates (Gielen et al. 2008, 2011a), although some sources exist with PAH emission and hence a carbon-rich component (Gielen et al. 2011b). The lack of objects with third dredge-up characteristics is possibly due to the large fraction of post-RGB stars in the sample.

To compute the total dilution of the post-AGB star by accreted gas, we will use the average difference between the oxygen abundance on the one hand, and the iron and titanium² abundances on the other hand. The dilution factor of the post-AGB star then becomes

$$\begin{array}{c}\hfill log({\mathrm{X}}_0/\mathrm{X})=[\mathrm{O}/\mathrm{H}]-([\mathrm{Fe}/\mathrm{H}]+[\mathrm{Ti}/\mathrm{H}])/2,\end{array}$$(3)

where [X/H] is the abundance of element X with respect to the solar value in units of dex, and X₀ is its initial photospheric value.

3. Sample stars

We have selected a sample of depleted post-AGB binaries based on several criteria. First, we required that the orbit of the candidate binary has been determined. The second criterion is that the chemical abundances must follow a plateau-type depletion pattern (see Paper I). This allows us to better constrain the dilution factor of the photospheres for our targets. Indeed, for saturated depletion patterns we can only derive a lower limit on the dilution factor, while for plateau-type objects we can restrict successful accretion models to those that match the location of the star in the log(X₀/X)−T_eff plane, as discussed in Sect. 4.2. The aforementioned criteria resulted in a selection of three stars: EP Lyr, RU Cen, and HD 46703. In what follows, we will discuss these individually. The observational properties are given in Table 1.

3.1. EP Lyr

EP Lyr was first classified as an RV Tauri star by Rosino (1951). Manick et al. (2017) analysed the variability of EP Lyr and found a fundamental pulsation period of 41.6 d, and a first overtone at 83.6 d. Based on this pulsation period, we find a luminosity of L_puls = 5150 L_⊙ using Eq. (1). This is in agreement with the luminosity determined from the Gaia parallax and SED fit given by $L_{\mathrm{SED}}={5500}_{-1400}^{+2300}{L}_{\odot }$ (Oomen et al. 2019). This luminosity corresponds to a core mass of 0.56 M_⊙.

Gonzalez et al. (1997) performed a spectroscopic analysis for EP Lyr and found an (average) effective temperature of 6200 K³. Furthermore, the chemical abundances for this star show that it is strongly depleted, and in Paper I we classified this object as having a plateau-type depletion pattern. Remarkably, the abundances of volatile elements such as S (−0.61 dex) and Zn (−0.70 dex) are significantly lower than the C, N, and O abundances (∼solar). Consequently, we derive the dilution factor from Eq. (3), with [Fe/H] = − 1.8 dex and [Ti/H] = − 2.0 dex.

The SED shows a rather weak infrared excess compared to other post-AGB stars (L_IR/L_* ≤ 3%, De Ruyter et al. 2006; Gielen et al. 2011a), suggesting a low dust mass in the disc. The excess peaks at mid-infrared wavelengths, implying that the inner rim of the (dusty) disc is located far from the central binary. The cold, low-mass disc could imply that the disc is highly evolved. Furthermore, despite the oxygen-rich photosphere of the star, the disc shows a mixed chemistry with clear signs of polycyclic aromatic hydrocarbon (PAH) emission (Gielen et al. 2011a).

The binary nature of EP Lyr was first discovered by Gonzalez et al. (1997). The system has a relatively wide orbit of P_orb = 1151 days and is also significantly eccentric (Manick et al. 2017; Oomen et al. 2018). The light curve of EP Lyr does not show any sign of periodic obscuration by a circumbinary disc (known as the RVb phenomenon, Percy 1993; Waelkens & Waters 1993; Pollard et al. 1996). This implies that the inclination angle of the system is less than about 60°. Furthermore, the Hα line shows a double-peaked emission profile, which we interpret as a signature of a circumstellar accretion disc. The Hα line has a strong blue-shifted absorption component at the orbital phase where the companion moves in front of the post-AGB star, at superior conjunction. This is indicative of a jet launched by the companion (Gorlova et al. 2012, 2015; Bollen et al. 2017, 2019). Because the typical opening angle of jets in post-AGB stars is ∼30° and systems with low inclination angles show a permanent absorption by the jet, we assume that the inclination angle of EP Lyr is larger than 30°. The orbital properties of the binary are derived based on the range of possible inclination angles of the system. The median values and 1σ confidence intervals of the resulting distributions are presented in the lower part of Table 1.

3.2. RU Cen

The binary orbit of RU Cen was first determined by Gielen et al. (2007) and was later updated by Oomen et al. (2018). It is a wide binary with an orbital period of 1489 days and is one of the most eccentric post-AGB binaries known with e = 0.62. Its SED shows only a modest IR excess with L_IR/L_* = 39%. Furthermore, the infrared excess peaks at relatively long wavelengths (∼10 μm), which implies relatively cold dust temperatures at the inner rim (Gielen et al. 2007), similar to EP Lyr. Moreover, a mixed disc chemistry is observed with both oxygen- and carbon-rich dust species (Arneson et al. 2017).

RU Cen has been classified as an RV Tauri star by O’Connell (1960). It shows no periodic obscuration by the circumbinary disc as a result of the binary motion (Pollard et al. 1996). There are two dominant pulsation periods: one at 32.30 days and a first overtone at 64.64 days (Maas et al. 2002). This relatively short pulsation period implies a low luminosity from the PLC relation equal to 1900 L_⊙. The luminosity determined from the SED fitting is lower than this value (L_SED = 1100 ± 200 L_⊙). Since the Gaia parallax is unreliable for binaries, we adopt the pulsation luminosity for this object. This corresponds to a post-RGB star with a mass of M₁ ≈ 0.44 M_⊙.

A detailed chemical analysis for RU Cen has been performed by Maas et al. (2002). The effective temperature is given as 6000 ± 250 K, although the formal error of 250 K might be underestimated given the large pulsations RU Cen undergoes. The photosphere shows signs of strong depletion, with [Fe/H] = − 1.9 dex and [Ti/H] = − 1.96 dex. The Zn abundance is also relatively low at −1 dex, and is significantly lower that the C, N, and O abundances at around −0.4 dex. For that reason, we derive the dilution factor as in Eq. (3), which is presented in Table 1.

Similar to EP Lyr, the light curve of RU Cen does not show the RVb phenomenon, suggesting an inclination angle of less than 60°. Because this is a southern target, we do not have a timeseries of optical spectra to investigate the Hα variability. However, the mass function in this system is relatively large. If the inclination angle of RU Cen would be less than 30°, the companion mass would be larger than 7 M_⊙ which is very hard to explain from an evolutionary point of view. Consequently, we apply the same range in inclination angles as EP Lyr. The resulting orbital properties are presented in Table 1.

3.3. HD 46703

HD 46703 was classified as a post-AGB star by Luck & Bond (1984) based on its high luminosity and low surface gravity, and a decade later its binary nature was discovered (Waelkens & Waters 1993). Its orbital period is about 600 days (Oomen et al. 2018), which is significantly smaller than that of EP Lyr and RU Cen. Furthermore, the orbit is significantly eccentric, with e = 0.30.

Hrivnak et al. (2008) performed a detailed photometric study of HD 46703. They classified it as a semi-regular variable, showing (relatively) small-amplitude variability (ΔV ≈ 0.4 mag) with a dominant pulsation period of about 29 days. Assuming that this oscillation is driven by the fundamental radial mode of the star, we derive a luminosity using Eq. (1) equal to 2450 L_⊙. The luminosity based on the Gaia distance is much higher, namely ${7500}_{-2100}^{+3500}{L}_{\odot }$. However, for this object, the Gaia parallax is very likely affected by the orbital motion.

With L_puls = 2450 L_⊙, HD 46703 is on the boundary of being a post-RGB or post-AGB object. The luminosity of this star is quite uncertain given the semi-regular pulsation behaviour and the discrepancy with the value based on Gaia. Therefore, we cannot differentiate between the post-RGB or post-AGB nature of the object. Using solar-metallicity MESA models, we find that the luminosity of 2450 L_⊙ is well-reproduced by a 0.46 M_⊙ core-mass RGB star close to igniting He in its core. Consequently, we take M₁ = 0.46 M_⊙ as the mass of the primary.

The depletion pattern of HD 46703 shows a relatively smooth plateau profile, with similar abundances of Fe and Ti (Hrivnak et al. 2008). The abundances of volatile elements such as S and Zn are about 1 dex higher than those of the refractory elements, but still show an underabundance with respect to the CNO elements. Here, the same argument applies in which the S and Zn abundances could be attenuated by the depletion process.

The SED shows that there is almost no IR excess present for HD 46703, similar to EP Lyr (L_IR/L_* ≤ 3%, De Ruyter et al. 2006; Gielen et al. 2011a). The cold dust temperature implies that there is little to no dust in the disc close to the star. This could mean that the disc is either old such that the inner rim has moved away from the binary, or the dust-to-gas ratio is very small (Gielen et al. 2011b). The former possibility is corroborated by the low activity of the Hα line, which implies very weak ongoing accretion. At superior conjunction, there is an (also weak) increased blue shifted absorption component. Moreover, no RVb signature is observed in the light curve. Consequently, we assume the same range of inclination angles of 30−60°.

4. Methods

4.1. Overview

The methodology in this paper can be divided into two main parts: the depletion analysis and the binary evolution calculations. The depletion analysis is similar to that of Paper I and consists of running a grid of MESA models (Paxton et al. 2011, 2013, 2015, 2018) along the post-AGB phase for each of the three stars in our sample. We include accretion of depleted gas onto the post-AGB star and follow the abundance changes at the stellar surface. In doing so, we assess what accretion rates and disc masses are required to reproduce the observed depletion pattern of the star.

We simulate the evolution of the binary independently from the single-star MESA models in which we reproduce the depletion pattern. For the MESA simulations that successfully describe the depletion pattern, we use information on the evolution of both star and disc as input for the binary evolution calculations. The disc evolution allows us to compute the effect of disc-binary interactions on the binary semi-major axis and the binary eccentricity. Finally, we assess whether at the end of the evolution, the current orbit of each post-AGB star can be adequately reproduced.

4.2. MESA modelling

We used version 10398 of the one-dimensional stellar evolution code MESA to model the onset of chemical depletion in post-AGB stars. We evolved the post-AGB stars as single stars, but included the effect of accretion of metal-poor gas from a circumbinary disc. This method is similar to that previously used in Paper I, hence we refer the reader to Sect. 4.2 of that paper for more information on the physics used in the models.

We prepared the MESA models by first evolving a 1.5 M_⊙ star along the RGB or AGB phase until the star reached the desired He core mass and by subsequently removing the envelope. For the post-RGB star RU Cen, the luminosity suggests a core mass of 0.44 M_⊙ (Sect. 3.2). In this case, we evolved the star along the RGB until the helium core reached 0.44 M_⊙, after which we increased the mass loss rate to 10⁻⁴ M_⊙ yr⁻¹ to remove the envelope. In order to allow for different starting points for post-RGB evolution, we remove the envelope until different amounts of envelope mass are left, or equivalently, until the star reaches different temperatures during the post-RGB phase. This method is repeated for EP Lyr and HD 46703, but with core masses of 0.56 M_⊙ and 0.46 M_⊙, respectively.

To investigate the depletion process, we evolved a set of models with a range of values for the accretion parameters Ṁ₀ and M_d, 0. These provide the time evolution of the accretion rate according to

$$\begin{array}{c}\hfill \dot{M}(t)={\dot{M}}_0{(1+\frac{2{\dot{M}}_{0}t}{M_{\mathrm{d},0}})}^{-3/2},\end{array}$$(4)

which is based on a viscous disc evolution model described in detail in Appendix A. For the initial accretion rate Ṁ₀, we used four different values, namely 10⁻⁷, 10^−6.6, 10^−6.3, and 10⁻⁶ in units of M_⊙ yr⁻¹. This range of accretion rates reproduces the global population of depleted post-AGB stars well (Oomen et al. 2019). For the initial disc mass M_d, 0, we used values of 10⁻³, 10^−2.5, 10⁻², 10^−1.5 in units of M_⊙, which covers the range of observed disc masses (Bujarrabal et al. 2013).

The parameters Ṁ₀ and M_d, 0 are directly related to the viscous evolution of the disc. The disc evolves faster for higher accretion rates paired with lower disc masses, and vice versa. Therefore, the combination of Ṁ₀ and M_d, 0 provide constraints on the viscosity parameter α in the disc evolution model (see Appendix A). This information is used in the evaluation of the disc-binary interactions, discussed below.

As was done in our previous work, we allow the post-AGB star to start its evolution at different temperatures (T₀). This is the temperature of the central star when the accretion from the circumbinary disc starts. This reflects the possibility that the star leaves its binary interaction with different remaining envelope masses, or equivalently, different radii. We expect the stellar radius at the start of the post-AGB phase to be determined by the Roche-lobe size. Consequently, by modelling a range of starting temperatures, we probe a range of initial Roche-lobe sizes and a range of initial semi-major axes. Possible values of T₀ range from the minimum temperature the star can reach on the RGB or AGB (∼3000 K) to its current effective temperature (∼6000 K). We tested different values for T₀ within this range and selected a region of interest for each star in the sample in order to fully probe the parameter space of the accretion model (see Sect. 5).

Next, we select those models that can reproduce the observed depletion pattern of each post-AGB star in the log(X₀/X)−T_eff plane. This in turn allows us to determine the evolutionary timescale of the post-AGB star from the start of accretion until the current point in evolution. Furthermore, the time evolution of the accretion rate (Eq. (4)), and hence of the disc mass, are also known. This information is used to evaluate the interactions between the disc and the binary which we describe in the next subsection.

4.3. Disc-binary interactions

There are several ways in which the disc can interact with the binary. Firstly, the binary can induce density waves at resonant locations in the disc, which carry energy and angular momentum from the binary to the disc (Lubow & Artymowicz 1996; Dermine et al. 2013; Vos et al. 2015). This mechanism decreases the orbital angular momentum of the binary. Secondly, accretion of high-angular momentum gas onto the binary acts to increase the total angular momentum of the binary (Muñoz et al. 2019, 2020; Moody et al. 2019; Duffell et al. 2020). In this case, gas inside the binary cavity exerts a positive gravitational torque on both stars.

4.3.1. Resonance torque

Smoothed-particle hydrodynamics (SPH) simulations by Artymowicz et al. (1991) show that the rotating gravitational potential of a binary can exert a significant torque on the disc by producing spiral density waves at resonant locations in the disc (Goldreich & Tremaine 1979). Its effect on eccentricity in binary evolution models has been investigated in detail by Dermine et al. (2013) and Vos et al. (2015), while Rafikov (2016a) includes the effects of disc evolution in his models. An important characteristic of the resonance torque is that it is independent of the strength and the width of the resonance area, but is completely determined by the global evolution of the disc (Lubow & Artymowicz 1996). In other words, one can write

$$\begin{array}{c}\hfill {\dot{J}}_{\mathrm{res}}=-J_{\mathrm{d}}/{\tau }_{\nu }=-J_{\mathrm{d}}\alpha {\left(\frac{H}{R}\right)}^2\mathrm{\Omega }(R),\end{array}$$(5)

where J_d is the total angular momentum of the circumbinary disc, ${\tau }_{\nu }={(\alpha {\left(\frac{H}{R}\right)}^2\mathrm{\Omega }(R))}^{-1}$ is the (global) viscous evolution timescale of the disc, (H/R) is the aspect ratio of the disc, and Ω(R) is the Keplerian rotation rate at the half-angular momentum radius R. The viscosity parameter α is constrained by the accretion parameters and the binary properties, as is described in Appendix A. Typical values of α found in literature for protoplanetary discs range from 10⁻³ to 10⁻¹ (e.g. Andrews et al. 2009; Rafikov 2017; Ansdell et al. 2018), which largely agrees with values we find from our depletion analysis (see final column in Tables C.2–C.4). In Eq. (5), we use the viscous timescale as defined by Lubow & Artymowicz (1996), since the resonance torque in their SPH simulations agrees well with this definition. We assume a thick disc with aspect ratio (H/R) = 0.2 (Kluska et al. 2018), which develops the high disc luminosities generally observed in post-AGB binaries.

On the basis of SPH simulations, it is shown that at small and moderate eccentricities the outer Lindblad resonance at l = 1 and m = 2 is most important (Artymowicz et al. 1991). This resonance is located at Ω = Ω_b/3 and is responsible for truncating the inner rim of the disc (Lubow & Artymowicz 1996). At small or moderate eccentricities (e ≤ 0.2), we describe the change in eccentricity from the (l, m) = (1, 2) resonance as (Dermine et al. 2013)

$$\begin{array}{c}\hfill {\dot{e}}_{\mathrm{res}}=\frac{l}{m}\frac{2}{{\tau }_{\nu }}\frac{J_{\mathrm{d}}}{J_{\mathrm{orb}}}\frac{1-e^2}{e+\frac{0.01\alpha }{e}}(\frac{1}{\sqrt{1-e^2}}-\frac{l}{m})·\end{array}$$(6)

Here, J_orb is the orbital angular momentum of the binary and τ_ν is as defined in Eq. (5). The change in eccentricity is strongest at low eccentricities and reaches its maximum when e ≃ 0.1α^1/2. However, below this critical point, the (l, m) = (1, 2) resonance can no longer support the inner rim of the disc, and ė rapidly vanishes to zero (Lubow & Artymowicz 1996). At eccentricities larger than this critical point, ė decreases as 1/e. This trend continues as the eccentricity increases and more resonances become excited. When the eccentricity reaches 0.5 − 0.7, these resonances start to cancel, effectively reducing ė to zero. In order to model the change in eccentricity in the regime 0.2 <  e ≤ 0.7, we follow Vos et al. (2015) by using

$$\begin{array}{c}\hfill {\dot{e}}_{\mathrm{res}}=\frac{l}{m}\frac{2}{{\tau }_{\nu }}\frac{J_{\mathrm{d}}}{J_{\mathrm{orb}}}(\frac{c_1}{e}+c_2)\end{array}$$(7)

where c₁ and c₂ are constants that are computed such that Eq. (7) smoothly connects to Eq. (6) at e = 0.2, while going to zero at e = 0.7. When the eccentricity increases beyond 0.7, we set ė_res equal to zero.

4.3.2. Accretion torque

Recent hydrodynamical simulations have shown that gas inside the circumbinary disc cavity exerts a torque on the binary (Miranda et al. 2017). Muñoz et al. (2019) show that the bulk of the torque arises from the extended circumstellar discs around both components of the binary. The overall torque that impacts the binary scales linearly with the mass accretion rate and the specific orbital angular momentum of the binary. Other recent studies performing 2D and 3D simulations of circumbinary discs also show that the binary gains angular momentum from disc interactions (Muñoz et al. 2019, 2020; Moody et al. 2019; Duffell et al. 2020). There is a remarkable agreement between these different works and different system parameters on the proportion of this relation, which is well-described by

$$\begin{array}{c}\hfill {\dot{J}}_{\mathrm{accr}}=0.7\dot{M}{a}_{\mathrm{b}}^2{\mathrm{\Omega }}_{\mathrm{b}},\end{array}$$(8)

where ${\dot{J}}_{\mathrm{accr}}$ is the torque due to accretion, Ṁ is the accretion rate onto the binary, a_b is the binary semi-major axis, and Ω_b is the orbital angular velocity.

Muñoz et al. (2019) also compute the change in eccentricity in their models by evaluating the forces exerted by each cell in the simulation on the components of the binary for four different eccentricities, namely 0.0, 0.1, 0.5, 0.6. Even though the torques are similar in each of these models, the effect on the eccentricity is different. For the circular case, the magnitude of ė is close to zero, while for e = 0.1 the eccentricity increases, and for the large-eccentricity models the eccentricity decreases. This suggests that the sign of ė flips at some point (near e ≈ 0.5). However, since there are only three relevant data points available, we cannot derive a reliable analytic formulation of ė as a result of accretion. In order to see whether this mechanism can produce the eccentric orbits observed in binary post-AGB stars, we will assume that the effect of accretion is independent of eccentricity and has the largest value found in the simulations of Muñoz et al. (2019) (at e = 0.1) which is

$$\begin{array}{c}\hfill {\dot{e}}_{\mathrm{accr}}=2.42\dot{M}{M}_{\mathrm{b}}^{-1}·\end{array}$$(9)

In the models by Muñoz et al. (2019, 2020), Moody et al. (2019), Duffell et al. (2020), the Lindblad resonances as described by Lubow & Artymowicz (1996) do not seem to be excited, or produce a much smaller torque than predicted by the model. Muñoz et al. (2019) show that the torque produced by the circumbinary disc is about a factor of 20 smaller than the torque produced by the gas inside the cavity. However, from an order of magnitude estimation (taking τ_ν ≈ 0.5 M_d/Ṁ, $J_{\text{orb}}\approx M_{\text{b}}{a}_{\text{b}}^2{\mathrm{\Omega }}_{\text{b}}$, and assuming the specific angular momentum of the disc is about twice that of the binary orbit), one finds that the torques produced by resonances in the disc should be of the same order or even larger than those from gas accretion. We interpret this as implying either that the resonances are much weaker in reality than predicted by the analytical model of Lubow & Artymowicz (1996), or that grid-based hydrodynamical codes cannot properly model the resonance torques, as opposed to the more ballistic approach of the SPH code of Artymowicz et al. (1991).

4.4. Binary evolution modelling

The depletion analysis using MESA modelling gives information on the combinations of starting temperature (T₀), initial disc mass (M_d, 0), and initial accretion rate (Ṁ₀) that can explain the depletion pattern. For each of the successful accretion models, we use this information in the binary evolution calculations. These binary evolution calculations are completely decoupled from the single-star MESA models discussed in Sect. 4.2 to maximise computational efficiency. We present the methodology of the binary evolution modelling below.

The first step is selecting an initial orbit for the binary. We always start from an almost circular orbit (e = 0.01). To select the initial semi-major axis a_b, 0 of the binary, we impose two boundary conditions. The first condition is that we require the star to remain inside its Roche lobe at all times in the evolution. Since we treat the depletion analysis in MESA and the binary evolution separately, we cannot include the effect of Roche-lobe overflow in our simulations. This condition is verified by comparing the Roche-lobe radius (Eq. (2)) to the stellar radius of the MESA model at each time step. If a model for a given a_b, 0 experiences Roche-lobe overflow at any point in its evolution, the initial orbit is likely too compact and the model is rejected.

The second condition we impose is that the initial semi-major axis must be small enough such that Roche-lobe overflow (RLOF) could have occurred at some point in its previous evolution. Post-RGB systems can only form if an external agent (i.e. the companion) removes its envelope on the RGB, since stellar winds are not strong enough at this point. We impose this condition by computing the maximum initial semi-major axis based on the maximum possible radius a star can reach for its given core mass while varying the envelope mass. HD 46703 has a core mass of about 0.46 M_⊙. The largest size a star with this core mass can have on the RGB is approximately 0.8 AU. From Eq. (2), the star would fill its Roche lobe in a circular orbit for a_b, 0 ≈ 2.7 AU (using the current masses of both components). Hence we impose a maximum a_b, 0 of 2.7 AU, such that the progenitor of HD 46703 can still have filled its Roche lobe in the past. Applying a similar approach to RU Cen, we find R_max ≈ 0.7 AU and a_max ≈ 3.0 AU, which is larger than for HD 46703 because its companion mass is higher. For EP Lyr, we have have a maximum radius on the AGB of R_max ≈ 1.3 AU and we find a_max ≈ 4.5 AU. Though this restriction is not strictly relevant for the progeny of AGB stars, we keep this assumption for the post-AGB models of EP Lyr since disc formation is often associated with the substantial overfilling of the donor star’s Roche lobe. Regardless, we find that the condition of a_b, 0 ≤ 4.5 AU does not restrict the parameter space of EP Lyr in our results below.

At the start of the binary evolution calculations, we initialise the initial disc parameters discussed in Appendix A. The initial inner and outer disc radii are determined based on the initial semi-major axis of the system. Next, the disc angular momentum is computed by integrating from the inner to the outer radius. Once the initial conditions for both disc and binary have been set, we evolve the system up to the current point in evolution, which is determined by the depletion analysis.

In each time step, we compute the relevant torques acting on the disc and binary and evaluate the changes in angular momentum, semi-major axis, and eccentricity. We use a simple forward Euler method to integrate the parameters (a_b, e, and J_d) over the evolution timescale defined by the depletion modelling. We divide the evolution timescale into equal time steps for the integration of the orbit. Our tests show that 10 000 time steps are necessary and sufficient to ensure a smooth and accurate evolution, while keeping the computational cost relatively low. Using only 10³ time steps gives irregularities for the slowly evolving post-RGB models due to the non-linear behaviour of the resonance and tidal torques, while increasing the number of time steps to 10⁵ produces almost identical results as 10⁴ time steps.

The initial semi-major axis is then optimised using a least-squares minimalisation to find the model for which the final orbit is as close as possible to the observed orbit (in e vs. a_b space). Even though the uncertainties in the observations in Table 1 do not have a Gaussian distribution, the least-squares method will provide a good indication of whether the disc-binary interactions can adequately reproduce the orbit.

Several mechanisms act to change the orbital parameters of the binary. The disc-binary interactions discussed in Sect. 4.3 produce torques which result in changes in eccentricity and semi-major axis. Additionally, an important binary process we consider is spin-orbit coupling due to tides (Zahn 1977). The differential gravitational force of the companion causes distortions in the equilibrium shape of the star. The companion exerts a torque on the post-AGB star due to the dissipation of energy in the tidal bulges, driving the rotational spin to synchronicity with the orbital frequency. Furthermore, the loss of mechanical energy due to tidal friction acts to reduce the eccentricity in most circumstances (Hut 1981). The details on how we handle tides in our models are described in Appendix B.

To evaluate the effect of the disc-binary interactions on the orbit, we investigate several scenarios, presented in Table 2, in which we start the evolution from a circular orbit (e = 0.01) and allow disc-binary interactions to pump the eccentricity. The first case is our standard scenario, in which we include the effects of resonances, accretion, and tides on the binary orbit. This will be discussed intensively in Sect. 6.1. In the next two cases, we simulate the evolution while ignoring the torque from accretion (case 2) and resonances (case 3). This allows us to better assess the effects of these mechanisms on the orbit. Finally, we investigate the possibility that the equilibrium tide model is not adequate for post-AGB stars and that tides are not effective. In this case, we remove the effect of tides on the binary.

Based on these different cases, we compute the total torque acting on the binary at each time step as

$$\begin{array}{c}\hfill {\dot{J}}_{\mathrm{orb}}={\dot{J}}_{\mathrm{accr}}+{\dot{J}}_{\mathrm{res}}+{\dot{J}}_{\mathrm{tides}},\end{array}$$(10)

where ${\dot{J}}_{\mathrm{accr}}$ is given by Eq. (8), and ${\dot{J}}_{\mathrm{res}}$ is given by Eq. (5). The torque originating from tides follows from spin-orbit coupling, hence ${\dot{J}}_{\mathrm{tides}}=-{\dot{J}}_{\mathrm{tides},\mathrm{s}}$ with ${\dot{J}}_{\mathrm{tides},\mathrm{s}}$ given by Eq. (B.4)⁴. For the disc, we can write ${\dot{J}}_{\mathrm{d}}=-{\dot{J}}_{\mathrm{res}}-{\dot{J}}_{\mathrm{accr}}$, which are the negative contributions of the disc-binary terms in Eq. (10). For each case, we only take into account the relevant terms as described in Table 2. We note that the torques contributing to Eq. (10) are averaged over the orbit, hence it is unnecessary to treat these equations as phase dependent.

Similarly, we can write the total effect on the eccentricity as

$$\begin{array}{c}\hfill \dot{e}={\dot{e}}_{\mathrm{accr}}+{\dot{e}}_{\mathrm{res}}+{\dot{e}}_{\mathrm{tides}},\end{array}$$(11)

where ė_accr is as computed from Eq. (9) and ė_res is as given in Sect. 4.3.1. The effect of tides is given in Appendix B and provides either a positive or negative contribution to ė, depending on the rotational velocity of the star.

Finally, we use Eqs. (10) and (11) to compute the change in semi-major axis by using

$$\begin{array}{c}\hfill {\dot{a}}_{\mathrm{b}}=2a_{\mathrm{b}}(\frac{{\dot{J}}_{\mathrm{orb}}}{J_{\mathrm{orb}}}+\frac{e\dot{e}}{(1-e^2)})·\end{array}$$(12)

Here we assume that the masses of both stars are conserved. This is a fair assumption, since M_b ≫ M_d.

5. Results: depletion modelling

5.1. HD 46703

We selected a range of starting temperatures (T₀) to cover the parameter space of the accretion model as well as possible. For HD 46703, this results in four different starting temperatures: 3250 K, 3500 K, 3750 K, and 4500 K. At T₀ <  3250 K, the models evolve too slowly and do not become depleted even for the highest disc masses and accretion rates, while at T₀ >  4500 K, our lowest accretion rates are still too high and the models are depleted too quickly. Models for different T₀ are depicted in four separate panels in Fig. 1. Because HD 46703 has a plateau-type depletion pattern, successful models are required to go through the observed location in the log(X₀/X)−T_eff plane. All of the models that do cross the 1σ boundary covered by the observations are used in the analysis of the disc binary interactions below and are presented in Table C.1.

Fig. 1. Depletion versus Teff plot for HD 46703 with starting temperatures of 3250 K, 3500 K, 3750 K, and 4500 K from top to bottom, respectively. The black error bar denotes the observed value of HD 46703. The legend is given at the top of the figure, in which different colours and line styles represent different initial accretion rates and initial disc masses, respectively. Ṁ0 is expressed in M⊙ yr−1 and Md, 0 is given in M⊙.

In Fig. 1, there is a clear anti-correlation between T₀ and initial disc mass. This trend is explained by the fact that the stellar atmosphere can be diluted most rapidly at high T_eff values, because at that point the outer convective envelope mass becomes very low. Consequently, what matters most is the rate of accretion at those high T_eff. In slowly evolving post-RGB stars, like HD 46703 and RU Cen, the evolution timescale is long, such that the main parameter determining the accretion rate is the initial disc mass (see Fig. 5 in Oomen et al. 2019).

We see that for the lowest T₀ of 3250 K (upper panel of Fig. 1), only the highest disc-mass models are capable of reproducing the depletion pattern of HD 46703. The star still has a substantial H-rich envelope at 3250 K of ∼0.03 M_⊙, hence it takes over 100 000 years before it reaches its current point in evolution. In order for the disc to still provide sufficient accretion at later stages in the evolution, the initial disc mass should be high.

5.2. RU Cen

For RU Cen, the full parameter space of the accretion model is well represented by models with T₀ values of 3400 K, 3700 K, 4000 K, and 4500 K (see Fig. 2 from top to bottom, respectively). Similar to HD 46703, RU Cen is a slowly evolving post-RGB star of low luminosity. Consequently, we observe the same trend in which lower values of T₀ require higher initial disc masses to become depleted. The collection of models that fit the observed location of RU Cen are presented in Table C.1.

Fig. 2. Depletion versus Teff plot for RU Cen with starting temperatures of 3400 K, 3700 K, 4000 K, and 4500 K from top to bottom, respectively. See legend at the top for definitions of colours and line styles.

5.3. EP Lyr

For EP Lyr, we selected three values for T₀: 3500 K, 4000 K, and 4500 K. These models are depicted in three separate panels in Fig. 3. Increasing T₀ beyond 4500 K produces similar results as those of T₀ = 4500 K. We therefore restrict our parameter space to the given range.

Fig. 3. Depletion versus Teff plots for EP Lyr with starting temperatures of 3500 K, 4000 K, and 4500 K from top to bottom, respectively. See legend at the top for definitions of colours and line styles.

For EP Lyr, the set of models that fit the observed parameters in Fig. 3 are more heterogeneous in terms of disc mass compared to the other two stars in our sample. This illustrates the degeneracy in the parameter space between Ṁ₀ and M_d, 0. However, we find that in none of the cases, models with low accretion rates of Ṁ₀ = 10⁻⁷ M_⊙ yr⁻¹ are able to reproduce the observed depletion pattern of EP Lyr. This is caused by the fast evolution of a 0.56 M_⊙ post-AGB star, as was shown in Paper I.

6. Results: orbital modelling

6.1. Case 1

In case 1, we investigate the combined effect of the accretion torque and the resonance torque on the evolution of the binary. The accretion torque is exerted by gas inside the disc cavity and leads to an increase in orbital angular momentum, while the resonance torque removes orbital angular momentum by launching spiral density waves at the inner rim of the circumbinary disc. As discussed in Sect. 4.3, both mechanisms pump the binary eccentricity and change the semi-major axis. Furthermore, we also take into account the tidal interaction, which tends to circularise the binary.

Figure 4 shows the evolution of the binary orbit in the e vs. a_b plane for each of the stars in the sample. In each of the panels, we show the best-fitting models (by least-squares analysis) for all the accretion models that reproduce the observed depletion value in Figs. 1–3. The models start from a circular orbit (e = 0.01) and gradually become eccentric due to the interaction with the disc. However, none of the models can reproduce the orbital properties of the stars in our sample. The initial and final orbital parameters for all the models are shown in Tables C.2–C.4.

Fig. 4. Best-fitting models in the e vs. ab plane for case 1 with both accretion and resonance torques. The black error bars show the location for each of the stars in the sample. The line styles used in this plot are the same as for the models that fit the observations in Figs. 1–3. The symbols at the end points represent models of different T0. The models start from a circular orbit with e = 0.01, which is shown by the thin black line. The grey dashed curves are lines of constant RL, in which the given temperature corresponds to the radius for this star’s luminosity such that R* = RL. Models for HD 46703 are shown in the top panel, for RU Cen in the middle panel, and for EP Lyr in the bottom panel.

The top panel of Fig. 4 shows the case-1 models for HD 46703. Several trends can be observed here. The low-T₀ models start and end in wider orbits. This is in the first place because the lower T₀ models are limited by their Roche lobes. The grey dashed curves in Fig. 4 show lines of constant R_L, such that models above those lines will fill their Roche lobes at that given temperature. As a star evolves, it becomes hotter and is able to evolve above those lines without filling its Roche lobe. Additionally, models with lower T_eff have a larger R_* and are subject to much stronger tidal forces, since ė_tides ∝ (R_∗/a_b)⁸. Consequently, these models need to start in wider orbits so as to reduce the tidal interaction. For the models with T₀ = 3250 K in the top panel of Fig. 4, the effect of tides is particularly clear, as the eccentricity decreases at later times in the evolution.

To evaluate the effect of the different mechanisms acting on the eccentricity, we show the evolution of ė for a model of HD 46703 with parameters T₀ = 3250 K, Ṁ = 10⁻⁶ M_⊙ yr⁻¹, and M_d, 0 = 10^−1.5 M_⊙ (solid orange line in top panel of Fig. 4) in the top panel of Fig. 5. The two most important effects on the orbital eccentricity are tides and resonances. Accretion only has a minute effect on the eccentricity. This can be understood from the small eccentricity pumping effect in Eq. (9). Since the total disc mass is in general less than 1% of the total binary mass, the maximum pumping of eccentricity will be a ∼0.02 increase.

Fig. 5. Time evolution of ė (top panel) and $\dot{J}$ (bottom panel) for the different mechanisms used in case 1 for HD 46703. The parameters used in this model are Md, 0 = 10−1.5 M⊙ and Ṁ0 = 10−6.0 M⊙ yr−1 with T0 = 3250 K, which corresponds to the solid orange line in the top panel of Fig. 4. The solid blue line represents the sum of the three other effects, i.e. resonances, accretion, and tides. The two most important effects on the orbital evolution are resonances and tides.

Resonances are most effective at low eccentricities, which results in a large ė_bin early in the evolution. However, as the eccentricity increases, tides become stronger since ė_tides scales linearly with eccentricity to first order. This effect is countered, however, by the contraction of the star. Since ė_tides ∝ (R_∗/a_b)⁸, the tides reduce in strength again as the star evolves, and almost vanish when the effective temperature of the star increases beyond 5000 K. This allows e to increase again slightly.

The evolution of the torque exerted on the binary system is shown in the bottom panel of Fig. 5. Overall, the resonance torque is about five times stronger than the accretion torque throughout the evolution of the disc. The net result is that the binary loses angular momentum, resulting in the decrease of the semi-major axis for HD 46703 in Fig. 4. Furthermore, it is clear from Fig. 5 that the spin-orbit coupling due to tides only weakly changes J_orb, but provides a positive contribution since the star rotates super-synchronously during its evolution (see Fig. B.1).

We show the case-1 models for RU Cen in the middle panel of Fig. 4. The maximum initial semi-major axis for RU Cen to have experienced Roche-lobe overflow is a_max = 3.0 AU, as described in Sect. 4.4. This semi-major axis is smaller than the currently observed semi-major axis of 3.8 AU (albeit in an eccentric orbit). This causes the most optimal models for RU Cen to start at the maximum a_b, 0. The final semi-major axes of the models are smaller than the observed value and the eccentricities reach values slightly higher than 0.2.

To reproduce the orbit of RU Cen, the semi-major axis should increase, while the eccentricity needs to be pumped substantially to its currently observed value of 0.62 ± 0.07. While the resonance torque is the most effective interaction to pump the eccentricity, this can only lead to a decrease of the orbital size. It is therefore impossible for this mechanism to produce the current orbital properties of RU Cen. Moreover, the large eccentricity of 0.62 is hard to explain for any orbital size using resonances, since these will start to cancel when the eccentricity reaches ∼0.5 and become very ineffective (Artymowicz et al. 1991).

The bottom panel of Fig. 4 displays the case-1 models of EP Lyr. The high luminosity of EP Lyr results in a faster evolution for this object, with timescales in the range of 1000−10 000 yr (see Table C.1). The disc-binary interactions do not have enough time to pump the eccentricity up to values above 0.15. Models with T₀ = 4000 K and T₀ = 4500 K have a final semi-major axis close to the observed value, while models with T₀ = 3500 K start in wider orbits to avoid Roche-lobe overflow throughout their evolution.

To better understand the effects of the individual disc-binary interaction mechanisms, we show the evolution of the models with only resonance torque (case 2) and only accretion torque (case 3) in the next subsection. We focus on HD 46703 since this star shows the greatest diversity in possible evolution scenarios. Finally, we will discuss models in which we ignore the equilibrium tide (case 4) in Sect. 6.3.

6.2. Case 2 and case 3

In case 2, we analyse the effect on the orbit when we only include the resonance torque described in Sect. 4.3.1, along with the tidal interaction. Lindblad resonances drive spiral density waves at the inner rim of the disc, which take away orbital energy and angular momentum from the binary. Figure 6 shows the evolution of the models for HD 46703 in the e vs. a_b plane. The model results are presented in Table C.2.

Fig. 6. Best-fitting models of HD 46703 in the e vs. ab plane for the effect of resonances in case 2. The symbols and line styles are the same as in the upper panel of Fig. 4.

It is clear that the results for case 2 are very similar to those of case 1. As was shown in Fig. 5, the resonance torque has a much larger impact on the orbit than the accretion torque. The Lindblad resonances are effective at pumping eccentricities, as was already shown by Dermine et al. (2013), Vos et al. (2015). However, for the case of HD 46703, the eccentricities only increase to ∼0.15, which is far below the observed value of 0.30 ± 0.02. Moreover, the final semi-major axis is too large for many of the low-T₀ models. Consequently, we cannot simultaneously reproduce the depletion pattern and the orbit of HD 46703.

In case 3, we investigate the effect of the accretion torque presented in Sect. 4.3.2. This exerts a positive torque on the binary including an eccentricity pumping effect. The orbital evolution is shown for the case of HD 46703 in Fig. 7 and presented in Table C.2. Accretion results in modest widening of the orbit for the low-T₀ models. All of the models remain circular, which confirms the small effect of accretion on the eccentricity, as shown in Fig. 5. We note that we assumed the maximum eccentricity increase found in the models of Muñoz et al. (2019); the actual effect is likely to be even smaller.

Fig. 7. Same as Fig. 6, but for the accretion torque described in Sect. 4.3.2.

6.3. Case 4

Case 4 explores the scenario in which tides are ignored, as a way of simulating the possibility that tides are not as strong as theory predicts. The equilibrium tide theory is well tested for the case of a red-giant star with a large, extended envelope. Even though post-AGB stars have an extended red-giant like structure, the envelope has a very low mass and density and convective damping may be less effective than for a normal red giant. Moreover, post-AGB envelopes are in a state of non-homologous contraction. In Fig. 8, we show the results of models with resonance and accretion torques, but without the effect of tides.

Fig. 8. Same as Fig. 4, but for models in which the effect of tides is ignored (case 4).

In the case of HD 46703, we observe two main differences in Fig. 8 with respect to the models that include tides in Fig. 4. Firstly, the models reach higher eccentricities, with the low-T₀ models coming close to the observed eccentricity. Secondly, we see a general shift of the models to lower a_b, 0 compared to case 1. The latter is due to the strong dependence of tides on the semi-major axis. Models for case 1 starting at larger semi-major axes are able to become more eccentric, hence provide a better fit to the observed value. The model with T₀ = 3250 K in the top panel of Fig. 8 reaches e ≈ 0.25, which is close to the observed value of 0.30. However, the final semi-major axis is much larger than the observed value. The other two models of T₀ = 3250 K with higher accretion rate (solid orange and red lines in Fig. 4) end up filling their Roche lobes for case 4, even when starting at a_max = 2.7 AU. The general result of case 4 for HD 46703 is that the models cannot simultaneously explain both orbital size and eccentricity.

For RU Cen and EP Lyr, the results of case 4 in Fig. 8 are similar to those of case 1. This is because tides play a less important role in the orbital evolution of these systems. RU Cen is still limited by the maximum initial semi-major axis of the system. For EP Lyr, the evolution timescale is too short for the eccentricity to be pumped to significant values by the disc-binary interactions. The parameters for all case-4 models are presented in Tables C.2–C.4.

6.4. Assumptions

In this work, we have made several assumptions within our models. Here, we investigate the effect of some of these assumptions on the results derived above.

6.4.1. Initial outer disc radius

In our disc model (Appendix A), we assumed that the initial outer disc radius R_out, 0 is situated at 10a_b in order to calculate the initial disc angular momentum. However, since both R_out, 0 and the initial disc angular momentum are unknown physical quantities, we investigate its effect on our results by changing its value for several models.

We arbitrarily chose to test this assumption for models of HD 46703 with T₀ = 3750 K. We recalculate the models for case 1, but change the value of R_out, 0 to 5a_b and 50a_b. The results are given at the bottom of Table C.2. We find that the final results are similar to those of the standard model. Consequently, the value of the initial outer disc radius has little impact on the orbital evolution, and our assumed value of 10a_b is reasonable.

The effect of this parameter can be understood from its interaction with the resonance torque, which has the largest effect on the binary properties in our models. Equation (5) gives the strength of the torque exerted by resonances as the ratio of J_d/τ_ν. Increasing the value of R_out, 0 increases the initial disc angular momentum via Eq. (A.5). Even though this also results in a larger viscosity parameter (Eq. (A.2)), the overall viscous timescale increases since the half-angular momentum radius scales linearly with R_out, 0 which significantly reduces Ω(R). Consequently, the overall effect of increasing the outer radius is a small reduction in the resonance torque, which does not significantly impact our results. On the other hand, changing R_out, 0 does result in changes in the viscosity parameter α. Since neither variable is directly constrained by observations for our objects, we cannot make conclusive statements on the values of these parameters based on our models.

6.4.2. Disc scale height

In our binary evolution models, we assumed a thick disc model with a large aspect ratio (H/R) = 0.2. This value is based on the large scale height inferred from post-AGB discs (Kluska et al. 2018). Furthermore, when comparing the viscous timescale from the resonance model of Lubow & Artymowicz (1996) to that of our disc model in Eq. (A.2), we find that these timescales are similar when taking (H/R) ≈ 0.2. As already noted in the previous subsection, the choice of R_out, 0 also influences the initial viscous timescale.

In the resonance model of Lubow & Artymowicz (1996), the global viscous timescale of the disc is related to the scale height via the unknown viscosity parameter α. By constraining the aspect ratio in our models, we therefore parametrise the viscous timescale via the viscosity parameter. Since we fix the viscosity parameter by our depletion analysis in MESA, one can see that a larger value for the aspect ratio represents a shorter viscous timescale of the disc in the resonance model, hence results in a stronger torque and more eccentricity pumping.

In an attempt to quantify the influence of the disc aspect ratio on our results, we evaluated the effect of a smaller value of this ratio, namely (H/R) = 0.1, which might be more representative for objects with a smaller infrared excess, such as EP Lyr and HD 46703. We present the results for models of HD 46703 with T₀ = 3750 K at the bottom of Table C.2. The effect is quite clear: a lower aspect ratio results in less eccentricity pumping. The final eccentricity of these models is about half of that of our standard models.

6.4.3. Stellar wind strength

Stellar winds in post-AGB stars still pose a major uncertainty in their evolution (Cranmer & Saar 2011; Miller Bertolami 2016). In Paper I we investigated the effect of stellar winds on the formation of depletion patterns; more precisely, we investigated its effect on the accretion properties required to reproduce the observed depletion patterns. We showed that weaker winds result in a much slower evolution of post-AGB stars. Therefore, higher disc masses are required for the star to become depleted.

In our standard model, we use a semi-empirical Reimers-type wind prescription as given by Eq. (4) of Schröder & Cuntz (2005) for the wind mass-loss rate in the post-AGB phase. To evaluate the effect of different mass-loss rates on our results, we produced models with a weaker stellar wind for the case of HD 46703. We reduced the wind mass-loss rate by a factor of 100 and performed the depletion modelling and orbital evolution modelling as before.

Figure 9 shows the MESA models for HD 46703 in the log(X₀/X)−T_eff plane for models with T₀ = 3750 K in the weak-wind limit. Comparing with the standard model in the bottom panel of Fig. 1, the initial disc mass needs to be higher in order to explain the depletion pattern of HD 46703. This is due to the much longer evolutionary timescale of the star in the weak-wind regime, as can be seen in Table C.1.

Fig. 9. Same as the bottom panel of Fig. 1 (with T0 = 3750 K), but for models with a stellar wind which is a factor of 100 weaker. In this case, more massive discs are required to explain the observed depletion pattern due to the longer evolution timescale.

For the successful models, we plot the binary evolution for case 1 in Fig. 10. Comparing these to the standard wind models, we find that the models require much larger initial semi-major axes. For the model with M_d, 0 = 10^−1.5 M_⊙ (solid orange line), we see the large effect of tides at reducing the eccentricity. Because of accretion, the model gained envelope mass and expanded, such that T_eff dropped to 3500 K. Furthermore, because of the longer evolution timescale, tides are much stronger since the star remains large for a very long time. Consequently, all models are limited by the tidal interaction before the eccentricity reaches the observed value.

Fig. 10. Eccentricity versus semi-major axis plot for case-1 models with T0 = 3750 K. The four transparent models on the left have standard wind parameters (bottom panel of Fig. 4), while the five models on the right are the successful models of Fig. 9 with a 100 times weaker wind.

To summarise, the effect of a weaker wind on our results is a much slower evolution, which means higher initial disc masses are required to produce our depletion pattern. The latter results in stronger disc-binary interactions, while the former gives more time to pump the eccentricity. Despite this, we find that the models are still unable to explain the orbit for the case of HD 46703. As the eccentricity increases, the equilibrium tide becomes stronger and starts to counteract the pumping mechanisms. We conclude that, even though stellar winds play an important role in our results, the observed orbits of post-AGB binaries still cannot be explained.

7. Discussion

The goal of this work is to investigate whether disc-binary interactions can produce the eccentric orbits in post-AGB stars and other post-interaction binaries. Figures 4–8 show that none of the evaluated interaction mechanisms provide a satisfactory solution for the orbits of EP Lyr, RU Cen, and HD 46703. In particular the accretion torque as given in Muñoz et al. (2019) does not significantly impact the eccentricity throughout post-AGB evolution. The effect of tides and resonances (if present) do play a relatively important role, but can only pump the eccentricity up to ∼0.2.

7.1. Orbital discrepancy between observations and models

A key factor in the failure to reproduce the observed orbits is the evolution timescale of the stars. Models that start at higher temperatures evolve very rapidly to their current point in evolution. The eccentricity-pumping mechanisms therefore have insufficient time to produce the observed orbital eccentricities. Low-T₀ models generally reach higher eccentricities. These models evolve slowly, because they start their evolution with a high-mass envelope. They require high disc masses to become depleted, which increases the strength of the interactions. For these models, tides play an important role, since the stellar radius remains large for a long time. We find that as the eccentricity of those models increases, the torques evolve to an equilibrium situation in which ė_tides ≈ −ė_res. Because the resonances start to cancel for increasing eccentricity while tides become much stronger, we expect this equilibrium to arise for all orbits as e reaches ∼0.2, unless the star spins up fast enough for tidal circularisation to become ineffective (Fig. B.1), or unless the stellar radius has sufficiently decreased such that the tidal interaction becomes weaker.

Another important issue for explaining the orbit arises with the requirement for the star to remain within its Roche lobe throughout post-AGB evolution. Given their current luminosities and orbits, the stars in our sample will fill their Roche lobes at temperatures below about 4500 K. Post-AGB stars spend a long time at low T_eff, but evolve quickly once the temperature increases beyond 4500 K. For that reason, the orbit does not change much when the star becomes hotter than 4500 K. In order to explain the orbit, one needs to pump the eccentricity during the slower stage in evolution, that is at low T_eff but large R_*. However, this conflicts with the requirement to avoid RLOF.

In nature, the stars could in principle evolve with an additional phase of RLOF. Even though we cannot include this in our current models, we expect this to change the models in two ways. Firstly, RLOF would manifest itself as an extra mass-loss agent for the post-AGB star, which speeds up the evolution of the star. Because the post-AGB phase is a phase of contraction and the companion is higher mass, mass transfer will always be stable. Secondly, the orbit will also change as a result of mass transfer. Phase-dependent RLOF could also have an eccentricity pumping effect (Soker 2000; Bonačić Marinović et al. 2008), but this would only have a small contribution because of the tiny amount of mass left in the envelope that can be transferred.

A further conundrum concerns the orbits of wide post-RGB systems like RU Cen. Post-RGB stars must have formed via mass transfer, since stellar winds alone are not expected to be strong enough to remove the envelope on the RGB. The present semi-major axis of RU Cen is too large to accommodate a Roche-lobe filling RGB star in a circular orbit. Therefore, we imposed a maximum initial semi-major axis of 3 AU. However, the disc-binary interactions we investigate primarily result in a decrease of a_b, making it impossible to reproduce the orbit in this way.

This result suggests that this system already left the preceding phase of strong interaction in an eccentric orbit. We show such models for the case of RU Cen in Fig. 11, where we allow for a non-zero initial eccentricity. The low-T₀ models are still somewhat limited by their Roche lobes. The current Roche-lobe size of RU Cen at periastron passage suggests that the star detached with T_eff ≈ 4500 K. Models that have T₀ = 4500 K evolve on fast enough timescales such that they are not affected by tidal evolution or disc-binary interactions, hence their orbits remain essentially unchanged. Consequently, Fig. 11 shows that the orbit of RU Cen, as well as its depletion pattern, can be reproduced if the initial orbit is allowed to be eccentric. This formation scenario requires the eccentricity of the binary to have been set by the unknown binary interaction mechanism that led to the expulsion of the progenitor’s envelope, while slightly widening the orbit. This can potentially be accomplished by mass loss and RLOF near periastron (Bonačić Marinović et al. 2008), especially when the mass transfer rate becomes very large and the mass-loss timescale becomes much shorter than the circularisation timescale (e.g. Kashi & Soker 2018).

Fig. 11. Eccentricity versus semi-major axis plot for case 1 models of RU Cen. These models are similar to the middle panel of Fig. 4, except the initial eccentricity is allowed to be non-zero. The star must remain inside its Roche lobe during its evolution. Starting points are given by a small circle, the end points are shown as open symbols. For T0 = 4500 K the orbital parameters hardly change and coincide with the observed point.

Finally, other eccentricity pumping mechanisms during the post-AGB phase could be at play that we have not considered. A notable possibility could be dynamical interactions with a tertiary component (Toonen et al. 2020), in which the Lidov-Kozai effect can induce eccentricity in the inner and outer orbits (see Naoz 2016, and references therein). Even though it is unlikely that the eccentric orbits of the whole population of post-AGB binaries are the result of triple interactions, it could provide a solution for a fraction of the systems, including special cases such as RU Cen.

7.2. Model uncertainties

It is important to note that the luminosity has a large impact on the results. The luminosity plays a dominant role in the evolution timescale of the star, since it determines both the nuclear burning rate and the stellar wind strength. Higher-luminosity stars have a much shorter evolution timescale, which is a key parameter in our models to pump the eccentricity. Moreover, luminosity determines the radius of a star at a given temperature. Lower-luminosity stars are therefore less constrained by their Roche-lobe sizes during their orbital evolution and experience a weaker tidal interaction.

In this work, we used an empirical PLC relation to estimate the luminosity. The pulsation period is generally well-determined and is related to the radius of the star. For RV Tauri pulsators, there is a lot of intrinsic scatter in the PLC relation (Ripepi et al. 2015). The large pulsation periods also hamper the spectroscopic determination of T_eff, which can vary over the orbital cycle by up to 1500 K. The changing colour of the star therefore induces additional uncertainty in deriving the luminosity from the PLC relation. Especially for HD 46703, the uncertainty is large due to the semi-regular nature of the pulsation period. It would therefore be preferable to use parallax-based luminosities.

Currently, the second data release of Gaia does not take into account the effect of binary orbits on the parallax determination, hence this is too unreliable to use in this work. In particular orbits with periods in the range of 100−700 days are sensitive to systematic errors (Pourbaix 2019). We look forward to the second part of Gaia DR3, which will improve on the distances (and luminosities) significantly. This will allow for better constraints on the stellar model in addition to those of the PLC relation.

In addition to the luminosity, there are several uncertainties in the depletion analysis that propagate to our models. The values we derived for the dilution factor are sensitive to the O abundance, for which there are few good lines to use for the abundance determination. It is therefore difficult to trace back the abundances in the stellar atmosphere before the depletion process initiated. However, it can be seen in Figs. 1–3 that the accretion models are not very sensitive to the dilution factor, since the dilution of the envelope in our models increases steeply around 6000 K. A small difference in observed dilution factor can therefore be modelled with similar accretion parameters. An accurate determination of T_eff for these stars is however important.

For the modelling of depletion patterns in MESA, the stellar wind strength plays an important role. It has a large impact on the evolution timescale, hence on our results as was shown in Sect. 6.4. Furthermore, the amount of mixing in the envelope of the star changes the rate at which the atmosphere becomes diluted, as discussed in Paper I. However, this effect is not expected to be important for temperatures below about 6000 K.

The disc model also introduces uncertainties in our results. In this work, we used a very simple model for the viscous evolution of the disc. As explained in Appendix A, our disc model assumes evolution without a central torque. However, we know that the disc-binary interactions investigated in this work exchange angular momentum between disc and binary. This impacts the temporal behaviour of the accretion rate, which might decrease more slowly with time than given by Eq. (4) (see Appendix A).

Furthermore, other disc evolution mechanisms could be at play other than viscosity. In particular disc mass loss can be effective in post-AGB binaries. A potentially important effect is photoevaporation, which is expected to be one of the leading mechanisms for the dispersal of protoplanetary discs (Alexander et al. 2014; Owen & Kollmeier 2019). High-energy photons, primarily produced by the accretion of gas onto the central stars, ionise the disc stripping it of gas and dust. This mechanism could be especially important for post-RGB models with longer evolutionary timescales, since accretion is expected to dominate early on in protoplanetary disc evolution (Alexander et al. 2014). Other possible mechanisms that could accelerate disc evolution are magnetohydrodynamical disc winds (Blandford & Payne 1982) and the formation of second-generation planets (Turner et al. 2014), although post-AGB discs are presumed to be stable against gravitational instabilities (Kluska et al. 2018).

7.3. Additional model constraints

For the complete sample of post-AGB binaries, it turns out that many orbits have an eccentricity in the range of 0.2−0.3 (see Fig. 5 of Oomen et al. 2018). This is only slightly larger than the eccentricities we find in our models. It is therefore possible that some post-AGB orbits can be explained with disc-binary interactions, assuming that the resonance torque is indeed as strong as is claimed by Lubow & Artymowicz (1996). It is still a concern that the strength of this torque was not reproduced in more recent works by for example Muñoz et al. (2019).

In future work, we will expand our sample to other post-AGB binaries with saturated depletion patterns. For these objects, we only have lower limits on the required accretion rates and disc masses required to reproduce the depletion (see Paper I for an extensive discussion). Although this removes important constraints on the disc model and post-AGB evolution, such an analysis will provide additional insights in the general reproducibility of the orbits of post-AGB binaries by disc-binary interactions. Moreover, the sample of saturated-type post-AGB binaries is much larger and more diverse.

Observed disc properties could provide important additional constraints on the accretion models derived from the depletion analysis in MESA. For now, there is much degeneracy in the parameter space of Ṁ₀, M_d, 0, and T₀. Deriving disc parameters, such as current disc mass, inner rim radius/temperature, and outer disc radius, could partly alleviate these degeneracies. However, this requires a detailed interferometric study of these systems. To probe the inner rim dust structure, near- or mid-infrared observations are required (e.g. Kluska et al. 2018; Ertel et al. 2019), while the outer gas structure can be probed at mm-wavelengths (e.g. Bujarrabal et al. 2018). Kluska et al. (2019) obtained near-IR data with PIONIER on the VLTI for RU Cen, but could not detect the dust inner rim which is probably too cold. For EP Lyr and HD 46703, currently no interferometric data is available.

The disc masses of post-AGB binaries are very uncertain, even for well-studied systems such as IRAS 08544−4431. From a radiative transfer model of the SED and interferometric data, it is possible to obtain a rough estimate of the total dust mass in the disc (Kluska et al. 2018). However, the total disc mass is still subject to the uncertain gas-to-dust ratio of the disc. On the other hand, mm-wavelength studies can trace the total mass of the gas component at larger radii (Bujarrabal et al. 2018). This assumes that the CO transition lines are optically thin, which is not always the case. Regardless, based on the small IR excesses in the SED of the objects in our sample, we estimate the current disc masses to be smaller than ∼10⁻³ M_⊙ for EP Lyr and HD 46703. This favours models with initial disc masses of 10⁻³ M_⊙, which reach the current point in evolution with disc masses in the range of 3 − 8 × 10⁻⁴ M_⊙. However, these models reach relatively low eccentricities < 0.05. Another possibility is that the discs start out more massive, but lose mass due to extra disc mass-loss mechanisms that we did not take into account.

In addition to the IR excess, a time series of optical spectra can give some information on the accretion properties. A jet launched by the companion as a result of its accretion produces a P-Cygni profile in the Hα line when the companion moves in front of the post-AGB star, at its superior conjunction. Based on the jet absorption, it is possible to derive the accretion rate onto the companion (Bollen et al. 2019).

8. Conclusions

We modelled the effect of accretion during the post-AGB phase with the MESA code to constrain the evolution of the star and circumbinary disc by reproducing the observed chemical depletion pattern. Using these constraints, we investigated several disc-binary interaction processes. We analysed the orbital evolution by including the effects of gas accretion inside the binary cavity, Lindblad resonances, and the equilibrium tide. We summarise our findings below.

We find that none of the models are able to simultaneously reproduce both depletion and orbit for the three stars in our sample (i.e. EP Lyr, RU Cen, and HD 46703). The different disc-binary interaction mechanisms we investigated cannot pump the eccentricity to the observed value within the evolutionary timescales of the stars. The torques that are produced by accretion streams and accretion discs inside the cavity of the circumbinary disc have very little effect on the binary eccentricity, and cause a small expansion of the orbit. The torque originating from the outer (l, m) = (1, 2) Lindblad resonance can efficiently pump the eccentricity to about 0.2, but quickly loses its potency at higher eccentricities. At that point, the equilibrium tide becomes of similar strength, inhibiting a further increase in e.

In addition to tidal interaction, the requirement to avoid RLOF is the main constraining factor. As the binary becomes more eccentric, the Roche-lobe radius becomes smaller around periastron, quickly leading to RLOF. In order to avoid RLOF during post-AGB evolution, the initial semi-major axis needs to be large enough such that the post-AGB star remains inside its Roche lobe at all times. Consequently, most models in which the eccentricity pumping starts at low T_eff and large R_* end up in too wide orbits compared to the observed semi-major axis.

Potentially successful models therefore need to start their post-AGB evolution at higher T_eff. Since the stellar radius is smaller, the model can avoid RLOF and is less affected by tidal interaction. However, the evolution timescale at high T_eff becomes increasingly short, and the disc-binary interaction mechanisms do not have enough time to pump the eccentricity. Consequently, the general trend in our results is that models that start accretion at low starting temperatures (T₀) become more eccentric but have too wide orbits since they had to avoid RLOF, while models that start at high T₀ have a semi-major axis close to the observed value, but end with low eccentricity.

Within this framework, the orbits of wide post-RGB systems, such as RU Cen, cannot be reproduced. Since post-RGB systems must be the result of mass transfer, during the interaction on the RGB the orbit must be small enough for RLOF and is expected to be circular when the star enters the post-RGB phase. Since the strongest disc-binary interactions (Lindblad resonances) have a tendency to decrease the semi-major axis, this can never lead to the wide, eccentric orbit of RU Cen.

Another possible solution to explain the orbits of these systems is that they leave their previous binary interaction phase in an eccentric orbit. If the stars in our sample became detached while in their current orbit, the evolution occurs too fast for tides to change the orbit, while accretion can still produce the observed depletion pattern. However, this would require the system to have gained its eccentricity during the previous phase of strong mass loss which led to the expulsion of the envelope.

We investigated the impact of several assumptions on our results and found that these do not impact our results in such a way that the models provide a better fit to the observed orbits. There could be other factors limiting our results, such as the poorly known stellar luminosities or the simple disc model used in this work. We expect the former to improve with the third data release of Gaia. The latter could still be improved by including a more detailed prescription of the effect of torques on disc evolution. Furthermore, effects such as photo-evaporation could become relevant for the long-lived post-RGB discs, although this would be an additional disc mass loss mechanism and result in weaker disc-binary interactions. Finally, we intend to expand our analysis to a larger sample of post-AGB stars in the future.

Acknowledgments

GMO acknowledges support of the Research Foundation – Flanders under contract G075916N and under grant number V434818N. HvW acknowledges support from the Research Council of the KU Leuven under grant number C14/17/082. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. This research has made use of NASA’s Astrophysics Data System Bibliographic Services. The authors would like to thank the anonymous referee for the constructive comments that have improved this paper. Furthermore, the authors thank Rob Izzard and Christoffel Waelkens for their additional feedback.

References

Appendix A: Disc model

The disc is described by a viscous evolution model. It can be shown that if the evolution of the disc is dominated by internal viscosity, a self-similar solution can be found for the evolution of the disc (Pringle 1991; Rafikov 2016b). This allows us to derive the time evolution of important disc properties, such as disc mass and accretion rate, based on their initial values.

In this work, we focus on the simple case in which viscosity is independent of the local surface density. A well-known solution occurs when there is no torque at the centre of the disc, such that the accretion rate can be described by (Rafikov 2016a)

$$\begin{array}{c}\hfill \dot{M}(t)=\frac{M_{\mathrm{d},0}}{2t_0}{(1+\frac{t}{t_0})}^{-\frac{3}{2}},\end{array}$$(A.1)

where M_d, 0 is the initial disc mass, and t₀ is the initial viscous timescale given by

$$\begin{array}{c}\hfill t_0=\frac{4}{3}\frac{\mu }{k_{\mathrm{B}}}\frac{a_{\mathrm{b}}}{\alpha }{\left(\frac{4\pi \sigma {(G{M}_{\mathrm{b}})}^2}{\zeta L_{\ast }}\right)}^{1/4}{\left(\frac{\eta }{I_L}\right)}^2,\end{array}$$(A.2)

with μ the mean atomic weight, k_B the Boltzmann constant, a_b the semi-major axis, α the viscosity parameter, σ the constant of Stefan-Boltzmann, M_b the total mass of the binary, L_* the luminosity of the post-AGB star, and ζ a constant that describes the geometrical effect of disc irradiation (∼0.1). The two final parameters in Eq. (A.2) are related to the angular momentum of the disc, with η giving the specific angular momentum of the disc in units of specific binary angular momentum, and I_L is a constant that describes the distribution of angular momentum (≈1.8, see Fig. 5 of Rafikov 2016b, for λ = 0.5).

The accretion rate given by Eq. (A.1) decreases over time as a result of the decreasing density in the disc. This is caused by two factors. Firstly, the total mass in the disc decreases as matter flows into the binary cavity. Secondly, as the disc evolves, the outer radius increases over time, thereby decreasing the overall density in the disc. Without a central torque, this evolution results in a time dependence of Ṁ ∝ t^−3/2. However, in this work, we aim to investigate the case of a non-zero central torque, which impacts the time-dependent behaviour of the disc (Rafikov 2016b,a). Recent 2D hydrodynamical models of finite discs suggest the accretion rate evolves as Ṁ ∝ t^−4/3 (Muñoz et al. 2020). However, in this work we will assume the simple form of accretion evolution with the formulation as given in Eq. (A.1).

In the model of Rafikov (2016b), the presence of a non-zero central torque requires the inner surface density distribution to scale with r^−3/2. Even though observations of post-AGB stars find that a double power-law distribution better fits the density distribution (Hillen et al. 2015; Kluska et al. 2018), these observations are based on the dust distribution, which might be different from the distribution of gas. Consequently, we will keep the simple power-law distribution for the surface density with Σ = Dr^−3/2 for the sake of consistency with our model.

The proportionality constant D for the surface density profile is computed by requiring that integrating the surface density over the whole disc gives the total disc mass:

$$\begin{array}{c}\hfill M_{\mathrm{d}}={\displaystyle {\int }_{R_{\mathrm{in}}}^{R_{\mathrm{out}}}\mathrm{\Sigma }(r)2\pi r\mathrm{d}r.}\end{array}$$(A.3)

To evaluate the integral in Eq. (A.3), we need to determine the inner and outer boundaries of the disc. The size of the cavity becomes larger when the eccentricity increases (Mösta et al. 2019) and depends on the structure of the disc (Artymowicz & Lubow 1994). The inner radius of the disc can be expressed as (Artymowicz & Lubow 1994; Dermine et al. 2013)

$$\begin{array}{c}\hfill R_{\mathrm{in}}=a_{\mathrm{b}}[1.7+\frac{3}{8}\sqrt{e}log\left({\alpha }^{-1}{\left(\frac{H}{R}\right)}^{-2}\right)],\end{array}$$(A.4)

where (H/R) is the aspect ratio of the disc which we assume to be equal to 0.2. The outer radius of the disc R_out increases over time as a result of the viscous evolution of the disc. The initial value of the outer disc radius is a free parameter and will determine the initial angular momentum of the disc J_d, 0 since

$$\begin{array}{c}\hfill J_{\mathrm{d}}={\displaystyle {\int }_{R_{\mathrm{in}}}^{R_{\mathrm{out}}}\mathrm{\Sigma }(r)\mathrm{\Omega }(r)2\pi r^3\mathrm{d}r.}\end{array}$$(A.5)

Here, we take R_out, 0 = 10a_b, though this initial value is arbitrary. As the disc evolves, the outer radius can be derived by solving Eqs. (A.3) and (A.5) for R_out(t). The integrals are evaluated numerically by dividing the radial coordinate into 1000 parts and applying the trapezium rule.

Similar to Paper I, we rewrite Eq. (A.1) as a function of the initial accretion rate Ṁ₀ and initial disc mass M_d, 0 by substituting t₀ = Ṁ_d,0/(2Ṁ₀), such that

$$\begin{array}{c}\hfill \dot{M}(t)={\dot{M}}_0{(1+\frac{2{\dot{M}}_{0}t}{M_{\mathrm{d},0}})}^{-3/2}.\end{array}$$(A.6)

In Eq. (A.6), Ṁ(t) is defined as the rate at which mass flows from the disc into the binary cavity. Since mass ratios in post-AGB binaries are in general larger than 0.3 (Oomen et al. 2018), both components of the binary system accrete roughly equal amounts of mass (Farris et al. 2014; Muñoz et al. 2020; Duffell et al. 2020). Consequently, we assume in our MESA models that the post-AGB star accretes half of the mass entering the cavity, or Ṁ_PAGB = 0.5 × Ṁ_bin.

Appendix B: Equilibrium tide

We include the effect of tides on the orbital evolution to properly investigate the evolution of the binary (Zahn 1977). We follow Hut (1981) for an (exact) formulation of the equilibrium tide. These tides are caused by the deformation of the star under the gravitational influence of the companion. This produces a tidal bulge that lags behind or precedes the motion of the secondary in the case of a rotational spin that is not synchronised to the orbit. More precisely, the bulge is offset with respect to the line connecting the centres of both stars, which produces a gravitational torque on the post-AGB star.

This simple model of the equilibrium tide has been derived by Hut (1981) and is given by

$$\begin{array}{cc}\hfill \dot{e}=& -27{\left(\frac{k}{T}\right)}_{\mathrm{c}}{q}^{-1}(1+q^{-1}){\left(\frac{R_{\ast }}{a_{\mathrm{b}}}\right)}^8\frac{e}{{(1-e^2)}^{13/2}}\hfill \\ \hfill & \times [f_3(e^2)-\frac{11}{18}{(1-e^2)}^{3/2}{f}_4(e^2)\frac{{\mathrm{\Omega }}_{\mathrm{spin}}}{{\mathrm{\Omega }}_{\mathrm{orb}}}],\hfill \end{array}$$(B.1)

where f₃ and f₄ are defined as in Hut (1981):

$$\begin{array}{cc}\hfill f_3(e^2)& =1+\frac{15}{4}{e}^2+\frac{15}{8}{e}^4+\frac{5}{64}{e}^6\hfill \end{array}$$(B.2)

$$\begin{array}{cc}\hfill f_4(e^2)& =1+\frac{3}{2}{e}^2+\frac{1}{8}{e}^4.\hfill \end{array}$$(B.3)

Here, q is the mass ratio defined as M₁/M₂ and the ratio ${\left(\frac{k}{T}\right)}_{\mathrm{c}}$ is computed using Eq. (30) in Hurley et al. (2002).

The final term in Eq. (B.1) includes the ratio of the spin frequency over the orbital frequency, which determines the position of the tidal bulge with respect to the line of centres. If ${\mathrm{\Omega }}_{\mathrm{spin}}>\frac{18}{11}{\mathrm{\Omega }}_{\mathrm{orb}}$ (in the quasi-circular limit), the torque exerted by the bulge leads to a pumping of the eccentricity, while if ${\mathrm{\Omega }}_{\mathrm{spin}}<\frac{18}{11}{\mathrm{\Omega }}_{\mathrm{orb}}$, the eccentricity tends to decrease.

For post-AGB stars, we expect the progenitor system to exit its phase of strong interaction in a synchronised fashion. Synchronicity does not necessarily persist during post-AGB evolution, hence we evaluate the change of rotation rate of the star as it evolves. Several processes act to change the rotation rate of the star. In what follows, we assume that the envelope is rigidly rotating and completely decoupled from the degenerate core.

First of all, the star shrinks and the envelope mass decreases as a result of stellar evolution. This causes the moment of inertia of the envelope to become smaller and the star to rotate faster. We compute the moment of inertia I of the envelope at each time step in our MESA model by taking the sum $\sum \frac{2}{3}{m}_i{r}_i^2$ from the base of the envelope to the surface, where m_i is the mass of cell i and r_i is the distance of that cell to the centre of the star.

In addition to changes to the structure, the star can gain angular momentum due to accretion of gas from the circumbinary disc. We assume that the accreted gas has the specific angular momentum of Keplerian rotation at the surface of the star. Since gas will be accreted at the equator, we can write ${\dot{J}}_{\mathrm{accr},\mathrm{s}}=\dot{M}\sqrt{G{M}_{\ast }{R}_{\ast }}$, where M_* and R_* are the mass and radius of the post-AGB star. The value of Ṁ(t) is known at each time step from the disc model, while R_* is given by our MESA model.

Even though accretion is taking place, stellar winds and nuclear burning cause the mass of the envelope to decrease over time. Stellar winds also take away angular momentum at the surface of the star. This is computed as ${\dot{J}}_{\mathrm{wind},\mathrm{s}}=\frac{2}{3}{\dot{M}}_{\mathrm{wind}}{R}_{\ast }^2{\mathrm{\Omega }}_{\mathrm{spin}}$, where Ṁ_wind is the mass lost from stellar winds and is a negative quantity. We assume that the wind is spherically symmetric, hence the factor of 2/3.

Finally, tides drive the post-AGB star towards co-rotation with the orbit. Similar to the eccentricity evolution of Eq. (B.1), we apply the formalism derived by Hut (1981) such that ${\dot{J}}_{\mathrm{tides},\mathrm{s}}$ is given by

$$\begin{array}{cc}\hfill {\dot{J}}_{\mathrm{tides},\mathrm{s}}=& 3{\left(\frac{k}{T}\right)}_{\mathrm{c}}\sqrt{G(M_1+M_2)a_{\mathrm{b}}}\frac{M_2^2}{M_1}{\left(\frac{R_{\ast }}{a_{\mathrm{b}}}\right)}^8\frac{1}{{(1-e^2)}^6}\hfill \\ \hfill & \times [f_2(e^2)-{(1-e^2)}^{3/2}{f}_5(e^2)\frac{{\mathrm{\Omega }}_{\mathrm{spin}}}{{\mathrm{\Omega }}_{\mathrm{orb}}}],\hfill \end{array}$$(B.4)

where f₂(e²) and f₅(e²) are defined as

$$\begin{array}{cc}\hfill f_2(e^2)& =1+\frac{15}{2}{e}^2+\frac{45}{8}{e}^4+\frac{5}{16}{e}^6\hfill \end{array}$$(B.5)

$$\begin{array}{cc}\hfill f_5(e^2)& =1+3e^2+\frac{3}{8}{e}^4.\hfill \end{array}$$(B.6)

Given these four processes, we compute the effect on the rotation of the post-AGB star with the initial condition of co-rotation with the orbit. The initial spin angular momentum of the envelope is given by J_spin, 0 = I₀Ω_orb, 0. The different torques acting on the envelope give

$$\begin{array}{c}\hfill {\dot{J}}_{\mathrm{spin}}={\dot{J}}_{\mathrm{accr},\mathrm{s}}+{\dot{J}}_{\mathrm{wind},\mathrm{s}}+{\dot{J}}_{\mathrm{tides},\mathrm{s}}.\end{array}$$(B.7)

The contraction of the star does not exert a torque, but changes the moment of inertia over time. At each time step, the new spin angular momentum and moment of inertia of the star determine the rotation rate as Ω_spin = J_spin/I. Because the angular frequency of rotation changes, the effect of tides change as well (Eqs. (B.1) and (B.4)). We show as an example the evolution of the rotation rate for a model of HD 46703 with T₀ = 3500 K in Fig. B.1.

Fig. B.1. Spin evolution of HD 46703 for a model with T0 = 3500 K, Ṁ0 = 10−6.3 M⊙ yr−1, and Md, 0 = 10−2 M⊙. Upper panel: torques acting on the rotation of the post-AGB star due to accretion (dashed orange line), stellar winds (dotted green line), and tides (dash-dotted red line). Lower panel: evolution of the surface rotational velocity at the equator (orange line). The blue line gives the ratio Ωspin/Ωorb, which starts at co-rotation (=1).

The upper panel of Fig. B.1 shows that early on in the evolution, accretion of gas increases the spin angular momentum, which quickly increases the angular velocity (lower panel). Because Ω_spin/Ω_orb can become larger than 18/11 in some cases, equilibrium tides can cause the eccentricity to increase rather than decrease.

In our implementation of tides in our simulations, we assume that tidal dissipation does not impact the envelope structure of the post-AGB star. Indeed, the energy dissipation rate is a small fraction of a solar luminosity, hence is several orders of magnitude smaller than the stellar luminosity.

Appendix C: Model results

In Table C.1, we present the parameters of the successful accretion models for each star from Figs. 1–3. We label each individual model in the second column. Models H1–H14 are the different models of HD 46703, while H15WW–H19WW are models where we assume a weaker wind. For RU Cen, the models are labelled R1–R15 and for EP Lyr these are labelled E1–E14. In the next three columns, we show the parameter combinations of the successful models, which are given by the starting temperature T₀, the initial accretion rate Ṁ₀ in units of M_⊙ yr⁻¹, and the initial disc mass M_d, 0 in units of M_⊙. Columns 5–7 give the evolution timescale in years and the disc properties at the end of the evolution (i.e. the predicted current values), with Ṁ_f being the final accretion rate onto the binary and M_d, f the final disc mass of the model.

In Tables C.2–C.4, we present the results of the binary evolution calculations. We have omitted models for which RLOF during the evolution could not be avoided. The first column shows the method of the calculations as presented in Table 2, while in the second column we give the labelled accretion models from Table C.1. In the remaining columns, we have from left to right the initial semi-major axis a_b, 0 in AU, the final semi-major axis a_b, f in AU, the final eccentricity e_f, the final outer disc radius R_out in AU, and the viscosity parameter α which was used in the model. The viscosity parameter was calculated from Eq. (A.2) by using the accretion properties and the initial semi-major axis.

All Tables

All Figures

Fig. 1. Depletion versus Teff plot for HD 46703 with starting temperatures of 3250 K, 3500 K, 3750 K, and 4500 K from top to bottom, respectively. The black error bar denotes the observed value of HD 46703. The legend is given at the top of the figure, in which different colours and line styles represent different initial accretion rates and initial disc masses, respectively. Ṁ0 is expressed in M⊙ yr−1 and Md, 0 is given in M⊙.

In the text

	Fig. 2. Depletion versus Teff plot for RU Cen with starting temperatures of 3400 K, 3700 K, 4000 K, and 4500 K from top to bottom, respectively. See legend at the top for definitions of colours and line styles.
In the text

	Fig. 3. Depletion versus Teff plots for EP Lyr with starting temperatures of 3500 K, 4000 K, and 4500 K from top to bottom, respectively. See legend at the top for definitions of colours and line styles.
In the text

Fig. 4. Best-fitting models in the e vs. ab plane for case 1 with both accretion and resonance torques. The black error bars show the location for each of the stars in the sample. The line styles used in this plot are the same as for the models that fit the observations in Figs. 1–3. The symbols at the end points represent models of different T0. The models start from a circular orbit with e = 0.01, which is shown by the thin black line. The grey dashed curves are lines of constant RL, in which the given temperature corresponds to the radius for this star’s luminosity such that R* = RL. Models for HD 46703 are shown in the top panel, for RU Cen in the middle panel, and for EP Lyr in the bottom panel.

In the text

Fig. 5. Time evolution of ė (top panel) and $\dot{J}$ (bottom panel) for the different mechanisms used in case 1 for HD 46703. The parameters used in this model are Md, 0 = 10−1.5 M⊙ and Ṁ0 = 10−6.0 M⊙ yr−1 with T0 = 3250 K, which corresponds to the solid orange line in the top panel of Fig. 4. The solid blue line represents the sum of the three other effects, i.e. resonances, accretion, and tides. The two most important effects on the orbital evolution are resonances and tides.

In the text

	Fig. 6. Best-fitting models of HD 46703 in the e vs. ab plane for the effect of resonances in case 2. The symbols and line styles are the same as in the upper panel of Fig. 4.
In the text

	Fig. 7. Same as Fig. 6, but for the accretion torque described in Sect. 4.3.2.
In the text

	Fig. 8. Same as Fig. 4, but for models in which the effect of tides is ignored (case 4).
In the text

	Fig. 9. Same as the bottom panel of Fig. 1 (with T0 = 3750 K), but for models with a stellar wind which is a factor of 100 weaker. In this case, more massive discs are required to explain the observed depletion pattern due to the longer evolution timescale.
In the text

	Fig. 10. Eccentricity versus semi-major axis plot for case-1 models with T0 = 3750 K. The four transparent models on the left have standard wind parameters (bottom panel of Fig. 4), while the five models on the right are the successful models of Fig. 9 with a 100 times weaker wind.
In the text

Fig. 11. Eccentricity versus semi-major axis plot for case 1 models of RU Cen. These models are similar to the middle panel of Fig. 4, except the initial eccentricity is allowed to be non-zero. The star must remain inside its Roche lobe during its evolution. Starting points are given by a small circle, the end points are shown as open symbols. For T0 = 4500 K the orbital parameters hardly change and coincide with the observed point.

In the text

Fig. B.1. Spin evolution of HD 46703 for a model with T0 = 3500 K, Ṁ0 = 10−6.3 M⊙ yr−1, and Md, 0 = 10−2 M⊙. Upper panel: torques acting on the rotation of the post-AGB star due to accretion (dashed orange line), stellar winds (dotted green line), and tides (dash-dotted red line). Lower panel: evolution of the surface rotational velocity at the equator (orange line). The blue line gives the ratio Ωspin/Ωorb, which starts at co-rotation (=1).

In the text