© The Authors 2023

1 Introduction

The light curve of an exoplanetary system may show transits, occultations, and phase-curve variations. The transit technique offers a unique opportunity to accurately determine the radii of transiting exoplanets. By complementing the photometric transit observations with radial velocity data, we can establish the planetary mass and mean density. The phase curves describe the scattering and reflecting properties of an atmosphere at different orbital phase. Available phase curves and occultations are considered to offer the best opportunity to study the three-dimensional (3D) structure of planetary atmospheres (Winn 2010; Parmentier & Crossfield 2018).

In such analyses of the transits as well as the phase curves, the following four problems may arise: (i) The stellar activity, stellar variability (including pulsation and granulation), and instrumental noise cause difficulties in finding and restoring the exact shape of the transit and phase curves (e.g. Oshagh 2018; Sulis et al. 2020; Wong et al. 2020). The transit depth can be also affected by stellar spots yielding the wrong planetary radii (Oshagh et al. 2013; Szabó et al. 2021). If sudden and discontinuous flux variations (jumps) occur in the flux measurements due to a cosmic ray hit or other kind of instrumental effect, then establishing the mean flux level of the host star is difficult. This can lead to further changes in the transit depth because the normalised flux levels are different before and after the jump (e.g. Mislis et al. 2010). (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 significant) must also be separated from the phase variation. The ellipsoidal and reflection effects have higher order harmonics of the orbital frequency. If it is not modelled carefully, a badly modelled ellipsoidal effect can affect the shape of the reflection curve, possibly leading to erroneous conclusions. We refer to Sect. 2 for a detailed explanation of the beaming, reflection, and ellipsoidal effects. (iv) In addition, the exact shape of the phase curve is not known without detailed a priori knowledge of the atmosphere (composition, scattering, and reflecting properties, scale height, clouds, particle sizes of the aerosols, and so on; see Garcia Munoz & Isaak 2015).

We carried out extensive tests on synthetic light curves to overcome the problem identified by point (i), namely, the stellar and instrumental noise sources were modelled by a wavelet transform. We investigated 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 will apply our method to KELT-9b (Csizmadia et al. 2023, hereafter Paper III). We also demonstrate the benefits of the wavelet-based noise handling approach in the precision and accuracy of the transit parameters via a direct comparison to a method where no noise treatment was implemented in Kálmán et al. (2023, hereafter Paper II).

One way to solve problems (ii-iv) is to use prescribed forms of the phase curve and to improve the description of the ellipsoidal effect. In Paper III, we attempted to fit single cosinusoidal and Lambertian, as well as Kane-Gelino- and Kopal-type phase curves (Kopal 1959; Madhusudhan & Burrows 2012; Kane & Gelino 2011; Luger et al. 2022) to the time-series data of KELT-9b obtained by the TESS space telescope. These four different phase functions yield reflection curves that exhibit significantly different shapes. In Sect. 2, we detail the models of the ellipsoidal, beaming, and reflection effects used for the detailed analysis of the KELT-9b light curve in Paper III. We also update the gravity-darkening model of TLCM in Sect. 3. The wavelet model and its validation are presented in Sect. 4. The summary of this study and our conclusions can be found in Sect. 5.

2 Modelling of out-of-transit variations

The out-of-transit 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). Then the out-of-transit variation is defined as the sum of the ellipsoidal, beaming, and reflection effects. Hereafter, we describe our modelling of these effects. There are other kinds of photometric variations that we consider part of our red noise model (see Sect. 4), such as stellar pulsation, stellar activity, possible additional eclipses caused by another star, instrumental, and other (non-white noise-like) effects. We note that we distinguish between the phase curve (only the reflection effect, without the beaming and ellipsoidal effects) and the phase function that contains the time-dependence of the phase curve.

2.1 Phase curve

We utilised the Transit and Light Curve Modeller (TLCM, Csizmadia 2020) code, which expresses the phase-curve variation (F_ph) according to the following form: $$\frac{F_{\text{ph}}}{F_{\text{star}}}=\left(\frac{I_{\text{planet}}}{I_{\text{star}}}{\left(\frac{R_{\text{planet}}}{R_{\text{star}}}\right)}^2+A_{\text{geometric}}{\left(\frac{R_{\text{planet}}}{R_{\text{star}}}\frac{R_{\text{star}}}{d}\right)}^2\text{Φ}\left(\alpha \right)\right).$$(1)

Here, F_ph and F_star are the reflected and stellar fluxes, respectively. The phase-angle α in the phase function Φ is: $$\cos \left(\alpha +\epsilon \right)=\cos \left(\omega +\upsilon \right)\sin i_{\text{p}},$$(2)

where ω is the argument of the periastron of the planetary orbit and v is the true anomaly. The inclination of the planetary orbit is denoted by i_p. In Eq. (1), I is the surface specific intensity in the passband of the observation, d is the mutual star-planet actual distance, R refers to the radii, and A_geometric is the so-called wavelength-dependent geometric albedo.

The angle ε takes into account that there can be a phase shift in the phase curve due to atmospheric circulation, namely, the brightest point of the planet can be shifted eastward or westward relative to the substellar point¹ (Parmentier & Crossfield 2018). The observed values of ε vary widely: 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, that is, by interchanging the indices appropriately and shifting the phase curve by 180°. This may be important for detached eclipsing binary stars, which TLCM can also model, taking into account all these effects for both stellar components (Csizmadia 2020).

2.2 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, cloud properties, and the single scattering albedo (e.g. Garcia Munoz & Isaak 2015). Establishing the exact form of the phase function requires a priori knowledge of the atmospheric properties, which is not always available. In addition, it requires complex and lengthy numerical calculations (Garcia Munoz & Isaak 2015). For the data analysis, analytically expressed approximations are often the best approach. In this paper series, we try the following four phase functions offered by TLCM on KELT-9b, as follows.

The first is a simple cosine function: $$\text{Φ}\left(\alpha \right)=\frac{1}{2}\left(1-\cos \alpha \right).$$(3)

The second is Lambertian, as per: $$\text{Φ}\left(\alpha \right)=\frac{\sin \alpha -\alpha \cdot \cos \alpha }{\pi }.$$(4)

The third is taken from Kane & Gelino (2011), where the phase angle must be measured in degrees: $$\text{Φ}\left(\alpha \right)={10.0}^{-0.4\cdot \left(0.09\left(\alpha /100°\right)\right)+2.39{\left(\alpha /100°\right)}^2)-0.65{\left(\alpha /100°\right)}^3)}.$$(5)

This formula is based on 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 atmospheric particle properties can change 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: $$A\prime \left(d\right)=\left(A_{\text{geometric}}+0.2\frac{{\text{e}}^{d-1}-{\text{e}}^{1-d}}{{\text{e}}^{d-1}+{\text{e}}^{1-d}}\right).$$(6)

Of course, we may question how well this Kane-Gelino phase function performs for hot Jupiters, where the response of the atmosphere can be quite different than in the case of the cooler Jupiter and the terrestrial Venusian atmosphere. We explore this more fully in Paper III, where we apply this phase function to KELT-9b.

The fourth and final phase function tried here is based on the theory of binary star phase function, which accounts for umbral and penumbral effects up to the fourth order in phase-angle. This is taken from Kopal (1959): $$\text{Φ}\left(\alpha \right)=C_0+C_1\cos \alpha +C_2\cos 2\alpha +C_4\cos 4\alpha ,$$(7)

where $$C_0=\frac{8}{3{\pi }^2}{\left(\frac{R_j}{d}\right)}^2+\frac{1}{16}{\left(\frac{R_j}{d}\right)}^3+\frac{2}{3{\pi }^2}{\left(\frac{R_j}{d}\right)}^4-\frac{2K_j}{{\pi }^2}\frac{R_j^2{R}_{3-j}^2}{d^4},$$(8) $$C_1=\frac{1}{3}{\left(\frac{R_j}{d}\right)}^2+\frac{1}{4}{\left(\frac{R_j}{d}\right)}^4,$$(9) $$C_2=\frac{16}{27{\pi }^2}{\left(\frac{R_j}{d}\right)}^2+\frac{3}{16}{\left(\frac{R_j}{d}\right)}^3+\frac{4}{15{\pi }^2}{\left(\frac{R_j}{d}\right)}^4+\frac{4K_j}{3{\pi }^2}\frac{R_j^2{R}_{3-j}^2}{d^4},$$(10) $$C_4=\frac{16}{675{\pi }^2}{\left(\frac{R_j}{d}\right)}^2-\frac{52}{105{\pi }^2}{\left(\frac{R_j}{d}\right)}^4+\frac{4K_j}{15{\pi }^2}\frac{R_j^2{R}_{3-j}^2}{d^4},$$(11)

and the limb darkening correction is given for the reflection in a linear form as: $$K_j=1-\frac{12}{5\pi }\frac{5+\left(\pi -5\right)u_1}{3-u_1}.$$(12)

According to Kopal (1959), C₃ is zero. For the planet, a firstorder approximation consists of taking the linear limb darkening coefficient, u₁ = 0. When we calculate the planet’s reflection effect, j = 2, where R₂ = R_planet; while if one calculates the star’s reflection effect, j = 1, with R₁ = R_star. We note that in the 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 otherwise negligible in star-planet systems but can be important for detached or even closer binary star systems.

The planet-to-star radius ratio (R_planet/R_star) and the scaled semi-major axis ratio (a/R_star) are known from analysis of the transit light curve, or (as in this study) they are fitted simultaneously with the phase curve parameters.

In eccentric orbits, the star-planet distance, d, varies as: $$d=\frac{a\left(1-e^2\right)}{1+\cos \upsilon },$$(13)

where e is the orbital eccentricity (not to be confused with the Euler number in Eq. (6)).

2.3 Dayside and nightside emission

Clearly, the dayside and the nightside phase curves are varying in the opposite phase and therefore we have $${\text{Φ}}_{\text{nightside}}\left(\alpha \right)=1-{\text{Φ}}_{\text{dayside}}\left(\alpha \right).$$(14)

The total planetary phase curve is the sum of the dayside and the nightside emission at a certain phase: $$F_{\text{ph}}=F_{\text{nightside}}\left(1-\text{Φ}\left(\alpha \right)\right)+F_{\text{dayside}}{\left(\frac{a}{d}\right)}^2\Phi \left(\alpha \right).$$(15)

Comparing Eq. (15) to Eq. (1), we can relate the measured quantities to the parameters we are searching for via: $$\frac{F_{\text{nightside}}}{F_{\text{star}}}=\frac{I_{\text{planet}}}{I_{\text{star}}}{\left(\frac{R_{\text{planet}}}{R_{\text{star}}}\right)}^2$$(16)

and $$\frac{F_{\text{dayside}}}{F_{\text{star}}}=\frac{F_{\text{nightside}}}{F_{\text{star}}}+A_{\text{geometric}}{\left(\frac{R_{\text{planet}}}{R_{\text{star}}}\frac{R_{\text{star}}}{a}\right)}^2.$$(17)

We note that I_planet/I_star, R_planet/R_star, A_geometric and the reciprocal of R_star/a are fitting parameters in TLCM and they can be measured from the simultaneous fit to the transit, occultation, and phase curves.

We note that we assumed the very same phase function for the reflected light and for the dayside and nightside emission in TLCM, which is, of course, just an approximation of reality. However, we used the simplest possible approach.

2.4 Ellipsoidal effect

The flux variation caused by the ellipsoidal shapes of the components is characterised following Kopal (1959; j = 1 for the star and j = 2 for the planet): $$\begin{array}{l}\frac{F_{\text{ellipsoidal,j}}}{F_{\text{star}}}\hfill \\ =f_j\frac{M_{3-j}}{M_j}{\displaystyle {\sum }_{l=2}^4{w}_{jl}}{\left(\frac{1+e\cos \upsilon }{1-e^2}\right)}^l{\left(\frac{R_j}{a}\right)}^{l+2}{p}_l\left(cos\left(\upsilon +\omega \right)\sin i_{\text{p}}\right),\hfill \end{array}$$(18)

with the gravity-darkening correction, $$\tau =\frac{4\times 14388.\text{0μm}\cdot K}{\lambda T_{\text{eff,}j}\left(1.0-e^{-14388.0/\lambda /T_{\text{eff,}j}}\right)}.$$(19)

Here, hc/k_B = 14388µm · K, the effective wavelength of the observation λ is given in microns, T_eff is in kelvins, and according to Kopal (1959): $$w_2=\frac{2\left(15+u_1\right)\times \left(1+k_2\right)}{5\left(6-2u_1-3u_2\right)\left(1+\frac{\tau }{4}\left(\frac{5}{1+k_2}-1\right)\right)},$$(20) $$w_3=\frac{\left(35u_1+48u_2\right)\times \left(1.0+k_3\right)}{7\left(6-2u_1-3u_2\right)\left(1.0+\frac{\tau }{10}\left(\frac{7}{1+k_3}-2.0\right)\right)},$$(21) $$w_4=\frac{9}{8}\frac{9\left(4u_1+7u_2-4\right)\times \left(1.0+k_4\right)}{8\left(6-2u_1-3u_2\right)\times \left(1.0+\frac{\tau }{18}\left(\frac{9.0}{1.0+k_4}-3.0\right)\right)}.$$(22)

The 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₁ = u_a + 2 · u_b and u₂ = −u_b. We note, of course, that 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, metallicity, radius, and evolutionary status. The k_i; values of the primary star or of the host star are fixed at their theoretically calculated values from Claret (2004). We note that the apsidal-motion constant is half of the Love number (Csizmadia et al. 2019). In addition, in Eq. (18), we have f₁ = 1 for the star and f2 = I_planet/I_star(R_planet/R_star)² for the planet (or secondary star) in the respective photometric passband (same as in Eq. (16)).

2.5 Beaming effect

The beaming effect is characterised as $$B_j=\pm f_j{\left(\frac{R_j}{a}\right)}^2{\text{Ω}}_j\left(T_{\text{eff,}j},\log g_j,Z_j\right)\frac{K_j}{c}\left(e\cos \omega +cos\left(\upsilon \left(t\right)+\omega \right)\right).$$(23)

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

3 Gravity darkening

3.1 Gravity-darkening model

The previous version of TLCM was able to model the gravity-darkening effect with some simplifications and, therefore, our former treatment was valid only for planets where R_planet/R_star < 0.2. That approach was first presented as part of Lendl et al. (2020). In a nutshell, that approach projected the intensity distribution of a gravity-darkened star to the spherical stellar surface. The transit event was calculated via the usual spherical star-spherical planet assumption (Mandel & Agol 2002) and the light loss was corrected for by modifying the light loss event with the intensity of the gravity-darkened star behind the planetary centre (Lendl et al. 2020).

In the newer version of TLCM, a simplification-free model is used for modelling the gravity-darkening phenomenon. We used a full Roche-geometry (Wilson 1979) to describe the exact shape of the star and to calculate the viewing angles and derivatives of the potential on the surface. The new approach is valid for any planet-to-star radius ratio.

The effect of gravity darkening is taken into account in the following way. The local surface effective temperature is calculated from: $$T_{\text{local}}=T_{*}{\left(\frac{\left|\nabla V\right|}{{\left|\nabla V\right|}_{\text{pole}}}\right)}^{\beta },$$(24)

where the local surface potential² is with q = M_planet/M_star: $$V=\frac{n^2{a}^3}{\left(1+q\right)r}+\frac{1}{2}{\omega }_{\text{rot}}^2{r}^2{\sin }^{2}b.$$(25)

The 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 is the stellar 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: $${\omega }_{\text{rot}}=\frac{{\left(V\sin I_{\text{star}}\right)}_{\text{sp}}}{R_{\text{star}}\sin I_{\text{star}}}.$$(26)

While I_star is a fitting parameter, 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. The stellar radius, R_star, is taken from isochrones at every iterational step in the following way: the effective temperature and metallicity of the star are 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): $${\rho }_{\text{star}}=\frac{3\pi }{G{P}^2\left(1+q\right)}{\left(\frac{a}{R_{\text{star}}}\right)}^3.$$(27)

Then, the R_star value is given by the corresponding isochrones which are obtained from ρ_star, T_eff, and metallicity, as described in Csizmadia (2020).

In this new approach, we calculate the total flux coming from the unobscured visible stellar surface that is used to normalise the light curve; then, we calculate numerically the light loss behind the visible disc of the planet with the help of a 2D Gauss-Legendre integration.

We fit two angles: the inclination of the stellar rotational vector and its longitude of the node. These values, along with the planet’s sky-projected position, define 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). We note that such conversions are also available for CoRoT, Kepler/K2, and CHEOPS in TLCM at present.

We note that longitude of node of the stellar rotational axis (denoted by Ω_star) is related to the angle between the projected view of the stellar rotational axis and the 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 we can fix one of them.³ Barnes et al. (2011) pointed out that photometry does not distinguish between prograde or retrograde rotation of the host star; therefore, there is a degeneracy in the modelling results. According to their analysis, the following scenarios are also possible if we obtain Ω_star as a solution: 360° − Ω_star, 180° − Ω_star, and 180° + Ω_star – as directly follows from the aforementioned expression of cos λ.

3.2 Validation of the gravity-darkening approach

The gravity-darkening part of TLCM was tested on the first 18.2 days of the Kepler light curve 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 for this check, which have a phase of ±0.14 around the primary eclipse. We obtained a total of 11,169 data points. We fit these data with a cosine-like baseline variation and we applied the same period, effective temperature (T_eff = 8600 K), and contamination value as in Szabó et al. (2020). We set V sin I_star = 76.96 km s⁻¹ (Johnson et al. 2014). All other parameters were free. We compare our results to those of Szabó et al. (2020) and Johnson et al. (2014):

i_star: Szabó et al. (2020) found, from Kepler and TESS extensive photometry, the stellar rotational axis inclination to be i_star = 102.5° ± 0.8°, while we have found i_star = 103.2° ± 2.7° with the old, simplified approach (Lendl et al. 2020) and i_star = 110.0° ± 2.1° with the present new, Roche-geometry approach.

λ: Johnson et al. (2014) found, from spectroscopic Rossiter-Mclaughlin measurements, the projected stellar obliquity to be λ = 58.6° ± 2.0°, while we have found λ = 55.9° ± 13.8° with the old, simplified approach (Lendl et al. 2020) and λ = 57.1° ± 9.1° with the present new, Roche-geometry approach.

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), however, we did not. Those authors used all available Kepler and TESS photometry, whereas 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; the TLCM approach is validated.

4 Wavelet-based method to remove stellar activity signals and noise-reduction

We used wavelets to remove any signal induced by stellar activity or variability and to reduce the noise level stemming from unknown instrumental effects. This method is also able to correct for 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 re-point the solar panels or stopping observations to download data or telescope re-alignment, and so on. For instance, such data download gaps can be seen in the middle of every light curve in each sector of TESS. The mean data level shift of the same target between sectors of TESS or quarters of Kepler may be due to different satellite rotation, pointing, and contamination levels; so, for our purposes, they can act as flux jumps once again.

The model of TLCM we used is based on Csizmadia (2020). It is a sum of the gravity-darkened transit + occultation + beaming + reflection + ellipsoidal variations + wavelet-based red noise model + radial velocity curve (if available). The parameters describing the different effects are fitted simultaneously.

The wavelet model is based on the work of Carter & Winn (2009). This needs only two parameters to characterise 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. We note that the latter is not related to the rms of the red noise component.

One 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 beaming, reflection, or ellipsoidal variation, then the model gives a normalised flux = 1.0 for all out-of-transit and out-of-occultation points. 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 normalisation point at phase 0.25. This one point is not enough to estimate σ_w and σ_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 residuals to the fit must be equal to the average uncertainties of the photometric points. Mathematically, this meant the following: the residual curve is defined as: $$r_i=O_i-M_i,$$(28)

where O_i and M_i are, respectively, the observed and model fluxes for the ith observations. The residual of the ith data point is denoted by r_i. Then we calculated the log-likelihood of the noise model. To do so, we transferred the noise parameters, σ_w, σ_r, and all r_i values to the routines and algorithm of Carter & Winn (2009). These routines return with the log-likelihood of the noise model (the definition of this likelihood can be found in Carter & Winn 2009). This log-likelihood is penalised as: $$-\log L_{\text{new}}=-\log L+0.5\times N_{\text{data}}\times {\left(\frac{S\left(r_{\text{RN}}\right)}{M\left({\sigma }_i\right)}-1.0\right)}^2,$$(29)

where σ_i is the photometric uncertainty of the ith observation and M(σ_i) is the mean of the individual photometric uncertainties. Here, log L_new is the log-likelihood to be minimised during the optimisation process and used for the error estimation in the MCMC analysis. Also, log L is the log-likelihood of the wavelet fit to the residuals given by the algorithm of Carter & Winn (2009); N_data is the number of data points; and 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): $$r_{\text{RN,i}}=O_i-M_i-R{N}_i.$$(30)

For the subsequent tests, the free parameters are 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₂, the reflection shift parameter, ε, period, P, epoch, T₀, 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 last approach was tested in the following way. We took 310 ten-day long segments of 1 min (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. The parameter range tested can be found in Table 1 and the effects occurred in the light curves and successfully tested here are listed in Table 2. For all tests in this paper, we selected random combinations of the input parameters from the ranges below (or we kept their values fixed as noted in Table 1). We note that more tests were done in the R_planet/R_star < 0.2 range than over this limit. The transit impact parameter was always checked to provide b < 1 + R_planet/R_star and if this was not fulfilled (condition of having a transit for circular orbits), a new, random b was selected until this condition holds. We assumed circular orbits for all tests but this does not limit the results only to non-eccentric orbits. We note the high range of K (RV amplitude), which is due to some of the large mass ratios we included in the tests. Of our test cases, 28% had mass ratios less than 0.005 and 50% of them had less than 0.03. We then modelled these light curves in the above-mentioned 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 for the S/N ratio: $$S/N_{\left(q\right)}=\frac{q{\left(\frac{R_{\text{star}}}{a}\right)}^3}{\sqrt{{\sigma }_{\text{w}}^2+{\sigma }_{\text{r}}^2}}\times \sqrt{N},$$(31)

where N is the number of points in the light curve. For other phase-curve parameters (K, ε, A₂), we used: $$S/N_{\left(A\right)}=\frac{A_2{\left(\frac{R_{\text{planet}}}{R_{\text{star}}}\right)}^2{\left(\frac{R_{\text{star}}}{a}\right)}^2}{\sqrt{{\sigma }_{\text{w}}^2+{\sigma }_{\text{r}}^2}}\times \sqrt{N},$$(32)

while for J we used: $${S/N}_{\left(J\right)}=\frac{J{\left(\frac{R_{\text{planet}}}{R_{\text{star}}}\right)}^2}{\sqrt{{\sigma }_{\text{w}}^2+{\sigma }_{\text{r}}^2}}\sqrt{\frac{P{R}_{\text{star}}}{\pi a{t}_{exp}}\sqrt{{\left(1+\frac{R_{\text{planet}}}{R_{\text{star}}}\right)}^2-b^2}}\sqrt{N_{\text{occultations}}},$$(33)

and for the transit parameters (a/R_star, R_planet/R_star, b, limb darkening coefficients), we used: $$S/N=\frac{{\left(\frac{R_{\text{planet}}}{R_{\text{star}}}\right)}^2\sqrt{\frac{P{R}_{\text{star}}}{\pi a{t}_{exp}}\sqrt{{\left(1+R_{\text{planet}}/R_{\text{star}}\right)}^2-b^2}}}{\sqrt{{\sigma }_{\text{w}}^2+{\sigma }_{\text{r}}^2}}\sqrt{N_{\text{transit}}},$$(34)

where we took into account that transit and occultation duration and, hence, the number of in-transit (in-occultation) points are different for different impact parameters.

In Figs. 1–4, we show some examples of the tests: the injected light curve, the convolved light curve which contains all the red noise effect included in the Kepler light curves, the comparison of the synthetic ‘observations’, and the modelling, along with the red-noise-corrected light curve and the fits. We give further examples of the performance of the simultaneous fit of the transit + phase curve model + wavelet-based noise model in Appendix A.

We also show the results of the tests in Figs. 5–14. From these figures we can read the minimum S/N ratio needed to get reasonable accuracy in the parameter retrieval. The accuracy hereafter is defined by dividing the difference between the injected and retrieved values of the parameter by their respective injected values. For instance, PLATO defines the desired accuracy in planet-to-star radius in this way (Rauer et al. 2014). In some cases, the accuracy is defined differently from this, as the absolute difference between injected and retrieved parameter (b, i, ε, A₂, J, and the limb-darkening coefficients). Such a distinction is necessary in the definition of accuracy because of the presence of cyclic variables (angles).

Our conclusions are as follows. First, 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 that is strongly related to this parameter (via transit duration). We can safely assume that in most of the cases, we can obtain the stellar radius a priori from the SED fit combined with the Gaia parallax, asteroseismology, or other methods (Fig. 5).

• b: if S/N > 40, the impact parameter can be retrieved with high accuracy. To determine the impact parameter, we need a precisely known stellar radius (3% or better). If the stellar radius is less precisely known (3–6%), then most of the solutions lie in a good range, but some outliers appear (gray dots in Fig. 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 (Fig. 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, R_planet/R_star is better than 5%. While Morris et al. (2020) found that this precision cannot be reached with PLATO, our work differs from theirs in two ways. Firstly, 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 made to vary between 6 and 24 – then the granulation is not averaged out but it acts similarly to pseudo-red noise (Chiavassa et al. 2017). Morris et al. (2020) fitted their simulated light curves using a Gaussian process model and the residuals were still found to be on the order of 100 ppm. We do not resolve here the question of whether the wavelet method can model the effects of granulation, but it remains a possibility. Secondly, we considered that the stellar radius is known to a precision of (at worst) 2% (as for PLATO’s primary sample, this will be provided SED fitting or by asteroseismology) and we used this prior in our fit, whereas Morris et al. (2020) did not. If the granulation is negligible or being a white noise it can be averaged out by many transit measurements. However, when we observe only two or three transits, granulation can play a role in radius ratio determination (Chiavassa et al. 2017). Then the use of the stellar radius constraint from asteroseismology or an SED fit results is needed, especially for small planets – as in the case of an Earth-Sun radius ratio (k ~ 0.009; Rauer et al. 2014, Fig. 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 accuracy achieved is just 20%. This is enough to validate a planet candidate and may even confirm the planetary mass measurement independently of RV measurements (Fig. 9).

• A₂: The geometric albedo of the planet can be retrieved with at least ±0.05 accuracy if S/N_(A) > 11 (Fig. 10). The accuracy increases steeply with increasing S/N ratio.

• ε: To get the value of the reflection shift with this wavelet-based filtering method with a ±4° accuracy, we need at least S/N_(A) = 11, while to get it with better accuracy than ±2°, S/N_(A) > 25 is required (Fig. 11).

• J: To recover the surface brightness ratio of the star and the planet – which is possible from occultations – we need S/N_(J) > 10 (Fig. 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 (Figs. 13 and 14).

These results cannot be achieved 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 note that 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 the fitting procedure. Otherwise, it must be modelled explicitly and not subsumed into the wavelet model.

Fig. 1 Example of performance of the wavelet-based light-curve fit when jumps are present in the light curve. a: Kepler Q1 light curve segment convolved with the injected model, used for the test (Kepler target: 001571088). b: injected model light curve. c: raw Kepler Q1 light curve segment (black dots) and the model+wavelet fit. d: the red noise-corrected light curve (raw flux - wavelet based red noise, black dots) and the model fit (red line).

Fig. 2 Example of performance of the wavelet-fit when stellar spot modulation is present in the light curve with the same layout as Fig. 1. Kepler target: 002556755 with injected test light curve.

Fig. 3 Example of performance of the wavelet-fit when stellar pulsation-like stellar variability is present in the light curve with the same layout as Fig. 1. Kepler target: 004044353 with injected test light curve.

Fig. 4 Example of performance of the wavelet-fit when stellar pulsation-like stellar variability is present in the light curve with the same layout as Fig. 1. Kepler target: 010285114 with injected test light curve.

5 Summary and conclusions

Using numerical tests and synthetic planetary light curve models, we have shown that in an ideal case, where the right model of physical reality is well known, wavelets are able to reconstruct and to filter out stellar variability and instrumental noise effects, such as jumps, cosmic ray hits, discontinuities, detector ageing, and so on. In Sect. 4, we have provided limits in terms of signal-to-noise ratio for the accuracy of the planet and system parameter retrieval. The wavelet approach has worked well on a wide variety of possible noise sources and stellar variability phenomena and it was able to manage even high-amplitude or rapid stellar variability and instrumental noise sources (see the examples in Figs. 1–4).

To reach this high level of performance, we needed a penalty function during the optimisation 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 that must be reached by the wavelet-based noise modelling. Without such a precondition, there is a danger of overfitting, namely, we fit everything with combinations of wavelets instead of extracting the system information. This is further illustrated by the example of stellar pulsational such as the variability of KELT-9 in Paper III: it was fully modelled with wavelets and the pulsation was not visible in the red noise-corrected flux-residuals. This example in Paper III sets a caveat: every unmodelled, unknown effect will be incorporated into the wavelets and the information is lost. Therefore, a good model must be selected for the fits when we are working with wavelet-based noise models.

However, data overfitting can occur 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. For Gaussian processes, the number of free parameters can be much larger and it can be that one selects an inappropriate kernel for that noise modelling approach. While the red noise factor of the wavelet-based noise model has no physical meaning, sometimes the Gaussian process parameters can be linked to some physical process.

We left the limb-darkening coefficients free in the test. We can assume that applying a good prior on the limb darkening may further increase the performance and we can get better results at even lower S/N ratios. We refer to Csizmadia et al. (2011) for details of how the impact parameter, scaled semi-major axis, and planet-to-star radius ratio are degenerate with limb-darkening coefficients. However, the present knowledge of limb darkening results in a preference for leaving the limb darkening coefficients as free parameters in the fit (Csizmadia et al. 2013; Espinoza & Jordan 2015; Agol et al. 2020).

We also validated the gravity-darkening treatment of TLCM for planets by modelling the Kepler light curve of Kepler-13Ab, a well-known object with asymmetric transits. We found results that 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 (i.e. within the error bars).

The latest version of TLCM with the updated ellipsoidal and reflection effects have been made available online⁴ upon publication.

Fig. 5 Result of the light curve test for the a/Rstar parameter. The ordinate is the S/N-ratio defined by Eq. (34). 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% accuracy.

Fig. 6 Result of the light curve test for the impact parameter b. The ordinate is the S/N-ratio defined by Eq. (34). The y-axis is the difference between the simulated and the retrieved impact parameter. We also plotted the 1σ error bar of the impact parameter for this figure. See the meaning behind the black and gray points in Fig. 5.

Fig. 7 Result of the light curve test for the inclination. The ordinate is the S/N-ratio defined by Eq. (34). The y-axis is the difference between the simulated and the retrieved inclination values. The vertical lines are the 1σ error bars. See the meaning of the black and gray points in Fig. 5.

Fig. 8 Result of the light curve test for the planet-to-star radius ratio parameter, k = Rplanet/Rstar. The ordinate is the S/N-ratio defined by Eq. (34). The y-axis is the difference between the simulated and the retrieved planet-to-star radius ratio values (k = Rplanet/Rstar). The red and blue points are the small super-Earths and Earths with 0.015 < k < 0.03 (red) and 0.005 < k < 0.015 (blue). We note that these radius ratios correspond to Super-Earth- and Earth-sized (red) as well as Neptune-sized (blue) planets around a solar-sized star. We also note that 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 with a 3% accuracy relative to the injected stellar radius, while gray points represent the cases where we had obtained them with a 3–6% accuray.

Fig. 9 Result of the light curve test for the mass ratio parameter q = Mplanet/Mstar. The ordinate is the S/N-ratio 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.

Fig. 10 Result of the light curve test for the albedo. The ordinate is the S/N-ratio defined by Eq. (32). The y-axis is the difference between the simulated and the retrieved planetary albedo values. The horizontal dashed lines denote ±0.05 absolute errors in albedo-determination. Note: in this figure, the x-axis has a logarithmic scale for better visibility. See the meaning of the black and gray points in Fig. 5.

Fig. 11 Result of the light curve test for the reflection shift parameter ε (cf. Eq. (2)). The ordinate is the S/N-ratio defined by Eq. (32). The abscissa is the difference between the simulated and the retrieved reflection shift values. The horizontal dashed lines denote ±2° relative errors in the radius ratio. Note: in this figure, the x-axis has a logarithmic scale for better visibility. See the meaning behind the black and gray points in Fig. 5.

Fig. 12 Result of the light curve test for the surface brightness ratio. The ordinate is the S/N-ratio defined by Eq. (33). The abcissa is the difference between the simulated and the retrieved planet-to-star surface brightness ratio values. Note: in this figure, the x-axis has a logarithmic scale for better visibility. See the meaning of the black and gray points in Fig. 5.

Fig. 13 Result of the light curve test for the limb darkening coefficient u+. The ordinate is the S/N-ratio defined by Eq. (34). The abscissa is the difference between the simulated and the retrieved u+ limb darkening coefficient combination values. See the meaning of the black and gray points in Fig. 5.

Fig. 14 Result of the light curve test for the limb darkening coefficient u−. The ordinate is the S/N-ratio defined by Eq. (34). The abscissa is the difference between the simulated and the retrieved u− limb darkening coefficient combination values. See the meaning of the black and gray points in Fig. 5.

Acknowledgements

The authors gratefully acknowledge the European Space Agency and the PLATO Mission Consortium, whose outstanding efforts have made these results possible. We thank DFG Research Unit 2440: “Matter Under Planetary Interior Conditions: High Pressure, Planetary, and Plasma Physics” for support. We also acknowledge support by DFG grants RA 714I14-1 within the DFG Schwerpunkt SPP 1992: “Exploring the Diversity of Extrasolar Planets”. K.Sz. acknowledges the support of the PRODEX Experiment Agreement No. 4000137122 between the ELTE Eötvös Loránd University and the European Space Agency (ESA-DISCI-LE-2021-0025) and the Lendület LP2018-7 2022 grant of the Hungarian Academy of Science. Project no. C1746651 has been implemented with the support provided by the Ministry of Culture and Innovation of Hungary from the National Research, Development and Innovation Fund, financed under the NVKDP-2021 funding scheme.

Appendix A Modelling test result examples

In Figures A.1–10, we give further examples of how the simultaneously fitted transit+phase curve model + wavelet-based model performs on the Kepler Q1 light curves with injected transits.

Fig. A.1 Kepler target 004851239, 1 min cadence data. Here (and in other Figures of the appendix) black dots are the SAP fluxes while the red line is the simultaneous fit of the transit+phase curve model and the noise model. Note: the transits and phase curve variations here and in other test cases in this paper were not in the original Kepler data: they were injected by us for test purposes.

Fig. A.2 Kepler target 004566474.

Fig. A.3 Kepler target 004464952.

Fig. A.4 Kepler target 002284679.

Fig. A.5 Kepler target 001719472.

Fig. A.6 Kepler target 012062443.

Fig. A.7 Kepler target 012020590.

Fig. A.8 Kepler target 011873252.

Fig. A.9 Kepler target 011128126.

Fig. A.10 Kepler target 011127190.

References

All Tables

All Figures

Fig. 1 Example of performance of the wavelet-based light-curve fit when jumps are present in the light curve. a: Kepler Q1 light curve segment convolved with the injected model, used for the test (Kepler target: 001571088). b: injected model light curve. c: raw Kepler Q1 light curve segment (black dots) and the model+wavelet fit. d: the red noise-corrected light curve (raw flux - wavelet based red noise, black dots) and the model fit (red line).

In the text

	Fig. 2 Example of performance of the wavelet-fit when stellar spot modulation is present in the light curve with the same layout as Fig. 1. Kepler target: 002556755 with injected test light curve.
In the text

	Fig. 3 Example of performance of the wavelet-fit when stellar pulsation-like stellar variability is present in the light curve with the same layout as Fig. 1. Kepler target: 004044353 with injected test light curve.
In the text

	Fig. 4 Example of performance of the wavelet-fit when stellar pulsation-like stellar variability is present in the light curve with the same layout as Fig. 1. Kepler target: 010285114 with injected test light curve.
In the text

Fig. 5 Result of the light curve test for the a/Rstar parameter. The ordinate is the S/N-ratio defined by Eq. (34). 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% accuracy.

In the text

	Fig. 6 Result of the light curve test for the impact parameter b. The ordinate is the S/N-ratio defined by Eq. (34). The y-axis is the difference between the simulated and the retrieved impact parameter. We also plotted the 1σ error bar of the impact parameter for this figure. See the meaning behind the black and gray points in Fig. 5.
In the text

	Fig. 7 Result of the light curve test for the inclination. The ordinate is the S/N-ratio defined by Eq. (34). The y-axis is the difference between the simulated and the retrieved inclination values. The vertical lines are the 1σ error bars. See the meaning of the black and gray points in Fig. 5.
In the text

Fig. 8 Result of the light curve test for the planet-to-star radius ratio parameter, k = Rplanet/Rstar. The ordinate is the S/N-ratio defined by Eq. (34). The y-axis is the difference between the simulated and the retrieved planet-to-star radius ratio values (k = Rplanet/Rstar). The red and blue points are the small super-Earths and Earths with 0.015 < k < 0.03 (red) and 0.005 < k < 0.015 (blue). We note that these radius ratios correspond to Super-Earth- and Earth-sized (red) as well as Neptune-sized (blue) planets around a solar-sized star. We also note that 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 with a 3% accuracy relative to the injected stellar radius, while gray points represent the cases where we had obtained them with a 3–6% accuray.

In the text

	Fig. 9 Result of the light curve test for the mass ratio parameter q = Mplanet/Mstar. The ordinate is the S/N-ratio 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.
In the text

Fig. 10 Result of the light curve test for the albedo. The ordinate is the S/N-ratio defined by Eq. (32). The y-axis is the difference between the simulated and the retrieved planetary albedo values. The horizontal dashed lines denote ±0.05 absolute errors in albedo-determination. Note: in this figure, the x-axis has a logarithmic scale for better visibility. See the meaning of the black and gray points in Fig. 5.

In the text

Fig. 11 Result of the light curve test for the reflection shift parameter ε (cf. Eq. (2)). The ordinate is the S/N-ratio defined by Eq. (32). The abscissa is the difference between the simulated and the retrieved reflection shift values. The horizontal dashed lines denote ±2° relative errors in the radius ratio. Note: in this figure, the x-axis has a logarithmic scale for better visibility. See the meaning behind the black and gray points in Fig. 5.

In the text

	Fig. 12 Result of the light curve test for the surface brightness ratio. The ordinate is the S/N-ratio defined by Eq. (33). The abcissa is the difference between the simulated and the retrieved planet-to-star surface brightness ratio values. Note: in this figure, the x-axis has a logarithmic scale for better visibility. See the meaning of the black and gray points in Fig. 5.
In the text

	Fig. 13 Result of the light curve test for the limb darkening coefficient u+. The ordinate is the S/N-ratio defined by Eq. (34). The abscissa is the difference between the simulated and the retrieved u+ limb darkening coefficient combination values. See the meaning of the black and gray points in Fig. 5.
In the text

	Fig. 14 Result of the light curve test for the limb darkening coefficient u−. The ordinate is the S/N-ratio defined by Eq. (34). The abscissa is the difference between the simulated and the retrieved u− limb darkening coefficient combination values. See the meaning of the black and gray points in Fig. 5.
In the text

	Fig. A.1 Kepler target 004851239, 1 min cadence data. Here (and in other Figures of the appendix) black dots are the SAP fluxes while the red line is the simultaneous fit of the transit+phase curve model and the noise model. Note: the transits and phase curve variations here and in other test cases in this paper were not in the original Kepler data: they were injected by us for test purposes.
In the text

	Fig. A.2 Kepler target 004566474.
In the text

	Fig. A.3 Kepler target 004464952.
In the text

	Fig. A.4 Kepler target 002284679.
In the text

	Fig. A.5 Kepler target 001719472.
In the text

	Fig. A.6 Kepler target 012062443.
In the text

	Fig. A.7 Kepler target 012020590.
In the text

	Fig. A.8 Kepler target 011873252.
In the text

	Fig. A.9 Kepler target 011128126.
In the text

	Fig. A.10 Kepler target 011127190.
In the text