The power of wavelets in analysis of transit and phase curves in presence of stellar variability and instrumental noise I. Method and validation

Stellar photometric variability and instrumental effects, like cosmic ray hits, data discontinuities, data leaks, instrument aging etc. cause difficulties in the characterization of exoplanets and have an impact on the accuracy and precision of the modelling and detectability of transits, occultations and phase curves. This paper aims to make an attempt to improve the transit, occultation and phase-curve modelling in the presence of strong stellar variability and instrumental noise. We invoke the wavelet-formulation to reach this goal. We explore the capabilities of the software package Transit and Light Curve Modeller (TLCM). It is able to perform a joint radial velocity and light curve fit or light curve fit only. It models the transit, occultation, beaming, ellipsoidal and reflection effects in the light curves (including the gravity darkening effect, too). The red-noise, the stellar variability and instrumental effects are modelled via wavelets. The wavelet-fit is constrained by prescribing that the final white noise level must be equal to the average of the uncertainties of the photometric data points. This helps to avoid the overfit and regularizes the noise model. The approach was tested by injecting synthetic light curves into Kepler's short cadence data and then modelling them. The method performs well over a certain signal-to-noise (S/N) ratio. In general a S/N ratio of 10 is needed to get good results but some parameters requires larger S/N, some others can be retrieved at lower S/Ns. We give limits in terms of signal-to-noise ratio for every studied system parameter which is needed to accurate parameter retrieval. The wavelet-approach is able to manage and to remove the impacts of data discontinuities, cosmic ray events, long-term stellar variability and instrument ageing, short term stellar variability and pulsation and flares among others. (...)


Introduction
The light curve of an exoplanetary system may show transits, occultations and phase-curve variations. The transit technique offers a unique opportunity to determine the accurate radii of transiting exoplanets. Complementing the photometric transit observations with radial velocity data, the planetary mass and mean density can be determined. The phase-curves describe the scattering and reflecting properties of an atmosphere at different orbital phase. Phase-curves and occultations are considered as the best opportunity to study the three-dimensional structure of planetary atmospheres (Parmentier & Crossfield 2018;Winn 2010).
In such transit-as well as in phase curve-analysis, the following four problems can arise. (i) The stellar activity, stellar variability -including pulsation and granulation, too -and instrumental noise cause difficulties to find and to restore the exact shape of the transit-and phase-curves (e.g. Oshagh 2018; Sulis et al. 2020). The transit depth can be also affected by stellar spots yielding wrong planetary radii. If sudden and discontinouos flux variations (jumps) occur in the flux measurements due to a cos-mic ray hit or other kind of instrumental effect, then there is a difficulty to establish the mean flux level of the host star. This can lead to further change in the transit depth because the normalized flux level is different and maybe not well fitted to each other before and after the jump. (ii) The beaming-effect might be degenerate with the reflection effect and care is needed to separate them from each other (Csizmadia 2020). (iii) The ellipsoidal effect, when it is significant, must be also separated from the phase-variation. The ellipsoidal and the reflection effects have higher order harmonics of the orbital frequency and if it is not modelled carefully, the badly modelled ellipsoidal effect can affect the shape of the reflection curve and this may lead to misconclusions. (iv) In addition, the exact shape of the phase curve is not known without a detailed a priori knowledge of the atmosphere (composition, scattering and reflecting properties, scale height, clouds, particle sizes of the aerosols etc, Garcia Munoz & Isaak 2015).
We make extensive tests on synthetic light curves to overcome problem (i), namely the stellar and instrumental noise sources are modelled by a wavelet-transform. We investigated Article number, page 1 of 11 arXiv:2108.11822v1 [astro-ph.EP] 26 Aug 2021 A&A proofs: manuscript no. pow how well the wavelet-method can filter out the stellar variability and instrumental noise effects. We show in the present study that the wavelet-transform is a powerful tool to model flux-variations of stellar and instrumental origin which increase the accuracy and precision of parameter retrieval. In a subsequent paper we apply our method to KELT-9b (Csizmadia et al, 2021, submitted, Paper II).
To solve problems (ii-iv) one way can be to use prescribed forms of phase-curves and to improve the description of the ellipsoidal effect. In Paper II we attempt to fit single cosinusoidal, Lambertian, Kane-Gelino-and Kopal-type phase curves to the time-series data of KELT-9b obtained by the TESS space telescope. These four different phase functions yield significantly different shaped reflection curves. In Section 2 we detail the model of the ellipsoidal, beaming and reflection effects used for the fine-analysis of the KELT-9b light curve in Paper II. We also update the gravity darkening model of TLCM in Section 3. The wavelet model and its test are presented on Section 4. The summary of this study and our conclusions can be found in Section 5.

Model of out-of transits variation
The out-of-transits variation is usually divided into components of reflection and ellipsoidal effects. When it became observable with space-based telescopes, this list was extended by the beaming-effect component (Zucker et al. 2007;Faigler & Mazeh 2011).
the sum of the ellipsoidal, beaming and reflection (phasecurve) effects. We describe our model of these effects hereafter. There are other kind of photometric variations which we consider part of our red noise model (see Section 4): stellar pulsation, stellar activity, possible additional eclipses caused by another star, instrumental and other -non-white noise-like -effects etc. Note that we distinguish between phase curve (only the reflection effect, without the beaming and ellipsoidal effects) and the phase function which contains the time-dependence of the phase curve.

Phase curve
We utilize Transit and Light Curve Modeller (TLCM, Csizmadia 2020) which expresses the phase curve variation (F ph ) in the following form: Here F ph and F star are the reflected and stellar fluxes, respectively. The phase-angle α in the phase function Φ is where ω is the argument of the periastron of the planetary orbit and v is the true anomaly. i p is the inclination of the planetary orbit. I is the surface specific intensity in the passband of the observation, d is the mutual star-planet actual distance. R denotes the radius. A geometric is the so-called wavelength-dependent geometric albedo. The angle ε takes into account that there can be a phaseshift in the phase curve due to atmospheric circulation, i.e. the brightest point of the planet can be shifted eastward or westward relative to the substellar point 1 (Parmentier & Crossfield 2018).
The observed values of ε vary between wide ranges, from -70 to +50 degrees (Parmentier et al. 2016;Bell et al. 2021) and references therein.
For the star the phase curve can be obtained in a similar way by interchanging the indices appropriately and shifting the phase curve by 180 • . This may be important for detached eclipsing binary stars where TLCM can also be used for modelling and takes all these effects of the two stellar component into account (Csizmadia 2020).

Time dependence of the phase function
The exact form of the phase-function Φ strongly depends on wavelength, the chemical composition of the atmosphere, the particle size in it, optical depth, clouds properties, single scattering albedo (e.g. Garcia Munoz & Isaak 2015). Establishing the exact form of the phase function needs a priori knowledge on the atmospheric properties which are not always available. In addition, it requires complex and lengthy numerical calculations (Garcia Munoz & Isaak 2015). For data analysis, one can try analytically expressed approximate formulae, too. In this paper series we probe the following four phase functions offered by TLCM on KELT-9b as follows.
First one is a simple cosine-like: The second one is a Lambertian: The third one is taken from Kane & Gelino (2011) where the phase angle must be measured in degrees: This formula is based on the observations of Venus and Jupiter and it takes into account that these planets have significant backward-scattering due to their clouds (Hilton 1992). On eccentric orbits the particle properties can change in the atmosphere due to the variable insolation. According to Kane & Gelino (2011), the geometric albedo in Eq. (1) must be replaced by A (d) for this case as Of course, one can ask how well this Kane-Gelino phase function can perform on hot Jupiters where the response of the atmosphere can be quite different than in the case of the cooler Jupiter and the terrestrial like venusian atmopshere. We try this phase function on the hot Jupiter KELT-9b in Paper II, too, and this will give an answer. The fourth and the last phase function tried here is based on the theory of binary star phase function which takes umbral and penumbral effects into account up to the fourth order of the phase-angle. This is taken from Kopal (1959): where and the limb darkening correction is given for the reflection in a linear form as C 3 is zero according to Kopal (1959). For the planet, one can take linear limb darkening coefficient as u 1 = 0 as first order approximation. When one calculate the planet's reflection effect, j = 2 where R 2 = R planet ; while if one calculates the star's reflection effect j = 1 with R 1 = R star . (In case of a detached binary system, these are the primary ( j = 1) and the secondary stars j = 2, respectively.) This phase function also takes back-warming effects into account which are negligible in star-planet system but it can be important for detached or even closer binary star system. The planet-to-star radius ratio (R planet /R star ) and the scaled semi-major axis ratio (a/R star ) are known from the transit and light curve analysis or -as in this study -fitted simultaneously with the phase-curve parameters.
In eccentric orbits the star-planet distance d varies as where e is eccentricity (not confuse with the Euler-number in Eq. 6).

Dayside and nightside emission
Clearly, the dayside and the nightside phase curves are varying in the opposite phase and therefore we have The total planetary phase-curve is the sum of the dayside and the nightside emission at a certain phase: Comparing Eq. (15) to Eq. (1) we can relate to the measured quantities to the parameters we search for via: and Note that I planet /I star , R planet /R star , A geometric and reciproc of R star /a are fitting parameters in TLCM and they can be measured from the simultaneous fit of transit, occultation and from phase-curve. We note that we assume the very same phase function for the reflected light and for the dayside/nightside emission in TLCM, which is, of course, just an approximation of reality. However, we use as simple approach as possible.

Ellipsoidal effect
The flux-variation caused by the ellipsoidal shapes of the components is characterized following Kopal (1959) ( j = 1 for the star and j = 2 for the planet): with the gravity darkening correction is Here hc/k B = 14388µm · K, effective wavelength of the observation λ is given in microns, T eff is in Kelvins, and according to Kopal (1959): Limb darkening coefficients u 1,2 of Kopal (1959) are from the transit fit. They are related to the limb darkening coefficients u a and u b of Claret (2004) via u 1 = u a + 2 · u b and u 2 = −u b . (Of course, the limb darkening coefficients can be different for the two objects in the system.) The apsidal motion constants k i are functions of stellar mass, metalicity, radius and evolutionary status. The k i s of the primary star or of the host star are fixed at their theoretically calculated values of Claret (2004). Note that the apsidal motion constant is half of the Love-number (Csizmadia et al. 2019). f 1 = 1 for the star and f 2 = I planet /I star (R planet /R star ) 2 for the planet (or secondary star) in the respective photometric passband (same as in Eq. 16).

Beaming effect
The beaming effect is characterized as Here B is the contribution of the beaming effect to the observed flux in units of stellar-flux, Ω j is the spectral index derived from theoretical stellar spectra of Munari et al. (2005), effective temperature, surface gravity log g, and metalicity Z. These spectra were convolved with the response function of used photometer (Csizmadia 2020). The star and the companion has the index j as before. K and c are the radial velocity amplitudes and speed of the light, respectively. The plus sign is valid for the companion and the minus sign stands for the star. For further details on the reflection, ellipsoidal and beaming effects see Csizmadia (2020).

Gravity darkening
TLCM is able to model the gravity darkening effect with some simplifications and therefore our treatment is valid only for planets where R planet /R star < 0.2. Our approach was first presented as part of the study in Lendl et al. (2020). Therefore we give only the following details.
Gravity darkening was invoked via a semi-small-planet approximation: the limb darkening and the transit event were calculated via the precise, analytic formulae of Mandel & Agol (2002). The effect of gravity darkening was taken into account in the following way. The local surface effective temperature was calculated from where the local surface potential 2 is with q = M planet /M star : Te polar temperature is not equal to the mean effective temperature of the star in the case of rotating stars (cf. Eq. 8 of Wilson 1979). The mean motion is denoted by n and b the astrographic latitude. The rotational angular velocity ω rot can be calculated from the known stellar radius, the measured or fitted V rot sin I star , and the fitted stellar inclination: The spectroscopically measured (V sin I star ) sp value can be kept fixed during the fit or it can be used as a Gaussian prior. I star is a fitting parameter. The stellar radius R star is taken from isochrones in every iterational steps in the following way: the effective temperature and metalicity of the star is known from spectroscopic measurements while the mean stellar density can be obtained from the scaled semi-major axis a/R star which is strongly related to the transit duration (Seager & Mallén-Ornelas 2003;Winn 2010;Csizmadia et al. 2015): Then, the R star value is given by the corresponding isochrones which are obtained from ρ star , T eff and metalicity as described in Csizmadia (2020).
We fit two angles: the inclination of the stellar rotational vector and its longitude of the node. These, and the planet's skyprojected position, define what is the local temperature behind the planetary disc. Then this temperature is converted via flux convolving the response function of TESS with the spectral library of Munari et al. (2005). (Such conversions are also available for CoRoT, Kepler/K2 and CHEOPS in TLCM now.) The light loss due to transits is given by the Mandel-Agol routines corrected for the gravity darkening by multiplying the limbdarkened intensity behind the planet's apparent center by the normalized gravity darkened fluxes.
Note that longitude of node of the stellar rotational axis -denoted by Ω star -is related to angle between the projected view of the stellar rotational axis and planetary orbit's angular momentum vector via cos λ = ± cos(Ω planet − Ω star ). Since we have set Ω planet = 90 • for sake of simplicity, we have λ = 90 • − Ω star . (The modelling is invariant against this transformation because only the difference between the longitudes of the nodes can be measured from photometry, so one can fix one of them. 3 Barnes et al. (2011) pointed out that photometry does not distinguish between prograde or retrograde rotation of the host stars therefore there is a degeneracy in the modelling results. According to their analysis, the following scenarios are also possible if one gets Ω star as solution: 360 • − Ω star , 180 • − Ω star , 180 • + Ω star as directly follows from the aforemementioned expression of cos λ.

Validation of the gravity darkening apporach
The gravity darkening part of TLCM was tested on the first 18.2 days long data set of Kepler-13Ab. This sanity check used the SAP-FLUX data of Kepler, which were cleaned by a floating median box-car filter. We selected the data points in the ±0.14 × P (P is the orbital period) for this check. We had in total 11,169 data points. We fitted these data with a cosine-like baseline variation (because the reflection effect is well developed in Kepler-13A) and we set the same period, effective temperature It is worth noting that Szabó et al. (2020) fixed the value of λ for their photometric fit at the value obtained by Johnson et al. (2014), we didn't. They used all available Kepler and TESS photometry, we used only 18.2 days of Kepler data for our check. These factors explain our larger error bars. The agreement between us and others are therefore reasonable and validates the TLCM-approach.

Wavelet-based method to remove stellar activity signals and noise-reduction
We used wavelets to remove any stellar activity/variability induced signal and to reduce the noise level stemming from unknown instrumental effects. This method is also able to manage the jumps in the light curve. These jumps or data discontinuities are sudden flux increases due to a cosmic ray impact event or due to instrumental properties after rotating the satellite to repoint the solar panels or stopping observations because of data download, or telescope re-alignment etc. For instance, such data download leaks can be seen in the middle of every light curves in each sector of TESS. The mean data level shift of the same target between sectors of TESS or quarters of Kepler may due to different satellite rotation, different pointing and different contamination level and for us act as flux-jumps again.
The model of TLCM we used is based on Csizmadia (2020). It is a sum of the gravity darkened variations+wavelet based red noise model + radial velocity curve if this latter one also is available. The parameters of the different effects are fitted simultaneously. The wavelet model is based on the work of Carter & Winn (2009). This needs only two parameters to characterize the red (or pink) noise present in the light curve: the white noise level (root mean square, rms) σ w and the red-noise factor σ r -this latter one is not related to the rms of the red noise component.
A difficulty in the application of the wavelets is that we do not know a priori the values of σ w and σ r . When only transits and occultations are present and the out-of-transit light curve part is free of any of beaming, reflection or ellipsoidal variation, then the model gives a normalized flux = 1.0 for all out-of transit and out-of-occultation point. Then the difference between the model and the observations can be used to estimate the wavelet parameters. However, we do not have any points where we know a priori the model flux parameters if out-of-transit variations are present, except the normalization point at phase 0.25. This one point is not enough for the estimation of σ w , σ r .
Therefore we fit the wavelet-parameters simultaneously with the free system parameters and we apply a penalty function (prior) in the fit. The penalty function was based on the requirement that the one-sigma scatter of the fit's residuals must be equal to the average uncertainties of the photometric points. Mathematically, this meant the followings. The residual curve is defined as where O i and M i are the observed and system model fluxes for the ith observations, respectively. r i is the residual of the ith data point. Then we calculate the loglikelihood of the noise model. To do that, we transfer the noise parameters σ w , σ r and all r i values to the routines and algorithm of Carter & Winn (2009). These routines return with the loglikelihood of the noise model (the definition of this likelihood can be found in Carter & Winn 2009). This loglikelihood is penaltized as Here − log L new is the minus loglikelihood to be minimized during the optimization process and used for the error estimation in the MCMC-analysis. − log L is the minus loglikelihood of the wavelet-fit to the residuals given by the algorithm of Carter & Winn (2009). N data is the number of data points. S (r RN ) is the standard deviation of the residuals after removing the red noise component RN i which is also provided by the routines of Carter & Winn (2009): σ i is the photometric uncertainty of the ith observation, while M(σ i ) is the mean of the uncertainties of all photometric individual uncertainties. For the subsequent tests, let the free parameters be the scaled semi-major axis a/R star , the planet-to-star radius ratio R planet /R star , the impact parameter b, the sum and the difference of the linear and quadratic limb darkening coefficients u + = u a + u b and u − = u a − u b , the mass ratio q = M planet /M star , the surface brightness ratio J, the geometric-albedo of the planet A planet , the reflection shift parameter ε, period P, epoch T 0 , and the wavelet-parameters σ w and σ r . We assume a circular orbit and thus we fix e = 0 for the tests. The radius of the star is assumed to be known and it is used as a prior in the way described in Csizmadia (2020).
This approach was tested in the following way. We took 310, 10-day-long segments of 1-minute (short cadence, SC) light curves from the Kepler Q1 database. We convolved these light curves with simulated systems which exhibit all previously mentioned effects: transit, occultation, beaming, ellipsoidal and reflection effects. Then we modelled them with the aforementioned way with TLCM. We plotted the difference between the simulated parameters and the retrieved ones as a function of the signal-to-noise (S/N) ratio. For the ellipsoidal effect (q) we used the following expression of the S/N-ratio: where N is the number of points in the light curve. For other phase-curve parameters (K, ε, A planet ) we have used and for the transit parameters (a/R star , R planet /R star , b, limb darkening coefficients) we have used (34) where we took into account that transit and occultation duration, and hence number of in-transit (in-occultation) points are different at different impact parameters.
In Figures 1-4 we show some examples of the test: the injected light curve, the convolved light curve which contains all the red-noise effect induced in the Kepler-light curves, the comparison of the 'synthetic observations' and the modelling, and finally the red noise corrected light curve and the fits.
We also show the results of the tests in Figures 5-14. From these figures we can read the minimum S/N-ratio needed to get reasonable accuracy in the parameter retrieval. We draw the following conclusions. When stellar variability or instrumental effects are present and they produce red noise in the light curve, we can set the following signal-to-noise ratio limits for the retrieval of the parameters with a wavelet+model fit with the following reasonable accuracies: a/R star : even at S /N ∼ 1 we can get good results (better than 6% relative error) and if S /N > 3 then we can get 2% or better accuracy in the scaled semi-major axis ratio. This is not surprising because we set a Gaussian prior on the stellar radius which is strongly related to this parameter (via transit duration). We can safely assume that in most of the cases we know the stellar radius a priori from SED-fit combined with Gaia-parallax, asteroseismology or from other methods (e.g. Csizmadia 2021, under review at Astronomical Journal) ( Figure 5). b: if S /N > 40 the impact parameter can be retrieved with high accuracy. The impact parameter determination needs a precisely known stellar radius (3% or better). If the stellar radius is less known (3-6%) then most of the solutions lie in a good range, but some outliers appears (gray dots in Figure 6). However, if we translate the impact parameter to inclination via cos i = b/(a/R star ), we find that the inclination values are always better than ±5 • when S /N > 40 (Figure 7). This causes very little difference in the planetary mass when the inclination value is used in the mass-function to determine the planetary mass. R planet /R star : the planet-to-star radius ratio is always better determined than ±2% if the S /N > 50. Even in the range of 10 < S /N < 50 is better than 5%. While Morris et al. (2020) found that this precision cannot be reached with PLATO, our work differs from their work at two points. First, we did not take the effect of granulation into account but Morris et al. (2020) did. We note that if the number of transits are small -in our example it was varied between 6 and 24 -then the granulation is not averaged out but it acts like pseudo-red noise (Chiavassa et al. 2017). Morris et al. (2020) fitted their simulated light curves with Gaussian process and the residuals was still found in the order of 100 ppm. We leave the issue open whether the wavelets can manage the effect of granulation but there is a possibility for that. Second, we considered that the stellar radius is known by at least 2% accuracy as a SED-fit or future asteroseismological part of PLATO will provide this for its primary sample and we used this prior in our fit -they did not. Then, if the granulation is negligible or it can be averaged out by many transit measurements being a white noise, the use of the asteroseismological or SED-fit based stellar radius constraints in the fit are able to provide radius ratio values what PLATO needs even for a Sun-Earth radius ratio (k ∼ 0.009) (Rauer et al. 2014, Figure 8). q: When S /N (q) > 20 then the approach is able to recover the mass ratio with a better accuracy than 10% with some rare exceptions when the reached accuracy is just 20%. This is enough to validate a planet candidate and it may confirm the planetary mass measurement independently of RV (Figure 9). A planet : The geometric albedo of the planet can be retrieved with at least ±0.05 accuracy if S /N (A) > 11( Figure 10).The accuracy increases fast as the signal-to-noise ratio increases. ε: To get the the value of the reflection shift with this wavelet-based filtering method with a ±4 • accuracy one needs at least S /N (A) = 11 while to get it with better accuracy than ±2 • S /N (A) < 25 is needed ( Figure 11). J: To recover the surface brightness ratio of the star and the planet -which is possible form occultations -one needs S /N (J) > 10 ( Figure 12). limb darkening: The limb darkening coefficient combinations u + and u − can be retrieved with ±0.01 accuracy of S /N > 100 but in some rare cases exceptions occur. (Figures 13 and 14).
We add that these results cannot be reached if we do not have a good prior on the stellar radius which helps to constrain the transit duration and thus the impact parameter.
We make the note that, of course, wavelets cannot replace the variable contamination effect. If the contamination is variable from sector to sector of TESS or variable from frame to frame as a consequence of the rotation of CHEOPS, for instance, then this contamination must be corrected before fitting procedure or it must be taken into account in the system's model and not in the wavelet-model. . The y-axis is the difference between the simulated and the retrieved scaled semi-major axis in percentage. Black dots denote the solutions where the stellar radius value was obtained to be with 3% accuracy relative to the injected stellar radius, while gray points represent the cases where we had obtained them with 3-6% accuray.

Summary & Conclusions
We have shown by numerical tests and injecting planetary light curve models that in an ideal case, where the right model of physical reality is well known, the wavelets are able to reconstruct and to filter out the stellar variability and instrumental noise effects, like jumps, cosmic ray hits, discontinuities, detector ageing etc. In Section 4 we have given limits in terms of signal-to-noise ratio for the accuracy of the planet and system parameter retrieval. The wavelet-approach worked well on a wide variety of possible noise sources and stellar variability phenomena, and it was able to manage even high amplitude or fast variable stellar variability and instrumental noise sources (see the examples in Figures 1-4 ).  . The y-axis is the difference between the simulated and the retrieved planet-to-star radius ratio values (k = R planet /R star ). The red and blue points are the small superearths and earths with 0.015 < k < 0.03 (red) and 0.005 < k < 0.015 (blue). Note that these radius ratios correspond to Super/Earthsized (red) and Neptune-sized (blue) planets around a solar-sized star. (The horizontal dashed lines denote ±2% relative errors in the radius ratio. For bigger companions, black dots denote the solutions where the stellar radius value was obtained to be with 3% accuracy relative to the injected stellar radius, while gray points represent the cases where we had obtained them with 3-6% accuray. To reach this high performance, we needed a penalty function during the optimization and uncertainty estimation process. The penalty function decreased the likelihood of the solution if the root mean square of the residuals of the system+wavelet fit deviated from the average uncertainty of the data points. In other words, we prescribed the white noise level what the waveletbased noise modelling must reach. Without such a precondition, there is a danger of overfit, i.e. we fit everything with combi- Fig. 9: The result of the light curve test. The ordinate is the S/Nratio defined by Eq. (31). The abcissa is the difference between the simulated and the retrieved planet-to-star mass ratio values. The horizontal dashed lines denote ±10% relative errors in the radius ratio. nations of wavelets instead of determining the noise properties of the light curve and to extract the system information. This is further illustrated by the example of KELT-9b stellar pulsational like variability in Paper II: it was fully modelled with wavelets and the pulsation was not visible in the red noise corrected fluxresiduals. This example in Paper II sets a caveat: every unmodelled, unkown effect will be incorporated into the wavelets and the information is lost. Therefore, a good model must be selected for the fits when one works with wavelet-based noise-models.
However, data overfit can be done with other methods as well, for example with Gaussian processes. In addition, the wavelet procedure of Carter & Winn (2009) used here, needs only two free parameters. In Gaussian processes, the number of free parameters can be much larger and it can be that one se-Article number, page 9 of 11 A&A proofs: manuscript no. pow  lects a non-appropriate kernel for that noise modelling approach. While the red noise factor of the wavelet-based noise model has no any physical meaning, sometimes the Gaussian process parameters can be linked to some physical process.
We left the limb darkening coefficients free in the test. One can imagine that applying a good prior on the limb darkening may further increase the performance and we can get better results at even lower signal-to-noise ratios. See (Csizmadia et al. 2011) how impact parameter, scaled semi-major axis, planet-tostar radius ratio is degenerated with limb darkening coefficients. However, the present knowledge of limb darkening prefers to leave the limb darkening coefficients free parameters in the fit  We also validated the gravity darkening treatment of TLCM for planets with modelling the Kepler light curve of Kepler-13Ab, a well-known object with asymmetric transits. We found results which are fully compatible with the spin-orbit angle λ obtained by Doppler-tomography results (Johnson et al. 2014) and with the stellar inclination value of Szabó et al. (2020) within 2 degrees which is within the error bars.
The latest version of TLCM with the updated ellipsoidal and reflection effects will be available at http://www.transits. hu once this paper is accepted.