© The Authors 2024

1. Introduction

The post-asymptotic giant branch (post-AGB) is a short phase of contraction after the envelope of a star has been stripped, either by a dusty wind on the asymptotic giant branch (e.g. Vassiliadis & Wood 1993; Habing & Olofsson 2004; Höfner & Olofsson 2018) or, as is more likely for binary systems, through strong binary interactions (e.g. Decin 2021). A post-AGB star continually contracts and heats up at a constant luminosity of ∼1000 − 10 000 L_⊙ (Van Winckel 2003), eventually reaching the white dwarf phase in ∼10 000 yr (Bertolami 2016). In this paper, we focus on a specific class of post-AGB stars, namely post-AGB binaries surrounded by a stable, dusty circumbinary disc (for a comprehensive review, see Van Winckel 2019). These systems show strong infrared excesses in their spectral energy distributions, indicative of a large dusty disc (Kluska et al. 2022). The binarity of these disc sources’ central stars is confirmed through radial velocity (RV) measurements (e.g. Oomen et al. 2018). This revealed orbital periods in the range of a few hundred to a few thousand days and eccentricities up to ≈0.65, while also showing that the secondary companions are likely main sequence (MS) stars. Meanwhile, Galactic sources resolved through both infrared and (molecular) sub-mm interferometry (e.g. Corporaal et al. 2023; Gallardo Cava et al. 2021), as well as single-telescope imaging (Ertel et al. 2019; Andrych et al. 2023), confirm the circumbinary nature of the dusty discs.

This type of circumbinary disc is likely formed from envelope material, stripped from the evolved primary and ejected during a poorly understood phase of binary interaction during the AGB or, even earlier, on the red giant branch (RGB; Van Winckel 2019). Examples of the latter population, typically distinguished from their post-AGB analogues by their systematically lower luminosities (e.g. Kamath et al. 2016), are called post-RGB binaries. Due to the lack of reliable luminosity estimates (see the note under Table 1) and for the sake of brevity, we will hereafter refer to circumbinary disc-bearing post-AGB/RGB systems simply as post-AGB binaries, despite the possibility that some of our sample stars might be post-RGB in nature. The physics described in this paper is applicable to both types of systems.

There is strong observational evidence for ongoing re-accretion streams from the dusty post-AGB circumbinary discs to the central stars. First, the primaries’ atmospheres often appear depleted in refractory elements, which is explained through the trapping of dust in the circumbinary disc and the subsequent infall of refractory-depleted gas on the primary. This dilutes its atmosphere, leaving behind a specific chemical signature (Oomen et al. 2019, and references therein). Second, thermal emission from a smaller accretion disc around the secondary star (i.e. a circumcompanion accretion disc) has been detected in near-infrared interferometry for several systems (Hillen et al. 2016; Anugu et al. 2023). Several more such detections are as of yet unpublished (N. Anugu, priv. comm.). Current estimates of the accretion rates in the circumcompanion discs indicate that re-accretion from the circumbinary disc is the most likely mechanism feeding them (Bollen et al. 2022). To be succinct during the rest of this paper, we will use the term accretion disc to refer to a post-AGB binary system’s circumcompanion accretion disc, while circumbinary disc is used explicitly to refer to the larger, dusty disc surrounding the entire system.

Using high-resolution spectra from a more than a decade long monitoring campaign (Van Winckel et al. 2009) with the HERMES spectrograph (Raskin et al. 2011), it was shown that the accretion discs in post-AGB binaries can launch fast, collimated outflows (e.g. Gorlova et al. 2012, 2015). Such outflows, colloquially called jets, are ubiquitously observed around a large variety of accreting astrophysical objects. Jet energetics range from extremely relativistic, as when launched from active galactic nuclei, to non-relativistic in both evolved and non-evolved systems. Examples of the latter include those launched from (proto)planetary nebulae ((P)PNe) and young stellar objects (YSOs). These jets are suspected to affect the properties of both their launching object and the ambient surrounding medium. Indeed, jets in (P)PNe are being studied as crucial nebular shaping agents (e.g. Lopez et al. 1993; Clairmont et al. 2022; Lora et al. 2023), while jets in YSOs are thought to be removing a significant fraction of the angular momentum in their accretion discs, promoting stellar growth and affecting the underlying disc structure (e.g. Ray & Ferreira 2021). The jets in post-AGB binaries are also non-relativistic, and thus share similarities to these systems in terms of jet-launching physics.

Post-AGB binary jets are detected through variable Hα line profiles. Due to the high contrast between the very luminous primary and the secondary in post-AGB binaries, the latter’s light cannot be detected. Outside of superior conjunction, the Hα profile consists of the primary’s photospheric absorption profile, superposed with two emission wings, which either follow the primary’s RV or stay fixed on the centre of mass (COM). The origin of these wings are as of yet unclear. At superior conjunction, a jet launched from the accretion disc can obscure the primary, causing a deep and extended absorption feature as background photons are scattered out of the line-of-sight (LOS). Parametric modelling of such Hα profiles in a series of papers constrained the geometry and kinematics of post-AGB binary jets, indeed revealing many similarities to YSO jets and again suggesting a MS nature for the secondaries (Bollen et al. 2017, 2019, 2021, 2022). Examples of typical jet-related Hα variations can be seen in the left column of Fig. 1.

Fig. 1. Dynamic spectra (see explanation in Sect. 3.1) for the target systems’ Hα lines. Left: observed spectral time-series after spectra selection. Middle: background spectrum used as input for the RT module. Right: best fitting disc wind model (parameters summarised in Table 4). A black dashed line marks superior conjunction, while a red selection box denotes the region used for the calculation of ${\chi }_{\nu }^2$.

Still, parametric prescriptions provide little connection to the actual underlying physics of the accretion disc and the accretion-ejection process. For this reason, Verhamme et al. (2024) (from here on out Paper I) introduced the use of self-similar magnetohydrodynamic (MHD) disc wind solutions (Jacquemin-Ide et al. 2019) to describe post-AGB binary jets. These solutions are based on the extended disc wind paradigm, which is deemed one of the main contenders for the dominant jet formation mechanism in YSOs, through both theoretical and observational considerations (Ferreira et al. 2006; Launhardt et al. 2022; Moscadelli et al. 2022). They consistently describe the accretion-ejection process, allowing for a direct connection to be made between the disc wind properties and the accretion disc. Nevertheless, Paper I only considered a parameter study using one MHD solution, which was heavily tailored to YSO jets, and applied it to one target system without a minimisation routine.

In this paper, we introduce an MHD disc wind solution fitting routine to match predictions for the Hα line profile, throughout the binary orbit, to the observed HERMES time-series. We applied this fitting routine to five different post-AGB binary targets using five different self-similar disc wind solutions as input to describe the jet geometry, density and kinematics. We aim to test the validity of the MHD disc wind paradigm for jet formation, and to reveal properties of the circumcompanion accretion discs and the re-accretion streams from the circumbinary discs feeding them. The self-similar MHD solutions are introduced in Sect. 2. Our methodology is described in Sect. 3, while our target systems are introduced in Sect. 4. Results from the fitting routine are presented in Sect. 5 and extensively discussed in terms of the underlying physics in Sect. 6. Our summary and conclusions are presented in Sect. 7.

2. Self-similar disc wind solutions

To describe our disc wind models, we use the super-Alfvénic, self-similar MHD solutions of Jacquemin-Ide et al. (2019). These accretion-ejection solutions, originally intended to describe YSO jets, are henceforth termed magnetised disc winds (MDWs). They assume a vertical magnetic field threading the midplane and driving the wind. Using a self-similar ansatz, they consistently solve the equations of resistive, single-fluid 2D MHD for the combined disc-jet system (e.g. Ferreira 1997; Casse & Ferreira 2000a). This gives the relative scaling of the wind geometry, density and kinematics. Absolute values for these quantities can then be calculated for any assumed combination of accretion rate and stellar mass (as is done in Sect. 3). We describe MDWs in general, as well as the specific MDWs used in this work, below.

2.1. MDW parameters

MDWs are fully described by a complex set of dimensionless parameters. Our MDW models refer only to the volume associated with the material outflowing from the disc, but their dynamics are always computed self-consistently with the underlying disc physics. Below, we give a brief, simplified description of the MDW parameters most relevant to us, all related to the disc structure. This includes the disc aspect ratio ϵ, disc magnetization μ, disc turbulence parameter α_m, the ejection efficiency ξ and the radiative efficiency η_rad. For more rigorous details, we refer the interested reader to Paper I and Jacquemin-Ide et al. (2019). Table 2 summarises these parameters for the five MDWs considered in this work. Fig. 2 shows an example wind structure calculated using #MDW 1 and 5, which we use to showcase the general effects of the parameters on the wind. An analogous figure for all five MDWs is presented in Appendix A.

Fig. 2. Example disc wind structures using two different MDWs. Top: using #MDW 1. Bottom: using #MDW 5. Left: hydrogen number density n. Middle: poloidal velocity vp. Right: toroidal/rotational velocity vϕ. Calculated assuming the HD 52961 system and i = 70°, rin = 5 R2, rout = 0.99 RRL, 2 and Ṁin = 10−3 M⊙/yr (see Sect. 3). We note the slight discontinuity at a fixed polar angle in the n and vϕ profile. This is an artefact from the leapfrog integration performed over the Alfvén surface in the MDW solutions (Jacquemin-Ide et al. 2019), in combination with our projection onto Cartesian coordinates. This introduces no uncertainty on our results.

The scale height of the disc is parametrised by h(r) = ϵr, with r the radius along the midplane and ϵ the disc aspect ratio. Higher ϵ values imply thicker accretion discs. #MDW 1 is a very thin disc at ϵ = 0.03, while all other MDWs we consider are considerably thicker at ϵ = 0.1. Under magneto-hydrostatic equilibrium, thicker discs rotate significantly slower (see Paper I, Eqs. 2 and 16), which weakens rotation of the wind as well. This is seen in the rightmost panels of Fig. 2.

The strength of the vertical magnetic field is parametrised by the disc magnetization μ = B²/(μ₀P), with B the vertical field strength, μ₀ the permeability of the vacuum and P the total pressure (gas + radiation), all evaluated at the disc midplane. The presence of such a vertical field leads to the generation of both radial and azimuthal field components, the latter providing a torque allowing for accretion and driving ejection (Ferreira & Pelletier 1995). Increasing μ increases the dynamical importance of the magnetic field. Our MDWs range from μ ∼ 10⁻¹ − 10⁻³, covering almost the entire μ-range of super-Alfvénic solutions found by Jacquemin-Ide et al. (2019).

Because the disc is magnetised, it is also assumed to be turbulent. This MHD turbulence is described using anomalous transport coefficients, namely a turbulent viscosity ν_ν, and magnetic diffusivity ν_m. These are assumed to follow a vertical Gaussian profile, and have respective midplane amplitudes ν_ν = α_νC_Sh and ν_m = α_mV_Ah, with α_ν the standard dimensionless Shakura-Sunyaev parameter for viscous angular momentum transport (Shakura & Sunyaev 1973), α_m an analogous parameter for the poloidal magnetic diffusivity, C_S the midplane sound speed and V_A the midplane Alfvén speed (so that μ = (V_A/C_S)²). These are related via the magnetic Prandtl number, 𝒫_m = ν_ν/ν_m. The diffusivity related to the toroidal magnetic field is then given by ${\nu }_m^{\prime }={\nu }_m/{\chi }_m$, where χ_m is the diffusivity anisotropy parameter. The MDWs we consider here assume 𝒫_m = 1 and χ_m ≈ 1, consistent with simulations of turbulence induced by the magneto-rotational instability (e.g. Fromang & Stone 2009; Lesur & Longaretti 2009; Zhu & Stone 2018). Setting α_m then sets all turbulence-related terms. We note that α_ν scales as α_ν = α_m𝒫_mμ^1/2 (Salvesen et al. 2016). Increasing μ thus increases both the turbulent and the laminar (wind) torques. The MDWs considered here assume either α_m = 1 or 2.

The disc accretion rate is parametrised as Ṁ_acc(r) ∝ r^ξ, with ξ the disc ejection efficiency. This is a crucial parameter, directly linking the accretion flow to the wind ejection. Conservation of mass then gives the mass lost to one of the lobes of the wind as Ṁ_wind = 0.5 · Ṁ_in((r_out/r_in)^ξ−1), with Ṁ_in the accretion rate at the inner disc rim, and r_out and r_in the outer and inner radius of the wind-launching disc (these are fitting parameters, see Sect. 3.3). Higher ξ values imply that more material is deviated to the wind, making the wind denser and increasing the wind absorption depth in Hα (see also Paper I, Eq. 3). This effect is seen clearly in the leftmost panels of Fig. 2. Throughout the rest of this paper, we will call MDWs with a large disc ejection efficiency ξ ‘efficient’ solutions.

These MDW solutions assume isothermal magnetic surfaces, meaning that the temperature of a field line remains fixed at the temperature of its anchoring point in the midplane. This simplifying assumption has no effect on the dynamics, since most of the wind power is carried by the magnetic field. For the purely magnetically and rotationally driven MDWs considered here, also known as ’cold’ magnetocentrifugal MDWs, the highest efficiencies are achieved at lower magnetization (see Fig. 7 in Jacquemin-Ide et al. 2019). Assuming that all energy of the magnetic structure is converted to kinetic energy, the asymptotic poloidal and toroidal velocities of the wind, along a field line anchored to the midplane at radius r₀, respectively follow $v_p^{\infty }=r_0{\mathrm{\Omega }}_{{\text{K}}_0}{(2\lambda -3)}^{1/2}$ and $r{v}_{\varphi }^{\infty }=r_0^2{\mathrm{\Omega }}_{{\text{K}}_0}\lambda$, with Ω_K0 the Keplerian angular velocity at r₀ and λ the magnetic lever arm (Blandford & Payne 1982). It can be shown under certain assumptions that λ ≈ 1 + 1/(2ξ) (Ferreira 1997). Deviations from this expression are small for all MDWs considered here (see Fig. 5 in Jacquemin-Ide et al. 2019). Increasing ξ thus also leads to a decrease in the lever arm λ, reducing the asymptotic velocities that the disc wind achieves. This effect is seen in the middle and rightmost panels of Fig. 2.

The energy budget of an accretion-ejection structure, settled between outer and inner radii r_out and r_in, is written as P_acc = 2P_wind + P_rad + P_adv, with P_acc the released accretion power, P_wind the power fed into one wind lobe, P_rad the power released as radiation from the disc and P_adv the thermal power carried further inward below r_in by the material accreting onto the central object (Ferreira 1997). Dividing this expression by P_acc leads to an expression for the disc radiative efficiency η_rad = P_rad/P_acc = 1 − 2η_wind − η_adv. If we were to consistently solve for the disc thermal equilibrium, namely by taking into account turbulent heating and radiative cooling (e.g. Combet & Ferreira 2008), both η_rad and the disc aspect ratio ϵ would be an outcome of the calculations. Here, since ϵ is a prescribed parameter of the MDW solutions, η_rad is found a posteriori using the expression above, since η_wind and η_adv are readily calculated for any MDW solution. The hotter and thicker the disc (the scale height h is related to the midplane temperature), the larger η_adv (η_adv ∝ ϵ²). One might thus expect discs with a higher aspect ratio ϵ to be less bright and have a lower η_rad. While this scaling remains correct for our MDWs, the variations in η_rad among them (see Table 2) are mostly due to variations in the fraction of power carried away in the wind lobes η_wind. The disc temperature is itself a direct function of the disc accretion rate and η_rad (see Eq. (2)).

2.2. Considered MDW solutions

We consider five different MDW solutions as the basis for our disc wind models. Due to the limited understanding of the properties of post-AGB binary jets, we chose solutions covering a wide range in MDW properties. #MDW 1 is a strongly magnetised and correspondingly inefficient solution with a thin disc. It also includes a small amount of additional heating at the disc surface (Casse & Ferreira 2000a). It has been proven to work well in describing YSO jets, and has been extensively employed to model their observations (Panoglou et al. 2012; Yvart et al. 2016). Moreover, it was used in the disc wind parameter study on HD 52961 presented in Paper I. Following Paper I’s advice, we also consider various more efficient MDWs with thicker discs. #MDW 2 and 3 are weakly magnetised, more efficient solutions with a thicker disc. The geometry of #MDW 2’s wind is particularly open compared to the other considered solutions (see Fig. A.1). #MDW 4 and 5 are also weakly magnetised, thicker disc solutions, but with even higher ejection efficiencies. They are the most efficient solutions out of the families of MDWs in Jacquemin-Ide et al. (2019) assuming either α_m = 1 or 2, respectively.

3. Methods

This section summarises the methodology for modelling and fitting synthetic Hα time-series using the MDW solutions described in Sect. 2. Sect. 3.1 lists the necessary input derived from the observed HERMES spectra. Sect. 3.2 briefly describes the radiative transfer (RT) module. Sect. 3.3 then presents our fitting routine, while Sect. 3.4 discusses the current omission of error calculations for the fitted parameters. Aside from the new fitting routine, our methods are based upon those developed in Bollen et al. (2019, 2020, development of the spectral input and RT routine) and Paper I (extension of the RT routine to accommodate MDW solutions), and refer the reader to these papers for the details of the associated modelling steps, which are only concisely described here.

3.1. Spectral input

Our modelling approach requires several inputs from the full spectral time-series of each target. We briefly describe these inputs here, but refer to Appendix B for details on how they are derived.

A selection of spectra for modelling was made from the full time-series. The resulting dynamic Hα spectra for our targets are shown in the left column of Fig. 1. A dynamic spectrum plots a time-series of spectra ordered in orbital phase ϕ_orb. The y-axis represents the orbital phase, while the x-axis represents wavelength in the RV space of the considered line (from the COM frame of reference). A colour gradient interpolated between spectra then represents the normalised flux. In addition, our routine requires the binary orbital parameters, based on SB1 orbital parameters fitted to the HERMES RVs, in order to position the binary stars (Table B.1). A background spectrum was constructed, representing the observed Hα spectra outside conjunction and unaffected by jet absorption. These background spectra are shown in the middle column of Fig. 1. Finally, an additional error term was calculated to take into account the inherent, non jet-related variability of the spectra (see Fig. B.1).

3.2. Radiative transfer

After using the orbital solution to place the stars in position according to the considered ϕ_orb, the RT routine draws rays from the primary’s surface through the accretion disc wind launched from the secondary and towards the observer. Using the previously constructed Hα background spectra as the initial radiation leaving the primary, the routine performs 1D RT along the rays to model the effect of disc wind absorption. The wind density and velocity profiles are set by the selected MDW solution, projected onto a Cartesian grid. Furthermore, the disc wind was assumed to be both oriented perpendicular to the orbit and in thermal equilibrium. Spatial resolutions of ≲0.05 AU and 0.1 AU, for the directions perpendicular to and along the wind axis respectively, were achieved for all models.

3.3. Fitting routine

For the sake of robustness and in order not to miss additional degenerate minima, we employed a grid-based approach to fitting the synthetic time-series, iterating over the remaining free parameters. The values assumed in the parameter grid are summarised in Table 3. The best-fitting model for every MDW was found by minimising the reduced chi-squared, ${\chi }_{\nu }^2$, within a box in the dynamic spectrum. This box was chosen via a range in ϕ_orb and RV, denoting the region of points used for the ${\chi }_{\nu }^2$ calculation (see Fig. 1). The parameters’ physical interpretations, effect on the dynamic spectra, units and limits are briefly summarised below (see Paper I for more details).

Inclination, i (°). Sets the angle between the normal to the orbital plane (i.e. the disc wind axis) and the LOS. A lower i value implies a wider binary orbit and a higher value of the secondary’s mass, M₂. This also increases the secondary’s radius, R₂ (calculated from the MS mass-radius relationship of Demircan & Kahraman 1991), and widens the coverage of the absorption feature in ϕ_orb as the LOS passes through higher up, geometrically broader regions of the jet. It also decreases the angle between the poloidal velocity and the LOS in well-collimated winds, causing the absorption to extend to higher absolute RV values. It lies between 0° and 90°.

Inner launching radius, r_in (R₂). The inner launching radius of the disc wind at the disc midplane. The field line anchored to this point is used as the inner boundary for the considered wind region. Increasing r_in cuts out the high velocity components of the wind close to the secondary. It is expressed in units of the secondary star’s radius, R₂, and is limited by the secondary’s stellar surface, r_in ≥ 1 R₂.

Outer launching radius, r_out (R_RL, 2). The outer launching radius of the disc wind. Its field line is used as the outer boundary of the wind. Increasing r_out adds material with low velocity to the outer wind regions, deepening the low RV absorption. This also makes the outer regions broader, making the absorption feature in the synthetic spectra wider in ϕ_orb coverage. It is expressed in units of the secondary’s Roche radius, R_RL, 2. It is limited by the inner launching radius and the Roche radius R_RL, 2 > r_out > r_in. r_out ≥ R_RL, 2 would namely imply that the circumcompanion accretion disc spills over onto the primary.

Inner rim accretion rate, log(Ṁ_in) (M_⊙/yr). The rate of material going through the inner edge of the accretion disc at r_in, and falling further onto the star. Disc wind densities increase proportionally to Ṁ_in (e.g. Jacquemin-Ide et al. 2019), and so does the total rate of mass lost to the wind lobes: Ṁ_wind = 0.5 · Ṁ_in((r_out/r_in)^ξ−1). This does not affect the wind geometry or kinematics, but causes deeper absorption signatures. The accretion rate at the outer disc edge then follows Ṁ_out = (r_out.r_in)^ξ Ṁ_in, expressing the feeding rate in the re-accretion stream from the circumbinary disc onto the accretion disc.

Disc wind temperature, T_wind (K). Homogeneous equilibrium temperature for the disc wind assumed during RT. Higher temperatures correspond to higher electron population densities in the first excited hydrogen state, implying stronger Hα absorption. This effect is limited by the wind’s thermal blackbody emission, which will eventually dominate over the increase in absorption as T_wind increases. Since the wind feature is seen in absorption, we are limited by T_wind < T_pAGB, with T_pAGB the temperature of the primary (see Table 1), as the feature would appear in emission otherwise.

Primary radius, R₁ (R_⊙). The radius of the primary post-AGB star. Increasing it increases the ϕ_orb coverage of the wind absorption feature, since the primary is obscured by the wind for a longer duration. In order to limit computational time, we fixed it to the values found from parametric modelling by Bollen et al. (2020, 2021, 2022) (R_1, par, see Table 1). R₁ is one of the more robust, reliably retrieved parameters found from the parametric modelling.

3.4. Omission of parameter errors

Remarkably, parametric modelling of the Hα time-series using Markov chain Monte Carlo (MCMC) resulted in very small formal errors on the fitting parameters (e.g. Bollen et al. 2022). The inclination, for example, had formal errors typically only of the order ∼0.1°. Given the expected correlations, this already seems strongly underestimated. Moreover, runs with relatively minor changes in the input spectra or parameter exploration range resulted in best fit inclinations that could easily differ by > 5°. In the case of BD+46°442, we found an equally well-fitting, degenerate minimum that differed by as much as 20°.

The current MCMC implementation (the emcee package by Foreman-Mackey et al. 2013) seems both overconfident and badly conditioned, and does not properly take degeneracies into account. Preliminary error determination for the disc wind model fits in this work, using a simple χ² criterion, was similarly overconfident. As a result, we refrained from giving formal errors on our fitted parameters. We suggest to use the maximum step size in the parameter grid as a lower limit for the parameter errors (Table 3). In following works, parameter degeneracies can be significantly lifted by assuming inclinations from interferometric modelling of the circumbinary disc (assuming co-planarity with the accretion disc). Though outside the scope of this paper, the use of MCMC algorithms more suitable to high-dimensional and multi-modal distributions, such as Hamiltonian Monte Carlo (e.g. Hoffman & Gelman 2011), could lead to more robust error determination.

4. Target stars

This section describes our target binaries and their observed Hα spectra. Table 1 summarises the properties of these systems used for their selection and modelling. The observed dynamic Hα spectra, after spectra selection from the full time-series, are shown in the left column of Fig. 1. Most Galactic post-AGB binaries (see Kluska et al. 2022 for a full catalogue) for which HERMES time-series were assembled indeed show jet absorption. For this paper, we selected five targets covering a good range in orbital properties, metallicity, refractory depletion strength, inclination from previous parametric modelling and visual morphology of the Hα spectra. As a result, these five systems form a small but representative sample of the diversity seen in jet-launching post-AGB binaries.

4.1. HD 52961

HD 52961 shows extreme refractory depletion (Waelkens et al. 1991), as reflected by its remarkably low metallicity of [Fe/H] = −4.8. It has a fairly long orbital period at P ≈ 1300 d (known values in the population range from ≈100 − 2500 d) and a substantially eccentric orbit at e = 0.23 (values range from ≈0.0 − 0.65). Previous studies detected infrared emission from UV-excited polycyclic aromatic hydrocarbons (PAHs) and C⁶⁰ fullerenes (Gielen et al. 2011). Most importantly, it was the target of the parameter study performed using #MDW 1 in Paper I. Its observed dynamic Hα spectrum presents a deep double-peaked jet absorption feature slightly offset from conjunction, covering ϕ_orb ≈ 0.10 − 0.45. The feature reaches down to RV ≈ −70 km/s. In addition, it shows a small amount of redshifted jet absorption up to RV ≈ 60 km/s at ϕ_orb ≈ 0.3, which suppresses the rightmost emission wing of the background spectrum. The background emission wings follow the RV curve of the primary.

HD52961’s position in the HR diagram is close to the region occupied by the RV Tauri stars, a group of radially pulsating post-AGBs in the population II Cepheid instability strip. RV Tauri lightcurves show successive shallow and deep minima due to the pulsations (Preston et al. 1963; Alcock et al. 1998). If this is the only source of photometric variability in a system, it is called an RVa variable. Post-AGB binaries can, however, also show a long-term brightness modulation that follows the orbital period. This modulation is caused by high system inclinations, causing periodic obscuration of the primary by the circumbinary disc as the LOS grazes its edge (Manick et al. 2017). Regardless of whether they show the characteristic brightness dips of RV Tauri stellar pulsations or not, these systems are called RVb variables. HD 52961 is an example of an RVb. It has a long term periodic brightness modulation but, due to it lying just outside the instability strip, its stellar pulsations do not follow the RV Tauri pattern of successive shallow and deep brightness dips (Van Winckel et al. 1999). Previous parametric jet modelling implied a high inclination of i_par ≈ 70°, which is consistent with its RVb classification.

4.2. BD+46°442

BD+46°442 is a non-depleted system with moderate metallicity (typical values range from ≈ − 0.3 − 0.0) and an almost circularised, very short period orbit at P ≈ 140 d and e ≈ 0.1. Previous parametric modelling implied an intermediate inclination of i_par ≈ 50°. It does not show strong pulsations. Its dynamic Hα spectrum shows a jet absorption feature down to RV ≈ −250 km/s, slightly offset from conjunction. The feature spans ϕ_orb ≈ 0.4 − 0.7. The background emission wings stay centred on the COM. It was the first of the HERMES-observed systems to which the hypothesis of a circumcompanion-launched jet was applied (Gorlova et al. 2012).

4.3. TW Cam

TW Cam is a non-depleted system with moderate metallicity. It has an intermediate orbital period at P ≈ 660 d and substantial eccentricity at e = 0.25. It is an RV Tauri pulsator, with strong pulsations of up to 20 km/s amplitude in radial velocity (see Fig. 9 in Manick et al. 2017). This results in strong variability of the observed Hα spectrum. We refrained from trying to clean the observed spectra (e.g. by removing shock-affected spectra). It was tentatively classified as an RVb by Manick et al. (2017) due to a photometric variability component of low amplitude (ΔV = 0.07 mag) following the orbital period, implying a high inclination under the circumbinary disc-grazing LOS hypothesis. Giridhar et al. (2000), however, identified it as an RVa variable (not disc-grazing), and parametric jet modelling also predicted a low inclination of i_par ≈ 26°. This is due to the fact there is some jet absorption present throughout the entire orbit, which occurs when the inclination is smaller than the jet opening angle. This, together with the low amplitude of the long-term photometric variability, suggests it might not be due to periodic obscuration by the circumbinary disc. TW Cam’s dynamic Hα spectrum shows wide, double-peaked jet absorption centred at conjunction, reaching down to RV ≈ −300 km/s. The feature covers almost the entire orbit, but becomes more prominent at ϕ_orb ≈ 0.35 − 1.0. The background emission wings stay centred on the COM.

4.4. IRAS 19135+3937

IRAS 19135+3937 has only a weak signature of refractory depletion, and has a slightly sub-solar metallicity. It has the third shortest orbital period out of all known post-AGB binary orbit solutions, with P ≈ 130 d (Oomen et al. 2018), and a moderate eccentricity at e = 0.13. It has been classified as a class d semi-regular variable. This variability has previously also been attributed to obscuration by the circumbinary disc (Gorlova et al. 2015), implying a high inclination. Correspondingly, previous parametric modelling fit a high inclination at i_par ≈ 72° . The primary’s photosphere seems relatively stable (Gorlova et al. 2015). It has an Hα jet absorption feature slightly offset from conjunction, ranging down to RV ≈ −360 km/s and spanning ϕ_orb ≈ 0.3 − 0.7. In addition, it shows a small amount of redshifted absorption around ϕ_orb ≈ 0.40, reaching up to RV ≈ 50 km/s. The background emission wings stay centred on the COM.

4.5. HP Lyr

HP Lyr shows a very strong depletion pattern, though not quite as strong as that of HD 52961 (Oomen et al. 2019). It has the second longest period out of all known post-AGB binary solutions at P ≈ 1800 d (Oomen et al. 2018), and a significant eccentricity at e = 0.20. It is classified as an RVa variable due to the lack of long-term brightness variations (Manick et al. 2017). This is consistent with the fact that parametric jet modelling resulted in an almost pole-on inclination of i_par ≈ 19°. The RV pulsations of the primary are on the order of 5 km/s (see Fig. 11 in Manick et al. 2017). The observed background Hα spectrum is thus variable, though less severely than in the case of TW Cam. It has a deep Hα jet absorption feature, reaching down to RV ≈ −120 km/s and spanning ϕ_orb ≈ 0.15 − 0.7. Approximately from ϕ_orb ≈ 0.42 − 0.70, the jet absorption feature detaches from the photospheric absorption line, letting the leftmost background emission wing shine through. The background emission wings stay centred on the COM.

5. Results

In this section, we describe the resulting Hα profiles of our fitted models and compare them to the observations. The fitting routine of Sect. 3.3 was run using each of the MDWs in Table 2, resulting in five disc wind model fits per target, one per MDW solution. The parameters and dynamic Hα spectra for the best fitting models are shown in Table 4 and the right column Fig. 1, respectively. For completeness, the results of the fitting procedure using all MDW solutions are summarised in Appendix C.

5.1. HD 52961

#MDW 4 provides the best fit for HD 52961’s Hα time-series at ${\chi }_{\nu }^2$ = 5.32. It reproduces both the general absorption depth, RV extent and the phase positioning of the high ϕ_orb absorption peak. The lower ϕ_orb peak in the model, in contradiction with the observations, is redshifted. As discussed by Paper I, this redshift is induced by the prograde rotation of the wind material, causing the material’s velocity to be pointed away from the observer as the LOS to the primary passes the wind’s leading edge. Correspondingly, this also blueshifts the lower ϕ_orb absorption peak, as the trailing edge passes the LOS. This results in a strong RV asymmetry of the modelled wind absorption. Instead of being placed at its observed position of ϕ_orb ≈ 0.20, the modelled low ϕ_orb absorption peak is located at ϕ_orb ≈ 0.30, where it tries to fit the small redshifted absorption feature seen in the observations.

Compared to the other targets, the rotation effect is remarkably strong in HD 52961’s best fit. This is due to its high fitted inclination of i = 72°, causing the LOS to pass through regions close to the wind’s base, where the rotation is strongest. Such a high fitted inclination is consistent with its RVb classification. None of the other models manage to suppress the wind rotation enough to blueshift the low ϕ_orb peak. #MDW 1, 2 and 3 perform significantly worse at ${\chi }_{\nu }^2$ ≤ 6.92 due to the strong effect of rotation (see Fig. C.1 and Table C.1). Only #MDW 5, the most efficient of the thicker disc solutions, can reduce the wind rotation more and reaches ${\chi }_{\nu }^2$ = 5.63, but only manages to shift the low ϕ_orb peak into the RV position of the photospheric absorption line. We note that #MDW 5 fits at i = 56°, which is significantly lower than for the other MDWs (i > 70°). Given #MDW 5’s high efficiency, the density of the wind is too high close to the base. In order to compensate, the fit tends towards a low inclination in order to move the LOS to less dense regions. However, i = 56° seems fairly low, given the system’s robust RVb classification.

Despite significant mismatches the ${\chi }_{\nu }^2$ values are still fairly low. This is a general feature of our fits due to the large selection boxes used, which include a lot of continuum points. These hold little statistical significance, but by construction of the background spectrum match up well to the observations. As a result, different MDW fits can thus sometimes fit the global Hα pattern very differently while still achieving similar ${\chi }_{\nu }^2$ values (e.g. IRAS 19135+3937 #MDW 3 and 4, see Fig. C.1). We thus stress that the choice of the best fitting MDW model was guided by the ${\chi }_{\nu }^2$ values, but ultimately depended on how well the specific MDW managed to recover the global Hα pattern.

5.2. BD+46°442

BD+46°442 is best fitted by #MDW 3 at ${\chi }_{\nu }^2$ = 6.44. This model reproduces the absorption depth, ϕ_orb coverage and negative RV extent of the observed absorption feature. Two glaring mismatches with the observations remain, however. The first is the rotation-induced redshifted absorption feature centred at ϕ_orb ≈ 0.41, which is not seen in the observations. This feature is also present in the #MDW 1 and 2 fits, which both fit more poorly at ${\chi }_{\nu }^2$ ≤ 6.73. Rotation also blueshifts the lower ϕ_orb peak. This detaches the wind absorption from the photospheric absorption at ϕ_orb ≈ 0.55 − 0.70. The high efficiency #MDW 4 and 5 models do manage to suppress wind rotation to a significant degree, but produce poorer fits at ${\chi }_{\nu }^2$ ≤ 6.85. Their modelled absorption features are too broad in ϕ_orb, and extend too much downwards into negative RV values. This is due to their lower fitted inclinations (i < 45° versus i > 55° for the other fits), which are assumed because the wind regions close to the base are too dense. In addition, close to superior conjunction, the leftmost background emission wing shines through all fits, implying a general lack of low RV absorption in all the models.

5.3. TW Cam

TW Cam is best fitted by #MDW 3 at ${\chi }_{\nu }^2$ = 5.42. This model’s wind absorption feature reproduces the extent in ϕ_orb and negative RV range of the main absorption peaks in the observations fairly well. The absorption depth is generally too shallow, most likely because the parameter grid is limited to Ṁ_in ≤ 10⁻³ M_⊙/yr. The fit assumes a low inclination of i = 24°. This is lower than the typical wind opening angle at the point where the outer wind edge and LOS rays cross, causing the primary to be obscured and the model to show wind absorption at all times. The low inclination implies we look through regions high above the wind’s base, so the effect of rotation is a lot less pronounced. The observed blueshifted absorption peaks are fairly asymmetrical, where the lower ϕ_orb peak at ϕ_orb ≈ 0.5 is slightly more narrow in ϕ_orb coverage and is concentrated at more negative RV values than the other peak at ϕ_orb ≈ 0.8. In comparison, the modelled absorption peaks are fairly symmetrical. This is most likely due to a neglect of a wind tilt in our current models (see e.g. Bollen et al. 2020), which can be dealt with by adding a jet tilt parameter. Since this is a second order effect, it is postponed for future work. The wind absorption away from superior conjunction (e.g. at ϕ_orb ≈ 0.0 − 0.3) detaches from the photospheric absorption due to high poloidal velocities, causing a mismatch with the observations. There is too much absorption in a broad RV band, and the background spectrum immediately blueward of the photospheric absorption shines through the model. The #MDW 1 fit shows a similar effect, and correspondingly fits almost equally well at ${\chi }_{\nu }^2$ = 5.59. The higher efficiency #MDW 4 and 5 indeed manage to reattach the wind absorption outside of superior conjunction to the photospheric absorption due to their inherently lower poloidal velocities. In addition, #MDW 5 reproduces the observed absorption depth much better. They however fail in reproducing the observed wind absorption’s width in ϕ_orb, and as a result fit significantly poorer at ${\chi }_{\nu }^2$ ≤ 6.44. The #MDW 2 fit is also poorer at ${\chi }_{\nu }^2$ = 6.57. We thus have been unable to find the right combination of kinematics and density profile to fully reproduce TW Cam’s observations.

We note that the fitted inclinations are consistently on the lower end at i < 44°. Parametric modelling also implied a low inclination at i_par = 25.7° (Table 1). Assuming a wind axis approximately perpendicular to the orbital plane, this strongly challenges the RVb classification of TW Cam by Manick et al. (2017), instead supporting the RVa classification by Giridhar et al. (2000).

5.4. IRAS 19135+3937

IRAS 19135+3937 is best fitted by #MDW 4 at ${\chi }_{\nu }^2$ = 3.30. This model matches the observations very well, reproducing both the general RV extent, ϕ_orb width, absorption depth and weak redshifted absorption indent at ϕ_orb ≈ 0.40. Similar to the BD+46°442 models, a lack of low RV absorption close to superior conjunction is indeed present. This is also the case for all other fits for this target. In the #MDW 4 and 5 models this effect is however almost completely reduced due to the inherently lower poloidal velocities. #MDW 1 and 2 are poorer fits at ${\chi }_{\nu }^2$ ≤ 3.96, but #MDW 3 and 5 fit almost as well as #MDW 4 at ${\chi }_{\nu }^2$ = 3.49 and 3.53, respectively. Nevertheless, the #MDW 3 model’s wind absorption is too shallow, and suffers the same mismatches of overestimated rotation and lack of low RV absorption as the BD+46°442 #MDW 3 model. Aside from the #MDW 2 model, fitted inclinations are relatively low at i < 48°. #MDW 5 fits an especially low inclination at i = 28°, while parametric modelling fitted a high inclination at i_par = 71.8°. The parametric prediction is more consistent with the identification of obscuration by the circumbinary disc as the cause for the observed photometric variability (Gorlova et al. 2015).

5.5. HP Lyr

HP Lyr is best fitted by #MDW 4 at ${\chi }_{\nu }^2$ = 4.91. The model reproduces the observation very well, matching the RV extent and ϕ_orb width of the wind feature. The absorption depth is only slightly too shallow. The wind absorption from ϕ_orb ≈ 0.42 − 0.70 does not detach as much from the photospheric absorption line as in the observations due to excessive suppression of rotation, as the #MDW 1 and 3 models do reproduce this feature. #MDW 1, 2 and 3 however fit much worse overall with ${\chi }_{\nu }^2$ ≤ 6.29. The #MDW 5 model is extremely similar to the #MDW 4 one, reaching ${\chi }_{\nu }^2$ = 5.05, though with a more shallow absorption depth due to its lower inclination of i = 17° (i ≥ 29° for all other MDWs). As was the case for HD 52961, #MDW 5 most likely assumes this lower inclination to avoid the densest regions close to the wind base. However, this implies a very large mass for the secondary star at M₂ = 4.3 M_⊙. In contrast, known SB1 orbital parameters imply M₂ follows a distribution with a mean of 1.09 M_⊙ and standard deviation of 0.62 M_⊙ (assuming random orbit orientations; Oomen et al. 2018). An upper limit on M₂ might serve as a good lower limit for i during future fitting runs.

6. Discussion

This paper expands on Paper I in the pursuit to test MDW-based models in the context of post-AGB binary jets. We have successfully introduced our fitting routine and have made quantitative comparisons between the Hα line predictions for five different MDWs and the observed HERMES time-series of our targets. In this section, we will proceed with an in-depth discussion on several aspect of our fitted models and the associated physics, using the results of Paper I as a comparative guideline.

6.1. Jet rotation

The rotation of jets is a natural consequence of their formation, and is confirmed by direct observations of rotation and/or winding of the magnetic field lines in YSO jets (e.g. Launhardt et al. 2022; Moscadelli et al. 2022). From the lack of rotation-induced redshift in the observed HD 52961 wind absorption feature, Paper I concluded that wind rotation in this system must be strongly suppressed compared to what is predicted by the YSO-tailored #MDW 1. As suggested by the authors, we considered more weakly magnetised, thicker and more efficient MDW solutions, which naturally lead to weaker jet rotation.

For TW Cam and HP Lyr, both systems fitted under a low inclination (i < 28° and 37°, respectively, excluding the poorly fitting #MDW 2), the rotation-induced redshift is relatively weak for all considered MDW solutions (see Fig. C.1). As shown in Paper I, it are mostly the strongly rotating, sub-Alfvénic regions close to the disc surface which are responsible for the redshifted absorption. Since the low inclinations cause the LOS to pass through higher up wind regions, the rotation is mostly avoided.

The good fits for IRAS 19135+3937, using #MDW 4 and 5, both reproduce a small but noticeable rotation-induced redshift in the absorption at ϕ_orb ≈ 0.40, matching a similar feature in the observations. The other, less efficient solutions predict redshifts that are too strong (see Fig. C.1). This implies that while rotation in the actual IRAS 19135+3937 jet is indeed relatively weak compared to typical YSO jets, it is still present to some degree.

Redshifted wind absorption is still strongly present in the best fit model for HD 52961 using #MDW 4. The higher ξ #MDW 4 and 5 clearly suppress the rotation-induced redshift more compared to the lower ξ #MDW 2 and 3, and the lower ξ and ϵ #MDW 1. Due to the high fitted inclinations, rotation can even dominate the velocity field when using #MDW 1, 2 and 3 (see Fig. C.1). Nevertheless, the rotation suppression is never sufficient to match the data.

Also in the case of BD+46°442, the best fit using #MDW 3 does not seem efficient enough to suppress the redshifted feature, which is completely absent in the observations (see Fig. 1). While the #MDW 4 and 5 models do manage to fully suppress the redshift, their absorption features are slightly too broad in ϕ_orb, and reach RV values that are too negative (see Fig. C.1). These models fit lower inclinations of i = 44° and 28°, respectively, compared to i > 55° for the other models. Since #MDW 3 already fitted the absorption depth seen in the observations well at i = 56°, these more efficient solutions are too dense close to the base of the wind, and the inclination is correspondingly decreased so the LOS avoids these dense regions. As a result, the absorption feature becomes broader in phase, and since the angle between the LOS and the poloidal velocities decreases, absorption extends to more negative RV.

For HD 52961 the rotation suppression is simply not strong enough using the currently considered MDW solutions, while for BD+46°442 they do not provide the correct combination of kinematics and density profile. New MDW solutions with thicker accretion discs at ϵ > 0.1 must be considered to solve these discrepancies. As pointed out by Paper I, however, purely magnetically driven cold disc winds start failing as the disc becomes too thick. Discs that rotate too slowly simply cannot provide the rotational push to launch the material (Casse & Ferreira 2000a). As we will discuss in Sect. 6.6, the inclusion of irradiation by the primary star might overcome this limitation.

6.2. Accretion rates

Paper I raised the issue of lifetimes for the circumbinary discs, which feed the jet-launching accretion discs. Their #MDW 1 model preferred a feeding rate of Ṁ_out ≈ 10⁻⁴ M_⊙/yr. Unlike in some YSOs, there is no natal envelope which can replenish the circumbinary disc. Given that they have a typical total mass of up to M_cb = 10⁻² M_⊙ (derived from radio observations, e.g. Bujarrabal et al. 2013, 2015, 2016, 2017, 2018; Gallardo Cava et al. 2021), the circumbinary discs should then dissipate on a timescale τ_cb = M_cb/Ṁ_out ∼ 100 yr. The post-AGB phase lifetime for a lower-mass star can however easily be τ_pAGB ≳ 10 000 yr, and the re-accretion onto the primary might extend this by a factor of a few (Bertolami 2016; Oomen et al. 2019). With a large fraction of ≈30% of all observed post-AGBs being circumbinary disc-bearing (Kamath et al. 2014, 2015; Kluska et al. 2022), they should not be able to dissipate so quickly, and we’d expect τ_cb ∼ τ_pAGB. To push down the required feeding rate and increase the τ_cb estimate, Paper I suggested the use of more efficient MDW solutions, as such solutions can divert a larger fraction of the accretion stream toward the disc wind, requiring less material to be fed in order to obtain the same Hα absorption depth.

Across our targets, #MDW 1 (ξ = 0.04) fits in the range Ṁ_out ∼ 10⁻⁴ − 10⁻³ M_⊙/yr (see Table C.1). #MDW 2 and 3 (ξ ≈ 0.10) manage to push the lower end down by an order of magnitude to Ṁ_out ∼ 10⁻⁵ − 10⁻³ M_⊙/yr. #MDW 4 (ξ = 0.215) fits Ṁ_out ∼ 10⁻⁵ − 10⁻⁴ M_⊙/yr, while #MDW 5 (ξ = 0.305) fits Ṁ_out ∼ 10⁻⁶ − 10⁻⁴ M_⊙/yr. More efficient MDWs thus generally do predict lower feeding rates, but even for #MDW 5, the most efficient of all known cold magnetocentrifugal solutions (Jacquemin-Ide et al. 2019), the high end of the range still predicts τ_cb ∼ 100 yr, which seems much too short. Evidently, the feeding rates must be pushed down even further, ideally by one to two orders of magnitude. We have, however, reached the limit of efficiency for purely magnetocentrifugal cold winds from an ϵ = 0.1 disc with #MDW 5 (Jacquemin-Ide et al. 2019). As was put forward as a possibility in the conclusions of Paper I, an additional driving force is needed. Irradiation from the bright post-AGB primaries, with luminosities easily over 1000 L_⊙ (Table 1), is bound to deposit heat onto the accretion discs’ surfaces, adding an extra thermal push to the wind material. Indeed, such ‘warm’ magnetothermal solutions have already been explored in the context of YSO jets with coronal heating (Casse & Ferreira 2000b), which can readily reach huge efficiencies of ξ > 0.5. Such solutions will need to be explored in future work to push down the fitted values for Ṁ_in and Ṁ_out, and solve the current circumbinary disc lifetime discrepancy. Lower Ṁ_in values will also lead to slightly colder accretion disc predictions (see Sect. 6.5) and lower secondary field strengths required to truncate the accretion discs through magnetospheric funnelling (see Sect. 6.3). It will also allow for thicker, more slowly rotating disc solutions at ϵ > 0.1 to be explored (see Sect. 6.1).

6.3. Disc truncation and magnetospheric funnelling

Considering the inefficient #MDW 1 solution, Paper I discussed the necessity of large inner launching radii ≳20 R₂ in order to match the HD 52961 observations. This would cut out the inner wind regions, which were too fast in their model. As there is no apparent physical reason for this inner launching radius to be different from the inner radius of the bulk accretion disc itself, this seemed to imply that the HD 52961 accretion disc had a large inner cavity and was heavily truncated. The required inner radii were, however, deemed too high for this truncation to be possible via magnetospheric funnelling, a mechanism often invoked to explain disc truncation around YSOs (e.g. Bouvier et al. 1999, 2003), where the stellar magnetic dipole disrupts the accretion flow and funnels it along the magnetospheric field lines (see e.g. Bessolaz et al. 2008). To relieve this tension, Paper I suggested the use of more efficient MDW solutions with inherently lower poloidal velocities, which we consider here.

Both HD 52961 and IRAS 19135+3937 have a best fit with r_in = 1 R₂, implying an inner disc rim that reaches the stellar surface. This is consistent with the fact that they are best fitted by #MDW 4, the second-most efficient MDW solution at ξ = 0.215. Contrary to the solutions with lower ξ, the inner jet regions do not need to be cut out to suppress absorption at high negative RVs (see Table C.1). BD+46°442 and TW Cam are best fitted by #MDW 3, which is more efficient than #MDW 1, but still less than #MDW 4, and correspondingly fit an inner rim slightly further away at r_in = 5 R₂. It is not clear if the discs are truly truncated in these two systems, since a more careful selection of the used MDW solution might result in best fits at r_in = 1 R₂. In any case, these fits are consistent with the discs either reaching the stellar surface, or at least being truncated out to only several stellar radii.

Following Paper I and considering magnetospheric funnelling as the truncation mechanism, the strength of the stellar dipole at the stellar equator required to truncate the disc can be written as (Bessolaz et al. 2008):

$$\begin{array}{c}\hfill B_2\sim 11\mathrm{kG}·{\left(\frac{r_{\mathrm{in}}}{R_2}\right)}^{\frac{7}{4}}{\left(\frac{{\dot{M}}_{\mathrm{in}}}{{10}^{-4}{\mathrm{M}}_{\odot }/\mathrm{yr}}\right)}^{\frac{1}{2}}{\left(\frac{M_2}{1{\mathrm{M}}_{\odot }}\right)}^{\frac{1}{4}}{\left(\frac{R_2}{1{\mathrm{R}}_{\odot }}\right)}^{-\frac{5}{4}}.\end{array}$$(1)

To start truncating the disc (r_in > 1 R₂), the reasonably well-fitting models for the four targets discussed above imply magnetic field strengths greater than several kG and up to ∼15 kG. For comparison, a typical heavily accreting YSO at say Ṁ_in ≈ 10⁻⁷ M_⊙/yr can achieve mean field strengths of up to a few kG (e.g. Johns-Krull et al. 1999), and at such accretion rates this is strong enough to truncate the accretion disc out to several stellar radii (Bessolaz et al. 2008). Our secondaries are mature main sequence stars, not entirely dissimilar to for example the Sun. The average equatorial solar dipole strength is only on the level of ∼1 G, implying that our secondaries’ inherent field strengths would be much too low for truncation at the current required accretion rates. However, if the secondary’s magnetic field is driven by a stellar dynamo (see e.g. Vidotto et al. 2014; Wurster et al. 2018), the strong accretion streams in our systems could be expected to increase stellar rotation and rejuvenate the stellar magnetic field. If, in addition, the required Ṁ_in rates are indeed much lower than our current models imply (see Sects. 6.2 and 6.6), the accretion discs might be magnetically truncated, though in any case still close to the stellar surface. In future modelling, it will thus be reasonable to assume that the accretion disc either extends down to the stellar surface or is truncated out to several stellar radii at most, meaning r_in ∼ 1 R₂. Constraints on the secondaries’ field strengths might provide more stringent limits on r_in, but the brightness contrast with the primaries unfortunately precludes direct determination of B₂ from spectra, using for instance Zeeman broadening.

While the use of more efficient MDW solutions does lead to much smaller fitted r_in values for four out of five of our targets (Table C.1), one exception sticks out. The wind absorption feature in HP Lyr only extends down to RV ≈ −120 km/s, despite its low inclination (compare to TW Cam, where it reaches RV ≈ −320 km/s). Correspondingly, HP Lyr is fitted with r_in ≥ 20 R₂ for all MDW solutions (see Table C.1). HP Lyr’s best fit would imply a field strength of B₂ ≈ 1.6 MG. The #MDW 5 fit (Table C.1) predicts the lowest field strength value for this target, at B₂ ≈ 800 kG. These values are extraordinarily high, and impossible to reach via a MS stellar dynamo (e.g. because of saturation; Vilhu 1984). Even if the accretion rate in our HP Lyr models is reduced to ∼10⁻⁸ M_⊙/yr, the required field strength would only be brought down to ∼16 kG. Magnetospheric funnelling simply cannot be responsible for truncating HP Lyr’s accretion disc out to this large a radius, and the accretion disc should extend further down towards the secondary.

The consistently large required r_in value for HP Lyr’s models could be an artefact from the strongly simplified thermal structure in the RT routine, where we assumed a homogeneous wind temperature. This assumption forces all material included in the model to be set to a fairly high temperature, causing the entire wind structure to absorb. If HP Lyr’s accretion disc wind has an inverse temperature gradient in reality, with colder inner regions and hotter outer regions (though still bellow T_pAGB), only the hydrogen in the slow outer regions would be hot enough to absorb significantly in Hα. In that case, our fitting routine accounts for this by simply cutting out the fast inner wind regions. An inverse temperature gradient might, for example, be caused by interaction between the outer disc wind regions and the fast stellar wind launched by the post-AGB primary (with mass loss rates up to 10⁻⁷ M_⊙/yr; Van Winckel 2003), if the latter is strong enough and our accretion rate estimates are reduced. Future modelling of HP Lyr, and possibly similar objects that are fitting at high r_in values despite the use of very efficient MDW solutions, warrants a more careful calculation and analysis of the thermal wind structure. This can be done either in detail by solving the energy equation along the field lines for specific disc wind models, or, at lesser computational expense, by including parametrised wind temperature profiles in the RT routine.

6.4. Missing low RV absorption

Paper I noticed a consistent lack of low RV absorption in their #MDW 1 model compared to the HD 52961 observations, occurring close to superior conjunction. As a result, the Hα emission component in the background spectrum shined through the models. The authors attributed this to the void cavity around the disc wind axis in the models, which is present by construction. Such cavities are indeed confirmed observationally in disc winds around YSOs (e.g. López-Vázquez et al. 2023). As seen in Fig. 1, the best fitting model for BD°46+442 shows this effect noticeably. As does the IRAS 19135+3937 best fit, to a lesser extent. This feature persists in the other fits for these systems as well (Fig. C.1). For BD°46+442 and IRAS 19135+3937, the current disc wind models alone are unable to provide enough low RV material to explain the observations, even with higher efficiency MDW solutions.

There are, however, physical reasons to expect extra outflow components inside the disc wind cavity (Ferreira et al. 2006), which might be able to provide the missing absorption. If magnetospheric funnelling onto the secondary star is indeed active close to the stellar surface, we would expect a thermally driven stellar jet to be launched from the secondary’s poles and be collimated by the star’s open field lines (Matt & Pudritz 2005). Even if the disc does extend down to the stellar surface, the collision of the accreting material with the equator might heat up the stellar atmosphere. Hence, a thermal wind from the equator, externally confined by the disc wind’s magnetic field, might be launched as well.

In Appendix D we derive Eq. (D.6), a simple estimate of the required extra mass-loss in the cavity, Ṁ_cav, needed to supply the missing absorption (assuming a homogeneous cylindrical flow at T = T_wind). If Ṁ_cav ≪ Ṁ_in, we know that the extra outflow can be fed by material passing the inner disc rim in our models. If Ṁ_cav > Ṁ_in, our current models are inconsistent with the addition of an extra outflow component. Applying Eq. (D.6) to BD+46°442’s best fit, we derive Ṁ_cav/Ṁ_in ∼ 0.1%. Analyses of other model fits that show an analogous lack of low RV absorption (see Fig. C.1) result in similar values. While Eq. (D.6) is derived under idealised conditions in terms of opacity and kinematics, this still shows that the inner rim accretion rate in our current models, in combination with an additional low-efficiency jet launching mechanism in the cavity, can provide the missing absorption. The implementation of such an additional outflow component in our models, however, currently falls outside the scope of this paper.

6.5. Accretion disc temperature

Following Paper I and assuming that the accretion disc is geometrically thin and radiates like an optically thick blackbody, its effective temperature profile follows (generalised version of the expression in e.g. Frank et al. 2002):

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill T_{\mathrm{eff}}(r)& ={({\eta }_{\mathrm{rad}}·\frac{G{M}_2{\dot{M}}_{\mathrm{acc}}(r)}{8\pi \sigma r^3})}^{\frac{1}{4}}\hfill \\ \hfill & \approx 36.4\mathrm{kK}·{\left(\frac{{\eta }_{\mathrm{rad}}{\dot{M}}_{\mathrm{acc}}(r)}{{10}^{-4}{\mathrm{M}}_{\odot }/\mathrm{yr}}\right)}^{\frac{1}{4}}{\left(\frac{M_2}{1{\mathrm{M}}_{\odot }}\right)}^{\frac{1}{4}}{\left(\frac{r}{1{\mathrm{R}}_{\odot }}\right)}^{-\frac{3}{4}},\hfill \end{array}\end{array}$$(2)

with σ the Stefan-Boltzmann constant. We note how sensitive this is to the radii to which the disc reaches.

We use two example model fits for our discussion: HP Lyr #MDW 3 and TW Cam #MDW 2 (see Table C.1). Their T_eff(r) profiles are shown in the left column of Fig. 3. HP Lyr’s #MDW 3 model has a high inner rim radius at 20 R₂. As a result, despite the high accretion rate of 10⁻³ M_⊙/yr, its profile only reaches up to ≈5000 K. TW Cam’s #MDW 2 model fits a disc that both reaches the stellar surface and has a high accretion rate of 10⁻³ M_⊙/yr. Correspondingly, it reaches very high temperatures in the inner disc regions of up to ≈50 000 K.

Fig. 3. Example effective disc temperature profiles and spectra. Top: HP Lyr’s #MDW 3 model. Bottom: TW Cam’s #MDW 2 model. Left: Teff(r) profile (Eq. (2)). Dotted black lines indicate the disc boundaries. Middle: the model disc spectrum (Eq. (3)), normalised to the peak flux, shown as a black line. A scaled blackbody fitted to the whole spectrum is drawn as a red dot-dashed line. The H band location is shown in gray. Right: zoom-in of the spectrum in H band. The values at the six PIONIER wavelength channels are shown as black symbols. The blue dashed line shows a scaled blackbody fitted to these points only. Fitted blackbody temperatures for the full spectrum ($T_{\mathrm{full}}^{\mathrm{BB}}$) and H band only ($T_{\mathrm{H}}^{\mathrm{BB}}$) are shown.

6.5.1. The IRAS 08544–4431 accretion disc temperature estimate

Using H band interferometric data from the Very Large Telescope Interferometer/Precision Integrated-Optics Near-infrared Imaging ExpeRiment instrument (VLTI/PIONIER), Hillen et al. (2016) singled out the spatially unresolved flux contribution from the pole-on accretion disc in IRAS 08544−4431. This is one of the most well-studied post-AGB binaries, in which all components have been spatially resolved. Fitting the spectral slope with a scaled blackbody, they derived an effective blackbody temperature T_disc = 4000 ± 2000 K. This is currently the only direct observational temperature estimate known for a post-AGB binary accretion disc. Assuming this value as representative for the general population, Paper I concluded that the inner regions of their #MDW 1 model were too hot, and that the accretion discs should be fairly cold in general. However, by applying a similar scaled blackbody fit to our models, we show below that the H band interferometric estimate is inherently biased, and unsuitable for comparison with disc wind models without broader wavelength coverage of the data.

6.5.2. Spectral blackbody temperatures

Ignoring possible circum- or interstellar reddening, the relative flux scaling of the accretion disc with wavelength follows:

$$\begin{array}{c}\hfill L_{\lambda }\propto {\displaystyle {\int }_{r_{\mathrm{in}}}^{r_{\mathrm{out}}}}r·B_{\lambda }(T_{\mathrm{eff}}(r))\mathrm{d}r,\end{array}$$(3)

with B_λ the Planck function. The resulting spectra for our example models are shown in the middle column of Fig. 3. Scaled blackbodies were fitted to the full spectra to derive a blackbody temperature, $T_{\mathrm{full}}^{\mathrm{BB}}$. We note that this fits poorly to the multi-blackbody model spectra, which generally have a flattened peak and broader wings. $T_{\mathrm{full}}^{\mathrm{BB}}$ just serves as a rough but useful mean spectral temperature for the full disc spectrum. For the HP Lyr #MDW 3 model $T_{\mathrm{full}}^{\mathrm{BB}}=3600\mathrm{K}$, while for the TW Cam #MDW 2 model $T_{\mathrm{full}}^{\mathrm{BB}}=34300\mathrm{K}$. Indeed, due to the influence of the colder and geometrically larger outer disc regions, this is lower than the maximum achieved disc temperatures. Nevertheless, the TW Cam #MDW 2 model seems too hot to comply with the IRAS 08544−4431 estimate of T_disc = 4000 ± 2000 K.

The analysis above, however, neglects a crucial aspect. The PIONIER data used by Hillen et al. (2016) only covers the H band. We again fitted a scaled blackbody, but now only to the PIONIER wavelength channels, deriving a new spectral temperature, $T_{\mathrm{H}}^{\mathrm{BB}}$. The resulting fits are shown in the right column of Fig. 3. For HP Lyr #MDW 3, we derive $T_{\mathrm{H}}^{\mathrm{BB}}=3000\mathrm{K}$, agreeing well with the $T_{\mathrm{full}}^{\mathrm{BB}}$ value. For TW Cam’s #MDW 2 model, however, we fit $T_{\mathrm{H}}^{\mathrm{BB}}=4300\mathrm{K}$. This is much lower than $T_{\mathrm{full}}^{\mathrm{BB}}$, and suddenly agrees very well with the observational estimate for IRAS 08544−4431. As seen from Fig. 4, this sudden drop in derived spectral temperature is consistent for all hot disc models, which typically extend close to the stellar surface. In comparison, HP Lyr’s models are all cold due to the high fitted inner radii (≥20 R₂), and their spectral backbody temperature estimates are consistent between using the full spectrum or the H band only. Looking at the lower plot in the middle column of Fig. 3, the reason for this drop from $T_{\mathrm{full}}^{\mathrm{BB}}$ to $T_{\mathrm{H}}^{\mathrm{BB}}$ for the hot discs is elucidated. The peak disc flux for TW Cam’s #MDW 2 model is strongly blueshifted due to the hot inner disc regions. The H band itself lies far along the long-wavelength wing of the spectrum, and is dominated by the cold, large outer disc regions. The spectral slope is only weakly sensitive to the hottest parts of the disc in this narrow wavelength range. As a result, fitting the wavelength dependency of a hot disc’s relative H band spectrum with a scaled blackbody results in a spectral blackbody temperature estimate biased to a few 1000 K. The method presented in Hillen et al. (2016) is not a good probe of the disc’s thermal structure.

Fig. 4. Blackbody temperatures derived from the accretion disc’s spectral slope, considering either the full spectrum ($T_{\mathrm{full}}^{\mathrm{BB}}$, red diamonds), or only the H band ($T_{\mathrm{H}}^{\mathrm{BB}}$, blue dots). The orange shaded area denotes the ±1σ region of the Tdisc estimate for IRAS 08544−4431 by Hillen et al. (2016), around which the $T_{\mathrm{H}}^{\mathrm{BB}}$ values cluster.

Unlike what was posited by Paper I, hot accretion disc models thus do not contradict current near-infrared interferometric observations. If $T_{\mathrm{full}}^{\mathrm{BB}}\gtrsim 10000\mathrm{K}$, the model disc spectrum sharply peaks in the UV at ∼200 nm, blueward of the Balmer jump. If the total flux contribution of the accretion disc is significant, this should result in observable evidence for UV radiation. Interestingly, there are indeed other observational indicators for a hot component emitting in the UV. Several circumbinary disc-bearing post-AGB binaries show signatures of UV-excited PAH emission and/or a compact photodissociation region (PDR; e.g. Gielen et al. 2011; Bujarrabal et al. 2023). In the case of one system, the Red Rectangle, this has previously been attributed to a hot jet-launching accretion disc (Witt et al. 2009). Our target HD 52961 also shows PAH emission (Gielen et al. 2011), and our best fit does indeed predict a hot disc at $T_{\mathrm{full}}^{\mathrm{BB}}\approx 18500\mathrm{K}$.

We argue that it is very likely that several post-AGB binaries do host a hot accretion disc extending down close to the secondary’s surface, which can act as a source of UV photons. The detection of PAH emission and/or PDRs can already prove suggestive. However, only combined multi-wavelength interferometric observations in the near-infrared and other wavelength regimes, for example combining VLTI/PIONIER and Center for High Angular Resolution Astronomy/Stellar Parameters and Images with a Cophased Array (CHARA/SPICA) observations, can directly and reliably probe the temperature structure of the accretion discs on a system-by-system basis. This, in turn, will place very strong constraints on our disc wind models.

6.6. Trends in preferred MDW solution

From Table 4 and our discussion of the results in Sect. 5, it is clear that the targets show a clear preference for high efficiency MDW solutions, with the #MDW 2, 3, 4 and 5 fits (ξ > 0.10), depending on the target, mitigating or solving multiple issues faced by #MDW 1 (ξ > 0.04) in matching the observed Hα profiles or the expected circumbinary disc lifetime. Such high efficiencies are only possible, at least for cold disc winds, at low magnetizations of μ ∼ 10⁻³ − 10⁻². In addition, there is a strong preference for thicker discs at ϵ = 0.1 in order to reduce the effect of jet rotation on the model spectra. Two targets imply the need for even thicker discs at say ϵ ≳ 0.3 (Sect. 6.1). At such thicknesses, the driving of a cold, purely magnetocentrifugal wind starts to fail (Casse & Ferreira 2000a), and a thermal push from the primary’s irradiation will indeed be necessary to create a bona fide disc wind. If the irradiation is also able to consistently heat the accretion disc midplane and puff up the disc, this would physically justify the high aspect ratios in the first place. In addition, such magnetothermal winds can reach ξ > 0.5 and might fully solve the persisting issue regarding the circumbinary disc lifetime. One might remark that our most efficient solution, #MDW 5, already faced problems, as it systematically fitted unrealistically low inclinations for three of our targets (see Sect. 5). The link between wind kinematics and density profiles will however be very different in an irradiated MDW solution.

7. Summary and conclusions

We have successfully developed a fitting routine to fit the predictions of disc wind models, launched from a circumcompanion accretion disc and based on self-similar MHD disc wind solutions, to time-resolved observations of the Hα jet features in post-AGB binaries. The five used MHD solutions cover the range of weakly magnetised and highly efficient cold disc winds, while the five target systems form a fairly representative sample of the known jet-launching post-AGB binaries. The models manage to reproduce many of the observed features fairly well, showing that MHD disc winds can indeed for a large part explain the Hα features seen in post-AGB binaries. Several issues and points of interest remain, however.

Many of the models overestimate the wind’s rotation. This manifests as a redshifted absorption feature which is not observed, implying the need for more slowly rotating solutions, such as those with thicker accretion discs. This issue is especially pertinent for two targets: HD 52961 and BD+46°442.

The efficiency range covered by known cold magnetocentrifugal disc winds, topping out at ξ ≈ 0.30 for discs of aspect ratio ϵ = 0.1, is insufficient. Even at an efficiency of ξ = 0.305, predicted accretion rates may be as high as 10⁻⁴ M_⊙/yr, and correspondingly, predicted lifetimes for the dusty circumbinary disc as short as 100 yr. Given the omnipresence of these circumbinary discs and a typical post-AGB lifetime on the order of 10 000 yr, these accretion rates need to be pushed down further by about two orders of magnitude. This would require disc wind solutions with ξ ≳ 0.5, which are unachievable using pure magnetocentrifugal driving, but can be reached when thermal effects help lift up disc material (Casse & Ferreira 2000b).

High efficiency MHD solutions predict inner rims closer to the secondary’s stellar surface, and four out of five targets are indeed best fitted with an inner rim at a radius of at most five times the secondary’s stellar radius. At the current high fitted accretion rates and assuming magnetospheric funnelling is responsible for truncating the accretion discs, the required magnetic field strengths of the secondaries should be greater than about several kG to 15 kG. This is much higher than what can be generated by an isolated mature main sequence star. At reduced accretion rates and under possible rejuvenation of the stellar dynamo due to accretion, magnetospheric truncation close to the stellar surface remains plausible. HP Lyr stands out as an exception, consistently requiring a fitted inner rim radius equal to or greater than twenty times the stellar radius. This is impossible to achieve via magnetospheric funneling, even for lower accretion rates, requiring unreasonably high field strengths. This large inner radius could be an artefact due to the strong assumption of homogeneous wind temperature. Future modelling attempts should include more realistic and physically inspired temperature structures.

All BD+46°442 models show a consistent lack of low RV absorption close to superior conjunction. This is shown to be, at least in principle, resolvable via the inclusion of an additional low efficiency outflow component inside the empty cavity currently present around the disc wind axis, such as a thermal stellar jet launched from the secondary’s poles (Matt & Pudritz 2005).

A previous method of estimating the accretion disc’s temperature from H band interferometry, performed on IRAS 08544−4431 by Hillen et al. (2016), is shown to be biased to several 1000 K. Hot accretion disc models that reach down to the stellar surface, and can achieve temperatures over 10 000 K in the inner regions, are thus not inconsistent with current H band interferometric constraints. Instead, we argue that several post-AGB binaries that show evidence of UV-excited PAH emission and/or PDRs, such as our target HD 52961, indeed host such a hot accretion disc, providing a source of UV photons. Interferometric data covering a wider wavelength range, for example combining VLTI/PIONIER and CHARA/SPICA, will be required to reliably probe the thermal structure of the accretion discs. This would provide strong constraints on our disc wind models.

There is a clear preference for all targets towards the efficient, thicker MHD solutions with ξ > 0.10 and ϵ = 0.10, which are only achieved at a low magnetization level of μ ∼ 10⁻³ − 10⁻². For future work, we advise the use of warm magnetothermal solutions, which include heat deposition via irradiation from the bright post-AGB primary at the disc surface. This should provide an extra push to the wind material, readily allowing the solutions to reach ξ ≳ 0.5 and for even thicker disc solutions at aspect ratios as high or even higher as ϵ = 0.3 to be explored. This would provide a physically consistent solution to the problems currently faced in regards to overestimated wind rotation and the high required accretion rates.

Acknowledgments

TDP acknowledges support of the Research Foundation – Flanders (FWO) under grant 11P6I24N (Aspirant Fellowship). HVW acknowledges support of the FWO under grant G097619N. OV acknowledges support from the KU Leuven Research Council (grant C16/17/007: MAESTRO). DK acknowledges the support of the Australian Research Council (ARC) Discovery Early Career Research Award (DECRA) grant (DE190100813). DK is also supported in part by the ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. JJ acknowledges support by the NSF AST-2009884 and NASA 80NSSC21K1746 grants. Based on observations obtained with the HERMES spectrograph, which is supported by the FWO, Belgium, the Research Council of KU Leuven, Belgium, the Fonds National de la Recherche Scientifique (F.R.S.-FNRS), Belgium, the Royal Observatory of Belgium, the Observatoire de Genève, Switzerland and the Thüringer Landessternwarte Tautenburg, Germany. The computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Centre), funded by the FWO and the Flemish Government – department EWI. Finally, we thank the anonymous referee for helping improve the clarity of the paper.

References

Appendix A: MDW example structures

Fig. A.1 shows five example disc wind structures, calculated entirely analogously except using the different MDW solutions as input. It showcases the general behaviour of the MDW solutions.

Fig. A.1. Model disc wind structures using different MDWs. Left: hydrogen number density n. Middle: poloidal velocity vp. Right: toroidal/rotational velocity vϕ. From top to bottom, we show #MDW 1, 2, 3, 4 and 5. Calculated assuming the HD 52961 system and i = 70°, rin = 5 R2, rout = 0.99 RRL, 2 and Ṁin = 10−3 M⊙/yr. We note the slight discontinuity at a fixed polar angle in the n and vϕ profiles. This is an artefact from the leapfrog integration performed over the Alfvén surface in the MDW solutions (Jacquemin-Ide et al. 2019), in combination with our projection onto Cartesian coordinates.

Appendix B: Spectral inputs from HERMES time-series

In this appendix, we describe in detail how the different required inputs for the modelling routine are derived from the observed HERMES time-series. We cover the pre-selection of observed spectra in Sect. B.1, the orbital parameters in Sect. B.2, and the background spectrum and its associated uncertainty in Sects. B.3 & B.4, respectively.

B.1. Observed spectra selection

We require a sub-selection of spectra from the full HERMES time-series for the object of interest. In order to limit both cycle-to-cycle variability and computational cost, we select several spectra lying within one to two orbital periods from each other and with a reasonable phase coverage across the orbit. If spectra lie too close to each other in orbital phase (i.e. within a difference of 0.01), only the spectrum with the highest signal-to-noise around H_α is retained. Remaining gaps in the orbital phase, ϕ_orb, are then filled by manually selecting additional spectra. The number of selected spectra for our targets, as well as their average signal-to-noise around H_α, are shown in Table 1. The resulting H_α dynamic spectra for our targets are shown in the left column of Fig. 1.

B.2. Orbital parameters

The fitting routine makes extensive use of the system’s SB1 orbital parameters in order to set te geometry of the binary. These parameters are derived from fitting the primary’s RV curve with a Keplerian orbit while iteratively filtering out the signal due to stellar pulsations (Oomen et al. 2018). These parameters are summarised for our target systems in Table B.1. This does not provide a full orbital solution for the binary. In addition, an assumption needs to be made for both the inclination and the primary’s mass. The former is taken to be a fitting parameter (see Sect. 3.3), while the latter is set to M₁ = 0.6 M_⊙ (typical for single white dwarfs in the Galactic field; e.g. Tremblay et al. 2016). This allows for the calculation of the secondary’s mass, M₂, from the spectroscopic mass function, and solves the full 3D binary orbit.

B.3. Background spectrum

A dynamic background spectrum for the H_α line, valid for the entire orbit, needs to be constructed as input for the RT module. This background needs to represent the observed, unobscured spectra out of superior conjunction, which are unaffected by jet absorption. The construction procedure takes several observed unobscured spectra, and subtracts the RV shifted photospheric absorption component of the primary star. The photospheric component is represented using a suitable stellar spectrum from Coelho (2014). This leaves only the emission component, which is averaged over the spectra to create a template for the emission. This template is then added back on top of the photospheric absorption. If the observed emission component moves with the primary’s RV curve in the observed time-series, the template is RV shifted according to ϕ_orb before addition. If the emission is centred on the systemic velocity, only the photospheric component is RV shifted before addition. This results in a suitable dynamic background spectrum, which is assumed to be fully localised on the primary’s surface within the RT module. The background spectra constructed for our targets are shown in the middle column of Fig. 1.

B.4. Additional background uncertainty

The resulting dynamic background spectra are, by construction, smooth. However, looking at the observations in the left column of Fig. 1, we see significant variability in spectra outside of superior conjunction, unrelated to variable jet obscuration. This variability is physical, as most post-AGB stars are photometrically variable (e.g. Van Winckel 2003). As such, there is an inherent uncertainty associated to the constructed background.

In order to account for this inherent background variability, we include an additional error term, σ_back(λ). This term should be higher in the inherently more variable regions of the spectra. Correspondingly, the model predictions will be punished less severely there. We take several spectra outside of superior conjunction, and calculate the difference with the constructed background spectrum over the considered wavelengths. This difference profile is then averaged over the spectra to give σ_back(λ), which is added in quadrature to the S/N error to give the total spectral error ${\sigma }_{\text{tot}}(\lambda )={({\sigma }_{\text{S/N}}^2+{\sigma }_{\text{back}}^2(\lambda ))}^{1/2}$. σ_tot is dominated by σ_back for all our targets. If the H_α emission component follows the primary’s RV, the σ_tot profile is shifted correspondingly. With this additional error, we typically get a reduced chi-squared ${\chi }_{\nu }^2$ ∼ 1 when comparing the background spectrum to the observed unobscured spectra. This shows that the combination of the background spectrum and σ_tot is a good statistical representation for the unobscured spectra. The total spectral error profiles for our targets are shown in Fig. B.1.

Fig. B.1. Total Hα spectral error profiles, accounting for the variability of the observed unobscured spectra.

Appendix C: All disc wind model fits

The fitted parameters, for all five MDW solutions and for every one of our five targets, are summarised in Table C.1. The corresponding model dynamic spectra are shown in Fig. C.1. These fits use the same ${\chi }_{\nu }^2$ calculation boxes that are shown in Fig. 1. We note that the selection of the best fitting MDW model is guided by the ${\chi }_{\nu }^2$ values, but ultimately depends on the ability of the MDW to recover the global H_α pattern.

Fig. C.1. Fitted model dynamic spectra for all targets (Table 1) and considered MDW solutions, ordered left to right by their #MDW number. The corresponding model parameters are summarised in Table C.1.

Appendix D: Mass-loss rate for extra outflow inside the disc wind cavity

In this appendix, we derive a simple expression for the mass-loss rate that an extra outflow inside the disc wind cavity would need in order to supply the missing low RV absorption occurring close to superior conjunction in some of our models. We assumed a simple cylindrical geometry for this extra outflow.

Consider a disc wind model close to superior conjunction, where a chosen LOS ray passes through both the left and right wing of the wind lobe facing the observer (Fig. D.1). In addition, assume that there is some extra component of material flow inside the wind cavity. Neglecting the wind emission, the intensity at wavelength λ leaving the disc wind is then given by:

$$\begin{array}{c}\hfill I_{\lambda }^{\mathit{out}}=I_{\lambda }^0·e^{-\mathrm{\Delta }{\tau }_{\lambda ,tot}}=I_{\lambda }^0·e^{-(\mathrm{\Delta }{\tau }_{\lambda ,left}+\mathrm{\Delta }{\tau }_{\lambda ,cav}+\mathrm{\Delta }{\tau }_{\lambda ,right})},\end{array}$$(D.1)

Fig. D.1. Simple geometrical model for an extra outflow component in the disc wind cavity. The dashed line indicates the stretch of the LOS inside the cavity. It has length s. The extra component is taken to be a homogeneous flow of uniform velocity v, pointed along the jet axis and flowing along a cylinder of radius r = (s/2)⋅sin i.

with $I_{\lambda }^0$ the background spectrum. We have split the total optical depth along the LOS ray into the separate contributions of the two jet wings and the hypothetical material in the cavity. Identifying $I_{\lambda }^{\mathit{out}}$ with the observed spectrum, I_λ, obs, and $I_{\lambda }^0\cdot e^{-(\mathrm{\Delta }{\tau }_{\lambda ,left}+\mathrm{\Delta }{\tau }_{\lambda ,right})}$ as the spectrum modelled by our MDW fits, I_λ, mod, this gives:

$$\begin{array}{c}\hfill I_{\lambda ,obs}=I_{\lambda ,mod}·e^{-\mathrm{\Delta }{\tau }_{\lambda ,cav}}\Rightarrow \mathrm{\Delta }{\tau }_{\lambda ,cav}=ln\left(\frac{1}{1-\frac{\mathrm{\Delta }{I}_{\lambda }}{I_{\lambda ,mod}}}\right),\end{array}$$(D.2)

where ΔI_λ = I_λ, mod − I_λ, obs. Eq. (D.2) thus identifies how much optical depth is missing at a certain wavelength in order to bring the emission peak shining through our models in line with the observations. Fig. D.2 visualises the different quantities involved in the calculation of Eq. (D.2).

Fig. D.2. Example spectrum at ϕorb = 0.49 of BD+46°442’s #MDW 3 fit. The model is shown as a solid red line, while the observations are shown as black dots. The curve of ΔIλ/Iλ, mod is shown as a dashed blue line. The central wavelength value at which absorption is missing (RV ≈ −60 km/s) is indicated by the black dot-dashed line.

With the estimate for the extra required optical depth at hand, we can construct a simple geometrical model for our extra outflow component. This model is visualised in Fig. D.1. Consider the length s of an LOS ray inside the cavity. The contribution of the cavity material to the optical depth is given by Δτ_λ, cav = κ_λ, cav ⋅ ρ_cav ⋅ s, with κ_λ, cav the material’s H_α line opacity profile and ρ_cav its mass density, assuming homogeneity. For κ_λ, cav, we choose the H_α line opacity of the disc wind, at the point where the LOS ray enters the cavity. This mainly assumes that the material in the cavity has the same temperature as the disc wind (T = T_wind), since the H_α opacity is strongly dependent on temperature, but only weakly pressure-dependent at the considered densities. Choosing a range of wavelengths Δλ, spanning the background emission peak that shines through the model, this can then be averaged over Δλ:

$$\begin{array}{cc}\hfill {\overline{\mathrm{\Delta }\tau }}_{\mathit{cav}}& ={\overline{\kappa }}_{\mathit{cav}}·{\rho }_{\mathit{cav}}·s\Rightarrow {\rho }_{\mathit{cav}}=\frac{{\overline{\mathrm{\Delta }\tau }}_{\mathit{cav}}}{{\overline{\kappa }}_{\mathit{cav}}·s},\hfill \end{array}$$(D.3)

$$\begin{array}{cc}\hfill {\overline{\kappa }}_{\mathit{cav}}& =\frac{{\int }_{\mathrm{\Delta }\lambda }{\kappa }_{\lambda ,cav}\mathrm{d}\lambda }{\mathrm{\Delta }\lambda },\hfill \end{array}$$(D.4)

$$\begin{array}{cc}\hfill {\overline{\mathrm{\Delta }\tau }}_{\mathit{cav}}& =\frac{{\int }_{\mathrm{\Delta }\lambda }ln\left(\frac{1}{1-\mathrm{\Delta }{I}_{\lambda }/I_{\lambda ,mod}}\right)\mathrm{d}\lambda }{\mathrm{\Delta }\lambda }·\hfill \end{array}$$(D.5)

Eq. (D.3) can thus be used to obtain a single estimate for the density of the cavity material needed to supply the missing absorption. A representative LOS ray needs to be chosen, for which we chose the ray leaving the centre of the post-AGB primary.

Assume that the extra component consists of an outflow along a cylinder of radius r = (s/2)⋅sin i, with a uniform velocity v along the jet axis (Fig. D.1). The velocity v is related to the observed central RV value of the missing absorption via v = |RV|/cos i (i.e. we assume all the missing absorption to be concentrated at the central RV value). Using Eq. (D.3) to calculate ρ_cav, we retrieve an estimate of the extra component’s mass-loss rate:

$$\begin{array}{c}\hfill {\dot{M}}_{\mathit{cav}}\sim \pi r^2·v·{\rho }_{\mathit{cav}}=\frac{\pi {sin}^{2}i}{4cosi}·\frac{\mathrm{RV}·\mathrm{s}}{{\overline{\kappa }}_{\mathit{cav}}}·{\overline{\mathrm{\Delta }\tau }}_{\mathit{cav}}.\end{array}$$(D.6)

Typical values for models that show a clear lack of wind absorption near superior conjunction, for example #MDW 3 for BD+46°442 and IRAS 19135+3937, are ${\overline{\kappa }}_{\mathit{cav}}=10-100{\mathrm{cm}}^2/\mathrm{g}$, ρ_cav ≈ 10⁻¹⁵ g/cm³, v ≈ 100 km/s and Ṁ_cav ≈ 10⁻⁸ M_⊙/yr. These are the physical conditions that an additional outflow inside the cavity, for example a thermal stellar jet, would need to achieve in order to provide the missing absorption.

All Tables

All Figures

Fig. 1. Dynamic spectra (see explanation in Sect. 3.1) for the target systems’ Hα lines. Left: observed spectral time-series after spectra selection. Middle: background spectrum used as input for the RT module. Right: best fitting disc wind model (parameters summarised in Table 4). A black dashed line marks superior conjunction, while a red selection box denotes the region used for the calculation of ${\chi }_{\nu }^2$.

In the text

Fig. 2. Example disc wind structures using two different MDWs. Top: using #MDW 1. Bottom: using #MDW 5. Left: hydrogen number density n. Middle: poloidal velocity vp. Right: toroidal/rotational velocity vϕ. Calculated assuming the HD 52961 system and i = 70°, rin = 5 R2, rout = 0.99 RRL, 2 and Ṁin = 10−3 M⊙/yr (see Sect. 3). We note the slight discontinuity at a fixed polar angle in the n and vϕ profile. This is an artefact from the leapfrog integration performed over the Alfvén surface in the MDW solutions (Jacquemin-Ide et al. 2019), in combination with our projection onto Cartesian coordinates. This introduces no uncertainty on our results.

In the text

Fig. 3. Example effective disc temperature profiles and spectra. Top: HP Lyr’s #MDW 3 model. Bottom: TW Cam’s #MDW 2 model. Left: Teff(r) profile (Eq. (2)). Dotted black lines indicate the disc boundaries. Middle: the model disc spectrum (Eq. (3)), normalised to the peak flux, shown as a black line. A scaled blackbody fitted to the whole spectrum is drawn as a red dot-dashed line. The H band location is shown in gray. Right: zoom-in of the spectrum in H band. The values at the six PIONIER wavelength channels are shown as black symbols. The blue dashed line shows a scaled blackbody fitted to these points only. Fitted blackbody temperatures for the full spectrum ($T_{\mathrm{full}}^{\mathrm{BB}}$) and H band only ($T_{\mathrm{H}}^{\mathrm{BB}}$) are shown.

In the text

Fig. 4. Blackbody temperatures derived from the accretion disc’s spectral slope, considering either the full spectrum ($T_{\mathrm{full}}^{\mathrm{BB}}$, red diamonds), or only the H band ($T_{\mathrm{H}}^{\mathrm{BB}}$, blue dots). The orange shaded area denotes the ±1σ region of the Tdisc estimate for IRAS 08544−4431 by Hillen et al. (2016), around which the $T_{\mathrm{H}}^{\mathrm{BB}}$ values cluster.

In the text

Fig. A.1. Model disc wind structures using different MDWs. Left: hydrogen number density n. Middle: poloidal velocity vp. Right: toroidal/rotational velocity vϕ. From top to bottom, we show #MDW 1, 2, 3, 4 and 5. Calculated assuming the HD 52961 system and i = 70°, rin = 5 R2, rout = 0.99 RRL, 2 and Ṁin = 10−3 M⊙/yr. We note the slight discontinuity at a fixed polar angle in the n and vϕ profiles. This is an artefact from the leapfrog integration performed over the Alfvén surface in the MDW solutions (Jacquemin-Ide et al. 2019), in combination with our projection onto Cartesian coordinates.

In the text

	Fig. B.1. Total Hα spectral error profiles, accounting for the variability of the observed unobscured spectra.
In the text

	Fig. C.1. Fitted model dynamic spectra for all targets (Table 1) and considered MDW solutions, ordered left to right by their #MDW number. The corresponding model parameters are summarised in Table C.1.
In the text

	Fig. D.1. Simple geometrical model for an extra outflow component in the disc wind cavity. The dashed line indicates the stretch of the LOS inside the cavity. It has length s. The extra component is taken to be a homogeneous flow of uniform velocity v, pointed along the jet axis and flowing along a cylinder of radius r = (s/2)⋅sin i.
In the text

	Fig. D.2. Example spectrum at ϕorb = 0.49 of BD+46°442’s #MDW 3 fit. The model is shown as a solid red line, while the observations are shown as black dots. The curve of ΔIλ/Iλ, mod is shown as a dashed blue line. The central wavelength value at which absorption is missing (RV ≈ −60 km/s) is indicated by the black dot-dashed line.
In the text