© The Authors 2023

1 Introduction

The large amount of exquisite quality photometric data furnished by the Kepler and Transiting Exoplanet Survey Satellite (TESS) space missions (Borucki et al. 2010; Ricker et al. 2015), following the pioneering CoRoT satellite (Auvergne et al. 2009), has led to the discovery of thousands of exoplanets via the transit method, later confirmed by ground-based observations using the radial velocity method. These combined data have also allowed for the determination of the absolute physical parameters of the planets, once the properties of the host stars become known with a sufficient level of accuracy. There is now a large database that has become available, with different exoplanet radii and masses and on various orbits, which has allowed for detailed statistical studies to be carried out regarding their distribution and architecture (see Zhu & Dong 2021, for a review).

Given that the inferences on planetary radii and masses require a characterisation of their host stars, at least in terms of mass and radius, there are several ongoing, dedicated systematic efforts to achieve this goal. Asteroseismology has proven particularly useful both in giant (e.g. Chaplin & Miglio 2013) and solar-type stars (e.g. Silva Aguirre et al. 2015) via analyses of solar-type oscillations, the regular frequency pattern of which allows for the derivation of asteroseismic masses and radii with a high level of precision. The case of other pulsating star types is not yet so promising. Studies of δ Scuti stars, which could be the next logical step in extending asteroseismic studies, are still challenging, despite the tremendous progress in terms of available data, also thanks to the space telescope missions mentioned above (Guzik 2021). Although recent progress has been achieved in the case of early type δ Scutis (Bedding et al. 2020) and also in some of their regular counterparts (Murphy et al. 2021), the majority still defy a definitive asteroseismic solution. They have irregular, complex frequency patterns, caused by overly large, uneven rotation splittings due to their fast rotation and the as-yet unexplained excitation ‘rules’ of their modes. Additionally, the pulsation amplitudes have a broad range, namely, from 370 ppm to approximately 30 000 ppm in relative flux (Uytterhoeven et al. 2011). Stars of spectral types A/F, to which δ Scutis also belong, are generally excluded from radial velocity exoplanet surveys (Howard et al. 2010; Grandjean et al. 2020; Rosenthal et al. 2021) because their heavily rota-tionally broadened spectral lines make it challenging to obtain precise RV measurements. However, although not particularly favored in dedicated transit surveys too, they may still show up in subsequent transit searches among space mission data. Exoplanets found around such hot stars have the potential of verifying the theory of photoevaporation of primary planetary atmospheres under the harsh light of their hot young stellar hosts (Owen & Wu 2013, 2017), conceived to explain the ‘radius gap’ in the distribution of Kepler planetary radii (Fulton et al. 2017). Therefore, exoplanet searches directly targeting δ Scuti stars (e.g. Hey et al. 2021) can be extremely fruitful.

There are also transiting systems showing asymmetric light curves. The first of this kind was Kepler-13 (Szabó et al. 2011, 2012), explained as a combined effect of gravity-darkened stellar surface due to fast rotation and a tilted planetary orbit. In recent years, exoplanet-hosting stars exhibiting δ Sct pulsations and gravity darkening (WASP-33; Collier Cameron et al. 2010; Herrero et al. 2011; von Essen et al. 2014; Dholakia et al. 2022; Kálmán et al. 2022) were also investigated. At the same time, intrinsic stellar activity, such as spots, flares, pulsation, or granulation may also yield asymmetric transits. Studies inspecting modelled high-precision photometric light curves of exoplanet systems with an active host star showed that the determined transit parameters – such as transit time, transit depth, and transit width – are influenced by the presence of stellar spots (e.g. Oshagh et al. 2013). Moreover, Oshagh et al. (2016) investigated the effect of the polluting signal of stellar activity on the spectroscopic determination of such quantities using the Rossiter-McLaughlin effect, showing that in some cases, it is capable of distorting the fit and falsely rendering an aligned system as misaligned. Pulsations could also constitute a problem, especially when their amplitudes are too low to be properly detected and are separated from the light curve. It is true that when multiple individual transits are combined into an average transit curve, the same process also tends to suppress any periodic signal that is not in resonance with the planetary motion. However, we note that there are frequencies close to an orbital resonance that may have an asymmetric footprint on the transit curve when folded with the orbital period. Another concern could be the modulation phenomenon of the non-radial pulsations subject to transit. This phenomenon can be very pronounced in ordinary eclipsing stars ((e.g. Kim et al. 2002; Rodríguez et al. 2004, 2010; Maceroni et al. 2014). In case of transits, the modulation is expected to be much lower due to the much smaller size of the eclipsing body, but this could still cause distortions if the effect of multiple modes adds up to the level capable of affecting parameters that depend very sensitively on the minute details of the transit curve. Parameters of this kind include the spin–orbit angles, for example.

To mitigate the effects of irregular stellar activity, these nuisance signals can be handled as a time-correlated noise, in a wavelet-based approach (Carter & Winn 2009), as in the Transit and Light Curve Modeller (TLCM; Csizmadia 2020; Csizmadia et al. 2023). On the other hand, this approach is not entirely guaranteed to work in the case of multi-periodic pulsation signals because the red component is not assumed to give rise to such asymmetries.

The goal of our paper is the exploration of the effect of small-amplitude, multi-periodic, pulsations on the determination of planetary system parameters from transit light curves. We consider two scenarios for each of them. The first scenario features single-period, radial-only pulsations. In the second scenario, we examine the effect of the modulation of the multi-periodic, nonradial pulsations on the parameters, handled either as radial mode harmonics subjected to eclipses or as red noise. For each case, we ran a full light-curve analysis, aiming to identify those configurations that are particularly vulnerable to such polluting signals. In Sect. 2, we describe the methods used in this investigation. Section 3 presents the results and their interpretations. In Sect. 4, we discuss them in a broader context, regarding the implications on future transit light curve analysis procedures. In Sect. 5, we summarise our findings.

2 Methods

We considered a hypothetical exoplanet-star system with a pulsating host star without rotational distortion, thus exhibiting uniform equilibrium surface brightness with no gravity darkening and no particular symmetry axis – along with various pulsation scenarios, both with single period radial and multi-periodic nonradial pulsations. Artificial light curves were synthesized using two separate tools: one for computing the transits of the time-independent, uniform stellar surface, and a second one for doing the same for the time-dependent pulsating surface-brightness component. The two variations were added together and supplemented with an artificial noise. Then we used a third modelling tool to fit both the gravitationally undarkened and darkened transit models on the data. The latter models also included the spin–orbit angles as fitting parameters. We then identified the cases leading to misaligned configurations. All the steps are presented in detail below.

2.1 Synthetic light curves

We studied the effects of a transit on the pulsational pattern in a generic model based on the pixel-level in-silico visualisation of the pulsating stellar surface while it was being 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 were able to identify which of these are the most sensitive to planetary transits, namely, the modes in which 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 Bíró & Nuspl 2011; Pál 2012), 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_s, and conjunction parameter, b: $$\omega =\frac{2\pi }{P}\frac{a}{R_{\text{s}}}\frac{1}{\sqrt{1-b^2}}.$$(1)

The photometric dataset was modelled with a quadratic limb-darkening law with coefficients of µ₁ = 0.209 and µ₂ = 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 as well. 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 precisely taken into account. The usual features of limb darkening, rotational distortion, and gravity darkening of the host star were 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; that is, each side of the icosa-hedron was divided into 15 segments, giving 15² = 225 tiles for each of the 20 faces and resulting in 20 · 15² = 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 nonlinear pulsation patterns, so the stellar surface actually requires no rotation at all. The same limb darkening used for the static light curve was applied to the pulsation surface patterns as well, before generating the combined pulsation flux of all modes.

The reason behind using one tool for synthesizing the static flux and another for the pulsating flux component is that currently there is no suitable tool available capable of computing both components with the same accuracy. In particular, fitsh is virtually exact for transit light curves, but cannot submit arbitrary time-dependent surface patterns to the same transit computation, while pulzem computes both, but its accuracy is intrinsically limited by the rasterisation 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 percent of the transit depth render these errors below the applied artificial uncorrelated noise component. Therefore, we deemed the use of this combined simulation method appropriate.

The first investigated scenario was a test system consisting of a spherical star and a planet with a relative size of R_p ~ 0.1 R_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, ~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° was used, corresponding to the WASP-33 system, which ultimately inspired this work.

2.2 Modelling the pulsations

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\left(\theta ,\phi \right)$, with a degree, ℓ, 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 non-radial modes of the pulsations. During eclipses (or transits, in exo-planetary 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 ℓ 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), and in extreme cases, can lead to the appearance of ‘hidden modes’ during 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 chose to investigate 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 multi-periodic, multi-mode, non-radial pulsations, resembling the ones expected in δ 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: $$f\left(t\right)=A\cdot \sin \left(2\pi \left(ft+\phi \right)\right).$$(2)

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

Fig. 1 Modulations of pulsations by the eclipse computed in a hypothetical binary system. Both stars are spherical, the orbit is circular, the inclination is 85°, and the secondary radius is half of the primary radius. The left panel shows with blue lines the eclipse footprint of the secondary on the circular disc of the primary. Only the footprints corresponding to every fifth point of the light curves on the right panels are shown for clarity. The right panels show the actual light changes for two selected modes: (0,0) on top and (2,1) on bottom.

2.2.1 Single-mode, radial pulsations

We used the case of single, radial pulsations to test whether handling them as red noise is appropriate, namely, whether the fitting procedure still yields parameters close to their input values. Red noise was introduced in TLCM primarily to cope with the signatures of stellar activity, for which it proved to be excellent (Kálmán et al. 2023).

2.2.2 Multi-periodic, non-radial pulsations

We also tested how much bias can be expected from non-radial pulsations in various geometrical configurations. It is evident that the leading parameter is the relative exoplanet radius: larger radii influence larger areas on the stellar surface and thus cause larger modulations. Therefore, an upper limit can be determined for the possible detection.

Generally, we wanted to map the variation of the light-curve modulations with the alignment angle of the exoplanet orbit. Therefore, we created a sequence of systems by varying the projected spin-axis angle λ between 0 and 180 degrees, in steps of 15 degrees. This also changes the projection of the surface pulsation patterns on the apparent stellar disc with respect to the transit trajectory and also the type of modulation as a consequence. We assume that the symmetry axis of the pulsations coincides with the axis of rotation, as usual.

We note that pulzem offers an option to compute the instantaneous amplitude and initial phase of each mode in every moment, then output them directly as amplitude and phase modulation data. These quantities are mode number-specific, but do not depend on other characteristics of the pulsation; therefore, we used one fiducial pulsation with fixed amplitude, frequence, and initial phase for all modes. We modelled different (ℓ, m) modes in each system, including the radial mode (0,0), which serves as the base modulation relative to which all the others are evaluated. Even the amplitude of a radial mode suffers a slight decrease due to the transit obscuring part of the disc, but no phase modulation takes place, as the whole surface oscillates coherently.

As we progressed with our study, we aimed to analyse a more realistic model containing multiple frequencies. The amplitudes attributed to the modes are small compared to typical δ Sct pulsations, because they were inspired by the WASP-33 system; for this work, we assumed some incidental (ℓ, m) non-radial mode numbers (there is no known mode identification presented in literature for this system). The summary of frequencies and their inclusions in the nXX model series in this paper can be found in Table 3. In our final set of simulations and analyses, we utilised the same frequencies as in the case of n02 but increased their amplitudes by a factor of 10. These simulated scenarios are labelled as pXX.

2.3 Performing fits with TLCM

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_C, the star-to-planet radius ratio, R_p/R_S, the scaled semi-major axis, a/R_S, and the impact parameter, b. Stellar limb-darkening was taken into account with a quadratic limb-darkening formula, described by µ₁ and µ₂. The code 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 characterised by σ_r (for the red component) and σ_w (for the white component). The wavelet-based method for handling red noise of TLCM was tested in Csizmadia et al. (2023); Kálmán et al. (2023) and was found to be consistent. In our case, the correlated noise is made up solely of stellar pulsations. Kálmán et al. (2022) reported that the wavelet-based approach of TLCM can handle the δ 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. 2023), 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 (υ sin I_*) 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 β. There are two additional fitting parameters: the inclination of the stellar rotational axis, I_* (measured from the line of sight), and the projected spin–orbit angle λ.

For every synthesised light curve, we performed at least two series of fits. The first series corresponds to the input model configuration: no gravity darkening and aligned star configuration (I_* = 90°, Ω_* = 90° or λ = 0°). A second series assumes nonzero gravity darkening using a temperature gravity-darkening coefficient of β = 0.25 and a fixed υ sin I_* = 86.5 km s⁻¹, with both stellar angles as fitting parameters. We note that with this choice, the amount of gravity darkening from the polar to the equatorial zone is fully 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 fit the conjunction parameter in order to assess its influence on the solution in particular cases. The scaled semi-major axes, a/R_S, star-to-planet radius, R_p/R_S, mid-transit time, t_C, and period, P, were always fitted in our analysis, as well as the red noise σ_r and white noise σ_w. The limb-darkening parameters u₊ = µ₁ + µ₂ and u₋ = µ₁ − µ₂, 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 $$\text{AIC\hspace{0.17em}=\hspace{0.17em}}2\cdot n_{\text{par}}+n_{\text{obs}}\cdot \log \left(\text{RSS}/n_{\text{obs}}\right),$$(3) $$\text{BIC = }{n}_{\text{par}}\cdot \log n_{\text{obs}}+n_{\text{obs}}\cdot \log \left(\text{RSS}/n_{\text{obs}}\right),$$(4)

where n_par is the number of parameters, n_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_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 ∆AIC or ∆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 a 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, we identified the frequencies using Period04 (Lenz & Breger 2005), and subtracted them from the dataset as pure harmonic signals. After having ensured that the pre-whitening process was satisfactory, the TLCM analysis was repeated.

3 Results

3.1 Single-mode, radial pulsation

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 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 supports 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. Figure 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. We note that the aligned and free axis models, shown in blue and orange, are virtually indistinguishable in this setup.

Fig. 2 Summary of the results for selected iXX models. Blue colours mark the aligned case with β = 0 and magenta colours mark the free angles cases with β = 0.25 in all graphs. Left panel: relative deviations of a/RS and Rp/RS fits from their input values. Middle panel: fitted angles Ω* and I* free angles case assuming β = 0.25. Right panel: computed differences, ∆AIC and ∆BIC, between the free angles and fixed angles cases. Filled circles mark the ‘yn’ scenarios (red noise considered) and empty circles mark the ‘nn’ scenarios (no red noise). The two grey shaded zones correspond to differences of 5 and 10, respectively, representing two thresholds of rejection.

Fig. 3 Sample fits and residuals of three simulated systems from the iXX series. Top panel: folded light curves with grey points, with simulated transit curve samples for the first orbital period, including noise and pulsation with black circles. There are also the fits of the two models: the fixed angle, β = 0 model with blue dashed line, and the free angles, β = 0.25 model with orange solid line. The two fits are virtually indistinguishable on these graphs. Middle panel: residuals in the red noise case. Bottom panel: residuals in the white noise case. The colours along the horizontal line at y = 0 encode the sign of the difference between the predicted fluxes of the fixed angles (λ = 0, I* = 90; at β = 0) and free angles (β = 0.25) models. Positive difference is shown in red, negative in blue; otherwise it is magenta.

3.2 Multi-mode, non-radial pulsations

Figure 4 illustrates the modulations experienced by the integrated flux of some non-radial 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), the difference is especially striking, because the amplitude variations due to the transits are not larger than about 2–5% 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%, 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 and, therefore, virtually undetectable.

For completeness, we also investigated the dependence of the modulations of the same non-radial modes on the λ projected spin-orbit angle, while keeping I_* 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 non-radial nature. As a downside, they are not useful in assessing any obliquity of the host star either. 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 explored this possibility in the nxx series of trial fits. Table A.1 and Fig. 5 summarise the investigated cases and their results.

The first case, n00, yielded the correct configuration for the gravity-darkened model. Figure 6 shows the achieved fits to the dataset. Howeve, 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’ ∆BIC differences are in all cases large enough to signal a significantly lower plausibility of the gravity-darkened models with respect to 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, with 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 · f_orb and therefore potentially liable for causing a distortion that survives the folding and binning procedure. To test our hypothesis, we built 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 A.1 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-angle models less plausible than the simple, aligned-angle models with no gravity darkening.

We note here that (I_*, Ω_*), (I_*, 180° ± Ω_*), and (−I_*, 360° − Ω_*) 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 Ω_*.

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 Table A.3, left-most 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 ${\text{Ω}}_{\star }=2_{-148}^{+147}$. The light-curve solution and the original data are shown in the left-most 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 Table A.3, the solution obtained using the gravity-darkening model exhibits an essentially aligned configuration ($I_{\star }={93}_{-20}^{+15}$; ${\text{Ω}}_{\star }={81}_{-15}^{+32}$). In any case, the AIC and BIC values also reassure that the simpler model is the correct solution. For completeness, the light curve solutions are displayed in the middle and right panels of Fig. 9 as well.

Fig. 4 Example for modulations of some non-radial pulsations during the exoplanet transit phases in a system with aligned orbit (i = 90°, λ = 0°). The four panels on the left show the residuals in the signals after the subtraction of mode (0,0). The panels on the right show the variations in the amplitudes (top) and phases (bottom) of all visible non-radial modes. Note: the small vertical scale of the phase modulation diagram implies that such a change will be barely noticeable.

4 Discussion

Our primary goal has been to investigate the possible disturbing effects of pulsations of the host star on the determination of the parameters of a transiting exoplanet. We carried out the investigation using trial fits on synthetically generated data, contaminated artificially with various types of pulsation signals, with amplitudes exaggerated deliberately to provide a clear study of the effect.

We found that the modulations of non-radial pulsations in the integrated flux caused by the transits are so small that they would be unnoticeable in TESS observations. Thus, the non-radial nature of pulsations can be safely ignored; a radial treatment involving only the diminishing effect of the transits on the amplitudes is sufficient for the transit analysis.

The presence of single frequency pulsations did not deter the fitting in any of the investigated cases. The gravity-darkened, free-angle models with β = 0.25 yielded the aligned configuration and, essentially, the input parameters within errors. The simple models with β = 0 models also restored the same parameters. In addition, the lower AIC and BIC values of the latter indicate that gravity darkening is an unnecessary assumption that is not supported by the data. These criteria become even more decisive when the red noise is considered as a component of the model, confirming that this approach is useful in handling sinusoidal contaminating signals too.

We found, however, a number of multi-periodic, small amplitude pulsation signal combinations which deceived the modelling in the gravity-darkened case (e.g. for n02; see Fig. 7 and Table A.2), and yielded side tilted solutions with Ω_* ~ 0°. Strangely, intermediate values between 0 and 90 degrees have never occurred. Also, the stellar inclination essentially stayed close to the plane of the sky, with I_* ~ 80–85°. In this configuration, the ~5–10° deviation of the rotation axis from the plane of the sky mainly affects the profile depth, while the side tilt accounts for small asymmetries in the transit curve. It seems that the pulsation signals mainly affect the latter feature of the transit.

The investigation of the resulted misaligned solutions indicates that they are associated with distortions in the transit light curve caused by the presence of frequencies close to orbital resonance. The other group of frequencies had no such influence. Their indirect resonance is caused by the finite length of the modelled dataset, so the folded and binned light contribution could be different based on another number of observations, containing either more or fewer transits. Therefore, this effect depends not only on the frequency of a mode, but on the volume of the dataset as well. That said, their removal from the dataset prior to modelling made the AIC/BIC-based plausibility evaluation more certain.

Nevertheless, when it comes to model selection, the ∆AIC and ∆BIC differences generally point to the model without gravity darkening as the more plausible one. This is especially true for the ‘yn’ variant of ∆BIC, because its lowest values consistently followed the correct solutions in all model scenarios investigated in our work. This resonates well with the fact that the wavelet method for estimating the correlated (red) noise component is an integral part of TLCM, therefore the ‘yn’ variants of both AIC and BIC are the proper choice for comparing the different models.

The simulated pulsations should be considered as the lower amplitude regime within typical δ Sct pulsations. Therefore, we also investigated the case of more typical δ Sct-like pulsations in terms of amplitude. Our findings suggest that the enhanced amplitude reduces the capabilities of TLCM due to the large proportion of red noise in data, making it essential to remove the problematic signals. In our case, all frequencies have been cleared as simple harmonics only, which resulted in larger residuals during the transit. Despite this fact, the wavelet algorithm of TLCM easily handled the situation, and the aligned axis model was accepted as the correct solution. We note that if other frequencies were present alongside the original ten frequencies, the additional residuals could cause problems, but addressing this issue is beyond the scope of this paper.

Some exoplanetary systems around δ Scuti type host stars, such as WASP-33b (Kálmán et al. 2022) and KELT-26b (Rodríguez Martínez et al. 2020), have highly misaligned configurations, some of them close to perpendicular. In case of KELT-26b a spurious frequency close to the 18th multiple of the orbital frequency was reported, however, it may be unrelated to our findings, given that the misalignment has been confirmed with Doppler tomography. The possibility that the pulsations, especially the non-radial ones, could also affect the Doppler tomography data in a similar way cannot be excluded, but such an investigation is clearly beyond the scope of the present paper.

Based on our findings, we propose extreme caution with a simultaneous eclipse or transit modelling that handles a suspected pulsation signal as red noise. Our suggestion is that transit analysis could benefit from a cleaning of the data of the possible problematic frequencies beforehand, in addition to preparing to conduct a comparison of the various models using their BIC numbers – even in cases when proper spectroscopy measurements and analysis are available for a better determination of the angles.

Fig. 5 Summary of the results for the nXX cases. The structure of the figure is similar to Fig. 2, but on the left and middle panel the marginal distribution of the corresponding parameter is depicted with so called violin plot. The circles and bars are belonging to the found median values and the uncertainty ranges.

Fig. 6 Model fits achieved for test case n00. The models are shown in the column titles. The top row shows the synthetic data with black symbols and the modelled curves with blue lines. The middle and bottom rows show the fitted red noise and the final residuals, respectively. The coloured areas show the differences between the fits of the actual model and the leftmost model (i.e. aligned, fixed angle). Red and blue correspond to positive and negative discrepancy, respectively, and the regions with no difference are coded with magenta.

Fig. 7 Model fits achieved for test case n02. The structure and notations are the same as in Fig. 6.

Fig. 8 Contributions of some of the sinusoidal signals to the total flux around the transit phases. Top left: contributions of F2 around the individual transit events with blue symbols and shifted vertically from each other for a better view, overlaid on the folded transit curve. Bottom left: contribution of F2 after folding with the orbital period and binning, shown on a scale of 100 ppm. The right panels show the contributions of F2, F5, F6, and F10, as well as their sum.

5 Conclusion

In this work, we performed a successful analysis investigating the possible degeneration of pulsations and gravity darkening in the photometric modelling of exoplanet transit light curves. We constructed dozens of synthetic photometric light curves of increasing complexity with pulzem and fitsh and performed the fits with the TLCM in several pre-defined runs. Our results show that the amplitude and phase modulations caused by the transit shadow on the pulsating disc are negligible in the light curve. We found that the noise handling algorithm of the TLCM generally copes well with datasets containing single-mode pulsations, aligned positions were preferred in the false assumption of gravity darkening as well. This conclusion remains generally valid in the presence of multiple frequency signals, with some notable exceptions, however. According to our findings, frequencies close to resonance with the orbital frequency may escape the noise handling, cause distortions in the data, and thus lead to false results, especially regarding the spin–orbit angles. This can be avoided by carrying out a cleaning of the problematic frequencies prior to transit analysis and comparing the plausibility of tilted models with fixed angles models. In this respect, we found that the most reliable indicator for choosing between the various fitted models is the BIC ‘yn’ number, which handles the red noise modelling as an integral part of the model. Due to the noise-handling wavelet algorithm, there is an upper limit of applicability in the amplitude of the signals. In such a case of higher amplitude pulsations, it is better to clean the pulsations beforehand. A similar conclusion has been reached by other studies investigating the case of asymmetric transit fits due to the interference of signal coming from stellar activity (Oshagh et al. 2016).

Fig. 9 Summary of the results for the pXX cases. The structure of the figure is similar to Fig. 7, except that the left panel corresponds to a data series having all the frequencies, while on the other two columns show the ten detected frequencies were removed as sine model. Therefore, the differences between the models are only given for the latter.

Fig. 10 Illustration of the pre-whitening process on the data in the pxx series. The top and bottom rows correspond to the data before and after pre-whitening, respectively. Left panels: Fourier transform of the full dataset, shown both with the inclusion and exclusion of the transit phases (brown and green colours, respectively). Note: the amplitude of the prewhitened spectrum (bottom left) is shown with a ten-fold increase. The frequencies detected by Period04 are shown by filled or empty circles, depending on whether the signal is present or not in the data. Right panels: corresponding folded synthetic dataset (light gray), as well as the data of the first orbital cycle (black).

Acknowledgements

We would like to thank the anonymous referee for their suggestions and comments to improve this paper. This work was supported by the PRODEX Experiment Agreement No. 4000137122 between the ELTE Eötvös Loránd University and the European Space Agency (ESA-D/SCI-LE-2021-0025). Support of the Lendület LP2021-9/2021 and the Lendület LP2018-7/2022 grants of the Hungarian Academy of Science is acknowledged. 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 Tables for multi-mode, non-radial pulsations

The tables below present the obtained values from the analysis of n00 (Table A.1), n02 (Table A.2), pxx (Table A.3), and all the nxx (Table A.4).

References

All Tables

All Figures

Fig. 1 Modulations of pulsations by the eclipse computed in a hypothetical binary system. Both stars are spherical, the orbit is circular, the inclination is 85°, and the secondary radius is half of the primary radius. The left panel shows with blue lines the eclipse footprint of the secondary on the circular disc of the primary. Only the footprints corresponding to every fifth point of the light curves on the right panels are shown for clarity. The right panels show the actual light changes for two selected modes: (0,0) on top and (2,1) on bottom.

In the text

Fig. 2 Summary of the results for selected iXX models. Blue colours mark the aligned case with β = 0 and magenta colours mark the free angles cases with β = 0.25 in all graphs. Left panel: relative deviations of a/RS and Rp/RS fits from their input values. Middle panel: fitted angles Ω* and I* free angles case assuming β = 0.25. Right panel: computed differences, ∆AIC and ∆BIC, between the free angles and fixed angles cases. Filled circles mark the ‘yn’ scenarios (red noise considered) and empty circles mark the ‘nn’ scenarios (no red noise). The two grey shaded zones correspond to differences of 5 and 10, respectively, representing two thresholds of rejection.

In the text

Fig. 3 Sample fits and residuals of three simulated systems from the iXX series. Top panel: folded light curves with grey points, with simulated transit curve samples for the first orbital period, including noise and pulsation with black circles. There are also the fits of the two models: the fixed angle, β = 0 model with blue dashed line, and the free angles, β = 0.25 model with orange solid line. The two fits are virtually indistinguishable on these graphs. Middle panel: residuals in the red noise case. Bottom panel: residuals in the white noise case. The colours along the horizontal line at y = 0 encode the sign of the difference between the predicted fluxes of the fixed angles (λ = 0, I* = 90; at β = 0) and free angles (β = 0.25) models. Positive difference is shown in red, negative in blue; otherwise it is magenta.

In the text

Fig. 4 Example for modulations of some non-radial pulsations during the exoplanet transit phases in a system with aligned orbit (i = 90°, λ = 0°). The four panels on the left show the residuals in the signals after the subtraction of mode (0,0). The panels on the right show the variations in the amplitudes (top) and phases (bottom) of all visible non-radial modes. Note: the small vertical scale of the phase modulation diagram implies that such a change will be barely noticeable.

In the text

	Fig. 5 Summary of the results for the nXX cases. The structure of the figure is similar to Fig. 2, but on the left and middle panel the marginal distribution of the corresponding parameter is depicted with so called violin plot. The circles and bars are belonging to the found median values and the uncertainty ranges.
In the text

Fig. 6 Model fits achieved for test case n00. The models are shown in the column titles. The top row shows the synthetic data with black symbols and the modelled curves with blue lines. The middle and bottom rows show the fitted red noise and the final residuals, respectively. The coloured areas show the differences between the fits of the actual model and the leftmost model (i.e. aligned, fixed angle). Red and blue correspond to positive and negative discrepancy, respectively, and the regions with no difference are coded with magenta.

In the text

	Fig. 7 Model fits achieved for test case n02. The structure and notations are the same as in Fig. 6.
In the text

Fig. 8 Contributions of some of the sinusoidal signals to the total flux around the transit phases. Top left: contributions of F2 around the individual transit events with blue symbols and shifted vertically from each other for a better view, overlaid on the folded transit curve. Bottom left: contribution of F2 after folding with the orbital period and binning, shown on a scale of 100 ppm. The right panels show the contributions of F2, F5, F6, and F10, as well as their sum.

In the text

	Fig. 9 Summary of the results for the pXX cases. The structure of the figure is similar to Fig. 7, except that the left panel corresponds to a data series having all the frequencies, while on the other two columns show the ten detected frequencies were removed as sine model. Therefore, the differences between the models are only given for the latter.
In the text

Fig. 10 Illustration of the pre-whitening process on the data in the pxx series. The top and bottom rows correspond to the data before and after pre-whitening, respectively. Left panels: Fourier transform of the full dataset, shown both with the inclusion and exclusion of the transit phases (brown and green colours, respectively). Note: the amplitude of the prewhitened spectrum (bottom left) is shown with a ten-fold increase. The frequencies detected by Period04 are shown by filled or empty circles, depending on whether the signal is present or not in the data. Right panels: corresponding folded synthetic dataset (light gray), as well as the data of the first orbital cycle (black).

In the text