# Jet formation in post-AGB binaries: Confronting cold magnetohydrodynamic disc wind wind models with observations

Toon de Prins, Hans van Winckel, Jonathan Ferreira, Olivier Verhamme, Devika Kamath, Nathan Zimniak, Jonathan Jacquemin-Ide

► To cite this version:

Toon de Prins, Hans van Winckel, Jonathan Ferreira, Olivier Verhamme, Devika Kamath, et al.. Jet formation in post-AGB binaries: Confronting cold magnetohydrodynamic disc wind wind models with observations. Astronomy & Astrophysics - A&A, 2024, 689, pp.A151. 10.1051/0004-6361/202449631. hal-04641851

HAL Id: hal-04641851

https://hal.science/hal-04641851v1

Submitted on 26 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Jet formation in post-AGB binaries

## Confronting cold magnetohydrodynamic disc wind models with observations*

T. De Prins<fnref n="1" />,**, D. H. Van Winckel<fnref n="1" />, D. J. Ferreira<fnref n="2" />, O. Verhamme<fnref n="1" />, D. Kamath<fnref n="3" />,<fnref n="4" />, N. Zimniak<fnref n="2" />, and J. Jacquemin-Ide<fnref n="5" />

Received 16 February 2024 / Accepted 12 June 2024

## ABSTRACT

**Context.** Jets are launched from many classes of astrophysical objects, including post-asymptotic giant branch (post-AGB) binaries with a circumbinary disc. Despite dozens of detections, the formation of these post-AGB binary jets and their connection to the inter-component interactions in their host systems remains poorly understood.

**Aims.** Building upon the previous paper in this series, we consider cold self-similar magnetohydrodynamic (MHD) disc wind solutions to describe jets that are launched from the circumcompanion accretion discs in post-AGB binaries. Resulting predictions are matched to observations. This both tests the physical validity of the MHD disc wind paradigm and reveals the accretion disc properties.

**Methods.** Five MHD solutions are used as input to synthesise spectral time-series of the Hα line for five different post-AGB binaries. A fitting routine over the remaining model parameters is developed to find the disc wind models that best fit the observed time-series.

**Results.** Many of the time-series' properties are reproduced well by the models, though systematic mismatches, such as overestimated rotation, remain. Four targets imply accretion discs that reach close to the secondary's stellar surface, while one is fitted with an unrealistically large inner radius at ≥20 stellar radii. Some fits imply inner disc temperatures over 10 000 K, seemingly discrepant with a previous observational estimate from H band interferometry. This estimate is, however, shown to be biased. Fitted mass-accretion rates range from ~10⁻⁶–10⁻<fnref n="3" /> M☉/yr. Relative to the jets launched from young stellar objects (YSOs), all targets prefer winds with higher ejection efficiencies, lower magnetizations and thicker discs.

**Conclusions.** Our models show that current cold MHD disc wind solutions can explain many of the jet-related Hα features seen in post-AGB binaries, though systematic discrepancies remain. This includes, but is not limited to, overestimated rotation and underestimated post-AGB circumbinary disc lifetimes. The consideration of thicker discs and the inclusion of irradiation from the post-AGB primary, leading to warm magnetothermal wind launching, might alleviate these.

**Key words.** accretion, accretion disks – stars: AGB and post-AGB – binaries: spectroscopic – circumstellar matter – ISM: jets and outflows

## 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

<footnote n="1">¹ Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, 3001 Leuven, Belgium</footnote>

<footnote n="2">² Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="3">³ Department of Physics & Astronomy, School of Mathematical and Physical Sciences, Macquarie University, Sydney, NSW 2109, Australia</footnote>

<footnote n="4">⁴ Astronomy, Astrophysics and Astrophotonics Research Centre, Macquarie University, Sydney, NSW 2109, Australia</footnote>

<footnote n="5">⁵ Center for Interdisciplinary Exploration & Research in Astrophysics (CIERA), Physics and Astronomy, Northwestern University, Evanston, IL 60202, USA</footnote>

<footnote>* Based on observations made with the *Mercator* Telescope, operated on the island of La Palma by the Flemish Community, at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias.</footnote>

<footnote>** Corresponding author; toon.deprins@kuleuven.be</footnote>
<table><caption>**Table 1.** Target system properties.</caption><thead><tr><th>Object parameter</th><th>HD 52961</th><th>BD+46°442</th><th>TW Cam</th><th>IRAS 19135+3937</th><th>HP Lyr</th></tr></thead><tbody><tr><td>T<sub>pAGB</sub> (K)</td><td>6000</td><td>6250</td><td>5250</td><td>6250</td><td>6000</td></tr><tr><td>[Fe/H]</td><td>-4.8</td><td>-0.8</td><td>-0.5</td><td>-1.0</td><td>-1.0</td></tr><tr><td>Depletion</td><td>Strong</td><td>None</td><td>None</td><td>None</td><td>Strong</td></tr><tr><td>P (d)</td><td>1288.6 ± 0.3</td><td>140.82 ± 0.02</td><td>662.2 ± 5.3</td><td>126.97 ± 0.08</td><td>1818 ± 80</td></tr><tr><td>e</td><td>0.23 ± 0.01</td><td>0.085 ± 0.005</td><td>0.25 ± 0.04</td><td>0.13 ± 0.03</td><td>0.20 ± 0.04</td></tr><tr><td>L<sub>Gaia</sub> (L<sub>☉</sub>)</td><td>3300<sup>+1600</sup><sub>-400</sub></td><td>2400<sup>+1100</sup><sub>-1100</sub></td><td>1500<sup>+400</sup><sub>-500</sub></td><td>3600<sup>+700</sup><sub>-1200</sub></td><td>11600<sup>+12100</sup><sub>-6200</sub></td></tr><tr><td>i<sub>par</sub> (°)</td><td>70.6 ± 0.6</td><td>49.6 ± 0.6</td><td>25.7 ± 0.2</td><td>71.8 ± 0.5</td><td>18.6 ± 0.7</td></tr><tr><td>R<sub>1,par</sub> (R<sub>☉</sub>)</td><td>131 ± 4</td><td>27.3 ± 0.6</td><td>83.0 ± 0.6</td><td>22.6 ± 0.4</td><td>75 ± 1</td></tr><tr><td>N<sub>obs</sub></td><td>42</td><td>40</td><td>36</td><td>32</td><td>40</td></tr><tr><td>(S/N)<sub>avg</sub></td><td>110</td><td>70</td><td>70</td><td>40</td><td>60</td></tr><tr><td>Refs.</td><td>1, 2, 6, 7, 10</td><td>3, 6, 7, 8</td><td>4, 6, 7, 9</td><td>2, 6, 7, 8</td><td>5, 6, 7, 9</td></tr></tbody></table>

**Notes.** Relevant properties of the target systems studied in this work, including the primary's effective temperature *T*<sub>pAGB</sub>, metallicity [Fe/H] and refractory depletion strength. In addition, we list the binary period *P* and eccentricity *e*, the *Gaia* DR3 luminosity *L*<sub>Gaia</sub>, and the inclination *i*<sub>par</sub> and primary star radius *R*<sub>1,par</sub> found from previous parametric modelling. Below that, the number of selected observed spectra *N*<sub>obs</sub>, and their average signal-to-noise ratio (S/N)<sub>avg</sub>, as well as the references used. *T*<sub>pAGB</sub> and [Fe/H] values (formal errors of 250 K and 0.3 dex, respectively) were based on (1) Waelkens et al. (1991), (2) Rao & Giridhar (2014), (3) Gorlova et al. (2012), (4) Giridhar et al. (2000) and (5) Giridhar et al. (2005). The refractory depletion strength classification, as well as *P* and *e* were taken from (6) Oomen et al. (2018). *L*<sub>Gaia</sub> was derived from (7) Gaia Collaboration (2023), *i*<sub>par</sub> and *R*<sub>1,par</sub> from (8) Bollen et al. (2020), (9) Bollen et al. (2021) and (10) Bollen et al. (2022). We note that the luminosities are not fully reliable, as none of these systems were flagged as astrometric binaries in *Gaia* DR3.

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.

Still, parametric prescriptions provide little connection to the actual underlying physics of the accretion disc and the
<figure><img src="image_2.png" /><figcaption>**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_v^2$.</figcaption></figure>
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\alpha$ 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 $\epsilon$, disc magnetization $\mu$, disc turbulence parameter $\alpha_m$, the ejection efficiency $\xi$ and the radiative efficiency $\eta_{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

<table><caption>**Table 2.** Dimensionless MDW parameters.</caption><thead><tr><th>#MDW</th><th>ε</th><th>μ</th><th>α<sub>m</sub></th><th>ξ</th><th>η<sub>rad</sub></th></tr></thead><tbody><tr><td>1</td><td>0.03</td><td>0.351</td><td>2</td><td>0.04</td><td>0.118</td></tr><tr><td>2</td><td>0.1</td><td>2.85 × 10<sup>-3</sup></td><td>1</td><td>0.103</td><td>0.545</td></tr><tr><td>3</td><td>0.1</td><td>1.53 × 10<sup>-2</sup></td><td>1</td><td>0.104</td><td>0.445</td></tr><tr><td>4</td><td>0.1</td><td>2.3 × 10<sup>-3</sup></td><td>1</td><td>0.215</td><td>0.364</td></tr><tr><td>5</td><td>0.1</td><td>2.2 × 10<sup>-3</sup></td><td>2</td><td>0.305</td><td>0.325</td></tr></tbody></table>

**Notes.** Relevant dimensionless parameters of the self-similar MDW solutions considered in this work. The solutions are labelled by number (\#MDW) and ordered by ascending ejection efficiency $\xi$. See text Sect. 2.1 for these parameters' definitions.

the general effects of the parameters on the wind. An analogous figure for all five MDWs is presented in Appendix A.

The scale height of the disc is parametrised by $h(r) = \epsilon r$, with $r$ the radius along the midplane and $\epsilon$ the disc aspect ratio. Higher $\epsilon$ values imply thicker accretion discs. \#MDW 1 is a very thin disc at $\epsilon = 0.03$, while all other MDWs we consider are considerably thicker at $\epsilon = 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 $\mu = B^2/(\mu_0 P)$, with $B$ the vertical field strength, $\mu_0$ 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 $\mu$ increases the dynamical importance of the magnetic field. Our MDWs range from $\mu \sim 10^{-1}-10^{-3}$, covering almost the entire $\mu$-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 $\nu_v$, and magnetic diffusivity $\nu_m$. These are assumed to follow a vertical Gaussian profile, and have respective midplane amplitudes $\nu_v = \alpha_v C_S h$ and $\nu_m = \alpha_m V_A h$, with $\alpha_v$ the standard dimensionless Shakura-Sunyaev parameter for viscous angular momentum transport (Shakura & Sunyaev 1973), $\alpha_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 $\mu = (V_A/C_S)^2$). These are related via the magnetic Prandtl number, $P_m = \nu_v/\nu_m$. The diffusivity related to the toroidal magnetic field is then given by $v'_m = \nu_m/\chi_m$, where $\chi_m$ is the diffusivity anisotropy parameter. The MDWs we consider here assume $P_m = 1$ and $\chi_m \approx 1$, consistent with simulations of turbulence induced by the magneto-rotational instability (e.g. Fromang & Stone 2009; Lesur & Longaretti 2009; Zhu & Stone 2018). Setting $\alpha_m$ then sets all turbulence-related terms. We note that $\alpha_v$ scales as $\alpha_v = \alpha_m P_m \mu^{1/2}$ (Salvesen et al. 2016). Increasing $\mu$ thus increases both the turbulent and the laminar (wind) torques. The MDWs considered here assume either $\alpha_m = 1$ or 2.

The disc accretion rate is parametrised as $\dot{M}_{acc}(r) \propto r^\xi$, with $\xi$ 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 $\dot{M}_{wind} = 0.5 \cdot \dot{M}_{in}((r_{out}/r_{in})^\xi - 1)$, with $\dot{M}_{in}$ the accretion rate at the inner disc rim, and $r_{out}$ and $r_{in}$ the outer and inner
<figure><img src="image_3.png" /><figcaption>**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 $v_p$. Right: toroidal/rotational velocity $v_\phi$. Calculated assuming the HD 52961 system and $i = 70^\circ$, $r_{in} = 5 R_2$, $r_{out} = 0.99 R_{\text{RL,2}}$ and $\dot{M}_{\text{in}} = 10^{-3} \text{ M}_\odot/\text{yr}$ (see Sect. 3). We note the slight discontinuity at a fixed polar angle in the $n$ and $v_\phi$ 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.</figcaption></figure>

radius of the wind-launching disc (these are fitting parameters, see Sect. 3.3). Higher $\xi$ values imply that more material is deviated to the wind, making the wind denser and increasing the wind absorption depth in $H\alpha$ (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 $\xi$ '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_0$, respectively follow $v_p^\infty = r_0\Omega_{K_0}(2\lambda - 3)^{1/2}$ and $rv_\phi^\infty = r_0^2\Omega_{K_0}\lambda$, with $\Omega_{K_0}$ the Keplerian angular velocity at $r_0$ and $\lambda$ the magnetic lever arm (Blandford & Payne 1982). It can be shown under certain assumptions that $\lambda \approx 1 + 1/(2\xi)$ (Ferreira 1997). Deviations from this expression are small for all MDWs considered here (see Fig. 5 in Jacquemin-Ide et al. 2019). Increasing $\xi$ thus also leads to a decrease in the lever arm $\lambda$, 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_{\text{out}}$ and $r_{\text{in}}$, is written as $P_{\text{acc}} = 2P_{\text{wind}} + P_{\text{rad}} + P_{\text{adv}}$, with $P_{\text{acc}}$ the released accretion power, $P_{\text{wind}}$ the power fed into one wind lobe, $P_{\text{rad}}$ the power released as radiation from the disc and $P_{\text{adv}}$ the thermal power

carried further inward below $r_{\text{in}}$ by the material accreting onto the central object (Ferreira 1997). Dividing this expression by $P_{\text{acc}}$ leads to an expression for the disc radiative efficiency $\eta_{\text{rad}} = P_{\text{rad}}/P_{\text{acc}} = 1 - 2\eta_{\text{wind}} - \eta_{\text{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 $\eta_{\text{rad}}$ and the disc aspect ratio $\epsilon$ would be an outcome of the calculations. Here, since $\epsilon$ is a prescribed parameter of the MDW solutions, $\eta_{\text{rad}}$ is found a posteriori using the expression above, since $\eta_{\text{wind}}$ and $\eta_{\text{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 $\eta_{\text{adv}}$ ($\eta_{\text{adv}} \propto \epsilon^2$). One might thus expect discs with a higher aspect ratio $\epsilon$ to be less bright and have a lower $\eta_{\text{rad}}$. While this scaling remains correct for our MDWs, the variations in $\eta_{\text{rad}}$ among them (see Table 2) are mostly due to variations in the fraction of power carried away in the wind lobes $\eta_{\text{wind}}$. The disc temperature is itself a direct function of the disc accretion rate and $\eta_{\text{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 $\alpha_m = 1$ or 2, respectively.

## 3. Methods

This section summarises the methodology for modelling and fitting synthetic $H\alpha$ 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\alpha$ 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 $\phi_{\text{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\alpha$ 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 $\phi_{\text{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\alpha$ 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 $\le 0.05$ AU and 0.1 AU, for the directions perpendicular to and along the wind axis respectively, were achieved for all models.

<table><caption>**Table 3.** Fitting routine parameter grid.</caption><thead><tr><th>Parameter</th><th>Grid values</th></tr></thead><tbody><tr><td><i>i</i> (°)</td><td><i>i</i><sub>0</sub><sup>(a)</sup>, <i>i</i><sub>0</sub> + 4, ..., <i>i</i><sub>0</sub> + 60</td></tr><tr><td><i>r</i><sub>in</sub> (<i>R</i><sub>2</sub>)</td><td>1, 5, 10, 15, 20, 25</td></tr><tr><td><i>r</i><sub>out</sub> (<i>R</i><sub>RL,2</sub>)</td><td>0.05, 0.1, 0.2, ..., 0.9</td></tr><tr><td>log (<i>M</i><sub>in</sub>) (M<sub>☉</sub>/yr)</td><td>-8, -7, -6, -5, -4, -3</td></tr><tr><td><i>T</i><sub>wind</sub> (K)</td><td><i>T</i><sub>pAGB</sub>, <i>T</i><sub>pAGB</sub> - 200, ..., <i>T</i><sub>pAGB</sub> - 2000</td></tr><tr><td><i>R</i><sub>1</sub> (R<sub>☉</sub>)</td><td>Fixed <sup>(b)</sup></td></tr></tbody></table>

**Notes.** An ellipsis (...) indicates that the stepsize is maintained for the rest of the grid values. <sup>(a)</sup> $i_0 = 20°$ for all targets except HP Lyr, for which $i_0 = 5°$ was chosen due to the low inclination implied from previous modelling (see Table 1). <sup>(b)</sup>Fixed to the values derived from previous modelling by Bollen et al. (2020, 2021, 2022) (see Table 1).

### 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^2$, within a box in the dynamic spectrum. This box was chosen via a range in $\phi_{\text{orb}}$ and RV, denoting the region of points used for the $\chi^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_2$. This also increases the secondary's radius, $R_2$ (calculated from the MS mass-radius relationship of Demircan & Kahraman 1991), and widens the coverage of the absorption feature in $\phi_{\text{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_{\text{in}}$ ($R_{\text{RL,2}}$).** 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_{\text{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_2$, and is limited by the secondary's stellar surface, $r_{\text{in}} \ge 1 R_2$.

**Outer launching radius, $r_{\text{out}}$ ($R_{\text{RL,2}}$).** The outer launching radius of the disc wind. Its field line is used as the outer boundary of the wind. Increasing $r_{\text{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 $\phi_{\text{orb}}$ coverage. It is expressed in units of the secondary's Roche radius, $R_{\text{RL,2}}$. It is limited by the inner launching radius and the Roche radius $R_{\text{RL,2}} > r_{\text{out}} > r_{\text{in}}$. $r_{\text{out}} \ge R_{\text{RL,2}}$ would namely imply that the circumcompanion accretion disc spills over onto the primary.

**Inner rim accretion rate, log ($\dot{M}_{\text{in}}$) (M<sub>☉</sub>/yr).** The rate of material going through the inner edge of the accretion disc at $r_{\text{in}}$, and falling further onto the star. Disc wind densities increase proportionally to $\dot{M}_{\text{in}}$ (e.g. Jacquemin-Ide et al. 2019), and so
does the total rate of mass lost to the wind lobes: $\dot{M}_{\text{wind}} = 0.5 \cdot \dot{M}_{\text{in}}((r_{\text{out}}/r_{\text{in}})^{\xi} - 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 $\dot{M}_{\text{out}} = (r_{\text{out}}/r_{\text{in}})^{\xi} \dot{M}_{\text{in}}$, expressing the feeding rate in the re-accretion stream from the circumbinary disc onto the accretion disc.

**Disc wind temperature, $T_{\text{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$\alpha$ absorption. This effect is limited by the wind's thermal blackbody emission, which will eventually dominate over the increase in absorption as $T_{\text{wind}}$ increases. Since the wind feature is seen in absorption, we are limited by $T_{\text{wind}} < T_{\text{pAGB}}$, with $T_{\text{pAGB}}$ the temperature of the primary (see Table 1), as the feature would appear in emission otherwise.

**Primary radius, $R_1$ (R$_{\odot}$).** The radius of the primary post-AGB star. Increasing it increases the $\phi_{\text{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,\text{par}}$, see Table 1). $R_1$ 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$\alpha$ 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 $\chi^2$ 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$\alpha$ spectra. Table 1 summarises the properties of these systems used for their selection and modelling. The observed dynamic H$\alpha$ 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$\alpha$ 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 \approx 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$^{60}$ 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$\alpha$ spectrum presents a deep double-peaked jet absorption feature slightly offset from conjunction, covering $\phi_{\text{orb}} \approx 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 $\phi_{\text{orb}} \approx 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_{\text{par}} \approx 70^\circ$, 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 \approx 140$ d and $e \approx 0.1$. Previous parametric modelling implied an intermediate inclination of $i_{\text{par}} \approx 50^\circ$. It does not show strong pulsations. Its dynamic H$\alpha$ spectrum shows a jet absorption feature down to RV ≈ -250 km/s, slightly offset from conjunction. The feature spans $\phi_{\text{orb}} \approx 0.4–0.7$. The background emission wings stay centered 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 \approx 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} \approx 26^\circ$. 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 $\approx -300$ km/s. The feature covers almost the entire orbit, but becomes more prominent at $\phi_{orb} \approx 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 \approx 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} \approx 72^\circ$. 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 $\approx -360$ km/s and spanning $\phi_{orb} \approx 0.3-0.7$. In addition, it shows a small amount of redshifted absorption around $\phi_{orb} \approx 0.40$, reaching up to RV $\approx 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 \approx 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} \approx 19^\circ$. 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 $\approx -120$ km/s and spanning $\phi_{orb} \approx 0.15-0.7$. Approximately from $\phi_{orb} \approx 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_v^2 = 5.32$. It reproduces both the general absorption depth, RV extent and the phase positioning of the high $\phi_{orb}$ absorption peak. The lower $\phi_{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 $\phi_{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 $\phi_{orb} \approx 0.20$, the modelled low $\phi_{orb}$ absorption peak is located at $\phi_{orb} \approx 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^\circ$, 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 $\phi_{orb}$ peak. \#MDW 1, 2 and 3 perform significantly worse at $\chi_v^2 \le 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_v^2 = 5.63$, but only manages to shift the low $\phi_{orb}$ peak into the RV position of the photospheric absorption line. We note that \#MDW 5 fits at $i = 56^\circ$, which is significantly lower than for the other MDWs ($i > 70^\circ$). 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^\circ$ seems fairly low, given the system's robust RVb classification.

Despite significant mismatches the $\chi_v^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_v^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_v^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_v^2 = 6.44$. This model reproduces the absorption depth, $\phi_{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 $\phi_{orb} \approx 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_v^2 \le 6.73$. Rotation also blueshifts the lower $\phi_{orb}$ peak. This detaches the wind absorption from the photospheric
<table><caption>**Table 4.** Fitting parameters and $\chi_v^2$ values for the best fitting disc wind models.</caption><thead><tr><th>Object</th><th>#MDW</th><th>i<br>(°)</th><th>r<sub>in</sub><br>(R<sub>2</sub>)</th><th>r<sub>out</sub><br>(R<sub>RL,2</sub>)</th><th>M&#x0307;<sub>in</sub><br>(M<sub>☉</sub>/yr)</th><th>T<sub>wind</sub><br>(K)</th><th>χ<sup>2</sup><sub>v</sub></th><th>M&#x0307;<sub>out</sub><br>(M<sub>☉</sub>/yr)</th><th>M&#x0307;<sub>wind</sub><sup>(a)</sup><br>(M<sub>☉</sub>/yr)</th></tr></thead><tbody><tr><td>HD 52961</td><td>4</td><td>72</td><td>1</td><td>0.9</td><td>10<sup>-4</sup></td><td>4200</td><td>5.32</td><td>3 × 10<sup>-4</sup></td><td>1 × 10<sup>-4</sup></td></tr><tr><td>BD+46°442</td><td>3</td><td>56</td><td>5</td><td>0.8</td><td>10<sup>-5</sup></td><td>5250</td><td>6.44</td><td>1.2 × 10<sup>-5</sup></td><td>1 × 10<sup>-6</sup></td></tr><tr><td>TW Cam</td><td>3</td><td>24</td><td>5</td><td>0.9</td><td>10<sup>-3</sup></td><td>4850</td><td>5.42</td><td>1.4 × 10<sup>-3</sup></td><td>2 × 10<sup>-4</sup></td></tr><tr><td>IRAS 19135+3937</td><td>4</td><td>48</td><td>1</td><td>0.9</td><td>10<sup>-5</sup></td><td>5050</td><td>3.30</td><td>2.4 × 10<sup>-5</sup></td><td>7 × 10<sup>-6</sup></td></tr><tr><td>HP Lyr</td><td>4</td><td>29</td><td>20</td><td>0.3</td><td>10<sup>-4</sup></td><td>5000</td><td>4.91</td><td>1.4 × 10<sup>-4</sup></td><td>2 × 10<sup>-5</sup></td></tr></tbody></table>

Notes. In addition to the fitted parameters, we also show the feeding rate at the outer disc rim and the wind mass-loss rate. <sup>(a)</sup> $\dot{M}_{\text{wind}}$ denotes the mass-loss rate of one of the wind lobes. Since the wind is bipolar, conservation of mass implies $M_{\text{out}} = M_{\text{in}} + 2M_{\text{wind}}$.

absorption at $\phi_{\text{orb}} \approx 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_v^2 \le 6.85$. Their modelled absorption features are too broad in $\phi_{\text{orb}}$, and extend too much downwards into negative RV values. This is due to their lower fitted inclinations ($i < 45^{\circ}$ versus $i > 55^{\circ}$ 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 left-most 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_v^2 = 5.42$. This model's wind absorption feature reproduces the extent in $\phi_{\text{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 $\dot{M}_{\text{in}} \le 10^{-3} \, \text{M}_{\odot}/\text{yr}$. The fit assumes a low inclination of $i = 24^{\circ}$. 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 $\phi_{\text{orb}}$ peak at $\phi_{\text{orb}} \approx 0.5$ is slightly more narrow in $\phi_{\text{orb}}$ coverage and is concentrated at more negative RV values than the other peak at $\phi_{\text{orb}} \approx 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 $\phi_{\text{orb}} \approx 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_v^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 $\phi_{\text{orb}}$, and as a result fit significantly poorer at $\chi_v^2 \le 6.44$. The #MDW 2 fit is also poorer at $\chi_v^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^{\circ}$. Parametric modelling also implied a low inclination at $i_{\text{par}} = 25.7^{\circ}$ (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_v^2 = 3.30$. This model matches the observations very well, reproducing both the general RV extent, $\phi_{\text{orb}}$ width, absorption depth and weak redshifted absorption indent at $\phi_{\text{orb}} \approx 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_v^2 \le 3.96$, but #MDW 3 and 5 fit almost as well as #MDW 4 at $\chi_v^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^{\circ}$. #MDW 5 fits an especially low inclination at $i = 28^{\circ}$, while parametric modelling fitted a high inclination at $i_{\text{par}} = 71.8^{\circ}$. 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_v^2 = 4.91$. The model reproduces the observation very well, matching the RV extent and $\phi_{\text{orb}}$ width of the wind feature. The absorption depth is only slightly too shallow. The wind absorption from $\phi_{\text{orb}} \approx 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_v^2 \le 6.29$. The #MDW 5 model is extremely similar to the #MDW 4 one, reaching $\chi_v^2 = 5.05$, though with a more shallow absorption depth due to its lower inclination of $i = 17^{\circ}$ ($i \ge 29^{\circ}$ for all other MWs). 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_2 = 4.3 \, \text{M}_{\odot}$. In contrast, known SB1 orbital parameters imply $M_2$ 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^\circ$ and $37^\circ$, 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 $\phi_{\text{orb}} \approx 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 $\xi$ #MDW 4 and 5 clearly suppress the rotation-induced redshift more compared to the lower $\xi$ #MDW 2 and 3, and the lower $\xi$ and $\epsilon$ #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 $\phi_{\text{orb}}$, and reach RV values that are too negative (see Fig. C.1). These models fit lower inclinations of $i = 44^\circ$ and $28^\circ$, respectively, compared to $i > 55^\circ$ for the other models. Since #MDW 3 already fitted the absorption depth seen in the observations well at $i = 56^\circ$, 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 $\epsilon > 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 $\dot{M}_{\text{out}} \approx 10^{-4} \text{ M}_{\odot}/\text{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_{\text{cb}} = 10^{-2} \text{ M}_{\odot}$ (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 $\tau_{\text{cb}} = M_{\text{cb}}/\dot{M}_{\text{out}} \sim 100 \text{ yr}$. The post-AGB phase lifetime for a lower-mass star can however easily be $\tau_{\text{pAGB}} \gtrsim 10\,000 \text{ 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 $\approx 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 $\tau_{\text{cb}} \sim \tau_{\text{pAGB}}$. To push down the required feeding rate and increase the $\tau_{\text{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 ($\xi = 0.04$) fits in the range $\dot{M}_{\text{out}} \sim 10^{-4}-10^{-3} \text{ M}_{\odot}/\text{yr}$ (see Table C.1). #MDW 2 and 3 ($\xi \approx 0.10$) manage to push the lower end down by an order of magnitude to $\dot{M}_{\text{out}} \sim 10^{-5}-10^{-3} \text{ M}_{\odot}/\text{yr}$. #MDW 4 ($\xi = 0.215$) fits $\dot{M}_{\text{out}} \sim 10^{-5}-10^{-4} \text{ M}_{\odot}/\text{yr}$, while #MDW 5 ($\xi = 0.305$) fits $\dot{M}_{\text{out}} \sim 10^{-6}-10^{-4} \text{ M}_{\odot}/\text{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 $\tau_{\text{cb}} \sim 100 \text{ 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 $\epsilon = 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$_{\odot}$ (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 $\xi > 0.5$. Such solutions will need to be explored in future work to push down the fitted values for $\dot{M}_{\text{in}}$ and $\dot{M}_{\text{out}}$, and solve the current circumbinary disc lifetime discrepancy. Lower $\dot{M}_{\text{in}}$ values will also lead to slightly colder accretion disc predictions (see
Sect. 6.5) and lower secondary field strengths required to trun-
cate the accretion discs through magnetospheric funnelling (see
Sect. 6.3). It will also allow for thicker, more slowly rotating disc
solutions at $\epsilon > 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 $\ge 20 R_2$ 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 dif-
ferent 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 possi-
ble 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_{\text{in}} = 1 R_2$, 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 $\xi = 0.215$. Contrary to the solutions with lower $\xi$, the inner jet regions do not need to be cut out to suppress absorption at high negative RVs (see Table C.1). BD+46$^{\circ}$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_{\text{in}} = 5 R_2$. 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_{\text{in}} = 1 R_2$. 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):

$$
B_2 \sim 11 \, \text{kG} \cdot \left(\frac{r_{\text{in}}}{R_2}\right)^{\frac{7}{4}} \left(\frac{\dot{M}_{\text{in}}}{10^{-4} \, \text{M}_{\odot}/\text{yr}}\right)^{\frac{1}{2}} \left(\frac{M_2}{1 \, \text{M}_{\odot}}\right)^{\frac{1}{4}} \left(\frac{R_2}{1 \, \text{R}_{\odot}}\right)^{-\frac{5}{4}}. \quad (1)
$$

To start truncating the disc ($r_{\text{in}} > 1 R_2$), the reasonably well-
fitting models for the four targets discussed above imply mag-
netic field strengths greater than several kG and up to ~15 kG.
For comparison, a typical heavily accreting YSO at say $\dot{M}_{\text{in}} \approx$
$10^{-7} \text{ M}_{\odot}/\text{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 $\dot{M}_{\text{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_{\text{in}} \sim 1 R_2$. Constraints
on the secondaries' field strengths might provide more stringent
limits on $r_{\text{in}}$, but the brightness contrast with the primaries unfor-
tunately precludes direct determination of $B_2$ from spectra, using
for instance Zeeman broadening.

While the use of more efficient MDW solutions does lead to much smaller fitted $r_{\text{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 $\approx -120$ km/s, despite its low inclination (compare to TW Cam, where it reaches RV $\approx -320$ km/s). Correspondingly, HP Lyr is fitted with $r_{\text{in}} \ge 20 R_2$ for all MDW solutions (see Table C.1). HP Lyr's best fit would imply a field strength of $B_2 \approx 1.6$ MG. The \#MDW 5 fit (Table C.1) predicts the lowest field strength value for this target, at $B_2 \approx 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 $\sim 10^{-8} \text{ M}_{\odot}/\text{yr}$, the required field strength would only be brought down to $\sim 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_{\text{in}}$ value for HP Lyr's mod-
els 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_{\text{PAGB}}$), only
the hydrogen in the slow outer regions would be hot enough
to absorb significantly in $H\alpha$. 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^{-7} \text{ M}_{\odot}/\text{yr}$; Van Winckel 2003), if the lat-
ter 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_{\text{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, occur-
ring close to superior conjunction. As a result, the $H\alpha$ emis-
sion 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 construc-
tion. 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
<figure><img src="image_4.png" /><figcaption>Fig. 3. Example effective disc temperature profiles and spectra. Top: HP Lyr's #MDW 3 model. Bottom: TW Cam's #MDW 2 model. Left: $T_{\text{eff}}(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_{\text{full}}^{\text{BB}}$) and H band only ($T_{\text{H}}^{\text{BB}}$) are shown.</figcaption></figure>

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, $\dot{M}_{\text{cav}}$, needed to supply the missing absorption (assuming a homogeneous cylindrical flow at $T = T_{\text{wind}}$). If $\dot{M}_{\text{cav}} \ll \dot{M}_{\text{in}}$, we know that the extra outflow can be fed by material passing the inner disc rim in our models. If $\dot{M}_{\text{cav}} > \dot{M}_{\text{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 $\dot{M}_{\text{cav}}/\dot{M}_{\text{in}} \sim 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):

$$ T_{\text{eff}}(r) = \left( \eta_{\text{rad}} \cdot \frac{GM_2 \dot{M}_{\text{acc}}(r)}{8\pi\sigma r^3} \right)^{\frac{1}{4}} \\ \approx 36.4 \, \text{kK} \cdot \left( \frac{\eta_{\text{rad}} \dot{M}_{\text{acc}}(r)}{10^{-4} \, \text{M}_{\odot}/\text{yr}} \right)^{\frac{1}{4}} \left( \frac{M_2}{1 \, \text{M}_{\odot}} \right)^{\frac{1}{4}} \left( \frac{r}{1 \, \text{R}_{\odot}} \right)^{-\frac{3}{4}}, \quad (2) $$

with $\sigma$ 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_{\text{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_2$. As a result, despite the high accretion rate of $10^{-3} \, \text{M}_{\odot}/\text{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^{-3} \, \text{M}_{\odot}/\text{yr}$. Correspondingly, it reaches very high temperatures in the inner disc regions of up to ≈50 000 K.

### 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/PIIONIER), 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 black-body temperature $T_{\text{disc}} = 4000 \pm 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:

$$L_{\lambda} \propto \int_{r_{\text{in}}}^{r_{\text{out}}} r \cdot B_{\lambda}(T_{\text{eff}}(r)) dr, \quad (3)$$

with $B_{\lambda}$ 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_{\text{full}}^{\text{BB}}$. We note that this fits poorly to the multi-blackbody model spectra, which generally have a flattened peak and broader wings. $T_{\text{full}}^{\text{BB}}$ just serves as a rough but useful mean spectral temperature for the full disc spectrum. For the HP Lyr #MDW 3 model $T_{\text{full}}^{\text{BB}} = 3600$ K, while for the TW Cam #MDW 2 model $T_{\text{full}}^{\text{BB}} = 34300$ 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_{\text{disc}} = 4000 \pm 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_{\text{H}}^{\text{BB}}$. The resulting fits are shown in the right column of Fig. 3. For HP Lyr #MDW 3, we derive $T_{\text{H}}^{\text{BB}} = 3000$ K, agreeing well with the $T_{\text{full}}^{\text{BB}}$ value. For TW Cam's #MDW 2 model, however, we fit $T_{\text{H}}^{\text{BB}} = 4300$ K. This is much lower than $T_{\text{full}}^{\text{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 ($\ge 20 R_2$), 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_{\text{full}}^{\text{BB}}$ to $T_{\text{H}}^{\text{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.

Unlike what was posited by Paper I, hot accretion disc models thus do not contradict current near-infrared interferometric observations. If $T_{\text{full}}^{\text{BB}} \ge 10000$ 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.

<figure><img src="image_5.png" /><figcaption>Fig. 4. Blackbody temperatures derived from the accretion disc's spectral slope, considering either the full spectrum ($T_{\text{full}}^{\text{BB}}$, red diamonds), or only the H band ($T_{\text{H}}^{\text{BB}}$, blue dots). The orange shaded area denotes the $\pm 1\sigma$ region of the $T_{\text{disc}}$ estimate for IRAS 08544–4431 by Hillen et al. (2016), around which the $T_{\text{H}}^{\text{BB}}$ values cluster.</figcaption></figure>

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_{\text{full}}^{\text{BB}} \approx 18500$ 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 ($\xi > 0.10$), depending on the target, mitigating or solving multiple issues faced by #MDW 1 ($\xi > 0.04$) in matching the observed $H\alpha$ profiles or the expected circumbinary disc lifetime. Such high
efficiencies are only possible, at least for cold disc winds, at low magnetizations of $\mu \sim 10^{-3}-10^{-2}$. In addition, there is a strong preference for thicker discs at $\epsilon = 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 $\epsilon \gtrsim 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 $\xi > 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 circumbin- ion accretion disc and based on self-similar MHD disc wind solutions, to time-resolved observations of the H$\alpha$ 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 sam- ple 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$\alpha$ 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 $\xi \approx 0.30$ for discs of aspect ratio $\epsilon = 0.1$, is insufficient. Even at an efficiency of $\xi = 0.305$, predicted accretion rates may be as high as $10^{-4} M_{\odot}/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 $\xi \gtrsim 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 $\xi > 0.10$ and $\epsilon = 0.10$, which are only achieved at a low magnetization level of $\mu \sim 10^{-3}-10^{-2}$. 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 $\xi \gtrsim 0.5$ and for even thicker disc solutions at aspect ratios as high or even higher as $\epsilon = 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.

**Acknowledgements.** 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

* Alcock, C., Allsman, R. A., Alves, D. R., et al. 1998, AJ, 115, 1921

* Andrych, K., Kamath, D., Kluska, J., et al. 2023, MNRAS, 524, 4168

* Anugu, N., Kluska, J., Gardner, T., et al. 2023, ApJ, 950, 149

* Bertolami, M. M. M. 2016, A&A, 588, A25

* Bessolaz, N., Zanni, C., Ferreira, J., Keppens, R., & Bouvier, J. 2008, A&A, 478, 155

* Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883
Bollen, D., Winckel, H. V., & Kamath, D. 2017, A&A, 607, A60

Bollen, D., Kamath, D., Van Winckel, H., & De Marco, O. 2019, A&A, 631, A53

Bollen, D., Kamath, D., De Marco, O., Van Winckel, H., & Wardle, M. 2020, A&A, 641, A175

Bollen, D., Kamath, D., Van Winckel, H., De Marco, O., & Wardle, M. 2021, MNRAS, 502, 445

Bollen, D., Kamath, D., Van Winckel, H., et al. 2022, A&A, 666, A40

Bouvier, J., Chelli, A., Allain, S., et al. 1999, A&A, 349, 619

Bouvier, J., Grankin, K. N., Alencar, S. H. P., et al. 2003, A&A, 409, 169

Bujarrabal, V., Alcolea, J., Van Winckel, H., Santander-García, M., & Castro-Carrizo, A. 2013, A&A, 557, A104

Bujarrabal, V., Castro-Carrizo, A., Alcolea, J., & Van Winckel, H. 2015, A&A, 575, L7

Bujarrabal, V., Castro-Carrizo, A., Alcolea, J., et al. 2016, A&A, 593, A92

Bujarrabal, V., Castro-Carrizo, A., Alcolea, J., et al. 2017, A&A, 597, L5

Bujarrabal, V., Castro-Carrizo, A., Van Winckel, H., et al. 2018, A&A, 614, A58

Bujarrabal, V., Alcolea, J., Castro-Carrizo, A., et al. 2023, A&A, 677, L18

Casse, F., & Ferreira, J. 2000a, A&A, 353, 1115

Casse, F., & Ferreira, J. 2000b, A&A, 361, 1178

Clairmont, R., Steffen, W., & Koning, N. 2022, MNRAS, 516, 2711

Coelho, P. R. T. 2014, MNRAS, 440, 1027

Combet, C., & Ferreira, J. 2008, A&A, 479, 481

Corporaal, A., Kluska, J., Van Winckel, H., Kamath, D., & Min, M. 2023, A&A, 671, A15

Decin, L. 2021, ARA&A, 59, 337

Demircan, O., & Kahraman, G. 1991, Ap&SS, 181, 313

Ertel, S., Kamath, D., Hillen, M., et al. 2019, AJ, 157, 110

Ferreira, J. 1997, A&A, 319, 340

Ferreira, J., & Pelletier, G. 1995, A&A, 295, 807

Ferreira, J., Dougados, C., & Cabrit, S. 2006, A&A, 453, 785

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Frank, J., King, A., & Raine, D. 2002, Accretion Power in Astrophysics (Cambridge: Cambridge University Press)

Fromang, S., & Stone, J. M. 2009, A&A, 507, 19

Gaia Collaboration (Vallenari, A., et al.) 2023, A&A, 674, A1

Gallardo Cava, I., Gómez-Garrido, M., Bujarrabal, V., et al. 2021, A&A, 648, A93

Gielen, C., Cami, J., Bouwman, J., Peeters, E., & Min, M. 2011, A&A, 536, A54

Giridhar, S., Lambert, D. L., & Gonzalez, G. 2000, ApJ, 531, 521

Giridhar, S., Lambert, D. L., Reddy, B. E., Gonzalez, G., & Yong, D. 2005, ApJ, 627, 432

Gorlova, N., Van Winckel, H., Gielen, C., et al. 2012, A&A, 542, A27

Gorlova, N., Van Winckel, H., Ikonnikova, N. P., et al. 2015, MNRAS, 451, 2462

Habing, H. J., & Olofsson, H. 2004, Asymptotic Giant Branch Stars (New York: Springer)

Hillen, M., Kluska, J., Le Bouquin, J. B., et al. 2016, A&A, 588, L1

Hoffman, M. D., & Gelman, A. 2011, ArXiv e-prints [arXiv:1111.4246]

Höfner, S., & Olofsson, H. 2018, A&ARv, 26, 1

Jacquemin-Ide, J., Ferreira, J., & Lesur, G. 2019, MNRAS, 490, 3112

Johns-Krull, C. M., Valenti, J. A., & Koresko, C. 1999, ApJ, 516, 900

Kamath, D., Wood, P. R., & Van Winckel, H. 2014, MNRAS, 439, 2211

Kamath, D., Wood, P. R., & Van Winckel, H. 2015, MNRAS, 454, 1468

Kamath, D., Wood, P. R., Van Winckel, H., & Nie, J. D. 2016, A&A, 586, L5

Kluska, J., Van Winckel, H., Coppée, Q., et al. 2022, A&A, 658, A36

Launhardt, R., Pavlyuchenkov, Ya. N., Akimkin, V. V., et al. A&A, 678, A135

Lesur, G., & Longaretti, P.-Y. 2009, A&A, 504, 309

Lopez, J. A., Meaburn, J., & Palmer, J. W. 1993, ApJ, 415, L135

Lora, V., Toalá, J. A., González-Carbajal, J. I., et al. 2023, ArXiv e-prints [arXiv:2304.08668]

López-Vázquez, J. A., Lee, C. F., Fernández-López, M., et al. 2023, ArXiv e-prints [arXiv:2312.03272]

Manick, R., Van Winckel, H., Kamath, D., Hillen, M., & Escorza, A. 2017, A&A, 597, A129

Matt, S., & Pudritz, R. E. 2005, MNRAS, 356, 167

Moscadelli, L., Sanna, A., Beuther, H., Oliva, A., & Kuiper, R. 2022, Nat. Astron., 6, 1068

Oomen, G.-M., Van Winckel, H., Pols, O., et al. 2018, A&A, 620, A85

Oomen, G.-M., Van Winckel, H., Pols, O., & Nelemans, G. 2019, ApJ, 629, A49

Panoglou, D., Cabrit, S., Pineau des Forêts, G., et al. 2012, A&A, 538, A2

Preston, G. W., Krzeminski, W., Smak, J., & Williams, J. A. 1963, ApJ, 137, 401

Rao, S. S., & Giridhar, S. 2014, Rev. Mex. Astron. Astrofis., 50, 49

Raskin, G., Van Winckel, H., Hensberge, H., et al. 2011, A&A, 526, A69

Ray, T., & Ferreira, J. 2021, New Astron. Rev., 93, 101615

Salvesen, G., Simon, J. B., Armitage, P. J., & Begelman, M. C. 2016, MNRAS, 457, 857

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Tremblay, P.-E., Cummings, J., Kalirai, J. S., et al. 2016, MNRAS, 461, 2100

Van Winckel, H. 2003, ARA&A, 41, 391

Van Winckel, H. 2019, The Impact of Binary Stars on Stellar Evolution (Cambridge: Cambridge University Press)

Van Winckel, H., Waelkens, C., Fernie, J. D., & Waters, L. B. F. M. 1999, A&A, 343, 202

Van Winckel, H., Lloyd Evans, T., Briquet, M., et al. 2009, A&A, 505, 1221

Vassiliadis, E., & Wood, P. R. 1993, ApJ, 413, 641

Verhamme, O., Kluska, J., Ferreira, J., et al. 2024, A&A, 684, A79

Vidotto, A. A., Gregory, S. G., Jardine, M., et al. 2014, MNRAS, 441, 2361

Vilhu, O. 1984, A&A, 133, 117

Waelkens, C., Van Winckel, H., Bogaert, E., & Trams, N. R. 1991, A&A, 251, 495

Witt, A. N., Vijh, U. P., Hobbs, L. M., et al. 2009, ApJ, 693, 1946

Wurster, J., Bate, M. R., & Price, D. J. 2018, MNRAS, 481, 2450

Yvart, W., Cabrit, S., Pineau des Forêts, G., & Ferreira, J. 2016, A&A, 585, A74

Zhu, Z., & Stone, J. M. 2018, ApJ, 857, 34
**Appendix A: MDW example structures**

<caption>Fig. A.1 shows five example disc wind structures, calculated entirely analogously except using the different MDW solutions</caption>

as input. It showcases the general behaviour of the MDW solutions.

<figure><img src="image_6.png" /><figcaption>**Fig. A.1.** Model disc wind structures using different MDWs. Left: hydrogen number density *n*. Middle: poloidal velocity $v_p$. Right: toroidal/rotational velocity $v_\phi$. From top to bottom, we show #MDW 1, 2, 3, 4 and 5. Calculated assuming the HD 52961 system and *i* = 70°, $r_{in} = 5R_2$, $r_{out} = 0.99 R_{RL,2}$ and $\dot{M}_{in} = 10^{-3} \text{ M}_\odot/\text{yr}$. We note the slight discontinuity at a fixed polar angle in the *n* and $v_\phi$ 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.</figcaption></figure>
# 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, $\phi_{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 the 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_1 = 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_2$, 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 $\phi_{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, $\sigma_{back}(\lambda)$. 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 $\sigma_{back}(\lambda)$, which is added in quadrature to the S/N error to give the total spectral error $\sigma_{tot}(\lambda) = (\sigma_{S/N}^2 + \sigma_{back}^2(\lambda))^{1/2}$. $\sigma_{tot}$ is dominated by $\sigma_{back}$ for all our targets. If the $H_α$ emission component follows the primary's RV, the $\sigma_{tot}$ profile is shifted correspondingly. With this additional error, we typically get a reduced chi-squared $\chi_v^2 \sim 1$ when comparing the background spectrum to the observed unobscured spectra. This shows that the combination of the background spectrum and $\sigma_{tot}$ is a good statistical representation for the unobscured spectra. The total spectral error profiles for our targets are shown in Fig. B.1.
<table><caption>**Table B.1.** Spectroscopic orbital parameters.</caption><thead><tr><th>Object<br>Parameter</th><th>HD 52961</th><th>BD+46°442</th><th>TW Cam</th><th>IRAS 19135+3937</th><th>HP Lyr</th></tr></thead><tbody><tr><td><i>P</i> (d)</td><td>1288.6 ± 0.3</td><td>140.82 ± 0.02</td><td>662.2 ± 5.3</td><td>126.97 ± 0.08</td><td>1818 ± 80</td></tr><tr><td><i>e</i></td><td>0.23 ± 0.01</td><td>0.085 ± 0.005</td><td>0.25 ± 0.04</td><td>0.13 ± 0.03</td><td>0.20 ± 0.04</td></tr><tr><td><i>T</i><sub>0</sub> − 24 · 10<sup>5</sup> (d)</td><td>7308 ± 24</td><td>55233.9 ± 1.3</td><td>55111 ± 18</td><td>54997.7 ± 1.0</td><td>56175 ± 61</td></tr><tr><td><i>ω</i> (°)</td><td>297.4 ± 5.9</td><td>95.8 ± 3.3</td><td>144.4 ± 10</td><td>66.0 ± 4.4</td><td>14.2 ± 13</td></tr><tr><td><i>K</i><sub>1</sub> (km s<sup>−1</sup>)</td><td>13.1 ± 0.3</td><td>23.8 ± 0.1</td><td>14.1 ± 0.6</td><td>18.0 ± 0.6</td><td>7.8 ± 0.2</td></tr><tr><td><i>γ</i> (km s<sup>−1</sup>)</td><td>6.2 ± 0.2</td><td>−98.13 ± 0.08</td><td>−49.8 ± 0.5</td><td>2.1 ± 0.4</td><td>−115.6 ± 0.2</td></tr><tr><td><i>a</i><sub>1</sub> sin <i>i</i> (AU)</td><td>1.507 ± 0.034</td><td>0.3074 ± 0.0014</td><td>0.83 ± 0.04</td><td>0.209 ± 0.008</td><td>1.27 ± 0.06</td></tr><tr><td><i>f</i>(<i>m</i>) (M<sub>☉</sub>)</td><td>0.274 ± 0.019</td><td>0.195 ± 0.003</td><td>0.174 ± 0.022</td><td>0.075 ± 0.008</td><td>0.083 ± 0.007</td></tr></tbody></table>

**Notes.** Values taken from Oomen et al. (2018). Parameters include the orbital period *P*, eccentricity *e*, time of periastron passage *T*<sub>0</sub>, argument of the periastron *ω*, primary semi-amplitude *K*<sub>1</sub>, systemic velocity *γ*, projected primary semi-major axis *a*<sub>1</sub> sin *i* and the mass function *f*(*m*).

<figure><img src="image_7.png" /><figcaption>**Fig. B.1.** Total *H*<sub>α</sub> spectral error profiles, accounting for the variability of the observed unobscured spectra.</figcaption></figure>
## 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^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^2$ values, but ultimately depends on the ability of the MDW to recover the global $H_\alpha$ pattern.

<table><caption>**Table C.1.** Parameters and $\chi^2$ values for all fitted disc wind models.</caption><thead><tr><th>Object</th><th>#MDW</th><th>i<br>(°)</th><th>r<sub>in</sub><br>(R<sub>2</sub>)</th><th>r<sub>out</sub><br>(R<sub>RL,2</sub>)</th><th>M&#775;<sub>in</sub><br>(M<sub>☉</sub>/yr)</th><th>T<sub>wind</sub><br>(K)</th><th>χ<sup>2</sup><sub>v</sub></th><th>M&#775;<sub>out</sub><br>(M<sub>☉</sub>/yr)</th><th>M&#775;<sub>wind</sub><br>(M<sub>☉</sub>/yr)</th><th>R<sub>RL,2</sub><br>(R<sub>2</sub>)</th><th>R<sub>2</sub><br>(R<sub>☉</sub>)</th></tr></thead><tbody><tr><td rowspan="5">HD 52961</td><td>1</td><td>80</td><td>25</td><td>0.9</td><td>10<sup>-4</sup></td><td>4200</td><td>9.36</td><td>1.1 &middot; 10<sup>-4</sup></td><td>5 &middot; 10<sup>-6</sup></td><td>232</td><td>0.89</td></tr><tr><td>2</td><td>76</td><td>1</td><td>0.9</td><td>10<sup>-4</sup></td><td>4200</td><td>6.92</td><td>1.8 &middot; 10<sup>-4</sup></td><td>4 &middot; 10<sup>-5</sup></td><td>231</td><td>0.91</td></tr><tr><td>3</td><td>72</td><td>5</td><td>0.9</td><td>10<sup>-4</sup></td><td>4200</td><td>7.40</td><td>1.4 &middot; 10<sup>-4</sup></td><td>2 &middot; 10<sup>-5</sup></td><td>231</td><td>0.91</td></tr><tr><td><b>4</b></td><td><b>72</b></td><td><b>1</b></td><td><b>0.9</b></td><td><b>10<sup>-4</sup></b></td><td><b>4200</b></td><td><b>5.32</b></td><td><b>3 &middot; 10<sup>-4</sup></b></td><td><b>1 &middot; 10<sup>-4</sup></b></td><td><b>228</b></td><td><b>0.93</b></td></tr><tr><td>5</td><td>56</td><td>5</td><td>0.9</td><td>10<sup>-5</sup></td><td>4400</td><td>5.63</td><td>3 &middot; 10<sup>-5</sup></td><td>1 &middot; 10<sup>-5</sup></td><td>212</td><td>1.10</td></tr><tr><td rowspan="5">BD+46°442</td><td>1</td><td>64</td><td>15</td><td>0.9</td><td>10<sup>-3</sup></td><td>5650</td><td>6.73</td><td>1.06 &middot; 10<sup>-3</sup></td><td>3 &middot; 10<sup>-5</sup></td><td>54</td><td>0.87</td></tr><tr><td>2</td><td>68</td><td>1</td><td>0.9</td><td>10<sup>-5</sup></td><td>5050</td><td>7.69</td><td>1.6 &middot; 10<sup>-5</sup></td><td>3 &middot; 10<sup>-6</sup></td><td>55</td><td>0.84</td></tr><tr><td><b>3</b></td><td><b>56</b></td><td><b>5</b></td><td><b>0.8</b></td><td><b>10<sup>-5</sup></b></td><td><b>5250</b></td><td><b>6.44</b></td><td><b>1.2 &middot; 10<sup>-5</sup></b></td><td><b>1 &middot; 10<sup>-6</sup></b></td><td><b>46</b></td><td><b>0.96</b></td></tr><tr><td>4</td><td>44</td><td>1</td><td>0.5</td><td>10<sup>-5</sup></td><td>5050</td><td>6.85</td><td>2.0 &middot; 10<sup>-5</sup></td><td>5 &middot; 10<sup>-6</sup></td><td>26</td><td>1.20</td></tr><tr><td>5</td><td>28</td><td>1</td><td>0.4</td><td>10<sup>-6</sup></td><td>5450</td><td>7.18</td><td>2.4 &middot; 10<sup>-6</sup></td><td>7 &middot; 10<sup>-7</sup></td><td>17</td><td>2.12</td></tr><tr><td rowspan="5">TW Cam</td><td>1</td><td>28</td><td>15</td><td>0.9</td><td>10<sup>-3</sup></td><td>4650</td><td>5.59</td><td>1.08 &middot; 10<sup>-3</sup></td><td>4 &middot; 10<sup>-5</sup></td><td>107</td><td>2.00</td></tr><tr><td>2</td><td>44</td><td>1</td><td>0.6</td><td>10<sup>-3</sup></td><td>4650</td><td>6.57</td><td>1.6 &middot; 10<sup>-3</sup></td><td>3 &middot; 10<sup>-4</sup></td><td>89</td><td>1.14</td></tr><tr><td><b>3</b></td><td><b>24</b></td><td><b>5</b></td><td><b>0.9</b></td><td><b>10<sup>-3</sup></b></td><td><b>4850</b></td><td><b>5.42</b></td><td><b>1.4 &middot; 10<sup>-3</sup></b></td><td><b>2 &middot; 10<sup>-4</sup></b></td><td><b>98</b></td><td><b>2.52</b></td></tr><tr><td>4</td><td>24</td><td>1</td><td>0.9</td><td>10<sup>-4</sup></td><td>4850</td><td>6.44</td><td>2.6 &middot; 10<sup>-4</sup></td><td>8 &middot; 10<sup>-5</sup></td><td>98</td><td>2.52</td></tr><tr><td>5</td><td>20</td><td>1</td><td>0.6</td><td>10<sup>-4</sup></td><td>4450</td><td>6.68</td><td>3 &middot; 10<sup>-4</sup></td><td>1 &middot; 10<sup>-4</sup></td><td>58</td><td>3.41</td></tr><tr><td rowspan="5">IRAS 19135+3937</td><td>1</td><td>56</td><td>10</td><td>0.9</td><td>10<sup>-3</sup></td><td>5850</td><td>4.28</td><td>1.08 &middot; 10<sup>-3</sup></td><td>4 &middot; 10<sup>-5</sup></td><td>57</td><td>0.67</td></tr><tr><td>2</td><td>68</td><td>5</td><td>0.9</td><td>10<sup>-3</sup></td><td>5650</td><td>3.96</td><td>1.2 &middot; 10<sup>-3</sup></td><td>1 &middot; 10<sup>-4</sup></td><td>60</td><td>0.59</td></tr><tr><td>3</td><td>48</td><td>5</td><td>0.9</td><td>10<sup>-5</sup></td><td>5450</td><td>3.49</td><td>1.2 &middot; 10<sup>-5</sup></td><td>1 &middot; 10<sup>-6</sup></td><td>54</td><td>0.75</td></tr><tr><td><b>4</b></td><td><b>48</b></td><td><b>1</b></td><td><b>0.9</b></td><td><b>10<sup>-5</sup></b></td><td><b>5050</b></td><td><b>3.30</b></td><td><b>2.4 &middot; 10<sup>-5</sup></b></td><td><b>7 &middot; 10<sup>-6</sup></b></td><td><b>54</b></td><td><b>0.75</b></td></tr><tr><td>5</td><td>28</td><td>1</td><td>0.9</td><td>10<sup>-5</sup></td><td>4650</td><td>3.53</td><td>3 &middot; 10<sup>-5</sup></td><td>1 &middot; 10<sup>-5</sup></td><td>42</td><td>1.32</td></tr><tr><td rowspan="5">HP Lyr<sup>a</sup></td><td>1</td><td>37</td><td>25</td><td>0.8</td><td>10<sup>-4</sup></td><td>4800</td><td>7.17</td><td>1.1 &middot; 10<sup>-4</sup></td><td>5 &middot; 10<sup>-6</sup></td><td>247</td><td>1.00</td></tr><tr><td>2</td><td>61</td><td>25</td><td>0.8</td><td>10<sup>-4</sup></td><td>5000</td><td>8.19</td><td>1.2 &middot; 10<sup>-4</sup></td><td>1 &middot; 10<sup>-5</sup></td><td>298</td><td>0.65</td></tr><tr><td>3</td><td>33</td><td>20</td><td>0.3</td><td>10<sup>-3</sup></td><td>4400</td><td>6.29</td><td>1.16 &middot; 10<sup>-3</sup></td><td>8 &middot; 10<sup>-5</sup></td><td>87</td><td>1.15</td></tr><tr><td><b>4</b></td><td><b>29</b></td><td><b>20</b></td><td><b>0.3</b></td><td><b>10<sup>-4</sup></b></td><td><b>5000</b></td><td><b>4.91</b></td><td><b>1.4 &middot; 10<sup>-4</sup></b></td><td><b>2 &middot; 10<sup>-5</sup></b></td><td><b>82</b></td><td><b>1.32</b></td></tr><tr><td>5</td><td>17</td><td>20</td><td>0.3</td><td>10<sup>-4</sup></td><td>5200</td><td>5.05</td><td>1.4 &middot; 10<sup>-4</sup></td><td>2 &middot; 10<sup>-5</sup></td><td>61</td><td>2.91</td></tr></tbody></table>

**Notes.** The best fitting model per target is marked in bold. In addition to the fitted parameters, we show the feeding rate at the outer disc rim and the wind mas-loss rate, and also provide $R_{RL,2}$ in units of $R_2$ and $R_2$ in solar radii so the extent of the disc can be calculated for future reference. <sup>(a)</sup>The values for HP Lyr #MDW 2 correspond to a secondary, high inclination minimum in the parameter space. The first minimum represented an unphysical solution, implying an unrealistically large secondary mass of $M_2 \approx 23 M_\odot$ and, unlike what is observed, a wind absorption feature spanning the entire orbit.
<figure><img src="image_8.png" /><figcaption>**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.</figcaption></figure>
## 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 $\lambda$ leaving the disc wind is then given by:

$$I_{\lambda}^{\text{out}} = I_{\lambda}^{0} \cdot e^{-\Delta\tau_{\lambda,\text{out}}} = I_{\lambda}^{0} \cdot e^{-(\Delta\tau_{\lambda,\text{left}}+\Delta\tau_{\lambda,\text{cav}}+\Delta\tau_{\lambda,\text{right}})}, \quad (\text{D.1})$$

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}^{\text{out}}$ with the observed spectrum, $I_{\lambda,\text{obs}}$, and $I_{\lambda}^{0} \cdot e^{-(\Delta\tau_{\lambda,\text{left}}+\Delta\tau_{\lambda,\text{right}})}$ as the spectrum modelled by our MDW fits, $I_{\lambda,\text{mod}}$, this gives:

$$I_{\lambda,\text{obs}} = I_{\lambda,\text{mod}} \cdot e^{-\Delta\tau_{\lambda,\text{cav}}} \Rightarrow \Delta\tau_{\lambda,\text{cav}} = \ln\left(\frac{1}{1 - \frac{\Delta\lambda}{I_{\lambda,\text{mod}}}\cdot\frac{1}{\Delta\tau_{\lambda,\text{cav}}}}\right), \quad (\text{D.2})$$

where $\Delta I_{\lambda} = I_{\lambda,\text{mod}} - I_{\lambda,\text{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).

<figure><img src="image_9.png" /><figcaption>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) \cdot \sin i$.</figcaption></figure>

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

<figure><img src="image_10.png" /><figcaption>Fig. D.2. Example spectrum at $\phi_{\text{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 $\Delta I_{\lambda}/I_{\lambda,\text{mod}}$ is shown as a dashed blue line. The central wavelength value at which absorption is missing (RV $\approx -60$ km/s) is indicated by the black dot-dashed line.</figcaption></figure>

$\Delta\tau_{\lambda,\text{cav}} = \kappa_{\lambda,\text{cav}} \cdot \rho_{\text{cav}} \cdot s$, with $\kappa_{\lambda,\text{cav}}$ the material's $H_{\alpha}$ line opacity profile and $\rho_{\text{cav}}$ its mass density, assuming homogeneity. For $\kappa_{\lambda,\text{cav}}$, we choose the $H_{\alpha}$ 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_{\text{wind}}$), since the $H_{\alpha}$ opacity is strongly dependent on temperature, but only weakly pressure-dependent at the considered densities. Choosing a range of wavelengths $\Delta\lambda$, spanning the background emission peak that shines through the model, this can then be averaged over $\Delta\lambda$:

$$\overline{\Delta\tau_{\lambda,\text{cav}}} = \overline{\kappa_{\lambda,\text{cav}} \cdot \rho_{\text{cav}} \cdot s} \Rightarrow \rho_{\text{cav}} = \frac{\overline{\Delta\tau_{\lambda,\text{cav}}}}{\overline{\kappa_{\lambda,\text{cav}} \cdot s}}, \quad (\text{D.3})$$

$$\overline{\kappa_{\lambda,\text{cav}}} = \frac{\int_{\Delta\lambda} \kappa_{\lambda,\text{cav}} d\lambda}{\Delta\lambda}, \quad (\text{D.4})$$

$$\overline{\Delta\tau_{\lambda,\text{cav}}} = \frac{\int_{\Delta\lambda} \ln\left(\frac{1}{1-\Delta I_{\lambda}/I_{\lambda,\text{mod}}}\right) d\lambda}{\Delta\lambda}. \quad (\text{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) \cdot \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 $\rho_{\text{cav}}$, we retrieve an estimate of the extra component's mass-loss rate:

$$\dot{M}_{\text{cav}} \sim \pi r^2 \cdot v \cdot \rho_{\text{cav}} = \frac{\pi \sin^2 i}{4 \cos i} \cdot \frac{RV \cdot s}{\overline{\kappa_{\text{cav}}} \cdot \overline{\Delta\tau_{\text{cav}}}}. \quad (\text{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_{\text{cav}}} = 10-100 \text{ cm}^2/\text{g}$, $\rho_{\text{cav}} \approx 10^{-15} \text{ g/cm}^3$, $v \approx 100 \text{ km/s}$ and $\dot{M}_{\text{cav}} \approx 10^{-8} \text{ M}_\odot/\text{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.