© A. Tkachenko et al. 2022

1 Introduction

Stellar astrophysics has received an impressive boost with the launch and successful in-orbit operations of space missions such as the Microvariability and Oscillations of Stars (MOST, Walker et al. 2003), Convection, Rotation, and planetary Transits (CoRoT, Auvergne et al. 2009), Kepler/K2 (Borucki et al. 2010; Howell et al. 2014), and Transiting Exoplanet Survey Satellite (TESS Ricker et al. 2015) missions. While the first three missions have now been retired, TESS remains currently operational and has the potential to continue its science operations until 2030–2035. Furthermore, the PLAnetary Transits and Oscillations of stars (PLATO, Rauer et al. 2014) European Space Agency mission is being planned for launch in 2026 and will continue the tradition of collecting large amounts of nearly uninterrupted, long time-based, ultra-high precision, and duty cycle space-based photometry of stars and stellar systems.

Space-based photometric observations of stars and their ensembles are often complemented with ground-based spectroscopic observations, enabling stellar astrophysics studies of an unprecedented level of precision. For example, it has recently been demonstrated with asteroseismology (Aerts et al. 2010) that models of stellar interior structure and evolution have a largely incomplete treatment of angular momentum transport for intermediate- to high-mass stars (e.g. Mosser et al. 2012; Beck et al. 2012; Cantiello et al. 2014; Kawaler 2015; Aerts et al. 2017, 2019; Aerts 2021). These findings have in turn driven new theoretical developments that concern the physics of internal gravity waves (e.g. Rogers et al. 2013; Rogers 2015; Fuller et al. 2015a) and internal magnetic fields (e.g. Fuller et al. 2015b, 2019).

Parallel to detailed asteroseismic studies of the interior structures of stars, the field of binary stars and multiple systems has also received an unprecedented boost thanks to the abovementioned space missions. Nowadays, binary stars of various orbital and stellar type configurations are being found and observed in large numbers, many of which also contain intrinsically variable stellar components (e.g. Slawson et al. 2011; Matijevič et al. 2012; Conroy et al. 2014; Kirk et al. 2016; Sekaran et al. 2020; IJspeert et al. 2021). The co-existence of intrinsic variability of stars and binarity enables powerful synergies where, in particular, asteroseismic studies benefit from knowledge of highly precise fundamental properties of stars that come from binary dynamics (e.g. Schmid et al. 2015; Schmid & Aerts 2016; Johnston et al. 2019; Sekaran et al. 2021), while new aspects of tidal evolution theory can be developed and observationally verified thanks to interactions between stellar oscillations, magnetism, and binarity. Theory of tidal excitations of stellar pulsations was recently tested and improved thanks to the discovery and detailed observational studies of tidally induced pulsations in close binary systems (e.g. Welsh et al. 2011; Fuller & Lai 2012; Fuller et al. 2013; Fuller 2017). At the same time, the effects of tides on self-excited stellar pulsations are also being actively studied and their interactions are being quantified (e.g. Bowman et al. 2019; Handler et al. 2020; Fuller et al. 2020; Fuller 2021). Furthermore, it has recently been shown in observational incidences that some 5–10% of global magnetic fields among intermediate- to high-mass isolated stars are presumably of fossil origin (e.g. Alecian et al. 2019; Mathys 2017). Yet, there is a dearth of magnetic intermediate- to high-mass stars in close, short-period binaries, which is now claimed to be due to the dissipation of the fossil magnetic fields by the turbulent magnetic diffusion induced by the saturated tidal flows (Vidal et al. 2019).

Besides serving as natural laboratories to study physics of interactions, binary stars are essential objects in galactic and stellar evolution, and in distance scale contexts. First and foremost, the percentage of stars found in binary and higher-order multiple systems ranges from some 40% for solar-like stars (Cohen et al. 2020a) up to 55–60% for intermediate- to high-mass stars (e.g. Cohen et al. 2020b; Bodensteiner et al. 2021; Banyard et al. 2022). Hence, any galactic-scale study that relies on the stellar component has to take the phenomenon of binarity into consideration. Furthermore, spectroscopic double-lined (SB2) eclipsing binaries are a prime source of precise and accurate fundamental properties of stars. With the stellar masses and radii being measured with an accuracy better than 3% (e.g. Torres et al. 2010; Tkachenko et al. 2014; Debosscher et al. 2013; Pavlovski et al. 2009, 2014, 2018; Southworth 2013, Southworth 2021), such binary systems provide one of the most stringent observational tests of stellar structure and evolution models. It has been demonstrated that low-mass stars in binaries tend to have radii that are inflated by some 5–15% compared to the respective non-magnetic model predictions (e.g. Ribas 2006). Inclusion of the magnetic fields in these models helps to resolve the issue (e.g. MacDonald & Mullan 2017), demonstrating the efficiency of binaries in probing and improving stellar models in the regime of low-mass (M-type) stars. Similarly, binarity is being exploited to test interior physics in the models of intermediate- to high-mass stars in a way that is independent of and complementary to the method of asteroseismology. Several studies offer probes into the level of core-boundary mixing in the form of the convective core overshooting (Guinan et al. 2000; Claret & Torres 2016, 2017, 2018, 2019; Martinet et al. 2021), while others go as far as to report the pertinent need for higher convective core masses in models of intermediate- to high-mass stars (e.g. Tkachenko et al. 2020; Johnston 2021). In addition, interferometric eclipsing SB2 systems provide the means of measuring their accurate dynamical parallaxes, thus offering an important reference for astrometric missions such as Hipparcos and Gaia (e.g. Gallenne et al. 2019; Pavlovski et al. 2022). Last but not least, eclipsing (binary) and transiting (exoplanet) systems offer a unique opportunity to study the effects of spin-orbit (mis)alignment through observations and interpretation of the Rossiter-McLaughlin (RM) effect (e.g. Rossiter 1924; McLaughlin 1924; Albrecht et al. 2012; Triaud 2018; Wang et al. 2018; Addison et al. 2018; Kamiaka et al. 2019).

With the amount and, more importantly, the unprecedented quality of (space-based) photometric and (ground-based) spectroscopic observations available nowadays, considerable time and effort are being spent to improve methods and techniques for the analysis of stars in binary and higher-order multiple systems. For example, following the development of binary light curve analysis algorithms such as Wilson-Devinney (Wilson & Devinney 1971), Phoebe (Prša & Zwitter 2005), and jktebop (Southworth et al. 2004) in the pre-space missions era, new methods such as those implemented in the ellc (Maxted 2016) and Phoebe 2.0 (Prša et al. 2013; Degroote et al. 2013; Prša 2018; Conroy et al. 2019, 2020, 2021) codes are being developed, which are both more efficient and often include a better treatment of physics. In spectroscopy, the original method of cross-correlation (e.g. Bracewell 1965) has commonly been used in the past, but it proved impractical in many cases due to its strong dependence on the underlying template (synthetic) spectrum and its inability to pick out line profile variations intrinsic to the star (unless they are exceptionally strong, e.g. in large amplitude radially pulsating stars). Yet, two-dimensional cross-correlation analysis (e.g. as implemented in the TODCOR software package, Mazeh & Zucker 1994) is often used for the RV determination in SB2 systems, despite the above-mentioned dependence on the underlying model and the number of degeneracies to which this method is subject. A considerable improvement has been achieved with the introduction of the broadening function method (Rucinski 1992), which is different from the original method of cross-correlation in that it uses a template that is not subject to spectral line broadening mechanisms such as stellar rotation and/or macro-turbulence. This makes the method less model-dependent in the sense that small deviations of the template parameters from the true properties of the star are less critical for the overall shape of the obtained broadening function. In addition, the use of a template that incorporates only the intrinsic (e.g. thermal) broadening of spectral lines allows us to resolve fine structures in the ensemble of lines in the observed stellar spectrum, and enables the detection and further analysis of line profile variations caused by the intrinsic variability of the star.

The need to achieve an extremely high signal-to-noise ratio (S/N) to be able to detect low-amplitude magnetic signatures in polarisation spectra has driven development of the least-squares deconvolution (LSD; Donati et al. 1997) multi-line averaging technique. The LSD profiles obtained with this method are also ideal for inference of radial velocities (RVs) of stars (including those in binary systems) and for studying temporal spectral changes due to the intrinsic variability of stars. In this paper, we present a generalisation of the original LSD method for binary stars. In particular, we focus on the determination of precise, orbital-phase-resolved radial velocities (RV) of the individual binary components with a careful treatment of the RM effect. In Sect. 2, we provide a short introduction to the original LSD method, as well as the literature overview with respect to the recent improvements and generalisations of the method. Our own generalisation of the method to spectroscopic double-lined binaries is described in detail in Sect. 3. The method is further tested and validated on artificial spectra of SB2 stars in Sect. 4, where we also touch upon its applicability range and limitations. We close the paper with conclusions and an outline of future prospects in Sect. 5.

2 The least-squares deconvolution (LSD) method

The LSD method was originally introduced by Donati et al. (1997) with the purpose of detecting and measuring weak surface magnetic fields from polarised spectra of stars. In its original formulation, the method relies on the following two fundamental assumptions: (1) all spectral lines in the considered spectral interval have a similar profile shape, and (2) all lines in the spectrum add up in a strictly linear fashion. Assuming both conditions are adequate for representing a stellar spectrum, the LSD problem can be formulated as a convolution of an a priori unknown average profile Z(υ), with a pre-computed line mask M: $$I=M*Z\left(\upsilon \right),$$(1)

where υ and I stand for the velocity and model spectrum, respectively, while the line mask M represents a set of delta functions (wavelength vs. predicted central line depths). The concept of representing observed spectrum as a convolution of an a priori unknown average profile with a list of delta functions allows the S/N to be significantly boosted. The exact amount of the S/N boost depends on the spectral type of the star and on the wavelength range employed for the calculation of the LSD profiles, and is typically anywhere between a factor of five to 50. Such a large increase in the S/N of the data in turn enables the detection of stellar surface magnetic fields as weak as a few tens of Gauss from the originally moderate S/N polarisation spectra (e.g. Blazère et al. 2016a,b; Kochukhov et al. 2019).

In practice, neither of the abovementioned fundamental assumptions is fulfilled. It is a well known fact that strong (close to resonance) lines differ substantially in shape from their weak counterparts, and that lines with pronounced damping wings (e.g. Balmer lines in the spectra of main-sequence intermediate-mass stars and/or magnesium lines in solar-like stars) exhibit shapes that are very distinct from any other type of spectral lines. Furthermore, it is fair to assume that individual spectral lines add up linearly when blending occurs due to a superposition of two stellar spectra (i.e. in the case of a double-lined spectroscopic binary) and/or due to a large rotational broadening. However, the above assumption is no longer valid when absorption coefficients of two or more neighbouring lines overlap in wavelength, which results in anon-linear addition of the lines, forming a blend. Both these limitations give rise to non-negligible uncertainties in the predicted intensities of spectral lines in the model spectrum I obtained by a convolution of the LSD profile Z(υ) with the input line mask M. This implies that LSD profiles computed from a time series of observed spectra can be studied with respect to the time-dependent variability they exhibit. However, each of these LSD profiles separately can hardly be interpreted as an isolated spectral line with properties averaged across the ensemble of spectral lines from which the LSD profile is computed (Kochukhov et al. 2010).

Several studies have suggested improvements to the original LSD technique with the ultimate goal of overcoming the major shortcomings of the method associated with its fundamental assumptions. Sennhauser et al. (2009) and Sennhauser & Berdyugina (2010) proposed using Minnaert’s analytical expression (Minnaert 1935) to account for non-linearity in spectral line blending when it comes to representation of strong, optically thick lines. The authors demonstrated that their method of non-linear deconvolution with deblending (NDD) offers a better representation of blended spectra than the original LSD method, which assumes linearity. Kochukhov et al. (2010) propose an improved method of least-squares deconvolution (iLSD), which includes the option of a multi-profile LSD that represents stellar spectrum as a superposition of an arbitrary number of scaled average profiles. The iLSD method offers the opportunity to group spectral lines according to their predicted depths and to represent each of the above groups with its own average profile. Given that weak and strong lines tend to show significant differences in their profile shapes, the approach proposed by Kochukhov et al. (2010) allows for partial compensation of one of the fundamental assumptions in the original formulation of the method, namely a shared profile shape between all spectral lines engaged in the calculation of the LSD profile. Furthermore, the authors demonstrate that the LSD profile cannot be interpreted as a real spectral line with average properties when it comes to Stokes I (intensity) and Stokes Q (linear polarisation) spectra. At the same time, the Stokes V average profiles do closely resemble the behaviour of an isolated line with average properties, provided that the global magnetic field is weaker than some 1 kG.

Tkachenko et al. (2013) elaborate further on the idea proposed by Kochukhov et al. (2010) and implement in their approach the possibility of including multiple line masks ${\displaystyle {\sum }_{i=1}^N{M}_i}$ along with a multi-profile LSD. In particular, such an approach allows one to compute composite average profiles from spectra of double-lined binary stars, taking into account differences in the atmospheric parameters and chemical compositions of the two binary components. Furthermore, the authors implement a line-strength correction algorithm that aims to improve upon the LSD representation of the stellar spectrum. Corrections are applied to the line strengths locally by minimising the differences between the observed and LSD-based model spectra. These corrections are purely mathematical and cannot be interpreted in terms of physics of stellar atmospheres, yet they allow for a significantly improved representation of the LSD-based model spectrum. This spectrum will have a much higher S/N than the original observed spectrum, and thus the procedure by Tkachenko et al. (2013) offers an efficient way of de-noising stellar spectra.

Lastly, Asensio Ramos & Petit (2015) implement the LSD algorithm under the Bayesian framework, which enables a fast calculation of the LSD average profile with a Gaussian prior. Strachan & Anglada-Escudé (2017) introduce a method of differential least-squares deconvolution (dLSD) that employs a high S/N combined spectrum of the star as a template and searches for a convolution kernel that needs to match the template to an observation of the star. The method is developed specifically for studying planet-star obliquity through the interpretation of the observed RM effect.

Fig. 1 LSDBinary algorithm flow chart. Bottom left panels show output of the LSDinit programme in terms of the LSD profiles computed from synthetic spectra (left) and the continuum flux ratio of the two binary components as a function of the wavelength (right). The output of the LSDBinary programme is illustrated in the two bottom right panels, showing LSD profiles of the individual binary components (left) and a comparison between the composite observed (black) and LSD-based model (red) spectra of the binary. See Sect. 3 for details.

3 LSDBinary: generalisation of the LSD method for spectroscopic double-lined binaries

The LSD method is, by design, sensitive to signals that are common to the majority of spectral lines in a stellar spectrum. Therefore, even under the assumption of a single-line mask, M, LSD remains a powerful method for the detection of SB2 binary systems. For example, the LSD profile computed from the composite spectrum of an SB2 system that employs a line mask corresponding to the primary component, will contain signatures of both binary components with a caveat that the contribution of the secondary component to the LSD profile will be somewhat diluted. Measuring RVs and/or studying line profile variations (LPVs) of individual binary components in these composite LSD profiles still proves nearly as difficult as with the original observed composite spectra because, despite the significantly enhanced S/N in the LSD profile, line blending remains the dominant source of uncertainty. In practice, this means that the original LSD method, as well as its numerous generalisations outlined in Sect. 2, despite being effective for the detection of SB2 systems, are not suitable for the calculation of LSD profiles of the individual binary components and inference of their precise RVs.

In this paper, we present a generalisation of the original LSD method applied to SB2 binary star systems. In particular, we focus on the separation of the spectral contributions of the individual binary components in velocity space to a level of precision that would enable inference of precise and accurate RVs of both stars and, if possible, enable their LPVs to be resolved. In doing so, we focus primarily on the difficult case of the in-eclipse orbital phases that are characterised by a high degree of blending, while we pay little attention to the out-of-eclipse phases where spectral contributions of the two stars can typically be easily separated from each other. This specific focus on the in-eclipse spectra is done solely for the purpose of demonstrating the algorithm’s performance in the most difficult circumstances, while the algorithm is designed to handle all orbital phases in the spectroscopic observations of eclipsing SB2 systems. As a baseline, we adopt the modified LSD method by Tkachenko et al. (2013), which employs a modified, fast version of the Levenberg-Marquardt algorithm (Marquardt 1963) implemented by Piskunov & Kochukhov (2002). The algorithm has the capabilities outlined in the previous section and serves as the ‘central engine’ of the LSDBinary software package presented here.

3.1 Pseudo-code

Figure 1 provides a flow chart of the LSDBinary software package with the three main components - input, central engine, and output - being indicated. Here, we provide a more detailed description of the entire process, in the form of a pseudo-code.

3.1.1 The LSDInit algorithm

The LSDinit algorithm comprises three (optionally four) major steps, whose purpose is to take care of all of the necessary preparations for the core calculations with the LSDBinary algorithm. In this module, we set up initial guesses for the individual LSD profiles and the wavelength-dependent flux ratio of the two stars, and we pre-compute local corrections to the LSD model spectrum to ensure its closest match to the observations. The three (optionally four) above-mentioned steps are described in detail below.

The SynthV radiative transfer code (Tsymbal 1996) is employed to compute synthetic spectra in an arbitrary wavelength range, with options for the variable micro-turbulent velocity and chemical composition of the star, including the possibility for their vertical stratification. Spectra are synthesised for different positions on the stellar disk to account for centre-to-limb intensity variations, also known as the limb darkening effect. The SynthV code is publicly available and its latest version is distributed as part of the LSDBinary package. However, SynthV can also be replaced with any other preferred radiative transfer code, provided its output is tuned to the input requirements of the LSDBinary code.

The Convolve code (Tsymbal 1996) performs the integration of specific intensities over the visible stellar disk and the convolution of the obtained spectrum with the Gaussian profile with the full width at half maximum (FWHM) corresponding to resolving power R of the instrument. In addition, the convolution is performed with the projected rotational velocity v sin i of the star, and optionally the macro-turbulent velocity. The macroturbulent broadening is implemented in the radial-tangential formalism, following Gray (1992).

The Rossiter-McLaughlin Effect (rme) algorithm can be used as an alternative to calculate the disk-integrated synthetic spectra with the SynthV and Convolve suit of codes, and is therefore an optional step. The rme algorithm, as the name suggests, allows us to compute time series of the in-eclipse synthetic spectra subject to distortion caused by the RM effect (Rossiter 1924; McLaughlin 1924) and the variable associated with the orbital-phase binary components’ flux ratio. RME is a geometric effect that occurs when the (stellar or planetary) companion transits across the disk of the rotating primary component and blocks part of its visible surface. In doing so, the companion also blocks part of either the blue- (ingress) or red- (egress) shifted (due to stellar rotation) light of the primary component, causing asymmetries in its observed line profiles. Our rme algorithm solves the problem analytically by computing the surface area of the eclipsed star that is being blocked by the companion and performs the disk integration over the remaining visible part of its disk. The rme algorithm relies on orbital configuration of the system as input, and delivers distorted line profiles and the orbital-phase-dependent components’ flux ratio, to be used for the calculation of the initial guess LSD profiles (see below) for the in-eclipse phases. Therefore, the option of using the rme algorithm is most suitable for binary systems whose orbital configuration is known in advance, though the algorithm does not require the orbital parameters to be known with particularly high precision and/or accuracy.

The LSDsynth algorithm is designed to compute theoretical LSD profiles for both the binary components from the corresponding disk-integrated synthetic spectra (with or without the RME taken into account) and line masks. To speed up calculations, as the initial guess for the LSD profile, the code employs either a Gaussian profile with the FWHM corresponding to the resolving power R of the instrument, or a rotational kernel for the cases of slow and moderate to rapid rotation, respectively. An example of the synthetic LSD profiles computed with the LSDsynth algorithm for both binary components is shown in Fig. 1 (bottom left panel). As discussed in Sect. 2, due to the natal assumptions and limitations of the original LSD method, the LSD-based model spectrum I, computed by means of the convolution of the LSD profile Z(v) with the line mask M, fails to closely reproduce the input stellar spectrum from which the LSD profile is computed. The LSDsynth algorithm developed by us overcomes that problem by computing and applying local corrections to line intensities in the LSD-based model spectrum in order to provide its closest match to the input synthetic stellar spectrum. These relative intensity corrections take the following form:(r_synthetic − r_model)/r_model, where r stands for the normalised flux, while ‘synthetic’ and ‘model’ refer to the input synthetic and LSD-based model stellar spectrum, respectively. Therefore, these local intensity corrections represent fractional differences between the input synthetic and LSD-based model spectra at each wavelength point, and ensure the closest match between the two when applied to the LSD-based model spectrum. Lastly, the LSDsynth algorithm computes the ratio of the binary components’ fluxes as a function of wavelength and approximates it with the second-degree polynomial. An example of the functional form of this wavelength-dependent flux ratio is shown in Fig. 1 (bottom left panel). In summary, taking synthetic spectra and line masks as input, the LSDsynth algorithm delivers: (i) initial guess synthetic spectra-based LSD profiles, (ii) local fractional intensity corrections to the LSD-based model spectrum, and (iii) best-fit polynomial coefficients that allow one to reproduce the wavelength-dependent functional form of the components’ flux ratio.

3.1.2 The LSDBinary algorithm

The LSDBinary algorithm represents the ‘core engine’ of the method and allows us to compute time series of the LSD profiles of both binary components separated in velocity space. The algorithm has three major components discussed below.

‘Input’: The algorithm relies on the output of the LSDinit module in terms of the initial guess for: (i) the synthetic spectrum-based LSD profiles of both binary components, (ii) local fractional intensity corrections to the LSD-based model spectra of both stars, and (iii) the functional form of the components’ flux ratio. In addition to the above, the LSDBinary algorithm also requires the initial guess for the RVs of both binary components, which can be provided either in the form of a table, RV versus orbital phase, or in the form of orbital parameters of the binary system. In the latter case, orbital-phase-resolved RVs are computed from the provided orbital elements within the LSDBinary module and are employed as the initial guess. Lastly, LSDBinary requires two line masks, one for each binary component, to be provided as well.

‘Core calculations’: LSDBinary utilises the LSD algorithm implemented in Tkachenko et al. (2013), in particular its functionality, to compute multi-component LSD profiles with the employment of multiple line masks. We modify the algorithm in order to account for relative contributions of the two binary components to the system’s composite spectrum. These contributions are defined by the components’ flux ratio, whose wavelength-dependent functional form is pre-computed within the LSDinit module and their radii ratio. The latter parameter can either be fixed (e.g. if it is known from the binary light curve solution) or used as a free parameter in LSDBinary. The problem of solving for the individual LSD profiles is approached in an iterative fashion with the employment of the modified version of the Levenberg-Marquardt algorithm (see above and Piskunov & Kochukhov 2002). In each iteration, the algorithm minimises the difference between the input observed composite spectrum of the binary system and its LSD-based composite model spectrum counterpart. The latter is computed as a linear sum of the individual binary components’ model spectra, which are subject to their local fractional intensity corrections before being added up.

‘Output’: The LSDBinary algorithm has Level-0 and Level-1 outputs. The former refers to the individual LSD profiles and the LSD-based model spectra. These model spectra are provided for: (i) both binary components separately, and are the product of the convolution of the primary and secondary LSD profiles with their respective line masks, and taking the individual fractional intensity corrections into account; and (ii) the binary system as a whole, where the composite model spectrum is represented by the linear sum of the primary and secondary LSD-based model spectra, with the components’ flux and radii ratio being taken into account. On the other hand, Level-1 output refers to the quantities inferred directly from the LSD profiles and the process of their calculation (i.e. from the Level-0 output). These quantities are the RVs of both binary components and, where it is set as a free parameter, the radii ratio of the two stars as a function of the orbital phase. Individual RVs are computed by matching the components’ LSD profiles to the respective initial guess synthetic spectrum-based LSD profiles with the employment of the golden search algorithm (Kiefer 1953). Therefore, RVs delivered by the LSDBinary algorithm are on the scale relative to the initial guess LSD profiles. An example of Level-0 output in terms of the individual LSD profiles and the LSD-based composite model spectrum overlaid on the input observed spectrum of the binary system is shown in Fig. 1 (bottom right panel).

Fig. 2 Graphical presentation of the LSDBinary algorithm. Left: initial guess synthetic LSD profiles computed with the LSDinit module, with (red stars) and without (black dots) the RM effect, along with the resulting ‘observed’ LSD profile (solid black line) computed with the LSDBinary module. Middle: in-eclipse phase-resolved RV curves, with (red stars) and without (black dots) the RM effect taken into account in the computation of the initial guess LSD profiles. The RVs measured relative to the respective initial guess LSD profiles are indicated with the arrows. Right: radii ratio as a function of orbital phase, with (red stars) and without (black dots) the RM effect taken into account in the computation of the initial guess synthetic spectrum-based LSD profiles. Arrows indicate the same data points as in the middle panel. See text for details.

3.2 Step-by-step algorithm demonstration

In order to provide a step-by-step demonstration of the LSDBinary algorithm performance, we simulate the in-eclipse spectra of an Algol-like binary system with a mid-A-type primary and a K-giant secondary as components. Time series of orbital-phase-resolved input spectra are simulated with our RME module for a representative S/N of ~100. In this ideal scenario, line masks are taken to fully correspond to the atmospheric parameters of the primary and secondary components. The initial guess and synthetic spectrum-based LSD profiles, as well as the functional form of the wavelength-dependent flux ratio are computed with the LSDInit module under the following assumptions: (i) LSD profiles of both stars are symmetric and their flux ratio is constant with orbital phase, and (ii) the LSD profile of the primary component is distorted due to the RM effect and the components’ flux ratio is orbital-phase dependent. In this particular example, as well as in the rest of this work, LSD profiles are computed from a 500 Å wide wavelength interval centred at 5250 Å. This wavelength interval is free of the Balmer lines and contains a large amount of metal lines in the spectra of A- and F-type stars that we are concerned with in this work.

The left panel in Fig. 2 shows a singleshot system’s LSD profile (solid black line) computed with the LSDBinary module at orbital phase ϕ = 0.03 (ϕ = 0.0 corresponds to the centre of the primary eclipse). The initial guess synthetic spectrum-based symmetric LSD profile (computed with the SynthV and Convolve suite of codes) and the initial guess LSD profile distorted by the RM effect (computed with the RME module), are shown with the black dots and red stars, respectively. We stress that, irrespective of the shape of the initial guess LSD profile, the output ‘observed’ LSD profile is highly asymmetric and has a shape that is remarkably similar in both cases. In other words, the Level-0 data products of the LSDBinary algorithm in the form of the system’s LSD profile and the corresponding LSD-based model spectrum are barely sensitive to the assumed initial guess. Radial velocities computed from this singleshot system’s LSD profile and measured relative to the respective initial guess synthetic spectrum-based LSD profiles (that are computed at laboratory wavelengths) are indicated with arrows in the middle panel in Fig. 2 (designated as ‘RV_noRME’ and ‘RV_RME’ for the cases of the symmetric and the distorted initial guess profiles, respectively). The middle panel also demonstrates the orbital-phase-resolved RV curves for clarity of comparison between the two initial guess cases. We note that the two RV curves differ from each other substantially, with one of them resembling the pure orbital motion of the primary component (red symbols), whereas the other one shows signatures of the line-profile distortions due to the RM effect (black symbols). Because RV is an integral quantity of the line profile, it is sensitive to both the global shift of the profile due to the orbital motion of the star and to the profile distortion, no matter what the true cause of that distortion is. In our example (see left panel in Fig. 2), the RV of the ‘observed’ LSD profile has a positive difference of some 15 km s⁻¹ with respect to the symmetric initial guess LSD profile. On the other hand, the same ‘observed’ LSD profile has a negative RV difference of some 10 km s⁻¹ with respect to the asymmetric initial guess LSD profile, which shares laboratory wavelengths with its symmetric version but now contains distortions due to the RM effect.

The effect of the initial guess LSD profile and the corresponding orbital-phase-resolved functional form of the components’ flux ratio on the inferred radii ratio parameter is illustrated in the right panel in Fig. 2. One can see that, in the case of a distorted initial guess LSD profile due to the RM effect initial guess, the radii ratio curve appears nearly flat (red stars). This result is indicative of LSDBinary finding a constant with an orbital-phase value of the radii ratio that, in this particular case, also converges to the assumed input value of R₂/R₁ = 1.2 in the simulations. This is because the flux contribution of the primary component varies over the course of its eclipse and this variability is taken into account in the calculation of the initial guesses with the RME module. Threfore, the primary’s spectrum is modelled correctly in the LSDBinary module. On the other hand, the assumption of a constant initial guess LSD profile and flux contribution of the primary component over the course of its eclipse leads to appreciable variations in the inferred radii ratio parameter (black dots). This result is indicative of an incorrectly assumed orbital-phase-independent flux contribution of the primary component. By assuming no flux variation for the primary component over the course of its obscuration by the secondary, we overestimate the primary’s contribution to the composite spectrum of the binary system. The amount of this overestimation is variable and depends strongly on the surface area of the primary that is being blocked by the companion star. In practical terms, the assumption of a constant with an orbital phase flux ratio leads to a significant overestimation of line depths for the primary and causes an appreciable mismatch between the observed and LSD-based model composite spectra of the binary system. In an attempt to compensate for this mismatch, the LSD-Binary module seeks to lower the contribution of the primary component to the total flux of the system, a goal that is most efficiently achieved by increasing the radii ratio parameter R₂/R₁. The effect is expected to be more pronounced at phases where a larger surface area of the primary is being blocked by the companion star, namely in the centre of the primary eclipse (just as we observe in the right panel in Fig. 2).

In conclusion, although the Level-0 data products of the LSDBinary algorithm (LSD profiles and LSD-based model spectra) are barely sensitive to the choice of the initial guess LSD profiles and constant versus the variable flux ratio, its Level-1 data products (RVs and the radii ratio parameter) are influenced significantly by the choice made. In the next section, we demonstrate how the sensitivity of the Level-1 data products to the choice of the initial guess can be exploited to improve upon the assumed orbital parameters of a binary system.

4 LSDBinary: method validation, applicability range, and limitations

In this section, we present a few key tests based on simulated time series of spectra, which are necessary to get a grasp on the applicability range of the LSDBinary algorithm. In these simulations, we consider two artificial binary systems that are representative of two particular classes of binary stars. Specifically, these are detached systems with similar stellar components and semi-detached Algol-type systems, with one component being significantly cooler and more evolved than the other. RZ Cas (e.g. Soydugan et al. 2006; Tkachenko et al. 2009) and KIC 11285625 (Debosscher et al. 2013) are used as reference systems in our calculations, and their relevant systemic and atmospheric parameters are summarised in Table 1. The RME module is used to synthesise the composite spectra of both binary systems, with a particular focus on the in-eclipse phases. These orbital phases are by far the most challenging because: (i) spectral-line contributions of the two binary components overlap in the velocity space, (ii) line profiles of the primary component are largely distorted due to the RM effect, and (iii) the flux contribution of the primary component is orbital-phase dependent.

4.1 Effects of the signal-to-noise ratio and resolving power

In the first instance, we test the ability of the LSDBinary method to handle spectra of variable S/N and resolving power R. To that end, we add Poisson noise to our synthesised time series to simulate spectra of S/N = 40, 80, and 120, for which we also consider three regimes of resolving power: low, medium, and high, at R = 5000, 25 000, and 60000, respectively. Figures 3 and 4 summarise the results obtained for both binary systems. One can see that the common feature to both systems is the ability of the LSDBinary algorithm to handle spectra in the entire range of S/N values: the obtained LSD profiles and the LSD-based model spectra have S/Ns that are (as expected) substantially higher than the input spectra, while the line-profile shapes do not display any significant dependency on the S/N value. This is not an unexpected result given that the LSD technique is, by design, most efficient for stars that display large number of lines in their spectra, and both our simulated binary systems with their A- and F-type primary components meet the above requirement. Furthermore, we do not observe any significant changes in the obtained Level-0 data products when degrading the resolving power from high (R = 60000) to medium (R = 25000) for the case of the RZ Cas system. However, we lose the ability to resolve the primary component’s spectral-line distortions due to the RM effect when degrading spectral resolution to its lowest value of R = 5000. Comparing the cases of RZ Cas and KIC 11285625, we find that the effect of spectral resolution is much more important in the latter case: local distortions of the LSD profile of the primary component in the KIC 11285625 system are already smeared out at a medium resolving power of R = 25 000. This is likely explained by the much smaller value of the projected rotational velocity of the star than in the case of the primary component of RZ Cas and, to a lesser extent, by relative contributions of the two stellar components of KIC 11285625 to its composite spectrum. Overall, the orbital and stellar configurations of the KIC 11285625 system suggest a much lower amplitude of the RM effect than in the case of RZ Cas, which in turn requires higher spectral resolution to resolve the effect.

Fig. 3 Results of the application of the LSDBinary algorithm to artificial in-eclipse spectra of the RZ Cas binary system. Top block: LSD profiles of the primary (red) and secondary (black) components. Bottom block: LSD-based model spectra (red) overlaid on the input ‘observed’ spectra (black) of the system. In each of the blocks, rows correspond to different values of the resolving power (from top to bottom, R = 60 000, 25 000, and 5000), while columns reflect the change in S/N of the spectrum (from left to right, S/N = 120, 80, and 40).

4.2 Effect of variable input atmospheric parameters

In the next test, we investigate the effect of changing the atmospheric parameters of either binary component on the Level-0 data products of the LSDBINARY algorithm. The effect of variable and/or uncertain atmospheric parameters is expected to propagate into the LSD profiles through the changing properties of the input line mask, as well as through the associated change in the components’ flux ratio. Since the effect is qualitatively similar for both simulated binary systems, here we present the results obtained for the RZ Cas system only. We also present the results for a single combination of R = 60 000 and a S/N = 120, since the effect of changing these characteristics has been discussed in the previous section.

Figures 5 and 6 illustrate the effect of changing the atmospheric parameters of the primary and secondary components, respectively, on the obtained LSD profiles. Effective temperature T_eff, surface gravity log g, and metallicity [M/H] are varied symmetrically with respect to their true values (see Table 1) and within typical spectroscopic uncertainties: ~5% in T_eff and 0.3–0.4 dex in log g and [M/H]. The LSD profile of the primary changes most notably when varying the component’s T_eff and [M/H] (top and bottom left panels in Fig. 5), while it appears to be barely sensitive to the variations in log g (middle left panel in Fig. 5). Lack of sensitivity to the latter parameter is likely explained by the limited effect of the star’s surface gravity on the overall shape and, to a lesser extent, strength of the metal lines in the T_eff regime of A- and F-type stars. At the same time, changing parameters of the primary component have little to no effect on the LSD profile of the secondary star (right column in Fig. 5).

Similarly, varying parameters of the secondary star has little to no effect on the LSD profile of the primary (left column in Fig. 6). At the same time, the LSD profiles of the secondary component appear to be most sensitive to changes in the T_eff (top right panel in Fig. 6) and log g (middle right panel in Fig. 6) of the star, and less sensitive to its metallicity [M/H] (bottom right panel in Fig. 6). We note that the secondary component in an Algol-type system is more evolved and about twice as cool as the primary, which explains a more appreciable sensitivity of its LSD profile to the log g parameter.

Ultimately, we stress that no significant changes in the overall shape of the LSD profile are observed for either of the binary components when their atmospheric parameters are varied within typical spectroscopic uncertainties. Instead, Figs. 5 and 6 reveal an appreciable effect on the global scaling of the component’s LSD profile in terms of its depth and strength. These findings are in line with the results of Tkachenko et al. (2013), who also demonstrate that ability of the LSD technique to reproduce an overall shape of the line profile is independent of the assumption of atmospheric parameters for the line mask calculation, provided they do not deviate from their true values by more than typical spectroscopic uncertainties.

Fig. 4 Same as for Fig. 3 but for the KIC 11285625 system.

Fig. 5 Effect of changing atmospheric parameters of the primary component of the RZ Cas system on the LSD profiles of the primary (left) and secondary (right). From top to bottom: Teff, log g, and [M/H] are varied between 8400 and 9200 K (step width of 100 K), 3.9 and 4.7 dex (step width of 0.1 dex), and −0.3 and 0.3 dex (step width of 0.1 dex), respectively. The LSD profiles shown in red correspond to the true parameters of the star listed in Table 1.

Fig. 6 Effect of changing atmospheric parameters of the secondary component of the RZ Cas system on the LSD profiles of the primary (left) and secondary (right). From top to bottom: Teff, log g, and [M/H] are varied between 4500 and 5100 K (step width of 100 K), 3.1 and 3.9 dex (step width of 0.1 dex), and −0.3 and 0.3 dex (step width of 0.1 dex), respectively. The LSD profiles shown in red correspond to the true parameters of the star listed in Table 1.

4.3 Effect of variable atmospheric parameters, radii ratio, and orbital inclination

In Sect. 3.1, we briefly touched upon how the LSDBINARY algorithm can be used in conjunction with the assumption of a symmetric initial guess LSD profile to measure the in-eclipse RV variations associated with both the orbital motion of the star and the geometric RM effect. Such RV curves can be used for follow-up studies, for instance to assess the degree of spin-orbit misalignment, if present in the system. At the same time, we showed that the assumption of an orbital-phase-independent flux ratio of the two stars (that comes with the assumption of a symmetric initial guess LSD profile) results in an appreciably variable with an orbital phase radii ratio R₂/R₁ delivered by LSDBinary. Obviously, the true radii of the stars (and hence their ratio) do not change in the course of the eclipse and the observed variation of R₂/R₁ is a pure effect of ignoring the orbital-phase dependency of the flux ratio in the initial guess calculations. In essence, the R₂/R₁ parameter is being exploited in this case as a global scaling factor for the light contribution of the primary component that varies over the course of its eclipse.

In this section, we investigate the effect of the variable (at the step of the initial guess calculation in the LSDInit module) stellar radii ratio and orbital inclination angle on the data products of the LSDBinary algorithm. Similar to Sect. 4.2, we focus on the case of an artificial Algol-type system resembling stellar and orbital properties of the RZ Cas binary system (see Table 1), where the artificial ‘observed’ spectra are characterised by R = 60 000 and a S/N = 120. To account for the RM effect in the test, the initial guess orbital-phase-resolved LSD profiles and wavelength-dependent functional form of the components’ flux ratio are computed with the RME module.

In the current test, the radii ratio R₂/R₁ and orbital inclination angle i parameters are varied one at the time while keeping the other one fixed. Thus, the two scenarios considered are: (i) the R₂/R₁ parameter is varied in the range from 0.8 to 1.6 (step width of 0.05) while keeping i = 82°, and (ii) the i parameter is varied in the range from 70° to 90° (step width of 2°) while keeping R₂/R₁ = 1.19. In total, we consider 28 R₂/R₁-i parameter combinations (17 and 11 in the cases (i) and (ii), respectively) when computing the initial guess LSD profiles and components’ flux ratio with the LSDInit module. The initial guess is then passed on to the LSDBinary algorithm that takes the above-described artificial ‘observed’ spectra as input and returns optimised LSD profiles of both binary components and their radii ratio parameter R₂/R₁.

The top row in Fig. 7 summarises the results obtained for the first case (i), where the initial guess radii ratio parameter R₂/R₁ is varied while keeping the inclination angle i fixed. The left and right panels demonstrate the orbital-phase-resolved R₂/R₁ and RV curves, respectively, where the parameter combination that represents the closest match to the true parameter values is shown with red stars. One can see that this particular solution favours a constant with the orbital-phase value of the radii ratio and the RV curve that closely resembles pure orbital motion of the primary component. At the same time, setting the initial guess R₂/R₁ to a larger or smaller parameter value (e.g. 1.4 or 1.1 as indicated with the filled triangles and squares in the top left panel in Fig. 7) results in the R₂/R₁ versus orbital-phase curve that deviates appreciably from the constant value of the respective parameter. The same conclusion holds true for the RV curve where one starts to see a residual signal from the RM effect, which increases progressively as the deviations of the initial guess R₂/R₁ parameter get larger from its true value. The former is the direct consequence of an incorrectly assumed initial guess components’ flux ratio due to variable R₂/R₁ and fixed effective temperatures, while the latter is explained by the variations in LSD profiles obtained with LSDBinary and associated with the changing RM effect geometry due to variations in R₂/R₁.

A similar effect is observed when the initial guess orbital inclination angle parameter i is varied while keeping R₂/R₁ fixed (see the bottom row in Fig. 7). Changing the orbital inclination angle i to a larger or smaller value (e.g. 88° or 74° as indicated in the bottom left panel with the filled triangles and squares, respectively) leads to a smaller or larger area of the primary component being visible during its eclipse phases. This in turn leads to a larger or smaller flux contribution from the primary star in the initial guess, the effect that is compensated by a larger or smaller value of the R₂/R₁ parameter as optimised in the LSD-Binary module. As can be seen in Fig. 7 (bottom left panel), the effect is progressive with orbital phase and reaches its point of extremum at the centre of the primary eclipse. Just as is the case with the variable initial guess R₂/R₁ parameter, the orbital inclination angle i has strong effect on the shape of the initial guess LSD profile, which leads to the presence of a residual RM effect signal in the RV curve when the value of i is not optimal (see bottom right panel in Fig. 7).

In the last set of tests on artificial data, we investigate how the in-eclipse orbital-phase-resolved RV curve of the eclipsed primary component is affected by uncertainty in the determination of its T_eff and the components’ radii ratio R₂/R₁. To start with, we assume the R2/R1 parameter is known by other means, such as from a light curve solution of the system. In this particular scenario, we fix R₂/R₁ to its true value of 1.19 (see Table 1, column ‘RZ Cas’) and investigate how well we can recover the orbital inclination angle i of the system when the T_eff of the primary is varied within its typical spectroscopic uncertainties of 3–5%, or ~300 K in the absolute value for the system considered here. Since the RM effect is accounted for in the calculation of the initial guess LSD profiles and components’ flux ratio, the optimal value of the inclination angle is the one that provides minimal deviation of the RV curve obtained with the LSDBinary algorithm from the component’s RV curve that resembles its pure orbital motion. The top row in Fig. 8 summarises the obtained results by considering two extreme values of the effective temperature of the primary component: 8500 K (left panel) and 9100 K (right panel). The results are illustrated for three values of the orbital inclination angle i: two extreme cases of 74^° and 90°, and the optimum value of 82° we arrived at. We do not employ any minimisation algorithm to search for an optimum value of the inclination angle, but instead perform a grid search in the range i ∈ [74°,90°] with a step width of 1°. One can see that the true value of the inclination angle i = 82° is recovered well in both cases considered in this test, whereas other solutions are characterised by non-negligible deviations of individual RV measurements from those expected for the pure orbital motion of the star. These results reinforce our conclusions made earlier in Sect. 4.2 that any variations in atmospheric parameters that are within their typical spectroscopic uncertainties do not have a significant effect on the shape of the obtained LSD profiles and, thus, not on the RVs inferred from them.

Similarly, we now consider a case where the R₂/R₁ parameter is estimated with the LSDBinary algorithm from the out-of-eclipse spectra instead of it being fixed to an a priori known value. Considering the same two extreme values of T_eff of the primary, namely 8500 K and 9100 K, we compute initial guess LSD profiles and the components’ flux ratio for the out-of-eclipse phases employing the SynthV and Convolve suite of codes. Therefore, we assume symmetric LSD profiles for both binary components and a constant flux ratio, a well-justified assumption for the out-of-eclipse phases. We then let LSD-Binary optimise for the R₂/R₁ parameter by minimising the difference between the ‘observed’ and LSD-based model composite spectra. We find that R₂/R₁ = 1.02 and 1.36 provide the best match between the model and ‘observations’ when the T_eff of the primary is set to 8500 K and 91 000 K, respectively. The difference of some 15% between the best fit R₂/R₁ and its true value is the direct consequence of an incorrectly assumed T_eff of the primary component. A lower (higher) initial guess value of the primary’s effective temperature leads to underestimation (overestimation) of its flux contribution to the total light of the system. The inconsistency arising between the model and ‘observations’ is efficiently remediable by varying the R₂/R₁ parameter that, along with the effective temperatures, controls individual flux contributions of both binary components. Therefore, a lower (higher) initial guess T_eff of the primary component is compensated by its larger (smaller) radius and leads to a decrease (an increase) in the R₂/R₁ parameter.

With the R₂/R₁ parameter being estimated from the out-of-eclipse-spectra, we now proceed with the modelling of the in-eclipse spectra. The initial guess LSD profiles and the components’ flux ratio are computed with the RME module for incorrectly assumed values of R₂/R₁ and T_eff of the primary: (i) T_eff = 8500 K, R₂/R₁ = 1.02, and (ii) T_eff = 9100 K, R₂/R₁ = 1.36. The bottom row in Fig. 8 shows the results obtained for the above cases (i) (left panel) and (ii) (right panel) and forthree different values of the inclination angle i: 74° (open circles), 82° (red stars), and 90° (open squares). Although the RV curve corresponding to i = 82° still provides the smallest deviations from the RV curve resembling the pure orbital motion of the primary component compared to all other values of i considered in this test, the residual distortions due to the RM effect are significant. In our previous test with the variable initial guess T_eff of the primary component (see above and top row in Fig. 8), we demonstrated that the effect of this parameter on the in-eclipse RV curve is minimal, subject to the correct assumption of the orbital inclination angle. Therefore, we conclude that the R₂/R₁ parameter has a significantly larger effect on Level-0 and Level-1 data products of the LSDBinary algorithm than effective temperature of the star, despite both of these parameters being important for a correct prediction of the binary components’ flux ratio. This conclusion is also largely intuitive: whereas the T_eff of the primary significantly impacts its flux contribution to the total light of the system, the R₂/R₁ parameter also has a large effect on the geometry of the RM effect.

Fig. 7 Effect of the variable initial guess radii ratio R2/R1 (top row) and the orbital inclination angle i (bottom row) parameters on the orbital-phase-resolved radii ratio (left column) and radial velocity (right column) curves obtained with the LSDBinary algorithm. Filled squares and triangles in the top left panel refer to the initial guess R2/R1 = 1.1 and 1.4, respectively. The same symbols in the bottom left panel refer to the initial guess inclination angle i values of 74° and 88°, respectively. Red stars indicate a solution corresponding to the parameter values of R2/R1 = 1.2 and i = 82°; a combination that represents the closest match to the true parameter values listed in Table 1. See text for details.

5 Discussion and conclusions

In this study, we have presented a generalisation of the least-squares deconvolution (LSD) technique for eclipsing, spectroscopic double-lined binary systems. The main focus is on the separation of spectral contributions of individual binary components in the velocity space so that their RVs can be measured from the obtained LSD profiles with high precision. In our tests, we put a large emphasis on the in-eclipse phases where spectral contributions of the binary components are heavily blended. To that end, we develop a dedicated RME code that allows one to compute the initial guess LSD profiles and components’ flux ratio in the presence of the RM effect. At the same time, we include an option of initial guess symmetric LSD profiles and a constant with the orbital-phase components’ flux ratio to provide the means of working with the out-of-eclipse phases as well. This latter option can also be applied to the in-eclipse orbital phases, in which case an RV curve that is distorted due to the RM effect can be measured directly from the observations.

The LSDBinary software package developed in this work is freely available on GitHub¹ and comprises: (i) the SynthV and Convolve suite of codes to compute line masks and initial guess symmetric LSD profiles; (ii) the RME code to account for the RM effect in the calculation of the orbital-phase-dependent initial guess LSD profiles and binary components’ flux ratio; (iii) the LSDInit module that serves as a wrapper around the above-mentioned codes and whose main purpose is to organise input required for the core calculations; and (iv) the LSDBinary module that performs core calculations and, in its most basic configuration, delivers LSD profiles, LSD-based model spectra, and RVs of both binary components.

In this work, we also provide an extensive test of the LSDBinary algorithm based on artificial spectra of two binary systems: an Algol-like system with a main-sequence A-type primary and an evolved K-type giant secondary as components, and a detached binary with two F-type components. The main results and conclusions of our study can be summarised as follows:

• The algorithm has low sensitivity to the S/N of the input observed spectra, provided a sufficient number of atomic lines are present in the stellar spectrum. In practice, this makes the LSDBinary algorithm most suitable for stars of spectral type A and later. The regime of coolest stars whose spectra are dominated by molecular lines should be avoided as well, for the reasons of the current lack of the algorithm validation in the respective parameter space;

• Since resolving details of the line profile variations is an important aspect in the study of the in-eclipse spectra and the RM effect, the resolving power of the instrument should be sufficiently high. In this study, we investigate the regimes of low (R = 5000), medium (R = 25 000), and high (R = 65 000) resolution and find that the low resoltion is largely insufficient, irrespective of the configuration of the binary system and the amplitude of the RM effect. Depending on how large the RM effect is, a medium resolving power of the instrument of R ≈ 25 000 might suffice, although we recommend to push for a higher resolution R so that all details of the line profile variations can be resolved;

• We found the shapes of the obtained LSD profiles of both binary components to be largely insensitive to the variations in their atmospheric parameters T_eff, log g, and [M/H], provided these variations are within typical spectroscopic uncertainties. The largest effect is observed on the global scaling of the LSD profiles, which leaves the RVs inferred from them unaffected. On the other hand, variations in T_eff are found to have large effect on the components’ radii ratio inferred with the LSDBINARY algorithm. This result is largely intuitive given that an incorrect assumption for the T_eff of the star leads to overestimation or underestimation of its flux contribution to the total light of the system, the effect that is efficiently compensated in the LSDBinary algorithm by variations in the stellar radius through the radii ratio parameter R₂/R₁. We note, however, that LSD profiles, and hence RVs inferred from them, are expected to be largely affected by non-uniform distributions of physical quantities (e.g. temperature, velocity, and/or chemical elements) over the stellar surface. The LSDBinary code presented in this work includes a multi-profile LSD capability that, with just a mild development effort, could be exploited for the calculation of multi-temperature and/or multi-chemical composition LSD profiles²;

• The orbital inclination angle i and radii ratio R₂/R₁ are predictably the two parameters that have largest effect on the shape of the LSD profiles of the eclipsed binary component, and hence the shape and amplitude of the RM effect. Assuming the orbital configuration of the binary system and one of the above parameters are known, with the other one being optimised, the tandem of the RME and LSDBINARY modules successfully reconstruct the pure orbital motion-driven RV curve of the eclipsed star. Should this not be the case, it is an indication of either incorrectly determined atmospheric parameters of either or both binary components, or of deficiencies in stellar atmosphere models that fail to predict correctly the components’ flux ratio.

In practice, the latter finding can be exploited to either perform a consistency check between the photometric and spectroscopic solutions of a binary system, or to fine-tune the orbital inclination angle and/or the components’ radii ratio parameters from spectroscopic data. To demonstrate that, we return to the artificial in-eclipse data of an RZ Cas-like binary system whose configuration and atmospheric parameters are assumed to be known, except for the orbital inclination angle i that needs fine-tuning (see Table 1 for the list of the assumed parameters). We employ the tandem of the RME and LSDBinary codes to compute LSD profiles and RVs of both binary components from time series of the simulated in-eclipse spectra, while keeping the radii ratio parameter R₂/R₁ free in the LSDBINARY module. The results are illustrated in the top and bottom panels in Fig. 9 for RV and R₂/R₁ as a function of orbital phase, respectively. Both curves are shown for three values of the inclination angle i: 76° (open circles), 82° (red stars), and 88° (open squares). One can see that the best-fit solution is achieved at the true value of the orbital inclination angle and corresponds to the RV curve resembling the pure orbital motion of the star and the constant value of the radii ratio. Thus, simultaneously minimising deviations of the in-eclipse RV curve from its pure orbital motion version and of the R₂/R₁ curve from a constant value, allows one to fine-tune (or perform a spectroscopic consistency check for) one of the parameters considered in this test.

We note that the LSDBinary software package has been successfully applied to spectroscopic time series of Algol-type systems with pulsating primary components, R CMa (Lehmann et al. 2018) and RZ Cas (Lehmann et al. 2020). In both of these studies, the algorithm was used to compute the separated LSD profiles of both binary components and to measure precise orbital-phase-resolved RVs from them. The separated LSD profiles of the primary component of the R CMa system were additionally used to study line profile variations caused by non-radial pulsations of the star (Lehmann et al. 2018).

Fig. 8 Effect of incorrectly assumed binary parameters on the RV curve of the primary component. Top row: effect of incorrectly assumed initial guess Teff of the primary component, with R2/R1 being fixed to its true value of 1.19 (see Table 1). Left and right panel shows the effect of underestimated and overestimated effective temperature values, i.e. Teff = 8500 K and 9100 K, respectively. Red stars show the optimal solution obtained with LSDBinary corresponding to the orbital inclination angle i = 82°. Two more solutions corresponding to i = 74° (open circles) and 90° (open squares) are shown for comparison. Bottom row: combined effect of incorrectly assumed initial guess Teff of the primary and components’ radii ratio R2/R1. Left and right panel shows the cases of underestimated and overestimated parameters: Teff = 8500 K and R2/R1 = 1.02, and Teff = 9100 K and R2/R1 = 1.36, respectively. See text for details.

Fig. 9 Example of fine-tuning the inclination angle parameter i by minimising deviations of the R2/R1 (bottom) and RV (top) measurements obtained with the LSDBinary algorithm from their respective constant value and pure orbital motion-driven curve. The open circles, red stars, and open squares show the results obtained for the orbital inclination angle i = 76°, 82°, and 88°, respectively. The solid black lines shows the orbital motion-driven RV curve of the primary and the constant systemic value of R2/R1 = 1.19 in the top and bottom panels, respectively.

Acknowledgements

The research leading to these results has (partially) received funding from the KU Leuven Research Council (grant C16/18/005: PARADISE) and from the BELgian federal Science Policy Office (BELSPO) through PRODEX grant PLATO. Part of this work was supported by the German Deutsche Forschungsgemeinschaft, DFG project number Ts 17/2-1 and LE1102/3-1.

References

All Tables

All Figures

Fig. 1 LSDBinary algorithm flow chart. Bottom left panels show output of the LSDinit programme in terms of the LSD profiles computed from synthetic spectra (left) and the continuum flux ratio of the two binary components as a function of the wavelength (right). The output of the LSDBinary programme is illustrated in the two bottom right panels, showing LSD profiles of the individual binary components (left) and a comparison between the composite observed (black) and LSD-based model (red) spectra of the binary. See Sect. 3 for details.

In the text

Fig. 2 Graphical presentation of the LSDBinary algorithm. Left: initial guess synthetic LSD profiles computed with the LSDinit module, with (red stars) and without (black dots) the RM effect, along with the resulting ‘observed’ LSD profile (solid black line) computed with the LSDBinary module. Middle: in-eclipse phase-resolved RV curves, with (red stars) and without (black dots) the RM effect taken into account in the computation of the initial guess LSD profiles. The RVs measured relative to the respective initial guess LSD profiles are indicated with the arrows. Right: radii ratio as a function of orbital phase, with (red stars) and without (black dots) the RM effect taken into account in the computation of the initial guess synthetic spectrum-based LSD profiles. Arrows indicate the same data points as in the middle panel. See text for details.

In the text

Fig. 3 Results of the application of the LSDBinary algorithm to artificial in-eclipse spectra of the RZ Cas binary system. Top block: LSD profiles of the primary (red) and secondary (black) components. Bottom block: LSD-based model spectra (red) overlaid on the input ‘observed’ spectra (black) of the system. In each of the blocks, rows correspond to different values of the resolving power (from top to bottom, R = 60 000, 25 000, and 5000), while columns reflect the change in S/N of the spectrum (from left to right, S/N = 120, 80, and 40).

In the text

	Fig. 4 Same as for Fig. 3 but for the KIC 11285625 system.
In the text

Fig. 5 Effect of changing atmospheric parameters of the primary component of the RZ Cas system on the LSD profiles of the primary (left) and secondary (right). From top to bottom: Teff, log g, and [M/H] are varied between 8400 and 9200 K (step width of 100 K), 3.9 and 4.7 dex (step width of 0.1 dex), and −0.3 and 0.3 dex (step width of 0.1 dex), respectively. The LSD profiles shown in red correspond to the true parameters of the star listed in Table 1.

In the text

Fig. 6 Effect of changing atmospheric parameters of the secondary component of the RZ Cas system on the LSD profiles of the primary (left) and secondary (right). From top to bottom: Teff, log g, and [M/H] are varied between 4500 and 5100 K (step width of 100 K), 3.1 and 3.9 dex (step width of 0.1 dex), and −0.3 and 0.3 dex (step width of 0.1 dex), respectively. The LSD profiles shown in red correspond to the true parameters of the star listed in Table 1.

In the text

Fig. 7 Effect of the variable initial guess radii ratio R2/R1 (top row) and the orbital inclination angle i (bottom row) parameters on the orbital-phase-resolved radii ratio (left column) and radial velocity (right column) curves obtained with the LSDBinary algorithm. Filled squares and triangles in the top left panel refer to the initial guess R2/R1 = 1.1 and 1.4, respectively. The same symbols in the bottom left panel refer to the initial guess inclination angle i values of 74° and 88°, respectively. Red stars indicate a solution corresponding to the parameter values of R2/R1 = 1.2 and i = 82°; a combination that represents the closest match to the true parameter values listed in Table 1. See text for details.

In the text

Fig. 8 Effect of incorrectly assumed binary parameters on the RV curve of the primary component. Top row: effect of incorrectly assumed initial guess Teff of the primary component, with R2/R1 being fixed to its true value of 1.19 (see Table 1). Left and right panel shows the effect of underestimated and overestimated effective temperature values, i.e. Teff = 8500 K and 9100 K, respectively. Red stars show the optimal solution obtained with LSDBinary corresponding to the orbital inclination angle i = 82°. Two more solutions corresponding to i = 74° (open circles) and 90° (open squares) are shown for comparison. Bottom row: combined effect of incorrectly assumed initial guess Teff of the primary and components’ radii ratio R2/R1. Left and right panel shows the cases of underestimated and overestimated parameters: Teff = 8500 K and R2/R1 = 1.02, and Teff = 9100 K and R2/R1 = 1.36, respectively. See text for details.

In the text

Fig. 9 Example of fine-tuning the inclination angle parameter i by minimising deviations of the R2/R1 (bottom) and RV (top) measurements obtained with the LSDBinary algorithm from their respective constant value and pure orbital motion-driven curve. The open circles, red stars, and open squares show the results obtained for the orbital inclination angle i = 76°, 82°, and 88°, respectively. The solid black lines shows the orbital motion-driven RV curve of the primary and the constant systemic value of R2/R1 = 1.19 in the top and bottom panels, respectively.

In the text