© The Authors 2024

1 Introduction

A star with an initial mass ranging from 0.8 to 8 M_⊙ will eventually evolve through the asymptotic giant branch (AGB) stage. During this phase, the star undergoes a process of significant mass loss (10⁻⁷ to 10⁻⁴ M_⊙ yr⁻¹) through a slow (5–20 km s⁻¹) wind that enriches the interstellar medium (Habing & Olofsson 2004). When leaving the AGB phase, its effective temperature increases and goes through a short-lived transition of post-AGB and pre-planetary nebula (lasting about a few thousand years) phase. The majority of planetary nebulae (PN) show a variety of aspherical (bi- or multi-polar) structures (Balick & Frank 2002).

Those bipolar structures are expected to arise from binary interaction with a close-by companion (Jones & Boffin 2017). In the preceding phase, the post-AGB phase, binarity has been widely detected and studied (Oomen et al. 2018 and references therein). Common building blocks of post-AGB binaries have been identified (e.g. Kluska et al. 2022), consisting of a central binary surrounded by a stable circumbinary disk of gas and dust. They often exhibit a bipolar outflow which is likely launched by an accretion disc around the companion (e.g. Bollen et al. 2022 and references therein). Thanks to the near-infrared (NIR) interferometry, these building blocks can also be spatially resolved (e.g. Kluska et al. 2019; Corporaal et al. 2023). While these post-AGB binaries represent about a third of all known Galactic post-AGB stars (Kluska et al. 2022), they are poorly connected to the AGB progenitor sample, as evidence for binaries among these progenitors is lacking. Due to the intrinsic properties of AGB stars, detecting binaries is not a straightforward task. Given their photosphere is bright and dynamic, the detection of orbital motion is possible only if it exceeds the amplitude of the intrinsic pulsation of the stellar envelope (Jorissen 2004).

The late-type carbon star V Hydrae (V Hya) is among the best candidates to study the impact of a companion on the circumstellar environment of AGB stars. Secular photometric monitoring indicates that the star exhibits two photometric periods: a long period of 6160 ± 400 days with an amplitude of 3 mag and a short period of 530 ± 30 days with an amplitude of 1.5 mag (Knapp et al. 1999). While the 530-day period is commonly associated with the Mira pulsation of the AGB star, the origin of the 17 yr period is still under debate. Several scenarios have been suggested, many of which invoke the presence of a binary companion: eclipse by a dusty cloud surrounding a companion (Knapp et al. 1999), extinction by the plasma ejection of the companion (Sahai et al. 2016), a triple system with an inner eccentric orbit of 8 yr (Salas et al. 2019), along with possible episodic dust formation (Lloyd Evans 1997). In addition, the presence of a binary companion has been inferred on the basis of the near-UV (NUV) excesses (Sahai et al. 2008), detections of blue emission lines in the spectrum (Lloyd Evans 1991, 2009), and the exceptional structure of the large-scale environment of the star, as we describe below.

Several radio observations in the CO bands highlight the presence of an environment shaped by different velocity regimes: a slow (15 km s⁻¹) dense equatorial wind, a medium-velocity (45 km s⁻¹) spherical wind, and a fast (150–200 km s⁻¹) bipolar outflow (Knapp et al. 1997). Kahane et al. (1996) mapped the fast wind component in the CO J = 2−1 transition as a bipolar flow whose blue-shifted lobe is pointing toward the east with an inclination axis of about 30°. The presence of an eastward high-velocity knotty jet was also detected in the [S II] emission by HST/STIS (Sahai et al. 2016) and its possible radio-equivalent by recent high-resolution ALMA observations together with the discovery of several concentric rings inside the dense equatorial structure (Sahai et al. 2022). All the reported observations confirm the large-scale complexity of the system, although no model (so far) has been able to assemble all the pieces of the puzzle and link the large-scale bipolar structure with the expected binary motion.

In this paper, we bring on new observations from the high-resolution HERMES spectrograph (Raskin et al. 2011). The long-term spectral monitoring of V Hya provides unprecedented information on the system dynamics and new constraints on its binary motion. We report radial-velocity and spectral-line variations in phase with the long-term photometric period and interpret them in the framework of a binary system with a period of 17 yr. The results are twofold: they provide the first spectroscopic evidence of a binary-AGB system and highlight the complex phenomena of the short-lived transition phase of V Hya, making it a case study for further investigation.

This paper is structured as follows. The observations and data reduction are presented in Sect. 2. The determination of the orbital parameters is described in Sect. 3. In Sect. 4, we compare the spectroscopic signal with the visual photometry. In Sect. 5, we analyze the temporal variation of selected spectral lines. The spatio-kinematic modeling of the Na I resonance doublet is introduced in Sect. 6. In Sect. 7, we discuss the possible origin of the long obscuration event and put our new spectroscopic results in the context of previous observations. In Sect. 8, we outline our main findings.

2 Observations

2.1 Spectroscopic data

V Hya has been monitored with the HERMES spectrograph mounted on the 1.2 m Mercator telescope at La Palma (Raskin et al. 2011). A total of 130 spectra were taken on different dates between 2011 and 2023, with an average exposure time of 900 s. The high-resolution mode was used, providing a resolution of R = 86000 over a wavelength range from 3800 to 9000 Å. The dedicated python-based pipeline was used for the data reduction. The radial velocity (RV) was obtained by cross-correlating the object spectrum with a carbon-star template and fitting a Gaussian function to the mean line profile, as represented by the resulting cross-correlation function (CCF). The mask used to compute the CCF is built from a HERMES spectrum of the R-type carbon star BD +02°4338, Doppler-corrected to zero velocity from a correlation with the master Arcturus template. All lines with a depth getting below the relative level of 0.8 with respect to the local continuum were kept in the mask. This mask, covering the wavelength range 476.9–654.7 nm, was selected because it provides the cleanest CCF for V Hya (deepest and least asymmetric). Thirty spectra were excluded from the orbital analysis because they bear the signature of shocks. Due to the presence of shocks linked to the stellar pulsation, the CCF profile at those epochs close to maximum light was distorted (asymmetric broadening or double-peaked; Alvarez et al. 2001). Table A.1 summarises the log of the observations, including the observing dates of each spectra, the computed RVs and their formal uncertainties (due to the quality of the CCF fit). An uncertainty of 0.07 km s⁻¹ expressing the long-term stability of the spectrograph is added quadratically to the formal uncertainty of the CCF fit. This long-term accuracy is estimated from the average RV standard dispersion of a sample of RV standard stars monitored along with the science targets (see Fig. 10 in Raskin et al. 2011).

2.2 Photometric data

For more than a century, V Hya has been monitored by variable-star observers. We use the archival visual photometric data from the American Association of Variable Stars Observers, AAVSO (Kloppenborg 2023). The overall dataset (9000 measurements) covers a time interval from 1913 to 2023, allowing us to monitor six 17-yr cycles.

2.3 Astrometric data

Astrometric measurements from both HIPPARCOS (Perryman et al. 1997) and Gaia (Gaia Collaboration 2022) ESA missions have been used to constrain additional parameters that are not accessible from spectroscopic data alone, such as the orbital inclination. The difference between the instantaneous proper motion computed by the two missions and the average proper motion between the two reference epochs provides an estimate of the acceleration in the sky plane. That acceleration being non-null for V Hya (Brandt 2021) is a signature of a long-trend motion and gives additional constraints on the orbital solution.

3 Radial-velocity analysis

3.1 Methods

In this section, we describe the method used to retrieve the orbit from the RV curve. In a first step (hereafter, step 1), the contribution associated to the pulsation was fitted together with a Keplerian orbit. During a second step (hereafter, step 2), the cleaned (pulsation removed) RV curve was combined with the described astrometric data and fitted with the ORVARA¹ code based on a Markov chain Monte Carlo (MCMC) method (Brandt et al. 2021).

Fig. 1 Radial velocities of V Hya. Top: RV curve as a function of time. The full line is the Keplerian orbit and the dashed line is the composite signal (orbit and pulsation). Middle: phase-folded cleaned orbit of V Hya. Bottom: residual values with the ±3σ from the mean is represented as the grey-shaded region.

3.1.1 Disentangling the pulsation (step 1)

The RV curve exhibits two periodic variations: a short-period signal corresponding to the Mira pulsation (following the 530 d periodicity) as well as a long-term trend (see Fig. 1). The long trend can be attributed to the 17 yr modulation observed in the light curve (Knapp et al. 1999) but, as this period is longer than the time spanned by the HERMES monitoring, an entire cycle is not visible. As mentioned in Sect. 2.1, the errors in the RV from HERMES consist of two terms: the error intrinsic to the CCF fit (photon noise and template mismatch) and a 0.07 km s⁻¹ uncertainty to account for a possible offset of the wavelength zero-point. To better constrain the fit, nine additional RV values obtained by Barnbaum et al. (1995) have been added to the RV curve obtained with the HERMES spectrograph. Fortunately, those data, falling on the previous orbital cycle, correspond to different orbital phases as compared to the HERMES data, allowing us to cover almost the whole orbit. To account for an eventual offset of the RV zero point reference between the two instruments, an additional 0.2 km s⁻¹ uncertainty is added to the RV data from HERMES and represents the minimum value required to increase the relative statistical weight of the oldest RV data from Barnbaum et al. (1995) to a significant level, given their higher uncertainties (± 1.5 km s⁻¹).

The RV curve was fitted by a combination of two functions: a Keplerian orbit and an asymmetric sine, sin_a. The asymmetric sine is expected to mimic the S-curve observed in pulsating stars (Jorissen 2004) and is mathematically defined as:

$$A{\sin }_a\left({\omega }_0\left(t+t_0\right)\right)=\frac{A}{\text{Γ}}\arctan \left(\frac{\text{Γ}\sin \left({\text{ω}}_0\left(t+t_0\right)\right)}{1-\Gamma \cos \left({\omega }_0\left(t+t_0\right)\right)}\right),$$(1)

where ω₀ = 2π/P is the angular frequency, A is the signal amplitude, t₀ the time reference and Γ the tilt parameter satisfying |Γ| ≤ 1. Γ = 0 corresponds to the usual sine function of period, P, whereas Γ = +1 (resp. −1) corresponds to the saw-tooth function tilted to the left (right, resp.).

In the fitting process, the asymmetric sine period was fixed at 530 days and the Keplerian period was forced to remain in the range 6310±1220 days. Those values were obtained by computing the Lomb–Scargle periodogram of the light curve and taking the 3σ uncertainty interval for the long-period signal (see Fig. B.1). In total, the fitting routine is composed of a nine-parameter function (3 for the asymmetric sine and 6 to describe the orbital motion).

After a first fitting, the contribution of the pulsation is subtracted from the RV and the cleaned signal is fitted with a Keplerian orbit alone using a χ² -minimisation. The method allows us to get precise orbital parameters at the price of imposing the period. The standard deviations of the orbital parameters are obtained by computing the square root of the diagonal elements of the covariance matrix.

3.1.2 Orbital parameters with ORVARA (step 2)

As a complementary analysis, meant to confirm our results and to obtain the full orbital solution (including the inclination), another independent method was used on the RV curve cleaned from the pulsation. The software package ORVARA, developed by Brandt et al. (2021), uses a parallel tempering MCMC method to simultaneously fit radial-velocity and astrometric measurements. Combining the measurements allows the companion mass and inclination for binary systems to be computed (Escorza & De Rosa 2023), which would otherwise be inaccessible with just the spectroscopic measurements. The code first fits the RV data then adds the absolute astrometric measurements: astrometric position (α, δ), proper motion (μ_δ, μ_δ), and Gaia DR3 parallax for each MCMC step. It allows us to simultaneously fit the full orbital parameters (period, P; eccentricity, e; semi-major axis, a; inclination, i; argument of the periastron of the visible star, ω; time of periastron passage, T₀; angle of the ascending node Ω); the mass of the two stars; and an additional radial-velocity jitter, σ_vr. Details on the numerical implementation can be found in the reference paper (Brandt et al. 2021).

For V Hya, the primary mass was set as a Gaussian prior of 1.9 ±1 M_⊙ and the Gaia DR3 parallax was used as prior, corresponding to a distance of 434± 21 pc. As the parallax uncertainties for AGB stars is often underestimated, different prior values were tested, corresponding to the recomputed distances from Andriantsaralaza et al. (2022): ${484}_{-177}^{+107}$ pc and ${529}_{-131}^{+248}$ pc. However, no significant difference was found in the final orbital parameters and, therefore, the Gaia DR3 parallax has been adopted in the remainder of the paper. We used 15 MCMC temperatures and for each temperature, with 100 walkers run with 5 x10⁶ steps per walker. After a first try, the MCMC was oscillating between two orbital solutions: one with a high eccentricity for a period of 35 yr and a mass ratio above 4 and another with a low eccentricity, a period of 17 yr, and a smaller mass ratio. Such a high mass ratio would imply a companion that is more massive than V Hya initial mass, which is inconsistent with its evolutionary stage. We therefore re-started the MCMC with an additional prior on the eccentricity (e < 0.15) to favour the orbital solution with the shorter period, compatible with the observed modulation in the light curve.

3.2 Results

The orbital parameters found with the two methods (step 1 and step 2) are listed in Table 1. Despite their different implementation and assumptions, both methods produce a low-eccentricity orbit (compatible with a circular nature) with a period consistent with the light-curve modulation. The ORVARA solution has larger uncertainties, but it is expected to be more robust and does not require to impose the period. Below, we discuss the values obtained and their sensitivity regarding the astrometric measurement.

3.2.1 Spectroscopic fit (step 1)

The pulsation was fitted with an asymmetric sine amplitude of 5.18 km s⁻¹. Integrating the curve over half a period to obtain the corresponding radial distance covered leads to a value of ΔR = 0.25 au. For a stellar radius of 2.5 au (see Sect. 6.1), this variation in radial distance corresponds to a variation of ± 10% of R_*. over a pulsation cycle. The mass function, f(m), obtained is equal to 0.18 M_⊙. This value gives partial indication of the companion mass with RV data only depending on the inclination angle and the mass of the primary. The inclination angle cannot be obtained with the spectroscopic signal alone. It requires additional information from astrometric measurements and an estimation of the primary mass.

3.2.2 Astrometric-spectrometric orbit (step 2)

The right panel of Table 1 gives the orbital parameters obtained when simultaneously fitting the RV curve and astrometric data. The ${\chi }_{\text{red}}^2$ is the overall reduced χ², taken as the sum of the χ² for the HIPPARCOS proper motion $({\chi }_{\text{H}}^2)$, the Gaia proper motion $({\chi }_{\text{G}}^2)$, and the long-term HIPPARCOS-Gaia proper motion $({\chi }_{\text{HG}}^2)$. The different χ² can be found in the central column of Table 2. Figure 2 displays the best fit for the RV curve and the proper motions from HIPPARCOS and Gaia. The corner plot of the MCMC chain is displayed in Fig. C.1. The best-fitting value for the inclination is 37 7 ± 2 2^°. This lies within the range of previous estimations: about 30° for Knapp et al. (1997), and 45° for Sahai et al. (2022).

3.3 Sensitivity with respect to the HIPPARCOS data

In Fig. 2, the data point corresponding to the proper motion in right ascension at HIPPARCOS epoch, ${\mu }_{\alpha *}^{\text{H}}$, appears as an outlier in the best-fitted model. The instantaneous values of the proper motion at HIPPARCOS epoch are computed by Brandt (2021), as the 60/40 weighted average between the instantaneous proper motion for the original and the second reduction (Perryman et al. 1997; van Leeuwen 2008). For V Hya, the values in right ascension between the two reductions (−14.21±1.41 masyr⁻¹ and −11.02±1.14 masyr⁻¹, resp.) are not compatible within their uncertainties (as it is for the declination component, 2.72±1.14 masyr⁻¹ and 2.29±1.11 masyr⁻¹ resp.). This discrepancy can be explained by the variability-induced mover (VIM) flag associated with V Hya for the first reduction. Using an improved chromaticity correction adapted for red stars, Pourbaix et al. (2003) did not confirm the VIM solution; whereas in the second reduction, its recomputed solution is a standard five-parameter solution that does not bear any VIM processing. It indicates the uncertain nature of the ${\mu }_{\alpha *}^{\text{H}}$ datapoint. We thus decided to keep the 60/40 mean value computed in the catalog of acceleration (Brandt 2021), but to assess the robustness of the obtained orbit regardless the ${\mu }_{\alpha *}^{\text{H}}$, a sensitivity analysis was performed. The ${\mu }_{\alpha *}^{\text{H}}$ value was artificially offset by +4, +2, −2, and −4 mas yr⁻¹ and, for each new value, ORVARA was run again to fit the RV curve and the new astrometric indicators, using the same MCMC strategy described in Sect. 3.1.2, but without any prior on the eccentricity. The results are the following: for the positive offset all orbital parameters stay unchanged, except for a change in the inclination sign. The solution with an inclination of i = 180 − 38 = 142° (and Ω = 180 + 20°) is favoured over the solution with i = 38° (and Ω = 180 − 20°). It corresponds to a µ_δ signal in anti-phase, compared to the best-fitting curve from Fig. 2 (right panel), resulting in a change in the orbit motion direction with respect to the previous solution. For the negative offset, the obtained orbital solution is bimodal: either the previous estimated solution (period of 17 yr with e < 0.1, i = 38° and a =11.2 au) or the solution with a period that is twice larger (and, thus, a mass ratio greater than 5), implying a companion mass higher than 4 M_⊙. This solution is more probable for higher shift of the µ_α*. parameter but corresponds to an important increase of the χ² (see Table 2). As discussed in Sect. 3.1.2, such a massive companion is considered as unlikely. The fit of the PM and RV curves for the four artificial ${\mu }_{\alpha *}^{\text{H}}$ are displayed in the appendix, together with their respective corner plot. We conclude that the obtained orbital parameters are robust against the unreliable ${\mu }_{\alpha *}^{\text{H}}$ value. Therefore, we took the orbital parameters (i, Ω, M₁/M₂, a) derived from the astrometry as indicative, rather than precise, and considered the orbital motion direction as undetermined: either counterclockwise rotation of the companion around the primary on the sky (i = 40^° solution) or clockwise rotation (i = 180 − 40° solution) Such a degeneracy could be lifted with additional constraints from astrometric data (forthcoming Gaia DR4 catalog) or through interferometric observations (Planquart et al. 2024).

Fig. 2 ORVARA results: RV curve (left), proper motions in right ascension (centre) and declination (right). In all plots, the black thick line represents the best-fitting orbit, while 40 other well-fitting orbits are included and colour-coded as a function of the companion mass.

3.4 Nature of the companion

The derived mass function (Table 1) allows us to constrain the nature of the companion, hereafter referred to as ‘V Hya B’. With the inclination derived at step 2, V Hya B is expected to be more massive than V Hya A. With the masses of the two components being highly correlated (see corner plot in Fig. C.1), no precise mass can be retrieved for the companion without imposing V Hya A mass. Taking a primary mass of 1 M_⊙ as lower limit, the companion mass would be of 1.9 M_⊙. Such a mass is too high to correspond to a white dwarf and therefore the companion is likely to belong to the main sequence phase. A main sequence companion with a mass above 2.0 M_⊙ would corresponds to an A-type star. The main-sequence nature of V Hya B was previously inferred from the measured UV excess by GALEX (Sahai et al. 2008). Such UV excess would require a hot main sequence star.

To test the nature of the companion obtained from our orbital analysis with the UV excess, the synthetic spectra of the two stars were superimposed on the system photometry from UV to near-IR (NIR). The V Hya B spectrum was modelled as an A0 star with a temperature of 9950 K from Kurucz models (Castelli & Kurucz 2003) and a radius 1.5 R_⊙ corresponding to a mass of 2 M_⊙. The AGB star was modelled by a black-body of 2700 K (inferred temperature from Knapp et al. 1999), as no synthetic spectra for carbon stars exist for λ < 300 nm. The SED is displayed in Fig. 3. To account for the photometric variability induced by the stellar pulsation, the uncertainties of the photometric data points have been set to 2 mag in the optical and 1 mag in the near-infrared. The spectra are reddened by combining the effect of two sources: the interstellar extinction and the extinction due to the circumstellar envelope (CSE) of the AGB star.

The visual interstellar extinction was estimated to be A_V = 0.09 mag (Lallement et al. 2019), and the extinction curve of Gordon et al. (2009) was adopted. The interstellar contribution is of negligible effect compared to the effect of the circumstel-lar dust. The effect of the dust extinction is obtained assuming a CSE composition made of amorphous carbon, AmC (Rouleau & Martin 1991). The dust grains are treated as spherical, with their radius a distributed following the standard MNR distribution (dn/da = a^−3.5, with a_min = 0.005 µm and a_max = 0.25 µm; Mathis et al. 1977). The total extinction curve is displayed in the inset of Fig. 3, and the total A_V is 3.90 mag. The additional 3.54 mag corresponds to a column density for AmC dust grains of the order of 10¹² cm⁻².

The cold nature of the primary does not make it possible to reproduce the observed flux in the UV, whose expected flux in the far-UV (FUV) filter lies several orders of magnitude below its measured value. On the contrary the A0 spectrum, dominating the overall flux for wavelengths below 300 nm, reproduces perfectly the GALEX photometric data in the NUV and FUV filters. The SED confirms the assumption of Sahai et al. (2008) who attributes the UV excess of V Hya to the presence of a hot companion. It has to be stressed that the black body radiation curve overestimates the flux at short wavelengths for red stars such as carbon stars, and explains the observed discrepancy between the model and the bluest optical point (Johnson B) in Fig. 3. Under these circumstances, we could expect to detect signatures of the hot companion even up to 450 nm, as suggested by Lloyd Evans (1991).

4 Comparison with the light-curve modulations

In this section, we compare the RV periodic motion with the light-curve variation.

4.1 The pulsation cycle

To highlight the RV signal associated to the pulsation, the contribution from the Keplerian orbit was removed. For the light-curve variation, the 17 yr signal was fitted by taking a 530-day rolling mean, thus, the signal associated to the Mira variation can be obtained by removing the obtained rolling mean.

In Fig. 4, the RV and light curve signals cleaned from the long variation are plotted as a function of the pulsation phase. The two signals are found to be in quadrature, with the RV signal being in advance with respect to the AAVSO light curve $(\text{Δ}\varphi =-\frac{\pi }{2})$. Such a quadrature delay is expected and reflects the change of the AGB photospheric radius: as the star undergoes expansion of its envelope, its brightness decreases (descending part of the light curve). The expansion leads to blue-shifted RVs for the observer (lower lobe of the RV signal). The same phase delay was already found for the CS Mira star R CMi (see Fig. 5 from Jorissen 2004).

Fig. 3 Spectral energy distribution of the V Hya system from UV to NIR (colored square symbols) and model spectra (black: V Hya as a black body, blue: companion as an A0-type star). The diamonds represent the expected FUV and NUV fluxes from the AGB star (black) and the hot companion (blue). The inset displays the adopted extinction curve from the dust (interstellar + circumstellar).

4.2 The 17-yr cycle

In Fig. 5, the RV solution (extrapolated to early cycles) is over-plotted over more than a century of photometric observations, covering six dimming events. Their duration and depth vary from cycle to cycle, with a mean amplitude of 3 mag and a mean duration of 30% of the cycle. The shape of the dimming events is symmetric: the ingress and egress have the same slope, with a value that varies from 0.001 to 0.003 mag day⁻¹, depending on the cycle.

The two long-period signals are again found to be in quadrature. But, in contrast with the short-period signal associated to the pulsation, the RV signal is phase-delayed $(\text{Δ}\varphi =-\frac{\pi }{2})$ with respect to the light curve, hinting at a different mechanism than intrinsic variations to explained this phase lag. This time, it can easily be explained in the framework of binary motion: the minimum light found by the AAVSO monitoring always corresponds to the superior conjunction, when the primary, namely, V Hya A, is behind its companion and its RVs switch from red-shifted to blue-shifted. A direct consequence is that the obscuration causing the deep minima in the light curve occurs during the companion passage in front of V Hya.

Due to the particular shape of the eclipse, with a long duration and a deep obscuration, and the system inclination, a simple eclipsing scenario by a stellar body has to be rejected. Instead, a more extended and opaque object attached to the companion is needed. The inclination angle of about 35–40^° (see Table 1) implies that the eclipsing object must have an important extension in the direction perpendicular to the orbital plane. A possible explanation could be that the primary star is eclipsed by a vertical dense outflow.

Fig. 4 RV curve cleaned from the Keplerian motion, the AAVSO light curve with the 530 days rolling-mean signal removed, and their best-fitting functions phase-folded with respect to the pulsation period. The red dots correspond to the RV data affected by shocks and not used in the analysis (see Sect. 2).

5 Time series analysis of spectral lines

As discussed in the introduction, bipolar outflows have been observed and studied in post-AGB binaries and their signature can be clearly detected in optical spectra. Bollen et al. (2022) showed that the Hα line in the spectra of post-AGB binaries with jets shows a variable profile whose time-series can be modeled to obtain different jet parameters. In this section, we report on the temporal variation of particular spectral lines: the sodium and potassium doublets, Hα, and the C₂ (0,0) vibra-tional band, and we connect them with the stellar and orbital star periodicity.

5.1 Sodium doublet

The most interesting features are found in the sodium doublet Na I D-lines. In Fig. 6, the Na I D-profiles at two different orbital phases are plotted, together with a reference spectrum. The phase 0.5 corresponds to the light minimum (superior conjunction, V Hya behind its companion) and the phase 0 corresponds to the inferior conjunction. A strong and variable blue-shifted absorption is present. This absorption is visible in both lines of the doublet and is most blue-shifted around phases 0 and less blue-shifted around phases 0.5. The profiles of the two sodium lines do not seem to differ.

To highlight the variation of the Na I D-line profile, the spectra are plotted as a function of the orbital phase and interpolated between the phases to obtain a smooth representation of the variation. This is shown in Fig. 7. A sample of spectra at different orbital phases can be found in Fig. D.1. Since the HERMES monitoring duration is shorter than the orbital phase, the phases between 0.05 and 0.35 are not yet covered and are left blank in the figures. The 17 yr periodicity of the Na I modulation was confirmed by prior spectroscopic monitoring of Lloyd Evans & Carter (1993). Indeed, a weakening of Nai absorption followed by an increasing of its emission core is reported at the beginning of the fading event of 1993 (see red rectangle in Fig. 5).

The large absorption band seen in Fig. 6 is also visible in Fig. 7 and reaches its maximum blue-shifted value, −150 km s⁻¹ at phase 0.5, and its minimum at phase 0. At phase 0.5, an emission (represented in red) is also present on the right side of the blue-shifted absorption. The combination of blue-shifted absorption and emission is reminiscent of P-Cygni profiles. The periodicity of the absorption pattern seems to be the same as the orbital one but it is phase-shifted. More precisely, they are in quadrature as the absorption band is at its maximum blue-shifted value at maximum light when the star velocity equals the centre-of-mass velocity. This is easily seen in Fig. 7, thanks to the sine black line representing the orbital motion. This phase shift is quite puzzling since one would expect it to be in phase if the absorption was due to the primary star or in anti-phase if it was due to the companion motion. Thus this feature cannot be related to the motion of the stars but must have a different origin. A possible explanation is the presence in the binary system of a conical outflow with a velocity gradient. As the system goes along its orbital motion, the cone obstructs V Hya light with a varying angle and therefore a varying absorption velocity.

We can see that even though the shapes are different (Na I D2 exhibits a double-peak emission whereas Na I D1 only a left-centred peak), their equivalent widths are quite similar. The ratio between D2 and D1 equivalent widths range from 0.8 to 1, for a spectral window of 2 Å. Those values being well below 2 (the expected values for the case of an optically-thin media, and ratio of the oscillator strength), indicating that the media is therefore optically thick.

Fig. 5 AAVSO lightcurve of V Hya recorded since 1913, displaying six evenly spaced dimming events. The blue curve is the extrapolated spectroscopic orbit. The red regions are the observing epochs of sodium absorption (see Sect. 5.1) from the HERMES monitoring (filled rectangle), and from the observing dates of Lloyd Evans & Carter (1993, dashed rectangle).

Fig. 6 Na I doublet near superior conjunction (phase 0.46, top) and near inferior conjunction (phase 0.98, bottom). The dashed line is the spectrum of the reference carbon star X Tra used as template for the modelling. The vertical dashed lines indicate the centres of the two components of the doublet.

Fig. 7 Dynamic spectra of the sodium doublet (top row), potassium doublet (middle row), C2(0,0, bottom left), and Hα (bottom right) of V Hya as a function of the orbital phase. In each panel, the observed phases are shown twice to guide the eyes, time runs down. The colour represents the pseudo-continuum-normalised fluxes taken as the median value of the spectral window. The black curve represents the primary motion while the vertical line marks the systemic velocity, and the horizontal lines on the left indicate the observed phases.

5.2 Potassium doublet

Peculiarities are also observed in another alkaline metal: the potassium doublet K I D-lines at λ = 7664.89 Å and 7698.96 Å.

A large blue-shifted absorption band is present and shares the same phase dependency as the sodium doublet, being most blue-shifted at phase 0.5. The absorption is narrower, especially at phase 0 where the bandwidth reaches about 20 km s⁻¹. Pseudoemission lines following the orbital variation are visible in the red-shifted part of the spectrum (in red in Fig. 7). Their evolution follows the orbital trend, indicating that those lines are of photospheric origin, contrary to the alkaline resonance lines whose variations are phase-shifted, hinting at a different origin above the photosphere.

5.3 Swan band

The Swan bands are vibrational bands of the C₂ molecule. Their emission is temperature-dependent and begins to develop above 2450 K (King 1948). The most intense one is the (0,0) transition at 5165 Å. As V Hya is a cool carbon star, only emission in the (0,0) band is observed in its spectrum. In Fig. 7, the C₂ profile seems to change with the orbital phase. The line profile around λ = 5165 Å is quite different between superior and inferior conjunctions. Around phase 0.5 a pseudo-emission seems to be present but this feature is absent around phase 0. Such an abrupt emission of C2 was also reported during the obscuration event of 1993 (Lloyd Evans & Carter 1993; Lloyd Evans 1997) confirming the relation between the C₂ emission and the star drop of brightness.

5.4 Balmer lines

No amplitude modulation correlated with the orbital motion is found in the Balmer lines. Only the Hα emission is strong enough to be observed with a high signal-to-noise ratio (bottom right panel of Fig. 7) and mainly shows radial-velocity variation that follows the orbital trend with a faint amplitude modulation due to the Mira pulsation. No blue-shifted absorption as for the Na I or K I lines is seen. Compared to other carbon stars with a Mira pulsation, the intensity of Hα for the star V Hya is quite reduced overall, even outside the eclipse phases, reaching at most five times the continuum. As previously reported by Lloyd Evans (1991), the Hγ line is found in absorption and could be associated to the hot companion signature. Indeed, no absorption is found in the Balmer lines of the reference carbon star X Tra.

6 Spatio-kinematic modelling of the outflow

As discussed in Sect. 5.1, the presence of varying blue-shifted absorption in the sodium doublet can be interpreted as the signature of a high-velocity outflow attached to the companion star. Bollen et al. (2021) studied the jet imprints on Balmer lines for a sample of post-AGB binaries for a variety of jet configurations. Their modelling approach consists in fitting the model spectral lines created by a parametric jet model to the phase-resolved Hα profile using an MCMC routine.

To deduce the jet structure (geometry, velocity, and density profile) we used the same spatio-kinematic modelling engine² and adapted it for the phase-resolved sodium line at 5895.92 Å. Hereafter, we summarise the overall fitting routine adopted to fit V Hya sodium spectral line over the range 0 to −250 km s⁻¹. We refer to the reference paper (Bollen et al. 2019) for a more general description of spatio-kinematic jet structures and their numerical implementation. Among the jet models available, we restrict ourselves to the stellar jet model, which is similar to the model developed by Thomas et al. (2011) for the Red Rectangle Nebula. It has been shown that the fitting routine provide similar results for different configurations, hence we choose the most generic model containing the smallest number of free parameters and assumptions about the launching site of the outflow. The stellar jet model and its different parameters are described in Appendix E.1. As the absorption feature (believed to be the jet imprint) is present at all observed phases which cover nearly 70% of the orbit, the jet opening angle is expected to be larger than the inclination angle of the system. Therefore, no observed spectrum can be used as template for the unattenuated (background) spectrum in the fitting process. Instead, we used the high-resolution spectrum of the carbon star X Tra as template for the unattenuated spectrum. The spectrum was modulated in RV to reproduce V Hya orbital motion and in flux at to reproduce the 530 days intrinsic variability. The dynamic spectrum is shown in Appendix E.2.

The procedure is meant to quantitatively model the absorption of the stellar contribution through the jet. However, in the case of V Hya, an extra emission pattern is seen in the sodium profile and needs to be modelled and implemented on top of the Bollen et al. (2019) method. The modulated emission pattern is not found in other carbon-star spectra, suggesting it to be the result of the star-jet interaction rather than of the star photosphere. We interpreted the blue-shifted emission as rescattering of V Hya light in our line of sight. Additional details on the implementation and normalisation of those two effects are described in Appendix E.3. In summary, the general fitting procedure is a three-step process: first, the synthetic model spectrum is used to fit the stellar-jet model to the phase-dependent observations and to obtain the best-fitting absorbing jet. Then, from the obtained spatio-kinematic parameters, the additional source term is modelled. Lastly, the jet parameters are fine-tuned by fitting the absorbing jet on an updated spectrum consisting in the synthetic model and the modulated source term.

6.1 Best-fit model

The jet parameters are listed in Table 3, together with their best-fitting value and their initial range. In the table, i represents the inclination of the orbital plane, α is the jet opening angle, θ_cav is the angle of the jet cavity, θ_tilt is the angle between the jet axis and the normal to the orbital plane, υ₀ is the gas velocity at the jet axis, υ_α is the velocity at the cone edge, p_v and p_d are the power indices governing the density and velocity laws (see Appendix E.1.), c_τ is the scaling factor (see Appendix E.1.), and R₁ is the radius of the primary star. Due to a strong correlation between the maximal jet opening angle, α, and the system’s minimal inclination, we constrained the inclination angle, i, to its value obtained from the orbital analysis (see Table 1). For the definitions of all the other parameters, we refer to Appendix E.1. The velocity and density structure across the jet is displayed in Fig. 8. The velocity follows a decreasing and nearly-linear trend starting from υ₀ at the jet axis to υ_α at the jet edge. The maximum deprojected velocity, υ₀, is below 200 km s⁻¹, although such a value is well above the velocity of the wind feeding the AGB circumstellar envelope.

The normalised density profile is maximum at the edge of the cone and minimum along its axis, with a cavity angle θ_cav. Such a hollow cone is also encountered for post-AGB binaries (Bollen et al. 2019) and reveals the presence of a cavity along the jet axis. The presence of a cavity in the V Hya system was already put forward to explain the launching mechanism of the [S II] jets observed by Sahai et al. (2016) and Scibelli et al. (2019). The observed and modelled dynamic spectra of the Na I D1 line is displayed in Fig. 9. It can be seen that the model reproduces the temporal modulation of the blue-shifted absorption quite well. The obtained jet geometry and orbital configuration are sketched in Fig. 10.

6.2 Limitation of the model

The above-described analysis allows us to conclude that a high-velocity outflow launched by a companion orbiting V Hya is able to satisfactorily reproduce the observed spectral modulation. The choice of such a modelling tool was motivated by previous large-scale observations of high-velocity outflows (Knapp et al. 1997; Hirano et al. 2004) and by the close proximity of V Hya evolutionary stage to the post-AGB phase where those structures are widely found. However, no direct evidence of the companion motion is seen in the spectrum as its flux contribution is expected to negligible above 450 nm (see Sect. 3.4), and the exact launch site of the jet remains unknown. Further analyses are needed to constrain the temperature and the absolute density structure through a full 3D radiative transfer code. Those ingredients and the putative properties of the companion (see Sect. 3.4), are needed to fully characterise the physical mechanism feeding and launching the jet.

Fig. 8 Best-fit model of the velocity and density laws as a function of the jet angle. The density at the edge is set to 1.

Fig. 9 Observed dynamic spectrum of the Na I D1 line (left) and the model spectrum obtained by the spatio-kinematic fitting of a conical jet (right).

7 Discussion

In this section, we discuss the possible origin of the long-period obscuration event in the framework of the binary model obtained. Alternative explanations not invoking a companion are then discussed. Finally, we confront the obtained model of the V Hya system with previous large-scale observations to build a multi-scale consistent picture of the system.

7.1 Obscuration by dust

As highlighted by the orbital analysis in Sect. 4, the dimming episodes of about 3 mag occurring every 17 yr are associated with the superior conjunction and the high-velocity absorption of alkaline resonance lines. These observations can be explained in the framework of an absorbing gaseous jet passing in front of the line of sight. However, a gaseous jet alone cannot be responsible for the visual dimming event, as such a broad-band attenuation requires particles with a gray absorbing power, such as dust grains. A light curve modulation induced by the variation of the dust grains column density through the line of sight with the orbital phase would require a non-uniform distribution of particles above the orbital plane, due to the inclination well below 90 degrees. For a mean magnitude drop of 3 mag, the observed modulation can only be recovered by assuming the highest dust density at the jet centre, with a dust column density that is three times denser at the superior conjunction when the line of sight crosses the axis of the cone. The dust distribution through the jet polar angle would be opposite to what was found for the gas density in the jet, better described by a hollow cone at its centre. It raises the question of the physical mechanism behind the dust and gas distributions. Only a better characterisation of the absolute gas density and temperature using a full 3D radiative transfer code would allow us to determine more precisely the dust formation limit, the mass-to-dust ratio, and the exact dust spatial distribution in the system and understand their possible coupling with the gas. However, such a complex investigation is out of the scope of the current paper, as we are focused on the system dynamics. Interestingly, some post-AGB systems also exhibit a long-term modulation of their light curve, but of a smaller amplitude, (<0.5 mag) associated with the radial-velocity and atomic-line variations. In RV-Tauri stars, those obscurations occur at the inferior conjunction, the opposite situation as encountered here, explained by the particular orbital configuration of the system where the line of sight grazes the circumbinary disc due to the high orbital inclination (Manick et al. 2019). A few other post-AGB stars with a bipolar morphology show small-amplitude obscuration events together with radial-velocity and Hα variation: For instance, V501 Lup exhibits obscurations of about 0.5 mag in the visible and near-infrared which share the orbital periodicity (2727±26 days) and could be the result of dust absorption and scattering by the nebula (Manick et al. 2021). For 89 Her, a long-term obscuration with a periodicity of 5000 d is reported together with modulation of Hα, Na I D, and Ba II lines. With the long period of the light curve being larger than the orbital one (289 days, Bollen et al. 2021), two scenarios have been suggested: an orbital wobbling of the circumbinary disc inducing a misalignement with the hour-glass structure or recurrent mass-ejection episodes (Gangi et al. 2021).

Fig. 10 Geometric representation to scale of the binary system where V Hya A is represented as an orange filled circle and V Hya B as a black dot (to be visible, its radius was increased). V Hya A mass was set to 2 M⊙, corresponding to a mass of 2.5 M⊙ for V Hya B. Upper left panel: system depicted edge-on at the superior conjunction. The black line is oriented along the line of sight. The green cross indicates the centre of mass. Lower left panel: system viewed face-on, with the Roche equipotential in black line and the Roche lobe of the two stars in dashed gray. Right panel: system projected on the sky-plane at the superior conjunction for the counterclockwise solution, assuming a parallax of 2.31 mas. The black line is the best projected orbit, the dashed black line represents the line of nodes. The effect of the mass-ratio uncertainty on the orbital separation is illustrated by the colored ellipse.

7.2 Origin of the long-period variability

For the sake of completeness, it ought to be mentioned that alternative scenarios not involving a binary system can be invoked to explain the obscuration of variable stars: a long-period secondary pulsation period that is twelve times longer than the fundamental tone of 530 days or an episodic dust condensation from the star. Below, we describe both of them and argue why we consider them less likely. For the first one, as discussed in Sect. 4, the phase shift between the radial velocity and light curves is found to be +90° for the long cycle. Therefore, the usual scenario of radial pulsation induced by the so-called kappa-mechanism (associated instead with a −90° phase shift) cannot be responsible for the 17 yr periodic fading. Besides, this long period would make V Hya an outlier in the usual period-luminosity diagram, presented by Wood (2000), being too long to belong to the D-sequence (“long-secondary periods”).

The second interpretation, namely, episodic dust formation, is the usual mechanism put forward to account for the abrupt fading by several magnitudes of R Coronae Borealis (R CrB) variables. Those stars exhibit light curves with important dimming events separated by several years (Clayton 1996). As these obscurations are correlated with intrinsic dust formation, their occurrence is far from periodic and their fading episodes show an abrupt decrease followed by a gradual restoration of the star brightness. The only exception is the sub-class of DY-Per objects, labeled as cold hydrogen-deficient R CrB stars, whose light curves exhibit a rather regular periodicity in their fading. Considering the regularity of V Hya light curve over a century of observations and previous hints in favour of the binary hypothesis (Hirano et al. 2004; Sahai et al. 2008, 2016) we consider intrinsic variability scenarios to be less likely to explain the long-term periodicity of the star.

Fig. 11 Sky-projected system orientation compared to radio observation. The contour levels is a stereogram view of the 12CO J = 2−1 high-velocity emission. The contour levels are drawn at 6, 8, and 10 mJy beam−1 for two maps equally shifted by ±100 km s−1 from the centre-of-mass velocity. The blue-shifted contours are displayed in blue, and the red-shifted in red. The dashed green line represents the orbital line of node, the black arrow the jet axis and the dashed black line sketches the cone edges.

7.3 Connection to large-scale bipolar outflows

The characterisation of the orbital parameters (including the system inclination) associated with the detection of the high-velocity jet allows us to link the different building blocks of the system; in particular, the jet that is supposed to be attached to the companion, with its sky-projected orientation. A sketch of the system is displayed in Fig. 10. As discussed, the inclination found is consistent with the previous inclination estimates from the bipolar outflow measurement. Additionally, the ascending node being close to 180° (see Table 1), the jet axis (perpendicular to the orbital plane) points towards the east. Its direction is consistent with previous measurements, at a larger scale, of the high-velocity blue-shifted [SII] emission and the high-velocity parabolic outflow found in the radio-emission of CO molecules (Sahai et al. 2016, 2022). The line of nodes and jet axis angle have an uncertainty of ±10° due to the inclination degeneracy (see Sect. 3.1.2).

Figure 11 superimposes the sky-projected jet axis and the system’s line of node on the high-velocity emission obtained from ALMA observations (PI: Sahai R., Project ID: 2018.1.01113.S). Both geometry and orientation between the two structures match. This result is to be underlined as the two sets of observations, having entirely independent approaches, provide a complementary set of evidence for the system morphology. On the one hand, in our study, the conical geometry and orientation of the jet is deduced by fitting phase-dependent blue-shifted absorption of optical resonance lines and compare them with the system orbit. On the other hand, the bipolar jet is mapped using a one-shot observation in the radio-emission of ¹²CO J = 2−1 at larger scale.

This compatibility strongly suggests that they may share a common origin: the modeled conical jet could be seen as the innermost part of a large-scale parabolic structure. This scenario does demonstrate the link between the dynamical binary interaction (as inferred from our spectral monitoring) and the morphology of the large-scale bipolar nebula.

8 Conclusion

Thanks to an unprecedented 12-yr-long spectroscopic monitoring of V Hya, we performed a temporal study of the dynamical processes at the origin of its two periodic variations. The main results obtained in this paper are two-fold.

In the first step, the radial-velocity curve revealed a short variation of 530 days superimposed on a long trend. Combining the RV data with secular visual photometry and astrometric acceleration, we disentangled the signal associated with the Mira-like pulsation of 530 days from the 17 yr cycle. The latter is interpreted as the signature of a Keplerian orbit for a period equal to the long photometric cycle of 6351 days, a low eccentricity (e = 0.02 ± 0.02), and an inclination of 37°±2°, compatible with a main sequence companion.

In a second step, the time-resolved spectroscopy unveiled the presence of a high-velocity gas stream in the system. The resonant lines of sodium and potassium display phase-dependent absorption, blue-shifted up to −150 km s⁻¹ at the epoch corresponding to the visual minimum, and superior conjunction. We interpret the gaseous jet as the imprint of an interacting binary system. A spatio-kinematic modeling of the jet launched from the companion orbiting around V Hya is fitted to reproduce the observed spectral-line modulation. The obtained jet model is described as a hollow cone with a wide opening angle and a velocity gradient along its radial direction. The highest velocities are oriented along the jet axis, whose projection on the sky plane is pointing toward the east. In the meantime, the jet geometry obtained (velocity and direction) is consistent with previous large-scale observations of the high-velocity bullets found in [SII] lines by Sahai et al. (2016), and high-velocity bipolar flow from the molecular lines of CO (Knapp et al. 1997; Hirano et al. 2004; Sahai et al. 2022).

It is, to our knowledge, the first case where the temporal variation of the light curve, high-velocity gas, and the AGB star velocity are reported to be correlated. Our proposed scenario – similar to what is found in many post-AGB binaries (e.g. Bollen et al. 2022 and references therein), where the evolved star orbits around the jet launched by its main sequence companion – explains the reported observations in a common framework. To conclude, this multi-epoch study reveals the complex behavior of V Hya, which can be seen as a benchmark case to better understand the shaping mechanism induced by a binary system at the origin of bipolar planetary nebulae.

Further studies reproducing the radiative transfer inside the circumstellar and circumbinary environments of V Hya system could help to constrain the connection between the binary interaction and the large-scale bipolar morphology. Additionally, enlarging the presented analysis to a larger sample of AGB stars with known companion or bipolar morphology, for instance, TU Tau (Gray et al. 2023), R Crt (Khouri et al. 2020), UV Aur (Herbig 2009), RR UMi (Dettmar & Gieseking 1983) would be needed to obtain a statistical distribution of the orbital parameters of such systems.

Acknowledgements

L.P. is FNRS research fellow. A.J. acknowledges support from the Fonds de la Recherche Fondamentale Collective (FNRS, F.R.F.C.) of Belgium through grant PDR T.0115.23, and from FNRS-F.R.S. through grant CDR J.0100.21 (‘ASTRO-HERMES’). Use was made of the Simbad database, operated at the CDS, Strasbourg, France, and of NASA’s Astrophysics Data System Bibliographic Services. This research made use of Numpy, Matplotlib, SciPy and Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2018). This work presents results from the European Space Agency (ESA) space mission Gαiα. Gαiα data are processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia MultiLateral Agreement (MLA). The Gaia mission website is https://www.cosmos.esa.int/gaia. The Gaia archive website is https://archives.esac.esa.int/gaia. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2018.1.01113.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. We acknowledge with thanks the variable star observations from the AAVSO International Database contributed by observers worldwide and used in this research.

Appendix A Radial velocities

Tables A.1 lists the radial-velocities from the HERMES spectrograph used in Sect. 3.

Appendix B Light-curve periodogram

Figure B.1 displays the Lomb-Scargle periodogram of the AAVSO light curve and its four prominent peaks (Scargle 1982). The two left-most peaks correspond to the observed 6319 d and 531 d periods with a full width at half maximum (FWHM) of 940 days and 6 days, respectively. The two right-most peaks are aliases located at beating frequencies between the 530 day signal, f₅₃₀, and the one-year signal (f₃₆₅, red continuous line): f_alias = f₃₆₅ ± f₅₃₀.

Fig. B.1 Lomb-Scargle periodogram of the AAVSO light curve

Appendix C orvara results

Appendix C.1 Convergence study

Figure C.1 displays the corner plot from ORVARA. A strong correlation is present between the mass of the two stars (M_pri, M_sec) and, consequently, their separation, a. The inclination posterior distribution is slightly bimodal, with a main peak located at i = 37°±2°. The secondary faint peak corresponds to 180° − 37°±2°. The convention is taken such as i < 90° corresponds to counterclockwise rotation of the secondary around the primary while i > 90° is the opposite situation.

The convergence of the chains was assessed by comparing the variance between chains with the variance within each single chain through the potential scale reduction factor (PSRF) estimator (Gelman & Rubin 1992). The PSRF was < 1.1 for each parameter.

Appendix C.2 Sensitivity study on the ${\mu }_{\alpha \ast }^H$ point

Figures C.2, C.3, C.4, and C.5 display the proper motion (PM) and radial-velocity (RV) curves, together with the corner plot obtained for four ${\mu }_{\alpha \ast }^H$ values, offset from the original value of the catalog of acceleration (Brandt 2021), by +4, +2, −2, and −4 mas/yr, respectively.

Fig. C.1 Corner plot of some derived orbital parameters: the mass of the two stars, semi-major axis, eccentricity, and system inclination.

Fig. C.2 PM curve in right-ascension and declination (top panels), RV curve (middle panel), and corner plot (bottom panel) for ${\mu }_{\alpha \ast }^H$ offset by +4 mas/yr.

Fig. C.3 Same as Fig. C.2 but for ${\mu }_{\alpha \ast }^H$ offset by +2 mas/yr.

Fig. C.4 Same as Fig. C.2 but for ${\mu }_{\alpha \ast }^H$ offset by −2 mas/yr.

Fig. C.5 Same as Fig. C.2 but for ${\mu }_{\alpha \ast }^H$ offset by −4 mas/yr.

Appendix D Representative line profiles

Figures D.1 and D.2 show the temporal behaviour of spectral lines presented: Hα, Na I, K I, and C₂ (0,0). The latter, being in the bluest part of the spectrum with a low S/N, has been smoothed by a Gaussian filter of FWHM = 2.8 km s⁻¹.

Fig. D.1 Sodium doublet (leftmost and middle panels) and Hα (rightmost panel) lines as a function of the orbital phase. The vertical dashed line marks the systemic velocity of −7.94 km s−1. The corresponding phase and observing date (expressed in JD-2400000) are displayed on the right-most axis.

Fig. D.2 Same as Fig. D.1 but for potassium doublet lines (left-most and middle panels) and the C2 (0,0) line (right-most panel).

Appendix E Spatio-kinematic modelling

Appendix E.1 Jet configuration

In this section we recall the stellar jet configuration used in this study, as introduced in Bollen et al. (2020).

Fig. E.1 Schematic geometry of the stellar jet model and associated parameters. Each point in the jet is described by a coordinate pair: (θ,z), where θ is the angle from the jet axis and z is the vertical distance from the orbital plane.

The velocity of the matter inside the jet is described by a latitudinal dependent power law:

$$\upsilon (\theta )=\left[{\upsilon }_0+\left({\upsilon }_{\alpha }-{\upsilon }_0\right)f_1(\theta )\right]1_r,$$(E.1)

with

$$f_1(\theta )=\frac{e^{-|\theta /\alpha |p_{\upsilon }}-e^{-p_{{}_{\upsilon }}}}{1-e^{-p_{\upsilon }}}$$(E.2)

Here, 1_r is the radial vector, θ the polar latitudinal coordinate, α the half-opening angle of the jet, p_υ is a power index to fit, υ₀ is the velocity at the centre, and υ_α the velocity at the cone edge.

The density profile is also latitude-dependent and decreases with the jet height z as follows:

$$\rho (\theta ,z)={\left(\frac{\theta }{\alpha }\right)}^{p_d}{z}^{-2},$$(E.3)

where p_d is a power index to fit and the z dependence is fixed to satisfy the continuity equation along the cone height.

The stellar jet model is composed of 10 free parameters: the inclination of the system, i, which we constrain by the astrometric results: the jet opening angle, α; the velocity at the edge, υ_α, and at the centre, υ₀; the power index of the velocity law, p_υ; the power index of the density, p_d; θ_cav, the angle of the cavity within the jet; θ_tilt is the angle misalignment with respect to the orbital axis (z coordinate); R₁ is the radius of the emitting surface; and c_υ, a general scaling factor for the optical depth.

Appendix E.2 Dynamic model spectrum

Archive spectra of carbon stars from the UVES Paranal Observatory Project (POP) database (Bagnulo et al. 2003) have been compared to select a template spectrum for V Hya; among them HD134453 (X Tra), HD92055 (U Hya), and HD32088. The three tested stars show a similar behavior around the sodium doublet lines, confirming the unusual nature of V Hya spectral profile. We kept X Tra as a template model. In Fig. 6, the spectrum of V Hya is superimposed to the one of X Tra for two opposite orbital phases.

Fig. E.2 Dynamic model spectrum of X Tra for the sodium doublet at 5895.92 Å (left). Dynamic model spectrum with an additional 530 d-modulation (right).

To obtain the dynamic synthetic spectrum at 5895.92 Å, the observed X Tra spectrum was Doppler-shifted according to V Hya orbital motion (see Fig. E.2, left). A 530-d intensity modulation is added to mimic the effect of the intrinsic scatter due to the Mira-like pulsation. The final dynamic spectrum used as template for the fitting routine is displayed in the right panel of Fig. E.2.

Appendix E.3 Radiative transfer through the jet

At phase ϕ, the total intensity at a specific wavelength, I_λ. is obtained by summing the contribution of two terms:

$$I_{\lambda }(s)=I_{\lambda }^0{e}^{-\tau (s)}+\beta S_{\lambda }^{scat}.$$(E.4)

The first term describes the attenuation of the photospheric light due to the absorption of photons along the line of sight by the neutral atoms present the jet. The region attenuating the incoming stellar light is delimited by a cylinder normal to the effective stellar surface. The situation is illustrated in Fig. E.3. $I_{\lambda }^0$ represents the unattenuated initial intensity and τ(s) is the optical depth at a specific position s along the line of sight in the jet. This term is proportional to the density at that point τ(s) ∝ c_υρ(s)Δs, with Δs the distance between two consecutive points on the line of sight and c_υ a scaling parameter to fit. We refer to the reference paper (Bollen et al. 2019) for additional details on the numerical implementation of the radiative transfer equation and the spatio-kinematic modelling.

The second term is a scattering source. Photons emitted from the star in a direction different from the line of sight may also intersect the cone and interact with gaseous media. The interaction is modelled by a change in the direction of propagation. A fraction of the photons is bounced back in the line of sight and contributes to the scattering source of the gaseous cone – $S_{\lambda }^{scat}$ (illustrated with the red arrow in Fig. E.3). Knowing the velocity field inside the cone, the source term is modelled by performing a first-order ray-tracing of a grid of scattering points evenly spaced in the cone. The final source function is obtained by summing the contribution of each scattered photon depending on their final projected velocity and weighted by their probability to be re-scattered in the direction of the observer. This probability only depends on the deflection angle between the photon initial direction and the line of sight.

Fig. E.3 Illustration of the stellar light propagating through the cone: the light passing through the gray shaded area is attenuated by the jet cone (see the blue line tracing the ray). The light emitted in another direction but passing through the rest of the cone can be re-scattered or re-emitted into the line of sight (red arrow).

The orbital modulation, β, accounts for the normalisation correction due to the different origins of two emitting sources unresolved by the instrument. The scattering source being emitted from an extended region offset from the star, it does not follow V Hya orbital motion and is therefore not affected by the periodic obscuration event. The spectra displayed in the analysis are always normalised by averaging over a given spectral window dominated by photospheric pseudo-emission. The obtained window-normalised spectra can therefore be compared, independently of their flux modulation. The side effect of such a normalisation is that any extrinsic emission, unaffected by V Hya obscuration, is artificially modulated by a factor inversely proportional to the stellar flux modulation: appearing brighter when the star fades and nearly vanishing at the opposite phase. In practice, such extrinsic emission will vary in strength with the orbital phase, ϕ, according to a correction factor, expected to be related to the light-curve modulation, equal to:

$${\beta }_{\varphi }={10}^{0.4\text{Δ}{m}_{\varphi }},$$(E.5)

where ∆m_ϕ is the visual light-curve modulation expressed in magnitude. Figure E.4 illustrates the effect of both types of normalisation on the relative strength of the Na I D1 line. The two spectra displayed are taken near light maximum regarding the pulsation cycle but at two different orbital periods. The magnitude difference in the visual between the two epochs is equal to 1.31, corresponding to a flux attenuation factor of about 3.31.

Fig. E.4 Normalisation effect on the Na I D1 line profile for spectra taken at two different phases, (a) The normalisation is obtained by dividing each spectrum by its mean over a 100 Å window. This procedure is incorrect (see text), (b) The normalisation between the two spectra takes into account the relative flux ratio of 0.33 between the two phases corresponding to a drop of 1.31 mag. The RV zero point is set to λ = 5895.92 Å.

References

All Tables

All Figures

	Fig. 1 Radial velocities of V Hya. Top: RV curve as a function of time. The full line is the Keplerian orbit and the dashed line is the composite signal (orbit and pulsation). Middle: phase-folded cleaned orbit of V Hya. Bottom: residual values with the ±3σ from the mean is represented as the grey-shaded region.
In the text

	Fig. 2 ORVARA results: RV curve (left), proper motions in right ascension (centre) and declination (right). In all plots, the black thick line represents the best-fitting orbit, while 40 other well-fitting orbits are included and colour-coded as a function of the companion mass.
In the text

	Fig. 3 Spectral energy distribution of the V Hya system from UV to NIR (colored square symbols) and model spectra (black: V Hya as a black body, blue: companion as an A0-type star). The diamonds represent the expected FUV and NUV fluxes from the AGB star (black) and the hot companion (blue). The inset displays the adopted extinction curve from the dust (interstellar + circumstellar).
In the text

	Fig. 4 RV curve cleaned from the Keplerian motion, the AAVSO light curve with the 530 days rolling-mean signal removed, and their best-fitting functions phase-folded with respect to the pulsation period. The red dots correspond to the RV data affected by shocks and not used in the analysis (see Sect. 2).
In the text

	Fig. 5 AAVSO lightcurve of V Hya recorded since 1913, displaying six evenly spaced dimming events. The blue curve is the extrapolated spectroscopic orbit. The red regions are the observing epochs of sodium absorption (see Sect. 5.1) from the HERMES monitoring (filled rectangle), and from the observing dates of Lloyd Evans & Carter (1993, dashed rectangle).
In the text

	Fig. 6 Na I doublet near superior conjunction (phase 0.46, top) and near inferior conjunction (phase 0.98, bottom). The dashed line is the spectrum of the reference carbon star X Tra used as template for the modelling. The vertical dashed lines indicate the centres of the two components of the doublet.
In the text

Fig. 7 Dynamic spectra of the sodium doublet (top row), potassium doublet (middle row), C2(0,0, bottom left), and Hα (bottom right) of V Hya as a function of the orbital phase. In each panel, the observed phases are shown twice to guide the eyes, time runs down. The colour represents the pseudo-continuum-normalised fluxes taken as the median value of the spectral window. The black curve represents the primary motion while the vertical line marks the systemic velocity, and the horizontal lines on the left indicate the observed phases.

In the text

	Fig. 8 Best-fit model of the velocity and density laws as a function of the jet angle. The density at the edge is set to 1.
In the text

	Fig. 9 Observed dynamic spectrum of the Na I D1 line (left) and the model spectrum obtained by the spatio-kinematic fitting of a conical jet (right).
In the text

Fig. 10 Geometric representation to scale of the binary system where V Hya A is represented as an orange filled circle and V Hya B as a black dot (to be visible, its radius was increased). V Hya A mass was set to 2 M⊙, corresponding to a mass of 2.5 M⊙ for V Hya B. Upper left panel: system depicted edge-on at the superior conjunction. The black line is oriented along the line of sight. The green cross indicates the centre of mass. Lower left panel: system viewed face-on, with the Roche equipotential in black line and the Roche lobe of the two stars in dashed gray. Right panel: system projected on the sky-plane at the superior conjunction for the counterclockwise solution, assuming a parallax of 2.31 mas. The black line is the best projected orbit, the dashed black line represents the line of nodes. The effect of the mass-ratio uncertainty on the orbital separation is illustrated by the colored ellipse.

In the text

Fig. 11 Sky-projected system orientation compared to radio observation. The contour levels is a stereogram view of the 12CO J = 2−1 high-velocity emission. The contour levels are drawn at 6, 8, and 10 mJy beam−1 for two maps equally shifted by ±100 km s−1 from the centre-of-mass velocity. The blue-shifted contours are displayed in blue, and the red-shifted in red. The dashed green line represents the orbital line of node, the black arrow the jet axis and the dashed black line sketches the cone edges.

In the text

	Fig. B.1 Lomb-Scargle periodogram of the AAVSO light curve
In the text

	Fig. C.1 Corner plot of some derived orbital parameters: the mass of the two stars, semi-major axis, eccentricity, and system inclination.
In the text

	Fig. C.2 PM curve in right-ascension and declination (top panels), RV curve (middle panel), and corner plot (bottom panel) for ${\mu }_{\alpha \ast }^H$ offset by +4 mas/yr.
In the text

	Fig. C.3 Same as Fig. C.2 but for ${\mu }_{\alpha \ast }^H$ offset by +2 mas/yr.
In the text

	Fig. C.4 Same as Fig. C.2 but for ${\mu }_{\alpha \ast }^H$ offset by −2 mas/yr.
In the text

	Fig. C.5 Same as Fig. C.2 but for ${\mu }_{\alpha \ast }^H$ offset by −4 mas/yr.
In the text

	Fig. D.1 Sodium doublet (leftmost and middle panels) and Hα (rightmost panel) lines as a function of the orbital phase. The vertical dashed line marks the systemic velocity of −7.94 km s−1. The corresponding phase and observing date (expressed in JD-2400000) are displayed on the right-most axis.
In the text

	Fig. D.2 Same as Fig. D.1 but for potassium doublet lines (left-most and middle panels) and the C2 (0,0) line (right-most panel).
In the text

	Fig. E.1 Schematic geometry of the stellar jet model and associated parameters. Each point in the jet is described by a coordinate pair: (θ,z), where θ is the angle from the jet axis and z is the vertical distance from the orbital plane.
In the text

	Fig. E.2 Dynamic model spectrum of X Tra for the sodium doublet at 5895.92 Å (left). Dynamic model spectrum with an additional 530 d-modulation (right).
In the text

	Fig. E.3 Illustration of the stellar light propagating through the cone: the light passing through the gray shaded area is attenuated by the jet cone (see the blue line tracing the ray). The light emitted in another direction but passing through the rest of the cone can be re-scattered or re-emitted into the line of sight (red arrow).
In the text

Fig. E.4 Normalisation effect on the Na I D1 line profile for spectra taken at two different phases, (a) The normalisation is obtained by dividing each spectrum by its mean over a 100 Å window. This procedure is incorrect (see text), (b) The normalisation between the two spectra takes into account the relative flux ratio of 0.33 between the two phases corresponding to a drop of 1.31 mag. The RV zero point is set to λ = 5895.92 Å.

In the text