Issue 
A&A
Volume 644, December 2020



Article Number  A121  
Number of page(s)  18  
Section  Stellar atmospheres  
DOI  https://doi.org/10.1051/00046361/202039355  
Published online  10 December 2020 
Spectroscopic longterm monitoring of RZ Cas
I. Basic stellar and system parameters^{★}
^{1}
Thüringer Landesternwarte Tautenburg,
Sternwarte 5,
07778
Tautenburg, Germany
email: lehm@tlstautenburg.de
^{2}
Department of Astronomy and Space Sciences, Erciyes University,
38039
Kayseri, Turkey
^{3}
Astronomy and Space Sciences Observatory and Research Center, Erciyes University,
38039
Kayseri, Turkey
^{4}
National Astronomical Research Institute of Thailand,
260 Moo 4, T. Donkaew, A. Maerim,
Chiangmai
50180, Thailand
^{5}
Institute of Astronomy,
KU Leuven, Celestijnenlaan 200D,
3001
Leuven, Belgium
^{6}
Institute of Astronomy, Russian Academy of Sciences,
119017,
Pyatnitskaya Str. 48,
Moscow,
Russia
Received:
7
September
2020
Accepted:
4
November
2020
Context. RZ Cas is a shortperiod Algoltype system showing episodes of mass transfer and δ Sctlike oscillations of its massgaining primary component. This system exhibits temporal changes in orbital period, v sin i, and the oscillation pattern of the primary component.
Aims. We analyse highresolution spectra of RZ Cas that we obtained during a spectroscopic longterm monitoring lasting from 2001 to 2017. In this first part we investigate the atmospheric parameters of the stellar components and the time variation of orbital period, v sin i, and radial velocities (RVs), searching for seasonal changes that could be related to episodes of mass exchange and to a possible activity cycle of the system triggered by the magnetic cycle of the cool companion.
Methods. We used spectrum synthesis to analyse the spectra of both components of RZ Cas. The study of variations of the orbital period is based on published times of primary minima. We used the leastsquares deconvolved (LSD) binary program to derive separated RVs and LSD profiles of the components. From the LSD profiles of the primary we determined its v sin i. Using Markov chain Monte Carlo simulations with the PHOEBE program, we modelled the RV variations of both components.
Results. Spectrum analysis resulted in precise atmospheric parameters of both components, in particular in surface abundances below solar values. We find that the variation of orbital period is semiregular and derive different characteristic timescales for different epochs of observation. We show that the RV variations with orbital phase can be modelled when including two cool spots on the surface of the secondary component. The modelling leads to very precise masses and separation of the components. The seasonal variation of several parameters, such as v sin i, rotationorbit synchronisation factor, strength of the spots on the cool companion, and orbital period, can be characterised by a common timescale of the order of nine years.
Conclusions. We interpret the timescale of nine years as the magnetic activity cycle of the cool companion. In particular the behaviour of the dark spots on the cool companion leads us to the interpretation that this timescale is based on an 18yr magnetic dynamo cycle. We conclude that the masstransfer rate is controlled by the variable depth of the Wilson depression in the magnetic spot around the Lagrangian point L1. In the result, based on available data, we observe a damped activity cycle of the star, starting with a high masstransfer episode around 2001 with a calculated masstransfer rate of 1.510^{−6} M_{⊙} yr^{−1}, followed by quiet periods in 2006 and 2009, slightly higher activity around 2013 and 2014, and again followed by quiet periods in 2015 and 2016. However, owing to missing data for years 2010 and 2011, we cannot exclude that a second high masstransfer episode occurred within this time span.
Key words: binaries: eclipsing / binaries: spectroscopic / binaries: close / stars: atmospheres
© 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 photometricspectroscopic 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 highresolution 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) modelindependent 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).
Algoltype systems (Algols, hereafter) are semidetached, interacting EBs consisting of a mainsequence star of spectral type BA (primary, hereafter) and an evolved FK type companion (secondary, hereafter). As a consequence of the abovementioned 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 massaccreting primary shows δ Sctlike oscillations. These stars are extraordinary objects for asteroseismic studies because they allow us to investigate shortterm dynamical stellar evolution during masstransfer 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 masstransfer episodes, the dynamical behaviour of the system, and the excitation of different nonradial pulsation (NRP) modes of the massgaining 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 socalled 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 masstransfer 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 nonaxisymmetric 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 socalled 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 leastsquares deconvolution (LSD) technique (Donati et al. 1997) and the pixelbypixel method of the FAMIAS program (Zima 2008), highldegree NRP modes can also be detected and identified from the signaltonoise 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 shortperiod (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 shortperiod light variability with a dominant oscillation mode of a frequency of 64.2 c d^{−1}. This finding was later confirmed from dedicated multisite 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 RossiterMcLaughlin 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 masstransfer effects.
After observing RZ Cas in 2001 and 2006, we took new time series of highresolution 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 lightcurve 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, highS/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 nonKeplerian 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 highfrequency 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.
Journal of observations listing the instrument, its spectral resolving power, year of observation, and mean Julian date.
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_{2} 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 T_{eff1} of 8700 K and T_{eff2} of 3800–5000 K. The best agreement with the flux ratio from photometry (red line) is obtained for T_{eff2} = 4400 K (black solid line). The upper green line represents T_{eff1} = 9000 K, T_{eff2} = 4500 K, the lower green line for T_{eff1} = 8400 K, T_{eff2} = 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_{2} = 3.7 (Rodríguez et al. 2004), and log g_{1} = 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_{2} ∕R_{1} = 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_{1} and s_{2} are the synthetic spectra of the components shifted for their RVs v_{1}, v_{2} (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_{2} and C_{1} 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 bestfitting radii ratio can then be determined from Eq. (3) using a simple leastsquares 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 χ^{2} 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.
Atmospheric parameters of primary and secondary of RZ Cas derived with multiple methods.
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 leastsquares 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 bestfitting synthetic composite spectrum. The second row compares the “observed” spectrum of the secondary with the bestfitting synthetic spectrum found for the secondary. The observed spectrum of the secondary was therefore computed by subtracting the bestfitting 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 OC 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), bestfitting 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 R_{2} ∕R_{1} = 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 R_{2}∕R_{1}=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 R_{2} ∕R_{1} = 1.17 (red) and on the adjusted ratio of R_{2}∕R_{1} = 0.87 (green). Bottom row: ratio of the photometric flux ratio to the flux ratio obtained from spectrum analysis based on R_{2} ∕R_{1} = 1.17. 
Elemental abundances.
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 shortperiod 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 outofeclipse 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 rotationbroadening 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 Algoltypical effects such as an inhomogeneous circumprimary gasdensity 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 pixelbypixel 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 pixelbypixel to build the undistorted profiles. We performed that for all seasons but 2008 for which we did not have enough data to apply the pixelbypixel method. The subtraction of the highamplitude, lowfrequency contributions found with FAMIAS (most of them are harmonics of the orbital frequency) are responsible for the cleaning, not the faint highfrequency 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 bestfitting 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 leastsquares 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 wellselected 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 lowl degree modes. For that, we cleaned the observed LSD profiles for all lowfrequency 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 χ^{2}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. outofeclipse phases measured from the LSD profiles observed in 2016 (black) and from the same profiles corrected for the LPV found with FAMIAS (red). 
Values of v sin i obtained from the three methods in multiple years.
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 lowfrequency contributions that we found with the pixelbypixel method of FAMIAS from the LSD profiles so that we removed theline broadening effects due to orbitalphase 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 longterm 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 rotationtoorbit synchronisation factor of the primary, F_{1}, 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 subsynchronously, 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 nonKeplerian 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^{1}, Bob Nelson’s Data Base of O–C Values^{2}, T_{Min} kindly provided by J. Kreiner^{3} (also see Kreiner 2004), and unpublished data from D. Mkrtichian. Crosschecking 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 leastsquares fit T_{Min} versus season number E, yielding P_{0} = 1.195250392(61) d and T_{0} = 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 RovithisLivaniou (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_{0}. 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 “bestselected” 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 longterm 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 noncyclicity. 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.
Orbital periods P_{RV} and times of minimum T_{RV} derived from the RVs from single seasons, and period changes △P_{phot} and periods P_{phot} derived from the photometric T_{Min}.
Fig. 9 Radial velocities of the primary in 2006 (black) and 2001 (red) vs. orbital phase. Phase is calculated from P and T_{min} as derived from the RVs in 2006. 
Fig. 10 O–C values calculated from the corrected T_{Min}, 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. 
Timescales and amplitudes of orbital period variation.
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 bestfitting 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 nonspherical 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 masstransfer 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: RossiterMcLaughlin effect in the RVs of the primary in different years (red) after subtracting the bestfitting 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 semidetached Algoltype binary system and its cool component fills its critical Roche lobe. Tidal distortions occur and lead to nonspherical shapes. According to Sybilski et al. (2013), nonnegligible effects on RVs occur for a < 20 (R_{1} + R_{2}), 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 quasiradial flow in the subadiabatic 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 shortperiod 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 wellknown WilsonDevinney 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 multiparameter 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 Phoebescripter 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 χ^{2} 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 ForemanMackey 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 ForemanMackey et al. (2013) that the iteration number should be at least ten times the autocorrelation time (i.e. typically 500 000 iterations per season).
Basic absolute values derived with PHOEBE.
7.2.2 Application
Spots are described in WD by two coordinates, colatitude δ 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^{2} T^{4} = 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 “semidetached 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 δ_{1} = 90°, λ_{1} = 0°. The location of Spot2 converged to δ_{2} = 90° without showing larger scatter, so we fixed δ_{2} as well.
7.2.3 Results
Free parameters were the synchronisation parameter for the primary, F_{1}, the radii of Spot1 and Spot2 on the secondary, and , the longitude of Spot2, , and the systemic velocity V_{γ}. The RV semiamplitudes 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 bestfitting solutions without spots, that is when only the nonspherical 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, nonsystematic 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 bestfitting PHOEBE solution for F_{1} = 0, that is for a nonrotating primary, from the observed RVs. Panels h and i in Fig. 15 show that the positive amplitudes reveal a similar behaviour as F_{1} or R_{2}, 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 bestfitting PHOEBE solutions from the input data (see Fig. A.2). All variations can be described by a sinusoid plus a longterm trend. Because the seasonal sampling was not sufficient to determine a second frequency describing the longterm trend, we fixed it to 10^{−6} c d^{−1} and determined only amplitudes and phases. Table 10 lists the timescales derived from the bestfitting sinusoids.
The bestfitting curves are shown in Fig. 15 by solid lines. We see that F_{1}, R_{2}, rms, and K_{p} vary almost in phase, whereas R_{1} varies in antiphase. 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 longterm 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.
Synchronisation factor of the primary, radii of both spots on the secondary, longitude of Spot2 on the secondary, and mean scatter of the RV residuals.
Fig. 14 Residuals of the RVs of the secondary after subtracting the bestfitting 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. 
Absolute values of positive and negative RV deviations due to RME.
Timescales in years determined from the seasonal variations of different parameters.
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 ppchain to CNOcycle 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 CNOcycled, carbondeficient 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 Algoltype 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 bestfitting PHOEBE solutions is strongly correlated with this finding.Thus, we conclude that no further masstransfer episode as strong as in or close to 2001 occurred in RZ Cas.
The modelling of the observed RVs of both components with MCMCPHOEBE 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 antiphase 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 massoutflow 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 18yr cycle of magnetic field change, including a reversal of the magnetic poles. The 9yr 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 18yr 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. Longliving 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 flipflop effect (see Yakobchuk & Berdyugina 2018, and references herein).
We believe that there is a direct analogy with sunspot activity. Sunspots show a saucershaped depression in the photosphere caused by the Lorentz force of the strong magnetic field within the spot, the socalled 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 Rochelobe 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 masstransfer rate. RZ Cas showed, in perfect agreement with the drafted scenario, high masstransfer rate in 2001 (and possibly in 2011, when observations are missing) when the Spot1 size was around 20 degrees and low activity state in 20062007 and in 20152016 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 masstransfer 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_{2} 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 masstransfer effects. This can be explained from Fig. 17, showing the gas density distribution around RZ Cas for two different masstransfer rates, calculated from 3Dhydrodynamic 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 masstransfer 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 gasstream 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 Algoltypical 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_{1} 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 subsynchronous rotation of the primary. Its outer layers are only accelerated to almost synchronousrotation during the active phase in 2001. The finding of subsynchronous 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 innerLagrangian point is accreted back onto the donor, enhancing orbital shrinkage. But these authors state that once the masstransfer rate decelerates and convection develops in the surface layers, tides should be effective enough to resynchronise the primary. To our knowledge, subsynchronous rotation was found so far in only one shortperiod 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 nineyear 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 masstransfer 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 masstransfer rate of 1.5 × 10^{−6} M_{⊙} yr^{−1}. This is a typical value observed for Algoltype stars (e.g. Hall & Neff 1976) and also agrees Hall et al. (1976), who calculate, from an O–C analysis, a mean masstransfer rate of RZ Cas of 1.0 × 10^{−6} M_{⊙} yr^{−1}, 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 masstransfer episode occurring close to that year. This then dropped down by about 5 km s^{−1} and it could be that its slight increase after 2009 is correlated with a slightly increased activity as indicated by the variation of F_{1} and rms.
Fig. 15 Time variations of different parameters (see text). The solid lines are calculated from the bestfitting periods and longterm trends. The dotted lines show possible sinusoidal variations with a period of 9 yr plus a longterm 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 10^{11} cm^{−3}) obtained from 3D hydrodynamic simulations of the RZ Cas system (Nazarenko & Mkrtichian, priv. comm.) for masstransfer 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 longterm monitoring of RZ Cas, we investigated highresolution 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 masstransfer phase down to its core so that CNOcycledmaterial 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} = 1.951(5) M_{⊙}, M_{2} = 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 9yr magnetic activity cycle of the cool companion of RZ Cas, caused by an 18yr 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 flipflop effect in their strengths and one spot that shows an 18yr 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 subsynchronous 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 highfrequency oscillations of the primary of RZ Cas will be presented in Part II of this article. In a third article, we will investigate the accretioninduced variability of He I lines detected in the spectra.
Acknowledgements
H.L. and F.P. acknowledge support from DFG grant LE 1102/31. V.T. acknowledges support from RFBR grant No. 1552 12371. A.D. is financially supported by the Croatian Science Foundation through grant IP 2014098656 and Erciyes University Scientific Research Projects Coordination Unit under grant number MAP20209749. 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 PHOEBEbased 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
 Aerts, C., de Pauw, M., & Waelkens, C. 1992, A&A, 266, 294 [NASA ADS] [Google Scholar]
 Aerts, C., Waelkens, C., DaszyńskaDaszkiewicz, J., et al. 2004, A&A, 415, 241 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Berdyugin, A., Piirola, V., Sakanoi, T., Kagitani, M., & Yoneda, M. 2018, A&A, 611, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Biermann, P., & Hall, D. S. 1973, A&A, 27, 249 [Google Scholar]
 Bíró, I. B., & Nuspl, J. 2011, MNRAS, 416, 1601 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Borkovits, T., Derekas, A., Fuller, J., et al. 2014, MNRAS, 443, 3068 [NASA ADS] [CrossRef] [Google Scholar]
 Budding, E., & Butland, R. 2011, MNRAS, 418, 1764 [NASA ADS] [CrossRef] [Google Scholar]
 Carroll, J. A. 1933, MNRAS, 93, 478 [NASA ADS] [CrossRef] [Google Scholar]
 Davis, P. J., Siess, L., & Deschamps, R. 2014, A&A, 570, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dervişoğlu, A., Pavlovski, K., Lehmann, H., Southworth, J., & Bewsher, D. 2018, MNRAS, 481, 5660 [CrossRef] [Google Scholar]
 Deschamps, R., Siess, L., Davis, P. J., & Jorissen, A. 2013, A&A, 557, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Deschamps, R., Braun, K., Jorissen, A., et al. 2015, A&A, 577, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Donati, J.F., Semel, M., Carter, B. D., Rees, D. E., & Collier Cameron, A. 1997, MNRAS, 291, 658 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
 Frandsen, S., Lehmann, H., Hekker, S., et al. 2013, A&A, 556, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gamarova, A. Y., Mkrtichian, D. E., Rodriguez, E., Costa, V., & LopezGonzalez, M. J. 2003, ASP Conf. Ser., 292, 369 [Google Scholar]
 Gokhale, M. H., & Zwaan, C. 1972, Sol. Phys., 26, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Gray, D. F. 2005, The Observation and Analysis of Stellar Photospheres (Cambridge: Cambridge University Press 2005) [Google Scholar]
 Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hadrava, P. 1995, A&AS, 114, 393 [NASA ADS] [Google Scholar]
 Hall, D. S., & Neff, S. G. 1976, IAU Symp., 73, 283 [Google Scholar]
 Hall, D. S., Keel, W. C., & Neuhaus, G. H. 1976, Acta Astron., 26, 239 [Google Scholar]
 Hatzes, A. P., & Zechmeister, M. 2007, ApJ, 670, L37 [NASA ADS] [CrossRef] [Google Scholar]
 Hatzes, A. P., & Zechmeister, M. 2008, J. Phys. Conf. Ser., 118, 012016 [NASA ADS] [CrossRef] [Google Scholar]
 Hatzes, A. P., Zechmeister, M., Matthews, J., et al. 2012, A&A, 543, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ibanoǧlu, C., Soydugan, F., Soydugan, E., & Dervişoğlu, A. 2006, MNRAS, 373, 435 [NASA ADS] [CrossRef] [Google Scholar]
 Ibanoǧlu, C., Dervişoğlu, A., Ćakırlı, Ö., Sipahi, E., & Yüce, K. 2012, MNRAS, 419, 1472 [NASA ADS] [CrossRef] [Google Scholar]
 Ilijic, S., Hensberge, H., & Pavlovski, K. 2001, Astrotomography: Indirect Imaging Methods in Observational Astronomy, eds. H. M. J. Boffin, D. Steeghs, & J. Cuypers (Berlin: Springer), 573, 269 [CrossRef] [Google Scholar]
 Jerzykiewicz, M., Handler, G., Shobbrook, R. R., et al. 2005, MNRAS, 360, 619 [NASA ADS] [CrossRef] [Google Scholar]
 Kalimeris, A., RovithisLivaniou, H., & Rovithis, P. 1994, A&A, 282, 775 [Google Scholar]
 Khalesseh, B., & Hill, G. 1992, A&A, 257, 199 [NASA ADS] [Google Scholar]
 Kolbas, V., Dervişoğlu, A., Pavlovski, K., & Southworth, J. 2014, MNRAS, 444, 3118 [NASA ADS] [CrossRef] [Google Scholar]
 Kreiner, J. M. 2004, Acta Astron., 54, 207 [NASA ADS] [Google Scholar]
 Kupka, F. G., Ryabchikova, T. A., Piskunov, N. E., Stempels, H. C., & Weiss, W. W. 2000, Balt. Astron., 9, 590 [Google Scholar]
 Lehmann, H., & Mkrtichian, D. E. 2004, A&A, 413, 293 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lehmann, H., & Mkrtichian, D. E. 2008, A&A, 480, 247 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lehmann, H., Tkachenko, A., & Mkrtichian, D. E. 2009, Commun. Asteroseismol., 159, 45 [CrossRef] [Google Scholar]
 Lehmann, H., Southworth, J., Tkachenko, A., & Pavlovski, K. 2013, A&A, 557, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lehmann, H., Tsymbal, V., Pertermann, F., et al. 2018, A&A, 615, A131 [CrossRef] [EDP Sciences] [Google Scholar]
 Lenz, P., & Breger, M. 2005, Commun. Asteroseismol., 146, 53 [NASA ADS] [CrossRef] [Google Scholar]
 Löptien, B., Lagg, A., van Noort, M., & Solanki, S. K. 2018, A&A, 619, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lovekin, C. C., & Deupree, R. G. 2008, ApJ, 679, 1499 [NASA ADS] [CrossRef] [Google Scholar]
 McLaughlin, D. B. 1924, ApJ, 60, 22 [Google Scholar]
 Mennekens, N., De Greve, J.P., van Rensbergen, W., & Yungelson, L. R. 2008, A&A, 486, 919 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mkrtichian, D. E., Kusakin, A. V., Gamarova, A. Y., & Nazarenko, V. 2002, in Astronomical Society of the Pacific Conference Series, Vol. 259, IAU Colloq. 185: Radial and Nonradial Pulsationsn as Probes of Stellar Physics, ed. C. Aerts, T. R. Bedding, & J. ChristensenDalsgaard, 96 [Google Scholar]
 Mkrtichian, D. E., Nazarenko, V., Gamarova, A. Y., et al. 2003, ASP Conf. Ser., 292, 113 [Google Scholar]
 Mkrtichian, D. E., Kusakin, A. V., Rodriguez, E., et al. 2004, A&A, 419, 1015 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mkrtichian, D. E., Kim, S. L., Rodríguez, E., et al. 2007, ASP Conf. Ser., 370, 194 [Google Scholar]
 Mkrtichian, D. E., Lehmann, H., Rodríguez, E., et al. 2018, MNRAS, 475, 4745 [NASA ADS] [CrossRef] [Google Scholar]
 Narusawa, S.Y., Nakamura, Y., & Yamasaki, A. 1994, AJ, 107, 1141 [NASA ADS] [CrossRef] [Google Scholar]
 Ohshima, O., Narusawa, S.Y., Akazawa, H., et al. 1998, Inf. Bull. Variable Stars, 4581, 1 [Google Scholar]
 Ohshima, O., Narusawa, S.y., Akazawa, H., et al. 2001, AJ, 122, 418 [NASA ADS] [CrossRef] [Google Scholar]
 Olson, E. C. 1982, ApJ, 259, 702 [NASA ADS] [CrossRef] [Google Scholar]
 Pamyatnykh, A. A., Dziembowski, W. A., Handler, G., & Pikall, H. 1998, A&A, 333, 141 [NASA ADS] [Google Scholar]
 Prša, A., & Zwitter, T. 2005, ApJ, 628, 426 [NASA ADS] [CrossRef] [Google Scholar]
 Przybilla, N., Firnstein, M., Nieva, M. F., Meynet, G., & Maeder, A. 2010, A&A, 517, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Raskin, G., van Winckel, H., Hensberge, H., et al. 2011, A&A, 526, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reed, M. D., Brondel, B. J., & Kawaler, S. D. 2005, ApJ, 634, 602 [NASA ADS] [CrossRef] [Google Scholar]
 Richards, M. T., Cocking, A. S., Fisher, J. G., & Conover, M. J. 2014, ApJ, 795, 160 [NASA ADS] [CrossRef] [Google Scholar]
 Rodríguez, E., García, J. M., Mkrtichian, D. E., et al. 2004, MNRAS, 347, 1317 [NASA ADS] [CrossRef] [Google Scholar]
 Rossiter, R. A. 1924, ApJ, 60, 15 [Google Scholar]
 RovithisLivaniou, H. 2001, Odessa Astron. Publ., 14, 91 [Google Scholar]
 Sablowski, D. P., Järvinen, S., & Weber, M. 2019, A&A, 623, A31 [CrossRef] [EDP Sciences] [Google Scholar]
 Saio, H. 1981, ApJ, 244, 299 [NASA ADS] [CrossRef] [Google Scholar]
 Schlesinger, F. 1910, Publications of the Allegheny Observatory of the University of Pittsburgh (Nabu Press), 1, 33 [NASA ADS] [Google Scholar]
 Shulyak, D., Tsymbal, V., Ryabchikova, T., Stütz, C., & Weiss, W. W. 2004, A&A, 428, 993 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Simon, K. P., & Sturm, E. 1994, A&A, 281, 286 [NASA ADS] [Google Scholar]
 Smith, M. A., & Gray, D. F. 1976, PASP, 88, 809 [NASA ADS] [CrossRef] [Google Scholar]
 Sterne, T. E. 1941, Proc. Natl. Acad. Sci., 27, 175 [NASA ADS] [CrossRef] [Google Scholar]
 Sybilski, P., Konacki, M., Kozłowski, S. K., & Hełminiak, K. G. 2013, MNRAS, 431, 2024 [NASA ADS] [CrossRef] [Google Scholar]
 Tkachenko, A. 2015, A&A, 581, A129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tkachenko, A., Lehmann, H., Tsymbal, V., & Mkrtichian, D. E. 2008, in 15th Proceedings of Young Scientists, eds. V. Y. Choliy & G. Ivashchenko, 33 [Google Scholar]
 Tkachenko, A., Lehmann, H., & Mkrtichian, D. E. 2009, A&A, 504, 991 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tkachenko, A., Lehmann, H., & Mkrtichian, D. 2010, AJ, 139, 1327 [NASA ADS] [CrossRef] [Google Scholar]
 Tkachenko, A., Van Reeth, T., Tsymbal, V., et al. 2013, A&A, 560, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tkachenko, A., Degroote, P., Aerts, C., et al. 2014, MNRAS, 438, 3093 [Google Scholar]
 Tsymbal, V. 1996, ASP Conf. Ser., 108, 198 [Google Scholar]
 Unno, W., Kiguchi, M., & Kitamura, M. 1994, PASJ, 46, 613 [Google Scholar]
 Van Hamme, W., & Wilson, R. E. 2005, Ap&SS, 296, 121 [CrossRef] [Google Scholar]
 Wilson, R. E., & Devinney, E. J. 1971, ApJ, 166, 605 [NASA ADS] [CrossRef] [Google Scholar]
 Wilson, A., & Maskelyne, N. 1774, Phil. Trans. R. Soc. London Ser. I, 64, 1 [Google Scholar]
 Yakobchuk, T. M., & Berdyugina, S. V. 2018, A&A, 613, A7 [CrossRef] [EDP Sciences] [Google Scholar]
 Zejda, M., Mikulášek, Z., & Wolf, M. 2008, A&A, 489, 321 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zima, W. 2008, Commun. Asteroseismol., 157, 387 [NASA ADS] [Google Scholar]
 Zucker, S., & Mazeh, T. 1994, ApJ, 420, 806 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
All Tables
Journal of observations listing the instrument, its spectral resolving power, year of observation, and mean Julian date.
Atmospheric parameters of primary and secondary of RZ Cas derived with multiple methods.
Orbital periods P_{RV} and times of minimum T_{RV} derived from the RVs from single seasons, and period changes △P_{phot} and periods P_{phot} derived from the photometric T_{Min}.
Synchronisation factor of the primary, radii of both spots on the secondary, longitude of Spot2 on the secondary, and mean scatter of the RV residuals.
Timescales in years determined from the seasonal variations of different parameters.
All Figures
Fig. 1 Continuum flux ratios between the components vs. wavelength. The dotted lines are calculated from synthetic continuum spectra with T_{eff1} of 8700 K and T_{eff2} of 3800–5000 K. The best agreement with the flux ratio from photometry (red line) is obtained for T_{eff2} = 4400 K (black solid line). The upper green line represents T_{eff1} = 9000 K, T_{eff2} = 4500 K, the lower green line for T_{eff1} = 8400 K, T_{eff2} = 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), bestfitting 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 R_{2} ∕R_{1} = 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 R_{2}∕R_{1}=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 R_{2} ∕R_{1} = 1.17 (red) and on the adjusted ratio of R_{2}∕R_{1} = 0.87 (green). Bottom row: ratio of the photometric flux ratio to the flux ratio obtained from spectrum analysis based on R_{2} ∕R_{1} = 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. outofeclipse 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 T_{min} as derived from the RVs in 2006. 

In the text 
Fig. 10 O–C values calculated from the corrected T_{Min}, 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: RossiterMcLaughlin effect in the RVs of the primary in different years (red) after subtracting the bestfitting 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 bestfitting 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 bestfitting periods and longterm trends. The dotted lines show possible sinusoidal variations with a period of 9 yr plus a longterm 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 10^{11} cm^{−3}) obtained from 3D hydrodynamic simulations of the RZ Cas system (Nazarenko & Mkrtichian, priv. comm.) for masstransfer 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.