Using Independent Component Analysis to Detect Exoplanet Reflection Spectrum from Composite Spectra of Exoplanetary Binary Systems

The model adopted in this paper to simulate the photometric variations of the stellar light reflected from a planet orbiting its host star is based on the discussion presented in Charbonneau et al. (1999) and reused by Martins et al. (2013). The phase-dependent flux ratio of the planet flux F_planet relatively to the star flux F_star is given by

$$\begin{eqnarray}&&\displaystyle \frac{{F}_{\mathrm{planet}}(\alpha )}{{F}_{\mathrm{star}}}={pg}(\alpha ){\left(\displaystyle \frac{{R}_{{\rm{p}}}}{a}\right)}^{2}\end{eqnarray} \tag{ 1 }$$

where p is the geometric albedo of the planet, R_p is the planet radius, a is the planet orbital distance and g(α) is the phase function. It holds for the case R_p ≪ R_star ≪ a. As described in Borra & Deschatelets (2018), the phase angle, α, is the angle between the star and the Earth when seen from the planet. It depends on the position of the planet in its orbit and takes a value in the range 0 ≤ α ≤ 180°. For an orbital inclination i and orbital phase ϕ, assuming a circular orbit, α can be derived from

$$\begin{eqnarray}&&\cos (\alpha )=-\sin (i)\cos (2\pi \phi ).\end{eqnarray} \tag{ 2 }$$

The orbital phase ϕ in Equation (2) traces the position of the planet on its orbit around the star. ϕ takes values between 0 and 1. In Charbonneau et al. (1999) it is assumed 0 at the time of maximum radial velocity, in other works (e.g., Borra & Deschatelets 2018), phase 0 corresponds to the time of the mid-point of the transit. In this paper we will use the latter definition, i.e., for us ϕ = 0 at the mid-point of the transit or of the inferior conjunction. For this particular choice, at a given time t, ϕ value can be calculated knowing the planet time of transit t₀ and its orbital period P_orb by

$$\begin{eqnarray}&&\phi =\displaystyle \frac{t-{t}_{0}}{{P}_{\mathrm{orb}}}.\end{eqnarray} \tag{ 3 }$$

The phase function g(α) in Equation (1) models the fraction of the maximum planet's reflected flux at full phase. It varies, depending on the planet position on its orbit, since we see a different portion of the planet illuminated hemisphere. As described in Langford et al. (2011), in the simplest case, we can assume an isotropic scattering over 2π sr (Lambert-law sphere) in which case g(α) has typically the form

$$\begin{eqnarray}&&g(\alpha )=\displaystyle \frac{[\sin (\alpha )+(\pi -\alpha )\cos (\alpha )]}{\pi }\end{eqnarray} \tag{ 4 }$$

and α is related to orbital phase ϕ as described by Equation (2).

The geometric albedo p, which appears in Equation (1), is the reflectivity of a planet measured at superior conjunction. It is a wavelength-dependent, dimensionless number, with values between 0 and 1 obtained as a ratio between the reflected and incident flux.

The radial velocity equations describing the system are

$$\begin{eqnarray}&&{\mathrm{RV}}_{\mathrm{star}}=\gamma +{K}_{\mathrm{star}}[\cos (\omega +f)+e\cos (\omega )]\end{eqnarray} \tag{ 5 }$$

and

$$\begin{eqnarray}&&{\mathrm{RV}}_{\mathrm{planet}}=\gamma -{K}_{\mathrm{planet}}[\cos (\omega +f)+e\cos (\omega )]\end{eqnarray} \tag{ 6 }$$

for the star and the planet respectively. γ corresponds to the system's barycenter radial velocity relatively to the Sun, ω is the argument of periastron, e is the eccentricity, and f is the true anomaly; K_star and K_planet are the radial velocity semi-amplitude of, respectively, the stellar and planetary orbits and can be computed from the orbital parameters as

$$\begin{eqnarray}&&{K}_{\mathrm{star}}=\displaystyle \frac{2\pi a}{{P}_{\mathrm{orb}}}\displaystyle \frac{{m}_{\mathrm{planet}}}{{m}_{\mathrm{star}}+{m}_{\mathrm{planet}}}\displaystyle \frac{\sin i}{\sqrt{1-{e}^{2}}}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{K}_{\mathrm{planet}}=\displaystyle \frac{2\pi a}{{P}_{\mathrm{orb}}}\displaystyle \frac{{m}_{\mathrm{star}}}{{m}_{\mathrm{star}}+{m}_{\mathrm{planet}}}\displaystyle \frac{\sin i}{\sqrt{1-{e}^{2}}}\end{eqnarray} \tag{ 8 }$$

where m_star and m_planet are the two masses of the star and planet, respectively. The true anomaly f, as a function of time t can be computed (see for example Mortier et al. 2016) from

$$\begin{eqnarray}&&\tan \displaystyle \frac{f}{2}=\sqrt{\displaystyle \frac{1+e}{1-e}}\tan \displaystyle \frac{E}{2}\end{eqnarray} \tag{ 9 }$$

where E, the eccentric anomaly, can be found solving the Kepler equation

$$\begin{eqnarray}&&E-e\sin E=2\pi \displaystyle \frac{t-{t}_{0}}{{P}_{\mathrm{orb}}}.\end{eqnarray} \tag{ 10 }$$

Knowing ϕ, the planet radial velocity relative to the star can be obtained as the difference between RV_planet and RV_star which could be also written as (see Charbonneau et al. 1999)

$$\begin{eqnarray}&&{\mathrm{RV}}_{\mathrm{planet},\mathrm{star}}=-{K}_{\mathrm{star}}\displaystyle \frac{{m}_{\mathrm{star}}+{m}_{\mathrm{planet}}}{{m}_{\mathrm{planet}}}\sin (2\pi \phi ).\end{eqnarray} \tag{ 11 }$$

This velocity can be used to compute the Doppler wavelength shift between the stellar and planet spectra at each ϕ in Equation (15).

The ICA technique was initially proposed to solve the BSS problem. In the BSS problem the task is to separate unknown individual signals (sources or also components) from mixtures of signals without a priori knowledge of the mixing process. In ICA, as the name implies, the basic goal is to find a linear transformation in which the underlying components are statistically as independent from each other as possible. ICA differs from the more common approach to component separation provided by Principal Component Analysis (PCA) which performs a linear transformation of the data to obtain mutually uncorrelated, orthogonal directions, called indeed principal components. If the hidden signals under investigation follow Gaussian distributions, uncorrelatedness is equivalent to mutual independence and algorithms such as PCA are able to separate them. However if the signals follow a non-Gaussian distribution (and most astronomically observed signals are predominantly non-Gaussian), it can be shown that uncorrelated signals are not necessarily mutually independent and hence methods like PCA fail to optimally separate individual components (Waldmann 2012). It is actually here where ICA-based separation methods could be of interest.

A comprehensive description of the ICA technique can be found in Hyvärinen et al. (2001). Hereafter we just recall that in ICA the model for the data can be expressed as

$$\begin{eqnarray}&&{\boldsymbol{X}}={\boldsymbol{AS}},\end{eqnarray} \tag{ 12 }$$

where ${\boldsymbol{X}}={({x}_{1},{x}_{2},\ldots ,{x}_{m})}^{T}$ are the observed m mixtures, ${\boldsymbol{S}}={({s}_{1},{s}_{2},\ldots ,{s}_{n})}^{T}$ are the latent variables (i.e., components) that cannot be observed and ${\boldsymbol{A}}$ is an unknown constant m × n matrix, called the mixing matrix. We use bold upper-case letters to denote matrices and lower-case letters (for example x₁) to denote, instead, vectors.

The challenge is to estimate the mixing matrix and its (pseudo) inverse de-mixing matrix, ${\boldsymbol{W}}={{\boldsymbol{A}}}^{-1}$, without any additional prior knowledge of either ${\boldsymbol{A}}$ and ${\boldsymbol{S}}$. The estimation of the matrix ${\boldsymbol{S}}$ with the knowledge of ${\boldsymbol{X}}$ is the linear source separation problem and can be achieved by maximizing some measure of independence. Several of such measures are used, like kurtosis, negentropy, or mutual information (Hyvärinen et al. 2001) and a large variety of algorithms implementing the above mentioned independence measures (e.g., Second Order Blind Identification, SOBI, Belouchrani et al. 1997; Joint Approximation Diagonalization of Eigenmatrices, JADE, Cardoso & Souloumiac 1993; fixed-point algorithm, FastICA, Hyvarinen 1999 and many others see Waldmann 2012) exists.

2.3.1. FastICA

In the present work we decided to use the FastICA algorithm since it is one of the fastest, most commonly used and it is well documented. It is also available in many programming languages and, in particular, we use the FastICA package for Matlab^TM. A detailed description of the algorithm can be found in Hyvarinen (1999). Here we just recall its main steps:

• 1.  
  data are centered (set them at mean zero);
• 2.  
  data are whitened (i.e., ${\boldsymbol{X}}$ data are transformed so that the components of the new vector $\tilde{{\boldsymbol{X}}}$ are uncorrelated and their variances equal to unity);
• 3.  
  the inverse de-mixing matrix ${\boldsymbol{W}}$ is estimated;
• 4.  
  the mixing matrix ${\boldsymbol{A}}$ and the component matrix ${\boldsymbol{S}}$ are computed.

The matrix ${\boldsymbol{W}}$ is estimated raw by raw with the following iterations scheme:

• (a)  
  choose initial (random) weight vector ${w}_{{\rm{i}}}$
• (b)  
  let ${w}_{i}^{+}=E[{\boldsymbol{X}}g({w}_{i}^{T}{\boldsymbol{X}})]-E[{g}^{{\prime} }({w}_{i}^{T}{\boldsymbol{X}})]{w}_{i}$, where g and g' are the derivatives of the chosen contrast function (see Equation (13))
• (c)  
  normalize w⁺_i
• (d)  
  repeat step (b) and (c) until convergence is achieved

where convergence means that the old and new values of w_i point in the same direction, i.e., their dot-product is (almost) equal to 1.

A contrast function G and its first and second derivatives, g and g', are generally used to approximate the negentropy of the system. In our FastICA package, the choice is restricted to the following functions:

$$\begin{eqnarray}\begin{array}{rcl}{G}_{1}(u) & = & \displaystyle \frac{1}{{a}_{1}}\mathrm{log}\,\cosh ({a}_{1}u)\\ {G}_{2}(u) & = & -\exp (-{u}^{2}/2)\\ {G}_{3}(u) & = & \displaystyle \frac{1}{4}{u}^{4}.\end{array}\end{eqnarray} \tag{ 13 }$$

The problem of most ICA algorithms is that they are based on methods related to gradient descent where the basic principle is to start in some initial point, and then make steps in a certain direction until a convergence criterion is met. In the case of the FastICA sequence the w_i initial vector of weights is generated at random and the algorithm stability is therefore not always deterministic. Moreover, it is known that Equation (12) led to two ambiguities:

• 1.  
  sign and variances of the independent components could not be determined;
• 2.  
  order of independent components could also not be determined.

Both ambiguities are directly related to the fact that both ${\boldsymbol{A}}$ and ${\boldsymbol{S}}$ are unknown. Any scalar multiplier in one of the components s_i could always be canceled by dividing the corresponding column a_i of ${\boldsymbol{A}}$ by the same scalar. FastICA assumes that each component has unit variance and correspondingly the matrix ${\boldsymbol{A}}$ is adapted to take into account this restriction; this solves point 1 apart for a sign ambiguity. Concerning the order, positions of independent components can be freely changed without affecting the Equation (12). To show this, it is enough to reshuffle Equation (12) by multiplying it with a permutation matrix ${\boldsymbol{P}}$ and its inverse as ${\boldsymbol{X}}={{\boldsymbol{AP}}}^{-1}{\boldsymbol{PS}}$. The matrix ${{\boldsymbol{AP}}}^{-1}$ is the new unknown mixing matrix with order of column position changed. A way to mitigate the FastICA instability and ambiguities is to run ICA algorithms many times and measure "somehow" the stability of the obtained components. This is the purpose of ICASSO, a software package aiming at investigating the relations among estimates from FastICA both programmatically as well as visually.

2.3.2. ICASSO

To represent a possible real case, we decided to simulate the 51 Peg + 51 Peg b planetary system, Indeed, this is also one of the first systems where reflected signal from an exoplanet has been detected at a significance of 3σ_noise (Martins et al. 2015). To properly characterize the system we use the parameters summarized in Table 1.

The model of the composite spectrum at different orbital phases is obtained, similarly as it was done in Borra & Deschatelets (2018), by adding to a stellar spectrum the same spectrum Doppler shifted and modulated by:

• 1.  
  a phase-dependent function
• 2.  
  the predicted geometric flux ratio of the planet relative to the star
• 3.  
  a reflecting albedo

as described in Sections 2.1 and 2.2.

In order to simulate the composite spectrum, we used as F_star(λ) the synthetic spectrum of 51 Peg computed by means of SPECTRUM v.276e (Gray & Corbally 1994) starting from an ATLAS12 atmosphere model (Kurucz 2005) at T_eff = 5787 K, log g = 4.45 dex, metallicity 0.15 dex and microturbulence 0.85 km s⁻¹ (data taken from Valenti & Fischer 2005). The synthetic spectrum was broadened to take into account the rotational velocity of 2.6 km s⁻¹, macroturbulence velocity of 3.95 km s⁻¹ and degraded at the resolution of ESO HARPS spectrograph (R ≈ 115,000). The same synthetic spectrum, but Doppler shifted, was used to compute the F_planet(λ). The planetary albedo p(λ) was assumed to be equal to the Neptune one given by Karkoschka (1998) and it is shown in Figure 1. More recent measurements by Madden & Kaltenegger (2018) are also available in the literature, but we decided to use the ones by Karkoschka (1998) due to their more extended wavelength coverage toward the blue. We decided to use the Neptune albedo because in the wavelength range of our simulation it shows prominent features with variations of the order of ±30%. The goodness of our method in extracting planet features can be in this way better assessed. Actually at the level of simulation what is important is to show the ability to recover the planet signal given in input; in principle, without any loss in generality, we could even not use any real, observed albedo, but simply adopt for p(λ) a monotonic function.

Zoom In Zoom Out Reset image size

In our simulation the spectra are limited to the 4900–5300 Å wavelength range and have a fixed step of 0.005 Å. The wavelength range has been chosen in a spectral region where there is a possibility of achieving good S/N in the case of real ground-base observations avoiding regions of telluric contamination and the presence of strong lines (e.g., hydrogen ones), yet keeping the computational time reasonable.

To create the synthetic observed mixture ${\boldsymbol{X}}$ (composed by different ${F}_{j}(\lambda )$), to be analyzed subsequently by ICA, one should note that it is not possible to use F_j(λ) for all different ϕ simultaneously. In fact, the contribution of the F_planet(λ, ϕ) in each F_j(λ), at different phase angles, would be seen by ICA as a different independent component (due to the different radial velocity shift). This will lead to an attempt to disentangle an observed mixture ${\boldsymbol{X}}$ composed, say, by mF_j(λ) in n = m + 1 components i.e., m F_planet(λ, ϕ) and one F_star, which has of course, no solution. We actually need to inject in ICA an observed mixture such that the dimensionality is correctly preserved. For our specific exoplanetary case this can be achieved, in the most simple case, by using a pair of spectra F_j(λ) where in one of them the exoplanet reflected spectrum is not present. The observed mixture ${\boldsymbol{X}}$ assumes in this case the form:

$$\begin{eqnarray}&&{\boldsymbol{X}}=\left(\begin{array}{c}{F}_{1}(\lambda )\\ {F}_{2}(\lambda )\end{array}\right)\end{eqnarray} \tag{ 14 }$$

where:

$$\begin{eqnarray}\begin{array}{rcl}{F}_{1}(\lambda ) & = & {F}_{\mathrm{star}}(\lambda )\\ {F}_{2}(\lambda ,\phi ) & = & {F}_{\mathrm{star}}(\lambda )+{F}_{\mathrm{planet}}(\lambda ,\phi )\\ & = & {F}_{\mathrm{star}}(\lambda )+p(\lambda ){\left[\displaystyle \frac{{R}_{\mathrm{planet}}}{a}\right]}^{2}g(\alpha ){F}_{\mathrm{star}}\\ & & \times \left(\lambda \left[1+\displaystyle \frac{{\mathrm{RV}}_{\mathrm{planet},\mathrm{star}}}{c}\right]\right).\end{array}\end{eqnarray} \tag{ 15 }$$

In Equation (15), ϕ ∈ [0, 1] is the usual orbital phase (see Equation (3)) and RV_planet,star is the radial velocity of the planet relative to that of the star (as derived in Equation (11)). F₁(λ) thus effectively represents a spectrum obtained at the planet occultation (ϕ = 0.5) if it is an eclipsing system, or at the inferior conjunction (ϕ = 0.0) if not, i.e., when no starlight is reflected by the planet toward the observer. F₂(λ), instead, contains also the contribution reflected by the planet at the corresponding phase ϕ. As discussed in Section 3 to check the independence of F_star and F_planet we computed $f({F}_{\mathrm{star}}| {F}_{\mathrm{planet}})$ and f(F_star) × f(F_planet) and we found an average difference, $\langle f({F}_{\mathrm{star}})| {F}_{\mathrm{planet}})-f({F}_{\mathrm{star}}\times f({F}_{\mathrm{planet}})\rangle$, equal to 2.7 × 10⁻⁶ ± 2.1 × 10⁻⁶. Such a low difference suggests that the two spectra are independent and thus we can be confident that ERLICA will be able to extract the two searched components, i.e., the stellar spectrum and the reflected light from the exoplanet properly modulated by the introduced albedo.

In order to simulate a realistic case, as discussed in Section 3, we added to both spectra realizations, by using Matlab^TM function normrnd, a Gaussian noise with mean zero and standard deviation F(λ)/S/N (where S/N is a parameter to specify the S/N of the mF_j(λ) spectra in the wavelength range of interest). Expressing the noise via S/N is convenient and allows us to estimate the validity of the obtained results as a function of varying S/N.

3.1.1. Extraction of ICA Components

In order to assess the behavior of ICA in disentangling the system individual components we started by applying ERLICA, i.e., FastICA+ICASSO, on the composite spectra obtained when 51 Peg and 51 Peg b are at the two phases presented in Figure 2. These configurations represent somehow the two competing "extreme" cases: the case of maximum available flux signal and still different planet and stellar radial velocities versus the maximum shift in radial velocity but a lower phase dependent flux ratio. The blue spectrum in the figure corresponds to the case of the exoplanet disk almost fully illuminated, ϕ = 0.45, but with a minimal separation in radial velocity (RV_planet,star ≈ 40.7 km s⁻¹). In this case, the phase dependent flux ratio (see Equation (1)) amounts to 2.9 × 10⁻⁴. The red spectrum, on the contrary, represents the case of maximum shift in radial velocity (ϕ = 0.25, RV_planet,star ≈ 131.7 km s⁻¹), but with a disk only partially illuminated leading to a phase dependent flux ratio of 9.7 × 10⁻⁵.

Zoom In Zoom Out Reset image size

As a starting point a S/N = 50,000 has been adopted to limit the influence of the noise and highlight the capabilities of the method.

In extracting the two independent components we have to face the two ambiguities described in Section 2.3. The first one, i.e., the sign of each component, can be solved taking into account that each F₂(λ, ϕ) must be the sum of two positive quantities. Therefore we determined the sign of each disentangled component by taking into account the sign of the corresponding element of the matrix ${\boldsymbol{A}}$.

The second ambiguity, i.e the order of the components, was resolved by comparing each of them with one of the two input spectra (each of them being dominated by the stellar signal) via cross-correlation. The output component with the highest value of the cross correlation is then labeled as component 1 and identified as the ICA estimate of F_star.

Setting up the simulation is however not enough: we need a way to judge and assess, with the results in hand, the "goodness" of ICA in extracting the two signals. We decided, therefore, to search for a quantitative, mathematical, estimator. A first attempt led us to calculate simply the standard deviation of the differences between the extracted components and the input F_star and F_planet. However, this estimator cannot be used in the case of real observations, since there will be no "a priori known" F_star and F_planet signals to be compared with. Therefore, we decided to use an approach similar to the one shown in Borra & Deschatelets (2018), i.e., to use, as a more reliable and feasible estimate, the peak intensity of the ACF of each obtained component. In fact, the higher this peak is, the higher the strength of the detected signal is as discussed in Section 4.3 of Borra & Deschatelets (2018) who adopted a similar approach to assess the detection of a planetary signal with their ACF method.

Therefore, we autocorrelate each of the two found components extracted by applying ERLICA on the composite spectra. The ACFs are computed, after resampling the components on a velocity scale, by relatively shifting them at steps, lag, using a constant velocity increment. The ACF profile is, by construction, always symmetrical and well centered at 0 position. Actually, correlating a signal by itself always produces a peak centered at 0 and equal to 1 after normalization. If there is a signal, i.e., a detection, it can be inferred from the ACF increase with respect to the one obtained if only noise is present at X values (shifts) equidistant and close to 0 both on the negative and positive sides. Therefore, to use the ACF as an estimator we removed the autocorrelation peak at 0, which is, for our purposes, not meaningful and substitute it with a value given by interpolation (usually with a Gaussian fitting) of the neighborhood values. The peak intensity of such a modified profile can be used as a robust estimator of the detection level when compared with the same quantity obtained in the case of an ACF computed when a detection is not possible. No-detection can be simulated by injecting in ERLICA two spectra F₁(λ) and F₂(λ) computed both at ϕ ≃ 0.0 thus using F_Planet = 0 in Equation (15). In the following, we use the ratio of the peak of the ACF of each i component over the no-detection ACF peak value to define the detection significance, D_i, and declare that there is a detection, if D_i is larger than 3.

In summary the main steps of the simulation are the following:

• (i)  
  create F₁(λ) and F₂(λ) according to Equation (15) for a chosen ϕ value;
• (ii)  
  add to F₁(λ) and F₂(λ) a Gaussian noise to mimic spectra with different S/N values;
• (iii)  
  create the observed mixture ${\boldsymbol{X}}$ and apply ERLICA;
• (iv)  
  fix signs and order of the extracted components;
• (v)  
  compute ACF of each component after re-sampling at steps of constant velocities;
• (vi)  
  compute detection significance D₁ and D₂ values.

3.1.2. Contrast Function Selection and S/N Effect

FastICA algorithm implementation allows the user to choose between three different contrast functions (see Equation (13)). The theoretical analysis that led to the adoption of the aforementioned contrasts function can be found in Hyvarinen (1999); here we want just to recall the main conclusions given there:

• 1.  
  G₁(u) is a good general-purpose contrast function;
• 2.  
  G₂(u) may be better when robustness is very important;
• 3.  
  G₃(u) is not recommended in the case of presence of outliers.

To chose the most appropriate contrast function for our scientific case we applied the ERLICA approach as described in Section 3.1 by using all the three contrast functions and varying S/N on the simulated 51 Peg + 51 Peg b planetary system with ϕ = 0.45. S/N has been varied in the range [1000–50,000] to span a noise amplitude interval from ten times to one fifth of the planet flux. For each value of S/N we computed the detection significance D_i of the disentangled components.

Figure 3 shows the results. Some points can be highlighted:

• 1.  
  detection level increases with increasing S/N, as expected, independently of the adopted contrast function;
• 2.  
  all the three contrast functions led to almost equivalent results, but for S/N = 5000. The G₂(u) "gauss" contrast function shows a non-monotonic behavior for S/N < 10,000 which seems to indicate a greater sensitivity to the FastICA initial conditions;
• 3.  
  in general a S/N > 5000 is required for a reliable detection (D > 3). This S/N limit corresponds to the case of a noise amplitude comparable to the planet signal one.

Zoom In Zoom Out Reset image size

Taking into account the main conclusions in Hyvarinen (1999) and considering that G₂(u) performs, sometimes, slightly better than the other contrast functions going toward low S/N we decided to use the G₂(u) "gauss" contrast function in our ICA analysis. It is worthwhile noticing that, as shown in our simulations and highlighted in Hyvarinen (1999), this choice of the contrast function is not really critical in the sense that any of the contrast estimators in the FastICA framework works well for (practically) any distributions of the independent components (contrary to what happens in other ICA algorithms).

3.1.3. Simulation Results

First of all, we apply ERLICA at the case ϕ = 0 where no detection of the second component is expected. The corresponding ACFs of the two component profiles are shown in Figure 4. In this case, ERLICA is clearly able to retrieve the first component, which has a well peaked ACF, whereas the ACF of the "second" one (actually ERLICA retrieves only noise) is flat and with a maximum value of about 0.01.

Zoom In Zoom Out Reset image size

Then, we repeat the analysis for the case of maximum available reflected flux (ϕ = 0.45) and for that of maximum shift in radial velocity (ϕ = 0.25) between F₁ and F₂. Figure 5 shows the first and second extracted components, and their ACFs, as obtained by our computation once FastICA is run with the Gaussian contrast function G₂(u) (see Equation (13)), S/N = 50,000, and ICASSO with M = 15 iterations. In both cases, as can be seen already by a visual comparison, ICA is able to disentangle the two components (see (a), (c) and (b), (d) panels, respectively).

Zoom In Zoom Out Reset image size

The comparison of the peaks of the ACF profiles, shown in panels (e) and (f), with the no-detection case shows that:

• 1.  
  in the case of the maximum available flux both components are retrieved and reconstructed; peak of the ACF profiles for the first component is close to 1, thus the detection can be considered optimal (D₁ ≃ 100); the detection of the second component is lower (D₂ ≃ 65), but still the albedo trend is clearly noticeable in panel (c);
• 2.  
  when the system is at ϕ = 0.25 the retrieved second component is much more noisy (see panel (d)) and this is reflected by a decrease in the peak of the corresponding ACF profile. Nevertheless, also in this case the detection of the planet signal is achieved, D₂ ≃ 20.

Figure 6 shows results obtained by using a S/N = 6500 where the detection of the second component could still be considered achieved (D₂ ≃ 3). As shown in panel (c), the peak of the ACF profile for the second component (green curve) indeed still slightly exceeds our assumed limit of 3 for a reliable detection and confirms what was shown also in Figure 3. It is worth highlighting that, also in this case of much lower S/N, the first component (panel (a)) is still optimally retrieved (D₁ ≃ 100). On the other hand no conclusion can be inferred visually for the albedo wavelength dependency (see panel (b)) without applying some procedure of noise filtering.

Zoom In Zoom Out Reset image size

In conclusion, our simulation results show that by using ERLICA we can effectively disentangle the individual components of a composite spectrum of an exoplanetary system like 51 Peg. For a composite spectrum with very high S/N (>6500) we demonstrate that we can be able to detect the planet signal in a system with a flux ratio on the order of 10⁻⁴. We also show that if the noise in the input spectra is of the order or smaller than the planet signal our method is also capable of providing the planet reflected spectrum i.e., the albedo wavelength dependence without further processing. For lower S/N the estimate of the wavelength dependence of the planet reflected spectrum would require some extra processing to increase its S/N.

In order to test our method on real data we decided to look at the case of eclipsing binary stars which has the advantage of a higher flux ratio than in the case of an exoplanetary system. It should be noted indeed that for our approach the differences between an exoplanetary (binary) system and an eclipsing star (binary) system are really minor:

• 1.  
  all the theoretical framework for radial velocity computation described in Section 2.2 applies, without losing in generality, to systems composed by two stars (basically, in the formulas it is enough to replace the subscript planet with star2);
• 2.  
  all the photometric variations due to phase dependent reflected light (Section 2.1) are, on the contrary, not relevant and can be neglected but the two spectra F_star1(λ) and F_star2(λ) can still be considered as two independent components.

The observed mixture ${\boldsymbol{X}}$ in this case is composed by the following two signals:

$$\begin{eqnarray*}\begin{array}{rcl}{F}_{1}(\lambda ) & = & {F}_{\mathrm{star}1}(\lambda )\\ {F}_{2}(\lambda ,\phi ) & = & {F}_{\mathrm{star}1}(\lambda )+{F}_{\mathrm{star}2}(\lambda ,\phi )\\ & = & {F}_{\mathrm{star}1}(\lambda )+{F}_{\mathrm{star}2}\left(\lambda \left[1+\displaystyle \frac{{\mathrm{RV}}_{\mathrm{star}2}(\phi )-{\mathrm{RV}}_{\mathrm{star}1}(\phi )}{c}\right]\right)\end{array}\end{eqnarray*}$$

where F₁(λ) is the spectrum of the system in secondary eclipse, 0.45 < ϕ < 0.55 (i.e., a spectrum taken when the secondary star is hidden behind the primary), and F₂(λ, ϕ) is the combined spectrum of the primary and of the secondary obtained at $0.1\lt | \phi | \lt 0.4$, thus avoiding the spectra taken during the primary eclipse where the problem of limb-darkening (see Winn 2010) and the Rossiter–McLaughlin (McLaughlin 1924; Rossiter 1924) effect modify the observed F_star1. We recall that the F_star2 contribution to F₂(λ, ϕ) at each phase is Doppler shifted in wavelength because of the radial velocity of the secondary star with respect to the primary one. The required statistical independence of the two components hidden in the composite spectra is guaranteed because of the shift in radial velocity of the two stellar spectra and even reinforced by the fact that, in most astronomical cases, the two stars belong to different spectral classes and, therefore, have quite different spectra.

Out of the several eclipsing binary systems described in the literature we decided to study R CMa. The eclipsing binary star R CMa is a short-period Algol-type system showing an extraordinary small mass ratio between its components. R CMa was known for a long time as the system of lowest total mass and as the prototype of a small group of stars called the R CMa-type stars, introduced by Kopal (1956) and characterized by low mass ratio, overluminosity of the primary, and oversized secondary. Budding & Butland (2011) give a comprehensive overview of the history of the investigations of R CMa. They performed a combined photometric, astrometric, and spectroscopic analysis of the R CMa system and end up with the stellar parameters given in Table 2.

Lehmann et al. (2018) used time series of high-resolution spectra and analyzed the decomposed spectra of the components together with the radial velocities obtained from decomposed, least-squares deconvolved mean line profiles (LSD profiles, see Donati et al. 1997). Their results confirm the values given by Budding & Butland (2011) for the masses and radii, and also for the T_eff of the secondary component, whereas the T_eff derived for the primary component is by 300 K lower (see Table 2). The authors did not find evidence of the presence of a third body in the system, as supposed by Radhakrishnan et al. (1984) or Ribas et al. (2002).

We have chosen R CMa as our test star because of its luminosity ratio,$\tfrac{{F}_{\mathrm{star}2}}{{F}_{\mathrm{star}}}\simeq 0.04$ in the visible, and because we can compare the results of our method with those obtained by Lehmann et al. (2018). The latter authors used high-resolution spectra of R CMa obtained with the HERMES spectrograph (Raskin et al. 2011). The spectra were reduced using the standard HERMES pipeline and, subsequently, normalized to the local continuum. Then the Fourier transformation-based KOREL program (Hadrava 1995, 2006) was used to disentangle the observed composite spectra. In fact, from a time series of spectra, the program delivers the decomposed spectra of the components, normalized to the common continuum of both stars, together with the optimum orbital elements, assuming pure Keplerian orbits. The spectra of the components, resulting from observations in all out-of-eclipse phases, were renormalized to the individual continua with the help of the wavelength-dependent continuum flux ratio which was derived from spectrum analysis.

In our analysis we used the same HERMES normalized spectra used by Lehmann et al. (2018) and we applied our method to each j possible pair of F₁(λ) and F₂(λ, ϕ) spectra. Thus, after every ERLICA run, i.e for each j pair of spectra, we checked the detection significance of the derived components, ${S}_{1}^{j}$ and ${S}_{2}^{j}$. Then we properly averaged them to build $\langle {S}_{1}\rangle$ and $\langle {S}_{2}\rangle$, i.e., our final estimates of F_star1 and F_star2, respectively. In making the average of the second component estimates we used the RV_star2,star1 from Lehmann et al. (2018) to put all the individual ${S}_{2}^{j}$ in the reference frame where RV_star2 = 0. The results are shown in Figure 7: in the top panel we plot $\langle {S}_{1}\rangle$ compared with the synthetic spectrum of the primary star computed using its atmospheric parameters given in (Lehmann et al. 2018, Table 1), in the middle panel $\langle {S}_{2}\rangle$ is compared with the corresponding synthetic spectrum, and in the bottom panel we show $\langle {S}_{1}\rangle$ and $\langle {S}_{2}\rangle$ ACFs. To evaluate the detection significance we show in the bottom panel also the ACF of the false $\langle {S}_{2}\rangle$ we obtained by using k pairs built with two spectra both taken during the secondary eclipse. As can be seen the $\langle {S}_{1}\rangle$ and $\langle {S}_{2}\rangle$ are in very good agreement with the corresponding synthetic spectra and their detection significance is D₁ ≃ 59 and D₂ = ≃54, respectively.

Zoom In Zoom Out Reset image size

A comparison of our results and those obtained by Lehmann et al. (2018) using KOREL is shown in Figure 8. As can be seen there is a very good agreement for the Primary spectrum (with an rms value of 0.01) and a satisfactory agreement for the Secondary spectrum (rms = 0.10). The obtained rms values are of the same order of those between the derived spectra and the corresponding synthetic ones. We recall that the validity of the synthetic spectra is limited by the physics included in the corresponding models and by the different program codes used to calculate them. Figure 8 demonstrates that our method provides estimates of the disentangled spectra as reliable as those obtainable by well proofed programs used in binary star analysis.

Zoom In Zoom Out Reset image size

In conclusion, combining the results of our simulations (see Section 3) and those obtained in the case of R CMa, we can say that the validity of our method is assessed.

As shown in Section 3.1.2, to reliably disentangle individual components of binary systems using our ICA-based method, it is mandatory to have at our disposal spectra with very high S/N. This is particularly demanding in the case of exoplanetary systems which are characterized by a very low flux ratio (of the order of 1 × 10⁻⁴ or less); simulations (see Section 3.1.3) show that, for a typical case, an S/N > 5000 is required. To reach such an S/N is probably outside the possibility of the current instrumentation and certainly not available in currently public available data.

A way to mitigate such limitations with the aim of also successfully applying our method to such demanding systems is to try to increase somehow the S/N. This can be achieved, for example, by:

• 1.  
  using averages of multiple spectra instead of single ones for the F₁(λ) and F₂(λ) in Equation (14) with the constraint that they should be taken at exactly the same ϕ values;
• 2.  
  using more than one pair of F₁(λ) and F₂(λ) and then averaging the retrieved components as done in Section 4.1.

Based on these considerations we decided to test our method on the real data of an exoplanetary system and in particular we decided to use 91 HARPS spectra of the system 51 Peg + 51 Peg b already used by Martins et al. (2015) and by Borra & Deschatelets (2018). The spectra were re-reduced by using HARPS DRS 3.4 and kindly made available to us by J.H.C. Martins. Martins et al. (2015), using the CCF method, and Borra & Deschatelets (2018), using the ACF method, detected from the analysis of these spectra the signal of 51 Peg b with a detection significance of 3.70 σ_noise and 5.52 σ_noise, respectively.

We limited our analysis to the wavelength region 5400 Å < λ < 6800 Å where the spectra have S/N > 200. Unfortunately this wavelength range is affected in several regions by the presence of telluric lines and removing them is not an easy task (see discussion in Smette et al. 2015). The method we applied is based on the fact that the HARPS spectra were obtained on different Julian days and, thus, they are affected by different heliocentric velocities. Therefore, after correcting all the spectra for the proper heliocentric velocity to put them in the wavelength laboratory rest frame, and after normalization, each spectrum shows the contamination of telluric lines at different wavelength positions (see upper plot in the upper panel of Figure 9). These wavelength positions can be easily identified by their anomalously large standard deviations with respect to the mean spectrum. Then, by sorting at each individual position the normalized flux of all the spectra, we separated the spectra which, in that specific point, have the higher signals from the others. The former are those less affected by telluric lines and their mean intensity value was used as an estimate of the telluric corrected flux value. Eventually, this value was adopted to substitute, at the corresponding wavelength point, the signal in the latter ones.

Zoom In Zoom Out Reset image size

An example of the correction procedure in one of the regions affected by telluric lines is shown in Figure 9 where:

• 1.  
  the upper panel shows three spectra and the differences from their mean before our telluric line correction;
• 2.  
  the central panel shows the normalized sorted flux values at one of the λ in the region, λ₀ = 6278.84 Å. This plot is used to select the "uncontaminated" spectra (green dots) and to compute an estimate of the telluric corrected flux value at λ₀(red line). The number of "uncontaminated" spectra, 12, is a good compromise obtained by visually inspecting several plots like that one shown;
• 3.  
  the bottom panel shows the same three spectra of the left panel after the correction with, at the bottom, the new differences from their mean.

As can be seen the adopted method works quite well even if some residual contamination (of the order of a few percent in normalized fluxes) still remains in very critical regions (see lower plot in the bottom panel of Figure 9).

To perform our ERLICA analysis the 91 HARPS spectra were Doppler shifted for the RV_star values and divided into two groups:

• 1.  
  Group 1: contains the 20 spectra that are near the inferior conjunction; they can be used, individually, as F₁(λ) in Equation (14);
• 2.  
  Group 2: contains the 71 spectra that are supposed to contain both F_star and F_planet since they have been observed out of the inferior conjunction; they can be used as F₂(λ) in Equation (14).

Now we are able to follow the steps outlined in Section 3.1. In particular:

• create the observed mixture ${\boldsymbol{X}}$ using for F₁(λ) a spectrum pertaining to group 1 and for F₂(λ) a spectrum of group 2 and apply ERLICA;
• fix signs and order of the extracted components S_star and S_planet to obtain estimates of F_star and F_planet;
• compute the ACF of each component after re-sampling at steps of constant velocities;
• compute detection significance D_star and D_planet values.

With the aim of increasing the S/N we paired each of the 20 F₁(λ) spectra individually with all the 71 spectra of group 2. Eventually we were able to apply ERLICA on 20 × 71 = 1420 mixtures X_j and for each j pair of F₁(λ) and F₂(λ) spectra, we obtained the ${S}_{\mathrm{star}}^{j}$ and ${S}_{\mathrm{planet}}^{j}$ estimates of F_star and F_planet. Moreover, to compute D_star and D_planet, we applied ERLICA on the 190 k pairs built using as F₂(λ) a spectrum of group 1 in order to check our F_planet estimate. In fact, in these cases we should not detect any signature from the planet and the analysis of the k obtained second ICA components can give estimates of possible "fake/spurious" S_planet introduced by ERLICA.

4.2.1. Results

Figure 10 shows the average ACFs of the ${S}_{\mathrm{star}}^{j}$, ${S}_{\mathrm{planet}}^{j}$, and "fake" ${S}_{\mathrm{planet}}^{k}$. As can be seen the mean ACF of ${S}_{\mathrm{star}}^{j}$ shows that our procedure derives a well defined S_star component which provides a very accurate estimate of F_star. The identification of S_star with F_star only is confirmed by the absence of any planetary signal if we apply to the ${S}_{\mathrm{star}}^{j}$ ACFs the same analysis used by Borra & Deschatelets (2018) to detect the signal from 51Peg b (see their Figure 7)

Zoom In Zoom Out Reset image size

As far as the average ACF of the ${S}_{\mathrm{planet}}^{j}$ is concerned the computation of its detection significance, as described in Section 3.1.1 requires removal of the central peak. Actually both the ACFs of the ${S}_{\mathrm{planet}}^{j}$ and of the "fake" ${S}_{\mathrm{planet}}^{k}$ show central peaks which extend from −2 to +2 km s⁻¹ suggesting the presence of some correlated noise. Thus we removed this velocity interval before computing the gaussian fits shown by the blue and yellow lines in Figure 10. The obtained D_planet = 2.3 value is below our adopted detection threshold (D = 3) and this result is qualitatively in agreement with those obtained from our simulation (see Figure 3). In fact, even if we assume that our F₁(λ), F₂(λ) pairs are uncorrelated (which is quite unlikely) the S/N of our spectra, ∼200, should be multiplied by the square root of the number of pairs, $\sqrt{1420}$, leading to a value of S/N ∼ 7500 which is slightly above the S/N limit from Figure 3. Thus, taking into account the not complete independence of our F₁(λ), F₂(λ) pairs and the different values of the S/N of the individual HARPS spectra, we can conclude that a low detection significance was expected. Furthermore the spectra coverage of the planetary orbit is neither complete nor homogeneous making the used HARPS data not an optimal set for our method (in particular we have very few spectra taken close to the superior conjunction or at phases near the maxima of abs(RV_planet,star)). In any case, it is worthwhile noting that if we had used an estimator of detection similar to those used in the CCF and ACF approaches (for example the ratio between the peak of the fitted ACF and the standard deviation of the ACF pixel intensity for abs(Δ vel) > 10 km s⁻¹) we would have obtained D_noise = 14.3 σ_noise, i.e., a value larger than those obtained by Martins et al. (2015) or by Borra & Deschatelets (2018).

Even if the detection significance, D_planet, of the second component is quite low we tried to derive its wavelength dependence in order to understand its nature, i.e., if it contains the F_planet signature or if it is mainly due to systematics. To do that the individual ${S}_{\mathrm{planet}}^{j}$ estimates were averaged after applying the proper Doppler shifts to put them in a reference system where RV_planet = 0.

Figure 11 illustrates the final results:

• 1.  
  in the top panel we plot in green the average estimate of ${S}_{\mathrm{star}}^{j}$, the synthetic spectrum of 51 Peg used in Section 3.1 in light blue, and, to show the regions affected by telluric line contamination we may have not corrected perfectly, a scaled (2×) telluric spectrum computed with SKYCALC⁵ in red;
• 2.  
  the bottom panel contains the average estimate of S_planet^j in black, the average estimate of ${S}_{\mathrm{star}}^{j}$ in green, and the SKYCALC spectrum in red. The black curve was smoothed by using a 200 points running average to reduce the noise and shows some relatively broad features.

Zoom In Zoom Out Reset image size

Whether the average estimate of ${S}_{\mathrm{planet}}^{j}$ is the real reflected signal from 51 Peg b is however questionable. In an attempt to answer to this question we compare in Figure 12 the average estimates of ${S}_{\mathrm{planet}}^{j}$ (black) with the "fake" ${S}_{\mathrm{planet}}^{k}$ (yellow), derived from spectra in the inferior conjunction where the reflected F_planet signature cannot be present. Due to the centering and whitening of the input data and of the ambiguities intrinsic in ICA described in Section 2.3 we did not try to recover the absolute value of the black and yellow spectra and, therefore, both have zero mean value (in Figure 12 the spectra were vertically shifted to increase the readability).

Zoom In Zoom Out Reset image size

As can be seen, there is a very low similarity between the star spectrum and the extracted and smoothed second component (see lower panel of Figure 11). In particular the strongest lines in the 51 Peg spectrum, i.e., the Na D doublet and Hα are not present in the black curves of Figure 12. This is in contradiction with the expected results (see Figure 5) where the reflected signal should consist of the stellar spectrum modulated by the planetary albedo. On the other hand, the presence of more evident features and the corresponding larger standard deviations of the black curves with respect to the yellow ones in Figure 12 seems to suggest that we indeed detect some signal from 51 Peg b. If this is true a possible explanation for the absence of the stellar lines in our average second component could be the presence of clouds in the atmosphere of 51 Peg b which may smooth and flatten the planet reflected spectrum (see discussion, for example in Gao et al. 2017). Another possibility is that the "detected" features are the remnant of not completely removed telluric lines since the strongest of them fall in the critical regions where the SKYCALC spectrum shows most of the telluric lines (see Figure 12).

In conclusion, we do not have a sound final answer about the nature of the features in the average extracted second component. The possible detection of the reflected spectrum of 51 Peg b to be confirmed would require us to repeat our analysis using new "ad hoc" obtained input data, i.e., at higher S/N, covering the whole range of orbital phases, in particular both conjunctions, with more than one spectrum at each phase, and with auxiliary spectra to be used for accurately removing the telluric contamination.