Minimal Data Fidelity for Detection of Stellar Features or Companions

Technological advances in instrumentation have led to an exponential increase in exoplanet detection and scrutiny of stellar features such as spots and faculae. While the spots and faculae enable us to understand the stellar dynamics, exoplanets provide us with a glimpse into stellar evolution. While the ubiquity of noise (e.g., telluric, instrumental, or photonic) is unavoidable, combining this with increased spectrographic resolution compounds technological challenges. To account for these noise sources and resolution issues, we use a temporal multifractal framework to study data from the SOAP 2.0 tool, which simulates a stellar spectrum in the presence of a spot, a facula or a planet. Given these controlled simulations, we vary the resolution as well as the signal-to-noise (S/N) ratio to obtain a lower limit on the resolution and S/N required to robustly detect features. We show 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, for any spectral resolution above 20,000, while a planet with a radial velocity (RV) of 10 m/s can be detected for a S/N (per pixel) of 600. Rather than viewing noise as an impediment, our approach uses noise as a source of information.


I. INTRODUCTION
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 [54,69,89].Theoretical studies substantially precede observational evidence [62,63,66,81,83,90] which, 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., 51, 53, 65, 82, 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 of the star and by masquerading as a transit event.
Detection of exoplanets is major accomplishment in physical science [52,74,88].The two most common methods of detecting exoplanets are the Radial Velocity Method and the Transit Method [e.g., 49,60,61].The radial velocity method measures the doppler shift in the stellar spectrum, which is related to the relative center of mass motion of the star-planet system.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 [59,79].Moreover, due to the high precision of the measurements required to detect the movement of the star, the method is very sensitive to noise [50,68,71,84].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 [56,67,70,76].
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 blue shift caused by the incoming half of the stellar surface is balanced by the red shift of the outgoing half.However, this balance can be broken due to the presence of spots or faculae, thereby producing a net radial velocity curve for the star.Because this radial velocity 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).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 datasets from analyses so that only clean datasets are examined.Here we quantitatively examine how noisy the data can be and yet still contain sufficient statistical in- formation 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.

II. DATA
We use simulated data from the Spot Oscillation And Planet (SOAP) 2.0 tool [58].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 and is 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 radial velocity.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 datasets with the following parameters: 1. Spot: A = 1%, 2. Facula: A = 1%, and 3. Planet:

Star
× 100% is the percent area of the stellar disk covered by the feature, K is the radial velocity, 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 8 times to produce 200 spectra in time [47].
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. [55], 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 six datasets, 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 A. We have also examined the original dataset to show that this averaging does not impact the results [47].

A. Multi-fractal Temporally Weighted Detrended Fluctuation Analysis
The essence of a multifractal system is the exhibition of different self-similar structure on multiple different spatial or temporal (or both) scales [e.g., 64].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.
We analyze all L time series using Multi-fractal Temporally Weighted Detrended Fluctuation Analysis (MF-TWDFA) [See Appendix A, 46, 47, 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 time scales, s.For intuition, one can think of the expectation value of F q (s) as the weighted sum of the auto-correlation function.The dominant time scales 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.

B. Principal Component Analysis
Another method that is widely used to extract the dynamics of a system is Principal Component Analysis (PCA) [77], 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 Σ = U DV 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 [72,73,75], and to detect exoplanet and stellar features [55,78,85].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.

IV. DISCUSSION
Figures 2 (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 [55], 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.
To examine how the resolution and noise impacts these detections, we vary the resolution from 20,000 to 200,000 and vary the S/N from 15 to 60 in the case of the facula and the spot and S/N from 50 to 800 for the planet.We define a statistic Ω, Ω = Number of wavelengths with a timescale of 25 units Total Wavelengths , (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 Ω T h = 0.5 implies that any feature with Ω ≥ 0.5 would correspond to a robust detection.Clearly, any value of Ω T h 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 F P C , in a manner similar to N P C of Davis et al. [55].While N P C 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 P C 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 F P C 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, F P C provides a more refined description of the variability in the data than does N P C , 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 P C , both F P C and Ω are defined on a scale of [0,1] thereby allowing us to perform a direct comparison between the two methods.
Similar to N P C , 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 pvalue < 0.001 [55].The definition of F P C is 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 and the j th ensemble member of the noisy data and the i th score of the ideal data, and For the same parameters as in Figure 3 we show N P C and the corresponding F P C for a spot, a facula and a planet in Figures 4 (a-f).Both Ω and F P C 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.
Figures 3 (a,c,e) show contour plots of Ω as a function of the resolution ρ and S/N.The robustness of our method is clear: We can detect a spot with S/N as low as 15 while for S/N ≥ 35, the results are independent of resolution.Because a facula exhibits noise behavior that differs from a spot, detecting it requires a higher S/N value, but still as low as 45.It is useful to compare this with Fig. 7 of Davis et al. [55] (their y-axis is S/N per resolution element), which is a PCA based study.
By comparing the lines of equal photon flux (S/N ∝ 1/ρ) with the contour lines in their Fig. 7, Davis et al. [55] demonstrated that in order to identify stellar signals in spectral data, having high resolution is more important than having high S/N.In contrast, we see that the red dashed lines in Figure 3 (a,c,e) are the lines of equal photon flux and run parallel to the contour lines of Ω, showing that with MF-TW-DFA low-resolution spectra can be successfully utilized to capture and study these stellar phenomena.This has significant implications for the analysis of both existing datasets and incipient missions, for which we can utilize the multifractal methodology described here to extract high quality information from what may have previously been viewed as noisy and/or low-resolution data.

B. Fluctuation Strength
Figures 3 (b,d,f) show the second moment of the fluctuation functions from Equation (A4) for all wavelengths.The main characteristics of these plots are: (i) Ostensibly all wavelengths are parallel to each other, with only a minute number of wavelengths deviating from the bulk (facula and planet) (see also Fig. 2), (ii) a substantial decrease in the fluctuation strength as we transition from spots to faculae to planets.The implication of this decrease in fluctuation strength is an increase in the sensitivity of fluctuation functions to noise, thereby requiring a higher S/N to robustly detect a feature.
Figures 2 and 3 (b,d,f) exhibit the essential property of the individual phenomenon discussed above that differentiates spots, faculae and planets.Namely, the abruptness and clarity of the crossover timescales (Fig. 2) and the strength of the fluctuations (Figs. 3 (b,d,f)).In the SOAP2.0code, a simulated facula is generated using a slight alteration of the spot spectrum, as compared to the intrinsic line-by-line variability between a sunspot spectra and the photosphere.Using PCA Davis et al. [55] found that the principal components of a spot are quite similar to those of a facula.Here we have shown that MF-TW-DFA is very efficient at detecting the small fluctuations for the facula (Fig. 2) and thus is able to differentiate between a spot and a facula.The fluctuation functions also differ with respect to the inter-wavelength spread, which is much smaller for the spot than for the planet.This spot-to-planet increase in inter-wavelength spread is an important quantitative detection diagnostic intrinsic to MF-TW-DFA.Indeed, the spectral line distinction between spots, plages, faculae and the presence of an exoplanet [57,86,87] underlie the effectiveness of our detection framework.This suite of characteristics can therefore be utilized to increase the robustness and fidelity of multi-feature detections in observations.

V. CONCLUSION
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 signalto-noise (S/N) ratio 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 radial velocity of 10 ms −1 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 fidelity sufficient to discern stellar features themselves from planets.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. of q the data is said to be a mono-fractal.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 [80].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 log 10 F 2 (s) changes slope with respect to log 10 s.We apply this method to each wavelength in the spectrum, giving us the dominant timescales of the complete spectrum.

Appendix B: Continuum Normalization
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.Fig. 5 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 Fig. 5 are in excellent agreement with those in Fig. 3.
FIG. 1. (a)Plan view of a planetary system orbiting the system's center of mass, i.e., both the planet and the star gravitationally interact, resulting in a Doppler effect from the light being observed from the star.(b) A spot on the stellar surface, as the star rotates about its axis.

FIG. 3 .
FIG. 3. Contour plots for the ratio Ω from Eqn. 1 and the corresponding fluctuation functions (no-noise and ρ = 200, 000) for two stellar features (a-d) and a stellar companion (e-f).(a, b) A spot with 1% coverage of the stellar disk; (c) A facula with 1% coverage of the stellar disk; (e) Planet with a Radial Velocity of 10 ms −1 .Whereas PCA[55]  requires a S/N ∼ 150 (per resolution element or pixel) to detect a spot of this size, our method can detect it for a S/N ∼ 15 (per pixel).This also shows how utilizing the statistical information of noise in the data rather than filtering it out refines the detection of such features.In particular, the second moment of the fluctuation functions show that the several orders of magnitude decrease (in the fluctuations) of spots to faculae to planets corresponds to an increased sensitivity to noise, thereby requiring a higher S/N for detection.Moreover, the inter-wavelength spread of these fluctuations as well as wavelength dependent fluctuations of faculae and planets can be used in differentiating between these features.The red dashed lines in (a,c,e) are lines of equal photon flux.These closely parallel to the contour lines of Ω, showing that MF-TW-DFA captures information at low-resolutions, which has important implications for studying existing and future datasets.

FIG. 4 .
FIG. 4. Contour plots for NP C computed as in Davis et al. [55] and FP C from Equation 2 for a spot (a,b), a facula (c,d) and a planet (e,f).Comparing the values of FP C here with Ω in Figure3shows that MF-TW-DFA provides information on a scale that is at least an order of magnitude finer than that provided by PCA and thus has a clear advantage in capturing stellar features and nonlinear dynamics.The dotted red line in the FP C plots in (d) and (f) is the contour line from the corresponding NP C plots in (c) and (e).Thus, while NP C gives a whole number, FP C provides a fine structure in the same data range.

FIG. 5 .
FIG. 5.The fluctuation functions (no-noise and ρ = 200, 000) for two stellar features (a, b) and a stellar companion (c).(a)A spot with 1% coverage of the stellar disk; (b) A facula with 1% coverage of the stellar disk; (c) Planet with a Radial Velocity of 10 ms −1 .Continuum normalization was performed on the data to account for changes in observational conditions.The normalization above was performed as the final step, after changing the resolution and addition of noise.In real world observations, normalization would be the final step in data processing as in (a, b, c) and these fluctuation functions are in excellent agreement with those in Fig.3.