© ESO 2018

1. Introduction

Mass transfer between close binary components has a profound impact on their evolution. Early models of the mass transfer by Kippenhahn & Weigert (1967) have explained the Algol paradox (Crawford 1955) and demonstrated that about 80% of mass transferred during the whole process is exchanged in less than 10% of its total duration (see also De Greve 1986, for recent, more sophisticated modeling). The character and outcome of the process, that is whether the mass transfer is conservative or whether some matter and angular momentum escapes from the system and forms a common envelope around the whole system (Kuiper 1941; Paczynski 1976), depend strongly on the properties of the binary in question before the beginning of mass exchange and on the actual mechanism of the mass transfer.

Kippenhahn & Weigert (1967) introduced an initial classification of systems undergoing mass exchange, depending on whether the mass-losing component overflows the Roche limit during the core hydrogen burning (case A), or shell hydrogen burning (case B). Later, the term “case AB” was suggested to distinguish between cases when the mass exchange starts as case A close to the exhaustion of hydrogen in the core and continues as case B later during the process. All such systems are progenitors of Algols, that is systems in the later stage of the mass exchange when the mass ratio had already been reversed. Especially for the more massive Algols, the mass exchange appears to be non-conservative as found by van Rensbergen et al. (2006) from the distribution of the mass ratios among the observed Algols. The actual mechanisms of the mass and angular momentum loss from the system are not well established and modeled as yet. Bisikalo et al. (2000) proposed a purely hydrodynamical mechanism, which assumes that the gas flow, after encircling the gainer, hits the original, denser flow from the mass-losing star and is deviated above and below the orbital plane in the form of bipolar jets. On the other hand, van Rensbergen et al. (2008) and Deschamps et al. (2013) have proposed that a radiative interaction (a hot spot) between the flow, accretion disk, and the gainer may be responsible for the mass loss in the equatorial plane.

One way to discriminate between different scenarios is to carry out detailed studies of systems undergoing a phase of rapid mass exchange and deduce the true distribution and kinematics of the circumstellar gas for them. The bright, well-known binary β Lyr A (HD 174638, HR 7106, HIP 92420) with a steadily increasing orbital period of 12ḍ94, is an archetype of such systems. The history of its investigation is more than two hundred years long, and here we refer only to a subset of the more recent studies relevant to the topic of this paper. Detailed reviews of previous studies can be found in Sahade (1980) and Harmanec (2002). Unless a clear distinction is needed, we shall simply use the name β Lyr to denote β Lyr A in the rest of the text.

β Lyr is currently in a phase of rapid mass exchange, although its initial mass ratio has already been reversed and is now q = m_g/m_d ≃ 4.50 (Harmanec & Scholz 1993), with the mass-losing component (donor) being the less massive of the two (m_d ≃ 4.50M_⊙, and m_g ≃ 13.3M_⊙). Its spectral type is B6-8 II and the effective temperature T_eff,d= 13 300 K (Balachandran et al. 1986). It is losing mass via a Roche-lobe overflow toward its more massive partner (gainer). The conservation of angular momentum (or equivalently, the Coriolis acceleration in the non-inertial corotating frame) forces the gas flow to encircle the gainer and to form an accretion disk around it (Huang 1963; Wilson 1974; Hubený & Plavec 1991). Being composed from hot, and mainly ionized material, it is optically (and also geometrically) thick in the continuum and because it is observed nearly edge-on, it obscures the gainer completely. It occupies almost the whole critical Roche lobe around the gainer in the equatorial plane and was found to have the temperature of its rim of 7000–9000 K (e.g., Linnell 2000; Mennickent & Djuraševć 2013). The geometry of the optically thick bodies (the donor and the disk) was reconstructed from near-infrared interferometric observations by Zhao et al. (2008). The current rate of the mass transfer between the binary components is high, ≈2.1 × 10⁻⁵ M_⊙ yr⁻¹ for a conservative transfer (Harmanec & Scholz 1993), and ≈2.9 × 10⁻⁵ M_⊙ yr⁻¹ for a non-conservative one (De Greve & Linnell 1994; van Rensbergen & De Greve 2016), as deduced from the large observed secular change of the orbital period of ≈19 s yr⁻¹ (Harmanec & Scholz 1993; Ak et al. 2007).

The presence of a hot spot has been advocated by Lomax et al. (2012) and by Mennickent & Djuraševć (2013). The latter authors also postulated a second “bright spot” that should arise from a spiral arm that is formed within the disk. Some part of the gas flow is also deflected in the direction perpendicular to the accretion disk and forms a pair of jets, whose existence was first proposed by Harmanec et al. (1996), and was confirmed by Hoffman et al. (1998), Ak et al. (2007), Ignace et al. (2008), and Bonneau et al. (2011) via different types of observations. Observational evidence of the mass loss from the system has been presented by Umana et al. (2000, 2002), who resolved a circumbinary nebula surrounding the system in radio emission and found that it extends along the direction of the bipolar jets. An attempt to image the optically thin medium in Hα has been carried out by Schmitt et al. (2009), but their observations lacked the spatial resolution needed to separate the individual structures.

The evolution of β Lyr from the initial to the current evolutionary stage was modeled by several authors. The early conservative model of mass exchange by Ziolkowski (1976) was found unrealistic by Packet & De Greve (1979), because it would lead to a contact system. The latest evolutionary tracks by Mennickent & Djuraševć (2013) and by van Rensbergen & De Greve (2016) are in a good mutual agreement, probably thanks to the fact that the former is based on evolutionary models by van Rensbergen et al. (2008). The latter is based on improved evolutionary models including tides, and predicts that the system undergoes a case AB mass transfer and that it originated from a detached binary with initial masses of 10.35 M _⊙ (donor) and 7 M _⊙, and a period of 2.36 d and its current age is 2.63 × 10⁷yr.

The aim of the present study is to analyze and use the very rich series of visible and infrared spectro-interferometric observations covering the whole orbit of β Lyr, complemented by series of standard UBVR, near-infrared and far-UV photometric observations. All these observations are compared to the predictions of several working models of optically thick components of the system, focusing on the size, shape and physical properties of the opaque accretion disk surrounding the gainer. To this end, we use an improved version of the modeling tool SHELLSPEC by Budaj & Richards (2004) developed for binaries embedded in a 3D moving circumstellar environment.

A continuation of this study (to be published later) will use differential spectro-interferometry and analyses of selected emission-line profiles to investigate the probable distribution and kinematics of optically thin parts of the circumstellar matter within the system.

2. Observations

Throughout this paper, reduced Julian dates RJD = HJD−2 400 000 are used. The quadratic orbital ephemeris by Ak et al. (2007):$$\begin{array}{l}{T}_{\text{min}.\mathrm{I}}(\text{HJD})=2\hspace{0.17em}408\hspace{0.17em}247.968(15)+12.913779(16)\cdot E\\ \hspace{1em}+3.87265(369)\times {10}^{-6}\cdot E^2,\end{array}$$(1) is used to compute orbital phases of β Lyr A. Phase 0 corresponds to superior conjunction of the mass donor.

Our study is based on the following sets of dedicated spectro-interferometric observations and multicolor photometric observations (the details on the observations and their reductions being provided in Appendices A and B).

2.1 Interferometry

Three spectro-interferometric instruments took part in a twelve nights long observational campaign aimed at β Lyr in 2013. We also have at our disposal all previous interferometric observations as detailed below.

• Navy Precision Optical Interferometer (NPOI; Armstrong et al. 1998): These observations were carried out in 16 spectral channels spread over wavelength region Δλ = 562 − 861 nm with two triplets of telescopes.

• Michigan InfraRed Combiner (MIRC; Monnier et al. 2004, 2006): These observations were acquired with six telescopes in H-band split into eight channels. Earlier observations in four-telescope mode have already been analyzed by Zhao et al. (2008) and qualitatively compared to a working model. We note also that these observations were acquired before the instrument was equipped with photometric channels (Che et al. 2010).

• Visible spEctroGraph and polArimeter (VEGA; Mourard et al. 2009; Mourard et al. 2011): These observations were taken in four spectral regions using medium spectral resolution R = 5000. In each of these regions two channels in continuum ≈10 − 15 nm wide were chosen. Either two or three telescopes were used.

The VEGA and MIRC instruments are mounted at Center for High Angular Resolution Astronomy (CHARA) interferometric array (Ten Brummelaar et al. 2005). On nine nights during the campaign, the MIRC and VEGA instruments were co-phased to record fringes simultaneously in the visible and infrared. Basic properties of the spectro-interferometric observations are listed in Table 1. The phase coverage of β Lyr A orbit with spectro-interferometric observations from individual instruments is shown in Fig. 1. A detailed overview of the observations is listed in Table A.1. The reduced quantities (visibilities, closure phases, and – for MIRC only – triple product amplitudes) are available at CDS in the OIFITS format (Pauls et al. 2005).

Fig. 1. Phase coverage of spectro-interferometric observations of β Lyr acquired by different instruments. δ denotes the relative declination (positive toward the north), and α the relative right ascension (positive toward the east). The black line shows the size and orientation of the β Lyr orbit in the sky, the blue dots show orbital phases corresponding to NPOI observations, the magenta dots to CHARA/VEGA observations, the green dots to CHARA/MIRC observations acquired in 2013, and the red dots to CHARA/MIRC observations acquired in 2006/2007. An arbitrary vertical shift of 0.2 mas is added to separate the various orbits.

2.2. Photometry

A new series of differential Johnson UBVR photometric observations were acquired at Hvar Observatory in 2013–2017, with earlier infrared photometry acquired by Jameson & Longmore (1976), and Taranova & Shenavrin (2005). Apart from that, we included NUV and FUV observations from OAO A-2, IUE, and Voyager satellites described in Kondo et al. (1994).

Moreover, there is a Johnson–Cousins differential BV(R)_C photometry acquired by PS at his private observatory in Brno, Czech Republic. The latter observations were not modeled, though, because of their limited phase coverage and a relatively simple standardization procedure. Therefore these observations served only as an independent check that the Hvar UBVR measurements do not miss an important light curve feature.

Journal of photometric observations is in Table 2. The UBVR observations acquired at Hvar observatory and BV(R)_C measurements collected by PS are available in Table B.1 at CDS.

3. Choosing the initial physical properties for detailed modeling

3.1. Distance estimates

3.1.1. Trigonometric parallax

Perryman & ESA (1997) published the HIPPARCOS parallax of β Lyr 0″̣00370±0″̣00052. van Leeuwen (2007a,b) carried out a new reduction of HIPPARCOS data to obtain a more accurate value of 0″̣00339±0″̣00017. These values translate to the following distances estimates and 1 − σ ranges:

270 pc; range 237–314 pc,

295 pc; range 281–311 pc,

for the original and improved HIPPARCOS parallax, respectively.

3.1.2 Dynamical parallax from the orbital solution and spatial resolution of the orbit

Zhao et al. (2008) used two different techniques of image reconstruction and a simple model with two uniformly illuminated ellipsoids and derived three different distance estimates from the model and two methods of image reconstruction, respectively:

314 pc; range 297–331 pc,

278 pc; range 254–302 pc,

274 pc; range 240–308 pc.

It is important to realize that these distances are based on the value of the projected value of the semimajor axis a sin i = (57.87 ± 0.62) R _⊙ defined by the Kepler’s Third Law and the binary masses, which were estimated from several previous studies. In this study, we shall use a similar approach to provide an independent distance estimate based on our new models of the continuum radiation of β Lyr.

3.1.3. β Lyr B as a distance indicator of β Lyr A

β Lyr B (HD 174664, BD+33°3224, STFA 39B) is the second brightest member of the β Lyr visual system. It is a main-sequence B5 V star. A somewhat puzzling is the fact that it is also an X-ray (Berghofer & Schmitt 1994). Abt et al. (1962) measured its RV on eleven blue photographic spectra and concluded that the star is a single-line spectroscopic binary with an orbital period of 43ḍ48 and an eccentric orbit. They also studied available astrometric observations for all visual members of the β Lyrs̃ystem and concluded that β Lyr B is gravitationally bound to β Lyr A. Abt & Levy (1976) obtained new photographic RVs of β Lyr B and carried out an error analysis of newly obtained and earlier RVs to conclude that there is no evidence for duplicity of β Lyr B. From spectral classification they concluded that β Lyr B is a normal B7V star close to the zero-age main sequence.

Here, we critically re-investigate these pieces of information to see whether β Lyr B can be used to another distance estimate of β Lyr A. We first study the kinematic and radiative properties of β Lyr B. Radial velocities (RVs) were measured on Ondřejov Reticon red spectra of β Lyr B secured from 1995 to 1996, and on six Ondřejov red CCD spectra from 2003 to 2016. The linear dispersion of all these spectra was 17.2 Å mm⁻¹, and their two-pixel resolution was 12 700. Additional details on the reduction of β Lyr B spectra are presented in Appendix C.

RVs of β Lyr B were determined using two methods:

1. An interactive comparison of direct and flipped line profiles on the computer screen until the best match is achieved in program SPEFO (Horn et al. 1996; Škoda 1996). RVs were measured independently on four spectral lines Si II 6347 Å, S II 6371 Å, Hα, and He I 6678 Å.

2. Via an automatic comparison of the observed and synthetic spectra in the PYTERPOL program (Nemravová et al. 2016)¹.

The latter method produced slightly less scattered RVs, hence only these are presented here. However, the RVs obtained by both methods are listed in Table C.1. At the same time, PYTERPOL also estimated radiative properties of β Lyr A : T_eff, log g, and the projected rotational velocity v sin i. The red spectra we used (Δλ ≃ 6200–6800 Å) contain numerous telluric lines, which would adversely affect the results of PYTERPOL. We handle the problem using two methods: (i) 1 Å intervals centered on each telluric spectral line were omitted from each spectrum, and (ii) telluric and stellar spectra were separated by spectral disentangling (Simon & Sturm 1994) in Fourier space by the KOREL program (Hadrava 1995, 1997). PYTERPOL was then used for both, individual observed spectra with telluric-line regions omitted, and to disentangled stellar spectrum, free of telluric lines. The synthetic spectra were taken from the BSTAR (Lanz & Hubený 2007) and AMBRE (Palacios et al. 2010) grids. The results were similar, though not exactly the same. We simply adopted their mean as a realistic estimate of radiative properties of β Lyr B, noting that both methods may introduce some systematic errors. The observed spectra still contained remnants of weaker telluric lines, while the disentangled spectra were slightly warped and had to be re-normalized. The formal errors derived from the two solutions are unrealistically small, as they do not reflect possible systematic errors. The results are summarized in Table 3 and a comparison of one observed profile, and disentangled profiles with the best-fitting synthetic spectra is shown in Fig. C.1. We assumed the solar composition in these fits guided by the results of Abt&Levy (1976) from the blue spectra since the red spectra at our disposal do not contain enough spectral lines to estimate the metallicity of β Lyr B reliably.

Do β Lyr A and β Lyr B indeed form a physical bounded system? All RV measurements of β Lyr B are plotted in Fig. 2. The rms errors are only available for Abt & Levy (1976) and for our new RVs and are shown in this plot. We explain the large scatter of the RVs from the photographic plates by the combination of four factors: moderate dispersion of the plates, variable quality of individual spectra, relatively large projected rotational velocity of β Lyr B, and line blending. It is our experience that due to blends with fainter spectral lines the apparent RVs of individual measured lines can differ systematically. Since different number of spectral lines were measured on different plates, their average RVs are then affected differently. This effect should be absent in PYTERPOL RVs, based on the comparison of whole spectral segments with interpolated synthetic spectra. Also our manual RV measurements in SPEFO should be basically free from these effects since one sees which part of the profile is measured. To check on this, we derived robust mean RVs for the four stronger lines seen in the Ondřejov spectra (using the algorithm published by Andrews 1972). They are summarized in Table 4 and confirm a very good agreement of all mean RVs of individual lines. The systemic velocity of β Lyr B based on PYTERPOL RVs is γ _B = −(18.1 ± 0.5) km s⁻¹. This number is a weighted mean of two estimates. The first one is the mean RV derived from RV measurements on individual spectra, and the latter is measured on the disentangled spectrum (both are presented in Table 3). For SPEFO RVs from individual lines one gets a 1-σ range from −17 to −21 km s⁻¹. These estimates are to be compared to the systemic RV of β Lyr A, γ _A, which was found to be in the range from about −17 to −20 km s⁻¹, depending on the spectrograph used (see for example, solution 5 in Table 10 of Harmanec et al. 1996). From this, one can conclude that the systemic RVs of β Lyr A and β Lyr B are identical within the range of their measuring errors. In spite of their large scatter, the RVs of β Lyr B measured on old photographic spectra by Petrie & Pearce (1961) and Abt et al. (1962), Abt & Levy (1976) do not show any obvious long-term RV trend.

Fig. 2. Available RVs of β Lyr B. The first three datasets are RVs measured on photographic spectra of moderate dispersion: green open circles denote the RVs measured by Petrie & Pearce (1961), the red open triangles come from Abt et al. (1962), and the black boxes are RVs from Abt & Levy (1976). The filled blue diamonds are the RVs measured on Ondřejov electronic spectra with PYTERPOL. The rms errors, available only for Abt & Levy (1976) and new RVs are also shown. We note that the most deviant RV of −45 km s−1 is based on a single line and comes probably from an underexposed spectrum.

The range of individual RVs based on our measurements with their corresponding rms errors could lead to the suspicion that β Lyr B could be a spectroscopic binary after all. To check on this, we run period searches down to 0ḍ5 for PYTERPOL RVs and also SPEFO RVs for individual lines. All the periodograms show a dense forest of comparably deep peaks but the frequencies found differ mutually, both for the individual lines and for the PYTERPOL RVs. From this, we reinforce the conclusion of Abt & Levy (1976) that β Lyr B is a single star.

The proper motions of β Lyr A and β Lyr B were also investigated. Measurements accessible through VizieR database (Ochsenbein et al. 2000) at CDS were downloaded. Only values of μ_α and μ_δ published after the HIPPARCOS mission were retained. Proper motions in the right ascension were corrected for the declination of both systems (see Appendix C.3 for details). Their weighted mean is given in Table 5. Individual measurements and the weighted average are shown in Fig. 3.

Fig. 3. Proper motions of β Lyr A and β Lyr B. The black points denote proper motion measurements of β Lyr A, the black triangles proper motions of β Lyr B, the red point the mean weighted proper motion of β Lyr A, and the red triangle the mean weighted proper motion of β Lyr B. Xµ, Yµ are Cartesian coordinates of the proper motion vector given by Eqs. (C.1) and (C.2). The former is north–south oriented, and the latter east–west oriented.

The close similarity of systemic RVs of β Lyr A and β Lyr B reinforce the hypothesis that the two systems formed in the same association or are even physically linked. The same cannot be said for the proper motions or for the respective tangential velocities, which differ from each other quite significantly. We conjecture that a part of this difference can be attributed to the mutual orbit – if β Lyr A and B are bound – or to the intrinsic velocity dispersion σ_v of the putative stellar association. If the A − B orbit is circular and seen more or less edge on, the Keplerian velocity $v_{\text{kepl}}=\sqrt{G{M}_{\text{tot}}/(\alpha d)}$ is on the order of 1.2 km s⁻¹, assuming the mass M_tot ≃ 20 M_⊙, the angular separation α ≃ 40 arcsec, and the distance d ≃ 325 pc. This should be compared with the relative tangential velocity v_t = |µ_A – µ_B| d ≃ 4.2 km s⁻¹. This value seems larger than v_kepl and might be on the order of typical σ_v.

Nonetheless, the spread of proper motion measurements of β Lyr B is large and the weighed mean is dominated by the latest Gaia observation. If proper motion measurements of β Lyr B are given equal weight, the mean proper motion is: µ_B = 3.4 ± 2.4 mas yr⁻¹, θ_B = 151 ± 43 deg, that is in agreement with measurements of β Lyr A. Therefore, we think that the proper motions of both systems do not contradict the conclusion that the systems formed within the same association. That alone allows us to use β Lyr B as a distance indicator for β Lyr A.

Gaia Collaboration (2016a) published the first parallax derived by the Gaia satellite (see Gaia Collaboration 2016b) for β Lyr B, 0″̣00310 ± 0″̣00036, which translates to a distance d = 323 pc; range 289–365 pc.

In April 2018, the second DR2 release of the Gaia catalog was made publicly available. The improved value of the β Lyr B parallax is 0″̣003006 ± 0″̣000054, which translates to a rather narrow distance range

d = 333 pc; range 327–339 pc.

We have also derived a spectroscopic distance of β Lyr B using its radiative properties (see Table 3, third Col.), and the mean all-sky UBV magnitudes based on 77 observations acquired at Hvar observatory during three seasons. The mean all-sky magnitudes:$$V=7\stackrel{\text{m}}{.}199,B-V=-0\stackrel{\text{m}}{.}089,U-B=-0\stackrel{\text{m}}{.}541$$were dereddened in a standard way, which resulted in:$$V_0=6\stackrel{\text{m}}{.}99,E(B-V)=0\stackrel{\text{m}}{.}070.$$Assuming that β Lyr B is a main-sequence star, its mass and radius can be estimated from relations of Harmanec (1988) to be $M_{\text{B}}={4.16}_{-0.09}^{+0.08}{M}_{\odot }$, and $R_{\text{B}}={2.53}_{-0.16}^{+0.18}{R}_{\odot }$, respectively. Adopting bolometric corrections of Popper (1980), the absolute magnitude of β Lyr B is $M_{\text{V}}=-{0.08}_{-0.16}^{+0.17}$ mag, and its spectroscopic distance is:$$d_{\text{B}}^{\text{Sp}.}=(259\pm 20)\text{pc}.$$This is a significantly lower value than what Gaia obtained but it is actually in qualitative agreement with recent findings by Stassun & Torres (2016) that Gaia DR1 parallaxes are smaller than photometrical ones for well studied eclipsing binaries. Their finding was confirmed also by Graczyk et al. (2017). Since we note that in many cases the DR2 parallaxes do agree with the DR1 ones within the quoted errors (significantly smaller in DR2), the above warning might also be relevant for the DR2 parallaxes.

As pointed out by the anonymous referee, possible duplicity of β Lyr B could increase its observed luminosity and, therefore, its distance from us. Without a direct imaging, one cannot exclude indeed the possibility that β Lyr B is composed of two similar B7V stars seen pole on. In that case the spectroscopic distance of such an object could be as large as ∼370 pc. While it would certainly be desirable to obtain new, high-dispersion and high-S/N spectra of β Lyr B and derive its more precise RVs to check on their constancy, we do not find the duplicity of β Lyr B as too probable. Given the fact that the object was detected as an X-ray source, it is more probable that any putative secondary would be a late-type star with a chromospheric emission. However, no lines of a cool secondary have ever been detected so such an object – if present – must be much fainter than the primary. If so, its presence would not affect our spectroscopic distance estimate significantly.

3.2. The mass ratio and orbital solutions

Because the RV curve of the star hidden in the disk is only defined by a pair of fainter S II absorption lines at 6347 and 6371 Å discovered by Sahade (1966) and Skulskii (1975), one should also consider some range of plausible mass ratios. Various relevant estimates are summarized in Table 6.

For the purpose of our modeling, we shall use the orbital solution of Harmanec & Scholz (1993) based on Si II lines for both stars (their solution 6 of Table 10), which was basically confirmed by more extended series of spectra by Harmanec et al. (1996). In particular, we shall adopt K₁ = 41.4 ± 1.3, K₂ = 186.30 ± 0.35, and a sin i = 58.19 R_⊙, which implies the mass ratio of 4.500.

4. Modeling the continuum radiation

Our modeling of β Lyr A and its application to spectro-interferometric and photometric observations is presented in the following section. Our models are based on the program SHELLSPEC, developed by Budaj & Richards (2004). This program was equipped with additional features that simplify modeling of binaries, computations of synthetic interferometric observables, and a solution of the inverse task.

First, we briefly outline the program and its new Python wrapper called Pyshellspec. Then, we proceed to the development of several alternative models of β Lyr in continuum and their comparison to observations.

4.1. Foundations of the model

The model is based on the existing program SHELLSPEC which was designed for the computation of synthetic light curves, spectra and images of stars and/or binaries surrounded by a moving 3D circumstellar medium by means of solving a one-dimensional LTE radiation transfer along some user-specified line-of-sight.

To provide an environment for solving the transfer, the modeled (non-stellar) objects cannot be represented with a 2D mesh covering only the photosphere of each object. Instead the grid has to sample the whole volume in 3D. The only two exceptions are models of stars, which are always opaque, hence their atmospheres form a boundary condition. In SHELLSPEC, each object is placed into an regularly sampled cuboid (divided into n_x × n_y × n_z cells) that represents “the Universe”. All embedded objects inherit the spatial sampling of the Universe. Every object has a simple geometric shape given by several parameters. Each cell occupied by an object is assigned its density (gas, electron, and dust), temperature and velocity.

We use the following setup: LTE level populations, LTE ionisation equilibrium, the line profile is determined by thermal, microturbulent, natural, Stark, Van der Waals broadenings, and the Doppler shift. The continuum opacity is caused by HI bound-free, HI free-free, H⁻ bound-free, and H⁻ free-free transitions. Moreover, we account for the line opacity of Hα, He I 6678, and He I 7065, for the future spectral line analysis of the VEGA data. Abundances are assumed to be solar. We use a small grid of synthetic spectra for the stars, 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, gravity darkening (in particular the Roche-filling donor), and the reflection effect concerning the heat redistribution over the surfaces.

On the other hand, we did not include the Thomson scattering on free electrons, the Rayleigh scattering on neutral hydrogen, because these are only implemented as optically thin (single) scattering in SHELLSPEC. It would be a much harder computational problem to account for multiple scattering in 3D moving optically thick medium, and would essentially prevent us to converge the model with many parameters. There is also neither irradiation nor reflection between the stars. As the disk is presumably hotter than the silicate condensation temperature, we account for neither Mie absorption on dust, Mie scattering, nor dust thermal emission.

The velocity field of an object is given either by a net velocity or by rotational velocity, or a combination of both. The observable quantities (flux and intensity) are computed for a line-of-sight, specified by two angles. The model is only kinematic, that is the radiation field has no effect on the state quantities of the circumstellar medium.

SHELLSPEC allows modeling of various structures, but we restrict ourselves to those that are relevant for the continuum model of β Lyr A. A very detailed description of SHELLSPEC is in Budaj & Richards (2004), and its latest improvements are described in Budaj (2011a,b).

4.2. Interface for automatic comparison and fitting

On output, SHELLSPEC computes the monochromatic flux F_v (in erg s⁻¹ cm⁻² Hz⁻¹), and intensity I_v(x, y) (in erg s⁻¹ cm⁻² Hz⁻¹) projected onto the plane perpendicular to the line-of-sight, where x and y are the Cartesian coordinates in this plane, but it does not carry out a comparison of these quantities to observed ones, or automatic optimization of the model parameters. Moreover, individual model components implemented in SHELLSPEC are not bound by any orbit, although some model parameters depend on orbital parameters. The components are almost independent and do not share parameters that are common. To overcome these limitations we wrapped SHELLSPEC in a Python interface².

4.2.1. Computation of synthetic magnitudes, squared visibilities and closure phases

In our approach, the passband flux (e.g., F_V) is computed from the monochromatic flux F_λ = F_vc/λ² for a single (effective) wavelength λ_eff, given by the transmission curve of the respective filter. Its relation to the passband magnitude is simply m_V = –2.5 log10 F_λ/F_calib, where F_calib denotes the calibration flux (also monochromatic), given for example by Johnson (1966), Kondo et al. (1994). This approximation is required because radiation transfer computations for many wavelengths would be very time consuming. Unfortunately, this may lead to a slight offset Δ_P between the observed and synthetic light curves. The offset is thus determined by a minimization of the following formula:$$S({\mathrm{\Delta }}_{\text{P}})={\displaystyle \sum _{i=1}^N\left|F_i^{\text{obs}}-F_i^{\text{syn}}+{\mathrm{\Delta }}_{\text{P}}\right|},$$(2) where F^obs is the observed flux, and F^syn the synthetic flux.

The FWHM of the passbands Δλ_eff of interferometric observations were as low as 20 nm in the visible and 50 nm in the infrared. Therefore the passband intensity I_p was calculated in a similar way from the monochromatic intensity I_v for the effective wavelength. This approximation is validated by the fact that the continuum does not change significantly through the narrow passbands. Images are normalized afterward, rendering any offset with respect to the emergent intensity insignificant.

The images produced by SHELLSPEC are centered on the primary component (the gainer), and their nodal line is aligned with the north–south direction. Hence, the image center has to be shifted to the system barycenter first and then rotated to a given longitude of the ascending node Ω. The complex visibility V(u, v) is computed by a two-dimensional discrete Fourier transform (DFT) of the image at spatial frequency b = (u, v) = (Bx, By)/λ, where Bx(By) denotes the projection of the baseline into east–west (north–south) direction. The triple product T3 is then computed as follows:$$T_3(b_1,b_2)=V(b_1)V(b_2)V^{*}(b_1+b_2),$$(3)where b₁ and b₂ denote spatial frequencies corresponding to a pair of baselines in a closing triangle. Performance of the two-dimensional fast Fourier transform (FFT) for computation of interferometric observables was also evaluated. The need for an extensive zero-padding of each image produced by SHELLSPEC to obtain a sufficient resolution, and interpolation within the two-dimensional array made it actually slower than DFT.

4.2.2. Interface for the solution of the inverse task

The optimization of model parameters was carried out through the minimization of the total χ² defined as the sum of the χ² of the different data sets. MIRC providing triple product quantities that are not totally independent from the V² estimates, we decided to use a specific weight on the MIRC V² and T₃ data.$${\chi }^2={\chi }_{\text{LC}}^2+{\chi }_{\text{IF}}^2,$$(4)where the contributions of individual types of observations are given by the following relations:$${\chi }_{\text{LC}}^2={\displaystyle \sum _{i=1}^{N_{\text{P}}}{\displaystyle \sum _{j=1}^{N_{\text{M}}}{\left(\frac{m_{i,j}^{\text{obs}}-{\tilde{m}}_{i,j}^{\text{syn}}}{{\sigma }_{i,j}}\right)}^2}},$$(5) $${\chi }_{\text{IF}}^2={\chi }_{{\text{IF}}_{\text{VEGA}}}^2+{\chi }_{{\text{IF}}_{\text{NPOI}}}^2+{\chi }_{{\text{IF}}_{\text{MIRC}}}^2,$$(6)with$$\begin{array}{c}{\chi }_{{\text{IF}}_{\text{VEGA}}}^2={\chi }_{{\text{V}}^{\text{2}}}^2,\\ {\chi }_{{\text{IF}}_{\text{NPOI}}}^2={\chi }_{{\text{V}}^{\text{2}}}^2+{\chi }_{\text{CP}}^2,\\ {\chi }_{{\text{IF}}_{\text{MIRC}}}^2=\frac{1}{2}\left({\chi }_{{\text{V}}^{\text{2}}}^2+{\chi }_{{\text{T}}_{\text{3}}}^2\right)+{\chi }_{\text{CP}}^2,\\ {\chi }_{{\text{V}}^{\text{2}}}^2={\displaystyle \sum _{i=1}^{N_{{\text{V}}^{\text{2}}}}{\left(\frac{{\left|V^{\text{obs}}\right|}_i^2-{\left|V^{\text{syn}}\right|}_i^2}{{\sigma }_i}\right)}^2},\\ {\chi }_{{\text{T}}_{\text{3}}}^2={\displaystyle \sum _{i=1}^{N_{{\text{T}}_{\text{3}}}}{\left(\frac{{\left|T_3^{\text{obs}}\right|}_i-{\left|T_3^{\text{syn}}\right|}_i}{{\sigma }_i}\right)}^2},\\ {\chi }_{\text{CP}}^2={\displaystyle \sum _{i=1}^{N_{{\text{T}}_{\text{3}}}}{\left(\frac{T_3{\varphi }_i^{\text{obs}}-T_3{\varphi }_i^{\text{syn}}}{{\sigma }_i}\right)}^2}.\end{array}$$(7) where m denotes the magnitude, m̃ the magnitude corrected for the offset given by Eq. (2), N_M the number of photometric observations for jth passband, N_p the number of fitted passbands, |V|² the squared visibility, |T₃| the modulus of the triple product, T₃ϕ ≡ arg{T₃} the closure phase, N_V² the number of squared visibility observations, N_T³ the number of triple product observations, and σ’s are the uncertainties of the corresponding observations. MIRC provides triple product quantities that are not entirely independent from the V² estimates and we thus decided to use a specific weight on the MIRC V² and T₃ data. No additional weighting of the data sets is considered but a detailed analysis of the convergence process is presented in Table 8.

Available engines for the global minimization of Eq. (4) are the differential evolution algorithm by Storn & Price (1997), implemented within the SciPy library, and the Simplex algorithm by Nelder & Mead (1965), implemented within the NLOPT library.

As the properties of some of the objects are linked to the orbital elements, we implemented an orbit binding the objects together for the global optimization. The orbit is given by the following elements: the period P at a reference epoch, the epoch of primary minimum T_min, the rate of the period increase Ṗ, the eccentricity e, the semimajor axis a, the mass ratio q, the inclination i, the argument of periastron ω, and the longitude of ascending node Ω. The orbit binds the two stars and supplies values of orbit-dependent parameters (the masses, Roche-lobe radii, distance between the primary and secondary, radial velocities, and the orbital phase for a given epoch). With this in hand, we then use SHELLSPEC to compute the model at different phases.

4.3. A model for β Lyr A

A model for individual components of β Lyr A is introduced here. While the model of the two binary components is straightforward, properties of the accretion disk surrounding the gainer remain uncertain. Hence several models of the accretion disk were constructed and tested.

4.3.1. A model of binary components

The two stellar components of β Lyr A were approximated with the following models:

• The donor has a Roche geometry and is completely filling the Roche limit and its rotation is synchronized with the orbital motion. Its shape and size are driven by the semimajor axis a, and the mass ratio q. The gravity darkening (von Zeipel law, α _GD = 0.25) and a linear limb darkening were assumed. Coefficients of the limb-darkening law were taken from van Hamme (1993) using the tri-linear interpolation scheme. For the interpolation, the polar temperature was used instead of the effective temperature, gravitational acceleration was approximated with that of a sphere with radius equal to the polar radius, and the solar metallicity was used. Parameters describing the donor are denoted by the index “d”.

• The gainer is approximated with a sphere, even though it likely rotates close to its critical velocity, v _crit = v_kepl(Rg), because of the ongoing accretion, and should thus have an ellipsoidal shape and significant gravity darkening. Alternatively, there can be a thin transition layer (Pringle 1981), if the gainer is not (yet) rotating critically. We also checked that the accretion rate Ṁ is low enough and that the associated radial velocities within the disk are much smaller than Keplerian, v_r ≪ v_kepl. Nevertheless, due to the presence of the accretion disk, the only radiation that may be able to penetrate the disk is that coming from polar regions, whose radiation and shape do not depart from a spherical star so much. The component is limb-darkened using the same tables and interpolation scheme as for donor. The radius was set to a value typical for B0.5 IV–V star, R _g = 6 R_⊙ (Harmanec 1988), which is also in agreement with radii adopted in earlier studies (e.g., Harmanec 1992; Linnell 2000). Parameters describing the gainer are denoted by the index “g”.

4.3.2. Models of the accretion disk

SHELLSPEC allows the user to set up several models of an accretion disk. They differ in shape, radial temperature, and density profiles as presented below. In Sect. 5.2, we give more details on the preferred geometries.

The disk plane of each model lies in the orbital plane of β Lyr; z–axis is perpendicular to this plane and goes through the center of each disk. Except for the envelope-shaped disk, the z–axis is also the axis of symmetry. The radius R is measured in the disk plane from the center of the disk. The shapes of accretion disks that we tested are (see also Fig. 4):

1. Slab (panel I in Fig. 4) is limited by two spherical surfaces with radii R_in, R_out, and two surfaces z = ±H. The vertical temperature profile T(z) as well as the density ρ(z) are constant. The velocity profile v(r) is Keplerian as in all other cases.

2. Wedge (panel II) is limited by two spherical surfaces with radii R_in, R_out, and two conical surfaces z = ±R sin ϑ, where ϑ is the half-opening angle of the accretion disk. H denotes the maximal height of the disk at its outer rim.

3. Lens (panel III) is limited by a spherical surface with radius R_in and an ellipsoidal surface whose semimajor axis is equal to R_out and semiminor axis to H.

4. Envelope (panel IV) is limited by a Roche equipotential whose shape is given by the semimajor axis a of the orbit, the mass ratio q, and the filling factor f_f. A synchronization between rotation and orbital motion is assumed. The vertical structure of the envelope is further limited by two surfaces z = ±H. The envelope can only be homogeneous and isothermal (ρ = const., T = const.).

5. Nebula is a standard disk with a Gaussian vertical density profile ρ(z) determined by the hydrostatic equilibrium; in other words, its scale height H is determined by the local temperature T(r). Additional parameters (h_inv, t_inv, h_wind) can be used to account for a temperature inversion in the disk atmosphere, or a non-zero constant density in the wind region. However, in order to prevent discretization artifacts (discussed separately in Appendix D), we assume the temperature inversion is not a sudden jump from T to t_invT, and in our SHELLSPEC T(z) changes linearly between h_invH and the outer limit. In order to account also for a non-hydrostatic disks we modified SHELLSPEC to include a multiplicative factor h_mul, which can be treated as a free parameter too.

Fig. 4. Geometric shapes of models of accretion disks. (I) Denotes the slab, (II) the wedge, (III) the lens, and (IV) the envelope. The latter is shown as if it was filling the Roche limit, that is Rout = RL1, where Rout = RL1 denotes point radius of the Roche lobe. Here it is shown this way to emphasize the geometry of the disk, but its shape may be more symmetric depending on its filling factor. See Sect. 4.3 for the description of individual disk shapes. The Nebula is not presented here, since no simple geometrical sketch can be given. It is illustrated, however, in Fig. 5.

Two types of radial temperature profiles T(R) were tested. The first one is a power-law:$$T(R) = T_0{\left(\frac{R}{R_{\text{in}}}\right)}^{{\alpha }_{\text{T}}},$$(8)where T₀ is the effective temperature at the inner rim of the accretion disk R_in, and α_T the exponent of the power-law. The second one corresponds to a steady accretion disk heated only by viscous dissipation (Shakura & Sunyaev 1973; Pringle 1981):$$T(R) = T_1{\left(\frac{R_{\text{g}}}{R}\right)}^{\frac{3}{4}}{\left(1 - \sqrt{\frac{R_{\text{g}}}{R}}\right)}^{\frac{1}{4}},$$(9) where T₁ is the characteristic temperature of the disk, which is actually never attained in the disk. The maximum temperature T_max = 0.488 T₁ is at the radius 49/36 R_g. The temperature along z–axis is constant.

The radial density profile of gas ρ(R) is always approximated by a power law:$$\rho \left(R\right) = {\rho }_0{\left(\frac{R}{R_{\text{in}}}\right)}^{{\alpha }_{\text{D}}},$$(10) where ρ₀ is the gas density at the inner rim of the accretion disk. The temperature in the whole accretion disk is too high ≳7000 K for condensation of dust grains, so we assumed that the dust density is zero. The electron density was computed assuming LTE and solar chemical composition. We simply assume the disk is stable over the time span of our observations, that is 1970–2017.

4.4. Realistic uncertainties of observational data

Our model represents a rather simplified view on β Lyr A, so it is no surprise that a straight preliminary comparison of the model and data gave reduced ${\chi }_R^2\gtrsim 20$. This would mean that our model is wrong, but our model is not very different from models presented in earlier studies of the system. We attributed this mismatch to systematic errors and we tried to compensate for some of them by taking the following steps:

• The observed light curves of β Lyr A, spread over more than 40 yr, exhibit small bumps (e.g., asymmetry near primary and secondary eclipses), or flickering (bottom of the primary minimum). They represent an intrinsic variability of β Lyr A that is beyond the capabilities of our model. The limited resolution of the model creates jags in the synthetic light curve. They are most pronounced at centers of both eclipses. The latter effect is higher, and the systematic signal was in general ≲0.05 mag (that is 3–5 times higher than true uncertainty of these measurements). Therefore the uncertainty of Hvar UVBR light curves was set to this value.

• The infrared light curves lacked uncertainty estimates. They were estimated by inspecting the scatter within intervals 0.05 wide in phase. This uncertainty estimate was highly phase-dependent, but their mean was ≃0.1 mag, hence we adopted this uncertainty for all infrared light curves.

• The squared visibility |V|² and closure phase T₃ϕ are commonly affected by systematic effects given by atmospheric fluctuation during observations. The CHARA/VEGA squared visibility observations obtained with one baseline over a short period of time (one block, ≈20 min) were frequently showing much larger spread than uncertainty of individual points. This spread is unlikely to be a product of the slow change of the projected baseline caused by the diurnal motion or produced by the intrinsic variability of β Lyr. Also our experience is that measurements of low squared visibility ≲0.05 have floor uncertainty ≈0.05. The uncertainties of all CHARA/VEGA measurements of |V|² were adjusted according to the following formula:$${\sigma }_{V^2}^{\text{new}}=(\begin{array}{ll}\text{max}\{{\sigma }_{V^2}^{\text{old}},{\sigma }_{V^2}^{\text{block}}\}\hfill & \text{for}\hspace{1em}V{|}^2>0.05,\hfill \\ \text{max}\{{\sigma }_{V^2}^{\text{old}},{\sigma }_{V^2}^{\text{block}},0.05\}\hfill & \text{for}\hspace{1em}V{|}^2<0.05,\hfill \end{array}$$(11)where ${\sigma }_{V^2}^{\text{old}}$ is the uncertainty estimated with the reduction pipeline (see Mourard et al. 2009, and references therein), and ${\sigma }_{V^2}^{\text{block}}$ is the standard deviation of points within one block of observations (see Appendix A for details).

• For CHARA/MIRC observations, the uncertainty of closure phase measurements that satisfy the condition S / N_|T3| ≈ 1 was adjusted following the formula by Monnier et al. (2012):$${\sigma }_{\text{CP}}^{\text{new}}=\text{max}\left\{{\sigma }_{\text{CP}}^{\text{old}},30\text{ deg}/S/N_{\left|{\text{T}}_3\right|},\frac{1}{5}{\left(\text{Δ}{T}_3\varphi \right)}_{\lambda }\right\},$$(12) where ${\sigma }_{\text{CP}}^{\text{old}}$ is the original value estimated by the reduction pipeline (see Monnier et al. 2004 and references therein), S / N_|T3| is the signal-to-noise ratio of the corresponding triple amplitude measurement, and (ΔT₃ϕ)_λ is the difference between the highest and the lowest closure phase measurement in a single passband for one block of CHARA/MIRC observations (see Appendix A for details). Similarly to Zhao et al. (2011) minimal uncertainty 1 deg was adopted for all CHARA/MIRC closure phase measurements. The uncertainties on the triple product amplitudes are also corrected to account for systematic effects as described by Monnier et al. (2012): we use an additive error of 1 × 10⁻⁵ and a multiplicative factor of 10%.

• A systematic offset was apparent for triple product amplitudes |T₃| obtained by the NPOI instrument. While the squared visibilities |V|² and closure phases T₃ϕ can be fitted with our model, there are series of |T₃| observations showing a clear trend with respect to the projected baseline B/λ, but with sudden increases of the amplitude. These are likely to be of an instrumental origin. Consequently, we decided to fit only the squared visibilities |V|² instead. This should not affect the fitting in a negative way, because |V|² observations should constrain the model anyway and the |T₃| data do not bring important constraints due to the short baselines of NPOI observations.

4.5. Simplifications reducing the computational time

The evaluation of χ² represented by Eqs. (4)–(6) for all available data turned out to be very demanding. The total time required for the computation of χ² exceeded two hours even for a moderate resolution 0.6 R_⊙ per pixel, and the grid size n_x × n_y × n_z = 251 × 251 × 126. Such a long computational time prevents an extensive search of the parametric space, so three approximations were introduced: (i) low resolution. All models were computed with the grid resolution 1 R_⊙ per pixel and grid size n_x × n_y × n_z = 161 × 161 × 81; (ii) several model parameters were fixed at values determined by previous investigators of β Lyr; and (iii) a “binning” of synthetic observable quantities was introduced. The last approximation differed for the magnitudes and interferometric observables.

Synthetic magnitudes were not computed for each observation time. Instead a synthetic light curve as a function of orbital phase was sampled equidistantly with one hundred points for each passband. Synthetic magnitudes for each observation time were then obtained by a cubic-spline interpolation of the synthetic light curve.

The binning in case of the interferometry limited the number of images I_v(x, y, t, λ_eff) computed to derive synthetic interferometric observables. The binning was introduced into the effective wavelength λ_eff and the orbital phase ϕ(t). The bin size was set to 100 nm for λ_eff and 0.001 for ϕ. The function I_v(x, y, t, λ_eff) was not sampled equidistantly as it was for the photometry. Instead the following simplification scheme was adopted:

1. List pairs (λ_eff, ϕ) for all observations.

2. Round the lists to the bin size (precision) given in the preceding paragraph.

3. Remove repeating items.

4. Compute I_v(x, y, t, λ_eff) only for the remaining items.

The image defined by the pair (λ_eff, ϕ) that is the nearest to the observation was chosen for the computation of interferometric observables. This means that error in the effective wavelength and orbital phase introduced by this approach cannot exceed half-width of their respective bin (Δλ_eff = 50nm, Δϕ = 0.0005 for this particular application). We note that for computations of the spatial frequency (B/λ) the exact value of effective wavelength for a given observation is used.

Using all these steps, the computational time was reduced by a factor of ten or more, and the evaluation takes about three minutes (on a 8-core processor). The parallelization was achieved in the Python wrapper by employing the standard multiprocessing module.

4.6. Modeling strategy

Advances achieved in the previous studies of β Lyr allowed us to significantly reduce the number of optimized parameters and focus more on the properties of accretion disk. The optimized parameters were: the inclination i, the longitude of ascending node Ω, the outer radius R_out of the disk, the semi-thickness H of the disk, the radial density ρ(R) and temperature T(R) profiles, and the distance d to the system. An attempt to include the semimajor axis a and the donor temperature T_eff,d among the optimized parameters has been done, but the former turned out to be completely correlated with the systemic distance, and the latter with the disk temperature. Therefore both were set to the values reported in earlier studies, even though they were likely correlated in these studies, too. Model parameters that were kept fixed during the optimization are listed in Table 7.

The initial set of parameters was based on models developed by Linnell (2000) and Ak et al. (2007). As a first step, photometry, visible interferometry, and infrared interferometry were each fitted independently of the remaining data. By this procedure, we established that all data are compliant with a similar model, although not exactly the same. Based on this information a generous range was selected for each optimized parameter and a search for global minimum of Eq. (4) was run over these ranges. Once the global fit converged we ran a local fit using the simplex algorithm to polish the result of the global method. This approach was repeated for each eligible combination of shape and radial temperature profiles.

Uncertainties of the optimal set of parameters were estimated from the convergence of the global fit. The resulting χ ² was scaled down to the ideal value, that is number of degrees of freedom. All solutions had their χ ² scaled down by the same factor as the optimal set of parameters. Solutions whose probability P_χ²(p) > 0.05, where P_χ² is the cumulative distribution function of the χ² probability density function, and p is the vector of optimized parameters, were accepted as possibly correct ones. The maximal differences between all accepted results and the optimal one were adopted as uncertainties.

4.7. Results

The following seven models of the accretion disk were optimized: (i) slab with power-law temperature and density, (ii) slab with steady-disk temperature and power-law density, (iii) wedge with power-law temperature and density, (iv) wedge with steady-disk temperature and power-law density, (v) lens with power-law temperature and density, (vi) envelope with homogeneous temperature and density profiles, (vii) nebula with power-law radial and exponential vertical density profiles.

Promising optimal models selected out of (i) to (vii) are shown in Fig. 5.

Fig. 5. Synthetic images of β Lyr A system for five different models: slab power-law, slab steady, wedge power-law, wedge steady, and nebula. The scale of grays corresponds to the monochromatic intensity Iv (in erg s−1 cm−2 sr−1 Hz−1). The wavelength is always λ = 545 mn and the orbital phase 0.25; α goes along east–west direction, and δ along north–south direction. It is worth noting that each model converged independently, but the outcomes are remarkably similar in terms of geometry.

Fig. 6. Synthetic images of β Lyr A for the best-fit nebula model, with total (not reduced) χ 2 = 103 233, shown for four different wavelengths: λ = 155 nm (FUV), 545 nm (V band), 1630 nm (H), and 4750 nm (M). The axes correspond to the right ascension α and declination δ (in mas), while the color scale corresponds to the monochromatic intensity Iv (in erg s−1 cm−2 sr−1 Hz−1). This is a small subset of all 2392 images (per one iteration) used to derive light curves, interferometric visibilities, closure phases, and triple product amplitudes.

Fig. 7. Observed and synthetic light curves of β Lyr A, shown for all 21 datasets (see their names in the right column). The light curves are phased according to 1. Vertical offsets are arbitrary. The best-fit model is again “nebula” with χ 2 = 103 233; the individual contribution arising from light curves comparison is ${\chi }_{LC}^2 = 6 918$. Synthetic data are denoted by yellow crosses, observed data by blue error bars, and residua by red lines (or circles). There are clear systematic differences especially for datasets iue.1250, iue.1365, oao2.1910, oao2.3320. At the same time, there are neighboring datasets matched relatively well. Sometimes, an intrinsic variability can be also seen (oao2.1430).

Fig. 8. Similar comparison as in Fig. 7, but for squared visibilities |V|2, with a contribution ${\chi }_{V^2}^2 = 54 137$. The |V|2 values are plotted against projected baseline B/λ (in cycles), and shifted vertically according to the dataset number. The are CHARA/MIRC data at the bottom, NPOI in the middle, and CHARA/VEGA at the top. Synthetic data are denoted by yellow crosses, observed data by blue error bars, and residua by red lines. A few outliers with large uncertainties, which do not contribute much to χ2 anyway, were purposely removed from the plot to prevent clutter. Even though there are some systematic differences for individual segments of data, overall trends seem to be correctly matched.

Fig. 9. Similar comparison as in Fig. 8, but for closure phases T3ϕ (left) and triple product amplitudes |T3| (right). Contributions to the total χ2 are ${\chi }_{CP}^2 = 29 153$, and ${\chi }_{T_3}^2 = 13 023$. As before, the values are plotted against projected baseline B/λ. T3ϕ measurements were available for NPOI (top half) and CHARA/MIRC (bottom half), while only |T3| from MIRC instrument were used.

The analysis led to the following findings:

• The fit of individual data subsets (photometry, visible and infrared interferometry) converged to similar, but not identical solutions. A comparison of these particular solutions against all available data have shown clear discrepancies for data that were not fitted. For example, a model fitting the infrared interferometry did not reproduce depths of minima of the visible light curve. A detailed uncertainty analysis was not carried out though, and the relative importance of these discrepancies was not quantitatively evaluated. We report this, because it might point to a possible discrepancy in our model and/or to a systematic effect that is still affecting our data and that was not suppressed with steps described in Sect. 4.4.

• The model with homogeneous envelope served primarily as a test whether the accretion disk around gainer could be homogeneous in temperature and density. This rather non-physical model tested our sensitivity to radial temperature and density profiles and it has shown that the sensitivity is indeed limited because the model with this profile shows only slightly higher total χ² than the other inhomogeneous models (by approximately 5%, which is statistically significant, though). The geometric size of the disk is similar to those obtained for the remaining models (see Table 9) – filling factor $f_{\text{f}} = {0.910}_{-0.073}^{+0.087}$, translates into outer (point) radius $R_{\text{out}} = {34.8}_{-2.8}^{+3.3}{R}_{\odot }$, semi-thickness $H = {9.2}_{-1.2}^{+1.1}{R}_{\odot }$, and temperature $T = 7 {230}_{-310}^{+480}$ K that is reached all over the surface, because no limb-darkening is present. It also leads to a higher orbital inclination $i = {95.4}_{-1.9}^{+0.5}$ deg. Nevertheless, the geometrical configuration and higher χ² make this model implausible.

• The model with a lens-shaped disk has been definitively ruled out. It was tested with both temperature and density radial profiles, but the χ² was significantly higher, by approximately 35%, than that obtained for the remaining shapes.

• Both slab and wedge disk shapes turned out to be plausible. Nevertheless, if we carefully compare the resulting χ² values summarized in Table 8, we can exclude the slab with a steady temperature profile, because its (not reduced) χ² = 114 981 is well above the 3 − σ level (that is 105 669) inferred from the best-fit model. The central region of the slab seems too cold, and especially the UV light curves do not match the observations (cf. its ${\chi }_{\text{LC}}^2$ contribution). In case of the wedge, the steady profile is compensated by the central star (gainer), which is partially visible in the opening. The resulting χ² range from 103 644 to 105 202, that is still within the 3 − σ level, or in terms of the reduced ${\chi }_{\text{R}}^2 = 3.80-3.86$.

• The nebula disk model provides the best-fit, with χ² = 103 233, which is equivalent to ${\chi }_{\text{R}}^2 = 3.78$. Even though it is not significantly better on its own, if we focus on a subset of observational data, namely the light curves (${\chi }_{\text{LC}}^2$ in Table 8), this fit is indeed significantly better than the others and it is thus our preferred model.

The optimal sets of parameters for the plausible models are listed in Table 9. The intervals in which the optimal solution was searched for with the global-minimization algorithm are given for each optimized parameter. Although the solutions were mostly equal from point of the total χ², only the nebula model is plotted against the data in Figs. 6–9.

5. Discussion

This section compares our results from Sect. 4 with theoretical predictions from earlier studies of β Lyr.

5.1. Distance and inclination of β Lyr A

The distance of β Lyr A inferred from our model is d = 319.7 ± 2.7 pc, which is a bit larger distance than our preliminary expectation based on β Lyr B. The final result is essentially based on a sin i = 58.19 R_⊙ value (see Table 7), which is ≲1% larger than the value accepted by Zhao et al. (2008) and is well within their uncertainty interval. Hence the choice of a slightly different s sin i cannot cause the difference between our distance and that of Zhao et al. (2008). Nevertheless, our estimate is based on a more correct physical model than that by Zhao et al. (2008), whose estimate closest to ours relied on a model with two uniform ellipses.

Our analysis reinforced the hypothesis that β Lyr A and β Lyr B share a common origin (see Sect. 3.1.3). In particular, Gaia distance d _B = (333 ± 6) pc of β Lyr B is in a fair agreement with our model distance estimate of β Lyr A. The spectroscopic distance of β Lyr B derived in this study is significantly lower (259 pc), but that may result from calibration uncertainties (see the discussion at the end of Sect. 3.1.3) or from possible, still unrecognized duplicity of β Lyr B.

Concerning the inclination of β Lyr, our estimates obtained for the wedge- and slab-shaped disks do not agree with each other. This is not very surprising, because the wedge shape tends to attenuate the radiation from the hot central parts of the disk (and the gainer). Because this radiation is in fact observed, the fit converges to a lower inclination to expose central parts, and to compensate for the attenuation. The inclination of the nebula model is just in between.

5.2. Properties of the accretion disk

A critical discussion of disk parameters is presented here:

• Radius of the accretion disk. Dense parts of the accretion disk actually fill the corresponding Roche lobe. The front, back, and side radii are R_limit,g = 37.4 R_⊙, 31.7 R_⊙, and 30.3 R_⊙ respectively. The disk that was obtained reaches up to the Roche limit, although the “hard” upper limit of our optimisation procedure was as high as 35 R_⊙; the value is thus constrained by our observations. The tidal cutoff radius is as low as 26.3 R_⊙ for the given mass ratio q, but this is not necessarily the edge of a viscous disk (Papaloizou & Pringle 1977). In earlier studies (e.g., Linnell 2000) the disk was modeled by a solid body with prescribed radiative properties. Outer radii of our disk models R_out listed in Table 9 cannot be thus directly compared to radii obtained in earlier studies, because our disks are not optically thick starting from their rim. To obtain a comparable radius, a pseudo-photosphere approximated by the optical depth τ = 2/3 has to be found. It was searched along lines of sight perpendicular to the disk rim (x, y = 0, z = 0); x ∈ [0, R _out] R_⊙. It was realized that the photosphere forms almost up to the geometric R_out for the slab-shaped disk, $R_{\text{out},\text{ph}.}^{\text{slab}}\doteq 30R_{\odot }$, and slightly less for the wedge-shaped disk $R_{\text{out},\text{ph}.}^{\text{wedge}}\doteq 29R_{\odot }$, where index “ph” stands for the photosphere. Hence, our photospheric disk radius is in excellent agreement with that obtained by Linnell (2000), R = 30 R_⊙, and also with Mennickent & Djuraševć (2013), R = (2.83 ± 0.3) R_⊙ (if uncertainty ≈1.5R_⊙ is taken into account). This is demonstrated in Fig. 10, where the physical position of the photosphere within the accretion disk is shown for the wavelength 1630 nm.

• Semi-thickness of the disk. For the slab-shaped model, the value of H is the same as semi-thickness of the disk photosphere. For the wedge-shaped model, we searched for the physical position of photosphere along the following lines of sight perpendicular to the disk rim: (x = 0, y = 0, z); z ∈ [0, H]. This is qualitatively demonstrated in Fig. 11. Up to a certain z the wedge is opaque; for intermediate z the ray pierces through the first lobe and the optical depth τ = 2/3 is reached in the second lobe. For high z the disk is optically thin. Hence the semi-thickness of the “opaque” wedge is about 30% smaller than the geometric semi-thickness H = R_out sin υ ≐ 7 R_⊙ derived from Table 9. Our semi-thickness of the opaque disk ≃5 R _⊙ is in agreement with the result obtained by Mennickent & Djuraševć (2013), H = (5.50 ± 0.10) R_⊙, and substantially lower than the value by Linnell (2000), H = 8.0 R_⊙.

• Shape of the disk. Unfortunately, it seems almost impossible to distinguish between slab, wedge and nebula shapes. Nevertheless, it is clear that their semi-thickness is so large that the disk cannot be in a vertical hydrostatic equilibrium. A hydrostatic disk with a constant vertical temperature profile, T(z) = const., would have an exponential density profile (e.g., Pringle 1981):$$\rho (r,z) = {\rho }_0(r)\text{exp}\left(-\frac{z^2}{2H_{\text{e}q}^2}\right),$$(13)with the characteristic scale H_eq given by the temperature profile T(r):$$H_{\text{e}q}(r) = \sqrt{\frac{ℛT}{\mu }}\frac{1}{{\text{Ω}}_{\text{k}}},$$(14) where ℛ denotes the ideal-gas constant, µ the mean molecular weight, and Ω_k the Keplerian angular velocity. As we also verified with SHELLSPEC (using the nebula model with h _mul = 1), the resulting H_eq for our range of temperatures (30 000–7000 K) is always low, 0.2–1.2 R_⊙, and the gainer would be always visible. As a consequence, H ≫ H_eq is a proof that the disk is non-equilibrium, and the flow starting from the donor must have a non-negligible vertical velocities within the accretion disk. It is a matter of dynamical models (or spectro-interferometry in individual lines) to constrain the velocity field.

• Radial density profile and disk mass. The continuum data do not allow us to see below the pseudo-photosphere. This is evident from the high correlation between ρ₀ and α_D resulting from individual models, corr(ρ₀, α_D) ≃ 05–0.8. Hence our model does not provide the true density radial profile, but rather the minimal profile required to “produce continuum” at the correct radius and height. The disk mass given by our radial profiles is then ≃10^–4–10^–3 M_⊙. This disk mass estimate is essentially the same as that obtained by Hubený & Plavec (1991), Hubený et al. (1994). If it was real, it would suggest a very efficient accretion, because the accretion time scale τ_acc = M / Ṁ ≃ 5–50 yr is much shorter than the expected duration of β Lyr A mass-transfer phase. For a steady disk (Pringle 1981), it is simply assumed the viscosity always adapts to the constant accretion rate and the viscous time scale is then equal to the accretion one, τ_visc = τ_acc.

• Disk-rim temperature. The rim temperature given by Eqs. (8) and (9) is not directly comparable to those obtained with solid surface models, because the photosphere does not coincide with the geometric rim. To overcome this contradiction, one can adopt the photospheric temperature, that is shown in Figs. 10 and 11. Another approach is to take the mean value of the intensity distribution over the photospheric surface 〈I_disk(x, y, λ)〉 for a given wavelength and find a temperature corresponding to the Planck law B_v(T, λ) for the mean intensity. Using this approach, we determined the following rim temperatures at three wavelengths (given by the upper index in nm; Table 10):

  The rim temperature of slab-shaped models is comparable to that obtained by Mennickent & Djurašević (2013), T_rim = (8200 ± 400) K. It is difficult to tell which radial temperature profile is correct. The power-law seems to give a temperature slightly above, and steady-disk gives a temperature that is slightly below the estimate by Mennickent & Djurašević (2013). Linnell (2000) generally gives larger rim temperatures. Most of the rim of his accretion disk has T_rim = 9000 K, with two strips having twice higher temperature at the top and bottom of the disk. Spectroscopic studies (e.g., Harmanec & Scholz 1993; Ak et al. 2007) do not provide an accurate estimate, because spectral types from an early F-type to late A-type were attributed to rim of the accretion disk.

• Radial temperature profile. It is almost impossible to determine the whole temperature profile, because the disk is opaque (the photosphere forms at most a few solar radii below the disk rim), and the orbital inclination is very close to 90°, which prevents us from seeing the disk face-on. The first issue manifests in itself for the power-law radial temperature profile by the extreme correlation between the inner rim temperature T₀ and the exponent of the power-law ρ(T₀, α_T) = –0.92 for both models with slab and wedge. A comparison of models with slab and different radial temperature profiles (given by Eqs. (8) and (9)), and the temperature profile obtained by Mennickent & Djurašević (2013) is shown in Fig. 12. The radial profiles agree with each other up to ≃ 10 R_⊙ below the disk rim, that is below the pseudo-photosphere. The radial profile given by Eq. (9) derived by Shakura & Sunyaev (1973) is heated only by the viscous dissipation, but the disk in β Lyr enshrouds a B0.5 V star that must considerably heat the disk, too. A problem is that asymptotically a passive irradiated disk has the same radial temperature dependence T ∼ R^−3/4 as a steady disk (Friedjung 1985; Hubený 1990). Hence it is difficult to discern the two radial temperature models from their behavior in outer part of the disk. Calvet et al. (1991) modeled proto-stellar disks and found that the central regions of their disks are significantly heated by the embedded star. It is interesting to note that the temperature T₀ of the power law behavior is quite close with the temperature of the central star. This is a good indication for the power law temperature behavior in the disk. From afar the proto-stellar accretion disks are very similar to that surrounding β Lyr – the accretion rate is high, and there is a star in its center. Hence there are reasons to believe that the power-law provides a better description of the radial temperature profile. Finally the best-fitting power-law α_T = −0.95 is steeper than the canonical value (−3/4). This may suggest a presence of a transition layer, where the temperature falls more steeply at the outer disk rim.

Fig. 10. Top panel: optical depth τ (computed for the wavelength λ = 1630 nm) of the β Lyr A disk, which is observed approximately edge-on (that is from the top). The coordinate z thus corresponds to the line of sight, while x (and y) to the sky plane. The center is empty because the gainer is a non-transparent object. In this case, the inner-rim density of the wedge power-law model is set low, ρ0 = 10−9 g cm−3, to demonstrate a separation from the outer rim. The dashed line shows the physical position of the photosphere (τ = 2/3). The lines-of-sight grid exhibits “steps” due to the limited spatial resolution of our model (1 R⊙), even though the integration of the radiation transfer (the contribution function) is internally performed on a much finer grid. Bottom panel: corresponding photospheric temperature Tphot of the disk, which varies along with the radial temperature profile T(R).

Fig. 11. Optical depth τ in the vertical (perpendicular) cross-section of the disk; other parameters are the same as in Fig. 10. The lines of sight are seemingly different from the wedge shape, but this only because they start either in the vacuum, or at the non-transparent object.

Fig. 12. Comparison of radial temperature profiles of the accretion disk surrounding the gainer. The dotted line is the profile obtained by Mennickent & Djurašević (2013), and the dashed line is mean temperature of the two-temperature model developed by Linnell (2000). The red line represents the best-fitting power-law radial temperature profile (solution slab/pl in Table 9), and the red belt all plausible solutions given by the uncertainty of the inner rim temperature and the exponent of the power-law, and the red point position of photosphere (τ = 2/3), where temperature Tph. = 7 192 K is reached. The blue line represents the best-fitting steady-disk radial temperature profile (solution slab/sd in Table 9), the blue belt all plausible solutions given by uncertainty in the inner rim temperature, and the blue point position of photosphere (τ = 2/3), where temperature Tph: = 7 192 K K is reached. The photospheric temperatures were computed for a line of sight in the plane z = 0 and piercing through center of the accretion disk

5.3. Presence of “a hot spot”

The existence of a region heated by an interaction of the incoming flow and the accretion disk is usually required by theoretical models (e.g., Lubow & Shu 1975). For β Lyr Lomax et al. (2012) and Mennickent & Djurašević (2013) used the hot spot to explain a presence of bumps (or irregularities) in the light curve and polarized flux. We carried out an attempt to confirm their findings in continuum. A spot represented by a homogeneous sphere has been added to the slab power-law model (see Fig. 13). Parameters defining its radial position r_s within the accretion disk, the position angle θ _s with respect to line joining centers of both binary components, radius R_s, density ρ _s, and temperature T_s were optimized using the differential evolution and simplex algorithms; the optimal values are listed in Table 11.

Fig. 13. Synthetic images of slab power-law model which was further improved by a “spot” (that is a spherical structure between the primary and secondary, partly hidden inside the slab). The total (not reduced) χ2 value was decreased from 105 202 (without the spot) down to 102 005, which is a significant improvement. Four wavelengths are shown: λ = 155 nm (FUV), 545 nm (V band), 1630 nm (H), and 4750 nm (M).

First, only spot parameters were converged. As the spot is a substantial non-axisymmetric feature, the fit converged quickly to a location just between the primary and the disk. The spot temperature T_s about 10 000 K is logically between those of the primary and the outer rim of the disk. Second, all parameters were set free and converged again, because the original parameters might have been affected by a systematic error of the model (namely the missing spot). The procedure helped to decrease the original (not reduced) χ² = 105 202 down to102 005, which is a statistically significant improvement and we thus may confirm the existence of the spot on the basis of continuum observations. We verified that adding a spot to other models leads to an improvement of the same order.

As illustrated by Fig. 13, the “spot” detected by our modeling is not a tiny structure corresponding to the area of interaction between the gas stream from the donor and the disk. It likely represents an illuminated part of the disk rim, where the reflection of donor light, and irradiation heating occur. We note that such lateral temperature gradient has been observationally proven for another object with an optically thick disk, ε Aur (cf., e.g., Hoard et al. 2012).

In principle, it should be possible to add a second spot to the model, and so on, and expect further improvements of the χ², but we prefer to keep a simple model as long as possible. Otherwise, systematic uncertainties among different types of observations (light curves, squared visibilities, closure phases, etc.) might be hidden by a complex model. Moreover, there are techniques (like spectro-interferometry in lines, or Doppler tomography) better suited to pinpoint the orbital position of such non-axisymmetric features.

5.4. Comparison of SED from models and observations

Although we did not attempt to fit the spectral energy distribution (SED) of the system, it would be an important check of the model. That is why we compare synthetic SEDs to the observations of Burnashev & Skulskii (1978; see Fig. 14). It seems inevitable that SEDs exhibit some systematic offsets. Also the resolution of synthetic spectra does not match the observations. Moreover, we cannot expect that emission lines will be computed correctly, because the model is mostly focused on opaque medium and continuum flux. However, the absolute fluxes and their spread among primary minimum, secondary minimum and out of the eclipses seem to at least roughly correct. The expected systematic uncertainties of the observations (5%) might be almost of the same order, especially in UV, where the extinction is strong and variable. We thus believe it should be possible to match the observed SEDs with a future version of our model.

Fig. 14. Comparison of synthetic spectral-energy distributions (thick lines) of β Lyr A and observed SED (thin lines) according to Burnashev & Skulskii (1978). The monochromatic flux Fλ (in erg s−1 cm−2 cm−1) was measured in the range of λ = 330–740 nm. The color scale correspond to the orbital phase. Observation uncertainties are on the order of 5%. The model (corresponding to a “nebula” with tinv = 4.53) was not converged with respect to these observations and the SEDs thus inevitably exhibit some systematic offsets.

6. Conclusion and outlook

The properties of opaque bodies within β Lyr A system were studied. Our analysis was primarily targeted on the properties of the accretion disk surrounding the mass-gaining component of this close interacting binary.

For the description of interacting binary systems, we created a tool based on the SHELLSPEC code. It permitted us to significantly improve the modeling of stellar systems and the computations of radiation transfer in the co-moving circumstellar medium by Budaj & Richards (2004). Apart from improvements suitable for an optically thick medium, we can also compute interferometric observables and proceed with both local- and global-optimisation methods.

We then constructed several disk models that differed in shape, density and temperature profiles. These models were fitted to series of spectro-interferometric and photometric observations, both sampling the whole orbit. We also compared our results to those obtained by earlier investigators of the system (especially Linnell 2000; Mennickent & Djuraševć 2013) and to theoretical models of accretion disks.

Our results indicate that the opaque parts of the accretion disk have the outer radius R_out = (30.0 ± 1.0) R_⊙, the semi-thickness H = (6.5 ± 1.0) R_⊙ (for slab and wedge shapes), or equivalently the scale-height multiplication factor h_mul = 4.3 ± 0.3 (for nebula model; see the overview in Fig. 4). But the true location of the disk pseudo-photosphere slightly depends on the wavelength. The minimum mass should be 10^–4–10^–3 M_⊙. Given the thickness, the disk clearly cannot be in a vertical hydrostatic equilibrium. We have also determined the orbital inclination i = (9.5 ± 1.0) deg (as an average and range for admissible models), the longitude of ascending node Ω = (253.7 ± 1.0) deg, and the probable distance to the β Lyr A system d = (319.7 ± 2.7) pc.

The power-law temperature profiles (and also the steady-disk for wedge) seem compatible with the observations, but the central values remain very unconstrained, because the disk continuum is formed only a few solar radii below the disk rim.

An addition of a hot spot to our model improved the χ², so that we can consider the existence of the spot to be confirmed in the continuum radiation, although it may be actually a compensation of missing reflection from the disk, or heating of the disk by the companion. Its position may also correspond to a flow of material from the primary (donor).

The radiative and kinematic properties of neighboring β Lyr B have been determined too. Even though we were unable to prove β Lyr A and B orbit each other, they both likely originate from the same association.

This study presents a springboard to forthcoming analyses of the optically thin circumstellar medium in β Lyr A – it is crucial to know the properties of the opaque material too. Using a series of spectroscopic and spectro-interferometric observations of strong emission lines we intend to resolve and describe the structure and kinematics of the optically thin medium within this remarkable system. Consequently, it should be possible to better determine the radial profiles of the disk atmosphere. At the same time, the mass of jets would provide an accurate estimate of the mass and angular momentum loss from the system, which would offer an invaluable test for models of mass transfer in binary systems.

Acknowledgments

The constructive criticism of an earlier version of this manuscript by the anonymous referee is appreciated. This research was supported by the grants P209/10/0715, GA15-2112S, and GA17-00871S of the Czech Science Foundation, by the grant no. 250015 of the Grant Agency of the Charles University in Prague. This work is based upon observations obtained with the Georgia State University Center for High Angular Resolution Astronomy Array at Mount Wilson Observatory. The CHARA Array is supported by the National Science Foundation under Grant No. AST-1211929, AST-1411654, AST-1636624, and AST-1715788. Institutional support has been provided from the GSU College of Arts and Sciences and the GSU Office of the Vice President for Research and Economic Development. We thank M. Zhao for participating in the MIRC observations, and B. Kloppenborg for his support in the earlier attempts to use the SIMTOI modeling tool. Two Reticon spectra of β Lyr B were obtained by Dr. P. Hadrava. HB acknowledges financial support from the Croatian Science Foundation under the project 6212 “Solar and Stellar Variability”. The work of JB was supported by the VEGA 2/0031/18 and APVV-15-0458 grants. We acknowledge the use of the electronic database from the CDS, Strasbourg, and the electronic bibliography maintained by the NASA/ADS system. This work has made use of data from the European Space Agency (ESA) mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, http://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular those participating in the Gaia Multilateral Agreement. In April 2012, our late friend and colleague Olivier Chesneau pointed out to some of us the publication by Lomax et al. (2012) and expressed these inspirational words: “… the hot spot… is quite extended and should show up in the interferometric observables”. This was the starting point of this study and our results are dedicated to him.

References

Appendix A:

Details on the interferometric data reductions

Additional details on the reduction process of interferometric observations are presented in this section. Detailed characteristics of all observations are presented in Table A.1. Individual interferometric observations are available in the OIFITS format through CDS. These files contain only reduced data, that is calibrated squared visibilities and closure phases. Raw data are made available upon request.

A.1. Details on the reduction of CHARA/VEGA observations

The observations were carried out during a 11 days long campaign in 2013. Both MIRC and VEGA were observing at the same time, sharing the visible and infrared photons. The majority of observations was taken in the three-telescope (3T) mode. The only two exceptions were the first (22nd Jun, 2013) and the last (1st Jul, 2013) nights. Observations from those two nights were taken in two-telescope mode. Only four CHARA telescopes (denoted E1, E2, W1, and W2) were used. The long baselines were in general east–west oriented and the short baselines north–south oriented, because the projection of β Lyr orbit on the sky is roughly east–west oriented (see Fig. A.1). The length of projected baselines ranged from ≃65 m to ≃250 m.

Fig. A.1. (u, v) coverage of all interferometric observations. Colors correspond to three different instruments: NPOI (blue), CHARA/MIRC (green), CHARA/VEGA (magenta).

Fig. C.1. Comparison of an observed spectrum of β Lyr B obtained on RJD = 49866.4348 and disentangled spectra of β Lyr B with the best fitting spectra corresponding to their respective solutions listed in Table 3, that is synthetic spectra compared to observed spectrum correspond to slightly different parameters than those compared to disentangled spectra. The blue line denotes disentangled spectrum, the black line the observed spectrum, the red line best-fitting synthetic spectrum, and the yellow line fit residuals. Only principal spectral lines (Si II 6347 Å, Si II 6371 Å, Hα, He I 6678 Å) that are present within the studied spectral region Δλ ≃ 6200–6800 Å and their surroundings are plotted.

The observations were carried out with two cameras in four passbands that were centered at following wavelengths λ _C ∈ {535, 656, 706.5, 815} nm. Resolution of the recorded spectra was R = 5000. Individual frames were recorded with frequency 100 Hz, and were grouped into blocks containing 2500 frames. An observation typically contained ≃20 blocks. Within these blocks the frames were coherently summed and raw squared visibility was determined for each block. The whole passband was not used, but two narrow channels were selected in each passband. The following channels were selected: Δ_CH ∈ {525–535, 535–545, 643–653, 658–667, 687–702, 709–724, 810–825, 825–840} nm. The channels avoid major spectral lines in these regions. The only exceptions are regions 643–653 nm, and 658–667 nm, which are partially affected by wings of Hα line, and 825–840 nm, which is affected by water-vapor lines. Narrow channels were chosen because of the limited coherence of the waves due to the atmospheric turbulence.

Four calibrators were observed in order to calibrate the instrumental visibilities. Due to variability of atmospheric conditions during night, a calibrator was observed before and after each observation of β Lyr. A list of calibrators and their properties are listed in Table A.2. The calibrators were chosen with help of tool SearchCal³ developed by Bonneau et al. (2006). The uniform-disk diameters of calibrators were taken from the JMMC catalog of stellar diameters by Lafrasse et al. (2010).

In order to avoid modeling of highly inaccurate observations, all blocks that have instrumental visibility with S/N < 2 were removed. Also we performed simple filtering based on residual optical path delay (OPD). Correct blocks of visibility measurements have very similar OPD. Hence if OPD of one (or more) blocks deviate significantly from a mean OPD based on all blocks acquired within one measurement, this block is very likely wrong. Therefore if a block OPD differed from the mean by more than two standard deviations, it was removed.

A.2. Details on the reduction of CHARA/MIRC observations

The instrument CHARA/MIRC was used to measure squared visibilities and closure phases. MIRC performs real-time group delay tracking in the H-band. The observations come from two observational runs: (i) 2006–2007 campaign that was already analyzed by Zhao et al. (2008) and whose description can be found there, and (ii) 2013 campaign, together with VEGA, whose description follows. All observations used here were taken across the near-IR H-band.

Since 2011, MIRC can combine light from six telescopes, so it is able to record 15 squared visibilities and 20 closure phase observations. The H-band is split into eight channels with absolute wavelength accuracy ±0.25%. A more thorough description of the instrument is in studies by Monnier et al. (2010), and by Che et al. (2010), Che et al. (2012). Using Fourier transform techniques, the visibilities are measured, averaged, and corrected for biases. The bispectrum is formed using the phases and amplitudes of three baselines that form a closed triangle (Monnier et al. 2007). Amplitude calibration was performed using real-time flux estimates derived through use of a beam splitter following spatial filtering for improved performance. Lastly, observations of reference calibrators (see Table A.2) throughout the night allowed for correction of time-variable factors such as atmospheric coherence time, vibrations, differential dispersion, and birefringence in the beam train. The uncertainties of closure phase measurements were adjusted (see Sect. 4.4). MIRC and VEGA shared the same calibrators for this observing campaign.

A.3. Details on the reduction of NPOI observations

All NPOI observations were taken with the six-beam combiner. Visibilities, complex triple amplitudes, and closure phases were recorded in 16 narrow-band channels between 5500 Å and 8500 Å. The calibrators are taken from a list of single stars maintained at NPOI with diameters estimated from V and (V–K) using the surface brightness relation by Mozurkewich et al. (2003) and van Belle et al. (2009) (see Table A.2). Values of E(B–V) were derived from comparison of the observed and theoretical colors as a function of spectral type by Schmidt-Kaler in Aller et al. (1982). Values for the extinction derived from E(B–V) were compared to estimates based on the maps by Drimmel et al. (2003), and used to correct V if they agreed within 0.5 mag. Even though the surface brightness relation based on (V–K) colors is to first order independent of the reddening, we included this small correction. For these observations, only one calibrator was used: HD 176437.

NPOI data and their reductions followed the procedure described by Armstrong et al. (1998) and Hummel et al. (2003). A pipeline written in GDL⁴ was used for the OYSTER⁵ NPOI data reduction package. The pipeline automatically edits the one-second averages produced by another pipeline directly from the raw frames, based on expected performance such as the variance of fringe tracker delay, photon count rates, and narrow-angle tracker offsets. Visibility bias corrections are derived as usual from the data recorded away from the stellar fringe packet. After averaging the data over the full length of an observation, the closure phases of the calibrators were automatically unwrapped so that their variation with time, as well as that of the visibility amplitude, could be interpolated for the observations of β Lyr. For the calibration of the visibilities, the pipeline used all calibrator stars observed during a night to obtain smooth averages of the amplitude and phase-transfer functions using a Gaussian kernel of 80 min in length. The residual scatter of the calibrator visibilities and phases around the average set the level of the calibration uncertainty and was added in quadrature to the intrinsic data errors. The amplitude calibration error of typically a few percent in the red channels up to 15% in the blue channels was added in quadrature to the intrinsic error of the visibilities. The phase calibration was good to about a couple of degrees.

Appendix B:

Details on the photometric data reductions

Hvar UBVR observations are differential observations, relative to χ Lyr (HD 176437), for which the following mean Hvar all-sky values from excellent nights were adopted:$$\begin{array}{c}V=3.253 \text{mag},\\ B-V=-0.0644 \text{mag},\\ U-B=-0.0308 \text{mag},\\ V-R=-0.0024 \text{mag}.\end{array}$$

The check star v² Lyr (HD 174602)$$\begin{array}{c}V=5.243 \text{mag},\\ B-V=0.0980 \text{mag},\\ U-B=0.1038 \text{mag},\\ V-R=0.1000 \text{mag},\end{array}$$was observed as frequently as the variable and β Lyr B was also observed as another check star on a number of nights. All Hvar UBVR observations were transformed to standard system through non-linear transformation formulæ using the HEC22 reduction program (see Harmanec et al. 1994; Harmanec & Horn 1998, for the observational strategy and data reduction)⁶. All observations were reduced with the latest HEC22 rel.18.2 version of the program, which allows the time variation of linear extinction coefficients to be modeled in the course of observing nights. The UBVR observations were reduced to Johnson bright standards, for which we derived robust mean UBVR values from individual observations published by Johnson (1966). The uncertainties of these observations estimated from measurements of the check star are σ_U,1 = 0.010 mag, σ_B,1 = 0.014 mag, σ_V,1 = 0.007 mag, and σ_R,1 = 0.013 mag. We underline that our R magnitudes were reduced to Johnson, not Cousins R values.

Differential Johnson–Cousins BV(R)_c observations were acquired at private observatory of Mr. Svoboda in Brno, Czech Republic using SBIG ST-7XME CCD camera mounted at 34 mm refractor. β Lyr and comparison star γ Lyr were observed simultaneously. The atmospheric extinction was assumed constant over the image. Seasonal transformation of the instrumental magnitudes into standard system was carried out using linear formulæ similar to Eq. (4) in Harmanec et al. (1977). Upper limit on the uncertainties of these observations are σ_B,2 = 0.013 mag, σ_V,2 = 0.013 mag, and σ_R,2 = 0.010 mag.

All calibrated UBVR observations acquired at Hvar observatory, and differential BV(R)_c observations acquired by PS are listed in Table B.1. This table is available at CDS.

Appendix C:

Supplementary material to analysis of β Lyr B

Supplementary material to analysis of β Lyr B (see Sect. 3.1.3) is presented here.

C.1. Details on the spectroscopic observations

All 13 electronic spectrograms were obtained in the coudé focus of the 2 m reflector and have linear dispersion of 17.2 Å mm^–1 and two-pixel resolution 12 600 (11–12 km s^–1 per pixel). The first 7 spectra (until RJD ≃ 50 235) were taken with a Reticon 1872RF linear detector and cover a spectral region from 6300 to 6730 Å. Complete reductions (bias subtraction, flat-fielding, extraction of 1 d spectrum, wavelength calibration, normalization) of these spectrograms were carried out by PH with the program SPEFO. The remaining spectra were secured with a SITe-5800 × 2000 CCD detector and cover wavelength interval from 6260 to 6760 Å. Their initial reductions (bias subtraction, flat-fielding, extraction of 1-d spectrum, and wavelength calibration) were carried out by MŠ in IRAF⁷ and their normalization by PH in SPEFO. In both cases the stellar continuum was approximated by Hermite polynomials that were fitted through several (suitably chosen) continuum points.

Photometric observations of β Lyr B were obtained at Hvar observatory and their reduction procedure is described in Appendix B.

C.2. Kinematic and radiative properties of β Lyr B

Additional details on the measuring of RVs and the modeling of observed spectra with synthetic ones follow:

• The RV measurements of β Lyr B obtained manually with SPEFO, and through comparison with synthetic spectra using PYTERPOL are listed in Table C.1. Four spectral lines were measured with the manual method on each spectrum. Instead of individual measurements, their mean and corresponding standard deviation are listed, because measurements on each spectral line did not differ systematically from RVs measured on the remaining spectral lines. The manually measured RVs give an impression that they slowly vary, but similar trend is not present in automatic measurements or RVs measured on photographic plates.

• Uncertainties of kinematic and radiative properties of β Lyr B obtained through modeling of its observed (or disentangled) spectra with synthetic spectra (solutions 1 and 2 in Table 3) were obtained through Markov chain Monte Carlo simulation implemented within emcee⁸ Python library by Foreman-Mackey et al. (2013). The posterior probability distribution of for each individual optimized was fitted with a Gaussian function. Standard deviation of the function was taken for uncertainty of the optimized parameter. Only statistical part of the total uncertainty was pinpointed by this approach. The wavelength ranges fitted for observed spectra were Δλ = {6337–6410, 6530–6600, 6660–6690} Å, and for disentangled spectra Δλ = {6338–6605, 6673–6724} Å. The region Δλ = 6337–6410 Å was not modeled for observed spectra, because it contains only few weak stellar lines and is densely polluted by telluric lines.

• Agreement between the observed (or disentangled) spectra and their best-fitting synthetic spectra (given by solutions 1 and 2 in Table 3) is demonstrated by Fig. C.1.

Only the major spectral lines and their vicinity are plotted. Only fit of one observed spectrum with high S/N is shown. Jags in the observed spectrum are remnants of telluric lines. Also we note that PYTERPOL does not require that the modeled spectra are equidistant. Hence it was not necessary to fill the gaps in spectra that emerged after the removal of telluric lines.

C.3. Proper motion of β Lyr A and B

Proper motions of β Lyr A and B were downloaded from the Vizier portal. Each 2-d vector had a component along the declination µ_δ and in the perpendicular direction along the right ascension µ_α. The latter coordinate was corrected for declination of both systems. The following coordinates were used:$$X_{\mu } = {\mu }_{\delta },$$(C.1) $$Y_{\mu } = {\mu }_{\alpha }\hspace{0.17em}\text{cos}\hspace{0.17em}\delta ,$$(C.2)where δ denotes declination. For β Lyr A δ _A = 33.362667 deg, and for β Lyr B δ _B = 33.351856 deg were adopted. All studied records in coordinates given by Eqs. (C.1), and (C.2) are listed in Table C.2.

Appendix D:

A note on disk model with a vertical temperature jump

The previous version of SHELLSPEC contained a “nebula” disk model with an exponential vertical density profile ρ(z), and a possible jump in the corresponding temperature profile T(z) to mimic a hotter disk atmosphere irradiated by the star. For this model, we realized there are bright spots in FUV close to the outer rim (see Fig. D.1). They actually helped to decrease the ${\chi }_{\text{LC}}^2$ contribution.

Fig. D.1. Synthetic image of β Lyr A for λ = 155 nm (OAO2 band). This model is based on a “nebula” with a vertical temperature jump, T(z) profile is determined by its scale height H and multiplication factor h mul = 5.27, together with temperature inversion scale h inv = 5.19 and factor t inv = 4.87. The resulting total χ 2 = 99430, with individual contributions as low as ${\chi }_{\text{LC}}^2 = 4 932$, ${\chi }_{V^2}^2 = 53 534$, ${\chi }_{\text{CP}}^2 = 28 282$, ${\chi }_{{\text{T}}_3}^2 = 12 684$, which is significantly better than the nominal nebula model presented in Table 9. However, the artifacts (bright spots close to the outer rim; black in this color scale) are clearly visible.

This is a well-known problem caused by a simple linear interpolation of both ρ and T quantities within the critical step of the optical depth, where the medium changes from thick to thin, and there is a large source S_v (and contribution) function in the middle. When the resolution of the model is increased twice, or four times, these spots subsequently disappear and ${\chi }_{\text{LC}}^2$ increases again. Users should be aware of these discretization artifacts, because they sometimes appear in the course of convergence (e.g., when the orbital inclination changes). The same would be true for models with overlapping optically thick and optically thin objects, like slab + slab, or wedge + flow. Nevertheless, this problem is an indication for us that FUV radiation and corresponding light curves should be described by a more complete model of the disk atmosphere.

All Tables

All Figures

Fig. 1. Phase coverage of spectro-interferometric observations of β Lyr acquired by different instruments. δ denotes the relative declination (positive toward the north), and α the relative right ascension (positive toward the east). The black line shows the size and orientation of the β Lyr orbit in the sky, the blue dots show orbital phases corresponding to NPOI observations, the magenta dots to CHARA/VEGA observations, the green dots to CHARA/MIRC observations acquired in 2013, and the red dots to CHARA/MIRC observations acquired in 2006/2007. An arbitrary vertical shift of 0.2 mas is added to separate the various orbits.

In the text

Fig. 2. Available RVs of β Lyr B. The first three datasets are RVs measured on photographic spectra of moderate dispersion: green open circles denote the RVs measured by Petrie & Pearce (1961), the red open triangles come from Abt et al. (1962), and the black boxes are RVs from Abt & Levy (1976). The filled blue diamonds are the RVs measured on Ondřejov electronic spectra with PYTERPOL. The rms errors, available only for Abt & Levy (1976) and new RVs are also shown. We note that the most deviant RV of −45 km s−1 is based on a single line and comes probably from an underexposed spectrum.

In the text

Fig. 3. Proper motions of β Lyr A and β Lyr B. The black points denote proper motion measurements of β Lyr A, the black triangles proper motions of β Lyr B, the red point the mean weighted proper motion of β Lyr A, and the red triangle the mean weighted proper motion of β Lyr B. Xµ, Yµ are Cartesian coordinates of the proper motion vector given by Eqs. (C.1) and (C.2). The former is north–south oriented, and the latter east–west oriented.

In the text

Fig. 4. Geometric shapes of models of accretion disks. (I) Denotes the slab, (II) the wedge, (III) the lens, and (IV) the envelope. The latter is shown as if it was filling the Roche limit, that is Rout = RL1, where Rout = RL1 denotes point radius of the Roche lobe. Here it is shown this way to emphasize the geometry of the disk, but its shape may be more symmetric depending on its filling factor. See Sect. 4.3 for the description of individual disk shapes. The Nebula is not presented here, since no simple geometrical sketch can be given. It is illustrated, however, in Fig. 5.

In the text

Fig. 5. Synthetic images of β Lyr A system for five different models: slab power-law, slab steady, wedge power-law, wedge steady, and nebula. The scale of grays corresponds to the monochromatic intensity Iv (in erg s−1 cm−2 sr−1 Hz−1). The wavelength is always λ = 545 mn and the orbital phase 0.25; α goes along east–west direction, and δ along north–south direction. It is worth noting that each model converged independently, but the outcomes are remarkably similar in terms of geometry.

In the text

Fig. 6. Synthetic images of β Lyr A for the best-fit nebula model, with total (not reduced) χ 2 = 103 233, shown for four different wavelengths: λ = 155 nm (FUV), 545 nm (V band), 1630 nm (H), and 4750 nm (M). The axes correspond to the right ascension α and declination δ (in mas), while the color scale corresponds to the monochromatic intensity Iv (in erg s−1 cm−2 sr−1 Hz−1). This is a small subset of all 2392 images (per one iteration) used to derive light curves, interferometric visibilities, closure phases, and triple product amplitudes.

In the text

Fig. 7. Observed and synthetic light curves of β Lyr A, shown for all 21 datasets (see their names in the right column). The light curves are phased according to 1. Vertical offsets are arbitrary. The best-fit model is again “nebula” with χ 2 = 103 233; the individual contribution arising from light curves comparison is ${\chi }_{LC}^2 = 6 918$. Synthetic data are denoted by yellow crosses, observed data by blue error bars, and residua by red lines (or circles). There are clear systematic differences especially for datasets iue.1250, iue.1365, oao2.1910, oao2.3320. At the same time, there are neighboring datasets matched relatively well. Sometimes, an intrinsic variability can be also seen (oao2.1430).

In the text

Fig. 8. Similar comparison as in Fig. 7, but for squared visibilities |V|2, with a contribution ${\chi }_{V^2}^2 = 54 137$. The |V|2 values are plotted against projected baseline B/λ (in cycles), and shifted vertically according to the dataset number. The are CHARA/MIRC data at the bottom, NPOI in the middle, and CHARA/VEGA at the top. Synthetic data are denoted by yellow crosses, observed data by blue error bars, and residua by red lines. A few outliers with large uncertainties, which do not contribute much to χ2 anyway, were purposely removed from the plot to prevent clutter. Even though there are some systematic differences for individual segments of data, overall trends seem to be correctly matched.

In the text

Fig. 9. Similar comparison as in Fig. 8, but for closure phases T3ϕ (left) and triple product amplitudes |T3| (right). Contributions to the total χ2 are ${\chi }_{CP}^2 = 29 153$, and ${\chi }_{T_3}^2 = 13 023$. As before, the values are plotted against projected baseline B/λ. T3ϕ measurements were available for NPOI (top half) and CHARA/MIRC (bottom half), while only |T3| from MIRC instrument were used.

In the text

Fig. 10. Top panel: optical depth τ (computed for the wavelength λ = 1630 nm) of the β Lyr A disk, which is observed approximately edge-on (that is from the top). The coordinate z thus corresponds to the line of sight, while x (and y) to the sky plane. The center is empty because the gainer is a non-transparent object. In this case, the inner-rim density of the wedge power-law model is set low, ρ0 = 10−9 g cm−3, to demonstrate a separation from the outer rim. The dashed line shows the physical position of the photosphere (τ = 2/3). The lines-of-sight grid exhibits “steps” due to the limited spatial resolution of our model (1 R⊙), even though the integration of the radiation transfer (the contribution function) is internally performed on a much finer grid. Bottom panel: corresponding photospheric temperature Tphot of the disk, which varies along with the radial temperature profile T(R).

In the text

	Fig. 11. Optical depth τ in the vertical (perpendicular) cross-section of the disk; other parameters are the same as in Fig. 10. The lines of sight are seemingly different from the wedge shape, but this only because they start either in the vacuum, or at the non-transparent object.
In the text

Fig. 12. Comparison of radial temperature profiles of the accretion disk surrounding the gainer. The dotted line is the profile obtained by Mennickent & Djurašević (2013), and the dashed line is mean temperature of the two-temperature model developed by Linnell (2000). The red line represents the best-fitting power-law radial temperature profile (solution slab/pl in Table 9), and the red belt all plausible solutions given by the uncertainty of the inner rim temperature and the exponent of the power-law, and the red point position of photosphere (τ = 2/3), where temperature Tph. = 7 192 K is reached. The blue line represents the best-fitting steady-disk radial temperature profile (solution slab/sd in Table 9), the blue belt all plausible solutions given by uncertainty in the inner rim temperature, and the blue point position of photosphere (τ = 2/3), where temperature Tph: = 7 192 K K is reached. The photospheric temperatures were computed for a line of sight in the plane z = 0 and piercing through center of the accretion disk

In the text

Fig. 13. Synthetic images of slab power-law model which was further improved by a “spot” (that is a spherical structure between the primary and secondary, partly hidden inside the slab). The total (not reduced) χ2 value was decreased from 105 202 (without the spot) down to 102 005, which is a significant improvement. Four wavelengths are shown: λ = 155 nm (FUV), 545 nm (V band), 1630 nm (H), and 4750 nm (M).

In the text

Fig. 14. Comparison of synthetic spectral-energy distributions (thick lines) of β Lyr A and observed SED (thin lines) according to Burnashev & Skulskii (1978). The monochromatic flux Fλ (in erg s−1 cm−2 cm−1) was measured in the range of λ = 330–740 nm. The color scale correspond to the orbital phase. Observation uncertainties are on the order of 5%. The model (corresponding to a “nebula” with tinv = 4.53) was not converged with respect to these observations and the SEDs thus inevitably exhibit some systematic offsets.

In the text

	Fig. A.1. (u, v) coverage of all interferometric observations. Colors correspond to three different instruments: NPOI (blue), CHARA/MIRC (green), CHARA/VEGA (magenta).
In the text

Fig. C.1. Comparison of an observed spectrum of β Lyr B obtained on RJD = 49866.4348 and disentangled spectra of β Lyr B with the best fitting spectra corresponding to their respective solutions listed in Table 3, that is synthetic spectra compared to observed spectrum correspond to slightly different parameters than those compared to disentangled spectra. The blue line denotes disentangled spectrum, the black line the observed spectrum, the red line best-fitting synthetic spectrum, and the yellow line fit residuals. Only principal spectral lines (Si II 6347 Å, Si II 6371 Å, Hα, He I 6678 Å) that are present within the studied spectral region Δλ ≃ 6200–6800 Å and their surroundings are plotted.

In the text

Fig. D.1. Synthetic image of β Lyr A for λ = 155 nm (OAO2 band). This model is based on a “nebula” with a vertical temperature jump, T(z) profile is determined by its scale height H and multiplication factor h mul = 5.27, together with temperature inversion scale h inv = 5.19 and factor t inv = 4.87. The resulting total χ 2 = 99430, with individual contributions as low as ${\chi }_{\text{LC}}^2 = 4 932$, ${\chi }_{V^2}^2 = 53 534$, ${\chi }_{\text{CP}}^2 = 28 282$, ${\chi }_{{\text{T}}_3}^2 = 12 684$, which is significantly better than the nominal nebula model presented in Table 9. However, the artifacts (bright spots close to the outer rim; black in this color scale) are clearly visible.

In the text