© ESO 2019

1. Introduction

HD 81357 (BD+58°1190, MWC 859, HIP 46377, 2MASS J09272389+5808342) is a little studied 79 star classified B8 in the HD catalogue. Merrill & Burwell (1950) included HD 81357 in the third edition of the Mount Wilson Catalogue of B and A emission stars and noted that on one-prism spectrograms of HD 81357 taken in January and December, 1948, Hα was a wide, bright line, possibly double, and superposed on a broad absorption. Allen (1973) obtained near-infrared (IR) photometry of emission-line stars. For HD 81357, he gives only H and K magnitudes and classifies it as “category X”, meaning a star for which the observed continuum distribution corresponds to that expected for its spectral type within the observational uncertainty of the observed magnitudes and the deduced interstellar reddening. Andrillat & Fehrenbach (1982) included HD 81357 in their catalogue of Be stars observed with a 512-pixel television equipment Multiphot¹ attached to the Echelec spectrograph of the Haute Provence Observatory 1.52 m reflector. A spectrum taken on February 13, 1981 showed double-peak Hα emission with I_V = 1.88, and I_R = 1.65. Buscombe (1984) confirmed that HD 81357 had a spectral type B8e. Halbedel (1996) determined its rotational velocity of 150 km s⁻¹ and spectral type B9e. She mentioned HD 81357 under the section entitled “Binaries” and remarked that “HD 81357 also shows a curious spectrum. It seems to be composite with an F/G secondary, though it is possible that it has a very rich assortment of shell lines”. To encourage further study of this object, she added that “it does seem to undergo minor velocity changes”. Koen & Eyer (2002) extracted new candidate periodic variables from the epoch photometry of the HIPPARCOS catalogue. The analysis of 129 measurements for HD 81357 yielded a frequency of 0.03884 d⁻¹ () with an amplitude of 00189, the type of the variability being denoted as U, meaning unsolved. Wheelwright et al. (2010) listed HD 81357 as an Herbig Ae/Be star binary without giving any reference. We were unable to find any report of the presence of a strong IR excess due to dust (a criterion to distinguish the classical and Herbig Be stars) in the bibliography of HD 81357, and note that the star is located outside the zone of the star formation, which is another defining characteristics of the Herbig Be/Ae stars. Koubský et al. (2012), motivated by the note of Halbedel (1996) and by the results of Koen & Eyer (2002), analysed ten spectra of HD 81357 taken in Hα and near-IR regions and measured the radial velocities (RVs) of a number of metallic lines, clearly belonging to a late spectral type (mainly Fe I and Ca II lines). They found that both – the RVs and HIPPARCOS photometry – could be folded with a period of 33.8 days, thus confirming the binary nature of the object. The corresponding semiamplitude of RV curve was 79 km s⁻¹. No lines corresponding to a B spectral type were reliably detected, with the exception of the H I lines, filled by emission. Here we present the first detailed study of the system, based primarily on the new spectral and photometric observations accumulated since Spring 2011 at four observatories.

2. Observations and reductions

Throughout this paper, we specify all times of observations using reduced heliocentric Julian dates

2.1. Spectroscopy

The star was observed at the Ondřejov (OND), and Dominion Astrophysical (DAO) observatories. The majority of the spectra cover the red spectral region around Hα, but we also obtained some spectra in the blue region around Hγ and infrared spectra in the region also used for the Gaia satellite. The journal of spectral observations is in Table 1.

The methods of spectra reductions and measurements were basically the same as in the previous paper 30 of this series devoted to BR CMi (Harmanec et al. 2015). See also Appendix A for details.

2.2. Photometry

The star was observed at Hvar over several observing seasons, first in the UBV, and later in the UBVR photometric system relative to HD 82861. The check star HD 81772 was observed as frequently as the variable. All data were corrected for differential extinction and carefully transformed to the standard systems via non-linear transformation formulae. We also used the HIPPARCOSH_p observations, transformed to the Johnson V magnitude after Harmanec (1998) and recent ASAS-SN Johnson V survey observations (Shappee et al. 2014; Kochanek et al. 2017). Journal of all observations is in Table 2 and the individual Hvar observations are tabulated in Appendix B, where details on individual data sets are also given.

3. Iterative approach to data analysis

The binary systems in phase of mass transfer between the components usually display rather complex spectra; besides the lines of both binary components there are also some spectral lines related to the circumstellar matter around the mass-gaining star (see e.g. Desmet et al. 2010). Consequently, a straightforward analysis of the spectra based on one specific tool (e.g. 2D cross-correlation or spectra disentangling) to obtain RVs and orbital elements cannot be applied. One has to combine several various techniques and data sets to obtain the most consistent solution. This is what we tried, guided by our experience from the previous paper of this series (Harmanec et al. 2015). We analysed the spectroscopic and photometric observations iteratively in the following steps, which we then discuss in the subsections below.

• 1.

  We derived and analysed RVs in several different ways to verify and demonstrate that HD 81357 is indeed a double-lined spectroscopic binary and ellipsoidal variable as first reported by Koubský et al. (2012).

• 2.

  We disentangled the spectra of both binary components in a blue spectral region free of emission lines with the program KOREL (Hadrava 1995, 1997, 2004a) to find out that – like in the case of BR CMi (Harmanec et al. 2015) – a rather wide range of mass ratios gives comparably good fits.

• 3.

  Using the disentangled component spectra as templates, we derived 2D cross-correlated RVs of both components with the help of the program asTODCOR (see Desmet et al. 2010; Harmanec et al. 2015, for the details) for the best and two other plausible values of the mass ratio. We verified that the asTODCOR RVs of the hotter component define a RV curve, which is in an almost exact anti-phase to that based on the lines of the cooler star. A disappointing finding was that the resulting RVs of the hot component differ for the three different KOREL templates derived for the plausible range of mass ratios.

• 4.

  We used the PYTERPOL program kindly provided by J.Nemravová to estimate the radiative properties of both binary components from the blue spectra, which are the least affected by circumstellar emission. The function of the program is described in detail in Nemravová et al. (2016) ².

• 5.

  Keeping the effective temperatures obtained with PYTERPOL fixed, and using all RVs for the cooler star from all four types of spectra together with all photometric observations we started to search for a plausible combined light-curve and orbital solution. To this end we used the latest version of the program PHOEBE 1.0 (Prša & Zwitter 2005, 2006) and tried to constrain the solution also by the accurate Gaia parallax of the binary.

• 6.

  As an additional check, we also used the program BINSYN (Linnell & Hubeny 1994) to an independent verification of our results.

3.1. The overview of available spectra

Examples of the red, infrared and blue spectra at our disposal are in Fig. 1. In all spectra numerous spectral lines typical for a late spectral type are present. In addition to it, one can see several H I lines, affected by emission and a few absorption lines belonging to the hotter component.

Fig. 1. Examples of available blue, red, and infrared spectra. From top to bottom: blue spectrum, red spectrum, enlarged part of a red spectrum near to Hα, and the infrared spectrum. All three regions contain numerous lines of the cool component.

3.2. Light and colour changes

All light curves at our disposal exhibit double-wave ellipsoidal variations with the orbital 338 period. Their amplitude is decreasing from the R to the U band, in which the changes are only barely visible. We show the light curves later, along with the model curves.

In Fig. 2 we compare the colour changes of HD 81357 in the U − B vs. B − V diagram with those known for several other well observed Be binary stars. One can see that its colour changes are remarkably similar to those for another ellipsoidal variable BR CMi while other objects shown are the representatives of the positive and inverse correlation between the light and emission strength as defined by Harmanec (1983) (see also Božić et al. 2013). In particular, it is seen that for KX And dereddened colours exhibit inverse correlation with the object moving along the main-sequence line in the colour-colour diagram from B1V to about B7V. On the other hand, CX Dra seems to exhibit a positive correlation after dereddening, moving from B3 V to B3 I.

Fig. 2. Variations of HD 81357 in the colour – colour diagram are compared to those known for some other Be stars observed at Hvar.

3.3. Reliable radial velocities

3.3.1. phdia RVs

Using the program phdia (see Appendix A), we measured RVs of a number of unblended metallic lines in all blue, red and infrared spectra at our disposal. A period analysis of these RVs confirmed and reinforced the preliminary result of Koubský et al. (2012) that these RVs follow a sinusoidal RV curve with a period of 3377 and a semi-amplitude of 81 km s⁻¹. As recently discussed in paper 30 (Harmanec et al. 2015), some caution must be exercised when one analyses binaries with clear signatures of the presence of circumstellar matter in the system. The experience shows that the RV curve of the Roche-lobe filling component is usually clean (with a possible presence of a Rossiter effect) and defines its true orbit quite well, while many of the absorption lines of the mass-gaining component are affected by the presence of circumstellar matter, and their RV curves do not describe the true orbital motion properly. It is therefore advisable to select suitable spectral lines in the blue spectral region, free of such effects.

3.3.2. KOREL maps

In the next step, we therefore derived a map of plausible solutions using a Python program kindly provided by J. A. Nemravová. The program employs the program KOREL, calculates true χ² values and maps the parameter space to find the optimal values of the semiamplitude K₂ and the mass ratio M₁/M₂. In our application to HD 81357, we omitted the Hγ line and used the spectral region 4373–4507 Å, which contains several He I and metallic lines of the B-type component 1. As in the case of a similar binary BR CMi (Harmanec et al. 2015) we found that disentangling is not the optimal technique how to derive the most accurate orbital solution since the interplay between a small RV amplitude of the broad-lined B primary and its disentangled line widths results in comparably good fits for a rather wide range of mass ratios. The lowest χ² was obtained for a mass ratio of 9.75, but there are two other shallow minima near to the mass ratios 7 and 16. The optimal value of K₂ remained stable near to 81–82 km s⁻¹.

3.3.3. Velocities derived via 2D cross-correlation

As mentioned earlier, the spectrum of the mass-gaining component is usually affected by the presence of some contribution from circumstellar matter, having a slightly lower temperature than the star itself (e.g. Desmet et al. 2010). This must have impact on KOREL, which disentangles the composite spectra on the premise that all observed spectroscopic variations arise solely from the orbital motion of two binary components. Although we selected a blue spectral region with the exclusion of the Hγ line (in the hope to minimize the effect of circumstellar matter), it is probable that KOREL solutions returned spectra, which – for the mass-gaining star – average the true stellar spectrum and a contribution from the circumstellar matter.

As an alternative, we decided to derive RVs, which we anyway needed to be able to combine photometry and spectroscopy in the PHOEBE program, with 2D cross-correlation. We used the asTODCOR program written by one of us (YF) to this goal. The software is based on the method outlined by Zucker & Mazeh (1994) and has already been applied in similar cases (see Desmet et al. 2010; Harmanec et al. 2015, for the details). It performs a 2D cross-correlation between the composite observed spectra and template spectra. The accuracy and precision of such measurements depend on both, the quality of the observations, and on how suitable templates are chosen to represent the contribution of the two stars.

In what follows, we used the observed blue spectra over the wavelength range from 4370 to 4503 Å, mean resolution 0.12 A mm⁻¹, and a luminosity ratio L₂/L₁ = 0.1. We adopted the spectra disentangled by KOREL for the optimal mass ratio 9.75 as the templates for the 2D cross-correlation. We attempted to investigate the effect of circumstellar matter on the RVs derived with asTODCOR for the primary using also the template spectra for the other two mass ratios derived with KOREL. The RVs are compared in Fig. 3. From it we estimate that, depending on the orbital phase, the systematic error due to the presence of circumstellar matter may vary from 0 to 3 km s⁻¹.

Fig. 3. Comparison of asTODCOR RVs of star 1 for the optimal mass ratio q = 9.75 (solid line) and for the two extreme mass-ratios KOREL templates (q = 16 as filled circles, and q = 7 as open circles). Typical errors of asTODCOR individual RVs are close to 1.0 km s−1 – see Table A.2.

The resulting asTODCOR RVs for the optimal mass ratio 9.75 are listed, with their corresponding random error bars, in Table A.2, while the SPEFO RVs of the cooler star are in Table A.3.

Assigning the three sets of asTODCOR RVs from the blue spectral region with the weights inversely proportional to the square of their rms errors, we derived orbital solutions for them. The solutions were derived separately for the hot component 1 and the cool component 2 to verify whether those for the hot star describe its true orbital motion properly. We used the program FOTEL (Hadrava 1990, 2004b). The results are summarised in Table 3. We note that the values of the superior conjunction of component 1 from the separate solutions for component 1 and component 2 agree within their estimated errors, but in all three cases the epochs from component 1 precede a bit those from component 2. This might be another indirect indication that the RVs of component 1 are not completely free of the effects of circumstellar matter.

A disappointing conclusion of this whole exercise is that there is no reliable way how to derive a unique mass ratio from the RVs. One has to find out some additional constraints.

3.4. Trial orbital solutions

In the next step, we derived another orbital solution based on 151 RVs (115 SPEFO RVs, 18 asTODCOR RVs of component 2 and 18 asTODCOR RVs of component 1). As the spectra are quite crowded with numerous lines of component 2, we were unable to use our usual practice of correcting the zero point of the velocity scale via measurements of suitable telluric lines (Horn et al. 1996). That is, why we allowed for the determination of individual systemic velocities for the four subsets of spectra defined in Table 1. All RVs were used with the weights inversely proportional to the square of their rms errors. This solution is in Table 4. There is a very good agreement in the systemic velocities from all individual data subsets, even for the phdia RVs from blue spectra, where only four spectral lines could be measured and averaged.

3.5. Radiative properties of binary components

To determine the radiative properties of the two binary components, we used the Python program PYTERPOL, which interpolates in a pre-calculated grid of synthetic spectra. Using a set of observed spectra, it tries to find the optimal fit between the observed and interpolated model spectra with the help of a simplex minimization technique. It returns the radiative properties of the system components such as T_eff, v sin i or log g, but also the relative luminosities of the stars and RVs of individual spectra.

In our particular application, two different grids of spectra were used: AMBRE grid computed by de Laverny et al. (2012) was used for component 2, and the Pollux database computed by Palacios et al. (2010) was used for component 1. We used the 18 file D blue spectra from OND, which contain enough spectral lines of component 1. The following two spectral segments, avoiding the region of the diffuse interstellar band near to 4430 Å, were modelled simultaneously:

Uncertainties of radiative properties were obtained through Markov chain Monte Carlo (MCMC) simulation implemented within emcee ³ Python library by Foreman-Mackey et al. (2013). They are summarised in Table 5 and an example of the fit is in Fig. 4. We note that PYTERPOL derives the RV from individual spectra without any assumption about orbital motion. It thus represents some test of the results obtained in the previous analysis based on asTODCOR RVs.

Fig. 4. Example of the comparison of an observed blue spectrum in two selected spectral regions with a combination of two synthetic spectra. The residuals in the sense observed minus synthetic are also shown on the same flux scale. To save space, the spectra in the bottom panel were linearly shifted from 1.0 to 0.3. The Hγ emission clearly stands out in the residuals in the first panel. See the text for details.

3.6. Combined light-curve and orbital solution in PHOEBE

To obtain the system properties and to derive the final ephemeris, we used the program PHOEBE 1 (Prša & Zwitter 2005, 2006) and applied it to all photometric observations listed in Table 2 and phdia and asTODCOR RVs for star 2. Since PHOEBE cannot treat different systemic velocities, we used actually RVs minus respective γ velocities from the solution listed in Table 4. For the OND blue spectra (file D of Table 1), we omitted the less accurate SPEFO RVs and used only asTODCOR RVs. Bolometric albedos for star 1 and 2 were estimated from Fig. 7 of Claret (2001) as 1.0 and 0.5, respectively. The coefficients of the gravity darkening β were similarly estimated as 1.0 and 0.6 from Fig. 7 of Claret (1998).

When we tried to model the light curves on the assumption that the secondary is detached from the Roche lobe, we were unable to model the light-curve amplitudes. We therefore conclude that HD 81357 is a semi-detached binary in the mass transfer phase between the components.

It is not possible to calculate the solution in the usual way. One parameter, which comes into the game, is the synchronicity parameter F, which is the ratio between the orbital and rotational period for each component. While it is safe to adopt F₂ = 1.0 for the Roche-lobe filling star 2, the synchronicity parameter F₁ must be re-calculated after each iteration as

(1)

where the equatorial radius is again in the nominal solar radius , the orbital period in days, and the projected rotational velocity in km s⁻¹. We adopt the value of 166.4 km s⁻¹ for v₁sini from the PYTERPOL solution.

It is usual that there is a very strong parameter degeneracy for an ellipsoidal variable. To treat the problem we fixed T_eff of both components obtained from the PYTERPOL solution and used the very accurate parallax of HD 81357 from the second data release of the Gaia satellite (Gaia Collaboration 2016, 2018) to restrict a plausible range of the solutions.

From a PHOEBE light-curve solution for the Hvar UBV photometry in the standard Johnson system it was possible to estimate the following UBV magnitudes of the binary at light maxima

We calculated a number of trial PHOEBE solutions, keeping the parameters of the solution fixed and mapping a large parameter space. For each such solution we used the resulting relative luminosities to derive the UBV magnitudes of the hot star 1, dereddened them in a standard way and interpolating the bolometric correction from the tabulation of Flower (1996) we derived the range of possible values for the mean radius R₁ for the Gaia parallax and its range of uncertainty from the formula

(2)

(Prša et al. 2016). This range of the radius was compared to the mean radius R₁ obtained from the corresponding PHOEBE solution. We found that the agreement between these two radius determinations could only be achieved for a very limited range of mass ratios, actually quite close to the mass ratio from the optimal KOREL solution.

The resulting PHOEBE solution is listed in Table 6 and defines the following linear ephemeris, which we shall adopt in the rest of this study

(3)

The fit of the individual light curves is shown in Fig. 5 while in Fig. 6 we compare the fit of the RV curve of component 2 and also the model RV curve of component 1 with the optimal asTODCOR RVs, which however were not used in the PHOEBE solution.

Fig. 5. UBVR light curves modelled with PHOEBE. Abscissa is labelled with orbital phases according to ephemeris (3).The same scale on ordinate was used to show the changes of the amplitude with passband. For all curves, also the rms errors of individual data points are shown.

Fig. 6. Top panels: RV curve of star 2 and the residuals from PHOEBE solution. The rms errors are comparable to the symbols’ size. Bottom panels: asTODCOR RV curve of star 1, not used in the PHOEBE solution, and its residuals from that solution.

The combined solutions of the light- and RV-curves demonstrated a strong degeneracy among individual parameters. In Table 7 we show the original and normalized χ² values for the individual data sets used. It is seen that the contribution of the photometry and RVs to the total sum of squares are comparable. A higher χ² for RVs might be related to the fact that we were unable to compensate small systematic differences in the zero points of RVs between individual spectrographs and/or spectral regions perfectly, having no control via telluric lines. The degeneracy of the parameter space is illustrated by the fact that over a large range of inclinations and tolerable range of mass ratios the program was always able to converge, with the total sum of squares differing by less than three per cent.

In passing we note that we also independently tested the results from PHOEBE using the BINSYN suite of programs (Linnell 1984; Linnell & Hubeny 1996; Linnell 2000) with steepest descent method to optimise the parameters of the binary system (Sudar et al. 2011). This basically confirmed the results obtained with PHOEBE.

4. Hα profiles

The strength of the Hα emission peaks ranges from 1.6 to 2.0 of the continuum level and varies with the orbital phase (cf. Fig. 8). For all spectra a central reversal in the emission is present. But in many cases the absorption structure is quite complicated.

Fig. 7. SPEFO RVs of the Hα emission wings with their estimated rms errors plotted vs. orbital phase.

Fig. 8. Dynamical spectra in the neighbourhood of the Hα profile in a grey-scale representation created in the program phdia. Left-hand panel: observed spectra; right-hand panel: difference spectra after the subtraction of a synthetic spectrum (including the Hα profile) of star 2. The phases shown on both sides of the figure, correspond to ephemeris (3). The spectra are repeated over two orbital cycles to show the continuity of the orbital phase changes.

The RV curve of the Hα emission wings is well defined, sinusoidal and has a phase shift of some 5 days relative to the RV curve based on the lines of component 1 (see Table 4 and Fig. 7). This must be caused by some asymmetry in the distribution of the circumstellar material producing the emission. In principle, however, one can conclude that the bulk of the material responsible for Hα emission moves in orbit with component 1 as it is seen from the semi-amplitude and systemic velocity of the Hα emission RV curve, which are similar to those of component 1. We note that in both, the magnitude and sense of the phase shift, this behaviour is remarkably similar to that of another emission-line semi-detached binary AU Mon (Desmet et al. 2010).

The orbital variation of the Hα emission-line profiles in the spectrum of HD 81357 is illustrated in the dynamical spectra created with the program phdia – see the left panel of Fig. 8. The lines of the cool component are seen both shortward and longward of the Hα emission. These lines can be used to trace the motion of Hα absorption originating from star 2. It is readily seen that the behaviour of the central absorption in the Hα emission is more complicated than what would correspond to the orbital motion of star 2. Three stronger absorption features can be distinguished: region 1 near to phase 0.0, region 2 visible from phase 0.4 to 0.5, and somewhat fainter feature 3 present between phases 0.65 and 0.85. The absorption in region 2 follows the motion of star 2, while the absorption line in region 3 moves in antiphase.

In the right panel of Fig. 8 we show the dynamical spectra of the difference profiles resulting after subtraction of an interpolated synthetic spectrum of star 2 (properly shifted in RV according to the orbital motion of star 2) from the observed Hα profiles. We note that the lines of star 2 disappeared, but otherwise no pronounced changes in the Hα profiles occurred in comparison to the original ones. Thus, the regions of enhanced absorption 1, 2, and 3, already seen on the original profiles are phase-locked and must be connected with the distribution of circumstellar matter in the system.

There are two principal possible geometries of the regions responsible for the Hα emission: Either an accretion disk, the radius of which would be limited by the dimension of the Roche lobe around the hot star, for our model; or a bipolar jet perpendicular to the orbital plane, known, for instance, for β Lyr (Harmanec et al. 1996), V356 Sgr, TT Hya, and RY Per (Peters & Polidan 2007) or V393 Sco (Mennickent et al. 2012), originating from the interaction of the gas stream encircling the hot star with the denser stream flowing from the Roche-lobe filling star 2 towards star 1 (Bisikalo et al. 2000). In passing we note that no secular change in the strength of the Hα emission over the interval covered by the data could be detected.

5. Stellar evolution of HD 81357 with mass exchange

Given a rather low mass and a large radius of the secondary, we were interested to test whether stellar evolution with mass exchange in a binary can produce a system similar to HD 81357. We use the binary stellar evolution program MESA (Paxton et al. 2011, 2015) in order to test a certain range of input parameters. We tried the initial masses in the intervals , , and the initial binary period P ∈ (2; 10) days. Hereinafter, we use the same notation as in the preceding text, so that M₁ is the original secondary, which becomes the primary during the process of mass exchange. The mass transfer was computed with Ritter (1988) explicit scheme, with the rate limited to , and magnetic breaking of Rappaport et al. (1983), with the exponent γ = 3. For simplicity, we assumed zero eccentricity, conservative mass transfer, no tidal interactions, and no irradiation. We used the standard time step controls.

An example for the initial masses , M₂  =  , the initial period P  =  2.4 d and the mass transfer beginning on the SGB is presented in Fig. 9. We obtained a result, which matches the observations surprisingly well, namely the final semimajor axis , which corresponds to the period P_syn  ≐  32.60 d (while P_obs  =  33.77 d), the final masses (), (), the maximum secondary radius (), with the exception of the primary radius (cf. ). Alternatively, solutions can be also found for different ratios of the initial masses M₁, M₂, and later phases of mass transfer (RGB), although they are sometimes preceded by an overflow. An advantage may be an even better fit of the final R₁, and a relatively longer duration of the inflated R₂ which makes such systems more likely to be observed.

Fig. 9. Long-term evolution of HD 81357 binary as computed by MESA (Paxton et al. 2015). The initial masses were M1 = 1.5 M⊙, M2 = 2.2 M⊙, and the initial period P = 2.4 d. Top: HR diagram with the (resulting) primary denoted as 1 (dashed black line), secondary as 2 (solid orange line) at the beginning of the evolutionary tracks (ZAMS). Middle: the radius R vs. mass M; the corresponding observed values are also indicated (filled circles). Bottom: semimajor axis a of the binary orbit (green) and the radii R1, R2 of the components vs. time t. The observed value a = 63.95 R⊙ and the time t ≐ 7.6 × 108 yr with the maximal R2 are indicated (thin dotted). Later evolutionary phases are also plotted (solid grey line).

Consequently, we interpret the secondary as a low-mass star with a still inflated envelope, close to the end of the mass transfer phase. We demonstrated that a binary with an ongoing mass transfer is a reasonable explanation for both components of the HD 81357 system. A more detailed modelling including an accretion disk, as carried out by Van Rensbergen & De Greve (2016) for other Algols, would be desirable.

6. Discussion

Quite many of similar systems with a hot mass-gaining component were found to exhibit cyclic long-term light and colour variations, with cycles an order of magnitude longer than the respective orbital periods (cf. a good recent review by Mennickent 2017), while others seem to have a constant brightness outside the eclipses. The later seems to be the case of HD 81357, for which we have not found any secular changes. It is not clear as yet what is the principal factor causing the presence of the cyclic secular changes (see also the discussion in Harmanec et al. 2015). A few yellow-band photoelectric observations of HD 81357 found in the literature all give 83 to 84. It is notable, however, that both magnitudes given in the HD catalogue are 79. Future monitoring of the system brightness thus seems desirable.

We note that we have not found any indication of a secular period change and also our modelling of the system evolution and its low mass ratio M₂/M₁ seem to indicate that HD 81357 is close to the end of mass exchange process. The same is also true about BR CMi studied by Harmanec et al. (2015). We can thus conjecture that the cyclic light variations mentioned above do not occur close to the end of the mass exchange.

Acknowledgments

We gratefully aknowledge the use of the latest publicly available versions of the programs FOTEL and KOREL written by P. Hadrava. Our sincere thanks are also due to A. Prša, who provided us with a modified version of the program PHOEBE 1.0 and frequent consultations on its usage, and to the late A. P. Linnell for his program suite BINSYN. We also thank J. A. Nemravová, who left the astronomical research in the meantime, for her contributions to this study and for the permission to use her Python programs, PYTERPOL and several auxiliary ones, and to C. S. Kochanek for the advice on the proper use of the ASAS-SN photometry. Suggestions and critical remarks by an anonymous referee helped us to re-think the paper, improve and extend the analyses, and clarify some formulations. A. Aret, Š. Dvořáková, R. Křiček, J. A. Nemravová, P. Rutsch, K. Šejnová, and P. Zasche obtained a few Ondřejov spectra, which we use. The research of PK was supported from the ESA PECS grant 98058. The research of PH and MB was supported by the grants GA15-02112S and GA19-01995S of the Czech Science Foundation. HB, DR, and DS acknowledge financial support from the Croatian Science Foundation under the project 6212 “Solar and Stellar Variability”, while the research of D. Korčáková is supported from the grant GA17-00871S of the Czech Science Foundation. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. We gratefully acknowledge the use of the electronic databases: SIMBAD at CDS, Strasbourg and NASA/ADS, USA.

References

Appendix A: Details of the spectral data reduction and measurements

The initial reduction of all OND and DAO spectra (bias subtraction, flat-fielding, creation of 1D spectra, and wavelength calibration) was carried out in IRAF. Optimal extraction was used. Normalization, removal of residual cosmics and flaws and RV measurements of Hα emission profiles were carried out with the program SPEFO (Horn et al. 1996; Škoda 1996), namely the latest version 2.63 developed by Mr. J. Krpata (Krpata 2008). SPEFO displays direct and flipped traces of the line profiles superimposed on the computer screen that the user can slide to achieve a precise overlapping of the parts of the profile of whose RV one wants to measure. RVs of a selection of stronger unblended lines of the cool component 2 (see Table A.1) covering the red and infrared spectral regions (available for all spectra), were measured using the program phdia written by LK. The program phdia is a web application which allows the creation of dynamical spectra and, using a similar principle as the program SPEFO, also RV measurements based on the comparison of direct and flipped line profiles. The input is a normalized digital spectrum. A robust mean RV with the rms error was derived for each spectrum to obtain the mean RV. In contrast to our usual approach, we have not measured a selection of good telluric lines to obtain an additional fine correction of the RV zero point of each spectrogram because of a strong blending of the telluric lines with the stellar lines of the cool component 2. Moreover, we measured the RVs of the steep wings of the Hα emission line repeatably at different times to obtain an estimate of the rms error for these measurements. We point out that although some broad and shallow lines of component 1 are seen in the spectra, their direct RV measurement is impossible because of numerous blends with the lines of component 2.

Appendix B: Details on the photometric data used

Below, we provide some details of the individual data sets and their reductions.

• Station 01 – Hvar: These differential UBV and later UBVR observations have been secured by HB, PH, and DR relative to HD 82861 (the check star HD 81772 being observed as frequently as the variable) and carefully transformed to the standard UBV(R) system via non-linear transformation formulæ using the HEC22 reduction program – see Harmanec et al. (1994) and Harmanec & Horn (1998) for the observational strategy and data reduction⁴. All observations were reduced with the latest HEC22 rel.18 program, which allows the time variation of linear extinction coefficients to be modelled in the course of observing nights. For the light-curve solutions we used normal points averaged over the typical observing sequence of 0055 and the corresponding rms errors.

• Station 61 – HIPPARCOS: These all-sky observations were reduced to the standard V magnitude via the transformation formulæ derived by Harmanec (1998) to verify that no secular light changes in the system were observed. However, for the light-curve solution in PHOEBE, we consider the HIPPARCOS transmission curve for the H_p magnitude and also use the original rms errors.

• Station 93 – ASAS-SN: These all-sky automated survey for supernovae V observations were adopted from the data server https://asas-sn.osu.edu/ (Shappee et al. 2014; Kochanek et al. 2017) and cleaned for some strongly deviating data points. We used only the observations from the bb camera, which were numerous enough. The rms errors provided by the on the fly calculator are unrealistically small (cf. Jayasinghe et al. 2019) and we used them only to assign relative weights to individual observations, inversely proportional to their square and estimated the mean rms scatter as 0.0065 mag.

Journal of all data sets is in Table 2.

All Tables

All Figures

	Fig. 1. Examples of available blue, red, and infrared spectra. From top to bottom: blue spectrum, red spectrum, enlarged part of a red spectrum near to Hα, and the infrared spectrum. All three regions contain numerous lines of the cool component.
In the text

	Fig. 2. Variations of HD 81357 in the colour – colour diagram are compared to those known for some other Be stars observed at Hvar.
In the text

	Fig. 3. Comparison of asTODCOR RVs of star 1 for the optimal mass ratio q = 9.75 (solid line) and for the two extreme mass-ratios KOREL templates (q = 16 as filled circles, and q = 7 as open circles). Typical errors of asTODCOR individual RVs are close to 1.0 km s−1 – see Table A.2.
In the text

Fig. 4. Example of the comparison of an observed blue spectrum in two selected spectral regions with a combination of two synthetic spectra. The residuals in the sense observed minus synthetic are also shown on the same flux scale. To save space, the spectra in the bottom panel were linearly shifted from 1.0 to 0.3. The Hγ emission clearly stands out in the residuals in the first panel. See the text for details.

In the text

	Fig. 5. UBVR light curves modelled with PHOEBE. Abscissa is labelled with orbital phases according to ephemeris (3).The same scale on ordinate was used to show the changes of the amplitude with passband. For all curves, also the rms errors of individual data points are shown.
In the text

	Fig. 6. Top panels: RV curve of star 2 and the residuals from PHOEBE solution. The rms errors are comparable to the symbols’ size. Bottom panels: asTODCOR RV curve of star 1, not used in the PHOEBE solution, and its residuals from that solution.
In the text

	Fig. 7. SPEFO RVs of the Hα emission wings with their estimated rms errors plotted vs. orbital phase.
In the text

Fig. 8. Dynamical spectra in the neighbourhood of the Hα profile in a grey-scale representation created in the program phdia. Left-hand panel: observed spectra; right-hand panel: difference spectra after the subtraction of a synthetic spectrum (including the Hα profile) of star 2. The phases shown on both sides of the figure, correspond to ephemeris (3). The spectra are repeated over two orbital cycles to show the continuity of the orbital phase changes.

In the text

Fig. 9. Long-term evolution of HD 81357 binary as computed by MESA (Paxton et al. 2015). The initial masses were M1 = 1.5 M⊙, M2 = 2.2 M⊙, and the initial period P = 2.4 d. Top: HR diagram with the (resulting) primary denoted as 1 (dashed black line), secondary as 2 (solid orange line) at the beginning of the evolutionary tracks (ZAMS). Middle: the radius R vs. mass M; the corresponding observed values are also indicated (filled circles). Bottom: semimajor axis a of the binary orbit (green) and the radii R1, R2 of the components vs. time t. The observed value a = 63.95 R⊙ and the time t ≐ 7.6 × 108 yr with the maximal R2 are indicated (thin dotted). Later evolutionary phases are also plotted (solid grey line).

In the text