© ESO 2019

1. Introduction

Astrophysical jets are frequently observed phenomena in the Universe, ranging from extremely energetic jets in active galactic nuclei (AGNs; Blandford et al. 2019) to non-relativistic stellar jets emerging from young stellar objects (YSOs; Frank et al. 2014) and planetary nebulae (PNe; Livio 1999). These jets are an important source of feedback to the interstellar medium and can influence the evolution of the launching engine and of the ambient medium (Soker 2016, and references therein).

Although a wide variety of jets are launched by distinct astrophysical objects, their general structure is very similar. In general, jets originate from a compact central object, they are two-sided, and have a certain degree of collimation, caused most likely by a magnetic field (de Gouveia Dal Pino 2005). Observations of YSOs clearly show two-sided or bipolar jets (Reipurth et al. 1997, 2000, 2002; Bally et al. 2002; Heathcote et al. 1996; Burrows et al. 1996). Observations of molecular outflows in YSOs also show two jet types: a wide-angle outflow and a highly collimated jet (Lee et al. 2000). The narrow, high-velocity jet propagates through the ambient medium and ends in a bow shock. These narrow jets are often accompanied by a slower, wide-angled outflow component which also carves a cavity in the ambient medium (Bally 2016; Melnikov et al. 2018).

Collimated outflows also exist in evolved objects such as PNe and have been extensively modelled (e.g., Greenhill et al. 1998, Orion Source I; Bacciotti et al. 2000, DG Tauri; Coffey et al. 2004, RW Aur, TH 28, and LkHα 321; Sánchez Contreras & Sahai 2001, He 3-1475). In this contribution, we focus on low- and intermediate-mass, binary post-asymptotic giant branch (AGB) stars recently discovered to have jets (Gorlova et al. 2012, 2015; Bollen et al. 2017). Our aim is to characterise these systems and determine the jet launching mechanism.

Binary interactions during the AGB phase have a significant impact on the evolution of the binary system. When the binary orbit is small and the expanding AGB star grows larger than its Roche lobe, mass transfer will ensue and the orbital elements will evolve. If unstable mass transfer takes place, the binary will go through a common envelope phase with substantial orbital shrinkage (Ivanova et al. 2013), leaving behind a close binary system (De Marco 2009). If the mass transfer is more stable, or if the initial orbital separation is large enough that Roche lobe overflow is not achieved, then the orbital separation might not decrease as much or even grow. A third possibility is an intermediate configuration with substantial mass transfer that however avoids a classical common envelope or delays it (Reichardt et al. 2019) and may result in a post-AGB binary system with periods of 100–2000 days (Van Winckel et al. 2009)¹.

Binary post-AGB systems display a distinctive near-infrared (near-IR) excess in their spectral energy distribution (SED) that has been shown to be a signature of a stable, Keplerian circumbinary disk (Van Winckel 2003; Bujarrabal et al. 2013; Kamath et al. 2014, 2015). This implies a connection between the binary interaction and the formation of the circumbinary disk.

The post-AGB binary system therefore comprises a circumbinary disk, an evolved post-AGB star (the primary component), and a companion, likely to be a main sequence star (Oomen et al. 2018). Using interferometric image reconstruction techniques, signatures of an additional component have been detected in the post-AGB binary system IRAS 08544−4431 (Hillen et al. 2016; Kluska et al. 2018), that are interpreted as an accretion disk around the companion, possibly fed by the primary or by reaccretion of circumbinary material.

High resolution, radial velocity (RV) monitoring of several post-AGB stars by Van Winckel et al. (2009) has detected variable absorption features interpreted as jets emanating from the companion (Gorlova et al. 2012, 2014, 2015; Bollen et al. 2017). During superior conjunction, when the companion is closest to the observer and the evolved star is behind, continuum photons from the primary are scattered out of the line-of-sight by the H-atoms in the jet, resulting in a P-Cygni line profile during these orbital phases. Since the light from the primary shines through different parts of the jet during the orbital motion of the binary system, these time-series provide us with a tomography of the jet that may uniquely constrain the spatial, velocity, and density structure of the jet.

Bollen et al. (2017) fitted the time-series spectra of the post-AGB binary BD+46°442 with a simple conical jet model and determined a rather wide jet opening angle (> 100°), where a low density fast central outflow is enclosed in a denser and slower outer layer.

We have now observed that a large fraction of well-monitored post-AGB binaries show the presence of similar features, indicating that jets are commonly found in these systems. In this study, we present a refined technique to determine improved parameters such as outflow velocities and density scales, which are crucial to establish the jet outflow momenta and determine mass-accretion rates as a function of orbital parameters. We use the Markov chain Monte Carlo (MCMC) technique to fit the observed Hα-line profile variations. The strength of this work and its application of MCMC is that we now have an automated model-fitting code, which provides more information on the complete kinematic, geometrical, and density structure of these jets. The ultimate aim of our study is to apply our fitting routine to all post-AGB binary systems with jets, so as to investigate the observed diversity of the jets, understand the jet launching mechanism, and the dynamics of the interaction that leads to the formation of post-AGB binaries.

This paper is structured as follows: in Sect. 2, we present the jet model that we implemented in this code. In Sect. 3, we elaborate on our model-fitting routine. In Sect. 4, we illustrate our jet model applied to IRAS 19135+3937, a post-AGB binary system (Gorlova et al. 2015), which is presented in Sect. 4. We discuss the results of our fitting routine in Sect. 5. We finish with our conclusions in Sect. 6.

2. The post-AGB jet model

The central stellar binary system in our model consists of an evolved star (the post-AGB star) and a companion (the low-mass main sequence star). We model the post-AGB star as a uniform disk that is always perpendicular to the line-of-sight, because the evolved star will be projected as a disk along the line-of-sight towards the observer. We divide this disk up into a Fibonacci grid, whose benefit is that each grid point is a surface with equal area, as can be seen in Fig. 1. We follow the light of the evolved star commencing from these grid points, along the line-of-sight towards the observer, to calculate the absorption by the jet.

Fig. 1. Mesh grid of the primary star component and the rays going through the jet. The primary component is sampled as a Fibonacci grid with 200 grid points. The absorption of the continuum light by the H-atoms in the jet is determined for each ray along the line-of-sight that departs from these points on the primary’s surface.

The luminosity of the companion star in these systems is often negligible relative to the primary’s luminosity. Hence, its flux contribution is too small to be observed in the spectra. For this reason, we consider the total flux to be from the primary and ignore the flux contribution of the companion star in our model.

In post-AGB binary systems, jets are launched from the accretion disk around the companion. The jets travel in opposite directions, perpendicular to the orbital plane of the binary system. During each orbital phase, we determine which rays from the evolved star pass through the jet cones in our model. We then follow these rays through the jet. Each ray through the jet is sampled by a number of grid points between the point at which the ray enters the jet and the point at which it exits. At these grid points, we calculate the local density, velocity, and velocity dispersion, from which the optical depth due to line scattering can be calculated over the whole ray and for all frequencies.

In the following subsections, we describe the geometrical structures of the jet that we implement in our model. We will then describe how we recreate the absorption feature in our model that is caused by the jet occultation by giving a detailed description of the radiative transfer of the AGB star light from the evolved star through the jet.

2.1. Jet geometry

We aim to determine the structure of jets launched in post-AGB binary systems. We choose a jet geometry based on jet launching mechanisms originally developed for YSOs, because of the similarity between the observational data between young stars and our evolved systems. Over the past decades, several jet launching theories have been developed in order to explain YSO jets. The most widely accepted ones are those of Matt & Pudritz (2005, Model A), the magneto-centrifugally accelerated disk wind model by Blandford & Payne (1982, Model B), and the X-wind model by Shu et al. (1994, Model C).

These jet models are applied by Kurosawa et al. (2006, 2011), Weigelt et al. (2011), and Tambovtseva et al. (2014) to determine the structure and physical properties of jets in YSOs. In these studies the authors determined the density and velocity structure of jets in classical T Tauri stars (CTTSs) by generating synthetic hydrogen and helium spectra and fitting those to their observations. Since these studies focused on the modelling of single star systems, the hydrogen and helium lines remain unchanged (neglecting changes due to changing accretion rates). Our work has extra complexity, which lies in the binary nature of the systems. This results in spectral lines that vary over the orbital period. We account for this by implementing the binary motion in our model, such that we can calculate the absorption caused by the jet during each orbital phase of the binary system.

In order to find which model best represents our diverse sample of post-AGB binaries, we implement all three and carry out fits to our data. Below, we give a detailed description of the velocity and density structure for each of these jet configurations.

2.1.1. Model A: “stellar jet” configuration

The first configuration is similar to that applied by Kurosawa et al. (2006), who modelled the Hα emission in CTTSs. They based this on the jet model by Matt & Pudritz (2005). According to Matt & Pudritz (2005) the star, that is surrounded by an accretion disk, shows that the spin-down torque of the star is most likely due to the stellar wind, travelling along the open field lines of the star. The stellar magnetic field is strong enough to truncate the accretion disk, which causes disk material to be channelled along the magnetic field lines, part accreting onto the star and part being ejected. Thus, the ejected matter originates from the disk, but is ejected along the open stellar magnetic field lines that are anchored on the star and therefore we call this mode the “stellar jet” model. The magnetic field lines of the companion star in our systems are assumed to be similar to those of YSOs².

In our model, we assume that the jet travels along the open magnetic field lines of the companion star. Hence, the matter in the jet is ejected radially away from the star, as can be seen in panel A of Fig. 2. We implement a latitudinally dependent velocity law for the matter inside the jet, as done by Thomas et al. (2013). The first one is:

$$\begin{array}{c}\hfill \mathit{v}=[v_0+(v_{\alpha }-v_0){\left(\frac{\theta }{\alpha }\right)}^2]\widehat{\mathit{r}},\end{array}$$(1)

Fig. 2. Upper panel: schematic sketch of the jet model, showing our geometric implementation of the three models. Lower panel: launching region and velocity field of the jet for the stellar jet (A), the disk wind (B), and the X-wind (C). Unlike the stellar jet and X-wind, the launching region for the disk wind is the region between the inner-disk radius din and the outer disk radius dout. This configuration is based on the disk wind configuration applied by Kurosawa et al. (2006).

where θ is the polar latitudinal coordinate, α the half-opening angle of the jet, v_α is the velocity at the jet edges and v₀ is the velocity along the jet axis. Hence, the jet velocities are independent of the radial component.

Just like the jet velocity, the density at each position in the jet has a latitudinal dependency, where the jet has a low density along the jet axis and reaches the highest densities along the jet edges. The density decreases as a function of jet height z as follows:

$$\begin{array}{c}\hfill \rho (\theta ,z)\propto \frac{{\theta }^p}{z^2}·\end{array}$$(2)

By choosing different values for p, the density dispersion in the jet can be increased or decreased in our model. The range of adopted values for the factor p is given in Table 1.

2.1.2. Model B: “disk-wind” configuration with inner jet

The magneto-centrifugal disk wind model, first introduced by Blandford & Payne (1982), is suggested to be the main launching mechanism for jets from accretion disks around stellar objects, such as black holes and young stars (Pudritz & Norman 1983). In this model, the magnetic field lines are anchored to the accretion disk and inclined with respect to the rotation axis of the accretion disk by at least 30°. The matter in the disk below and above the mid-plane will be accelerated due to the magneto-centrifugal force and follow these magnetic field lines like beads on a wire.

In our model, we implement this disk wind as our second jet configuration. The launching region of the jet in this configuration is the circum-companion accretion disk. The inner and outer disk radii are determined from the inner and outer jet angles in the model, and have set limits: the inner disk radius must be larger than twice the radius of the companion star and the outer disk radius must be smaller than the Roche lobe. Since we do not know the exact size of the companion star, we assume a radius given by the stellar mass-radius relation by Demircan & Kahraman (1991):

$$\begin{array}{c}\hfill R_{\star }=1.01R_{\odot }{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{0.724}.\end{array}$$(3)

To determine the Roche radius of the companion, we use the approximation by Eggleton (1983):

$$\begin{array}{c}\hfill r_L=\frac{0.49q^{2/3}}{0.6q^{2/3}+ln(1+q^{1/3})},\end{array}$$(4)

where q = M₁/M₂ is the mass ratio of the binary system. The jet configuration has two outflowing components: the inner jet component and the disk wind component. The inner jet component has a similar velocity law as the stellar jet discussed in Sect. 2.1.1. The disk wind originates from the circum-companion disk, in which the particles are ejected latitudinally from the region between the inner and outer disk edges at angles between α_in and α_out, as shown in panel B of Fig. 2. The disk wind velocity is related to the Keplerian velocity of the disk. This means that the disk wind velocity of the particles ejected from a certain radius in the disk equals the Keplerian velocity at that position in the disk, multiplied with a scaling factor c_v. This scaling factor is a free parameter in our model and has a value between 0 and 1. The angle-dependent velocity law is defined as follows:

$$\begin{array}{c}\hfill v=\{\begin{array}{cc}{c}_{\mathrm{v}}·\sqrt{\frac{\mathrm{G}{M}_2}{d}}=v_{\alpha },out\sqrt{\frac{tan{\alpha }_{\mathrm{out}}}{tan\theta }};\hfill & \mathrm{if}\theta >{\alpha }_{\mathrm{in}}\hfill \\ [v_0+(v_{\alpha },in-v_0){\left(\frac{\theta }{{\alpha }_{\mathrm{in}}}\right)}^2];\hfill & \mathrm{if}\theta <{\alpha }_{\mathrm{in}},\hfill \end{array}\end{array}$$(5)

with c_v the scaling factor, d the radial distance between the launch site in the disk and the centre of the companion star, M₂ the mass of the companion, and $v_{\alpha ,\mathrm{in}}=c_{\mathrm{v}}·\sqrt{G{M}_2/d_{\mathrm{in}}}$ the Keplerian velocity at the inner disk rim d_in.

The density profile in the jet for this system is determined as follows:

$$\begin{array}{c}\hfill \rho (\theta ,z)\propto \{\begin{array}{cc}{\left(\frac{\theta }{{\alpha }_{\mathrm{in}}}\right)}^{p_{\mathrm{out}}}{z}^{-2};\hfill & \mathrm{if}\theta >{\alpha }_{\mathrm{in}}\hfill \\ \hfill \\ {\left(\frac{\theta }{{\alpha }_{\mathrm{in}}}\right)}^{p_{\mathrm{in}}}{z}^{-2};\hfill & \mathrm{if}\theta <{\alpha }_{\mathrm{in}}.\hfill \end{array}\end{array}$$(6)

2.1.3. Model C: “X-wind” configuration with inner jet

The X-wind theory for YSOs was introduced by Shu et al. (1994) for a scenario where the interaction between the magnetosphere of rapidly rotating YSOs and the surrounding accretion disks can give rise to the formation of so-called X-winds. In this model, the central star is surrounded by an accretion disk that extends inwards towards the star, down to an inner disk radius of about twice the stellar radius of the central star. The stellar magnetic field is not able to penetrate the accretion disk smoothly, however. This causes the magnetic field to be squeezed together in the inner disk region, named the X-region. A fraction of the inflowing matter in the disk will be launched from the equatorial plane by the magnetic field in the X-region, thus creating an X-shaped outflow. This wind efficiently carries away angular momentum from the disk. In other words, this outflow originates from the inner disk edge, where the disk’s magnetic field and the co-rotating stellar field meet. A more detailed description of the X-wind can be found in Shu et al. (1997). Shang et al. (2007) pointed out that the disk wind and the X-wind are intrinsically of the same nature: both are magneto-centrifugally accelerated winds, with the only difference being the truncation of the magnetic field in the accretion disk. Where the magnetic field of a disk wind can penetrate through the whole accretion disk, the field of an X-wind will only truncate the inner disk region.

In our model, the X-wind configuration comprises two main jet components, each with a distinct velocity profile. The true launching region of the X-wind is about twice the stellar radius of the main sequence star. As this is relatively small compared to the scale of the binary system, we consider the launching point in our model to be located at the position of the main sequence star, similar as the stellar jet model. The velocity structure in this model is defined as follows:

$$\begin{array}{c}\hfill \mathit{v}=\{\begin{array}{cc}[v_{\alpha ,\mathrm{out}}+(v_{\mathrm{M}}-v_{\alpha ,\mathrm{out}})·cos\left(\frac{\pi }{2}\frac{\theta }{{\alpha }_{\mathrm{out}}}\right)]\widehat{r};\hfill & \mathrm{if}\theta >{\alpha }_{\mathrm{in}}\hfill \\ [v_{\alpha ,\mathrm{out}}+(v_{\mathrm{M}}-v_{\alpha ,\mathrm{out}})·cos\left(\frac{\pi }{2}\frac{\theta }{{\alpha }_{\mathrm{out}}}\right)\hfill \\ +(v_0-v_{\mathrm{M}})·cos\left(\frac{\pi }{2}\frac{\theta }{{\alpha }_{\mathrm{in}}}\right)]\widehat{r};\hfill & \mathrm{if}\theta <{\alpha }_{\mathrm{in}},\hfill \end{array}\end{array}$$(7)

where α_out and α_in are the outer and inner jet angles at the boundary between the two velocity laws, respectively. The variables v_α, out and v_α, in are the radial jet velocities at the outer jet angle and inner boundary angle, respectively. Finally, v_M is defined as:

$$\begin{array}{c}\hfill v_{\mathrm{M}}=\frac{v_{\alpha },in-v_{\alpha },out}{cos\left(\frac{\pi }{2}\frac{{\alpha }_{\mathrm{in}}}{{\alpha }_{\mathrm{out}}}\right)}+v_{\alpha },out.\end{array}$$(8)

This model is based on a velocity profile similar to that used by Federrath et al. (2014) who used this nested velocity profile in order to model the momentum transfer and feedback by jets and outflows launched from protostellar disks.

The density profile in this model is equal to the density profile of the disk-wind model, given in Eq. (6).

With hindsight, it would have been best to leave a cavity in both the disk wind and X-wind models, as the original theoretical models in fact prescribe. However, the density law that we use implies a low density outflow close to the centre for these two models, which results in a central outflow.

2.2. Radiative transfer through the jet

To model the amount of absorption by the jet, we determine the loss of intensity of the photospheric contribution of the primary due to scattering of continuum photons from the primary by the hydrogen atoms in the jet. This causes the observed absorption feature in the Hα-profiles of the spectra during superior conjunction. Since the gas in the jet has high outflow velocities, the absorption by the jet lobe that is pointed towards us will be mainly blue-shifted. The jet emission is assumed negligible and we do not include this in our model. Hence, the specific intensity along a ray aligned with the line-of-sight becomes:

$$\begin{array}{c}\hfill I_{\nu }(n,s)=I_{\nu }^0(n){\mathrm{e}}^{-{\tau }_{\nu }(n,s)},\end{array}$$(9)

where $I_{\nu }^0$ is the initial intensity, τ_ν is the optical depth at position s along the nth line-of-sight from the primary component. The photospheric intensity along each line-of-sight from the primary is equally weighted. Since there is a large velocity gradient along the line-of-sight, we implement the Sobolev approximation to determine the optical depth in our medium. Hence, the optical depth will be proportional to:

$$\begin{array}{c}\hfill {\tau }_{\nu }(n,s)\propto \frac{\chi (n,s)\mathrm{\Delta }s}{|\mathrm{d}{v}_s/\mathrm{d}s|},\end{array}$$(10)

where Δs is the length of the line element at position s, $v_s=\mathit{v}·\widehat{\mathit{n}}$ is the projected velocity along the line-of-sight, and |dv_s/ds| is the velocity gradient at position s. The optical depth is determined for each time step along each line-of-sight from the primary component towards the observer. The extinction coefficient χ(n, s) is

$$\begin{array}{c}\hfill \chi (n,s)=\frac{\pi e^2}{m_{\mathrm{e}}c}{f}_{\mathrm{lu}}(n_{\mathrm{l}}-\frac{g_{\mathrm{l}}}{g_{\mathrm{u}}}{n}_{\mathrm{u}}),\end{array}$$(11)

where m_e is the electron mass, f_lu is the oscillator strength, and where l and u denote the lower and upper levels of the transition. We assume hydrogen to be mainly neutral (HI), i.e. n_HI ≈ n_H. According to the Boltzmann equation,

$$\begin{array}{c}\hfill \frac{n_i}{n_{\mathrm{HI}}}\approx \frac{n_i}{n_1}=\frac{g_i}{g_1}exp(-\frac{h\nu }{k_{\mathrm{B}}T}),\end{array}$$(12)

where k_B is the Boltzmann constant, g is the statistical weight that takes into account the degeneracy of the energy states, i ∈ {l, u}, and where we assume that most of the hydrogen is in the ground state (n_HI ≈ n₁). Since the jet consists mainly of hydrogen and helium, we can write the mass density as

$$\begin{array}{cc}\hfill \rho & \approx n_{\mathrm{H}}{m}_{\mathrm{H}}+n_{\mathrm{He}}{m}_{\mathrm{He}}=n_{\mathrm{H}}{m}_{\mathrm{H}}+Y_{\mathrm{He}}4m_{\mathrm{H}}\hfill \\ \hfill & =n_{\mathrm{H}}{m}_{\mathrm{H}}(1+4Y_{\mathrm{He}}),\hfill \end{array}$$(13)

where n_H is the hydrogen number density, m_H is the mass of a hydrogen atom, n_He is the helium number density, m_He is the mass of a helium atom, and Y_He = n_He/n_H. By assuming the jet to be isothermal and in local thermodynamical equilibrium, we can use Eqs. (12) and (13) in Eq. (11) to obtain

$$\begin{array}{c}\hfill \chi (n,s)\propto \rho (n,s).\end{array}$$(14)

From this follows that

$$\begin{array}{c}\hfill {\tau }_{\nu }(n,s)\propto \frac{\rho (n,s)\mathrm{\Delta }s}{|\mathrm{d}{v}_s/\mathrm{d}s|}·\end{array}$$(15)

The resulting intensity will be the contribution of the intensity passing through the jet along each line-of-sight:

$$\begin{array}{c}\hfill I_{\nu }={\displaystyle \sum _n{I}_{\nu }^0(n)\prod _{s}exp(-\frac{c\rho (n,s)\mathrm{\Delta }s}{|\mathrm{d}{v}_s/\mathrm{d}s|}),}\end{array}$$(16)

where c is a scaling factor that is multiplied to the optical depth, such that the resulting intensity in our model can be scaled to the observed intensity in our spectra. The reason for the implementation of the scaling factor c is that this method does not allow us to determine the absolute jet density In the following section we explain how we implement the fitting of each model to the observations so as to decide the best jet model.

3. MCMC-fitting routine

In order to fit our model to the data, we use a MCMC routine. This was done with the emcee package (Foreman-Mackey et al. 2013), a pure-Python implementation of Goodman & Weare’s “Affine MCMC ensemble sampler” (Goodman & Weare 2010).

The free parameters of this fitting routine in our three jet configurations are: inclination angle of the binary orbit, i, jet half-opening angle, θ_out, scaling factor for the optical depth, c, jet velocity along the jet axis, v₀, jet velocity at the edges, v_out, and the radius of the primary, R_prim (in units of the semi-major axis a_prim of the primary). Depending on the adopted model, there are additional free parameters: the inner jet angle at the boundary between two velocity regions in the jet (θ_in for model B and C), the velocity at this boundary (v_in for model B), and the disk-wind scaling factor for the jet velocity (c_v for model C).

Before starting the MCMC routine, we choose 128 initial combinations of the free parameters, where each combination is referred to as a “walker”. These walkers are given a random initial distribution between predetermined boundaries. Some additional constraints are imposed on the walkers, i.e., the outflow velocity in the jet decreases for increasing polar angle, v_out <  v_in <  v₀. Secondly, the primary radius R_prim is kept smaller than 70% of the Roche radius of the primary component. This constraint is derived from the observations of our test object (see Sect. 4.1 for more details). The full list of model parameters and corresponding constraints is given in Table 2.

The walkers will start exploring parameter space from this initial position during the MCMC routine. The acceptance or rejection of a new position during each iteration depends on the log-likelihood of that new position. The log-likelihood is given by:

$$\begin{array}{c}\hfill ln\mathcal{L}=-\frac{1}{2}{\displaystyle \sum _i[\frac{{(y_i-m_i)}^2}{{\sigma }_{\mathrm{tot},i}^2}+ln\left(2\pi {\sigma }_{\mathrm{tot},i}^2\right)],}\end{array}$$(17)

with y_i the data, m_i the model, and σ_tot, i = σ_S/R, i + σ_bs, i the error on the data. The total uncertainty σ_tot, i in our data is the sum of the error on the spectra determined from the signal-to-noise, σ_S/R, i, and the uncertainty that comes from our determination of the background spectra σ_bs, i (see Sect. 4.2). The sum i is over all sampled wavelength points of our spectra. We run the MCMC chains until the walkers converge to a certain value for each parameter. The best-fitting parameters for each MCMC routine is chosen by fitting a Gaussian profile to the posterior probability distribution, which yields a mean and standard deviation for each parameter. We only use the iterations following the burn-in phase to determine the mean and standard deviations from the posterior distribution.

4. Using observations of a post-AGB star with jets to constrain our jet model

4.1. The post-AGB binary IRAS 19135+3937

To test this fitting methodology for jets launched by post-AGB binaries, we start with the system IRAS 19135+3937. the data for this object are part of a large radial velocity monitoring programme of post-AGB binaries by Van Winckel et al. (2009), which started in 2009 with the HERMES spectrograph, mounted on the 1.2 m Mercator telescope, La Palma, Spain (Raskin et al. 2011). We choose post-AGB binary IRAS 19135+3937 as test object for our fitting code, because we have a well-sampled orbital cycle for this object. Since the start of the monitoring campaign, we have obtained 90 spectra for this system, providing us with a good phase coverage of the orbital period. This post-AGB binary has a relatively short orbital period of 126.97 ± 0.08 d and an eccentricity of 0.13 ± 0.03 (Gorlova et al. 2015; Oomen et al. 2018). The orbital parameters of the system are listed in Table 3. Ideally, we would be able to derive the radius of the post-AGB star by combining the distance to the object with its photometric flux. However, in the case of IRAS 19135+3937, the complications are twofold. First, the distance from Gaia DR2 (Gaia Collaboration 2016, 2018) for IRAS 19135+3937 is an approximate estimate of the true distance, because the parallaxes in Gaia DR2 are determined using single-star models. Additionally, its orbital period (∼ 127 days, see Table 3) implies that the Gaia DR2 parallax for this star will have a similar amplitude as its projected orbit, thereby reducing the accuracy of the resulting distance estimate. Second, this object is a semi-regular variable, which implies that its flux shows significant variability. This makes it impossible to accurately determine the stellar radius of the post-AGB star. Therefore, for the modelling of this system, we keep the stellar radius as a free parameter between set boundaries.

The boundaries for the stellar radius are constrained based on its Roche-lobe filling factor³. A star that nearly fills its Roche lobe would be tidally distorted, causing ellipsoidal variations in its light curve (Wilson & Sofia 1976), which we do not observe for IRAS 19135+3937. This implies that the radius of the post-AGB star must be smaller than its Roche lobe. In order to set an upper-limit for this parameter, we look at the Roche-lobe filling factor of sequence E red giants, which are binary systems that show ellipsoidal variations in their photometric data. From the sample of sequence E red giants in the study of Nicholls et al. (2010), we derived an average Roche-lobe filling factor of 0.82  ±  0.12 R_RL. Based on these data, we set the upper boundary for the radius of the post-AGB star at 70% of the Roche radius.

Figure 3 shows the dynamic spectra of IRAS 19135+3937 made by 22 out of the 90 spectra. We choose to use these 22 spectra since they were taken during one cycle of 129.98 d. The dynamic spectra are created by showing the observed Hα line as a function of orbital phase, where we interpolated along the orbital phase for regions where we do not have observations. The absorption feature caused by the scattering of the evolved star’s light in the jet is clearly visible between orbital phases θ = 0.3 − 0.7. The depth, extent and duration of this absorption feature is dependent on the orbit of the binary system, the size of the evolved star, and the size, density, and velocity of the jet.

Fig. 3. Interpolated dynamic spectra for the Hα-line of IRAS 19135+3937, based on the 22 spectra (see Sect. 4.1). The dashed grey line indicates the systemic velocity (γ) of the binary system and the solid line indicates the RV of the post-AGB star.

4.2. Determination of the unobstructed Hα spectral profile

Our model calculates the amount of absorption of background light by the jet at each wavelength and returns the resulting synthetic spectrum. Hence, we have to provide our model with an initial background spectrum that does not include any absorption by the jet. Determining this initial Hα spectrum is not straightforward, because the background spectrum is more complex than just the photospheric Hα profile.

The first component is the photospheric Hα absorption line from the post-AGB stellar spectrum. A synthetic stellar spectrum from Coelho et al. (2005) with stellar parameters T_eff = 6000 K, log g = 1.0, and [M/H] = − 1.0 was used to replicate this photospheric component (Gorlova et al. 2015). The absorption component is shifted at each orbital phase, according to the RV of the primary component.

In addition to the photospheric absorption, we observe an emission feature that is present throughout the whole orbital phase and is superimposed on the photospheric absorption from the AGB star. This feature varies between a double-peaked emission profile during inferior conjunction and a single, red-shifted emission peak during superior conjunction. During superior conjunction, the blue-shifted emission peak disappears at least in part due to the absorption caused by the jet. The emission component is not centred on the radial velocity of the post-AGB star nor that of the companion (see Fig. 3) so it is unclear if the emission is from the circum-companion accretion disk, from the base of the jet itself, or from yet another source.

If the emitting source is located behind the jet during the whole orbital phase, it must be included in the initial background spectra that we use to calculate the absorption. If the emitting source is not located behind the jet, on the other hand, the double-peaked emission profile should not be included in the background spectra, but should rather be added up to the resulting spectra after the absorption by the jet is determined. In addition, we need to determine if this emission feature is constant in strength or if it varies throughout the orbital period.

We tested three possible situations, illustrated in Fig. 4: the emitting source is never blocked by the jet and is constant, the emitting source is blocked by the jet during superior conjunction and is constant, and the emitting source is blocked by the jet during superior conjunction and varies throughout the orbital phase.

Fig. 4. Scenarios to determine the contribution of the emission component to the spectra. Upper row: case I, only the photospheric Hα spectrum from the primary is absorbed by the jet (denoted by “$\stackrel{\text{jet}}{\to }$”). The double-peaked emission component is constant and independent and is thus added (denoted by “+”) to the resulting spectrum. Middle row: case II, the background spectrum is the combined photospheric Hα line plus the emission component. This light is then absorbed by the jet. Lower row: case III, This is the same as case II, except for the fact that the emission component varies as a function of orbital phase. The emission component shown here is taken at phase ϕ = 0.5, and thus weaker than the emission component at phase ϕ = 0 (as can also be seen in Fig. 5).

Case I: The emitting source is not blocked by the jet and is constant in strength. There are three components in each observed spectra: the photospheric absorption line, the double-peaked emission feature, and the absorption by the jet. The emission feature remains constant in this scenario, since it is not blocked by the jet. As shown in Fig. 4, the background spectrum is the photospheric Hα line of the post-AGB star. The absorption by the jet will thus only absorb light from this spectrum. Since the emitting source is not blocked by the jet, the resulting observed spectrum will be the sum of the absorbed spectrum and the emission component.

However, when we look at the absorbed spectrum for case I in Fig. 4, we notice that the flux goes below zero at certain wavelengths. This is not possible, since the only source of light absorbed by the jet is the primary’s photospheric light, the jet cannot absorb more than this light coming from the photospheric line. Hence, this scenario where the emission feature is not blocked by the jet is unphysical and we conclude that the emission region must be completely or partially behind the jet during superior conjunction.

Case II: The emitting source is blocked by the jet and constant in strength. In this scenario, we assume the double-peaked emission feature remains constant in flux, but some of its light is absorbed by the jet. The variation in strength that we observe in the dynamic spectra of Fig. 3 is thus only caused by the absorption by the jet. We tested this scenario by combining the photospheric Hα line with the strongest double-peaked emission component that we observed in the spectra, as is shown in the middle row of Fig. 4. The fitting routine will then calculate the amount of absorption by the jet and fit it to the data. This gives a very poor fit to our observations. The best-fitting jet solution is a jet that has an extremely high optical depth to account for the high amount of absorption needed in this case, which is not realistic since the jet is optically thin. We conclude that we are not able to obtain good fits with this set-up.

Case III: The emitting source is blocked or partially blocked by the jet and varies in strength. We assume an intermediate situation, where the observed intensity of the emitting source varies throughout the orbital period. In order to determine the background spectrum at every orbital phase, we decided to interpolate the region between orbital phase ϕ ∈ [0.15 − 0.8] in the dynamic spectra, as is shown in Fig. 5. In this way, we assume that the background spectrum varies smoothly between phase 0.15, when the emitting source is presumably unobstructed, and phase 0.8, when the jet has completed its passage in front of the source. Since we do not know exactly how much emission the spectra show between these phases, we used two methods to interpolate this region: a linear interpolation (left plot of Fig. 5) and a cubic interpolation (right plot of Fig. 5). Both of these phase dependent spectral series are used in our model-fitting code as background spectra from which we subtract photons.

Fig. 5. Example of the determination of the spectra to be used as background light absorbed by the jet (case III). The regions between orbital phase ϕ ∈ [0.15 − 0.8] are substituted by values interpolated between adjacent spectra using a linear interpolation (left) or a cubic interpolation (right).

Fig. 6. Strongest emission spectrum for case II and background spectrum used for case III (upper). The uncertainty on the determination of the background flux is the difference between the flux of the strongest emission spectrum and the flux of the spectrum used for case III, where the emission varies throughout orbital phase (lower).

We include an extra uncertainty from the determination of the background spectrum, σ_bs, in the total uncertainty used in Eq. (17). This uncertainty is related to the strength of the emission component at each wavelength bin. The uncertainty in the flux will be larger for wavelength bins where this emission component is stronger. We quantify this by using the difference between the strongest observed emission flux and the flux of the background spectrum in case III as an extra uncertainty, as is shown in Fig. 4. The bottom panel of Fig. 4 shows how the uncertainty is larger for wavelength bins with a stronger emission flux.

5. Spatio-kinematic modelling of the jet

We recreate the time-series of the Hα spectral profile for IRAS 19135+3937 using the three jet configurations described in Sect. 2 in our model, in order to determine the spatio-kinematic structure of the jet. Several fits are performed for different values of the factor p in the angle-dependent density profile in the stellar jet model, and the factors p_in and p_out in the X-wind and disk-wind model. These combinations are shown in Table 1. We created synthetic spectra that are used as background spectra for our fitting (see Sect. 4.2). The synthetic Hα spectra contain the photospheric absorption line from the evolved star plus the double-peaked emission feature interpolated between orbital phase ϕ ∈ [0.15 − 0.8] using a linear and a cubic interpolation, as described in Sect. 4.2. We eventually used these two dynamic background spectral series for our model fitting. By comparing the goodness of fit between the models that use the linearly interpolated background spectra and those using background spectra obtained by cubic interpolation, it is clear that the latter give better results for all runs. Hence, we only focus on the models that use the background spectra obtained by cubic interpolation. In the following, we further discuss our best-fitting synthetic spectra for the three jet configurations.

5.1. The best-fitting jet model

All results of our MCMC model-fitting routine are listed in Table A.1, including their density structure, best-fitting parameters, and goodness of fit (reduced χ² and Bayesian information criterion). The values of the model parameters are determined from the posterior density distribution of our MCMC-fitting routine. In order to extract these values and their corresponding standard deviations, we fit a Gaussian curve to the posterior density distribution. The mean and standard deviation of this Gaussian curve are reported alongside the best fit values in Table A.1.

Since we are comparing the fitting results of different configurations with a different set of parameters, we evaluate the relative goodness of fit of each model by comparing their Bayesian information criterion (BIC), which is defined as:

$$\begin{array}{c}\hfill \mathrm{BIC}=ln(n)k-2ln\mathcal{L},\end{array}$$(18)

with n the number of data points and k the number of model parameters. The BIC takes into account the number of parameters in a model. Hence a higher number of parameters gets penalised by a higher BIC value.

The best-fitting model is the model using the stellar jet configuration with the background spectrum obtained by cubic interpolation and p-factor p = 8 (see Table 2 for corresponding best-fitting parameters). The posterior density distribution of the model is represented in the corner plot of Fig. 7. The individual spectra and dynamic spectra of this model are shown in Figs. 8 and 9, respectively. The jet has a very wide half-opening angle of θ_in = 76°.

Fig. 7. Visualisation of the one- and two-dimensional projection of the posterior density distribution of our model parameters for the best-fitting model using the corner software by Foreman-Mackey (2016). The red lines indicate the values for the best-fitting parameters (see Table 2).

Fig. 8. Observations (blue) vs. best-fitting model spectra for model A (red). The panels increase in orbital phase from top to bottom and left to right.

Fig. 9. Dynamic spectra created from the observed spectra (left) and the model spectra (right). The model dynamic spectra where created using the best-fitting model with the stellar jet configuration. The encircled region in the model spectra highlights the extra red-shifted absorption feature that is created by the model, but which we do not observe, as explained in Sect. 5.1.

The jet for this model reaches a maximum velocity of v₀ = 1210 km s⁻¹ along the jet axis. Due to the large inclination angle of i = 79°, the projected jet velocities reach blue-shifted values of ≈400 km s⁻¹, as is the case for the observations. In order to determine the nature of the companion, we could compare these jet velocities with typical escape velocities of stellar objects. For a 1 M_⊙ white dwarf, the escape velocity is about 6500 km s⁻¹. For a 1 M_⊙ main sequence star, this is about 618 km s⁻¹. Since the observed jet velocities are closer to the escape velocity of a main sequence star, we conclude that the companion in this system is most likely a main sequence star. This conclusion is also supported by Oomen et al. (2018), who found that based on the initial mass of the companion it is likely a 1 M_⊙ main sequence star.

The absorption feature caused by the jet covers a large part of the orbital phase (between ϕ ∈ [0.25 − 0.75]). This suggests that the difference between orbital inclination and jet half-opening angle should be rather small. Since the absorption is not seen throughout the whole orbital phase, the jet angle is most likely smaller than the inclination angle. This is also the case In the results of our best-fitting model, where the jet has a half-opening angle of θ_out = 76° and the inclination angle is slightly larger (i = 79°). Figures 8 and 9 show a fairly good match between the model spectra and the observations. An interesting feature that appears in the model spectra, but not in the observed spectra is the red-shifted absorption feature, indicated in Fig. 9 by the dashed circle. This absorption feature in the model is due to absorption by the jet in the region where the light rays enter the jet. The jet particles in this part of the jet have a positive velocity when projected on to the line-of-sight through the jet. Hence, the absorption will be red-shifted. This absorption is not as significant in the observations as in the model calculations.

The velocity and density structure of this model, represented in Figs. 10 and 11, show that the jet is very dense at the edges and has an extremely low density along the jet axis.

Fig. 10. Velocity and density structure as a function of jet angle for the best-fitting model. The density is normalised such that ρedge = 1.

Fig. 11. Velocity field and density structure for the best-fitting model of the jet in IRAS 19135+3937. The colour of the density structure is represented in a logarithmic scale.

5.2. Comparison of the three jet models

Here, we investigate which of the three models gives the best representation of the jets in our post-AGB binary systems. In order to differentiate which model fits the data best, we need to look at the difference between the BIC values of the models, where a lower BIC value translates into a better fit. We find the goodness of fit for the best-fitting model of the stellar jet (model A), the disk wind (model B), and the X-wind (model C) to be BIC_A = −24 878, BIC_B = −24 867, and BIC_C = −24 805, respectively. According to Kass & Raftery (1995), a difference of |BIC₁−BIC₂| > 6 is a strong condition to determine that the model with the lowest BIC value fits the data better. This only works when the absolute BIC values are completely correct. In our case, the relative BIC values are reliable, but the absolute BIC values are not. The BIC value is highly dependent on the errors in the χ² calculations. Our error only takes into account the error on the spectra (determined from the signal-to-noise ratio) and the uncertainty from the determination of the initial background spectra, which is given as input to the model.

Thus, in order to compare the different models, we compare both their BIC values and the reduced χ² values in order to determine the best-fitting models. By doing so, it can be seen that model A fits best to the observations. This model only gives a slightly better fitting result than models B and C, however. The reduced χ² for model B is even slightly better than that of model A (${\chi }_{\text{red},\text{B}}^2=0.2428$ and ${\chi }_{\text{red},\text{A}}^2=0.2434$). If we would neglect that model B has more model parameters than model A (8 compared to 6), this would be the better fitting model. However, since we do want to penalise for extra model parameters, we conclude that model A has a slightly better fit than model B (ΔBIC = 11).

The BIC of the best-fitting X-wind configuration (model C) is higher than model A with ΔBIC = 73. For this model, the model parameters are also very similar to model A. The main difference between model C and model A is a smaller inclination angle and jet angle. These two angles are smaller by about 3° for model C.

In Fig. 12, we show the velocity profile and density profile in the jet for best-fitting models A, B, and C. This shows again that all three models tend to converge to a similar structure. The inner jet boundary of model B and C in particular tends to lie very close to the outer jet edge. Hence, for these model results, the jet is mainly governed by the inner velocity and density law. The “outer jet” region is a relatively small, negligible region in these models.

Fig. 12. Velocity and density structure as a function of jet angle for the best-fitting model for the stellar jet model (A), the disk wind model (B), and the X-wind model (C). The density is normalised such that ρedge = 1.

The half-opening angle of the jet lies between 72° and 77° for all three models. The velocity along the edges of the jet is about 10 − 20 km s⁻¹ and reaches maximal velocities of 950 − 1350 km s⁻¹ at the centre. The best fitting density law is the same for models A and B. Both model B and model C converge towards a similar simple structure as the stellar jet of model A. We expect that for some post-AGB binary systems, such as IRAS 19135+3937, the three models will converge to a similar structure, whilst for other systems, the three models will be easier to differentiate results. This is due to the large diversity in terms of orbital parameters and system sizes in these systems.

The best-fitting model for the stellar jet configuration has a density parameter p = 8. This configuration assumes a simple conical outflow with an increasing density from the pole towards the jet edge. The best-fitting model for the disk wind configurations has a similar density law in the inner jet (p_in = 8). The outer density law in model B and model C is p_out = −2. The inclination for all three configurations converges to ∼75 − 80°, which implies that the system is observed nearly edge-on. The study by Gorlova et al. (2015) has suggested this as well. Their conclusion was based on the semi-regular variability of the system, first observed by Sallman & Droege (2004). Gorlova et al. (2015) showed that the light variations are most likely caused by the partial obscuration of the evolved star by the puffed-up circumbinary disk during inferior conjunction. This obscuration can only occur if the inclination of the system is rather high, which is what we find.

Similarly to the inclination, the jet velocity at the edge and the jet angle for each of the three configurations converges to the same values. This means that the density and velocity structure in the different configurations will be alike as well, as can be seen in Fig. 12.

By knowing both the orbital parameters and the inclination of the orbital system, we can calculate the radius of the evolved star and the mass of the companion. Given the mass function and the inclination found in the best-fitting model (i = 79°), we calculate the mass of the companion star to be 0.44 M_⊙. The exact radius of the evolved star was a model parameter in our code, given in units of the semi-major axis of the primary. The result for radius of the evolved star in the best-fitting model is 0.66 ± 0.04 ⋅ a_prim, with a_prim the semi-major axis. Hence, given the projected semi-major axis a_prim ⋅ sin i (Table 3) and the inclination found from our model-fitting routine, the radius of the evolved star will be R_prim = 30 ± 3 R_⊙.

5.3. A possible cavity in the jet

We find that the density in the inner part of the jet is extremely low, compared to the outer part of the jet (see Fig. 10). This leads us to question whether the inner region contributes at all to the absorption caused by the jet. We test this by including a cavity (region of zero density) in the best-fitting jet model. We increase the jet cavity from 0° up to 50°, and compared the BIC values of these tests with the best-fitting model. By increasing the cavity to have opening angles from 0° to 30°, the synthetic spectra, and thus the goodness of fit, does not change. The density in this region in the best-fitting model is just too low in order to contribute to the absorption feature in the spectra. By further increasing the cavity in the jet from 30°, the BIC improves slightly and reaches the same value again when the cavity is increased up to ∼40°.

By including this cavity in our model, we show that the inner region of the jet (< 40°) has extremely low density. In our best-fitting models for both the X-wind and the disk wind, the bulk of the matter is launched at an angle > 40°. This is in accordance with the disk wind theory by Blandford & Payne (1982), where the magnetic field line and the jet axis must make an angle of > 30°. This does not mean that the inner regions are completely void of matter, but rather that the density is too low to have any effect on the Hα line profile. This consequently impacts the maximum outflow velocity in the jet. Assuming that the inner region does not contribute to the observed absorption feature, since it is too low a density, we can conclude that the maximum velocity in the jet is reached at an angle of about 40° (the innermost region of the jet that launches matter). For the best-fitting stellar jet model, this means that the jet reaches velocities of about 870 km s⁻¹, instead of 1210 km s⁻¹.

Interestingly, the existence of a cavity seems to be implied by the data. In future work, we will alter the density laws, such that the disk wind and X-wind models are closer to their initial conception launching gas in a cone geometry without a central outflow.

6. Summary and future work

We developed a fitting procedure that allows us to determine the spatial, velocity, and density structure of jet-creating post-AGB binaries. Our procedure implements three separated jet models: a stellar jet configuration (Model A), a disk wind configuration (Model B), and an X-wind configuration (Model C); which are based on jet launching mechanisms designed to explain YSO jets. Our procedure calculates synthetic spectral lines for these systems by calculating how much of the evolved star’s light is scattered by the jet as it passes between the post-AGB star and the observer. These synthetic spectral lines were then fitted to the observations. This shows that the companion is likely a main sequence star.

We fitted synthetic lines to the data of the post-AGB binary IRAS 19135+3937 for each of the three jet configurations. The best-fitting model is the stellar jet configuration, where the jet is governed by the density law: ρ ∝ θ⁸. The jet in this model is a very wide outflow of 76°, with a slow velocity component along the edges (v_out = 11 km s⁻¹) and a high-velocity component in the inner regions. The inner region in the jet has extremely low density, making the contribution of this region to the scattering in our model insignificant. The highest velocity in the jet is about 870 km s⁻¹. There is no significant difference between the results of the three jet configurations. Both the X-wind and the disk wind model converge to a solution similar to the stellar jet configuration. We conclude that our model succeeds in fitting the dynamic spectra well.

Unexpectedly, we detected an emission component that is of unknown origin. Our fits indicate that this source varies over the orbit and that this light is absorbed when the companion and its jet are at inferior conjunction.

In future work, we will apply our fitting procedure to our target sample of about 15, jet-creating, post-AGB binaries. We will apply this generic model to determine the configurations of all objects and study if we can infer additional information on the jet launching mechanisms for these systems. We also expect to be able to measure the actual jet density more accurately, something that will provide us with the accretion rates. Ultimately, we would like to be able to connect accretion rates, AGB star abundances, and circumbinary disk longevity in one narrative that will shed light on the origin of post-AGB binaries.

Acknowledgments

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This work was performed on the OzSTAR national facility at Swinburne University of Technology. OzSTAR is funded by Swinburne University of Technology and the National Collaborative Research Infrastructure Strategy (NCRIS). DK acknowledges the support of the Macquarie University New Staff funding. HVW acknowledges support from the Research Council of the KU Leuven under grant number C14/17/082. Based on observations obtained with the HERMES spectrograph, which is supported by the Research Foundation – Flanders (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.

References

Appendix A: Table with all fitting results for the three models.

In Table A.1 we present the best fitting parameters with corresponding reduced chi-squared and BIC for the fitting of the three jet models to IRAS 19135+3937. ΔBIC denotes the difference of the BIC with the BIC of the best-fitting model. The jet models are model A (the stellar jet model), model B (the disk-wind model), and model C (the X-wind model).

All Tables

All Figures

	Fig. 1. Mesh grid of the primary star component and the rays going through the jet. The primary component is sampled as a Fibonacci grid with 200 grid points. The absorption of the continuum light by the H-atoms in the jet is determined for each ray along the line-of-sight that departs from these points on the primary’s surface.
In the text

Fig. 2. Upper panel: schematic sketch of the jet model, showing our geometric implementation of the three models. Lower panel: launching region and velocity field of the jet for the stellar jet (A), the disk wind (B), and the X-wind (C). Unlike the stellar jet and X-wind, the launching region for the disk wind is the region between the inner-disk radius din and the outer disk radius dout. This configuration is based on the disk wind configuration applied by Kurosawa et al. (2006).

In the text

	Fig. 3. Interpolated dynamic spectra for the Hα-line of IRAS 19135+3937, based on the 22 spectra (see Sect. 4.1). The dashed grey line indicates the systemic velocity (γ) of the binary system and the solid line indicates the RV of the post-AGB star.
In the text

Fig. 4. Scenarios to determine the contribution of the emission component to the spectra. Upper row: case I, only the photospheric Hα spectrum from the primary is absorbed by the jet (denoted by “$\stackrel{\text{jet}}{\to }$”). The double-peaked emission component is constant and independent and is thus added (denoted by “+”) to the resulting spectrum. Middle row: case II, the background spectrum is the combined photospheric Hα line plus the emission component. This light is then absorbed by the jet. Lower row: case III, This is the same as case II, except for the fact that the emission component varies as a function of orbital phase. The emission component shown here is taken at phase ϕ = 0.5, and thus weaker than the emission component at phase ϕ = 0 (as can also be seen in Fig. 5).

In the text

	Fig. 5. Example of the determination of the spectra to be used as background light absorbed by the jet (case III). The regions between orbital phase ϕ ∈ [0.15 − 0.8] are substituted by values interpolated between adjacent spectra using a linear interpolation (left) or a cubic interpolation (right).
In the text

	Fig. 6. Strongest emission spectrum for case II and background spectrum used for case III (upper). The uncertainty on the determination of the background flux is the difference between the flux of the strongest emission spectrum and the flux of the spectrum used for case III, where the emission varies throughout orbital phase (lower).
In the text

	Fig. 7. Visualisation of the one- and two-dimensional projection of the posterior density distribution of our model parameters for the best-fitting model using the corner software by Foreman-Mackey (2016). The red lines indicate the values for the best-fitting parameters (see Table 2).
In the text

	Fig. 8. Observations (blue) vs. best-fitting model spectra for model A (red). The panels increase in orbital phase from top to bottom and left to right.
In the text

	Fig. 9. Dynamic spectra created from the observed spectra (left) and the model spectra (right). The model dynamic spectra where created using the best-fitting model with the stellar jet configuration. The encircled region in the model spectra highlights the extra red-shifted absorption feature that is created by the model, but which we do not observe, as explained in Sect. 5.1.
In the text

	Fig. 10. Velocity and density structure as a function of jet angle for the best-fitting model. The density is normalised such that ρedge = 1.
In the text

	Fig. 11. Velocity field and density structure for the best-fitting model of the jet in IRAS 19135+3937. The colour of the density structure is represented in a logarithmic scale.
In the text

	Fig. 12. Velocity and density structure as a function of jet angle for the best-fitting model for the stellar jet model (A), the disk wind model (B), and the X-wind model (C). The density is normalised such that ρedge = 1.
In the text