© ESO 2021

1. Introduction

β Lyr A (HR 7106, HD 174638) is an archetype of a semi-detached binary in a rather rapid phase of mass transfer (of the order of 2 ⋅ 10⁻⁵ M_⊙ yr⁻¹) between binary components. Its orbital period has been increasing by a high rate of 19 sec per year. While during the interferometric campaign in 2013 the value of the period was 12$\stackrel{\text{d}}{.}$9427, today in 2020 it is already 12$\stackrel{\text{d}}{.}$9440. The gainer (primary) is an early B star hidden in an optically thick accretion disc. The donor (secondary) is a late B star filling its Roche lobe. For a summary of numerous investigations of this object, we refer the reader to Sahade (1966), Harmanec (2002), and Skulskii (2020).

In Mourard et al. (2018), we summarised more recent studies and attempted to model optically thick matter within the system. The model of β Lyr A was constrained by wide-band light curves (far UV to far IR, hereafter FUV and FIR, respectively) and continuum interferometric measurements. It was thus sensitive to the properties of the Roche-filling secondary, the primary, and its opaque accretion disc. The following values were adopted from previous studies: semi-amplitudes of the radial-velocity (RV) curves K₁ = 41.4 km s⁻¹, K₂ = 186.3 km s⁻¹, implying the inverse mass ratio q = 0.223; the projected semi-major axis asini = 58.19 R_⊙; and the masses M₂ = 13.048 M_⊙, M₁ = 2.910 M_⊙. The model led to an estimated distance to the system of ∼320 pc and to the finding that the accretion disc fills the whole available space of the Roche lobe in the orbital plane.

This study represents an extension of that work to optically thin parts of circumstellar matter within the system. To this end, we used additional observational data, in particular spectral-energy distribution (SED) data, high-resolution spectroscopy, and differential interferometry to measure absolute fluxes, emission line profiles, and wavelength-dependent brightness distribution at the same time. This allowed us to model the properties of the disc atmosphere, jets, or possible shell-like structures. For this we used a geometrically constrained model, described by a limited set of geometrical objects and a limited number of parameters. This method was preferred because an image reconstruction from limited spatial frequencies of interferometric measurements was not possible. On the other hand, any geometrical model uses numerous assumptions; for example, Roche geometry is used for stellar surfaces, some symmetries, or a-priori knowledge.

2. Observational data

All observational data used in our previous study (Mourard et al. 2018) remain the same, that is, the light curves and optical interferometric data. We thus refer the reader to this work for a detailed description. We just recall that the interferometric measurements give access to information on the brightness spatial distribution of the source. The squared visibilities sample the Fourier transform of the distribution at a spatial frequency defined by the baseline vector projected on the plane of sky divided by the central wavelength of the observed band, B/λ. Closure phase and triple product amplitudes are self-calibrated estimators based on the interferometric data considered on a triplet of telescopes. Hereinafter, we describe only the additional data used.

2.1. Spectral-energy distribution (SED)

To constrain the absolute flux of β Lyr A, we used data from Burnashev & Skulskii (1978). This low-resolution absolute spectrophotometry covers the wavelength range from 3300 to 7400 Å: including the Hα, Hβ, Hγ, Hδ, the Balmer jump, as well as the He I 5876 and 6678 lines. The effective bands are 25 Å wide, which is not enough to resolve the spectral lines. These are well represented by our high-resolution spectra described below. The fluxes were calibrated on Vega (α Lyr), based on its absolute calibration by Tereshchenko & Kharitonov (1972). When interpreting these fluxes, one should be more careful in the near UV (NUV) region, where the calibration is generally more difficult.

We performed a dereddening of the absolute fluxes to account for interstellar extinction. For the galactic coordinates l = 63.1876°, b = 14.7835° and the distance moduli μ = 5log₁₀[d]_pc − 5 ≐ 7.3 to 7.6, we would expect a value of at most E(B − V)=0.020 (Green et al. 2015). Using standard relations for A_V = 3.1 E(B − V) and A_λ/A_V (Schlafly & Finkbeiner 2011), we increased the observed absolute fluxes F_λ accordingly.

2.2. High-resolution spectroscopy (SPE)

We have at our disposal 72 Ondřejov CCD spectra from 2013 to 2015 (JDs 2456450.37 to 2457294.30). The spectra have a linear dispersion of 17.2 Å mm⁻¹ and a two-pixel resolution of 12700. They cover the wavelength region of approximately 6300 to 6730 Å. Their initial reductions (flat-fielding, wavelength calibrations and creation of 1D spectra) were carried out by MŠ in IRAF. Normalisation and measurement of a selection of telluric lines to be used for a fine correction of the radial-velocity (RV) zero point were carried out by PH in SPEFO (Horn et al. 1996; Krpata 2008). For modelling, we used a selection of 11 spectra, covering different phases of the orbital period.

Given the nature of the β Lyr A system, we performed a Doppler-tomography analysis. This would enable us to resolve a 3D structure and velocity fields of the circumstellar matter.

2.3. Differential interferometry (VAMP, VPHI)

In Mourard et al. (2018), we presented an extensive interferometric dataset recorded during a coordinated campaign in 2013. Data were obtained on the NPOI array (Armstrong et al. 1998) and on the CHARA array (ten Brummelaar et al. 2005) with the MIRC (Monnier et al. 2004) and VEGA (Mourard et al. 2009, 2011) instruments. This first paper was dedicated to the study of the opaque accretion disc and made use of all the interferometric measurements in the continuum bands from 525 to 861 nm and in the H band.

For the work presented here, the VEGA data were used to measure differential complex visibilities in the Hα (6562 Å) and He I (6678 and 7065 Å) lines. Differential complex visibilities were estimated through the cross-spectrum of the interferometric data between a first wide reference spectral band and a second narrow analysis band crossing the first one. Knowing the shape of the object from the squared visibilities in the reference band, the differential data made it possible to extract information on the variation of the shape of the object at high spectral and spatial resolution over a small band. The amplitude gives information on the chromatic dependence of the size of the object, whereas the phase provides valuable information on the position on the sky.

Exactly 202 measurements (87 for Hα, 87 for He I 6678, and 28 for He I 7065) are available with both differential amplitude and phase, with a standard deviation of phase in the continuum smaller than 15°. In the case of a larger standard deviation of the differential phase, we obtained 498 additional measurements of the amplitude of the differential visibility (195 for Hα, 204 for He I 6678, and 99 for He I 7065) with a signal-to-noise ratio better than 5 in the continuum. The details of the observations are presented in Mourard et al. (2018).

The standard differential processing (Mourard et al. 2009) of the VEGA data has been modified to avoid the underestimation of the uncertainties of the differential quantities. For this work, we replaced them with the standard deviation of the measurements (both amplitude and phase) computed in the continuum part and multiplied by a factor equal to the square root of the flux of each narrow-band channel in order to correctly match the behaviour of the photon noise.

One example of an individual measurement is presented in Fig. 1. It should be noted that the amplitude of the differential visibility is normalised to 1 in the continuum and that the phase is arbitrarily set to a mean value of 0 in the continuum. It is also important to note that, in some cases, phase jumps may occur as the differential phase is defined only modulo 2π. The (u, v) coverage is shown in Fig. 2.

Fig. 1. Example of VEGA differential visibility measurement during the night of 27 Jun 2013 with the 140 m E2W2 baseline. Top: Hα line (continuum normalised to 1). Middle: amplitude of the differential visibility (normalised to 1 in the continuum). Bottom: phase of the differential visibility in degrees.

Fig. 2. (u, v) coverage of differential visibility and phase measurements from the CHARA/VEGA instrument. A subset of Hα data used for modelling is shown (black), together with additional He I 6678 (green) and 7065 (blue) data. The orientation is v >  0 north and u >  0 east.

3. Pyshellspec model

To account for all types of observational data, we had to significantly extend and improve our modelling tool, Pyshellspec¹. Its purpose is to calculate radiative transfer through the volume surrounding the binary. We used a joint χ² metric as follows:

$$\begin{array}{c}\hfill {\chi }^2={\chi }_{\mathrm{lc}}^2+{\chi }_{\mathrm{vis}}^2+{\chi }_{\mathrm{clo}}^2+{\chi }_{\mathrm{t}3}^2+{\chi }_{\mathrm{sed}}^2+{\chi }_{\mathrm{spe}}^2+{\chi }_{\mathrm{vamp}}^2+{\chi }_{\mathrm{vphi}}^2,\end{array}$$(1)

with the following individual contributions:

$$\begin{array}{c}\hfill {\chi }_{\mathrm{lc}}^2={\displaystyle \sum _{k=1}^{N_{\mathrm{band}}}\sum _{i=1}^{N_{\mathrm{lc}k}}}{\left(\frac{m_{\mathit{ki}}^{\mathrm{obs}}-m_{\mathit{ki}}^{\mathrm{syn}}}{{\sigma }_{\mathit{ki}}}\right)}^2,\end{array}$$(2)

$$\begin{array}{c}\hfill {\chi }_{\mathrm{vis}}^2={\displaystyle \sum _{i=1}^{N_{\mathrm{vis}}}}{\left(\frac{|V_i^{\mathrm{obs}}{|}^2-{|V_i^{\mathrm{syn}}|}^2}{{\sigma }_i}\right)}^2,\end{array}$$(3)

$$\begin{array}{c}\hfill {\chi }_{\mathrm{clo}}^2={\displaystyle \sum _{i=1}^{N_{\mathrm{clo}}}}{\left(\frac{arg{T}_{3i}^{\mathrm{obs}}-arg{T}_{3i}^{\mathrm{syn}}}{{\sigma }_i}\right)}^2,\end{array}$$(4)

$$\begin{array}{c}\hfill {\chi }_{\mathrm{t}3}^2={\displaystyle \sum _{i=1}^{N_{\mathrm{clo}}}}{\left(\frac{|T_{3i}{|}^{\mathrm{obs}}-{|T_{3i}|}^{\mathrm{syn}}}{{\sigma }_i}\right)}^2,\end{array}$$(5)

$$\begin{array}{c}\hfill {\chi }_{\mathrm{sed}}^2={\displaystyle \sum _{i=1}^{N_{\mathrm{sed}}}}{\left(\frac{F_{\lambda i}^{\mathrm{obs}}-F_{\lambda i}^{\mathrm{syn}}}{{\sigma }_i}\right)}^2,\end{array}$$(6)

$$\begin{array}{c}\hfill {\chi }_{\mathrm{spe}}^2={\displaystyle \sum _{i=1}^{N_{\mathrm{spe}}}}{\left(\frac{I_{\lambda i}^{\mathrm{obs}}-I_{\lambda i}^{\mathrm{syn}}}{{\sigma }_i}\right)}^2,\end{array}$$(7)

$$\begin{array}{c}\hfill {\chi }_{\mathrm{vamp}}^2={\displaystyle \sum _{k=1}^{N_{\mathrm{set}}}\sum _{i=1}^{N_{\mathrm{vamp}k}}}{\left(\frac{V_i^{\mathrm{obs}}-V_i^{\mathrm{syn}}{f}_k}{{\sigma }_i}\right)}^2,\end{array}$$(8)

$$\begin{array}{c}\hfill {\chi }_{\mathrm{vphi}}^2={\displaystyle \sum _{k=1}^{N_{\mathrm{set}}}\sum _{i=1}^{N_{\mathrm{vphi}k}}}{\left(\frac{arg{V}_i^{\mathrm{obs}}-arg{V}_i^{\mathrm{syn}}+g_k+h_k}{{\sigma }_i}\right)}^2,\end{array}$$(9)

where m denotes magnitudes in given passbands, |V|² squared visibility, argT₃ closure phase, |T₃| triple product amplitude, F_λ absolute monochromatic flux, I_λ normalised monochromatic flux, |V| differential visibility amplitude, and argV differential visibility phase. The latter two interferometric quantities are modified by a multiplicative factor f_k, an additive offset g_k, and a phase slip h_k (±360°) to correctly match the way these quantities are estimated as explained in Sect. 2.3.

Apart from new observables, more objects (jet, flow, shell) were implemented in Python, and correspondingly more free parameters. We performed some corrections necessary for the high-resolution spectroscopy, in particular velocity fields of all object are propagated to Shellspec. Phoenix absolute spectra (Husser et al. 2013), used as boundary conditions at the stellar surfaces of our 3D model, were converted from vacuum to air wavelengths and a higher resolution 0.1 Å was used. Fittings of factors, offsets, and slips per each interferometric dataset were included to minimise the difference between observed and synthetic differential visibilities and phases. As an option, we could use the subspace-searching simplex algorithm (or subplex; Rowan 1990) for the χ² minimisation, which is sometimes very useful. The Openmp (threads) parallelisation was applied per wavelength (for LC, VIS, VAMP, VPHI, etc. datasets) or per phase (for SED and SPE). To model extended optically thin structures, we had to extend the grid (usually 80 × 80 R_⊙), and optionally used a lower resolution (2 R_⊙ instead of 1 R_⊙). The majority of rays are in a non-empty space even with these approximations, and thus computationally the task is substantially more demanding than before. For our extensive dataset, we need 3564 synthetic images per iteration, and the number of iterations is about 10³ to achieve a convergence.

Some improvements of the original Shellspec (Budaj & Richards 2004; Budaj et al. 2005; Budaj 2011) were also implemented in Fortran. This includes a radial velocity field in the disc (added on top of the Keplerian field), simple shadowing with prescribed scale height H, which makes it possible to switch on scattering in the disc atmosphere, which is a variable step in the optical depth to prevent integration artefacts. We modified priorities of overlapping objects (jet priority is higher than that of nebulas, and envelope priority is lower than that of nebulas). There is a possibility to use two embedded grids, with a lower resolution for extended structures and a higher resolution in the centre. We also tested asymmetric jets, or temperature gradients in shells.

Nevertheless, we recall all physical properties of our model. We assumed local thermodynamical equilibrium (LTE) level populations and an LTE ionisation equilibrium. The line profile was determined by thermal, microturbulent, natural, Stark, Van der Waals broadenings and the Doppler shift. The continuum opacity is caused by H I bound-free, H I free-free, H⁻ bound-free, H⁻ free-free transitions, the Thomson scattering on free electrons, and the Rayleigh scattering on neutral hydrogen. The scattering processes are implemented only for optically thin environments (single scattering process), with the shadowing mentioned above. The scattering is non-isotropic and is described by the dipole phase function. We also account for the line opacity of Hα, He I, Si II, Ne I, and C II. Abundances were assumed to be either solar, increased up to three times (0.5 dex), sub-solar (in C; Sect. 5.8), or He-rich (Sect. 5.9). We used a small grid of synthetic spectra for the stars, which was generated by Pyterpol (Nemravová et al. 2016) from Phoenix, BSTAR, and OSTAR grids (Husser et al. 2013; Lanz & Hubený 2007, 2003). The stars are subject to the Roche geometry, limb darkening, and gravity darkening (in particular the Roche-filling donor).

On the other hand, we do not include optical irradiation of stars, reflection (because the hot primary is mostly hidden in the disc), Mie absorption on dust, Mie scattering, or dust thermal emission. We consider these missing opacity sources negligible, because temperatures in the system are too high for dust condensation.

As in the previous study, we used the quadratic ephemeris by Ak et al. (2007),

$$\begin{array}{cc}& T_{\min .\mathrm{I}}(\mathrm{HJD})=2408247.968(15)+12.913779(16)·E\hfill \\ \hfill & +3.87265(369)\times {10}^{-6}·E^2,\hfill \end{array}$$(10)

which corresponds to the primary minimum of the optical light curve and in our particular case of β Lyr A the donor (secondary) is behind the gainer (primary; hidden in its opaque disc). Initial conditions for further convergence generally correspond to our previous model based on the optically thick medium (Mourard et al. 2018), although this model did not produce sufficient emission in lines.

Parameter relations. For an easier interpretation of results, we review some of the parameter relations (i.e. the geometrical constraints), as they are implemented in the current version of Shellspec. The disc (also known as nebula) object is described in cylindrical coordinates (R, z; see also Table 2 for basic length scales):

$$\begin{array}{c}\hfill H(R)=h_{\mathrm{cnb}}\sqrt{\frac{\gamma k_{\mathrm{B}}T}{\mu m_{\mathrm{u}}}}\frac{1}{{\mathrm{\Omega }}_{\mathrm{k}}},\end{array}$$(11)

$$\begin{array}{c}\hfill \mathrm{\Sigma }(R)={\mathrm{\Sigma }}_{\mathrm{nb}}{\left(\frac{R}{R_{\mathrm{innb}}}\right)}^{e_{\mathrm{densnb}}},\end{array}$$(12)

$$\begin{array}{c}\hfill \rho (R,0)=\frac{\mathrm{\Sigma }}{\sqrt{2\pi }H},\end{array}$$(13)

$$\begin{array}{c}\hfill \rho (R,z)=\rho (R,0)exp[-min(\frac{z^2}{2H^2};\frac{h_{\mathrm{windnb}}^2}{2})],\end{array}$$(14)

$$\begin{array}{c}\hfill T(R,0)=T_{\mathrm{nb}}{\left(\frac{R}{R_{\mathrm{innb}}}\right)}^{e_{\mathrm{tmpnb}}},\end{array}$$(15)

$$\begin{array}{c}\hfill T(R,z)=T(R,0)max(1;1+(t_{\mathrm{invnb}}-1)\frac{|z|-h_{\mathrm{invnb}}H}{a_{\mathrm{neb}}H-h_{\mathrm{invnb}}H}),\end{array}$$(16)

$$\begin{array}{c}\hfill v_r(R)=\mathcal{H}(|z|-h_{\mathrm{velnb}}H)v_{\mathrm{nb}}{(1-\frac{R_{\mathrm{innb}}}{R})}^{e_{\mathrm{velnb}}},\end{array}$$(17)

$$\begin{array}{c}\hfill v_{\varphi }(R)=\sqrt{\frac{G{M}_{\star }}{R}},\end{array}$$(18)

where H is the scale height, γ the adiabatic exponent, k_B the Boltzmann constant, μ the mean molecular weight, m_u the atomic mass unit, G the gravitational constant, Ω_k = v_ϕ/R the Keplerian angular velocity, Σ the surface density, ρ the volumetric density, T the temperature, v_r the radial velocity, and v_ϕ the azimuthal velocity; ℋ(x) denotes the Heaviside step function.

The jet has a double cone shape with the opening angle a_jet and is described in spherical coordinates (R, θ):

$$\begin{array}{c}\hfill \rho (R)={\rho }_{\mathrm{jt}}{\left(\frac{R_{\mathrm{injt}}}{R}\right)}^2\frac{v_r(R_{\mathrm{injt}})}{v_r(R)}(1\pm a_{\mathrm{symjt}}),\end{array}$$(19)

$$\begin{array}{c}\hfill T(R)=T_{\mathrm{jt}}{\left(\frac{R}{R_{\mathrm{injt}}}\right)}^{e_{\mathrm{tmpjt}}},\end{array}$$(20)

$$\begin{array}{c}\hfill v_r(R)=v_{\mathrm{jt}}{(1-\frac{R_{\mathrm{cjt}}}{R})}^{e_{\mathrm{veljt}}}.\end{array}$$(21)

In our case, the base plane corresponds to the orbital plane, and the cone position is determined by the radial offset R_poljt and the polar angle α_jet. Similarly, the shell is also described in spherical coordinates:

$$\begin{array}{c}\hfill \rho (R)={\rho }_{\mathrm{sh}}{\left(\frac{R_{\mathrm{insh}}}{R}\right)}^2\frac{v_r(R_{\mathrm{insh}})}{v_r(R)},\end{array}$$(22)

$$\begin{array}{c}\hfill T(R)=T_{\mathrm{sh}}{\left(\frac{R}{R_{\mathrm{insh}}}\right)}^{e_{\mathrm{tmpsh}}},\end{array}$$(23)

$$\begin{array}{c}\hfill v_r(R)=v_{\mathrm{sh}}{(1-\frac{R_{\mathrm{csh}}}{R})}^{e_{\mathrm{velsh}}}.\end{array}$$(24)

We assume the spherical shell is centred on the primary and also encompasses other objects. It represents circumstellar matter that escaped farther away from the binary. All remaining parameters are explained in the caption of Table 1.

3.1. Doppler tomography

In order to create model spectra with enough emission, we started with two Ondřejov spectra taken at 0.288 and 0.785 phases, that is, out of eclipses, and converged our new model with additional optically thin objects. This simplified Doppler tomography was carried out to verify that the model is indeed capable of fitting the Hα profile. There is always a question of which objects should be included in the model and which should not. If the model were too simplistic, the objects would be distorted; if it were too complex, the objects could be unconstrained. After some preliminary tests, we used five objects (primary, secondary, disc, jet, shell) out of eight (spot, envelope, flow)².

Given the observed Hα profiles and their overall width, our model should include a large positive velocity with respect to the line of sight or gradient of v. Because H I, Si II, as well as He I, and Ne I are excited, we expect both low and high temperatures, or a gradient of T in the circumstellar medium. Apart from absorption lines arising in stellar atmospheres, the model is capable of creating a P Cygni profile due to winds, either in a disc atmosphere or in a surrounding shell. Alternatively, line profiles may be formed by overlapping velocity fields, in accord with the priorities of objects. We have to look for suitable net velocities of whole objects, but also for turbulent velocities, which significantly affect the optical depth along the line of sight; one should converge both at the same time.

Results of the first two-spectra model are shown in Fig. 3. The model easily created enough emission, and the EW of Hα is fitted very well. There are relatively minor systematics in the Hα profile, but major systematics can be seen for other spectral lines, especially He I 6678. The Si II 6347, 6371, and Ne I 6402 model lines have lower EW and depth than the observed ones, because there is a tension between the overall emission and the respective absorption. The region between 6500 and 6550 Å contains telluric lines that are not included in our model, but they should not affect the convergence in a negative way. Although we varied the chemical composition, there might be non-LTE effects (for He I) or some unaccounted for temperature gradients.

Fig. 3. Normalised spectra of β Lyr A for two out-of-eclipse phases (0.288 and 0.785). Observed data with uncertainties are indicated (blue), as well as synthetic data (yellow), residua (red), and large χ2 contributions (red circles). This 1st model was only fitted to these two spectra. The Hα emission profile and its EW are well described, but there are systematic differences for the He I 6678 line; the synthetic Si II 6347, 6371 and Ne I 6402 lines have low EW.

A more representative set of 11 spectra covering a representative range of orbital phases was fitted in the second step. The results are shown in Fig. 4. All optically thin objects (disc atmosphere, jet, shell) contribute substantially to the Hα emission flux. Because the spectra are normalised to the continuum flux, which is larger outside eclipses, the emission appears to vary in strength twice each orbit. The synthetic Hα profiles for this second model exhibit a variability similar to the observed ones, although at several phases there are systematic differences (both positive and negative). Comments related to He I, Si II, Ne I lines remain essentially the same.

Fig. 4. Normalised spectra for 11 phases covering the whole light curve, including primary and secondary eclipses. Colours in this figure, as well as in others, have the same meaning as in Fig. 3. The synthetic Hα profiles for the second model exhibit a variability similar to the observed ones, although at several phases there are systematic differences (both + and −).

It might seem easy to improve the fit further, but it is not the case for a geometrically constrained model, where all parameters have either geometrical or physical limits. We are practically sure the convergence works and that it is not a matter of one local minimum of χ². In order to improve the fit, we may need to relax some of our assumptions (e.g. the axial symmetry of the disc, the vertical symmetry of the jets, or the radial symmetry of the shell). We also did not account for any intrinsic variability of the source. However, we preferred to keep our model as simple as possible, at least at this stage.

3.2. Differential interferometry

Interestingly, all interferometric data indicate a decrease of the visibility amplitude |dV| when scanning across the Hα profile (see Fig. 5), and it seems to be almost independent on baseline length and orientation (Fig. 2). Such a general finding means that the core of the Hα emitting region is clearly resolved for all baseline lengths between 50 and 200 m. Moreover, the respective velocities must be large enough to occur in the wings of Hα.

Fig. 5. Observed differential visibility amplitude |dV| versus wavelength λ, normalised to 1 in the continuum (blue), and its decrease across the Hα profile. Uncertainties of |dV| are also plotted (grey). Synthetic visibilities (yellow) are shown for the two “extreme” values of the shell’s outer radius Routsh = 40 R⊙ (top) and 120 R⊙ (bottom).

Only an extended symmetric shell may cause |dV| to decrease. On the contrary, a disc (nebula) or jets emit usually from small (hot) areas, and they are both asymmetric, which would force |dV| to increase, at least for the shortest baselines. This is not observed. Our preliminary tests thus demonstrate the need to include the shell in our model.

Regarding the differential phase argdV measurements, it should be noted that some phase wrappings (±360°) are present in the data (cf. e.g. datasets 1, 11, and 12), which should not be a problem as we account for them in the model. The differential phases are obtained on three baselines, with two of them (E1E2 and E2W2) being oriented almost perpendicularly to the orbital plane, whereas the third one (W1W2) is very close in orientation to the orbital plane. Interestingly, the worst fits to the model are obtained systematically for this last orientation. This general finding is in agreement with the fact that the differential phases in Hα are dominated by the jets, but we should also conclude that our geometrically-constrained model is not flexible enough to explain all the observed features.

3.3. A joint “compromise” model

Our complete dataset is very heterogeneous. On one hand, this is an advantage that allows us to construct a very robust model of β Lyr A when everything is fitted together. On the other hand, when some types of measurements exhibit systematics, as mentioned above, the χ² contributions of the joint model would be worse, while observation-specific models would be better (see Sect. 3.4).

After 2761 iterations (Fig. 6), we obtained a model with reduced ${\chi }_{\mathrm{R}}^2={\chi }_{\mathrm{N}}^2$ values, which are summarised in Table 1 (first column). Previously used datasets are fitted only slightly worse than in Mourard et al. (2018), with ${\chi }_{\mathrm{lc}}^2=3.1$, ${\chi }_{\mathrm{vic}}^2=4.0$, ${\chi }_{\mathrm{clo}}^2=3.4$, ${\chi }_{\mathrm{t}3}^2=5.6$; the new datasets resulted in ${\chi }_{\mathrm{sed}}^2=4.7$, ${\chi }_{\mathrm{spe}}^2=44.1$, ${\chi }_{\mathrm{vamp}}^2=4.5$, ${\chi }_{\mathrm{vphi}}^2=44.0$, and the total ${\chi }_{\mathrm{R}}^2=17.0$.

Fig. 6. χ2 convergence (red) for joint model; individual contributions (LC, VIS, CLO, T3, SED, SPE, VAMP, VPHI) are also indicated. The model successfully converges to a local minimum. Some datasets have a substantially larger number of observations, (i.e. effectively a larger weight). The χ2 values are different from Table 1 because the model was re-converged several times, and uncertainties of some datasets were modified.

We have to explain why some contributions are large. For ${\chi }_{\mathrm{sed}}^2$, there is a difference in NUV below the Balmer jump, which is caused by a relatively low number of measurements compared to other datasets and thus a low weight. Lines are not fitted, because the flux is computed at monochromatic wavelengths; we would have to use a wavelength resolution about five times higher and perform a convolution with the instrumental profile of Burnashev & Skulskii (1978). Nevertheless, the line-to-continuum flux ratio is fully described by our SPE dataset. In the case of ${\chi }_{\mathrm{spe}}^2$, there are systematic differences at certain phases (0.288, 0.825), although others are good fits (0.454, 0.548, 0.712). A substantial contribution arises from high-temperature He I and Ne I lines and also from telluric lines, which are not fitted by our model. For ${\chi }_{\mathrm{vamp}}^2$, synthetic |dV| sometimes exhibit a narrower decrease (cf. dataset 1), or a peak in the middle of Hα (2, 3), although others are almost perfect fits (4, 5, 6, …). Finally, ${\chi }_{\mathrm{vphi}}^2$ is substantially increased because synthetic argdV are sometimes smoother (4, 6), there are possibly remaining phase slips (1, 11, 12), or mirroring of phases (2, 3); these numerous measurements have relatively high weight.

A visual comparison of all observed and synthetic datasets is shown in Figs. 7–14. The resulting geometrical model of β Lyr A in the continuum is shown in Fig. 15. In particular, we see optically thick objects – the primary, the secondary, and the disc – and partly also the jets, but not the tenuous shell. The same model for the wavelength range of Hα is shown in Fig. 16. We can clearly see optically-thin circumstellar matter emitting in Hα, including the velocity field. In the following, we describe individual components of our model, as inferred from the observations:

Fig. 7. Comparison of observed and synthetic phased light curves from FIR to FUV, computed for joint (“compromise”) model. The names of datasets are shown in the right column. The light curves were arbitrarily shifted in the vertical direction.

Fig. 8. Squared visibility |V2| (shifted by dataset number) for different baselines B/λ (in cycles). The CHARA/MIRC observations are at the bottom, NPOI in the middle, and CHARA/VEGA at the top (blue); synthetic data for the joint model are plotted for comparison (yellow). Overall trends (|V2|(B/λ)) seem to be correctly described, although there are some systematics for small groups of data.

Fig. 9. Similar to Fig. 8, for the closure phase argT3 vs. B/λ.

Fig. 10. Similar to Fig. 8, for the triple product |T3| vs. B/λ.

Fig. 11. Spectral-energy distribution (SED) overplotted for ten phases. Synthetic SEDs systematically differ from observations, especially in NUV and in the vicinity of the Balmer jump. This is not the best model in terms of ${\chi }_{\mathrm{sed}}^2$ (however, cf. Fig. A.5).

Fig. 12. Normalised spectra I(λ) for 11 phases. Both observed (thick) and synthetic (thin) spectra are plotted for comparison. Colours correspond to the phases (see also, values on the left). Uncertainties are not plotted, but they are of the order of 0.01. The vertical lines correspond to (from left to right): Si II 6347, Si II 6371, Ne I 6402, Hα, and He I 6678 wavelengths.

Fig. 13. Differential visibilities |dV| vs. λ normalised in continuum. Synthetic (yellow) visibilities for the joint model exhibit a similar decrease across the Hα profile to the observed ones (blue). However, some of the synthetic |dV|s are systematically lower in the blue wing and higher in the red wing.

Fig. 14. Similar to Fig. 13, for the differential phases argdV vs. λ. On average, phase changes are comparable in both observed and synthetic data. There are some remaining phase slips that are also fitted for.

Fig. 15. Continuum synthetic images of β Lyr A for joint model, computed for four monochromatic wavelengths (from top to bottom): 155 nm (FUV), 545 nm (V), 1630 nm (H), 4750 nm (FIR). The orbital phase is always 0.25. The appearance of the system changes substantially. Optically thin components (disc atmosphere, jet, shell) are not seen very well in continuum radiation.

Fig. 16. Line-profile synthetic images of β Lyr A computed for wavelength range of Hα (i.e. from 655.1 to 657.8 nm, at a fixed phase 0.25). Optically thin components are clearly visible. The disc atmosphere appears first, because its Keplerian velocities close to the inner rim produce the blue-shifted wing. The jet inclined towards the observer appears second, with high velocities being projected to the line of sight. Finally, there is the spherical shell, with relatively low velocities spanning the core of Hα, which obscures other small-scale structures. For λ > 6563 Å, all objects disappear in a reverse order. An animated version is avaiable online and at https://sirrah.troja.mff.cuni.cz/~mira/betalyr/.

Primary. The gainer is an object for which we fixed several parameters (Table 2). It is mostly hidden in the disc, but as one can see in the figures, its polar region is visible. It is the source of hot radiation that is also scattered by the circumstellar medium (CSM) towards the observer.

Secondary. The donor is filling the Roche lobe with limb and gravity darkening. The limb darkening coefficient is interpolated for given λ from van Hamme (1993) tables. The gravity darkening parameter was set to a value suitable for non-convective atmospheres of stars (0.25). One can see that the regions near the L1 point are indeed dimmer because of it. Its polar temperature is ∼14 000 K and the inferred polar radius about 14.1 R_⊙.

Disc. This is an axially symmetric accretion disc centred on the gainer. It has an outer radius of 31.5 R_⊙ and almost fills the Roche lobe. The density profile decreases with radius, and its slope (−0.57) is slightly less steep than in some Algols, where it attains −1.0 (Budaj et al. 2005; Atwood-Stone et al. 2012). The temperature profile also decreases, and its slope (−0.73) is slightly steeper than a typical profile due to irradiation (−0.5), but it is in surprisingly good agreement with the theoretical temperature profile of the steady viscous accretion discs (−0.75; Pringle 1981). The temperature at the inner rim of the disc reaches ∼30 000 K, which is a very reasonable value, comparable with the temperature of the gainer. The temperature inversion reaches 1.5, which means that the temperature increases in the vertical direction by this factor. This is most probably caused by the irradiation of the disc atmosphere. For this reason, we see that the disc is brighter on the top and bottom and dimmer in the middle. Moreover, the atmosphere scatters the radiation from the gainer.

The parameter h_cnb = 3.8 means that the vertical scale height is multiplied by this factor and is more extended than the expected equilibrium value. This may be due to non-negligible hydrodynamic flows within the disc. There is a significant radial velocity component in the surface layers (v_nb = 112 km s⁻¹), which might be due to stellar wind or radiative acceleration. The turbulence is relatively low (v_trbnb = 11 km s⁻¹), which means that Keplerian and radial components describe the velocity field very well.

Jet. There are two conical jets perpendicular to the orbital plane. They do not seem to be associated with the polar regions of the gainer (cf. R_poljt). The “net” velocity v_poljt assigned to the jets was treated as a free parameter, and we thus have to discuss whether its value 10 km s⁻¹ is reasonable or not. Because the projected orbital velocities of the primary and secondary are K₁ = 41 km s⁻¹ and K₂ = −186 km s⁻¹ (at the phase 0.25), we consider it to be reasonable, although this indicates that the jets may not follow Keplerian velocities at the disc’s rim.

The terminal (expansion) velocity v_jt is almost 700 km s⁻¹, and the respective exponent (e_eveljt = 1.27) is slightly smaller than that of the shell or disc. Turbulence is about ten times smaller than the terminal velocity, which indicates that the velocity field is fitted reasonably well. This object is optically thin in continuum, so it is constrained mainly by observations in the Hα line. Without this object, the reduced ${\chi }_{\mathrm{R}}^2$ would increase up to 59, which justifies its role in our model (although a re-convergence might decrease it again).

Shell. This spherical object extends to more than 70 R_⊙. The respective net velocity is low (−5 km s⁻¹); the shell may not be exactly centred and co-moving with the primary. Its terminal velocity is very low (only 79 km s⁻¹), but it is interesting that its velocity exponent is similar to that of the disc (∼1.9). On the other hand, turbulence is very high (v_trbsh = 102 km s⁻¹), which indicates that the velocity field is not well described by our formulation (Eq. (24)). This is probably not surprising given that it fills a broad spherical region in the vicinity of the orbiting stars where gravitational potential is far from being isotropic and radial. What is surprising is that our data indicate only a small radial temperature gradient (cf. e_tmpsh). It is also optically thin and constrained mainly by interferometry in Hα. If this object were excluded from the model, the reduced ${\chi }_{\mathrm{R}}^2$ would increase to 62, which also justifies its role.

Generally, the new model seems to be compatible with our previous model (Mourard et al. 2018), but we should emphasise that the major difference is the firm detection of previously conjectured structures (the jets and the shell), which was possible thanks to spectro-interferometric and spectral observations in the Hα region. There are minor differences, however: the disc’s outer rim radius is larger (R_outnb = 31.5 vs. 30 R_⊙ previously), the thickness slightly smaller (h_cnb = 3.8 vs. 4.3), and the orbital inclination also larger (i = 96° vs 93.5°). All these differences may be enforced by the need for emission in the Hα line, which is enhanced if the disc is more extended and more inclined. For jets perpendicular to the disc, larger i leads to larger line-of-sight velocities and also to larger asymmetry due to the obscuration by the disc. This is required in order to explain the asymmetric Hα profile. In this particular model, the asymmetries in the Hα profile arise mostly from overlapping velocity fields of objects with increasing priorities (shell  →  disc  →  jet).

3.4. Observation-specific models

Starting from the “compromise” model above, we converged the model again to fit individual datasets to understand the trends and potential disagreements. From the χ² convergence (Fig. 17), it is evident that our model is indeed capable of fitting individual datasets better. For example, the reduced ${\chi }_{\mathrm{lc}}^2$ can easily reach 2.2 (instead of 3.1). For other χ² values, see Table 1 (last row). Consequently, we think there are either systematic differences between datasets, or our (complex) model is still not complete. We may miss additional objects, some asymmetries, or a temporal variability. Alternatively, we may modify weights of individual datasets (and use e.g. w_lc = 10), but this does not “solve” the problem, of course.

Fig. 17. χ2 convergence for observation-specific models, when starting from the best-fit joint model. After performing up to 103 iterations, substantial improvements were achieved for some datasets (LC, T3, SED, VAMP, etc.). It confirms systematic differences between observational datasets.

The results of all observation-specific models are summarised in Table 1 (columns LC to VPHI). For the V band, it is also possible to compare the models visually (Fig. 18). The differences are demonstrated as the disc thickness and the intensity of its outer edge, which is proportional to the temperature profile T(r). In most models, the primary is directly visible, but the LC dataset tends to produce a continuum emission from a more extended hot area.

Fig. 18. Continuum synthetic images for observation-specific models (datasets LC, VIS, CLO, T3, SED, SPE, VAMP, and VPHI) for the wavelength 545 nm (V). The apparent differences (e.g. the thickness of the disc, the appearance of the primary) demonstrate systematics between datasets. Alternately, some datasets (e.g. VPHI) do not constrain certain parameters.

Given these results, uncertainties of model parameters (cf. Table 1, last column) were determined as the maximum differences between the joint model and relevant observation-specific models, because they are almost certainly dominated by systematics, not by the extent of local χ² minima, not even by the global one. To be more specific, σs for all velocities can only be constrained by the joint model and the SPE, VAMP, and VPHI datasets; similarly, σs for flux-related quantities (𝜚, T, R) can hardly be constrained by relative measurements.

4. Stellar evolution modelling

To understand the evolutionary stage of β Lyr A binary and its relation to the observations, we performed a simplified 1D modelling with the MESA stellar evolution program (Paxton et al. 2011, 2015). In Fig. 19, we present a nominal evolution of a binary with the initial masses M₁ = 10 M_⊙, M₂ = 6 M_⊙, and the orbital period P = 7 d. We assumed the solar composition. We used an explicit Ritter scheme and we restricted the mass accretion rate Ṁ up to 10⁻³ M_⊙ yr⁻¹. Although we performed a small survey of parameters (P = 5 to 50 d, M₁ = 8 to 10 M_⊙, while keeping M₁ + M₂ = const.), we restrict our discussion to the nominal case, because other values lead to binaries incompatible with β Lyr A, or to a common-envelope phase that is difficult to describe in 1D. The mass ratio q is already reversed, so index 1 corresponds to the observed secondary (donor), and 2 to the primary (gainer).

Fig. 19. Stellar evolution of a binary corresponding to β Lyr A. The component radii R1, R2 vs. the primary mass M1 are plotted (top), with observed values indicated by full circles, together with the distance a and R1, R2 vs. time t (middle); see horizontal dashed lines, and surface abundances X, Y, AC12, AO16 vs. time t (bottom). At t = 19.613 × 106 yr, the mass ratio q corresponds to the observed one.

A comparison with the observed values of a, R₁, R₂ shows some differences: an offset 10³ yr between the time of best-fit for a and the best-fit for M₁, M₂ (see dotted vertical lines). The synthetic radii are a factor of 1.2 and 2 larger at this moment. However, we recall that it is only a 1D model, without an accretion disc. Consequently, we consider these differences to be acceptable. The model of van Rensbergen & De Greve (2016) produced a qualitatively similar HRD for the gainer, but their initial period was shorter, P = 2.36 d.

At this stage, the surface chemical composition is already modified. At t = 19.612 × 10⁶ yr, there is a low C abundance (by a factor of 10²) and a high N abundance (by a factor of 5). At t = 19.613 × 10⁶ yr, even He abundance is increased up to 0.36 (and H is correspondingly decreased). For our modelling, it means we should also test models with a substantially modified chemical composition. This is in accord with Balachandran et al. (1986), who suggested He enrichment N(H) = 0.4, N(He) = 0.6 (by number), and also N to be overabundant, and C, O to be underabundant, namely C/N ≤ 0.11, O/N ≤ 0.25.

Interestingly, further evolution would lead to a detached system with a stripped He dwarf (secondary). A hot sub-dwarf of the sdB or sdO type is expected (Heber 2009; Lei et al. 2018). ϕ Per binary might be just in this (late) evolutionary stage (Mourard et al. 2015) and a dedicated comparative study might be very useful.

5. Alternative models

5.1. High-resolution model

For a spatial resolution increased twice to 1 R_⊙, some datasets are fitted even better: for example, ${\chi }_{\mathrm{lc}}^2$ (not reduced) decreased from 7083 to 6720 because the primary and the disc rim are better resolved and thus contribute more to FUV and NUV fluxes. On contrary, ${\chi }_{\mathrm{sed}}^2$ increased from 8578 to 9832, due to the same reasons. The most sensitive term seems to be ${\chi }_{\mathrm{spe}}^2$, which increased substantially from 588454 to 682016, because Hα line profiles are slightly “sharper” and the He I 6678 emission is enhanced; although the profiles remain qualitatively very similar. Other contributions are slightly decreased. This conclusion is preliminary, though, without a repeated convergence. In principle, even the high-resolution model could be converged again, which would decrease the increased χ². We conclude the model is resolution dependent. However, this is not necessarily a bad thing; a low-resolution model may simply represent shallower gradients, or less sharp transitions between objects.

5.2. Distance fixed to 294 pc

First, we compare our distance with that of Bastian (2019), who discovered the Gaia 8 cluster, with parallaxes around 3.4 mas or distance 294 pc, and the intrinsic dispersion of only 0.06 mas (i.e. 5 pc radial, 1° spatial). β Lyr A is located in the middle (spatially); if it is also in the centre of mass, then our photo-spectro-interferometric value d = (328 ± 7) pc is substantially larger.

If we fix the distance in our model to d = 294 pc instead, and converge the model, we obtain different parameters, of course. In particular, in Table 4 (column 294 pc) we can see that secondary temperature T_cp = 13 512 K is lower, and the disc’s outer radius R_outnb = 35.2 R_⊙ is larger. While the overall fit seems better, ${\chi }_{\mathrm{R}}^2=16.9$ (as compared to 17.0), mostly because ${\chi }_{\mathrm{spe}}^2$ and ${\chi }_{\mathrm{vamp}}^2$ contributions were improved, it is at the expense of other terms being much worse. Notably, ${\chi }_{\mathrm{sed}}^2$ was increased three times, which is unacceptable for us, and ${\chi }_{\mathrm{vic}}^2$ and ${\chi }_{\mathrm{clo}}^2$ were also increased. Moreover, the disc’s outer radius is too large and overshoots not only the tidal cut-off radius 26.3 R_⊙ (possible in principle as the disc is not an isolated system during ongoing mass transfer), but also its Roche lobe. This is why we still prefer the original model.

5.3. The mass ratio

We cannot make the primary M_⋆ and secondary mass M_cp free when we keep asini and the period P (as well as Ṗ) fixed. Nevertheless, if we free the mass ratio q (and also all other parameters), it mainly affects the size of the donor, which is the major source of light. We obtained parameters shown in Table 4 (column QRATIO). While the majority of them remained close to the previous (local) minimum, the value of q = 0.2177 is lower. Because we introduced one more parameter, it is logical that the overall fit is better, with ${\chi }_{\mathrm{R}}^2=16.5$. However, we do not consider the respective changes of parameters to be substantial.

5.4. Spot-like asymmetry

Similarly to Mourard et al. (2018), we introduced a spot to our model, which represents an additional spherical object. We converged not only spot parameters, but also all other parameters, to be sure that all objects can adapt to new geometrical constraints. The results are listed in Table 3 (for the spot itself) as well as in Table 4 (column SPOT). We see a minor improvement of the light curve, spectra, and differential interferometry, however, this is at the expense of other datasets (the reduced ${\chi }_{\mathrm{R}}^2$ decreased to 16.8). The position of the spot converged close to the donor–gainer line and the distance corresponds to the outer radius of the disc. Consequently, such a spot may represent either a part of the flow from the donor, the base of the jets, or an asymmetry of the disc rim. We consider the existence of this spot likely, although not as prominent as before (cf. Mourard et al. 2018), because there are additional objects in our model.

5.5. Asymmetric jets

Our stringent geometrical constraints may be partly relaxed by using the parameter a_symjt, which allows for an asymmetry of the jets (Eq. (19)). The resulting model is shown in Table 4, column ASYMJT. As before, we introduced one more parameter, and it is not surprising that the fit is better, with ${\chi }_{\mathrm{R}}^2=16.5$. However, the resulting value a_symjt = 0.02 is not far from zero. It seems that this additional parameter actually allowed for tiny adjustments of other parameters, and we cannot conclude that jets are asymmetric.

5.6. Jet temperature gradients

The fit of the He I 6678 line is far from being perfect. To improve it, we used a model with a substantial temperature gradient T(r) in the jets. We obtained a significantly better fit (${\chi }_{\mathrm{R}}^2=16.4$; see Table 4, column ETMPJT) with the slope e_tmpjt = −0.55 and the temperature at the base of the jets up to T_jt ≐ 30 000 K. It seems to improve both Hα and He I line profiles. Although the systematics in the core of He I line remained qualitatively the same, there are improvements in the wings of both lines. The parameters of these jets are reasonable because their temperature corresponds very well to the temperature of the gainer and the exponent to the heating by irradiation from the gainer. Consequently, it may be considered as our preferred model.

5.7. Shell temperature gradients

We also tried to enforce the temperature gradient in the shell by initially decreasing the slope e_tmpsh = −0.5 and adjusting the temperature T_sh  =  20 000 K accordingly (finishing with 6000 K at R_outsh). All parameters were converged again and the result is shown in Table 4 (column ETMPSH). We can see the gradient is preserved, but the fit is worse, with ${\chi }_{\mathrm{R}}^2=18.1$, especially the ${\chi }_{\mathrm{spe}}^2$ term. The problem is that the gradient also affects all other lines (Hα, Si II, Ne I) and creates an excess emission that prevents further convergence.

5.8. Low C abundance

We tried to use higher abundances for Si and Ne in order to explain the depth of the respective lines. On the other hand, according to Sect. 4, there might be a 100 times lower abundance of C in surface layers and consequently in the CSM. Our spectra do contain the region of C II 6578 and 6583 lines, but they are too weak. Nevertheless, we checked synthetic spectra with these transitions and solar abundance (3.31 × 10⁻⁴ by the number of atoms; Grevesse & Sauval 1998). It turned out it would create such strong C II emission, with a similar width as Hα, that the abundance must be low (see Fig. 20). The abundance 10⁻² of the solar value is fully compatible with observations; the upper limit is about 10⁻¹. Subsequently, these low abundances were also applied in the “joint” model.

Fig. 20. Synthetic spectra for solar abundances (yellow) and 10−2 lower abundance of C (grey); observed spectra (blue) are plotted for comparison. For the solar composition, C II 6578 and 6583 emission is too strong.

5.9. He-rich abundance

If we wish to use the non-solar abundances of Balachandran et al. (1986), the situation is much more complicated. We must not use standard atmospheric models, because the abundances N(H) = 0.4 and N(He) = 0.6 change their hydrostatic profiles, and the emerging spectra must be different.

As a preliminary check, we computed several NLTE models with the Tlusty program (Hubený & Lanz 1995, 2017). For the modified chemical composition, it was necessary to improve the convergence by adjusting several parameters. ³ The output of Tlusty was then used as an input for the Synspec program (Lanz & Hubený 2007; Hubený & Lanz 2017) to obtain a detailed synthetic spectrum with the line list of Kurucz. We computed spectra only for the effective temperature 14 000 K and a limited wavelength range of 6330 to 6700 Å (see Fig. 21). The level of continuum is higher by 3% for He-rich abundances. The (non-rotated) Hα line profile shows Lorentzian wings deeper by 6%, and the He I 6678 line exhibits a significant central absorption. Luckily, after rotational broadening, the profiles will not be so different from standard ones, and we proceed can thus without an extensive computation of a new grid of synthetic spectra. Moreover, one of our stars (primary, gainer) is partly hidden within the disc, and its radiation is reprocessed by circumstellar matter.

Fig. 21. Non-rotated NLTE synthetic spectra computed for stellar atmospheres with He-rich (N(H = 0.4), N(He) = 0.6; black) and solar (orange) compositions. The effective temperature was 14 000 K in both cases. For the He-rich case, the profile exhibits a central absorption.

For simplicity, we thus only used He-rich abundances for the CSM. The result is shown in Table 4 (column H0.4_He0.6). The model does converge, but the overall ${\chi }_{\mathrm{R}}^2=19.5$ remained high, with ${\chi }_{\mathrm{spe}}^2$ being the dominant term. On the basis of our modelling, we thus cannot confirm that abundances are He-rich.

5.10. Differential visibility in He I 6678

There are additional interferometric datasets, namely differential visibilities in the He I 6678 line. We performed a comparison only, not convergence. According to Figs. 22 and 23, there are systematic differences, with synthetic |dV|’s being often flatter than observed |dV|s. The problem is likely the same as in Sect. 5.7 (i.e. unidentified temperature gradients, which would create extended hot emission regions of He I but not of Hα). The comparison of these two figures with the similar ones for Hα (Figs. 13 and 14) is however very instructive. It is easily seen that the He I data are firstly less resolved in the core of the line than in Hα, and secondly that almost no phase signal is detected. This result is in favour of the dominating behaviour of the shell structure, with the jets playing no role here.

Fig. 22. Differential visibility amplitude |dV|. for He I 6678 line. The joint model was used, but not the convergence of the respective He I VAMP dataset.

Fig. 23. Differential visibility phase argdV for the He I line. Again, the joint model, but no convergence.

5.11. Additional high-resolution spectroscopy

There are also additional spectroscopic datasets, particularly the one obtained at the Ritter Observatory and used by Ignace et al. (2018). It has a higher resolution than our spectra (λ/Δλ = 26 000) and a very good phase coverage. The spectra were acquired in different seasons (1996 to 2000). When we performed a comparison, not convergence, of 11 representative spectra, there were systematic differences, mainly in the observed emission peak around the primary eclipse (i.e. the donor eclipsed, not the gainer) which is substantially higher than in our model. The observed Hα profiles also do contain smaller features that are not reproduced by our model. Nevertheless, the model can easily be adapted to the overall emission (e.g. by adjusting densities ρ_nb, ρ_jt, ρ_sh, or by moving the jet along with α_jt ≃ −140° where a second local minimum of ${\chi }_{\mathrm{spe}}^2$ is located). It may be an indication that the distribution of optically-thin CSM had been evolving on the time scale of >10 years (or >300 orbits). Several investigators, most recently Rucinski et al. (2019), also noted cycle-to-cycle changes in the shape of the binary light curve. If true, our model should be treated as a “time-averaged” snapshot of the system.

6. Conclusions (and problems)

We presented a geometrically constrained model of β Lyr A that takes into account all types of available observational data. It contains the primary, the Roche-filling secondary, the optically thick disc, its optically thin atmosphere, the jets, and the shell. We determined absolute sizes of all the components, physical properties (ρ, T, v profiles), and the distance to the system. They are summarised in Table 1.

Some of the parameters of the joint model are close to their maximum or minimum values, in particular T_cp, R_outnb, e_velsh. This is an independent indication that our model is not yet complete and that we may miss some features. For example, the outer radius R_outnb of the disc (nebula) almost touches the Roche lobe. This may induce perturbations and a precession of the disc. These instabilities are not accounted for in our model. Additionally, the outer rim may not be in an exact equilibrium, because of the ongoing mass transfer, and the secondary may induce spiral arms, that is, azimuthal variations in the disc (Panoglou et al. 2019).

The mass loss rate from jets is substantial. Given the surface area at the beginning of the cone, $S=2\pi R_{\mathrm{injt}}^2(1-\cos a_{\mathrm{jet}})$, and the respective expansion velocity, we get Ṁ_jt = 2Sv_rρ_jt ≃8.3 × 10⁻⁷ M_⊙ yr⁻¹, which is about 4% of the mass transfer rate $\dot{M}\simeq 2\times {10}^{-5}{M}_{\odot }{\mathrm{yr}}^{-1}$. Consequently, the mass transfer is not conservative, but it is not far from being so. The timescale related to the jets is τ ≃ R_outjt/v_jt ≃ 0.5 d, which is shorter than the orbital period. The jets are continuously replenished as they follow the orbital motion. It is interesting to integrate the mass loss over long timescales and to observationally check where the (expanded and cooled-down) CSM is located. According to measurements of the radio emission by Umana et al. (2000), the CSM is very extended (145 mas =̂ 10⁴ R_⊙ at our d), and the integrated mass is M ≃ 0.015 M_⊙. Consequently, the timescale would be of the order of τ ≃ M/Ṁ_jt ≃ 10⁵ yr, which agrees with the binary evolution timescale (van Rensbergen & De Greve 2016).

We kept masses of both components more or less fixed, in accord with previous spectroscopic analyses of individual absorption lines (Si II, Ne I). The distance is then determined mainly from interferometry and the SED, but spectroscopy (SPE) is also affected, because Keplerian velocity fields are determined by central masses. We recall that our preferred distance d = (328 ± 7) pc is larger than the 294 pc inferred by Bastian (2019; see Sect. 5.2).

Looking at the synthetic profiles of the He I 6678 line in detail, it is in a broad emission (due to the inner hot edge of the disc and Keplerian broadening), but the observed profile is steeper, double peaked (with a red peak being stronger), and has a blue-shifted central absorption (i.e. very similar to Hα). Consequently, the inner edge should be even more visible and should exhibit some absorption due to winds. Our model cannot easily adapt to this, because the emission in Hα would be immediately increased, and we already match its EW. Moreover, the differences in Hα and He I are of the same order, which is an indication of a compromise.

On contrary, the synthetic Si II 6347 and 6371 lines tend to produce a broad emission and a weak absorption, even at an increased metallicity, but the observed Si IIs exhibit deeper absorptions. Similarly, the synthetic Ne I 6402 line is in a broad (disc) emission, but the observed profile is rather flat.

According to a standard stellar evolution, the metallicity of the stellar surfaces – as well as of the CSM – can be substantially different from normal (solar) if the mass transfer have reached chemically modified layers. For the parameters corresponding to our model, we expect a low abundance of C (by a factor of 10²). A likely final outcome would be a detached system with a He-rich dwarf (similar to ϕ Per; Mourard et al. 2015).

From a technical point of view, our model is somewhat resolution dependent. The peak densities or temperatures are not necessarily well resolved; in other words, the profiles are effectively shallower/smoother. A higher resolution may be also needed to obtain smooth P Cygni profiles in thin layers with velocity gradients, such as in expanding atmospheres.

In order to improve the convergence of our model, it may be useful to use the least correlated parameters. For example, the total mass of the shell (instead of R_insh, R_routsh, and ρ_sh), or a suitable reference radius (between R_innb, R_outnb; inclusive), in accord with observational datasets, which are sensitive either to outer, or inner radii. In our case, outer radii would seem more appropriate.

Nevertheless, in spite of the remaining limitations discussed above, it is encouraging that our modelling of an extended set of different types of observational data led to a generally consistent quantitative picture of the system. We note that the principal physical properties obtained from several specific considered models are numerically quite stable and not too different from one model to another. Our models nicely confirm the conjecture that the so-called B spectrum (introduced and discussed in the classical early studies of β Lyr A) originates mainly in the jets, the disc atmosphere, or other circumstellar matter above/below the orbital plane, including their blueshift of about 50–100 km s⁻¹ (Harmanec 1992; Harmanec et al. 1996; Bonneau et al. 2011). Additionally, the presence of an extended shell, observationally detected by Ak et al. (2007), is required by our model and data. Another important result of this study is a convincing confirmation of the strong carbon underabundance, in accord with the models of the large-scale mass exchange in binaries. We thus believe that the present study constitutes a good starting point for even more sophisticated modelling of the system.

Movie

Movie of Fig. 16 Access here

Acknowledgments

We thank J. A. Nemravová for the initial development of the Pyshellspec. We also thank P. Koubský, R. Kříček, D. Korčáková as well as J. A. N., who obtained some of the spectra used in this study. We thank R. Ignace for providing us with his reduced spectra from the Ritter Observatory in phase space and for illuminating comments on his study. We also acknowledge the Ritter Observatory for making the archive of spectra public. M. B. and P. H. were supported by the Czech Science Foundation grant 19-01995S. J. B. was supported by the VEGA 2/0031/18 grant and by the Slovak Research and Development Agency under the contract No. APVV-15-0458. H. B. acknowledges financial support from the Croatian Science Foundation under the project 6212 “Solar and Stellar Variability”. The CHARA Array is supported by the National Science Foundation under Grant No. AST-1636624 and AST-1715788.

References

Appendix A: Figures for observation-specific models

The observation-specific models converged for a given dataset (LC, VIS, CLO, T3, SED, SPE, VAMP, or VPHI) are presented in Figs. A.1–A.8. Generally, the individual χ² contributions are smaller than those for the joint model (see Table 1). We also take these results as a confirmation that our model converges and is capable of fitting these datasets optimally.

Fig. A.1. Same as Fig. 7 for the observation-specific model.

Fig. A.2. Same as Fig. 8 for the observation-specific model.

Fig. A.3. Same as Fig. 9 for the observation-specific model.

Fig. A.4. Same as Fig. 10 for the observation-specific model.

Fig. A.5. Same as Fig. 11 for the observation-specific model.

Fig. A.6. Same as Fig. 12 for the observation-specific model.

Fig. A.7. Same as Fig. 13 for the observation-specific model.

Fig. A.8. Same as Fig. 14 for the observation-specific model.

All Tables

All Figures

	Fig. 1. Example of VEGA differential visibility measurement during the night of 27 Jun 2013 with the 140 m E2W2 baseline. Top: Hα line (continuum normalised to 1). Middle: amplitude of the differential visibility (normalised to 1 in the continuum). Bottom: phase of the differential visibility in degrees.
In the text

	Fig. 2. (u, v) coverage of differential visibility and phase measurements from the CHARA/VEGA instrument. A subset of Hα data used for modelling is shown (black), together with additional He I 6678 (green) and 7065 (blue) data. The orientation is v >  0 north and u >  0 east.
In the text

Fig. 3. Normalised spectra of β Lyr A for two out-of-eclipse phases (0.288 and 0.785). Observed data with uncertainties are indicated (blue), as well as synthetic data (yellow), residua (red), and large χ2 contributions (red circles). This 1st model was only fitted to these two spectra. The Hα emission profile and its EW are well described, but there are systematic differences for the He I 6678 line; the synthetic Si II 6347, 6371 and Ne I 6402 lines have low EW.

In the text

	Fig. 4. Normalised spectra for 11 phases covering the whole light curve, including primary and secondary eclipses. Colours in this figure, as well as in others, have the same meaning as in Fig. 3. The synthetic Hα profiles for the second model exhibit a variability similar to the observed ones, although at several phases there are systematic differences (both + and −).
In the text

	Fig. 5. Observed differential visibility amplitude |dV| versus wavelength λ, normalised to 1 in the continuum (blue), and its decrease across the Hα profile. Uncertainties of |dV| are also plotted (grey). Synthetic visibilities (yellow) are shown for the two “extreme” values of the shell’s outer radius Routsh = 40 R⊙ (top) and 120 R⊙ (bottom).
In the text

Fig. 6. χ2 convergence (red) for joint model; individual contributions (LC, VIS, CLO, T3, SED, SPE, VAMP, VPHI) are also indicated. The model successfully converges to a local minimum. Some datasets have a substantially larger number of observations, (i.e. effectively a larger weight). The χ2 values are different from Table 1 because the model was re-converged several times, and uncertainties of some datasets were modified.

In the text

	Fig. 7. Comparison of observed and synthetic phased light curves from FIR to FUV, computed for joint (“compromise”) model. The names of datasets are shown in the right column. The light curves were arbitrarily shifted in the vertical direction.
In the text

	Fig. 8. Squared visibility |V2| (shifted by dataset number) for different baselines B/λ (in cycles). The CHARA/MIRC observations are at the bottom, NPOI in the middle, and CHARA/VEGA at the top (blue); synthetic data for the joint model are plotted for comparison (yellow). Overall trends (|V2|(B/λ)) seem to be correctly described, although there are some systematics for small groups of data.
In the text

	Fig. 9. Similar to Fig. 8, for the closure phase argT3 vs. B/λ.
In the text

	Fig. 10. Similar to Fig. 8, for the triple product |T3| vs. B/λ.
In the text

	Fig. 11. Spectral-energy distribution (SED) overplotted for ten phases. Synthetic SEDs systematically differ from observations, especially in NUV and in the vicinity of the Balmer jump. This is not the best model in terms of ${\chi }_{\mathrm{sed}}^2$ (however, cf. Fig. A.5).
In the text

	Fig. 12. Normalised spectra I(λ) for 11 phases. Both observed (thick) and synthetic (thin) spectra are plotted for comparison. Colours correspond to the phases (see also, values on the left). Uncertainties are not plotted, but they are of the order of 0.01. The vertical lines correspond to (from left to right): Si II 6347, Si II 6371, Ne I 6402, Hα, and He I 6678 wavelengths.
In the text

	Fig. 13. Differential visibilities |dV| vs. λ normalised in continuum. Synthetic (yellow) visibilities for the joint model exhibit a similar decrease across the Hα profile to the observed ones (blue). However, some of the synthetic |dV|s are systematically lower in the blue wing and higher in the red wing.
In the text

	Fig. 14. Similar to Fig. 13, for the differential phases argdV vs. λ. On average, phase changes are comparable in both observed and synthetic data. There are some remaining phase slips that are also fitted for.
In the text

	Fig. 15. Continuum synthetic images of β Lyr A for joint model, computed for four monochromatic wavelengths (from top to bottom): 155 nm (FUV), 545 nm (V), 1630 nm (H), 4750 nm (FIR). The orbital phase is always 0.25. The appearance of the system changes substantially. Optically thin components (disc atmosphere, jet, shell) are not seen very well in continuum radiation.
In the text

Fig. 16. Line-profile synthetic images of β Lyr A computed for wavelength range of Hα (i.e. from 655.1 to 657.8 nm, at a fixed phase 0.25). Optically thin components are clearly visible. The disc atmosphere appears first, because its Keplerian velocities close to the inner rim produce the blue-shifted wing. The jet inclined towards the observer appears second, with high velocities being projected to the line of sight. Finally, there is the spherical shell, with relatively low velocities spanning the core of Hα, which obscures other small-scale structures. For λ > 6563 Å, all objects disappear in a reverse order. An animated version is avaiable online and at https://sirrah.troja.mff.cuni.cz/~mira/betalyr/.

In the text

	Fig. 17. χ2 convergence for observation-specific models, when starting from the best-fit joint model. After performing up to 103 iterations, substantial improvements were achieved for some datasets (LC, T3, SED, VAMP, etc.). It confirms systematic differences between observational datasets.
In the text

	Fig. 18. Continuum synthetic images for observation-specific models (datasets LC, VIS, CLO, T3, SED, SPE, VAMP, and VPHI) for the wavelength 545 nm (V). The apparent differences (e.g. the thickness of the disc, the appearance of the primary) demonstrate systematics between datasets. Alternately, some datasets (e.g. VPHI) do not constrain certain parameters.
In the text

Fig. 19. Stellar evolution of a binary corresponding to β Lyr A. The component radii R1, R2 vs. the primary mass M1 are plotted (top), with observed values indicated by full circles, together with the distance a and R1, R2 vs. time t (middle); see horizontal dashed lines, and surface abundances X, Y, AC12, AO16 vs. time t (bottom). At t = 19.613 × 106 yr, the mass ratio q corresponds to the observed one.

In the text

	Fig. 20. Synthetic spectra for solar abundances (yellow) and 10−2 lower abundance of C (grey); observed spectra (blue) are plotted for comparison. For the solar composition, C II 6578 and 6583 emission is too strong.
In the text

	Fig. 21. Non-rotated NLTE synthetic spectra computed for stellar atmospheres with He-rich (N(H = 0.4), N(He) = 0.6; black) and solar (orange) compositions. The effective temperature was 14 000 K in both cases. For the He-rich case, the profile exhibits a central absorption.
In the text

	Fig. 22. Differential visibility amplitude |dV|. for He I 6678 line. The joint model was used, but not the convergence of the respective He I VAMP dataset.
In the text

	Fig. 23. Differential visibility phase argdV for the He I line. Again, the joint model, but no convergence.
In the text

	Fig. A.1. Same as Fig. 7 for the observation-specific model.
In the text

	Fig. A.2. Same as Fig. 8 for the observation-specific model.
In the text

	Fig. A.3. Same as Fig. 9 for the observation-specific model.
In the text

	Fig. A.4. Same as Fig. 10 for the observation-specific model.
In the text

	Fig. A.5. Same as Fig. 11 for the observation-specific model.
In the text

	Fig. A.6. Same as Fig. 12 for the observation-specific model.
In the text

	Fig. A.7. Same as Fig. 13 for the observation-specific model.
In the text

	Fig. A.8. Same as Fig. 14 for the observation-specific model.
In the text