Stellar pulsations interfering with the transit light curve: configurations with false positive misalignment

We studied the effects of a transit on the pulsational pattern in a generic model that was based on the pixel-level in-silico visualisation of the pulsating stellar surface while it was occulted by the planet. The local intensities were summed up pixel-by-pixel, leading to a synthetic brightness of the star, when it exhibited a certain pulsational pattern that was partially occulted by the planet at a certain position along the transit chord. This technique enabled an extensive analysis of the different pulsation modes, and we could identify which of these are the most sensitive to planetary transits, i.e. which are the modes where planetary transits most commonly cause virtual amplitude and phase modulations due to a transiting companion.

We wanted to avoid fallacies caused by so-called inverse crimes, when the same model used both for simulating the data and then again for fitting them results in a much better recovery of the input parameters than in real applications, leading to overoptimistic conclusions. Therefore we used one set of tools for creating the synthetic light curves (fitsh/lfit; pulzem, see Pál, 2012; Bíró & Nuspl, 2011), and another tool for the modelling (TLCM, see Csizmadia, 2020).

The light variation due to the transits of the static surface brightness component was modelled using the fitsh package, which uses analytical formulae (Mandel & Agol, 2002) to compute a resolution-independent light curve. Its lack of capability to model rotationally distorted, gravity-darkened stars did not pose any disadvantage because our input configuration was just a slowly rotating, undistorted star with no gravity darkening. One of the input parameters of fitsh is a combination of the period $P$, relative semi-major axis $a/R_{\text{s}}$ and conjunction parameter $b$:

The photometric dataset was modelled with a quadratic limb-darkening law with coefficients $\mu_{1}=0.209$ and $\mu_{2}=0.216$ chosen during preliminary tests and kept fixed throughout the investigation.

The light variation due to the pulsations was generated with the forward modelling utilities of the pulzem package (Bíró & Nuspl, 2011). Created to perform eclipse mapping of pulsation patterns in eclipsing binaries, it properly takes into account the partial occultation of the surface pulsation patterns, and handles planet-sized eclipsing objects too. The geodesic grid scheme based on Hendry & Mochnacki (1992) was implemented, involving a sub-triangulation of the faces of an icosahedron inscribed into the stellar volume and then radially projecting the resulting grid onto the surface. The procedure yields a nearly uniform surface resolution with the least number of triangular tiles. In addition, the intersection of the eclipser’s shadow - a circular disc in our case - with the triangular tiles is exactly taken into account. The usual features of limb darkening, rotational distortion and gravity darkening of the host star are also properly implemented – although in our case the latter two were not required in the forward modelling. We gradually increased the grid resolution until the change in the normalised flux of the transit curve decreased below the level of the uncorrelated error component. The condition was fulfilled at the value of $15$ for the sub-triangulation parameter, i.e. each side of the icosahedron is divided into 15 segments, giving $15^{2}=225$ tiles for each of the $20$ faces, and resulting in $20\cdot 15^{2}=4500$ tiles for the whole stellar surface, half of which is actually visible for a stellar disc. Rotation of the static pattern is not required, while it only requires a slightly modified frequency for the non-linear pulsation patterns, so the stellar surface actually needs no rotation at all. The same limb darkening used for the static light curve was applied to the pulsation surface patterns too, before generating the combined pulsation flux of all modes.

The reason behind using one tool for synthesizing the static flux and another one for the pulsating flux component is that there is no suitable tool available yet capable of computing both components with the same accuracy. fitsh is virtually exact for transit light curves, but cannot submit arbitrary time-dependent surface patterns to the same transit computation. pulzem computes both, but its accuracy is intrinsically limited by the rasterization approach. Static transit light curves of fitsh and pulzem showed differences that could not be diminished by increasing the surface resolution beyond all limits. Although the pulsation patterns suffer from the same relative inaccuracies, their small amplitudes explored in this study, of a few percents of the transit depth, render these errors below the applied artificial uncorrelated noise component. Therefore the use of this combined simulation method is appropriate.

The first investigated scenario was a test system consisting of a spherical star and a planet with a relative size $R_{\text{p}}\sim$ 0.1 $R_{\text{s}}$, orbiting its host star with a period of about 2 days. This is a typical hot Jupiter configuration. The orbit was completely circular and seen edge-on (inclination is 90 degrees). Uncorrelated random noise was added to the light curve. Its level was chosen to match the error of an average target of  10-11 mag observed with the TESS mission, $\sim 0.0002$ in the normalised flux. The modelled transits had a time duration of 0.1 orbital phases, and the depth of the transit in was approximately 0.013. Table 1 summarises the input parameters of the two tools used for the simulation. Each parameter set was used for different investigation as written in the following subsections.

For the second scenario, we used an identical system with the only difference that an inclination value of $89.89^{\circ}$ was used, corresponding to the WASP-33 system which ultimately inspired this work.

The surface intensity pattern of the various small amplitude non-radial pulsations in stars with no significant rotational and tidal distortion can be described by spherical harmonics $Y_{\ell}^{m}(\theta,\varphi)$, with degree $\ell$ and azimuthal order $m$ (Aerts et al., 2010). The integrated flux coming from an uneclipsed, fully visible stellar disc naturally does not contain any information on the nonradial modes of the pulsations. During eclipses (or transits, in exoplanetary parlance) an area moving through the stellar disc is excluded from the integration, leading to a modulation of the periodic signals. The instant amplitudes and initial phases will show variation patterns that depend on the geometric configuration and also on the particular non-radial mode, labelled by the horizontal mode numbers $\ell$ and $m$ in the context of the spherical harmonics approximation. For eclipsing binary systems with similar stellar radii this can be significant (see Fig. 1 for examples), in extreme cases leading to the appearance of ”hidden modes” during the eclipses. As a consequence, in the absence of mode identification, without information on the modulation patterns, the proper disentangling of pulsations and eclipse flux variation becomes a challenge. Iterative procedures are involved to tackle the problem (e.g. Prša & Zwitter, 2005; Conroy et al., 2020). They are lengthy and cumbersome, but rewarding in the correct determination of the absolute parameters. For planetary transits with much smaller eclipsed regions, the modulations themselves are not expected to be of concern. Nevertheless we have investigated their magnitudes.

For the first of the scenarios mentioned above we assumed the presence of single period, radial pulsations. For the second scenario we assumed various multiperiodic, multimode, non-radial pulsations, resembling the ones expected in $\delta$ Scuti type variables. The corresponding synthetic time series were computed in all cases using pulzem, as described in the previous subsection.

Table 2 lists the properties of the pulsation signals used in the single radial mode case. These refer to the contribution of each pulsation mode to the integrated flux outside the transits, according to the formula

The forward modelling tool of pulzem scales the surface patterns and shifts their initial phases to yield the given parameters, before subjecting them to the eclipses.

We made use of the Transit and Light Curve Modeller code (TLCM; Csizmadia, 2020) to solve the combined synthetic light curves of the simulated stellar pulsations and the transits. The transit modelling in TLCM is done via the Mandel & Agol (2002) model, which is parameterized by the orbital period, $P$, the time of mid-transit, $t_{\text{C}}$, the star-to-planet radius ratio, $R_{\text{p}}/R_{\text{S}}$, the scaled semi-major axis, $a/R_{\text{S}}$, and the impact parameter, $b$. Stellar limb-darkening was taken into account with a quadratic limb-darkening formula, described by $\mu_{1}$ and $\mu_{2}$. TLCM also incorporates the wavelet-based routines of Carter & Winn (2009) for handling the time-correlated noise which can be made up of both instrumental and astrophysical effects (including flares, granulation, stellar pulsations, etc.). The noise is characterized by $\sigma_{r}$ (for the red component) and $\sigma_{w}$ (for the white component). The wavelet-based method for handling red noise of TLCM was tested in Csizmadia et al. (2021); Kálmán et al. (2022b) and was found to be consistent. In our case, the correlated noise is made up solely of stellar pulsations. Kálmán et al. (2022a) reported that the wavelet-based approach of TLCM can handle the $\delta$ Scuti type pulsations of WASP-33.

The uneven surface brightness of the host star due to the gravity darkening effect induced by its rapid rotation may yield asymmetric light curves in certain configurations (Barnes, 2009; Barnes et al., 2011; Ahlers et al., 2019). We used a slightly modified version of the gravity darkening model of TLCM (Lendl et al., 2020; Csizmadia et al., 2021), incorporating a full Roche-geometry (Wilson, 1979) which is used for the calculation of the exact shape of the stellar surface. The modelling takes the projected stellar rotational velocity ($v\sin I_{\star}$) into account when calculating the stellar surface potential. This is then used to describe the surface temperature variations in combination with the gravity darkening exponent $\beta$. There are two additional fitting parameters: the inclination of the stellar rotational axis $I_{\star}$ (measured from the line of sight) and the projected spin-orbit angle $\lambda$.

For every synthesized light curve we performed at least two series of fits. The first series corresponds to the input model configuration: no gravity darkening, aligned star configuration ($I_{\star}=90^{\circ}$, $\Omega_{\star}=90^{\circ}$ or $\lambda=0^{\circ}$). A second series assumed nonzero gravity darkening using a temperature gravity darkening coefficient $\beta=0.25$ and a fixed $v\sin I_{\star}=$ 86.5 km/s, and both stellar angles being fitting parameters. We note that with this choice the amount of gravity darkening from the polar to the equatorial zone is completely determined and thus its effect on the transit light curve is fully characterised by the two angles and the conjunction parameter $b$. Sporadically we also fitted the conjunction parameter in order to assess its influence on the solution in particular cases. The scaled semi-major axes, $a/R_{\text{S}}$, star-to-planet radius $R_{\text{p}}/R_{\text{S}}$, mid-transit time $t_{C}$, period $P$ were always fitted in our analysis as well as the red noise $\sigma_{r}$ and white noise $\sigma_{w}$. The limb-darkening parameters $u_{+}=\mu_{1}+\mu_{2}$ and $u_{-}=\mu_{1}-\mu_{2}$, as required by TLCM, were kept fixed.

For an objective selection of the best model we made use of the Akaike and the Bayesian Information Criteria ($AIC$ and $BIC$, respectively; Cavanaugh & Neath, 2019), defined as

where $n_{\text{par}}$ is the number of parameters, $n_{\text{obs}}$ is the number of observations, and $RSS$ is the sum of squared residuals. The only difference between the two formulae is that the penalty factor for the number of model parameters increases steeper for the BIC due to the $\log n_{\text{obs}}$ factor. The model with the lowest $AIC$/$BIC$ values is the most plausible, while the relevance of the other models is measured in terms of the difference of their given $AIC$/$BIC$ from the lowest one. If this difference $\Delta AIC$ or $\Delta BIC$ is smaller than 3, then the empirical support for the model is substantial. For differences between 4 and 7 the plausibility is considerably smaller, and above 10 the model is implausible.

Due to the pulsations being fitted in the form of red noise, we computed two sets for both criteria. In one of them the red noise is accounted for in the model, which is labelled ’yn’ (i.e. ’yes noise’). In the other one, the red noise is excluded from the model, instead it is considered just another addition to the residuals, hence it is labelled ’nn’ (i.e. ’no noise’).

In the case of pXX series the enhanced pulsation amplitude might imply a red noise component so large that it approaches the limit of analytic capability of TLCM. Whenever this happened, whe identified the frequencies using Period04 (Lenz & Breger, 2005), and subtracted them from the datased as pure harmonic signals. After having ensured that the prewhitening process is satisfactory, the TLCM analysis was repeated.

The presence of single frequency pulsations did not pose any difficulty for the wavelet algorithm of TLCM responsible for handling the noise. The original input parameters were successfully recovered within errors in all the 30 systems, as it is shown in Fig. 2 for a selection of cases. The increased scatter in the data translates into larger uncertainties in all the fitted parameters, as usual (left and middle panels of Fig. 2. The free axis models also yield essentially the aligned solution in the majority of cases, although with generally larger $AIC$ and always larger $BIC$ values, which reassures the higher plausibility of the simpler model. There are also small deviations in the obtained relative planet radii (left panel of Fig. 2), which however still fall well within their uncertainties obtained for the pulsation-free case. Fig. 3 shows the achieved fits to the data for three of them: i00, i10 and i11, together with the residuals shown separately for the red noise and the residuals respectively. Note that the aligned and free axis models, shown in blue and orange, are virtually indistinguishable in this setup.

Fig. 4 illustrates the modulations experienced by the integrated flux of some nonradial modes during the transits in our artificial system. The simulated light curves reveal very small apparent modulations of the ”transited” pulsations for all cases. When compared to ordinary eclipsing binaries (Fig. 1 in Sec. 2), the difference is especially striking, because the amplitude variations due to the transits are not larger than about 2-5 percents of the unperturbed value, and the phase variations are also of a mere 2 degrees at most. By contrast, typical amplitude modulations in eclipsing binaries stars are 20-50 percents, while the phase variations commonly reach 30-50 degrees. When this phenomenon is applied to realistic scenarios of transiting exoplanets, the modulations are much smaller than the measurement errors, therefore virtually undetectable.

For completeness we also investigated the dependence of the modulations of the same nonradial modes on the $\lambda$ projected spin-orbit angle while keeping $I_{\star}$ at 90^∘. We obtained an almost unnoticeable amplification of the modulations. Although we did not consider all the possible oblique configurations, the small modulations suggest that they would not be significant for any other case either.

Therefore it seems that the distortions in the modelled parameters due to pulsations do not depend on their radial or nonradial nature. As a downside, neither are they useful in assessing any obliquity of the host star. As outlined in the previous section, single-mode pulsations are not a problem. Multiple pulsations on the other hand may still lead to a notable distortion of the transit curve profile if they add up under an unfortunate constellation. We have explored this possibility in the nxx series of trial fits. Table 4 and Fig. 5 summarise the investigated cases and their results.

The first case, n00, yielded the correct configuration for the gravity-darkened model. Fig. 6 shows the achieved fits to the dataset. But a second case, n01, with similar randomly picked modes, gave a side tilted configuration, as did case n02 too (see Fig. 7), where all the modes have been included. However the associated ’yn’ $\Delta BIC$ differences are in all cases large enough to signal a significantly lower plausibility of the gravity-darkened models with respect and in favour of the simpler model with fixed angles.

Nevertheless, it is still instructive to explore the nature of distortions induced by the pulsation signals that governed the solution away from the aligned configuration. To do so, we made a series of trial runs in order to identify modes that contribute significantly to the related asymmetry in the transit curve. Four suspect modes were identified, corresponding to frequencies 2, 5, 6 and 10 of Table 3. They have been selected based on the resemblance of their combined contribution to the difference in the fitted curves of the gravity-darkened and the undarkened models, that discrepancy being responsible for the tilted configuration achieved by the former model. Details of the individual as well as the combined contributions of these frequencies are shown in Fig. 8.

The other suspicious candidate is F1, for the obvious reason of being very close to $30\cdot f_{\text{orb}}$ and therefore potentially liable for causing a distortion that survives the folding and binning procedure. To test our hypothesis, we composed three test cases, omitting either the four suspect modes (n03), or the resonant mode (n05), or both (n04). As presented in Fig. 5, case n03 showed no improvement regarding the model selection, the other two, n04 and n05 became decisive in the $BIC$ values. While the AIC values, which are within the treshold of rejection for the case n04, may suggest that the free axis model has a slightly higher probability of being the real solution, the BIC values and the discovered aligned axis positions indicate that this solution is not stable.

A comparison of Table 4 and Fig. 5 reveals that side configurations arise only when the resonant mode F1 is present in the data. We also emphasise that the ’yn’ BIC numbers render for all cases the gravity-darkened, free angles models less plausible than the simple, aligned angle models with no gravity darkening.

We note here that $(I_{\star},\Omega_{\star})$, $(I_{\star},180^{\circ}\pm\Omega_{\star})$ and $(-I_{\star},360^{\circ}-\Omega_{\star})$ are all equivalent solutions, because photometry does not distinguish between prograde and retrograde rotations (Barnes et al., 2011). These equivalent configurations clearly show up on Fig. 5 as the double peaks in the ”violin plots” illustrating the marginal distributions of the angle $\Omega_{\star}$.

Regarding the pxx model with enhanced pulsation , we found that the limit of capability for the red noise algorithm in TLCM is reached, as the resulting solution exhibits significantly larger uncertainties compared to nxx (see Tab 6, leftmost column). The ratio of the uncertainties of pxx to the typical ones from previous runs is around 10; for example, the value and uncertainty of $\Omega_{\star}$ = $2\begin{subarray}{c}+147\\ -148\end{subarray}$. The lightcurve solution and the original data is shown in the leftmost panels of Fig 9.

Due to this issue, we applied the prewhitening process with Period04, and after identifying the significant frequencies, we subtracted them as a simple harmonic from the dataset, as can be seen in Fig. 10. We note that there are some residuals after removal due to the simple model.

As presented in the middle and right columns of Tab. 6, the solution obtained using the gravity darkening model exhibits esentially aligned configuration ($I_{\star}$ = $93\begin{subarray}{c}+15\\ -20\end{subarray}$; $\Omega_{\star}$ = $81\begin{subarray}{c}+32\\ -15\end{subarray}$). In any case, the AIC and BIC values also reassure that the simpler model is the correct solution. For completeness, the lightcurve solutions are displayed in the middle and right panels of Fig. 9 as well.