© ESO 2020

1 Introduction

Various types of oscillating stars reside in binary systems and therefore a precise determination of their system and stellar parameters can be performed. In particular the combined photometric-spectroscopic analysis of eclipsing binaries (EBs) leads to the determination of absolute masses and radii of the components in a direct way (e.g. Lehmann et al. 2013; Frandsen et al. 2013). Thanks to spectral disentangling techniques (e.g. Simon & Sturm 1994; Hadrava 1995, Ilijic et al. 2001; Sablowski et al. 2019), faint components giving rise to a contribution of a few per cent can be detected in high-resolution spectra (e.g. Lehmann et al. 2013; Tkachenko et al. 2014; Borkovits et al. 2014). Such techniques ensure the determination of precise (up to one per cent) model-independent dynamical masses of stars, which can further be confronted with asteroseismic values provided at least one of the binary components pulsates. This also means that the results can be used to calibrate the asteroseismic mass determination that becomes more and more important, in particular for planet hosting stars (e.g. Hatzes & Zechmeister 2007, 2008; Hatzes et al. 2012).

Algol-type systems (Algols, hereafter) are semi-detached, interacting EBs consisting of a main-sequence star of spectral type B-A (primary, hereafter) and an evolved F-K type companion (secondary, hereafter). As a consequence of the above-mentioned facts, studying pulsating Algols from complementary spectroscopic and photometric data provides a valuable test of stellar evolutionary models. The group of oscillating eclipsing Algol stars (oEAs, hereafter; Mkrtichian et al. 2002) consists of eclipsing Algols with mass transfer where the mass-accreting primary shows δ Sct-like oscillations. These stars are extraordinary objects for asteroseismic studies because they allow us to investigate short-term dynamical stellar evolution during mass-transfer episodes, most probably caused by the magnetic activity cycle of the less massive secondary. Basic principles of the interaction between the magnetic cycle of the cool secondary, the occurrence of rapid mass-transfer episodes, the dynamical behaviour of the system, and the excitation of different non-radial pulsation (NRP) modes of the mass-gaining primary can be studied in great detail. Mkrtichian et al. (2018) showed for the first time that mass transfer and accretion influences amplitudes and frequencies of NRP modes of the oEA star RZ Cas, where amplitude changes are caused by the sensitivity of the mode selection mechanism to conditions in the outer envelope (e.g. Pamyatnykh et al. 1998)and frequency variations by the acceleration of the outer layers by mass and angular momentum transfer.

The details of angular momentum exchange in Algols determining their evolution are still not completely understood. This is in particularvalid for systems like the so-called R CMa stars (e.g. Budding & Butland 2011; Lehmann et al. 2013, 2018) that include companions of extremely low masses. But also for the oEA star RZ Cas, calculations showed that its actual configuration cannot be explained when assuming a purely conservative mass-transfer scenario in the past (Mennekens et al. 2008). In several Algols, such as RZ Cas (Lehmann & Mkrtichian 2008) or TW Dra (Zejda et al. 2008; Tkachenko et al. 2010), orbital period variations were observed that can be attributed to periods of rapid mass exchange. The gas stream hits the equatorial zones of the atmosphere of the pulsating star, transfers an essential amount of angular momentum, and forces the acceleration of its outermost surface layers, thus causing strong differential rotation. While rotation alters the frequencies of the individual axisymmetric modes (e.g. Saio 1981; Lovekin & Deupree 2008), it also lifts the degeneracy in frequency for the non-axisymmetric modes as observed for some β Cep variables (e.g. Aerts et al. 2004; Jerzykiewicz et al. 2005). This means that NRP modes are sensitive to the acceleration of the surface layers and that it is possible to probe the acceleration via the rotational splitting effect, as suggested by Mkrtichian et al. (2018) for RZ Cas. The study of the corresponding frequency shifts, together with the direct measurement of changes in the projected equatorial velocity (v sin i, hereafter) of the primary from its spectral line profiles leads to an estimation of the amount of matter transferred to the primary. The results can be compared to those obtained from 3D hydrodynamical calculations, assuming different rates of mass transfer (e.g. Mkrtichian et al. 2007), and finally explain the kind of mass and angular momentum transfer of the Algols.

One further advantage of oEAs is that the so-called spatial filtration or eclipse mapping effect, occurring as a result of the obscuration of parts of the stellar disc of the oscillating primary by the secondary during eclipses, simplifies the identification of the pulsation modes. The effect was predicted and found in photometric (e.g. Gamarova et al. 2003; Mkrtichian et al. 2004; Reed et al. 2005) and spectroscopic (Lehmann et al. 2018) observations. The dynamic eclipse mapping method was introduced by Bíró & Nuspl (2011) with the aim of NRP mode identification by reconstructing the surface intensity patterns on EBs. With modern methods such as the least-squares deconvolution (LSD) technique (Donati et al. 1997) and the pixel-by-pixel method of the FAMIAS program (Zima 2008), high-l-degree NRP modes can also be detected and identified from the signal-to-noise ratio (S/N) enhanced line profiles; a unique identification however is difficult (e.g. Lehmann et al. 2009).

RZ Cas (spectral type A3 V + K0 IV) is a short-period (P = 1.1953 d) Algol and one of the best studied oEA stars. A partial eclipse is observed during primary minimum (Narusawa et al. 1994). The primary was found by Ohshima et al. (1998, 2001) to exhibit short-period light variability with a dominant oscillation mode of a frequency of 64.2 c d⁻¹. This finding was later confirmed from dedicated multi-site photometric campaigns by both Mkrtichian et al. (2003) and Rodríguez et al. (2004). The latter authors obtained simultaneous Strömgren uvby light curves of RZ Cas. These authors present a detailed photometric analysis for both binarity and pulsation, deriving WD (Wilson & Devinney 1971) solutions in the four mentioned passbands as well as absolute parameters of both components, and confirming the pulsational behaviour of the primary component as found by Ohshima et al. (1998). We use the radii and flux ratios between the components derived by Rodríguez et al. (2004) as a reference in our spectroscopic investigation. A comprehensive overview on the observations and analysis of RZ Cas can be found in Mkrtichian et al. (2018).

Starting our spectroscopic investigation of RZ Cas in 2001, we were the first to detect rapid oscillations in its spectra (Lehmann & Mkrtichian 2004), and later on in spectra taken in 2006 (Lehmann & Mkrtichian 2008; Tkachenko et al. 2009). From the different amplitudes of the Rossiter-McLaughlin effect (RME hereafter; Rossiter 1924; McLaughlin 1924) observed during the primary eclipses in different seasons and the modelling of line profiles over the full orbital cycle using the Shellspec07_inverse program, we deduced that RZ Cas was in an active phase of mass transfer in 2001, whereas in 2006 it was in a quiet state. To model the surface intensity distribution of the secondary of RZ Cas, we had to include a large cool spot facing the primary, presumably originating from a cooling mechanism by the enthalpy transport via the inner Lagrangian point as suggested by Unno et al. (1994). Comparing the rapid oscillations found in the radial velocities (RVs) from 2001 and 2006 also with those derived from the light curves of RZ Cas taken over many years (see Mkrtichian et al. 2018), we found that the NRP pattern of RZ Cas changed from season to season. Different NRP modeshave been excited with different amplitudes in different years and also frequency variations of single modes were observed. Mkrtichian et al. (2018) suggested that these frequency variations could be caused by a temporary acceleration of the outer layers of the primary owing to angular momentum exchange by mass-transfer effects.

After observing RZ Cas in 2001 and 2006, we took new time series of high-resolution spectra in 2008 and 2009. The fact that we found a typical timescale of about nine years from the behaviour of the pulsation amplitudes but also from light-curve analysis (see Mkrtichian et al. 2018) forced us to start a spectroscopic monitoring of the star covering the years 2013 to 2017. We now investigate the complete data set with the aim to use the spectra and observed variations in RVs and line profile variations (LPV) to deduce stellar and system parameters, check for timely variable NRP pulsation patterns, and try to correlate all observed variations with the occurrence of quiet and active phases of the Algol system. Observations are described in Sect. 2. The spectra taken with the HERMES spectrograph are used for a detailed analysis with the aim to derive precise atmospheric parameters for both components of RZ Cas (Sect. 3). The extraction of RVs and calculation of mean, high-S/N line profiles with the newly developed LSDbinary program is described in Sect. 4. Using different methods, we measure the projected equatorial rotation velocity of the primary (Sect. 5). In Sect. 6 we use the O–C values collected from literature to compute the orbital period variations of RZ Cas over the last decades. We check for non-Keplerian effects in the orbital RV curves in Sect. 7.2 and try to model them using the PHOEBE (Prša & Zwitter 2005) program. The results are discussed in Sect. 8 followed by concluding remarks in Sect. 9.

Our investigation of NRP is based on high-frequency oscillations in the RVs and in LPV. Applied methods and results will be described in a forthcoming article (Paper II) and discussed together with the results presented in this work.

2 Observations

Spectra were taken over a total time span of 16 yr with the TCES spectrograph at the 2 m Alfred Jensch Telescope of the Thüringer Landessternwarte (TLS) Tautenburg and over one year with the HERMES spectrograph (Raskin et al. 2011) at the 1.25 m Mercator Telescope on La Palma. The TCES instrument is an echelle spectrograph in Coude focus and covers the wavelength range 4450–7550 Å. A 2015 upgrade resulted in improved efficiency so that this spectrograph could be used with higher spectral resolution (narrower entrance slit of one instead of two arcsec width). Table 1 gives the journal of observations.

The TCES spectra were reduced using standard MIDAS packages for Echelle spectrum reduction and HERMES spectra were reduced using the HERMES spectrum reduction pipeline. We used our own routines for the normalisation of the spectra to the local continuum. Instrumental shifts were corrected using an additional calibration based on a large number of telluric O₂ absorption lines.

3 Spectrum analysis

3.1 Spectra

The HERMES spectra taken in 2009 comprise the H_ɛ to H_α range, whereas the TLS spectra only include H_β and H_α. Moreover, the HERMES spectra provide the higher spectral resolution. That is why we decided to use them for spectrum analysis. We were able to decompose the spectra into the spectra of the components using the KOREL program (Hadrava 1995). Analysing the decomposed spectra, we faced problems with the continuum normalisation. As a result of the very faint contribution of the secondary, small deviations in the continuum of the composite spectrum from the true continuum leads to large deviations in the continuum of the decomposed spectrum of the secondary. A first normalisation was applied during spectrum reduction by fitting higher polynomials to nodes assumed to represent the local continuum in each of the extracted echelle orders. This approach is not accurate enough when the contribution of the secondary is as faint as in the present case. Corrections to the local continuum have to be applied during spectrum analysis by a comparison with the continua of synthetic spectra. The spectra are then renormalised by multiplying them with the local ratio of the continua of the synthetic spectrum to that of the reduced spectrum. The spectrum decomposition with KOREL, on the other hand, is based on Fourier transform and introduces undulations in the continuum that are additive and have to be subtracted. This introduces an ambiguity with the multiplicative corrections described before that we could not solve to reach the required accuracy. Instead, we used the average of nine composite HERMES spectra of RZ Cas taken at maximum separation of the components and having about the same RV. In this case, we have to apply multiplicative continuum corrections only, for which we used spline functions. Because the signal of the secondary gets fainter for smaller wavelengths and because of the occurrence of telluric lines in the red part of the spectra we limited the spectral range to 4000–5550 Å, including the three Balmer lines H_β, H_γ, and H_δ.

Fig. 1 Continuum flux ratios between the components vs. wavelength. The dotted lines are calculated from synthetic continuum spectra with Teff1 of 8700 K and Teff2 of 3800–5000 K. The best agreement with the flux ratio from photometry (red line) is obtained for Teff2 = 4400 K (black solid line). The upper green line represents Teff1 = 9000 K, Teff2 = 4500 K, the lower green line for Teff1 = 8400 K, Teff2 = 4300 K.

3.2 Methods

We used the GSSP program (Tkachenko 2015) to derive the atmospheric parameters of the components of RZ Cas. It is based on the spectrum synthesis method and performs a grid search in stellar parameters. Synthetic spectra were computed with SynthV (Tsymbal 1996), based on a library of atmosphere models computed with LLmodels (Shulyak et al. 2004) for the hot primary andon MARCS models (Gustafsson et al. 2008) for the cool secondary. Atomic data were taken from the VALD database (Kupka et al. 2000).

Besides the atmospheric parameters, GSSP also has to solve for the a priori unknown flux ratio of the components. In the following, we use the continuum flux ratio of the secondary (component 2) to primary (component 1) for comparison. It is interpolated from the uvby luminosities provided by Rodríguez et al. (2004) and shown by the red line in Fig. 1. To get an impression of the influence of the T_eff of primary andsecondary on the resulting flux ratios with wavelength, we computed synthetic continuum spectra for different T_eff of the components, assuming log g₂ = 3.7 (Rodríguez et al. 2004), and log g₁ = 4.4 and [Fe/H] = −0.42 as derived in Sect. 3.3. Synthetic flux ratios were computed from the ratio of the continuum spectra assuming a radii ratio of R₂ ∕R₁ = 1.17, also taken from Rodríguez et al. (2004). The best representative theoretical curve is shown by the black line in Fig. 1 and corresponds to T_eff1 = 8700 K, T_eff2 = 4400 K. From the green lines, we see that a variation of T_eff2 by 100 K requires a variation of T_eff1 by 300 K to approximately fit the photometric flux ratio.

We applied three different methods (labelled a, b, and c in what follows) to determine the atmospheric parameters together with the radii ratio and continuum flux ratio of the secondary to primary.

(a) We use the standard method of the GSSP program that takes the wavelength dependence of the flux ratio from the ratio of the synthetic continuum spectra calculated with SynthV and scales it with the square of the radii ratio. The latter is obtained from comparing the observed with the synthetic normalised spectra on a grid of atmospheric parameters. To obtain the atmospheric parameters together with the radii ratio, we minimise (1)

where o(λ) and s(λ) are the observed and synthetic composite spectra, respectively, and σ is the estimated mean error of o(λ). The synthetic spectrum is computed from (2)

where s₁ and s₂ are the synthetic spectra of the components shifted for their RVs v₁, v₂ (all spectra in line depths), and the flux ratio (3)

is the product of the squared radii ratio and the ratio of the continuum fluxes C₂ and C₁ per unit surface, determined with SynthV.

(b) We use Eqs. (1) to (3) as before, but fix the radii ratio to 1.17 as obtained from the uvby photometry.

(c) We replace the flux ratio V_F computed so far from Eq. (3) by the flux ratio obtained from the uvby photometry.No synthetic spectra are needed in this case. The best-fitting radii ratio can then be determined from Eq. (3) using a simple least-squares fit.

Because of the known degeneracy between various free parameters, we reduced their number. The spectrum of the secondary with its small contribution to total light suffers from low S/N and we fixed its log g to 3.7, as derived from light curve analysis (Rodríguez et al. 2004). Furthermore, we used the same values for its elemental abundances as derived for the primary, except for the iron and carbon abundances, which we determined separately. For the primary, we basically know the photometric value of log g = 4.33(3) from Rodríguez et al. (2004). Because of the good S/N of the spectra and the fact that the Balmer lines shapes are very sensitive to log g in this temperature range, we decided to use logg of the primary as a free parameter to compare the results with the photometric results.

The GSSP allows us to iterate atmospheric parameters such as T_eff, log g, and microturbulent velocity v_turb together with v sin i and the surface abundance of one chemical element at the same time. We started, in this way, to optimise [Fe/H] together with the other parameters, first for the primary, then for the secondary, and repeated until the change in parameters and χ² became marginal. In a next step, we fixed all parameters to the values obtained so far and optimised the elemental abundances of all other chemical elements for which significant contributions in the spectra could be found.

3.3 Results

Tables 2 and 3 list the results. in Table 2 is the radii ratio obtained from the analysis using methods a, b, or c. is the adjusted radii ratio that follows when we apply the least-squares fit based on Eq. (3), as used in method c for the two other methods as well. This means that when keeping all other parameters obtained with the various methods, the radii ratio has to be changed from to to give the best agreement with the photometric flux ratio. For method c, both radii ratios are identical. The quantity is used in the following as a measure for the agreement of the spectroscopic with the photometric results. It should be equal to 1.17 in the optimum case.

Figures 2–4 compare the results from the three methods. These figures show in the first row the observed spectrum together with the best-fitting synthetic composite spectrum. The second row compares the “observed” spectrum of the secondary with the best-fitting synthetic spectrum found for the secondary. The observed spectrum of the secondary was therefore computed by subtracting the best-fitting spectrum of the primary from the observed composite spectrum. These spectra are rescaled according to the obtained flux ratio and the observed spectrum of the secondary is correspondingly noisy. The bottom row compares the obtained flux ratio as a function of wavelength with that obtained from photometry. The differences seen by eye are marginal.

From Table 2 and Figs. 2–4 we conclude as follows: First, there are only small differences in the atmospheric parameters and elemental abundances derived with the different methods for the primary. Second, methods a and c give T_eff of the secondary that agree with each other within the 1 σ error bars, but both are distinctly higher than the value expected from photometry. The T_eff derived from method b) on the other hand, is in agreement with that value. Third, the iron and carbon abundances of the secondary cannot be distinguished from those of the primary within the 1σ errors (Table 2), but all three methods yield a distinctly lower [C/H] of the primary component compared to the solar value. Fourth, all three methods deliver flux ratios higher than the photometric ratio, as can be seen from the bottom panels in Figs. 2–4. Based on the atmosphericparameters derived with methods a and c, RZ Cas should have much smaller radii ratios (the in Table 2) than computed. From the results of light curve analysis by Rodríguez et al. (2004), however, we expect a radii ratio of 1.17 and can exclude that the radius of the secondary is larger than that of the primary. From method b, on the other hand, we obtain a value of that is larger than unity and not so far from 1.17.

There is a large degeneracy between various parameters such as T_eff, log g, v_turb, [Fe/H] (both components), [Mg/H] (cool component), and radii ratio. This explains why we end up in both methods a and c with extraordinary small radii ratios on costs of higher T_eff of the secondary (much too high when comparing with Fig. 1). We assume that the number of degrees of freedom is too high when trying to optimise the radii ratio together with the other parameters; the signal from the faint companion in our composite spectra is simply too low for that. Thus we observe, although it does not give the smallest root mean square (rms) of the O-C residuals, that the most reliable results are obtained when fixing the flux ratio to the photometric ratio as we did in method b.

Fig. 2 Results of spectrum analysis using method a (cf. Sect. 3.3, Table 2), shown for the Hδ -Hγ region (left) and the region around the Mg I b triplet (right). First row: observed, continuum adjusted composite spectrum (black), best-fitting synthetic spectrum (red), and shifted difference spectrum (green). The inverse of the applied continuum correction is shown in blue. Second row: same for the secondary. Third row:flux ratio of the secondary to primary from photometry (black) from spectrum analysis based on R2 ∕R1 = 0.926 (red) and assuming the adjusted ratio of 0.78 (green).

Fig. 3 As Fig. 2, but using method b. Third row: flux ratio from photometry (black) and from spectrum analysis based on R2∕R1=1.17 (red) and on the calculated ratio of 1.065 (green).

Fig. 4 AsFig. 2, but using method c. Third row: flux ratio from photometry (black) from spectrum analysis based on the photometric flux ratio assuming R2 ∕R1 = 1.17 (red) and on the adjusted ratio of R2∕R1 = 0.87 (green). Bottom row: ratio of the photometric flux ratio to the flux ratio obtained from spectrum analysis based on R2 ∕R1 = 1.17.

4 Radial velocities and LSD profiles with LSDbinary

The classical method of LSD (Donati et al. 1997) is based on using one line mask as a template. In the result, a strong deconvolved line profile is obtained for the star that best matches this template, whereas the contribution of the other component is more or less suppressed. Tkachenko et al. (2013) generalised the method so that it can simultaneously compute an arbitrary number of LSD profiles from an arbitrary number of line masks. In this work, we used the LSDbinary program written by V. Tsymbal, which computes separated LSD profiles for the two stellar components using as templates two synthetic spectra that are based on two different atmosphere models. Theprogram also delivers optimised values of the RVs of the components and their radii ratio. A first successful application of the program to the short-period Algol R CMa is presented in Lehmann et al. (2018); this work also includes a short description of the program algorithms and shows the advantages of LSDbinary against TODCOR (Zucker & Mazeh 1994) in the case of very small flux ratios between the components of binary stars.

We applied LSDbinary to the spectra of RZ Cas to obtain the separated LSD profiles and RVs of its components. Figure 5 shows LSD profiles computed from RZ Cas spectra taken during the primary eclipse in 2016 as an example. The RME can clearly be seen in the profiles of the primary, as well as the RV variation due to orbital motion in the profiles of both components. We will use the obtained RVs and LSD profiles for a detailed investigation of the stellar and system parameters of RZ Cas in this first part and of the pulsations of its primary component in Paper II.

Fig. 5 LSD profiles computed with LSDbinary from spectra taken during primary eclipse for the primary (left) and secondary (right) of RZ Cas. Phase zero corresponds to Min I.

5 Rotation velocity of the primary

We applied three different methods to determine v sin i of the primary of RZ Cas. In all cases, we only use spectra taken in out-of-eclipse phases.

5.1 Fourier method

First, we applied the Fourier method (Carroll 1933; Smith & Gray 1976; Gray 2005) to the calculated LSD profiles. This method (DFT hereafter) is based on the determination of the rotation-broadening related zero points in the Fourier power spectra of the profiles. We assumed a linear limb darkening law with a limb darkening coefficient β = 1.5 (or b = 0.6 with b = β∕(1 + β)). Varying the value of β leads to slight systematic effects in v sin i, but changes were minor compared to the differences to the results from the two other methods described below. The selection of zero points was based on a 3σ clipping, comparing the v sin i from single zero points with the mean values from all zero points.

In the result, we obtained v sin i values strongly varying with orbital phase. Figure 6 shows an example. We assume that this variation comes from Algol-typical effects such as an inhomogeneous circum-primary gas-density distribution and/or surface structures caused by the influence of mass transfer and gas stream from the secondary. In a next step, we removed all the variations found in the line profiles with the pixel-by-pixel method of FAMIAS (Zima 2008) by subtracting all frequency contributions computed with the mentioned program from the profiles. This method will be explained in detail in Paper II. For each season, we considered a certain pixel of all LSD profiles as part of one time series from which we subtracted the found contributions and did this pixel-by-pixel to build the undistorted profiles. We performed that for all seasons but 2008 for which we did not have enough data to apply the pixel-by-pixel method. The subtraction of the high-amplitude, low-frequency contributions found with FAMIAS (most of them are harmonics of the orbital frequency) are responsible for the cleaning, not the faint high-frequency oscillations due to pulsation. The resulting v sin i are shown in Fig. 6 in red. The distribution is now much flatter and we used it to calculate the v sin i of the corresponding season as its arithmetic mean. Results are listed in the third column of Table 4.

5.2 Single spectral line

Next, we looked for stronger spectral lines of the primary in the composite spectra that are free of blends from the cool companion. We found only one line consisting of the Fe II doublet 5316.61/5316.78 Å that fulfils that condition. For each season, we fitted the line profiles by synthetic spectra computed with the parameters listed in Table 2, method b to determine the best-fitting v sin i and its error. This means that we were fixing all atmospheric parameters except for v sin i to the solution determined from the spectra observed in 2009 when the star was in a relatively quiet phase. However, this approach does not account for the effects in active phases of RZ Cas such as the attenuation of light by circumstellar material as found for the year 2001 for example (Tkachenko et al. 2009). Thus we used two more free parameters in our fit, correction factor a counting for different line depths caused by the presumed effects, and factor b to correct the continuum of the observed spectrum in the vicinity of the Fe II line. Both were determined from a least-squares fit (4)

where P_s is the synthetic and P_o the observed line profile. Finally, we built the mean v sin i per season from the weighted mean of well-selected data points using 3σ clipping. Figure 7 shows two examples: one for 2001 in which the star was in an active phase, and one for 2014 in which the v sin i shows a much smoother behaviour with orbital phase. The results are listed in the fourth column of Table 4.

5.3 Using FAMIAS

In Paper II, we will use the moment methods (Aerts et al. 1992) of the FAMIAS program for a mode identification of low-l degree modes. For that, we cleaned the observed LSD profiles for all low-frequency distributions in the same way as described in Sect. 5.1. One free parameter in applying the moment method is the v sin i of the primary, which optimum value and 1σ error we obtained from the resulting χ²-distribution. The values are listed in the last column of Table 4. The number of observations in 2008 were not sufficient to apply this method.

Fig. 6 Values of v sin i obtained from DFT vs. out-of-eclipse phases measured from the LSD profiles observed in 2016 (black) and from the same profiles corrected for the LPV found with FAMIAS (red).

Fig. 7 Values of v sin i determined from the FeII 5317 Å line vs. orbital phase, shown for 2001 and 2014. The mean values were built from the values indicated in green.

Fig. 8 Values of v sin i measured from (a) DFT (black crosses), (b) the Fe II line (red squares), and (c) FAMIAS (green plus signs). The solid and dotted lines show possible sinusoidal variations (see text).

5.4 Comparison

The v sin i determined with the three different methods are shown in Fig. 8. The values obtained from DFT and FAMIAS agree in most cases within the 1σ error bars; there is a larger difference for 2006. A systematic offset can be observed for the v sin i values obtained from the Fe II line. This is in particular the case for 2001 when RZ Cas was in an active phase. We assume that the offset comes from the cleaning procedures that we applied when using the other two methods. In these cases we removed the low-frequency contributions that we found with the pixel-by-pixel method of FAMIAS from the LSD profiles so that we removed theline broadening effects due to orbital-phase dependent dilution effects by circumstellar material in this way.

Finally, we built a spline fit through the averaged values obtained from DFT and FAMIAS and interpret this as the typical variation of v sin i over the epoch of our observations. We see that v sin i was distinctly larger in 2001 compared to the other seasons. We assume that its value increased during the active phase of RZ Cas near to 2001 because of the acceleration of the outer layers of the primary by angular momentum transport via mass transfer from the cool component (see Sect. 8 for a more detailed discussion). A period search in the values averaged from DFT and FAMIAS gives a period of 9.6 yr, overlaid on a long-term trend (the solid line in Fig. 8). A fit based on the 9.0 yr period as discussed in Sect. 7.2.3 is shown by the dotted line. The difference is marginal. The last column in Table 4 lists the rotation-to-orbit synchronisation factor of the primary, F₁, calculated from its radius, taken from Rodríguez et al. (2004) as 1.67 ± 0.02 R_⊙, and the arithmetic mean of the v sin i measured with DFT and FAMIAS. The result shows that the primary rotates sub-synchronously, only in 2001 it reaches almost synchronous rotation velocity.

6 Orbital period changes

It is hard tosearch for orbital period changes in the range of a few seconds when the orbital RV curves are heavily distorted by non-Keplerian effects such as in the case of RZ Cas. Table 5 lists the orbital periods and times of minimum (T_Min hereafter) derived from the RVs in different seasons using the method of differential corrections (Schlesinger 1910; Sterne 1941). The 1σ errors of the period are of the order of two seconds or more, which is larger than the expected period changes. Also the inclusion of effects such as Roche geometry of the components and spots on the stellar surfaces into the PHOEBE calculations (Sect. 7.2) did not help us reach the desired accuracy in orbital period. The T_Min, on the other hand, could be very precisely determined.

When we plot the RVs versus orbital phase where the latter is calculated from T_min and P of one single season, we see a clear phase shift in RV compared to other seasons (Fig. 9). Since we do not have any information about the behaviour of the system in between, we cannot deduce period shifts from our RVs alone, however. To fix the problem, we used the photometric T_Min from literature. We collected data from the O–C Gateway¹, Bob Nelson’s Data Base of O–C Values², T_Min kindly provided by J. Kreiner³ (also see Kreiner 2004), and unpublished data from D. Mkrtichian. Cross-checking for duplicates and rejecting all visual observations, we ended up with 605 T_Min covering the time span from 1896 to 2019, to which we added the eight T_Min derived from our spectra. We converted all dates given as HJD to BJD based on terrestial time TT. Then we computed the overall best fitting period from a least-squares fit T_Min versus season number E, yielding P₀ = 1.195250392(61) d and T₀ = 2 453 775.89453(79). The resulting values (5)

are shown in Fig. 10.

There exist different approaches to determine the local period from an O–C diagram. The classic method is to fit segments by linear or parabolic functions to calculate constant periods or linearly changing periods per segment, respectively. But there are also approaches that assume a continuous change of the orbital period like that of Kalimeris et al. (1994), see Rovithis-Livaniou (2001) for an overview on other methods. We applied a similar procedure as used in Mkrtichian et al. (2018), interpolating the O–C data to a grid of step width one in season E, using spline fits together with 3σ-clipping and computing the period change from the local slope of the resulting fit. It is (6)

where T_E −T_E−1 is the local period and thus Eq. (6) describes the local difference △P = P−P₀. A necessary precondition for our method is that the observed data points are sufficiently dense so that the spline fit gives a reliable prediction of the behaviour between these points. For that reason, we only considered all T_Min observed after BJD 2 437 000.

Figure 11 shows in its top panel the O–C values together with the spline fit. The T_Min obtained from our RVs (T_RV in Table 5) are shown in green and fit very well. The bottom panel shows the resulting period changes by the solid line, where we assumed a “best-selected” smoothness parameter for the underlying spline fit. To give an impression of the influence of the smoothness parameter, we also show (dotted line) the period changes resulting from a much larger parameter, leading to a more relaxed spline fit. The positions of our spectroscopic observationsare indicated. Obtained period changes and periods are listed in Table 5 as △P_LC and P_LC.

Mkrtichian et al. (2018) derived typical timescales of 4.8, 6.1, and 9.2 yr from the O–C variations. We did a frequency search in the obtained period changes using the PERIOD04 program (Lenz & Breger 2005) and found the six periods listed in Table 6. We do not assume that the observed variation can be described as strictly cyclic but consider the found periods as typical timescales that describe the behaviour in certain seasons. Figure 12 illustrates this. It shows the calculated period changes together with the best fit of the six periods and also single contributions from four of the six periods. The longest period of 52.7 yr describes a long-term trend in the data. The 6.3 yr period describes the variations before BJD 2 442 560 (segment A), and the 8.6 yr period the behaviour between 2 445 240 and 2 456 070 (segment B). The 14.5 yr period is responsible for an amplitude variation over a longer timescale. The remaining two periods of 6.9 and 10.8 yr cannot directly be linked to the variations in this way, but improve the fit by counting for their non-cyclicity. When doing a separate frequency search in segments A and B, we find the dominant periods as 5.9 and 8.4 yr, respectively. This shows that we can only estimate the order of the time scales underlying the orbital period variations but cannot use the errors obtained from frequency search as a measure for the accuracy of the obtained values.

Fig. 9 Radial velocities of the primary in 2006 (black) and 2001 (red) vs. orbital phase. Phase is calculated from P and Tmin as derived from the RVs in 2006.

Fig. 10 O–C values calculated from the corrected TMin, full range inJD.

Fig. 11 O–C values and derived period changes. Top: O–C values (filled black circles) fitted by splines (red line). Values considered as outliers are shown by open circles and those derived from our spectra in green. Bottom: period changes calculated from the local slope of the spline fit (see text). Open circlesindicate the seasons of our spectroscopic observations.

Fig. 12 Calculated period changes (solid line) fitted by six frequencies (dashed line). The dotted lines show (from top to bottom) thecontributions of the 52.7, 6.3, 8.6, and 14.5 yr periods.

7 Radial velocity variations with orbital phase

7.1 Keplerian approach

Figure 13 shows the RVs folded with the orbital period after subtracting the best-fitting Keplerian orbits computed with the method of differential corrections as mentioned in Sect. 6. We use the observations from 2006 when RZ Cas was in a quiet state for comparison, shown by black dots. The deviations from a straight line (pure Keplerian motion) of the RVs of the secondary seen in 2006 are due to its non-spherical shape and inhomogeneous surface intensity distribution as discussed in Tkachenko et al. (2009), which we investigate in detail in Sect. 7.2. More or less strong deviations from the behaviour in 2006 can be seen in different years.

Looking at the behaviour of the RVs of the secondary, we find the strongest deviations in 2001 where we see large variations with orbital phase. Almost no differences to 2006 are observed in 2009, 2015, and 2016, whereas moderate differences occur in 2013 and 2014. Interpreting the strength of the deviations as activity indicator, we conclude that we observed RZCas in or just after a mass-transfer episode in 2001 and in a quiet state in 2006 (as already stated in Lehmann & Mkrtichian 2008, and Tkachenko et al. 2009), and that this quiet phase continued until 2009, followed by a slightly more active phase around 2013 and 2014, and falling back into a quiet state in 2015 and 2016.

Looking at the behaviour of the RVs of the primary, in particular at the amplitude and shape of the RME, we see a different picture. Theamplitude of the RME in 2001 is much larger and its shape is strongly asymmetric. In all the other years, however, we see almost no difference to 2006, except for three nights of observation covering the ingress of the eclipses: one in 2009, one in 2013, and one in 2015.

Fig. 13 Radial velocity residuals. Left: Rossiter-McLaughlin effect in the RVs of the primary in different years (red) after subtracting the best-fitting Keplerian orbits. For comparison, the RVs from 2006 are shown in black. In 2013, the RVs are divided into BJD < 2 456 600 (2013a, red) and BJD > 2 456 600 (2013b, green), and in 2015 into BJD < 2 457 250 (2015a, red) and BJD > 2 457 250 (2015b, green). Right: same for the RVs of the secondary, shown over a larger range in orbital phase. Phase zero corresponds to Min I in each case.

7.2 Analysis with PHOEBE

RZ Cas is a semi-detached Algol-type binary system and its cool component fills its critical Roche lobe. Tidal distortions occur and lead to non-spherical shapes. According to Sybilski et al. (2013), non-negligible effects on RVs occur for a < 20 (R₁ + R₂), that is when the separation of the components is smaller than 20 times the sum of their radii. Thus, the a priori assumption of Keplerian orbit, which considers the mass centre of the stellar disc coinciding with its intensity centre, is no longer valid. As already mentioned in previous sections, RZ Cas is undergoing stages of active mass transfer. A hot region around the equatorial belt has long been known (Olson 1982). Moreover, the recent comprehensive tomographic study of Richards et al. (2014) revealed clear indicators of mass stream activity in several short period Algols, including RZ Cas. Unno et al. (1994) found that the secondary of RZ Cas can be modelled only when assuming an unusually large value of the gravity darkening exponent of 0.53 and inferred that dark spots are present on the front and back sides of the secondary with respect to the primary. These authors supposed that quasi-radial flow in the sub-adiabatic stellar envelope from the deep interior is the cause of darkening. Tkachenko et al. (2008, 2009) confirmed the finding of the two spots, resuming the interpretation. Dervişoğlu et al. (2018) and Berdyugin et al. (2018), on the other hand, give an alternative explanation, showing for the short-period Algols δ Lib and λ Tau that the mass stream can produce a light scattering cloud in front of the surface of the Roche lobe filling secondary facing the primary.

7.2.1 Method

To account for all these effects, we need to use sophisticated models to simulate RV changes during the whole orbital motion via integrating both components surface intensities. For this purpose, we used the well-known Wilson-Devinney code (Wilson & Devinney 1971; Van Hamme & Wilson 2005) through the PHOEBE interface (Prša & Zwitter 2005). The WD code consists of light and/or RV curve synthesiser (LC) and parameter optimiser (DC) for fitting purpose. As we describe later, our attempts to optimise the multi-parameter fit failed with DC, which is expected for such complex configurations. Thus we used LC via our Markov chain Monte Carlo (MCMC) optimiser written in Phoebe-scripter extension, to both optimise the parameters and explore the parameter space (see Dervişoğlu et al. 2018).

In MCMC runs, we start from a random initial parameter set and add accepted parameters to the Markov chain. We need to run a lot of chains (or “walkers”) to avoid the issue of one Markov chain sticking in a local χ² valley. That is why it is advisable to start the simulation with at least two times the number of walkers of the number of parameters to be optimised (see Foreman-Mackey et al. 2013). For our last sets of simulation, we set the number of walkers >20 and the lowest iteration numbers are dynamically increased to fulfil the condition by Foreman-Mackey et al. (2013) that the iteration number should be at least ten times the auto-correlation time (i.e. typically 500 000 iterations per season).

7.2.2 Application

Spots are described in WD by two coordinates, co-latitude δ and longitude λ, and two physical parameters, temperature ratio compared to the normal surface T and radius R. In our initial run, we started to model the RVs of the cool component by taking one spot (Spot1) on the cool secondary facing the primary into account to mimic a scattering cloud or diffusive material between the components. All of our models converged to spot position δ ≈ 90°, λ ≈ 0°, that is the region that faces the primary. However, we found a strong correlation between spot size and temperature ratio of the form of R² T⁴ = constant, balancing a lower temperature ratio by a more extended spot. Thus, we could only fit the “strength” or “contrast” of the spot with respect to the stellar surface, fixing the temperature ratio to a reasonable lower limit of 0.76. As mentioned at the beginning of Sect. 7.2, Unno et al. (1994) suggested dark spots at the front and back sides of the cool secondary towards the primary caused by mass transfer. We therefore implemented a second spot (Spot2) on the opposite side of the surface of the secondary to check if this offers further improvement.

To check if we can detect a variation of the filling factor between active and inactive phases of RZ Cas, we used the “unconstrained mode” of WD and the surface potential of the secondary as a free parameter. In all of our runs, however, the filling factor of the cool component converged to unity. So for the final analysis, we fixed it to unity and used the “semi-detached mode”. Eccentricity was set to zero and the orbital inclination to 82°. The temperature ratio was fixed to 0.76 for both spots on the secondary and, as already mentioned, the position of Spot1 to δ₁ = 90°, λ₁ = 0°. The location of Spot2 converged to δ₂ = 90° without showing larger scatter, so we fixed δ₂ as well.

7.2.3 Results

Free parameters were the synchronisation parameter for the primary, F₁, the radii of Spot1 and Spot2 on the secondary, and , the longitude of Spot2, , and the systemic velocity V_γ. The RV semi-amplitudes were allowed to vary within the error bars derived in previous attempts and led to the values as listed in Table 7.

In Fig. A.1, we show the corner plot for season 2009. This shows that MCMC not only gives the most optimal parameter set, it also provides an error estimation on each parameter and possible correlations between the parameters. Finally obtained values are listed in Table 8, in which we also included the mean scatter calculated from the RV residuals after subtracting the PHOEBE solutions obtained for each season. These residuals are shown in Fig. A.2.

Figure 14 shows the influence of the inclusion of spots on the secondary into the model. We selected three seasons as examples. Including spots on the secondary did not effect the modelled RVs of the primary and vice versa. Thus, we only show the residuals of the RVs of the secondary. The black dots corresponds to the best-fitting solutions without spots, that is when only the non-spherical shapes of the components (in particular of the secondary) are taken into account. The RVs show a systematic deviation around Min II that is distinctly reduced when including Spot1 that faces the primary, resulting in the red dots. The RVs show smaller, non-systematic deviations around Min I, which is clearly present in the data from 2013a, small in 2016, and almost not visible in 2015b. Adding the second spot on the oppositeside reduces this scatter, as shown by the green dots, but in some cases the reduction is marginal.

One explanation for the asymmetry observed in the RME in the years 2009, 2013b, and 2015a (see Fig. 13) could be theimpact of a hot spot on the primary on its RVs. For comparison, we therefore additionally included a hot spot on the primary for these seasons, taking three more free parameters (temperature ratio, radius, and longitude) into account. The fit did not improve the solutions, however.

Table 9 lists absolute values of positive and negative amplitudes of the RME (maximum deviations in RV during ingress and egress, respectively) obtained after subtracting the best-fitting PHOEBE solution for F₁ = 0, that is for a non-rotating primary, from the observed RVs. Panels h and i in Fig. 15 show that the positive amplitudes reveal a similar behaviour as F₁ or R₂, whereas the negative amplitudes show a systematically decreasing trend with time (see next section for a discussion).

Figure 15 shows the seasonal variations of all investigated parameters. In panel a we added the period change derived from the O–C values (cf. Sect. 6), in panel b the v sin i values from the fit determined in Sect. 5.4, and in panel g the mean scatter of the RV residuals after subtracting the best-fitting PHOEBE solutions from the input data (see Fig. A.2). All variations can be described by a sinusoid plus a long-term trend. Because the seasonal sampling was not sufficient to determine a second frequency describing the long-term trend, we fixed it to 10⁻⁶ c d⁻¹ and determined only amplitudes and phases. Table 10 lists the timescales derived from the best-fitting sinusoids.

The best-fitting curves are shown in Fig. 15 by solid lines. We see that F₁, R₂, rms, and K_p vary almost in phase, whereas R₁ varies in anti-phase. Moreover, v sin i and K_p show almostidentical shapes of variation. Searching for a possible common period that explains the variations of all parameters, we found a period of 9 yr that best fits, together with the mentioned long-term trend, the variations of all parameters except for that can be fitted by a period of 18 yr, twice the period of 9 yr. For dP we had to add a third (optimised) period of 6.3 yr. The resulting fits are shown in Fig. 15 by dotted lines. In all cases, the quality is comparable to that obtained from the optimised values listed in Table 10.

Fig. 14 Residuals of the RVs of the secondary after subtracting the best-fitting PHOEBE solutions without including spots (black), including one spot (red), and including two spots (green) on the secondary, obtained (from top to bottom) for seasons 2013a, 2015b, and 2016.

8 Discussion

Our spectral analysis yields atmospheric parameters of the components that are in agreement with the results of light curve analysis by Rodríguez et al. (2004). For the first time, we determine the elemental surface abundances. For both components, we find [Fe/H] close to − 0.4, as for Ca, Cr, Mn, and Ni of the primary, whereas O, Mg, Si, Sc, Ti, and V show relative abundances between − 0.2 and − 0.1. For the carbon abundance of the primary, we find [C/H] = , which is remarkably different from the other elements. We think that the difference is significant. First, the considered spectral range shows a sufficiently large number of C I lines of the primary (SynthV computes line depths >5% for 32 rotationally unbroadened C I lines). Second, differences between the fits with [C/H] of − 0.4 and − 0.8 can clearly be seen by eye. Figure 16 shows this for the strongest carbon lines in the Hβ region of the spectrum of the primary (observed composite spectrum after subtracting the synthetic spectrum of the secondary). In the spectrum of the cool secondary, carbon is mainly present in form of the CH molecule bands. Compared to the primary, the signal is very weak and we can only say that [C/H] is below solar abundance, between − 0.3 and − 1.0 within the 1σ error bars.

Mennekens et al. (2008) tried to model the binary evolution of RZ Cas and found consistent solutions only for an initial mass ratio of q ≈ 3, which is about the inverse of the actual ratio. From the large initial mass of the donor in that case we can assume that the primary has switched during its evolution from pp-chain to CNO-cycle hydrogen burning, resulting in depleted carbon and enhanced nitrogen abundances in its core (e.g. Przybilla et al. 2010). From the reversal of mass ratio we conclude that the donor was stripped by mass loss down to its core so that the gainer accreted CNO-cycled, carbon-deficient material in the late phase of its evolution. This material has then mixed with the surface layers of the gainer, leading to the observed carbon abundance. In that case, the surface carbon abundance of the gainer cannot be lower than that of the donor star, however. Multiple authors have used the sketched scenario to explain the surface abundance anomalies observed in several other Algol-type stars; these include Ibanoǧlu et al. (2012) for GT Cep, AU Mon, and TU Mon, Kolbas et al. (2014) for u Her, and Dervişoğlu et al. (2018) for δ Lib.

Our further investigation was based on the separated LSD profiles and RVs of the two components of RZ Cas computed with LSDbinary, and on the times of minimum (O–C values) taken from literature. The change in orbital period computed from the slope of the O–C diagram cannot be characterised by one single timescale, as already found by Mkrtichian et al. (2018). These authors derived timescales of 4.8, 6.1, and 9. 2 yr. We observe from our slightly extended data set different timescales of variation in different seasons, such as 6.3 yr around JD 2 440 000, and 8.6 yr from 2 450 000 to 2 455 000 (which is about the time span of our spectroscopic observations), overlaid by longer periods. From our previous analysis using the Shellspec07_inverse code (Tkachenko et al. 2009), we know that RZ Cas was in an active state of mass transfer in 2001 and in a quiet state in 2006. Comparison of the RV residuals after subtracting a pure Keplerian orbit with those from 2006 gave us initial hints pointing to further activity periods (Fig. 13). From the RVs of the primary we see that the RME was distinctly enhanced only in 2001. A comparison of the RV residuals of the secondary with those from 2006 shows the largest deviation in 2001, weaker deviations in 2013 and 2014, and almost no deviations in 2009, 2015, and 2016. The amount of scatter found in the RV residuals after subtracting the best-fitting PHOEBE solutions is strongly correlated with this finding.Thus, we conclude that no further mass-transfer episode as strong as in or close to 2001 occurred in RZ Cas.

The modelling of the observed RVs of both components with MCMC-PHOEBE gave the best results when adding two dark spots on thesurface of the cool companion: Spot1 facing the primary, Spot2 on the opposite side. Spot1 was already found from the RVs from 2001 and 2006 by Tkachenko et al. (2009), whereas Unno et al. (1994) included two spots into their model, explaining the observed anomalous gravity darkening by a cooling mechanism by enthalpy transport due to mass outflow that leads to a reverse process of gravitational contraction.

We monitored both spots for the first time over decades. We find that Spot1 always exactly points towards the primary, with radii (a synonym for strength or contrast as mentioned in before) between 23° and 43°. Spot2, on the other hand, is much weaker, shows a variation in its position, and induces only a second order improvement in some of the observed seasons (cf. Fig. 14). The main findings are that the strengths of the two spots vary in anti-phase and that Spot2 shows different positions in longitude, varying around the longitude of the L2 point of 180°.

The fact that the strength of Spot1 is largest when the star is between 2006 and 2009 in a quiet phase (Fig. 15) speaks against the cooling by mass-outflow mechanism. Instead, we argue that the variations of the strengths of both spots can be explained by magnetic activity of the cool companion, assuming an activity cycle of 9 yr, based on an 18-yr cycle of magnetic field change, including a reversal of the magnetic poles. The 9-yr cycle can be found in the variations of all investigated parameters except for the longitude of Spot2, as we showed in the last section, and was also found by Mkrtichian et al. (2018). Spot2 shows a migration in longitude, returning after an 18-yr cycle to its position from 2001 (cf. Fig. 15). We assume that we observe similar surface structures on the cool secondary of RZ Cas such as for cool, rapidly rotating RS CVn binaries or single rapidly rotating variables of FK Com and BY Dra type. Long-living active regions were observed on opposite sides of these stars. The spots are of different intensity and also show switching of activity from one spot to the other on timescales of years or decades, which is known as the flip-flop effect (see Yakobchuk & Berdyugina 2018, and references herein).

We believe that there is a direct analogy with sunspot activity. Sunspots show a saucer-shaped depression in the photosphere caused by the Lorentz force of the strong magnetic field within the spot, the so-called Wilson depression (Wilson & Maskelyne 1774). This means that inside the spot the level of τ = 1 is located below the level of the photosphere outside the spot. For the Sun, the geometric depth of the depression is of the order of 600 km (e.g. Gokhale & Zwaan 1972; Löptien et al. 2018). For the Roche-lobe filling donor, the existence of such atmospheric depression close to the L1 point means that L1 is “fed” by atmospheric layers of lower density and thus the mass transfer is the more suppressed the deeper the depression is. The local magnetic field strength controls the strength (depression) of Spot1 and the height of atmospheric layers feeding the L1 point and in this way the mass-transfer rate. RZ Cas showed, in perfect agreement with the drafted scenario, high mass-transfer rate in 2001 (and possibly in 2011, when observations are missing) when the Spot1 size was around 20 degrees and low activity state in 2006-2007 and in 2015-2016 when the spot size was about 40 degrees.

The explanation for the asymmetry and different amplitudes observed in the RME in the years 2001, 2009, 2013b, and 2015a (see Fig. 13) could be the impact of the combined effect of acceleration of the photospheric layers in 2001 caused by the high mass-transfer rate and by screening the surface of the primary by the dense gas stream (Lehmann & Mkrtichian 2008). The amplitudes of the RME in RZ Cas in general show very different behaviour during ingress (K_p) and egress (K_n) of primary eclipse. The variation of K_p is similar to those of R₂ and rms, whereas the shape of the K_n variation perfectly matches that of v sin i (see Fig. 15). We assume that K_n is related to the true v sin i seen outside the eclipses and K_p is strongly influenced by mass-transfer effects. This can be explained from Fig. 17, showing the gas density distribution around RZ Cas for two different mass-transfer rates, calculated from 3D-hydrodynamic simulations by Nazarenko & Mkrtichian (priv. comm.). It can be seen that the equatorial zone of the surface of the primary is masked during ingress (at orbital phases ϕ = − 0.1..0) by the optical thick gas stream from the secondary, whereas we directly see the complete disc of the primary during egress (ϕ = 0..0.1), only very slightly hampered by optical thin circumbinary matter. In consequence, K_p is affected by the seasonal varying density of the gas stream (or mass-transfer rate) and K_n is correlated with the surface rotation velocity of the primary. All the variability seen in RME in the years from 2006 on is mainly related to the variations ofK_p, forced by changes in gas-stream density and variable screening and attenuation effects, while the rotation speed was nearly constant.

The synchronisation factor derived with PHOEBE is mainly based on the shape of the RME and this shape is strongly influenced by timely varying Algol-typical effects (see paragraph below). Our MCMC simulation did not count for the observed asymmetry in the RME either (i.e. the fact that the amplitudes K_p and K_n are different from each other) and so the synchronisation factor F₁ has to be considered as some kind of mean value. The shape of its time variation resembles the mean of the shapes of the K_p and K_n variations (cf. Fig. 15). The results obtained for the synchronisation factor in Sect. 5.4 from radius and v sin i of the primary,on the other hand, give sub-synchronous rotation of the primary. Its outer layers are only accelerated to almost synchronousrotation during the active phase in 2001. The finding of sub-synchronous rotation is surprising. It was suggested by Davis et al. (2014) to occur during the rapid phase of mass exchange in Algols, when the donor star is spun down on a timescale shorter than the tidal synchronisation timescale and material leaving the inner-Lagrangian point is accreted back onto the donor, enhancing orbital shrinkage. But these authors state that once the mass-transfer rate decelerates and convection develops in the surface layers, tides should be effective enough to re-synchronise the primary. To our knowledge, sub-synchronous rotation was found so far in only one short-period Algol, namely TV Cas, by Khalesseh & Hill (1992), in which the authors found synchronisation factors of 0.85 for the primary and 0.65 for the secondary component.

The RV residuals after subtracting the solutions obtained with our MCMC simulations (Fig. A.2) finally show that our model distinctly reduces the scatter compared to the residuals obtained from pure Keplerian motion. On the other hand, we can still recognise many of the signs of activity that we discussed in connection with Fig. 13. We assume that these still unexplained features are due to the distribution and density of circumbinary matter along the line of sight in different orbital phases, varying between different seasons according to the varying activity of the star. Finally, we can conclude from Fig. A.2 in the same way as from Fig. 13 that RZ Cas showed an extraordinary phase of activity around 2001 followed by a quiet phase in 2006 to 2009, slightly enhanced activity in 2013/2014, followed by a quiet phase again. The calculated fits on the variation of the radii of the two spots and the RME amplitude K_p based on the nine-year cycle, on the other hand, show extrema of the same amplitude such as in 2001 for the time around 2010–2011. Therefore, because of missing data in this period, we cannot exclude the possibility that a second mass-transfer episode of comparable strength like in 2001 occurred close to 2010–2011.

The orbital period was at maximum shortly after 2001, dropped down steeply until 2006, and was then more or less continuously rising. From Eq. (8) in Biermann & Hall (1973), assuming conservative mass transfer and conservation of orbital angular momentum, we obtain a mass-transfer rate of 1.5 × 10⁻⁶ M_⊙ yr⁻¹. This is a typical value observed for Algol-type stars (e.g. Hall & Neff 1976) and also agrees Hall et al. (1976), who calculate, from an O–C analysis, a mean mass-transfer rate of RZ Cas of 1.0 × 10⁻⁶ M_⊙ yr⁻¹, averaged over several episodes of period change. As mentioned in the introduction, evolutionary scenarios for RZ Cas have to count for mass loss from the system (Mennekens et al. 2008), which may be the case for Algols in earlier stages of their evolution in general (e.g. Ibanoǧlu et al. 2006; Deschamps et al. 2013, 2015). Thus, conservatism is not necessarily justified and the assumption of mass conservation can only lead to a raw estimation of the amount of mass transferred during the active phase of RZ Cas in 2001.

We assume that the enhanced value of v sin i in 2001 points to an acceleration of the outer layers of the primary of RZ Cas due to the mass-transfer episode occurring close to that year. This then dropped down by about 5 km s⁻¹ and it could be that its slight increase after 2009 is correlated with a slightly increased activity as indicated by the variation of F₁ and rms.

Fig. 15 Time variations of different parameters (see text). The solid lines are calculated from the best-fitting periods and long-term trends. The dotted lines show possible sinusoidal variations with a period of 9 yr plus a long-term trend, except for panel f, which shows an 18 yr variation. Values from 2013b and 2015a were considered as outliers.

Fig. 16 Fit of the observed spectrum of the primary (black) in the Hβ region by synthetic spectra with [C/H] of 0.0 (blue), − 0.4 (green), and − 0.8 (red).

Fig. 17 Logarithmic plots of gas density distributions (in 1011 cm−3) obtained from 3D hydrodynamic simulations of the RZ Cas system (Nazarenko & Mkrtichian, priv. comm.) for mass-transfer rates of 1 × 10−9 M⊙ y−1 (left) and 6 × 10−8 M⊙ y−1 (right). The viewing angles at first (ϕ = − 0.1) and last (ϕ = + 0.1) contact are indicated.

9 Conclusions

In this first of three articles related to spectroscopic long-term monitoring of RZ Cas, we investigated high-resolution spectra with respect to stellar and system parameters based on the RVs and LSD profiles of its components calculated with the LSDbinary program. The main goal was to search for further episodes of enhanced mass transfer occurring after 2001 and for a general timescale of variations possibly caused by the magnetic cycle of the cool companion.

From spectrum analysis we determined precise atmospheric parameters, among them low metal surface abundances, in particular [Fe/H] of − 0.42 and [C/H] of − 0.80 for the primary and [Fe/H] of the same order for the secondary. The carbon deficiency observed for the primary gives evidence that the outer layers of the cool secondary have been stripped in the fast mass-transfer phase down to its core so that CNO-cycledmaterial was transferred to the outer layers of the primary in later stages of evolution. The derived T_eff of the components and the logg of the primary agree withinthe 1σ error bars with the results from LC analysis by Rodríguez et al. (2004). From the RV analysis with PHOEBE, we derived very precise masses and separation of the components of M₁ = 1.951(5) M_⊙, M₂ = 0.684(1) M_⊙, and a = 6.546(5) R_⊙.

From several of the investigated parameters that show seasonal variations, such as orbital period, v sin i, strength of the spots on the secondary, synchronisation factor calculated with PHOEBE, and rms of the RV residuals after subtracting the PHOEBE solutions, we deduce a common time scale of the order of 9 yr. The variation of the orbital period is complex and can be described in detail only when adding further periods. We conclude that we see the effects of a 9-yr magnetic activity cycle of the cool companion of RZ Cas, caused by an 18-yr dynamo cycle that includes a reversal of polarity. This conclusion is strongly supported by the behaviour of the two dark spots on the surface of the secondary that show the flip-flop effect in their strengths and one spot that shows an 18-yr periodicity in longitudinal migration.

From the variations of orbital period and v sin i around 2001, we conclude that the determined v sin i is the projected equatorial rotation velocity of the outer layers of the primary accelerated by transferred matter and does not stand for the rotation velocity of the star as a whole. In all other seasons, the measured v sin i point to a sub-synchronous rotation of the gainer. At the present stage, we cannot give a physical explanation for this new and interesting finding, however.

Based on our available data, we conclude that RZ Cas was undergoing an episode of high mass transfer in 2001, in a quiet phase in 2006 and 2009, followed by a slightly more active phase in 2013 and 2014, and again in a quiet phase in 2015 and 2016. Because we did not observe the star in 2010 and 2011, we cannot exclude that a second episode of high mass transfer occurred in these years, which would agree with the derived magnetic activity cycle of nine years.

The results of our investigation of high-frequency oscillations of the primary of RZ Cas will be presented in Part II of this article. In a third article, we will investigate the accretion-induced variability of He I lines detected in the spectra.

Acknowledgements

H.L. and F.P. acknowledge support from DFG grant LE 1102/3-1. V.T. acknowledges support from RFBR grant No. 15-52- 12371. A.D. is financially supported by the Croatian Science Foundation through grant IP 2014-09-8656 and Erciyes University Scientific Research Projects Coordination Unit under grant number MAP-2020-9749. The research leading to these results has (partially) received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement N°670519: MAMSIE), from the KU Leuven Research Council (grant C16/18/005: PARADISE), from the Research Foundation Flanders (FWO) under grant agreement G0H5416N (ERC Runner Up Project), as well as from the BELgian federal Science Policy Office (BELSPO) through PRODEX grant PLATO. Results are partly based on observations obtained with the HERMES spectrograph, which is supported by the Research Foundation – Flanders (FWO), Belgium, the Research Council of KU Leuven, Belgium, the Fonds National de la Recherche Scientifique (F.R.S.-FNRS), Belgium, the Royal Observatory of Belgium, the Observatoire de Genève, Switzerland and the Thüringer Landessternwarte Tautenburg, Germany.

Appendix A Analysis of radial velocity variations with PHOEBE

Figure A.1 shows an example of a corner plot obtained with the PHOEBE-based MCMC method. Figure A.2 shows the RV residuals after subtracting the PHOEBE solutions (cf. Sect. 7.2.3).

Fig. A.1 Corner plot for the year 2009, showing the χ2 distributions for systemic velocity, synchronisation factor of the primary, radius of the spot on the secondary facing the primary, and longitude and radius of the opposite spot on the secondary.

Fig. A.2 Residuals after subtracting the PHOEBE solutions from the RVs of the primary (left) and the secondary (right) component of RZ Cas.

References

All Tables

All Figures

Fig. 1 Continuum flux ratios between the components vs. wavelength. The dotted lines are calculated from synthetic continuum spectra with Teff1 of 8700 K and Teff2 of 3800–5000 K. The best agreement with the flux ratio from photometry (red line) is obtained for Teff2 = 4400 K (black solid line). The upper green line represents Teff1 = 9000 K, Teff2 = 4500 K, the lower green line for Teff1 = 8400 K, Teff2 = 4300 K.

In the text

Fig. 2 Results of spectrum analysis using method a (cf. Sect. 3.3, Table 2), shown for the Hδ -Hγ region (left) and the region around the Mg I b triplet (right). First row: observed, continuum adjusted composite spectrum (black), best-fitting synthetic spectrum (red), and shifted difference spectrum (green). The inverse of the applied continuum correction is shown in blue. Second row: same for the secondary. Third row:flux ratio of the secondary to primary from photometry (black) from spectrum analysis based on R2 ∕R1 = 0.926 (red) and assuming the adjusted ratio of 0.78 (green).

In the text

	Fig. 3 As Fig. 2, but using method b. Third row: flux ratio from photometry (black) and from spectrum analysis based on R2∕R1=1.17 (red) and on the calculated ratio of 1.065 (green).
In the text

	Fig. 4 AsFig. 2, but using method c. Third row: flux ratio from photometry (black) from spectrum analysis based on the photometric flux ratio assuming R2 ∕R1 = 1.17 (red) and on the adjusted ratio of R2∕R1 = 0.87 (green). Bottom row: ratio of the photometric flux ratio to the flux ratio obtained from spectrum analysis based on R2 ∕R1 = 1.17.
In the text

	Fig. 5 LSD profiles computed with LSDbinary from spectra taken during primary eclipse for the primary (left) and secondary (right) of RZ Cas. Phase zero corresponds to Min I.
In the text

	Fig. 6 Values of v sin i obtained from DFT vs. out-of-eclipse phases measured from the LSD profiles observed in 2016 (black) and from the same profiles corrected for the LPV found with FAMIAS (red).
In the text

	Fig. 7 Values of v sin i determined from the FeII 5317 Å line vs. orbital phase, shown for 2001 and 2014. The mean values were built from the values indicated in green.
In the text

	Fig. 8 Values of v sin i measured from (a) DFT (black crosses), (b) the Fe II line (red squares), and (c) FAMIAS (green plus signs). The solid and dotted lines show possible sinusoidal variations (see text).
In the text

	Fig. 9 Radial velocities of the primary in 2006 (black) and 2001 (red) vs. orbital phase. Phase is calculated from P and Tmin as derived from the RVs in 2006.
In the text

	Fig. 10 O–C values calculated from the corrected TMin, full range inJD.
In the text

	Fig. 11 O–C values and derived period changes. Top: O–C values (filled black circles) fitted by splines (red line). Values considered as outliers are shown by open circles and those derived from our spectra in green. Bottom: period changes calculated from the local slope of the spline fit (see text). Open circlesindicate the seasons of our spectroscopic observations.
In the text

	Fig. 12 Calculated period changes (solid line) fitted by six frequencies (dashed line). The dotted lines show (from top to bottom) thecontributions of the 52.7, 6.3, 8.6, and 14.5 yr periods.
In the text

Fig. 13 Radial velocity residuals. Left: Rossiter-McLaughlin effect in the RVs of the primary in different years (red) after subtracting the best-fitting Keplerian orbits. For comparison, the RVs from 2006 are shown in black. In 2013, the RVs are divided into BJD < 2 456 600 (2013a, red) and BJD > 2 456 600 (2013b, green), and in 2015 into BJD < 2 457 250 (2015a, red) and BJD > 2 457 250 (2015b, green). Right: same for the RVs of the secondary, shown over a larger range in orbital phase. Phase zero corresponds to Min I in each case.

In the text

	Fig. 14 Residuals of the RVs of the secondary after subtracting the best-fitting PHOEBE solutions without including spots (black), including one spot (red), and including two spots (green) on the secondary, obtained (from top to bottom) for seasons 2013a, 2015b, and 2016.
In the text

	Fig. 15 Time variations of different parameters (see text). The solid lines are calculated from the best-fitting periods and long-term trends. The dotted lines show possible sinusoidal variations with a period of 9 yr plus a long-term trend, except for panel f, which shows an 18 yr variation. Values from 2013b and 2015a were considered as outliers.
In the text

	Fig. 16 Fit of the observed spectrum of the primary (black) in the Hβ region by synthetic spectra with [C/H] of 0.0 (blue), − 0.4 (green), and − 0.8 (red).
In the text

	Fig. 17 Logarithmic plots of gas density distributions (in 1011 cm−3) obtained from 3D hydrodynamic simulations of the RZ Cas system (Nazarenko & Mkrtichian, priv. comm.) for mass-transfer rates of 1 × 10−9 M⊙ y−1 (left) and 6 × 10−8 M⊙ y−1 (right). The viewing angles at first (ϕ = − 0.1) and last (ϕ = + 0.1) contact are indicated.
In the text

	Fig. A.1 Corner plot for the year 2009, showing the χ2 distributions for systemic velocity, synchronisation factor of the primary, radius of the spot on the secondary facing the primary, and longitude and radius of the opposite spot on the secondary.
In the text

	Fig. A.2 Residuals after subtracting the PHOEBE solutions from the RVs of the primary (left) and the secondary (right) component of RZ Cas.
In the text