Minimal Data Fidelity for Stellar Feature and Companion Detection

The study of spots and faculae on Sun-like stars underlies a key feature of exoplanet research, which has seen substantial growth in the past decade, motivated by fundamental biological and physical questions (Lanza et al. 2008; Czesla et al. 2009; Zellem et al. 2017). Theoretical studies substantially precede observational evidence (Schatten et al. 1978; Zwaan 1978; Hoyt & Schatten 1993; Fröhlich & Lean 2004; Foukal et al. 2006; Shapiro et al. 2016) that, due to the advent of large telescopes and satellites, has systematically emerged during the past several decades. Understanding the properties of a star along with its variability (and companions if present), is foundational to such studies.

The most easily observed surface features of stellar variability are spots (dark regions) or faculae (light regions). The general origin of these features is associated with magnetohydrodynamic and thermal convective processes, although their precise origin is the subject of ongoing study (see, e.g., Berdyugina 2005; Choudhuri 2016; Hanasoge et al. 2016; Schumacher & Sreenivasan 2020, and references therein). Although these stellar features have enabled us to study many processes, such as differential rotation and stellar evolution, with respect to the detection of exoplanets, they have also been viewed as obstructive noise. This latter perspective derives from the fact that spots and faculae can mimic the presence of an exoplanet, due to both changes in the radial velocity (RV) of the star and by masquerading as a transit event.

Detection of exoplanets is a major accomplishment in physical science (Wolszczan & Frail 1992; Mayor & Queloz 1995; Borucki et al. 2010). The two most common methods of detecting exoplanets are the Radial Velocity Method and the Transit Method (e.g., Aigrain et al. 2004; Fischer et al. 2014, 2016). The RV method measures the Doppler shift in the stellar spectrum, which is related to the relative center of mass motion of the star–planet system (see Figure 1(a)). This is a practical approach only when a Jupiter sized planet has an orbit at the distance of Mercury, and hence such exoplanets are typically called hot Jupiters. Although in principle the method is simple, false detections commonly occur due to overfitting (Dumusque et al. 2012; Rajpaul et al. 2016). Moreover, due to the high precision of the measurements required to detect the movement of the star, the method is very sensitive to noise (Struve 1952; Lovis & Fischer 2010; Lanza et al. 2011; Aigrain et al. 2012). Therefore, the approach requires the removal of noise from the "raw" signal, the sources of noise including stellar surface activity (granulation, stellar spots, faculae), long-term stellar activity, instrumental noise, and atmospheric/telluric noise. The transit method is based upon the decrease in the intensity flux coming from a host star as a planet passes between the star and the observer. Because this decrease in flux is typically about 0.5%–2%, this method is also very sensitive to noise. Additionally, intrinsic physical phenomenon, in particular spots and faculae, can easily mimic a planet, leading to a high false-detection rate (see Figure 1(b); Desort et al. 2007; Oshagh et al. 2014; Korhonen et al. 2015; Ligi et al. 2015).

The rotation of a "clean" star, i.e., a star without surface spots or faculae, will have a net zero shift in the spectrum. Namely, the blueshift caused by the incoming half of the stellar surface is balanced by the redshift of the outgoing half. However, this balance can be broken due to the presence of spots or faculae, thereby producing a net RV curve for the star. Because this RV can mimic the presence of a planet, it is important to be able to distinguish between the two signals.

Regardless of whether one is studying stellar features or detecting exoplanets, one seeks data with high resolution and large signal-to-noise ratio (S/N). We use the Spot Oscillation And Planet (SOAP) 2.0 tool (Dumusque et al. 2014), described in Section 2, to get clean simulated data of a star. This allows us to simulate different noise strengths that may be present in real data and analyze the data similar to a real data set. Data degradation can arise either due to unavoidable instrumental problems or noise sources, such as telluric contamination or light from a nearby star. This has often led to removal of data sets from analyses so that only clean data sets are examined. Here we quantitatively examine how noisy the data can be and yet still contain sufficient statistical information about the stellar features present. However, rather than removing the noise from the signal, we treat it as a source of information. Central to our approach is to utilize the separation of timescales of the noise arising from different physical phenomena. This approach based on a temporal multifractal framework is described in Section 3.1.

We use simulated data from the SOAP 2.0 tool (Dumusque et al. 2014). This tool simulates a rotating star using a grid of cells filled with a clean solar spectrum from the National Solar Observatory. A spot or a facula is simulated by replacing the spectrum in one or more of the grid cells (depending on the size of the feature) by a spectrum of the spot or facula (with flux scaled in the spot spectrum according to the contrast ratio between a facula and the photosphere), which then follows the rotation of the star. The final spectrum is obtained by summing all the spectra in the grid. These spectra are noise-free (S/N ∼1000) and have very high resolution (the spectra constitute ∼500,000 data points, and the spectral resolution of a spectrograph refers to the precision with which it distinguishes wavelengths in the spectrum, defined as ρ = λ/Δλ, where Δλ is the smallest difference in all resolved wavelengths in a spectrum. As resolution increases, so do the number of wavelengths resolved). Thus we obtain three different types of spectra (spot, facula, planet) from this tool, where we can change the size of the feature or the planet induced RV. For simplicity, the spot and facula are present on the equator of the star and the axis of rotation of the star and that of the planetary orbit (circular) are kept fixed at 90°. We analyze three data sets with the following parameters:

• 1.  
  Spot: A = 1%,
• 2.  
  Facula: A = 1%, and
• 3.  
  Planet: K = 10 m s⁻¹,

where $A=\tfrac{{R}_{{\rm{S}}/{\rm{F}}}^{2}}{{R}_{\mathrm{Star}}^{2}}\times 100 \%$ is the percent area of the stellar disk covered by the feature, K is the RV, R_S/F is the radius of the spot or facula and R_Star is the radius of the stellar disk. These spectra have a period of 25 time units, which are then stacked eight times to produce 200 spectra in time (Agarwal et al. 2017). The spot and facula coverage parameters above are typical for an "active" star observed via RV detection methods (Fischer et al. 2016; Giguere et al. 2016). For the planet RV, while the current state-of-the-art RV precision is at 1 m s⁻¹ (Fischer et al. 2016), many stars have been observed to be distributed between 1 and 10 m s⁻¹, with the S/N required to be very high at low RV. Figure 2 shows the normalized intensity of a star with a planet, with RV 10 m s⁻¹.

In reality, we never have clean spectra and it has been a technological challenge in instrument design to build spectrographs with extremely high resolutions. Therefore, following the methods of Davis et al. (2017), we (a) decrease the resolution of these clean spectra and (b) add noise to them to model real observations. The resolution is decreased by convolving the original spectrum with a Gaussian having a full width at half maximum equal to λ/ρ, where λ is a wavelength in the spectrum and ρ is the desired resolution. Gaussian white noise per pixel is then added in the spectrum depending on the desired S/N. Steps (a) and (b) are performed for all three data sets, with the aim to obtain the lower limits of S/N and ρ for such features to be detectable. For ease of computation, the final spectrum is then averaged to obtain a wavelength resolution of 1 Å. We have also examined the original data set to show that this averaging does not impact the results (Agarwal et al. 2017).

3.1. Multifractal Temporally-weighted Detrended Fluctuation Analysis

The essence of a multifractal system is the exhibition of different self-similar structures on multiple different spatial or temporal (or both) scales (e.g., Halsey et al. 1986). In highly nonlinear coupled dynamical systems, the intrinsic dynamics and the extrinsic effects, such as noise, nearly always operate on multiple spatial and temporal scales, rather than exhibiting a single scaling. This underlies the use of a multifractal method.

Many methods used to detect stellar features are influenced by the noise present in the data and therefore include a noise reduction strategy. To remove the bias induced by the noise removal strategy, or to define what qualifies as noise, we instead use a method that utilizes the statistical information present in the noisy data to characterize a stellar feature. Given m evenly spaced in time spectra, with each spectra spanning L wavelengths, we construct L time series of length m for each of the wavelengths. Here, each wavelength is analyzed separately, providing robustness to the results, as well as allowing us to distinguish between features that are present across all wavelengths versus those that are wavelength specific.

Zoom In Zoom Out Reset image size

We analyze all L time series using multifractal temporally-weighted detrended fluctuation analysis (see Appendix A; Agarwal et al. 2012, 2017, and references therein), which does not a priori assume anything about the structure of the data, save for its multifractality. Thus, we begin with an agnostic view about the system timescales and their corresponding dynamics. Moreover, rather than "removing" noise from the signal, we use it as a source of information on the dynamics of the system on all timescales present.

The approach has four stages, which we describe in Appendix A. We produce a statistical measure called the fluctuation function, F_q (s), each moment of which, q, is assessed on multiple timescales, s. For intuition, one can think of the expectation value of F_q (s) as the weighted sum of the autocorrelation function. The dominant timescales in a system are the those where F_q (s) versus s changes slope and the individual slopes are associated with the statistical dynamics of a system.

3.2. Principal Component Analysis

Another method that is widely used to extract the dynamics of a system is principal component analysis (PCA; Pearson 1901), which has a history that predates the discovery of multifractals. In PCA one transforms the underlying data into a set of linearly independent orthogonal vectors, with their directions set by the directions of maximum variance. Given a set of observations, one computes the covariance matrix of the original data, Σ, upon which one performs singular value decomposition. The resulting diagonal matrix D, such that Σ = UDV^T , gives the ordered set of eigenvalues denoting the amount of variance in the corresponding eigenvectors that are the columns of U. PCA has been very successful in analyzing data to, for example, to separate the dominant climate modes and long-term global temperature patterns on Earth (Mann et al.1998; Mantua & Hare 2002; Moon et al. 2018), and to detect exoplanet and stellar features (Thatte et al. 2010; Poppenhaeger & Schmitt 2011; Davis et al. 2017). PCA has also been used widely as a noise reduction technique in which the complete system is projected onto the eigenvectors corresponding to a threshold (low) variability of the system. But this may also result in losing information about the system if the processes corresponding to "noise" actually represent the low-variability dynamics of the system. Importantly, PCA also assumes that the system is governed by a Gaussian distribution and that the dynamics corresponding to each of the different modes is linearly independent. These assumptions are obviously not met a priori in a highly nonlinear system such as Earth's climate or the dynamics of a star–planet system.

Zoom In Zoom Out Reset image size

Figures 3(a)–(c) show the crossover timescales from the fluctuation functions for the original data set with no noise for a spot, facula, and planet, respectively. The main characteristics of these plots are: (i) almost all wavelengths show the same crossover timescale of 25 units, which is the rotational period of the star in the case of a spot and facula, and the orbital period of the star in the case of the planet; (ii) the spot has a clean timescale of 25 units while the facula shows a wavelength specific pattern, since not all wavelengths are affected by it (Davis et al. 2017); and (iii) the planet has some scatter around the timescale of 25 units. This provides a test bed for detecting such features as we degrade the simulation data. We have also tested our method by performing continuum normalization of the spectra to account for changes in observational conditions, which we describe in Appendix B. We show that normalizing the spectra only acts as a filter on the data, and thus does not change our analysis and results, as described below.

Zoom In Zoom Out Reset image size

To examine how the resolution and noise impacts these detections, we vary the resolution from 20,000 to 200,000, the S/N from 15 to 60 in the case of the facula, spot, and S/N from 50 to 800 for the planet. We define a statistic Ω,

$$\begin{eqnarray}&&{\rm{\Omega }}=\displaystyle \frac{\mathrm{Number}\ \mathrm{of}\ \mathrm{wavelengths}\ \mathrm{with}\ {\rm{a}}\ \mathrm{timescale}\ \mathrm{of}\ 25\ \mathrm{units}}{\mathrm{Total}\ \mathrm{Wavelengths}},\end{eqnarray} \tag{ 1 }$$

which is the ratio of number of wavelengths with a timescale of approximately 25 units to the total number of wavelengths. Given some S/N and resolution of the underlying spectra, Ω is the probability of detecting a stellar feature or a planet in the data. Thus, depending on the goals of a particular application, one can define a threshold probability for detecting a feature, for example defining Ω_Th = 0.5 implies that any feature with Ω ≥ 0.5 would correspond to a robust detection. Clearly, any value of Ω_Th has an important trade-off in terms of balancing the false positives and false negatives.

To compare MF-TW-DFA with PCA, we construct a metric ${{ \mathcal F }}_{\mathrm{PC}}$, in a manner similar to N_PC of Davis et al. (2017). While N_PC gives approximately the number of principal components in the data that can be used to extract useful information about the system, it does not quantify that information. N_PC is a whole number, computed by rounding up the average correlations between the scores (defined as the projection of the data matrix onto its eigenvectors) of the noisy data and the ideal data. Therefore, we have constructed a metric ${{ \mathcal F }}_{\mathrm{PC}}$ that measures the fraction of the total variability in these aforementioned principal components, and thus the amount of information in the data as measured by PCA. Hence, ${{ \mathcal F }}_{\mathrm{PC}}$ provides a more refined description of the variability in the data than does N_PC, which loses such detail due to the rounding process. Since Ω is a measure of probability in wavelength space, it should also be thought of as the information content in the data as measured by MF-TW-DFA. Irrespective of the value of N_PC, both ${{ \mathcal F }}_{\mathrm{PC}}$ and Ω are defined on a scale of [0, 1] thereby allowing us to perform a direct comparison between the two methods.

Similar to N_PC, we first compute the correlation between the first 10 scores of the noisy data (given S/N and resolution) and the ideal data (original S/N and resolution). We then compute the p value and choose those eigenvectors/eigenvalues for which the correlation has a p value < 0.001 (Davis et al. 2017). The definition of ${{ \mathcal F }}_{\mathrm{PC}}$ is

$$\begin{eqnarray}&&{{ \mathcal F }}_{\mathrm{PC}}\equiv \sum _{i=1}^{10}\sum _{j=1}^{50}{{ \mathcal I }}_{{R}_{{ij}}}\displaystyle \frac{{f}_{{ij}}}{50},\end{eqnarray} \tag{ 2 }$$

where f_ij is the fraction of total variability given by the i_th eigenvector of the j_th ensemble member, R_ij is the absolute correlation coefficient between the i_th score, the j_th ensemble member of the noisy data, and the i_th score of the ideal data, and

$$\begin{eqnarray}{{ \mathcal I }}_{{R}_{{ij}}}=\left\{\begin{array}{cl}1 & \mathrm{if}\ \mathrm{the}\ p \mbox{-} \mathrm{value}\ \mathrm{for}\ {R}_{{ij}}\lt 0.001\\ 0 & \mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ 3 }$$

For the same parameters as in Figure 4 we show N_PC and the corresponding ${{ \mathcal F }}_{\mathrm{PC}}$ for a spot, a facula and a planet in Figures 5(a)–(f). Both Ω and ${{ \mathcal F }}_{\mathrm{PC}}$ provide high-resolution measures of the amount of information in the noisy data. Clearly, at all S/N values and resolutions, MF-TW-DFA performs at least an order of magnitude better than PCA. The superior performance MF-TW-DFA is principally due to the fact that the method (a) is not restricted to capturing linear dependencies in the data as is PCA, and (b) uses "noise" as a source of information.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.1. S/N versus Resolution

4.2. Fluctuation Strength

4.3. Differentiating Spots, Faculae, and Planets

We have used simulation data of a stellar spectrum from the SOAP 2.0 tool in the presence of a spot, a facula or a planet to demonstrate the fidelity of a multifractal methodology for detecting exoplanets. The motivation is to provide a framework that naturally deals with unavoidable noise in observational systems. The approach is unique conceptually in that it uses noise as a source of information rather than treating it as something to be filtered out. By using controlled simulation data we systematically varied both the resolution and the S/N to determine a lower limit on the resolution and S/N required to robustly detect features using our multifractal method. For any resolution above 20,000, we found that a spot and a facula with a 1% coverage of the stellar disk can be robustly detected for a S/N (per pixel) of 35 and 60, respectively, whereas a planet with a RV of 10 m s⁻¹ can be detected for a S/N (per pixel) of 600. These results are compared to the standard PCA method, which requires a S/N 100 times larger than our multifractal method to detect the same information. A key aspect of these results, of relevance to considering which methodologies are most appropriate for observational data, is whether they have sufficient fidelity to discern stellar features themselves from planets. The framework described in this paper can be applied to data from previous generation instruments for more in-depth analysis as well as in designing and deploying the next-gen instruments and satellites for maximum efficiency and fidelity. Finally, the method allows for a systematic examination both of intrinsic stellar dynamical processes and practical noise sources, such as telluric or instrumental noise, to help refine observational schema.

S.A. and J.S.W. thank Debra Fischer for a series of extremely helpful conversations and they acknowledge NASA Grant NNH13ZDA001N-CRYO for support. J.S.W. acknowledges the Swedish Research Council grant No. 638-2013-9243, and a Royal Society Wolfson Research Merit Award for support.

For a given set of evenly spaced spectra with L wavelengths, we construct L time series and analyze each time series independently. We analyze the time series for each wavelength using MF-TW-DFA (Agarwal et al. 2012, and references therein). This method has previously been used to study the temporal structure of the Arctic sea ice extent (Agarwal et al. 2012) and velocity fields (Agarwal & Wettlaufer 2017), exoplanet detection without the use of model fitting (Agarwal et al. 2017), and multidecadal global climate dynamics modes on Earth (Moon et al. 2018). Multifractal behavior has been studied using a variant of the current algorithm, known as Multi-Fractal Detrended Moving Average (Gu & Zhou 2010) in other astrophysical phenomena as well such as to study stellar dynamics (de Freitas et al. 2016, 2017a, 2017b, 2019a, 2019b, 2021; Belete et al. 2018; de Franciscis et al. 2018, 2019).

MF-TW-DFA involves four steps, described here:

• 1.  
  Given the original time series X_i , construct a nonstationary profile Y(i) as

  $$\begin{eqnarray}&&Y(i)\equiv \sum _{k=1}^{i}\left({X}_{k}-\overline{X}\,\right),\ \ \mathrm{where}\ \ i=1,...,N.\end{eqnarray} \tag{ A1 }$$

  The profile is the cumulative sum of the time series and $\overline{X}$ is the average of the time series X_i .
• 2.  
  This nonstationary profile is divided into N_s = int(N/s) nonoverlapping segments of equal length s, where s is an integer and varies in the interval 1 < s ≤ N/2. Each value of s represents a timescale s × Δt, where Δt is the temporal resolution of the time series. The time series has a length that is rarely an exact multiple of s, which is handled by repeating the procedure from the end of the profile and returning to the beginning, thereby creating 2N_s segments.
• 3.  
  A point-by-point approximation ${\hat{y}}_{\nu }(i)$ of the profile is made using a moving window, smaller than s and weighted by separation between the points j to the point i in the time series such that ∣i − j∣ ≤ s. Note that in regular MF-DFA, as opposed to MF-TW-DFA, nth order polynomials y_ν (i) are used to approximate Y(i) within a fixed window, without reference to points in the profile outside that window.A larger (or smaller) weight w_ij is given to ${\hat{y}}_{\nu }(i)$ according to whether ∣i − j∣ is small (large; Agarwal et al. 2012). This approximated profile is then used to compute the variance spanning up (ν = 1,...,N_s ) and down (ν = N_s + 1,...,2N_s ) the profile as

  $$\begin{eqnarray}&&\begin{array}{l}\mathrm{Var}(\nu ,s)\equiv \displaystyle \frac{1}{s}{\sum }_{i=1}^{s}\{Y([\nu -1]s+i)\\ \qquad \ -\ \hat{y}([\nu -1]s+i)\}{}^{2}\\ \qquad \ \mathrm{for}\ \nu =1,\ldots ,{N}_{s}\ \mathrm{and}\\ \mathrm{Var}(\nu ,s)\equiv \displaystyle \frac{1}{s}{\sum }_{i=1}^{s}\{Y(N-[\nu -{N}_{s}]s+i)\\ \qquad \ -\ \hat{y}(N-[\nu -{N}_{s}]s+i)\}{}^{2}\\ \qquad \ \mathrm{for}\ \nu ={N}_{s}+1,\ldots ,2{N}_{s}.\end{array}\end{eqnarray} \tag{ A2 }$$
• 4.  
  Finally, a generalized fluctuation function is obtained and written as

  $$\begin{eqnarray}&&{F}_{q}(s)\equiv {\left[\displaystyle \frac{1}{2{N}_{s}}{\sum }_{\nu =1}^{2{N}_{s}}\{\mathrm{Var}(\nu ,s)\}{}^{q/2}\right]}^{1/q}.\end{eqnarray} \tag{ A3 }$$

We study the behavior of F_q (s) with respect to the timescale s. In particular, the generalized fluctuation function scales with the timescale s for a given moment q, with the scaling exponents given by the generalized Hurst exponents h(q),

$$\begin{eqnarray}&&{F}_{q}(s)\propto {s}^{h(q)}.\end{eqnarray} \tag{ A4 }$$

These Hurst exponents convey the statistical information in the data at some timescale. If h(q) is independent of q the data is said to be a monofractal. The second moment of the fluctuation function is generally studied more than others, since it can be related to more common power spectrum. For power spectrum S(f) ∝ f^−β , where f is the frequency and β is its decay exponent, h(2) = (1 + β)/2 (Rangarajan & Ding 2000). Therefore, since for white noise β = 0, it gives h(2) = 1/2. Similarly for pink noise β = 1, which gives h(2) = 1, for red noise β = 2, giving h(2) = 3/2, and so forth. This allows us to characterize the data using noise characteristics on multiple timescales, where the dominant timescales in the data set are the points where the fluctuation function ${\mathrm{log}}_{10}{F}_{2}(s)$ changes slope with respect to ${\mathrm{log}}_{10}s$. We apply this method to each wavelength in the spectrum, giving us the dominant timescales of the complete spectrum.

Due to change in observational conditions such as variation in cloud cover, the observed flux across wavelengths in spectral data is normalized to reduce such observation effects. This normalization acts as an additional filter on the magnitude of the observed flux. But it should also be noted that this normalization may interact with the stellar feature under study, since the timescale over which normalization is performed may be larger than the timescale of the stellar feature itself and thus mask any such features, e.g., stellar spots and faculae. Thus, while for a planet orbiting a star the normalization would mask out only the observational noise, for stellar spots and faculae, depending on their dynamic timescales, normalization can also act to mask the intrinsic dynamics. Figure 6 shows the fluctuation functions for a spot (a), facula (b), and a planetary-RV signal (c). Here, the normalization was performed as the final step in data processing, as would be done in actual observations. We also see that these fluctuation functions in Figure 6 are in excellent agreement with those in Figure 4.

Zoom In Zoom Out Reset image size