Free Access
Issue
A&A
Volume 623, March 2019
Article Number A31
Number of page(s) 10
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/201834836
Published online 27 February 2019

© ESO 2019

1. Introduction

Spectral disentangling is a technique used to separate component spectra of spectroscopic binaries (SB2) and multiples. In comparison to separation techniques (e.g. Bagnuolo 1983) which need prior knowledge of radial velocities (RVs), it only requires a set of spectra spread over the orbital period at known times. This is made possible by connecting the disentangling algorithm with a global optimization algorithm, working on the orbital elements (or individual radial velocities). The solution is a disentangled spectrum of the two components (in the case of SB2), which can be interpreted as an extracted average of all individual single spectra from the observations. In addition, the orbital elements best fitting the data are derived. However, many binaries are variable and the influence of the changes in spectral behaviour on the output of the disentangling needs to be investigated and properly understood.

This paper is dedicated to a discussion of procedures dealing with different kinds of variability and how to extract the information of variability from the result of the disentanglement.

Spectangular was published by Sablowski & Weber (2017; hereafter Paper I); we point the reader to that article for the description of the program. The first application to highly variable Wolf-Rayet (WR) stars was presented in Shenar et al. (2017). Spectangular is based on singular value decomposition (SVD) and runs in a downhill simplex optimization. It differs significantly from another disentangling code, KOREL, which works in the frequency domain. Application of KOREL to variable spectra was presented in Harmanec et al. (2004) and Hadrava et al. (2009).

In Sect. 2 we summarize improvements to the code since it was first published. The SVD leads to results of least-squares nature, hence random variability will be strongly suppressed in the disentangled spectra. This is shown in Sect. 2.2, where we discuss the influence and a way to correct for telluric features by two-step disentangling. Section 3 concentrates on specific cases of variability of stellar origin. Section 3.2 shows a way to search for periodic line profile variations using the information in the spectrum-to-spectrum differences between the disentangled spectra and the individual observations. Continuum brightening has an effect on the extracted spectra as well, and is discussed in Sect. 3.3. Since eclipsing binaries play an important role in stellar physics, Sect. 3.4 is dedicated to demonstrating the applicability of the code to such systems. Surface spots imprint another periodic variability on the line profiles, and we investigate it in Sect. 3.5. The main results are summarized in Sect. 4.

2. Code modifications and procedures

For data preparation, we use the in-house developed code CroCo (see also Paper I). This code is used to re-sample the observations on a logarithmic wavelength scale and for two-dimensional cross-correlation (CC). The step size blog for the logarithmic wavelength scale was set to be the minimum step size, blog = min[blog, i], i = 0…n, found from all n spectra (multiplied by a user-set increment), where blog, i,  i = 0…n. is the minimum step size in spectrum i (which usually refers to the bluest pixel). Since that can lead to a step size of zero length (data errors), we implemented a user definable step size. Alternatively, the median, blog = median[blog, i], or the arithmetic mean, blog = mean[blog, i], of all these minimum step sizes can be used. Furthermore, the interpolation of the spectra can now be performed with a cubed spline (or a linear) interpolation. The spline can lead to better approximations at larger step sizes, corresponding to less data on the logarithmic scale.

When running the optimization with Spectangular on the individual RVs, there can be spectra for which the code does not find proper RVs. There can be three reasons for this: (i) the large number of variables, e.g. if 10 spectra of a SB2 system are used, we already have to deal with 20 variables; (ii) the noise level within the spectra; (iii) the selection of the initial values. In the case of non-convergence of optimization the disentangling fails. The SVD extracts a mean spectrum, best adapted to the data by the determined RVs. If the RVs are not realistic, the disentangled spectra are distorted. It is necessary now to identify these spectra and either remove them from the data set or change the initial RVs. This is usually best done during a first optimization run. The code writes some statistical information of the residual spectra (sum, mean, square sum) into a file (diffstatisic.dat) that can be used to identify the worst offenders. However, in most cases they can already be identified by visual inspection of the plotted differences in the GUI window. These differences are calculated by shifting the disentangled spectra by the determined RVs of each observation and subtracting these shifted disentangled spectra from the observations. As an example, we show the differences for such a case in Fig. 1 where we can clearly identify two spectra with outlying RVs. User interaction is likely to be necessary.

thumbnail Fig. 1.

Example of the differences between the disentangled spectra and 19 individual observations of Capella. The optimization has failed to find proper RV values for two spectra, which can be seen by the strong differences (black crosses). The spectra have been shifted by 0.3 on the y-axis for better visibility.

A further way to increase the convergence to the global minimum is to change the transformation coefficients of the optimizer (see Paper I). Barton & Ivey (1991), among others, found that a higher value of the contraction coefficient (e.g. 0.75 instead of 0.5) can be advantageous. This leads to a slower contraction of the simplex and an increase in the possibility to escape from local minima.

The program can also be run without the GUI, in a cluster environment for example. Only some minor changes need to be made in the source code, which are described in the manual of the software1.

2.1. Handling broad wavelength ranges

Since SVD is very time consuming, Spectangular is used in narrow wavelength ranges. In high-resolution spectroscopy, however, we have spectra with data points of the order of 100 000. Therefore, we have implemented a routine in CroCo to create whole data sequences. This routine creates several short re-sampled spectrum pieces written in separate files. From our experience with echelle spectra, it can be advantageous to constrain the spectrum pieces by the length of a spectral order. This prevents problems rising from order merging. The advantage of using echelle orders separately is also known from two-dimensional cross-correlation (see e.g. Zucker et al. 2003). Furthermore, for example in the case of hot stars, long pieces of pure continuum can be excluded since the disentangling relies on shifts of spectral lines.

In Spectangular, we implemented a function to disentangle such a sequence. To achieve this, the code needs a criterion to stop the optimization and to continue with the next wavelength range. The criterion for the optimization algorithm is the summed squares of differences, r, between the disentangled spectra and the observations (see Fig. 1). After each iteration of the optimization, the mean m and the standard deviation s of the simplex are calculated (see Paper I). At the end of an optimization run, the simplex has contracted to a small region around a single point. In that case, m (and s) do not change any more from one iteration to the next. This stage is defined here as stagnation of the optimization. When this stage is reached, the re-initiation function implemented will re-initiate the optimization. Hence, a complete optimization needs multiple re-initiations of the simplex. The end of the whole optimization is reached if the optimization stagnates and the minimal r corresponds to the initial parameter set. The auto-stop function implemented will stop the optimizer if that stage is reached. In the case where a sequence is specified, the code will load the next data set and start again with the optimization procedure.

2.2. Influence and correction of telluric lines

The tellurics can be removed by using the differences between the disentangled spectra and observations. An example is shown in Fig. 2 using the Hα region of Capella. The telluric absorption appears as emission in the residuals (black = observation - RV-shifted disentangled spectra). Here, we use a spline to reproduce these features and correct the individual observations by that spline. In the same figure the uncorrected spectrum (blue) and the spectrum corrected from terrestrial absorption (purple) are shown. To avoid removing any variability by fitting the differences, special care should be taken in the line profiles under consideration to make sure that only known terrestrial lines were fitted and removed. To identify these lines we used the HITRAN data base (Gordon et al. 2017). However, when the disentanglement is repeated on these cleaned spectra we did not observe any differences between the resulting line profiles. This is due to the strong variability of the terrestrial line strengths from spectrum to spectrum (in many spectra they are hardly visible) such that their influence is negligible.

thumbnail Fig. 2.

Example of the remaining differences (black at bottom) from the disentanglement around the Hα region used to remove the terrestrial absorption from the observation (blue at 1.0) to get a cleaned spectrum (purple, shifted to 1.25 for better visibility).

To support this statement, we show in Fig. 3 a zoomed-in image of an overplot of disentangled Hα profiles of the primary of Capella (G8III). The purple line is the result from the cleaned spectrum and the blue from the raw spectrum. The black line at the bottom shows the difference between these two profiles. Hence, such a telluric-check (TC) should always be performed to prevent any influence on the disentangled profile from tellurics.

thumbnail Fig. 3.

Zoomed-in image of the disentangled Hα profile of the primary of Capella when using spectra cleaned from terrestrial (purple) and non-cleaned (blue). Black is the difference between the two profiles.

If these lines do not show such a strong variability (also in the case of the strong terrestrial O2 absorption bands) the procedure will fail. In this case, a more classical approach is necessary (by observed tellurics within a spectrum of a hot star or simulated spectrum) to remove these lines before the disentangling is performed.

In addition, a third static component can be taken into account by the code. However, this can further complicate the optimization and will increase the time for computation significantly. Finally, for all of these procedures to work, it is necessary to use non-heliocentric corrected spectra. Otherwise, we have to deal with a third moving and variable component.

2.3. Signal-to-noise ratio and flux ratio

A major advantage of the disentanglement (especially in contrast to simple separation methods) is the enhancement of the S/N in the resulting spectra (see Table 1 in Paper I). However, for this to work the data has to be of equal quality in the sense of S/N and normalization. Assuming that this is fulfilled by the measurements, the S/N in the results can be estimated by simple statistics,

SNR A , B = SNR f A , B N , $$ \begin{aligned} \mathrm{SNR}_{A,B} = \mathrm{SNR} f_{A,B}\sqrt{N}, \end{aligned} $$(1)

Table 1.

Extracted period of pulsation for artificial SB2 data with one pulsating component.

where N is the number of observations,

f A , B = { k / ( 1 + k ) for component A 1 / ( 1 + k ) for component B , $$ \begin{aligned} f_{A,B} = {\left\{ \begin{array}{ll} k/(1+k)&\quad \mathrm{for\,component\,A } \\ 1/(1+k)&\quad \mathrm{for\,component\,B }, \\ \end{array}\right.} \end{aligned} $$(2)

the strengths of component A and B in normalized spectra (fA + fB ≡ 1), SNRA, B the estimated S/N of the disentangled spectra for component A and B, and k is the flux-ratio. Hence, to conserve the S/N of the observations, at least N=1/min[ f A,B 2 ] $ N = 1/\min [f_{A,B}^2] $ spectra are required. For example, if k = 1 the two components have identical strength, fB = fA = 0.5, and at least N = 4 spectra are required to maintain the S/N in the disentangled spectra.

3. Disentangling variable spectra

3.1. General remarks

As stated in the previous section, the disentanglement results in a spectrum of each component which can be interpreted as an extracted mean from the individual composite observations. Hence, there are two general ways to investigate variability. First, the differences between the disentangled spectra and the individual observations (O-D) are affected by variability such as non-periodic line profile variability and telluric features with large variations from observation to observation. This variability can be reconstructed within the O-D spectra, and in the case of telluric lines corrected for in the observations (see Sect. 2.2). This approach is applied to periodic line variability in Sect. 3.2. Second, if one component can be treated as stable compared to a more variable component, the disentangled spectrum of that stable component can be subtracted from the individual observations. This will result in a time series of individual spectra of the variable component and makes more detailed line-profile analysis possible. This approach is used in Sect. 3.5.

However, especially in active binaries or systems with complex temporal changes, both components may show variability. Hence, the line-profile analysis need to be performed in the rest frame of the star of interest. Therefore, computing line moments, the interval of integration needs to be shifted according to the orbital motion of the component. The subtraction algorithm implemented uses the final RVs to shift the disentangled spectrum to be subtracted and, according to user preference, can also shift the resulting time series by the RVs of the non-subtracted component. This results in a time series of spectra with no (or constant) shift with respect to the rest frame.

In that context, in some binaries where the orbital period differs significantly from the rotational period, the following procedure may be applied if rotational effects (e.g. surface spots) are in the focus. Let us consider a binary with components A and B. The orbital period is P and the periods of rotation for the primary and secondary are PA = P and PB ≪ PA, respectively. The spectra of the components can only be extracted if the orbital motion is well covered by observations at different orbital phases. For the analysis of the line profiles, however, it is advantageous to use spectra in quadrature (maximum separation of the two components in RVs). If we are interested in rotationally modulated effects of component B, two data sets are necessary, one set that covers the orbital period P (set 1, see Fig. 4a) and another set that covers the period of rotation of component B (set 2, see Fig. 4b), preferably around quadrature.

thumbnail Fig. 4.

Panel a: data set 1 (see text) spread over the whole orbital period to be used for disentangling. Panel b: data set 2 (see text) at orbital quadrature spread over the rotation period of the secondary used for rotational modulated line profile analysis.

For rotation modulated analysis of the line profile of the secondary, the spectra of the components are disentangled by the use of data set 1 (data spread over orbital period) first. The disentangled spectra can then be used as templates to perform a two-dimensional CC with the spectra of set 2 covering the rotational period of component B. The obtained RVs from the CC can be used to subtract the disentangled spectrum of component A from data set 2. This method is applied in Sect. 3.5 to artificial data to create surface temperature maps.

3.2. Radial pulsations

As an example of a binary with a pulsating component, we generated 20 artificial SB2 spectra of a non-variable sharp-lined component (vsini = 5 km s−1) and a periodic pulsating (change in FWHM of the lines) broad-lined (vsini = 20 km s−1) component with different amplitudes of the pulsation and a S/N of 50. The period of pulsation was one-sixth of the orbital phase and we used a Mizar-like orbit to generate the data (see Behr et al. 2011).

The resulting spectra of the disentangling for this artificial data is shown in Fig. 5. Purple is a composite spectrum at orbital phase 0.05, blue is the disentangled mean spectrum of the pulsating component, and green shows the disentangled spectrum of the non-pulsating component (we note the strong noise suppression in the disentangled spectra). To detect the variability, we subtract the disentangled spectra from the composite spectra and summed up the differences over wavelength for each spectrum. These differences versus orbital phase are plotted in Fig. 6a (crosses) for the data with 2 km s−1 pulsational amplitude. The green curve is a sinusoidal fit to the data. Extracted periods of the pulsation are summarized in Table 1. Columns 1 and 2 list the amplitude of the pulsation in km s−1 and in per cent of the rotational velocity of the pulsating star (20 km s−1), respectively. The extracted periods are inverse periods of the pulsation in fraction of the orbital period, that is, number of pulsation within the orbital period. As mentioned earlier, this value was set to six. The standard error results from a sin-fit to the data in least-squares minimization. As expected, the error increases with decreasing amplitude.

thumbnail Fig. 5.

Result from disentangling of the artificial SB2 with a broad-lined pulsating component. Blue: Disentangled mean spectrum of the pulsating component. Green: Same, but for the non-pulsating component. Purple: Artificial spectrum at orbital phase 0.05 for comparison.

thumbnail Fig. 6.

Summed differences (crosses) and sin-fit (green line) vs orbital phase for the artificial SB2 data. Panel a: with primary amplitude of pulsation of 2 km s−1 to extract the pulsational period. Panel b: with amplitude of pulsation for the primary and secondary of 2 km s−1 and 1 km s−1, respectively.

In the case of a similar variability of the two components, the differences can be used in an equivalent way. However, in any case, the signatures in these differences need to be tracked in RVs in order to identify the component to which the variability is connected. With a single pulsation period of both components, we show the differences and a combined fit in Fig. 6b. Again, we used six pulsations within the orbital phase for the primary; we used four for the secondary, instead of none as in previous test set-up. The extracted number of pulsations from the differences for the primary and secondary are 6.028 ± 0.022 and 3.799 ± 0.054, respectively.

3.3. Flares

In contrast to the variability by pulsations, flares are sudden non-periodic events. Hence, it is important to investigate how such changes in the spectra influence the extracted mean spectrum from the disentanglement. A flare also imprints changes in the radiated continuum of the star (e.g. Machado et al. 1980).

We generated a series of spectra (orbital elements, as used in Sect. 3.2) with one showing a flare as a brightening of the continuum of one component. This can be seen by a significant increase in the squared summed differences for this spectrum. Such an event introduces a slope to the continuum of the resulting spectra. This slope further depends on the spectrum on which the flare event happened. We show three examples in Fig. 7. The disentangled spectra at the bottom correspond to a series of 20 spectra and the flare happened in spectrum number 5 (phase 0.25). In the middle graph, the flare happened in spectrum number 9 (phase 0.50) and at the top in spectrum number 15 (phase 0.75). If in the second data set (flare at phase 0.50) we exchange the flare spectrum with that at phase 0.00, we obtain the same result. Therefore, the tilt does not depend on the order of the spectra within the data set, but within the orbital motion. This also shows the importance of equally well-normalized spectra. Hence, in the case of strongly slanted disentangled spectra, the normalization of the observations needs to be checked and corrected.

thumbnail Fig. 7.

Results from disentangled artificial spectra (time series of 20) showing a flare at phase 0.25 (bottom), 0.50 (centre), and 0.75 (top). The sign of the slope introduced in the results changes at the middle of such an equidistantly spread data set.

The resulting slope cannot be corrected by re-normalization of the disentangled spectra by division of the continuum. It is a difference of fluxes between the two components. Therefore, the slope has to be corrected by subtracting the difference. Let the continuum be c(λ) and the disentangled spectrum s(λ). The disentangled spectra are assumed to be calculated with the correct flux ratio. The corrected disentangled spectrum is then

s ( λ ) = s ( λ ) ( c ( λ ) 1 ) . $$ \begin{aligned} s^{\prime }(\lambda ) = s(\lambda )-(c(\lambda )-1). \end{aligned} $$(3)

To conserve the flux in the resulting disentangled spectra, the continuum function c(λ) is the same for both spectra,

s A , B ( λ ) = s ( λ ) ( c ( λ ) 1 ) , $$ \begin{aligned} s^{\prime }_{A,B}(\lambda ) = s(\lambda ) \mp (c(\lambda )-1), \end{aligned} $$(4)

such that what is subtracted in A has to be added to B.

As long as the difference in normalization is constant with wavelength, the slope is represented by a linear function, c(λ)=mλ + b. Hence, we have implemented a function to correct the resulting disentangled spectra from such a tilt.

3.4. Eclipsing binaries

The contribution of star A and B to the flux of the composite spectra of an eclipsing binary varies with phase. Hence, it is necessary for the disentanglement of such systems to know the light curve and a sufficient model.

3.4.1. Light curve model

In this section a simple spherical geometry model (no eccentricity and limb darkening) to generate a light curve is used to demonstrate how the disentanglement behaves with the data.

Let r be the vector distance between the two stars. From classical mechanics (see e.g. Goldstein et al. 2000) we know that the distance to the centre of mass is

r A = r M B / ( M A + M B ) = r β A , r B = r M A / ( M A + M B ) = r β B . $$ \begin{aligned} r_{A}&= - r~M_{B}/(M_{A}+M_{B})=-r~ \beta _A, \nonumber \\ r_{B}&= r~M_{A}/(M_{A}+M_{B}) = r~ \beta _B. \end{aligned} $$(5)

The geometry of the situation is sketched in Fig. 8a. From that figure we find

x = r sin ( ϕ ) , y = r cos ( i ) cos ( ϕ ) , and z = r sin ( i ) cos ( ϕ ) , $$ \begin{aligned}&x = r \sin (\phi ), ~y = r \cos (i)\cos (\phi ),\nonumber \\&{\mathrm{and} }~z = r \sin (i)\cos (\phi ), \end{aligned} $$(6)

thumbnail Fig. 8.

Panel a: geometry of the binary model with spherical stars and orbit. The z-axis points towards the observer. Panel b: circle cross section of the stars during eclipse.

where ϕ is the phase angle and i the inclination. Substituting this into Eq. (5) yields the coordinates in the centre-of-mass system

x A = r sin ( ϕ ) β A , y A = r cos ( i ) cos ( ϕ ) β A and z A = r sin ( i ) cos ( ϕ ) β A , $$ \begin{aligned}&x_A = -r \sin (\phi )\beta _A,~y_A = -r\cos (i)\cos (\phi )\beta _A \nonumber \\&{\mathrm{and} }~z_A = -r \sin (i)\cos (\phi )\beta _A, \end{aligned} $$(7)

for component A and

x B = r sin ( ϕ ) β B , y B = r cos ( i ) cos ( ϕ ) β B , and z B = r sin ( i ) cos ( ϕ ) β B , $$ \begin{aligned}&x_B = r \sin (\phi )\beta _B,~y_B = r\cos (i)\cos (\phi )\beta _B, \nonumber \\&{\mathrm{and} }~z_B =-r \sin (i)\cos (\phi )\beta _B, \end{aligned} $$(8)

for component B. In addition, let be RA and RB the radii, AA and AB the unobscured areas, and LA and LB the luminosities of the stars. Hence, we can write

I ( ϕ ) = n ( L A A A / R A 2 + L B A B / R B 2 ) / ( 4 π ) = I A ( ϕ ) + I B ( ϕ ) , $$ \begin{aligned} \begin{aligned} I(\phi )&= n (L_A A_A / R^2_A + L_B A_B /R^2_B)/(4\pi )\\&=I_A(\phi ) + I_B(\phi ), \end{aligned} \end{aligned} $$(9)

for the brightness of the system, where n is a constant to normalize the light curve. The projected distance between the two stars in the plane perpendicular to the line of sight is

d AB = ( x B x A ) 2 + ( y B y A ) 2 . $$ \begin{aligned} d_{AB} = \sqrt{(x_B-x_A)^2 + (y_B - y_A)^2}. \end{aligned} $$(10)

During an eclipse, the common cross section of the stars is described by circular segments of area d A A,B = R A,B 2 ( α A,B sin α A,B )/2 $ d{A_{A,B}} = R_{A,B}^2({\alpha _{A,B}} - \sin {\alpha _{A,B}})/2 $. The angle α, the segment extents (see Fig. 8b), is found by the projected distance in the plane perpendicular to the line of sight by α A,B =2arccos[(( R B 2 R A 2 )+ d AB 2 )/(2 R A,B d AB )] $ {\alpha _{A,B}} = 2\arccos [( \mp (R_B^2 - R_A^2) + d_{AB}^2)/(2{R_{A,B}}{d_{AB}})] $. For example, in the case that zA >  zB star A is in front of star B. Hence the areas are A A =π R A 2 $ {A_A} = \pi R_A^2 $ and A B =π R B 2 (d A A +d A B ) $ {A_B} = \pi R_B^2 - (d{A_A} + d{A_B}) $.

An example light curve of this simple model is shown in Fig. 9.

thumbnail Fig. 9.

Example of a light curve modelled by the equations presented here. I is the (measurable) combined brightness of IA and IB and Ratio is the brightness (or light) ratio between the two components, IB/IA.

3.4.2. Artificial eclipsing binary

The values used for the model are summarized in Table 2. This model was used to create a set of 20 spectra equidistantly spread over the orbital phase (orbital parameters as in Sect. 3.2). From an observational point, the disentangling of such a set of data would be possible by using the spectra at phase 0.5 since the system shows a total eclipse; that is, only the signal of the primary is present in that spectrum. In wide eclipsing binaries, the disentangling can be performed with out-of-eclipse spectra. However, in close binaries the eclipses can cover a significant fraction of the orbital period (see Fig. 9). Furthermore, it is necessary for the disentanglement to use spectra that cover the orbital motion well. In the case of close binaries, spectra taken only out-of-eclipse may be insufficient for the disentangling. Hence, it is necessary to also use the spectra obtained during eclipse.

Table 2.

Model values for the light curve shown in Fig. 9.

The light ratio necessary for the disentangling can be computed by

k ( ϕ ) = 1 / f B ( ϕ ) 1 = f A ( ϕ ) / [ 1 f A ( ϕ ) ] , $$ \begin{aligned} k(\phi ) = 1/f_B(\phi )-1=f_A(\phi )/[1-f_A(\phi )], \end{aligned} $$(11)

where

f A ( ϕ ) = I A ( ϕ ) / I ( ϕ ) and f B ( ϕ ) = I B ( ϕ ) / I ( ϕ ) . $$ \begin{aligned} f_A(\phi )=I_A(\phi )/I(\phi )~{\mathrm{and} }~f_B(\phi )=I_B(\phi )/I(\phi ). \end{aligned} $$(12)

From Sect. 3.3, we expect to see a slope in the continuum of the disentangled spectra if the phase dependent flux ratios are not taken into account. The result, where a constant ratio of k(ϕ)=1 is applied to this data, is shown in Fig. 10a. The appearance of such a slope gives a hint either to a variable flux of the system or to incorrect normalization. Therefore, if the case of bad normalization can be excluded, additional data (photometry) is necessary for a proper disentangling of the system. On the other hand, if we can exclude variability of the system, the normalization of the data needs to be corrected. In the same figure, we can also see that the line strengths are wrong (negative values for component A). It is therefore obvious that using the correct flux ratio (also if constant) is of major importance for any further analysis of the disentangled spectra.

thumbnail Fig. 10.

Disentangled spectra of the artificial eclipsing binary. Panel a: phase dependent flux ratio has not been taken into account. As shown in the previous section, this introduces a slope to the disentangled spectra. Shown are the sharp-lined component A (purple) and B (green). Panel b: disentangled spectra with correct flux ratios used in comparison with the templates used to create the set of spectra. Template A and B are shifted by 0.1 and −0.1 on the y-axis for visibility. Panel c: as panel b, but with flux ratios from the optimization with start ratio of k = 1 for all phases.

The disentangled spectra with the correct flux ratios are shown in Fig. 10b together with the template spectra for both stars. They are shifted by 0.1 for A and −0.1 for B in the vertical direction for better visibility. In addition to the excellent agreement between disentangled spectra and templates, we can again see the enhancement of the signal-to-noise ratio.

3.4.3. Extracting the flux ratios

As mentioned in Sect. 3.4.2, knowing the flux ratios is essential for disentanglement in eclipsing systems. Since the extracted spectra is the mean spectrum present in all observations, the differences can also be used to optimize the flux ratio k(ϕ) in each spectrum. Therefore, we implemented a downhill simplex algorithm to find these flux ratios. This algorithm provides two targets for minimization: the summed squares of differences (as used for the spectral disentanglement) and the summed square of the slopes, m A 2 + m B 2 $ \sqrt{m^2_A+m^2_B} $, of linear functions, cA, B = mA, Bλ + bA, B, representing the continuum of each disentangled spectrum (see also Sect. 3.3). The latter turned out, in our test runs, to be too insensitive. However, the minimization of the residuals with an initial constant k(ϕ)=1 lead to flux ratios plotted in Fig. 11 together with the values from the model, as described in Sect. 3.4.1. Furthermore, the resulting spectra in comparison with the templates used to create the data set are shown in Fig. 10c.

thumbnail Fig. 11.

Comparison of the component strengths fA and fB from the model used (Model A and B) and after optimization (Optimized A and B) on flux ratio k(ϕ), Eq. (11), of the two components A and B. Also plotted are the differences of Model – Optimized.

Knowing the correct RVs was assumed for that test run. They can be obtained, in a first step, by a CC. An iterative procedure alternating the optimization on flux ratios and RVs (or orbital elements) could be applied to such data with no information from photometry.

Knowing the flux ratio is not sufficient to get the light curve. However, if I(ϕ) is measured we can obtain the individual brightnesses of the two components from Eq. (12). Therefore, the combination with photometry yields more information about the system, i.e. a separate light curve for both components.

3.5. Surface spots

Spots appear in the line-profiles as moving “bumps” following the rotation of the star. Therefore, we constructed an artificial spotted star (A; see Fig. 15a) and calculated the spectra from this configuration (the same orbital parameters as used in Sect. 3.2) with two spot features at almost opposite sides of the stellar surface, and put that star in orbit with an unspotted component (B). The resulting line-profiles are illustrated in Fig. 12a. Then the spectra are disentangled and the component (B) without spots is subtracted from the composite spectra using the RVs from a CC with the disentangled spectra as templates (see Fig. 12b). The resulting spectra of component A are corrected for RV to have the line-profile at the same wavelength, which ideally would be the rest wavelength. Especially when subtracting a sharp-lined component, small mismatches between the disentangled and the composite cause artefacts in the resulting “cleaned” spectrum. These artefacts are visible in Fig. 12b and follow the relative motion in wavelength between the two components. We show in Fig. 13a a zoom-in of a small spectral region where we subtracted the sharp-lined component and plot over the input spectrum at phase 0.05. In Fig. 13b is a further zoom-in of a single line, but for the whole spectral sequence to visualize the line-profile asymmetries caused by the spots.

thumbnail Fig. 12.

Panel a: artificial data set of SB2 system with a spot on the broad-lined component. Disentangled profiles renormalized to one are shown at the top in blue and purple. Panel b: disentangled spectrum of the sharp-lined component is subtracted from the composite spectra and the result is shifted to the rest wavelength. These spectra were used to create the Doppler maps shown in Fig. 15c.

thumbnail Fig. 13.

Panel a: zoom-in of three lines of the spectra shown in Fig. 12b. Purple is the input spectra of the spotted component, black the spectra of the spotted component after subtraction of the disentangled non-spotted component used to calculate the map in Fig. 15c. Panel b: zoom-in of the leftmost line from panel a, but for all the spectra used to construct the maps.

As shown in Sect. 3.2, the residual spectra can be used to analyse variability and to search for periods. Therefore, we plot the summed differences in Fig. 14 together with a fit accounting for two spots crossing the surface. Hence, these differences can be used to search for the rotation period of the stars.

thumbnail Fig. 14.

Summed differences of the residuals from the artificial spot data (crosses). The solid line is a fit for two spots.

The stellar surface of the test star is reconstructed using the iMap code (for details, see Carroll et al. 2007, 2009, 2012). We note here that the two-temperature approximation used in the input map (Fig. 15a) is not ideal for inversions as shown in Fig. 15b, but is sufficient to investigate the differences between ideal spectra and disentangled spectra.

In the resulting maps shown in Fig. 15, we can see a small suppression of the temperature contrast. In the case of the ideal data (Fig. 15b), the coolest temperature recovered is 5121 K when the input value is 5000 K, whereas the disentangled spectra lead to the coolest temperature of 5210 K, see Fig. 15c. The shift of the spectra to the same wavelength after subtracting one component is crucial to recover the temperatures. This is demonstrated by the map in Fig. 15d. This map was computed from spectra shifted by RVs obtained from a CC without fitting the maximum of the CC function, but by manually adjusting the spectra. Therefore, this data set shows residual RV shifts between the spectra of the order of 0.05 pix (≈ 100 m s−1), and the obtained coolest temperature (5257 K) is a bit higher than in the case of ideal RV shifts. We also note the presence of a very weak spot at phase 0.5 in these maps. We attribute the remaining difference in temperature to residuals of the subtraction and re-sampling effects.

thumbnail Fig. 15.

Panel a: input map used to create artificial spot spectra shown from four different viewing angles. Panel b: surface map using the artificial spot spectra without a second component. Panel c: resulting surface map from the artificial spot data after subtraction of disentangled spectra. Panel d: as panel c, but with RVs to align the spectra obtained by CC and manually adjusted.

Regardless of the small discrepancies in temperatures, the resulting Doppler maps in Fig. 15 display the two spots at correct positions. Hence, a proper analysis of a systematic suppression of the temperature difference of spot and stellar surface can be performed for an individual target.

3.6. Impact of variability on disentangling results

Since the actual errors on the orbital parameters determined by the optimization algorithm depend, in first order, on the data quality and the amplitude of the variability, only an idea of their magnitudes can be given here. In addition, in real-world applications we most likely have to deal with a combination of multiple effects.

In the case of eclipsing binaries, it was suggested in Sect. 3.4.3 to iteratively alter the optimization on flux ratios and orbital elements. To demonstrate this procedure, the data from Sect. 3.4.2 were used with an initial assumption of a constant flux ratio of k(ϕ)=1 and the orbital parameters as listed in Table 3 second column. The orbital elements used to create the data are listed in the first column. In runs 2 and 4 the flux ratios were optimized (no change in the orbital elements, therefore not listed). Since the period is often well determined, we fixed the period, P, and the time of periastron passage, T0. In addition, for this example we fixed the systemic velocity γ as well.

Table 3.

Orbital elements for the eclipsing binary after iterative optimization on orbital elements and flux ratios.

From the results, we can see that the amplitudes KA and KB show the strongest deviation. These amplitudes may be extremely important when it comes to the determination of masses. However, the overall solution after four runs is close to the ideal values especially if we consider the resolution of 1.7 km s−1 of the spectra. Another optimization run on the orbital elements did not improved the results further.

Nevertheless, as was shown in Fig. 10a, incorrect flux ratios have a strong influence on the extracted spectra in the sense of line-depth. Hence, if spectral analyses are to be performed on the disentangled spectra, special care needs to be taken on the flux ratios.

The final precision of the flux ratios should be limited primarily by the S/N of the spectra. The rms (between the ideal flux ratio used to create the data and the value obtained from optimization) from the 20 flux ratios obtained in Sect. 3.4.3 and obtained here from run 4 are 0.0047 and 0.0046, respectively. From the 20 spectra with a S/N of 50, we get 1 / ( 50 2 0 ) = 0.0045 $ 1/(50\sqrt20)=0.0045 $, which is in agreement with the rms values of the obtained flux ratios.

From this result we expect that a single flare spectrum has a negligible effect on the orbital parameters as long as the number of spectra without flares is much higher (> 10). Figure 16 shows a comparison of the line-profiles of the disentangled sprectrum with the template spectrum (as in Sect. 3.3). Noise was added to the template to construct the data. The continuum was a factor of three brighter during the flare. However, the disentangled spectra agrees very well with the template.

thumbnail Fig. 16.

Comparison of line-profiles between disentangled and template spectra using the same data as in Sect. 3.3.

In Fig. 17 the equivalent width (EW) of the line at 6721.73 Å from Fig. 13 is shown for the data used in Sect. 3.5. The ideal data (green x) shows the periodicity as reflected by the summed differences shown in Fig. 14. Therefore, both measures could be used to search for the stellar rotation periods. Purple crosses are used for the EW of the line after the disentangled spectrum of the constant component was subtracted (Fig. 13b) and the blue square is the data point from the disentangled spectrum of the spotted component. The subtraction leads to scatter of the data points of the EW measured.

thumbnail Fig. 17.

Equivalent width (EW) [Å] of the line at 6721.73 Å from the artificial spot data of Sect. 3.5. Green is for the ideal single-component spectra used to create the composite binary spectra, purple (cleaned) is after subtraction of the constant component spectrum from the disentanglement, and the blue square corresponds to the EW of that line of the disentangled spectrum.

4. Summary

In this work we have investigated the influence of variability on the result of the disentanglement. This included tellurics and their correction, pulsations (change in FWHM of the lines), flare as continuum brightening, eclipsing binaries, and line profile variation by surface spots. The correction of tellurics is possible as long as they are highly variable in strength. In this case, the telluric features are strongly suppressed in the disentangled spectra. Furthermore, these lines can be reconstructed for each observation from the residuals by a spline approximation. In the case of strong variability from spectrum to spectrum, these tellurics do not influence the disentangled spectra.

Two ways of treating the variability were presented. First, the analysis of the residuals between the disentangled spectra and the individual observations can be used to search for periodic changes. Second, subtracting the disentangled spectra of component B, for example, from the observations leads to a time series of spectra for component A which can be used for further analysis (e.g. Doppler imaging, line-profile analysis).

We note that the code cannot account for wavelength dependent flux ratios. Therefore, temperature effects introducing a slope-change of the continuum may further affect the disentangling. However, since the wavelength ranges are rather narrow (in high-resolution spectra) the slope in such a small region may have a negligible effect. Since it also depends on the temperature difference, object-specific simulations should be performed. A wavelength dependent flux ratio to account for such effects (also from irradiation/heating) is a subject of future improvements to the code.

We also presented the disentanglement on artificial data of an eclipsing binary. The code allows the user to specify the contribution to the flux of each spectrum and each component. As long as the correct flux ratios are used, the disentangled spectra are in very good agreement with the templates. Furthermore, it is now possible to run an optimization on the flux ratios of each observed spectrum. This now provides a method for obtaining information about the flux variation. However, we assumed that the radial velocities are known, for example by an initial cross-correlation. Hence, the disentangling approach is also usable to extract light curve information from the spectroscopic data. A combination of light curve modelling, spectral disentangling, and solving for the flux ratios is beyond of the scope of this paper.

Additionally, since the flux-ratios are not bound to a light curve model, the variation in the brightness of the system from other effects (e.g. non-symmetry) can be accounted for. In this context, we note that a simultaneous solution from light- and RV-curves can help to further improve the results, for example as was implemented in FOTEL (see Hadrava 1990). This coupling is beyond of the scope of this work and is the subject of future improvements of the code. Effects like apsidal motion, non-constant orbital period, and a third spectroscopically non-visible component can be taken into account by optimizing on the individual RVs.

Doppler imaging showed a sensitivity of the temperature difference between spot and surface to the alignment in wavelength of the time series of the spectra. This was shown by a comparison between maps obtained with spectra shifted by RVs from a CC and a manually fine adjustment compared to a map with the RVs computed to create the binary spectra. In a science application a fit of the CC maximum can be performed, and an iterative approach using line profiles derived from the Doppler maps as templates can be used to improve the accuracy of the RVs.

It is important to note that a proper normalization of the observations needs to be done with special care. This can introduce a slope to the resulting spectra. This is also true for a brightening of the continuum (e.g. by flares). Therefore, we finally note that a proper preparation of the data for the disentanglement is of great importance since it can mimic variability.


Acknowledgments

We thank the anonymous referee for the constructive comments that helped to improve the content and readability of this work. Based in part on data obtained with the STELLA robotic telescope in Tenerife, an AIP facility jointly operated by AIP and IAC.

References

  1. Bagnuolo, Jr., W. G. 1983, Lowell Obs. Bull., 9, 180 [NASA ADS] [Google Scholar]
  2. Barton, R. R., & Ivey, J. S. 1991, Winter Simulation Conference Proceedings, 945 [Google Scholar]
  3. Behr, B. B., Cenko, A. T., Hajian, A. R., et al. 2011, AJ, 142, 6 [NASA ADS] [CrossRef] [Google Scholar]
  4. Carroll, T. A., Kopf, M., Ilyin, I., & Strassmeier, K. G. 2007, Astron. Nachr., 328, 1043 [NASA ADS] [CrossRef] [Google Scholar]
  5. Carroll, T. A., Kopf, M., Strassmeier, K. G., & Ilyin, I. 2009, in IAU Symposium, eds. K. G. Strassmeier, A. G. Kosovichev, & J. E. Beckman, 259, 633 [NASA ADS] [Google Scholar]
  6. Carroll, T. A., Strassmeier, K. G., Rice, J. B., & Künstler, A. 2012, A&A, 548, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Goldstein, H., Poole, C., & Safko, J. 2000, Classical Mechanics, 3rd edn. (Addison Wesley), 71 [Google Scholar]
  8. Gordon, I., Rothman, L., Hill, C., et al. 2017, J. Quant. Spectr. Rad.Transf., 203, 3 hITRAN2016 Special Issue [NASA ADS] [CrossRef] [Google Scholar]
  9. Hadrava, P. 1990, Contrib. Astron. Obs. Skalnate Pleso, 20, 23 [Google Scholar]
  10. Hadrava, P., Slechta, M., & Skoda, P. 2009, A&A, 507, 397 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Harmanec, P., Uytterhoeven, K., & Aerts, C. 2004, A&A, 422, 1013 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Machado, M. E., Avrett, E. H., Vernazza, J. E., & Noyes, R. W. 1980, ApJ, 242, 336 [NASA ADS] [CrossRef] [Google Scholar]
  13. Sablowski, D. P., & Weber, M. 2017, A&A, 597, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Shenar, T., Richardson, N. D., Sablowski, D. P., et al. 2017, A&A, 598, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Zucker, S., Mazeh, T., Santos, N. C., Udry, S., & Mayor, M. 2003, A&A, 404, 775 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

All Tables

Table 1.

Extracted period of pulsation for artificial SB2 data with one pulsating component.

Table 2.

Model values for the light curve shown in Fig. 9.

Table 3.

Orbital elements for the eclipsing binary after iterative optimization on orbital elements and flux ratios.

All Figures

thumbnail Fig. 1.

Example of the differences between the disentangled spectra and 19 individual observations of Capella. The optimization has failed to find proper RV values for two spectra, which can be seen by the strong differences (black crosses). The spectra have been shifted by 0.3 on the y-axis for better visibility.

In the text
thumbnail Fig. 2.

Example of the remaining differences (black at bottom) from the disentanglement around the Hα region used to remove the terrestrial absorption from the observation (blue at 1.0) to get a cleaned spectrum (purple, shifted to 1.25 for better visibility).

In the text
thumbnail Fig. 3.

Zoomed-in image of the disentangled Hα profile of the primary of Capella when using spectra cleaned from terrestrial (purple) and non-cleaned (blue). Black is the difference between the two profiles.

In the text
thumbnail Fig. 4.

Panel a: data set 1 (see text) spread over the whole orbital period to be used for disentangling. Panel b: data set 2 (see text) at orbital quadrature spread over the rotation period of the secondary used for rotational modulated line profile analysis.

In the text
thumbnail Fig. 5.

Result from disentangling of the artificial SB2 with a broad-lined pulsating component. Blue: Disentangled mean spectrum of the pulsating component. Green: Same, but for the non-pulsating component. Purple: Artificial spectrum at orbital phase 0.05 for comparison.

In the text
thumbnail Fig. 6.

Summed differences (crosses) and sin-fit (green line) vs orbital phase for the artificial SB2 data. Panel a: with primary amplitude of pulsation of 2 km s−1 to extract the pulsational period. Panel b: with amplitude of pulsation for the primary and secondary of 2 km s−1 and 1 km s−1, respectively.

In the text
thumbnail Fig. 7.

Results from disentangled artificial spectra (time series of 20) showing a flare at phase 0.25 (bottom), 0.50 (centre), and 0.75 (top). The sign of the slope introduced in the results changes at the middle of such an equidistantly spread data set.

In the text
thumbnail Fig. 8.

Panel a: geometry of the binary model with spherical stars and orbit. The z-axis points towards the observer. Panel b: circle cross section of the stars during eclipse.

In the text
thumbnail Fig. 9.

Example of a light curve modelled by the equations presented here. I is the (measurable) combined brightness of IA and IB and Ratio is the brightness (or light) ratio between the two components, IB/IA.

In the text
thumbnail Fig. 10.

Disentangled spectra of the artificial eclipsing binary. Panel a: phase dependent flux ratio has not been taken into account. As shown in the previous section, this introduces a slope to the disentangled spectra. Shown are the sharp-lined component A (purple) and B (green). Panel b: disentangled spectra with correct flux ratios used in comparison with the templates used to create the set of spectra. Template A and B are shifted by 0.1 and −0.1 on the y-axis for visibility. Panel c: as panel b, but with flux ratios from the optimization with start ratio of k = 1 for all phases.

In the text
thumbnail Fig. 11.

Comparison of the component strengths fA and fB from the model used (Model A and B) and after optimization (Optimized A and B) on flux ratio k(ϕ), Eq. (11), of the two components A and B. Also plotted are the differences of Model – Optimized.

In the text
thumbnail Fig. 12.

Panel a: artificial data set of SB2 system with a spot on the broad-lined component. Disentangled profiles renormalized to one are shown at the top in blue and purple. Panel b: disentangled spectrum of the sharp-lined component is subtracted from the composite spectra and the result is shifted to the rest wavelength. These spectra were used to create the Doppler maps shown in Fig. 15c.

In the text
thumbnail Fig. 13.

Panel a: zoom-in of three lines of the spectra shown in Fig. 12b. Purple is the input spectra of the spotted component, black the spectra of the spotted component after subtraction of the disentangled non-spotted component used to calculate the map in Fig. 15c. Panel b: zoom-in of the leftmost line from panel a, but for all the spectra used to construct the maps.

In the text
thumbnail Fig. 14.

Summed differences of the residuals from the artificial spot data (crosses). The solid line is a fit for two spots.

In the text
thumbnail Fig. 15.

Panel a: input map used to create artificial spot spectra shown from four different viewing angles. Panel b: surface map using the artificial spot spectra without a second component. Panel c: resulting surface map from the artificial spot data after subtraction of disentangled spectra. Panel d: as panel c, but with RVs to align the spectra obtained by CC and manually adjusted.

In the text
thumbnail Fig. 16.

Comparison of line-profiles between disentangled and template spectra using the same data as in Sect. 3.3.

In the text
thumbnail Fig. 17.

Equivalent width (EW) [Å] of the line at 6721.73 Å from the artificial spot data of Sect. 3.5. Green is for the ideal single-component spectra used to create the composite binary spectra, purple (cleaned) is after subtraction of the constant component spectrum from the disentanglement, and the blue square corresponds to the EW of that line of the disentangled spectrum.

In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.