Synthetic Light Curve Design for Pulsating Binary Stars to Compare the Efficiency in the Detection of Periodicities

1 Instituto de Astrofísica de La Plata, La Plata B1900FWA, Argentina; lydia@fcaglp.fcaglp.unlp.edu.ar 2 Instituto de Estadística, Universidad de Valparaíso, Valparaíso 2340000, Chile; gunther.avila@postgrado.uv.cl (G.F.A.M.); alejandra.christen@uv.cl (A.C.) 3 Facultad de Ciencias Astronómicas y Geofísicas, La Plata B1900FWA, Argentina * Correspondence: aldialb@fcaglp.unlp.edu.ar † Current address: Observatorio Astronómico, Paseo del Bosque, La Plata B1900FWA, Argentina.


Introduction
Observational surveys of massive stars of spectral type O show that most are forming binary systems [1][2][3]. Studies carried out in the region of 30 Dor also show that the binarity frequencies among dwarf and giant B-type stars agree with those of the O-type stars [4,5]. However, in contrast with the high binary incidence among dwarf massive stars, binary systems with B supergiant components are very rare. Barbá et al. [6] reported a noticeable sharp drop in the number of spectroscopic binaries between the supergiants of spectral types O9.7 and B0. These authors attributed this peculiarity to observational bias but suggested also discussing the action of possible evolutive scenarios such as binary mergers, binary disruptions, etc., to explain this effect.
On the other hand, massive stars pulsate in regular and quasi-regular (strange modes) oscillations (cf. [7]). They may show β Cephei type modes (low-order p and g modes) and Slowly Pulsating B (SPB) type variability (high-order g modes) [8,9]. Strange modes might also be excited with long periods, of the order of hours or days [10][11][12][13]. The pulsation mechanism is, in general, well understood (e.g., the κ mechanism, mechanism, or strange-mode instabilities) but for other groups, such as the γ Dor stars, roAp stars, and S Dor stars, no consensus has yet been reached [14]. These pulsations are revealed by photometric variations that can be identified by asteroseismological methods, while in the stellar spectrum, they generate variability in the broadening of the line profiles [15,16]. However, when pulsating stars are found in binary systems, their luminosity variations are modulated by the orbital motion [17]. This is the light-time effect, equivalent to a periodic Doppler shift of the pulsation frequency. Therefore, these effects generate complex patterns in the light curves that could make it difficult to detect or estimate periods or even conceal a binary detection.
The importance behind the detection of pulsating stars in massive binary systems lies in their impact on the study of the mass discrepancy problem, the theory of stellar interiors, and the verification of stellar evolution models [18], as they allow the masses of each component to be measured accurately and independently.
The search for periodicities is one of the most basic tasks in time series data analysis. It is a crucial feature for classifying variable stars and deriving stellar parameters. It is usually solved using an estimation function called a periodogram. Various types of periodograms are used in practice; for example, the Lomb-Scargle periodogram [19,20] is a known algorithm for detecting periodic signals in unevenly and regularly sampled time series. When periodic or quasi-periodic fluctuations arise in an interrupted or transient event [21], the wavelet transform is well suited for detecting changes in the parameters of signals due to its characteristic of focusing on a limited time interval of the data [22]. Nevertheless, the wavelet analysis is appropriate for time series with regular sampling. Today this problem has been partially addressed since satellite missions (Kepler, TESS, among others) provide photometric data of pulsating stars and binary systems with almost fixed cadences of a few seconds or minutes. Hence, using techniques for regularly sampled data is pertinent in this context [21,[23][24][25].
On the other hand, it is well-known that power spectra can often be misleading. The performance of many different methods for period detection and estimation in light curves has been tested in the literature, for example by Graham et al. [26]. Those authors examined a huge number of sources from three datasets (CRTS, MACHO, ASAS) and ran the data through eleven different algorithms such as the Fourier transform, Lomb-Scargle, phase dispersion minimisation, among others. The study of different observed light curves showed that, at best, period-finding algorithms could retrieve the period of a regularly periodic object with a reasonable degree of accuracy in only about 50% of cases. However, the situation is much worse for stellar objects with semi-periodic, quasi-periodic or multi-periodic variability; typically, only around 10-20% of the cases are the periods successfully recovered, evidencing the low level of accuracy for most of the methods, as shown in Figure 19 given by Graham et al. [26]. Although the authors performed an excellent comparative study among many methods for period searching, it still has the disadvantage that the authentic periods of the signal could be miss-detected because, even when they are accurate, they might be inexact. Based on this weakness, it is necessary to explore simulation designs for which the true values of the periods are known but compare the performance of a method not analysed in the above work, such as wavelet analysis, with the revisited Fourier analysis. This is the primary goal of the work.
In addition, one of the usual drawbacks when analysing observed data is the noise level of the signal, which generates consequently noisy periodograms. Thus, specific peaks in a periodogram might be spurious and not due to the presence of any real periodic phenomena [27,28]. Other drawbacks can also arise, for example, when the sample data are not evenly spaced in time. This situation is very frequent in the case of astronomical observations, where the data are out of phase in time due to various natural phenomena. It is well known that these non-equispaced data produce the so-called alias [29], i.e., false periods, where it is difficult to discern which peak of the periodogram is the real one and which is its alias. Unfortunately, simple and rigorous methods are unavailable to solve the alias discrimination task. In practice, we need to have prior information about the phenomenon under investigation to discuss the number and reliability of the possible periods. For example, for binary systems, the half period is often the most significant peak in a periodogram [26].
In summary, for a thorough analysis of the light curves of massive stars, it is necessary to have (i) a valid mathematical model to describe the time series, (ii) a period analysis tool that supports the detection of variable frequencies, and (iii) a homogeneous time distribution of the data to ensure a stable solution. Although the Fourier transformation provides a confident signal frequency spectrum, we often lose information about when this frequency happens. An alternative method is the wavelet analysis [30][31][32] that allows searching for local periods in a given window of time to obtain both frequency and time resolution.
In this work, we proposed to use a simulation design developed to study how different periodic or quasi-periodic simulated phenomena are reproduced in a wavelet analysis scalogram and a Fourier periodogram. For comparison purposes, we generated multiperiodic synthetic light curves (with pre-fixed periods) for a binary system with a pulsating companion. These synthetic models, unlike others, consider a statistical noise, a gap, and a slightly irregular sampling. In addition, phenomena such as beating and light time effects were also emulated. We obtained surprising and interesting results that will be useful in future studies of the periodicities of light curves and which have possible implications for the low incidence rate of binary systems among B supergiants. Section 2 introduces the simulation design by detailing the variables and their features. Section 3 briefly describes the time series analysis methods used. Section 4 shows the results obtained from the wavelet and Fourier analyses. Section 5 discusses our results in the context of the expected objectives and future work. Finally, Appendix A shows the complete frequency list for each model and all the periodograms resulting from the Fourier analysis. Appendix B shows the list of periods detected for the synthetic light curves with the wavelet analysis and the power spectrum for each one.

Simulation Design
In statistics, the area of experimental designs is concerned with finding the design that best helps to determine and study the factors that influence a given variable. There are many different types of designs: the factorial, fractional factorial, block, Latin square, and nested designs, among others [33,34]. Factorial designs are those in which experiments corresponding to all possible combinations of the "levels" or values of the factors or measurements of interest are carried out. The 2 k designs are factorial designs with k factors or variables for which only two levels or values are considered for each. These designs with few values for each factor are often used when it is desired to determine which factors are affecting a certain process or phenomenon. Extreme values of quantitative variables (such as temperature or concentration) or extreme categories of qualitative variables (such as the presence-absence of an attribute, among other possibilities) are considered for the levels. This type of design makes it possible to determine which factors or variables really influence the change or variation of the process or phenomenon. They are used in the first stage of the study when many variables are known about a process, and it is desired to determine its most influential factors. Once the relevant factors have been determined, it is possible to move on to a more exhaustive design, which considers more levels or values for each factor or even a more complex design. Each power-of-2 design considers 2 k different runs, which can be a considerable amount depending on the number k. Because it is not always feasible to perform 2 k experiments or runs, one can consider fractional factorial designs (e.g., 2 k−p ), which consist of running only a fraction of the total amount of runs given by 2 k /2 p . Therefore, it is important to optimally select this fraction. For this purpose, the fraction of the design is chosen in a way that allows obtaining a complete factorial design when p factors are discarded, for instance, because they are considered negligible after running the experiment. In this case, the number of factors p to be removed agrees with the power p in the denominator of 2 k /2 p . This methodology is based on the belief that if one considers a large number of variables, only a smaller number is relevant, and it does not produce significant variations in the phenomenon. In a later step, the remaining executions could be carried out to obtain a complete 2 k design.
One of the most important uses of fractional factorial experiments is in factorscreening 1 [33,34]. Furthermore, one of the benefits of using a statistical methodology is to be able to have hypothesis tests that will give us ways to quantify the significance of the findings. In this work, a fractional factorial design was used to establish the models to be simulated and underlined through Fourier and wavelet analysis. The design is detailed in the following section.

Designs of Synthetic Light Curves
In this stage, instead of performing a real experiment, we considered using numerical simulations of light curves to know in advance the result that we should find and thus be able to make a diagnosis of the goodness of each of the methods, Fourier analysis and wavelet analysis, to be used to detect periodicity. In the context of this research, we performed the simulations of synthetic light curves of a pulsating star in an eclipsing binary system. To build them, we used the photometric data of a known binary, HD 19,356, obtained by the TESS (Transiting Exoplanet Survey Satellite) mission [35] with a temporal resolution of 120 s and a gap of about 4 days. Then, to obtain the synthetic light curve, we modelled the data with the software PHOEBE (PHysics Of Eclipsing BinariEs) [36,37], adopting the orbital period of 2.87 days corresponding to the selected binary system [38]. The resulting synthetic curve, referred to as model No. 0, is shown in Figure 1. This model will be the basis for the complete simulation design.  As the time series of many physical processes show stochastic and auto-correlated properties [39,40], we contaminated the synthetic curves with an ARMA (1, 1) noise [41]. Simulations of light curves with ARMA noise have already been used by authors such as Caceres et al. [40,42] to study different phenomena, such as exoplanetary transits. According to these authors, autoregressive modelling is more effective after minimising systematic variations due to instrumental or atmospheric conditions. Alberici [37] performed numerical simulations of light curves with three different types of stochastic processes to model noise: two stationary, a Gaussian process (white noise), an ARMA (1, 1) process (dependent of the past), and a non-stationary GARCH (1, 1) process (also dependent of the past). This author independently tested the Fourier and wavelet analysis and found no difference among the diverse types of noise models when they were analyzed with the same techniques.
Therefore, according to Alberici [37] and Caceres et al. [40,42], we opted for using an ARMA (1, 1) noise model, which has some dependence on the prior values, for the random fluctuation of the star luminosity [40]. The ARMA processes were used in this work to model the synthetic light curves, considered without systematic effects. The magnitude of the variability of the ARMA process was determined at a low value for simplicity, and it is left for future work to perform a study with higher variability of the noise process (i.e., where the variability of the luminosity could affect the detection of the binary star eclipse or pulsations, for example).
Therefore, a simulation design was planned with six variables and two categories (factors) each, defined as follows: 1.
Beating is the periodic variation in amplitude at a certain wave point due to the superposition of two waves having slightly different frequencies f i and f j , producing a time disturbance between constructive and destructive interference. For simplicity, we considered two waves of equal amplitude, A, and wave number, k, but with slightly different angular frequencies: The resulting wave is the sum of y i + y j , which gives where Simulations with and without beating phenomena were considered.
In this variable, we explored the effects on the light curve considering regular pulsation cycles with a period longer or shorter than the binary's orbital period.

3.
Amplitude of the pulsation with respect to the depth of the primary eclipse (X 3 ).
Two categories were considered for this variable. The first is that the pulsation amplitude is 20% of the depth of the primary eclipse, which we call a large amplitude. The second is that the pulsation amplitude is 4% of the depth of the primary eclipse, referred to as a small amplitude. The values of the amplitudes of the pulsations were chosen in this way because they agree with those observed in TESS light curves of some pulsating binary stars [43].
Scenarios where the light curve is affected by simultaneously one or two pulsation periods were considered.
This scenario considers a star in a binary system that pulsates sinusoidally with a single frequency. Its luminosity varies with time as a consequence of this pulsation. However, the orbital motion of the star leads to a periodic variation of the distance between the observer and the star, so the phase of the observed luminosity variation also varies with the orbital period [17]. The light-time effect leads to frequency multiplets in the Fourier spectra for pulsating binary stars. To model this effect in our synthetic light curves, we relied on the simplest case of a pulsating binary star with a circular orbital motion. Furthermore, to emulate the splitting of the pulsation period, we proceeded to affect the curve with three sinusoidal pulsations: the intrinsic pulsation frequency of the star and two additional oscillation frequencies due to the light-time effect, where the spacing between each of them and the intrinsic frequency of the pulsation is equal to the orbital frequency (for more details see [17]).
Let P o be the orbital period of the star and P p the period of the pulsation, then the multiple components of the luminosity variations will have the following frequencies: Furthermore, the amplitudes and phases of the multiple components can be used to derive meaningful information about the orbital elements and the radial velocity curve. Scenarios with and without the light-time effect were considered in this work.
Let us now consider a set of harmonic oscillators. Periodic motion means that there must be some time τ at which all the different oscillators have completed an integer number of cycles, i.e.: This means that all n k are integers; therefore, the quotient ω i /ω j must be a rational number. If, by contrast, it is irrational, any oscillator i or j might not complete an integer number of cycles in the time τ, giving rise to quasi-periodic pulsations (QPP). This kind of behaviour is a common feature in the solar atmosphere and flare stars where the flaring energy is released by magnetic reconnection (cf. [44]). Quasi-periodic pulsations are also observed in the α Cyg variables and are attributed to strange-mode oscillations [45,46]. A star oscillating with some periodic phenomenon would appear in a power spectrum as a peak at precisely that frequency (Dirac-δ function). Otherwise, it will appear as a broader (sometimes Lorentz-shaped) peak.
Simulations with and without quasi-regular periods were considered. Since we have defined X i variables (with i = 1, · · · , 6) and each has two levels or categories, the total design is 2 6 = 64 models. Since the design is optimal, we can consider only a fraction of the total design without losing information about the variables involved. Thus, in this work, we analysed 2 6 /2 1 = 32 models. The design parameter model is defined in Table 1. Table 1. Simulation design features. The table lists the phenomena, properties, and periods fixed to simulate the pulsations in each model. The light-time effect and the beating phenomenon generate periods that are labeled with the indices "t" and "b", respectively. In all cases, these effects are present only in the lowest pulsation period.

Time Series Analysis Methods
Signals were differentiated into stationary and non-stationary. The first ones were localised in the time since their frequencies did not vary. The spectral analysis of this type of signal was carried out using the Fourier transform, which allows its decomposition into infinite sinusoidal terms, transforming the signal from the time base to the frequency base and back again. When switching from the frequency domain to the time domain, valuable information is lost, which, due to the stationary nature of the signal, is irrelevant. However, in non-stationary signals, i.e., those whose frequencies vary with time, the loss of information becomes very relevant, and it is difficult to determine when a frequency change occurs [47]. This is why wavelet analysis is presented as a tool to obtain a detailed decomposition and reconstruction of signals with abrupt changes, using a multi-resolution analysis with windows of a variable length adapted to the change in frequency of the signal. This is the main property of the wavelet, as it minimises the computational cost of the analysis by not needing infinite terms. Moreover, it improves the detail of the detection by allowing the use of long time intervals in those segments where greater precision is required at low frequencies and smaller intervals where information is required at high frequencies. Unlike Fourier, where the basis functions are sines and cosines of infinite length, in the wavelet analysis, the basis functions are localised functions in frequency and time [48]. Thus, they are a suitable method for studying quasi-regular, transient and discontinuous phenomena, such as light curves of binary stars.
The software used to carry out data analysis are detailed below.

The Period04 Tool
The Period04 tool [49] is a hybrid program, written in Java and C++, dedicated to the statistical analysis of large astronomical time series which contain gaps or missing data over a significant time interval. The program allows the extraction of individual frequencies from the multi-periodic content of the time series through an analysis based on the Discrete Fourier Transform (DFT). The light curves were fitted using the formula: where Z is the zero-point, A i is the amplitude, f i is the frequency, and φ i is the corresponding phase.

The WaveletComp Package
WaveletComp is an R package [50][51][52] for continuous time series analysis based on the Morlet wavelet family [32]. This family of wavelets leads to a wavelet transform of the time series. The transform can be separated into real and imaginary parts, thus providing information about the local amplitude and instantaneous phase of any periodic process over time.
The Morlet "mother" wavelet [32], in the version implemented by WaveletComp, is: In this frame, the angular frequency ω is set to 6 [51], therefore the relation between the wavelet scale, a obtained with wavelet analysis, and the "equivalent Fourier period" [53] is given by: This correction is given directly by the WaveletComp package. The Morlet wavelet transform of a time series f (t) is defined as the convolution of the series with a set of "daughter wavelets" generated by the mother wavelet through translations in time, τ, and dilatations in frequency, a, being: where φ * denotes the complex conjugate function. The local amplitude of any periodic component of the series can be expressed as: Then, the square of the amplitude can be interpreted as the energy density of the wave in the time-frequency domain and is called the wavelet power spectrum: This power spectrum is represented by a scalogram, see Figure 2, where a color code, from blue (low power) to red (high power), is used to indicate the power spectrum range. Maximum probability values are represented with black lines. White lines delimit probability regions with a 0.1 significance level against a white noise null. The WaveletComp package returns by default a scalogram whose axes are period vs index (ordinate and abscissa, respectively), where index refers to the amount of data in the time series and is directly related to the time axis of the series. More information can be obtained from the average wavelet power over time, which will depend on how many cycles a certain period completes in the series and not only on the magnitude of the wavelet power. The red curve represents values with a significance level less than or equal to 0.05, all of which differ from white noise. The values with the highest average powers indicate the most relevant periods for the signal. A practical example of wavelet analysis for multi-periodic time series with periods 30, 38, and 80 is shown in Figure 2. Periods 30 and 38 generate a beating effect in the series due to their proximity. All periods with higher power are detected with solid black lines in the scalogram, and the significance decreases towards nearby periods. The color code on the right-hand sidebar indicates a higher significance for the period 30 and decreases until the period 80 (see central panel). This is due to the number of cycles each period completes in the time range of the series. In the right panel, on the average wavelet power plot, periods 30 and 38 are more difficult to identify individually, unlike period 80, although they show a striking feature in the scalogram over the entire timeline identifying the presence of the beating. The results show the mixing effect produced in power obtained with the wavelet analysis for nearby periods. In agreement with the power spectrum, the maximum average wavelet power is around period 30 and decreases for longer periods. Estimating the spectral power using the continuous wavelet transform of a nonperiodic time series can generate false spectral structures at the beginning and the end of the wavelet spectrum [54]. This is a known problem in the periodogram analysis caused by replacing the integrals over an infinite continuous signal by finite sums [55]. The "zero padding" method is used to solve this problem, which completes with zero values on both edges of the time series after a sample mean is removed. Then, the wavelet spectrum is estimated [53,54]. One of the drawbacks of applying this method is that it produces discontinuities at the edges of the wavelet spectrum. The region of the wavelet spectrum where padding effects become important is known as the cone of influence (COI) [53]. Therefore, the portions of the wavelet spectrum between the time axis and the COI should be considered unreliable because statistical significance tests in that region cannot be trusted [53,56]. Additionally, wavelet analysis performed with the WaveletComp package can only produce accurate results if the time series is equispaced.

Period Detection Criteria
To identify the periods detected with both tools in the synthetic light curves, the following selection criteria were established: • For the wavelet analysis, each light-curve model was separated into two parts because they present a discontinuity (or a gap) in the time sequence (see Figure 3). For all periods that showed a maximum average wavelet power greater than 10, the mean value of each period across time and its standard deviation were calculated as suggested by Roesch and Schmidbauer [51]. The latter was used as a measurement uncertainty of the period. They were considered more or less significant according to the color code of the scalogram, even those outside the cone of influence, since there is a significant improvement in the reconstruction of the light curves when considering those periods as well [37]. The period's harmonics were identified for all independent periods. For the harmonics of the orbital period, we searched for the values included in the intervals with endpoints given by P h ± FWHM/2 (see Table A2), where P h is an harmonic in the Model No. 0 and FWHM/2 = 2 log(2) * σ P h , being FWHM the full width at half maximum of a normal distribution with standard deviation σ P h at any expectation value. • In the Fourier analysis, frequencies with a signal/noise ratio ≥ 5 were selected [57]. Then, the uncertainties of the frequencies were calculated following a nonlinear least squares fitting procedure available in the software [49]. According to Bognár, Zs. et al. [58] and Lenz and Breger [49], we accepted a peak as a combination of frequencies if the amplitudes of the main frequencies were greater than that of their presumed combination term and the difference between the observed and predicted frequency was not greater than the Rayleigh resolution criterion of the data sample. Once the independent frequencies were selected, we reported the corresponding periods.

Results
Each synthetic light curve generated with the design model of Table 1 is shown in Figure 3. These curves were analysed using the wavelet and Fourier analysis. To give a detailed description of the results obtained and how they are interpreted, we selected a group of four models with specific characteristics in common and compared them with Model No. 0 (shown in Figure 1). For the remaining 28 models, we will comment on the general results obtained, following the same analysis criteria. The results for the complete set of models are in Appendices A and B. Table 2 lists the independent periods detected for Models No. 0, 1, 5, 17, and 21 with each method and their corresponding errors. Among these five models, the magnitude of the pulsation amplitude, the presence or absence of the beating phenomenon, and the values of the pulsation periods vary. The respective wavelet power spectrum, average wavelet power, and periodograms are illustrated in Figures 4 and 5.   Table 1. · · · 2.859 ± 0.001 · · · 1.106 ± 0.002 No. 21 4 · · · 7.924 ± 0.062 2.803 ± 0.053 · · · 10.144 ± 0.018 · · · 2.878 ± 0.004 1.103 ± 0.001 For Model No. 0, the wavelet and Fourier analysis recognises the binary's orbital period as a fundamental frequency and finds its harmonics. The wavelet scalogram for both parts of the curve shows a structure (seen as "faces") at the position of the minima of the light curve, also reproducing the eclipses of the binary (compare the timeline and structures of Figures 1 and 4). However, the period detected by the wavelet tool is systematically slightly shorter than the original orbital period. This shift could be related to the fact that the time scale, adopted from a TESS light curve observation, is not completely equispaced.
When examining the wavelet power diagrams for Model No. 1, we noticed some differences from Model No. 0, particularly in the features that identify the eclipses. For example, it loses detail in finding harmonics with periods 1.43 and 0.95 days, maybe due to the size of pulsation amplitude (20%). Even so, it finds the binary 2.87 days period with an error of 0.09, although in the second part of the curve and with much lower power than in Model No. 0. Additionally, a period of about 1.15 days is found, and its trace shows a sinusoidal modulation in the scalogram. Although it could be a spurious value, we associated this period as an independent value, likely with a harmonic of 8.96 or 10 days (8 × 1.15 = 9.2), where the value 8.96 days were artificially introduced to obtain the beating. However, this shift could be due to three near pulsation periods: 3e, 8.96, and 10 days. Furthermore, a period of 7.852 days is also present with high power and could be related to the 3e days period. Even though the wavelet analysis does not find those exact values, the software likely detects the presence of significant periods in that range.
Model No. 5 is affected by the same pulsation periods as Model No. 1 but with a smaller amplitude. Therefore, the model's scalogram and power diagram show a lower significance, but the same periods and scalogram features are maintained. As shown in Figure 2, the wavelet analysis performed with the WaveletComp package presents complications in differentiating nearby periods, particularly cases with a beating, such as Models No. 1 and 5, resulting in an intermediate value.  When analysing Model No. 17, we observed that the structures present in the scalogram are very similar to those of Model No. 1, both with a pulsation amplitude of 20% with respect to the depth of the primary eclipse. Since they share two periods, the 7.9 days period and the harmonic of 1.15 days, they could be related to the periods of 3e and 10 days.
Finally, the scalogram and average wavelet power of Model No. 21 is very similar to that found in Model No. 5 since they share two pulsations with amplitudes of 4%. However, in this case, it detects the binary period and another with a 7.92 days period, but with lower power than in Model No. 5.
The periodograms resulting from the Fourier analysis show remarkable similarities between the models affected by low-amplitude pulsations regardless of the periods involved. This is particularly seen in models No. 5 and 21. At the same time, slight variations are observed with respect to model No. 0. The pattern changes markedly for models No. 17 and 21, as seen in Figure 5). A similar situation occurs with models No. 6-8.
Quasi-periodic oscillations decrease the power of detection of the Fourier Transform method, tending not to detect them. The Period04 tool accurately sees the fundamental frequency of the binary orbital motion (P orb = 2.87 ± 0.01 days) and its harmonics in all cases. However, it observes two pulsation periods with a significant error close to P = 10 days, detected as 9.26 days and 12.9 days in Model No. 1 and 9.47 days in Model No. 5. The 9.26 days (or 9.47 days) period could be likely related to the beating phenomenon with P = 8.96 or the pulsation with P = 10 days, or a mix of both. This is supported because, for example, the period 9.26 days also appears in Model No. 2 (see Appendix A) that has regular oscillations (with P = 7 days and P = 10 days, and a beating of 7.7 days), but in this case, the period 7 days is inexactly seen as P = 6.657 days (see also On the other hand, the wavelet analysis found the binary's period satisfactorily, although it underestimated the value with an average of 2.81 days (of 31 models), which suggests a period shift of 0.06 days. Concerning the rest of the periods, this method found an equal or smaller number of periods than the Fourier analysis. Although the mean value has a systematic error, its variance is similar to that obtained with the Fourier method. It detected a few periods, mainly in models with low pulsation amplitude and where the period values coincide with a harmonic of the binary's orbital period. However, it found all periods in four models (Models No. 17, 18, 26 and 28). These have a high pulsation amplitude (20%) and absence of beating phenomenon.
A comparison of the methods is shown in Figure 6. The boxplots of the hit ratios of all the prefixed periods obtained by wavelet and Fourier analysis are in the left and right panels, respectively. We see that the median percentages (thick black line inside the box) are much lower for the wavelet analysis (<0.4) than for the Fourier analysis (about 0.8). Figure 7 presents a boxplot comparison for both techniques, of the hit ratio segmented by the categories of five of the variables included in the simulation design: beating phenomenon (X 1 ), period of the pulsation (X 2 ), amplitude of the pulsation (X 3 ), number of pulsations periods (X 4 ), quasi-regular pulsations (X 6 ), calculated with a wavelet and Fourier analysis, respectively. Since we were dealing with the half-model design, the result does not provide information about one variable.
Regarding the results obtained for the wavelet analysis, we can observe that the variable that influences its performance the most is the pulsation amplitude. At low amplitude, it is difficult to detect the intrinsic period of the model. The presence of the beating phenomenon and the fact that the period of the pulsation is quasi-regular also have some influence. The hit ratio increases when the pulsation period is greater than the orbital period and when there is more than one pulsation (see Figure 7). The analysis of the variance of the percentage of hits of a quasi-periodic oscillation with a large pulsation amplitude (i.e., 20%) finds the period only the 55% of the time. While for a regular oscillation, it detects the period 85% of the time. Instead, if the pulsation amplitude is low, it finds the period about 30% of the time in both cases. Quasi-regular or regular period detection is also affected by their interaction with the orbital period, depending if they are shorter or longer than it. The regular period detection rate varies from 60% (for shorter periods) to 50% (for longer periods) and for quasi-regular periods from 65% to 30%.
Unlike wavelet analysis, Fourier analysis is unaffected by the pulsation amplitude. The proportion of detected periods (an approximate value) increases if the period pulsation is longer than the orbital one, but they are more accurate when it is shorter (see Figure 7). The presence of beating and quasi-period decreases the detection capability a little. In both cases, the medians (thick black line in the box) is about 75%, which changes to 100% without the phenomenon. From an analysis of the variance of the hit ratio, we obtain that if the pulsation period is shorter than the orbital period, the presence of the beat phenomenon improves period detection (70% vs 80%). Otherwise, if the pulsation period is longer, the presence of this effect worsens it (95% vs 75%). Quasi-regular or regular period detection is also affected by their interaction with the orbital period, depending if they are shorter or longer than it. The period detection rate for regular oscillations varies from 100% (for shorter periods) to 80% (for longer periods) and for quasi-regular periods from 50% to 85%. For pulsation periods longer than orbital periods, the performance is similar.

Discussion and Conclusions
It is widely known that power spectra can be misleading and frequently lead to inaccurate conclusions. For this reason, many works present a comparative analysis of popular algorithms applicable to period searching. Sometimes the study of one particular light curve with various numerical methods resulted in a high dispersion or undetected periods as reported by Graham et al. [26] and Distefano et al. [59].
The objective of the current work was to evaluate the response of two methods, based on the Morlet wavelet function and continuous Fourier transform, in the production of bias, alias, or shifts in period detection using a simulation design with pre-fixed and known periods. For this purpose, we performed a simulation design of synthetic light curves of a pulsating star in an eclipsing binary system with an orbital period of 2.87 days. It is a 2 k−p -statistical design with k = 6 variables, in which the beating phenomenon and the light-time effect were considered. Of the total design method, only one-half was executed. The 32 time series models were analysed using the WaveletComp and Period04 tools. Although many time series analysis tools are available, the application of such a large simulation design to more than two techniques would be unattainable.
The obtained synthetic light curves, shown in Figure 3, present similarities with the light curves of pulsating binary stars obtained with the TESS mission (see Figure 8)   Comparing both applied methods, we observed that the orbital period and its harmonics were well detected in most models. Even though the Fourier method provided more accurate period detection, the wavelet analysis found it more times. When estimating the orbital period, the WaveletComp package trends to provide a low period value or period shift of approximately 0.06 days in all the models, but this has a similar dispersion to the orbital period obtained from the Fourier analysis. In agreement with this finding, Foster [22] also reported through numerical simulations that for a simple sinusoid, the maximum frequency obtained with wavelet analysis has an offset with respect to the signal frequency, which can be compensated. This author also commented that under irregular sampling, the detection of the period by the wavelet analysis varies with time, making it difficult to give an exact value of the same, especially with the increase in missing data in the sampling, which can reach up to 20% at high frequencies or find spurious periods at low frequencies. Although we considered "almost" a regular sampling, it may be that some of the still sampling irregularity present in the timeline scale selected from TESS observations has affected the performance of the WaveletComp package.
As shown in the previous section, pulsation amplitude is a very powerful variable in period detection. In the case of the wavelet analysis, more satisfactory results are found recovering the periods of the pulsations with large amplitudes. On the other hand, the Fourier analysis did not show difficulties in recognising the periods when the amplitude is low, but they were detected with significantly lower power than models with large amplitudes.
Neither tool was able to distinguish accurately the periods leading to a beating phenomenon when they were longer than the orbital period, resulting in both cases in an intermediate value. However, when the beating was formed by periodic pulsations shorter than the orbital period, the Fourier analysis found it in all cases, but quasi-regular periods remained unresolved. Concerning the beating period √ 2, it is necessary to remember that this period coincided (unfortunately, but possible in a real scenario) with the first harmonic of the selected orbital period of 1.43 days. Therefore, when a value of 1.4 days was detected, it was considered the harmonic of the orbital period. If the values set in the simulation design were unknown, it would be impossible to distinguish between them.
The same holds for the 0.94 days period. By analogy, wavelet analysis does not detect these periods either.
Finally, the wavelet power spectrum revealed in all models a particular structure located at the position of the primary eclipses of the binary in the synthetic light curve, even showing the instant of time in which they occur. Then, if a similar structure is found in a real data set, it could indicate binarity. In addition, as mentioned before, several harmonics of the orbital period were detected by both tools, which is also an indicator of binarity.
In general, the Fourier analysis by means Period04 tool had a higher proportion of periods correctly and accurately detected than the analysis performed with the WaveletComp package. However, the discrepancy could be attributed to the less information available in the latter when performing the period analysis since splitting the light curve into two parts was necessary. Despite the high cadence, the TESS data outside the large gap also present a slightly irregular sampling, which affects the analysis. Furthermore, this package uses a somewhat coarse resolution level, preventing it from accurately determining the periods encountered. Foster [22] also found some kind of asymmetry in period detections, with more values at low frequencies and suppression at high frequencies.
Considering the period hit rate analysed using synthetic light curves for eclipsing binaries, a bias in detecting some pulsation periods (beating phenomenon and light time effects) in multi-periodic light curves would be expected. The hit rate depends on the amplitude and frequency of the pulsation with respect to eclipse depth and orbital period. The situation would be more unfavourable for purely elliptical double stars, partial eclipses, or systems influenced by the gravitational attraction of a third body that might behave similarly to a sinusoidal pulsation. In some cases, only a light modulation could likely be acting because the binary orbit is not aligned with the observer's direction. In these cases, the period hit rate would be lower, masking the binary detection, for example, among B supergiants.
Although the recovery of all pre-fixed periods was not feasible using WaveletComp and Period04, valuable insights were obtained regarding how each technique responds to the various modelled situations. For instance, we observed that WaveletComp can only detect an intermediate value when confronted with two close periods due to the poor spectral resolution of this package. Additionally, we have noted that Period04 fails to recognise quasi-regular periods, but has good accuracy in detecting high-frequency oscillations. Furthermore, we found with both tools structures that indicate binarity (in the periodograms and scalograms), information that will be very useful when performing these analyses on the light curves of binary star candidates, among other findings.
Considering the advantages and disadvantages of the techniques analysed, we recommend using at least two diagnosis tools to conduct a detailed analysis of time series data to obtain confident results. Additionally, a second criterion should be considered to test the accuracy of the found frequencies or periods since the methods might find a harmonic, inaccurate period, or spurious value. Particularly, the wavelet method might predict only a harmonic of the true period. Since approximated values were detected with at least one method, a fine-tuning of trial periods around a certain percentage of its value, by applying phase diagrams, would help to recover accurate values.
Evaluating the behaviour of other tools based on wavelets that do not have this weakness (e.g., Weighted wavelet Z-transform (WWZ) [22]) is one of the goals to follow. This method also enables an improvement in the resolution of the period searches, allowing a high degree of accuracy. In addition, considering even more realistic models or observed light curves that present a gap and time-varying periods would allow us to take full advantage of the main feature of wavelet analysis.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

TESS
Transiting Exoplanet Survey Satellite PHOEBE PHysics Of Eclipsing BinariEs

Appendix A. Fourier Analysis Results
For each synthetic light curve, the frequencies with a signal/noise ratio ≥ 5 and their errors, are listed in Table A1. The fixed listed periods of the design models are listed in column 2. It does not include the period of the binary system (2.87 days) present in all the models, except for model No. 0, where is the only period present. Columns 3 and 4 list the frequency and the corresponding period in days. Column 5 indicates the independent frequencies found and the linear relationship among them. The respective periodograms are shown in Figures A1 and A2. No.