Based on TESS observations of a delta Scuti variable WASP33 obtained in 2019 and 2022 we thoroughly investigate the power spectrum of this target photometric flux and construct a statistically exhaustive model for its variability. This model contains 30 robustly justified harmonics detected in the both TESS sectors simultaneously 13 less robust harmonics detected in a single TESS sector a red noise and a quasiperiodic noise terms. This allowed us to greatly improve the accuracy of the exoplanet WASP33 b transit timings reducing the TTV residuals r.m.s. drastically by a factor of 3.5 from 63 s to 18 s. Finally our analysis does not confirm existence of a detectable orbital phase variation claimed by von Essen et al. 2020 based on WASP33 TESS photometry of 2019.
['Astrophysics - Earth and Planetary Astrophysics', 'Astrophysics - Instrumentation and Methods for Astrophysics', 'Astrophysics - Solar and Stellar Astrophysics']
Cleaning WASP33 b transits from the host star photometric variability analysis of TESS data from two sectors
{'Cleaning WASP-33 b transits from the host star photometric variability: analysis of TESS data from two sectors': 'Roman V. Baluev 1 , 2 /star and Eugene N. Sokov 3  \n- 1 Saint Petersburg State University, 7-9 Universitetskaya Emb., St Petersburg 199034, Russia\n- 2 Special Astrophysical Observatory, Russian Academy of Sciences, Nizhnii Arkhyz, 369167, Russia\n- 3 Central Astronomical Observatory at Pulkovo of Russian Academy of Sciences, Pulkovskoje sh. 65/1, St Petersburg 196140, Russia  \nAccepted 2024 .... Received 2024 ...; in original form 2024 ...', 'ABSTRACT': 'Based on TESS observations of a δ Scuti variable WASP-33 obtained in 2019 and 2022, we thoroughly investigate the power spectrum of this target photometric flux and construct a statistically exhaustive model for its variability. This model contains 30 robustly justified harmonics detected in the both TESS sectors simultaneously, 13 less robust harmonics detected in a single TESS sector, a red noise and a quasiperiodic noise terms. This allowed us to greatly improve the accuracy of the exoplanet WASP33 b transit timings, reducing the TTV residuals r.m.s. drastically, by a factor of 3 . 5, from 63 s to 18 s. Finally, our analysis does not confirm existence of a detectable orbital phase variation claimed by von Essen et al. (2020) based on WASP-33 TESS photometry of 2019.  \nKey words: techniques: photometric, stars: variables: Scuti, stars: individual: WASP33, planets and satellites: detection', '1 INTRODUCTION': "WASP-33 (HD 15082) is a bright 8.14 mag fast-rotating star with v sin i ∼ 90 km/s, located in Andromeda. Its the only known exoplanet WASP-33 b was detected by SuperWASP (Herrero et al. 2011). Its orbital period is estimated by P b = (1 . 219867 ± 4 . 2 · 10 -5 ) d (von Essen et al. 2020), mass by M b = (2 . 10 ± 0 . 14) M Jup , and radius by r b = (1 . 59 ± 0 . 07) R Jup (Chakrabarty & Sengupta 2019). It is remarkable that this exoplanet was the first one ever known to orbit a δ Scuti variable.  \nThe variability of the prototype star itself, δ Scuti, was discovered in the beginning of XX century (Campbell & Wright 1900), and then it was confirmed by spectral and photometric observations (Colacevich 1935; Fath 1935). Furthermore, Eggen (1956) combined this and a few similar stars together in a standalone class of variables, and so the δ Scuti type of variable stars appeared.  \nThese variables are pulsating stars belonging to spectral classes A0-F5 III-IV, located in the instability strip of the HR diagram (Baglin et al. 1973; Breger & Stockenhuber 1983). Their m V variation amplitudes range from 0 . 003 to 0 . 9 (typically a few hundreds) and periods from 0 . 01 d to 0 . 2 d. Most of the δ Scuti stars belong to type I population (Breger 1979), with a typical mass of ∼ 2 M /circledot (Milligan & Carson 1992). The power spectrum of δ Scuti  \nvariables may reveal tens of periods (Breger et al. 1999) ranging from hours to minutes. Such stars often exhibit significant changes in the form of their lightcurve, in their variability period and amplitude. Radial as well as non-radial pulsations can be observed. Nonradial pulsations indicate concurrent inward and outward surface moves, contrary to radial ones which reflect expansions and contractions of the star as a whole. A key role in the pulsation generating mechanism belongs to helium which is highly abundant in the atmosphere of these stars. When helium is heated and ionized it gains opacity, causing the star brightness be at minimum. However, the radiation absorbed makes helium to heat further and expand, causing its recombination, transparency increase and the star brightness to reach the maximum. Then helium starts to cool down and the cycle restarts again (Cox 1963).  \nLarge amount of data regarding δ Scuti variables were brought by (i) OGLE and MACHO surveys which detected about 3000 such variables in the Large Magellanic Cloud (Poleski et al. 2010; Garg et al. 2010), (ii) Kepler spacecraft (Borucki et al. 2010) that found about 750 candidate variables of this type, and (iii) TESS (Transiting Exoplanet Survey Satellite) space mission, for example (Antoci et al. 2019) analyzed a sample of 117 variables belonging to δ Scuti and γ Doradus class.  \nThe intrinsic variability of WASP-33 represents an obvious nuisance factor for analysing its exoplanet b transit lightcurves. The effect of this variability on the estimated  \nplanetary orbital parameters was previously investigated using out-of-transit observations from ground-based telescopes located in Germany and Spain (von Essen et al. 2014). Their investigation of the WASP-33 pulsations power spectrum revealed eight significant frequencies. It was concluded that cleaning the transit lightcurves from stellar pulsations practically does not affect estimated planetary parameters, but their uncertainties get decreased. Besides, von Essen et al. (2014) noticed that WASP-33 pulsations phases change over time. In a later analysis of TESS photometry of 2019, von Essen et al. (2020) claimed the detection of 29 pulsation frequencies. In addition to the power spectrum analysis, they also reported statistically significant detection of WASP-33 b secondary eclipse and of the orbital phase curve variation (period of P b , amplitude of (100 . 4 ± 13 . 1) ppm). Also, they took into account the 'ellipsoidal variation' that appears because of planetary gravitational effect on the star shape (period of P b / 2, theoretically predicted amplitude of 27 . 4 ppm).  \nA more recent analysis of WASP-33 variability was performed using data from TESS (K'alm'an et al. 2022). They investigated the Fourier spectrum of stellar oscillations reconstructed with TLCM (Transit and Light Curve Modeller, Csizmadia et al. 2023). This spectrum revealed three peaks with S/N > 8 that were close to a commensurability with the planetary orbital period (1 : 3, 1 : 12, and 1 : 25). It was argued that the distribution of 17 frequencies with the largest amplitudes tend to shift right from these subharmonics. As K'alm'an et al. (2022) noticed, this effect can be caused by the tidal perturbations from the planet.  \nOur primary goal here is to analyse WASP-33 pulsations spectrum based on the up-to-date TESS data (observations of 2019 and 2022), and to investigate whether (and how much) this can help us to improve the transit timings accuracy for the exoplanet WASP-33 b .", '2 TESS PHOTOMETRY OF WASP-33': "TESS observed WASP-33 in sector 18 (Nov 2019) and sector 58 (Nov 2022), lasting about 1 month each. In the DVT data downloaded from the MAST (Mikulski Archive for Space Telescopes) we found 14796 and 18951 photometric measurements for these two sectors, respectively. We identified 17 transits in sector 18 and 21 transits in sector 58, and for each transit we extracted a ± 4 h piece of the lightcurve about the predicted midtransit epoch. This subset formed our 'in-transit' dataset (ITD). The 'out-of-transit' dataset (OTD) included all observations located at least ± 2 h away from the centers of transits and secondary eclipses. Notice that ITD and OTD slightly overlap.  \nIn this work we will deal mainly with the OTD, since we aim to construct a model for the intrinsic star variability. The sector 18 portion of the OTD will be hereafter refered to as TESS 18 , and it contains 10684 datapoints. The sector 58 OTD portion (TESS 58 ) contains 13625 datapoints. ITD contains 8914 datapoints in total, with ∼ 200 observations per transit.  \nTESS data look close to being even with a constant timestep of ∆ t = 2 min. However, this is not strictly fulfilled even within just a single sector, because in the middle of each sector there was a gap, after which the 2 min grid  \npattern got shifted. Also, the OTD set contains holes in place of transits and secondary eclipses. But nevertheless the spectral window functions of TESS 18 , 58 data suggest that these data can be treated as nearly even, in particular there is a classic aliasing pattern suggesting the existence of a 'quasiNyquist' frequency f Ny = 1 / (2∆ t ) = 360 d -1 (see Fig. 2).", '3 BUILDING THE STATISTICALLY EXHAUSTIVE MODEL OF WASP-33 VARIABILITY': "Here our goal is to construct a statistically 'exhaustive' WASP-33 model that represents its OTD so that the residuals do not contain statistically significant periodic (or quasiperiodic) variations. In this attempt we aimed to apply a periodogram-based approach somewhat reminiscent of the CLEANest algorithm by Foster (1995). Namely, we represent the photometric variation as a sum of sinusoidal harmonics:  \nµ n ( t ) = A 0 + n ∑ k =1 A k cos (2 πf k ( t -T 0 ) + ϕ k ) , (1)  \nwhere A k , ϕ k , and f k = 1 /P k are free parameters.  \nHowever, if the model (1) was used in the white noise (WN hereafter) treatment then we would likely need hundreds of harmonics to represent the OTD exhaustively. After extraction of some initial (the most significant) harmonics, the residual periodogram did not reveal a single clearly dominating peak, looking more similar to a spectrum of nonwhite (banded) noise.  \nAlone such an observation is not yet decisive to claim that the data physically involve non-white noise (NWN hereafter). Two alternative models 'a lot of harmonics + WN' and 'fewer harmonics + NWN' may appear formally indistinguishable, i.e. statistically equivalent in how well they describe the data. But even in this case the NWN model may offer a great reduction in the number of model components (and parameters), so we can profit from a mathematically more simple µ .  \nTo take the NWN into account, we adopted an approach of Gaussian process (GP) fitting, introduced in (Baluev 2011). In this approach we should define a model for the NWN correlation function, which we now set to:  \nk n ' (∆ t ij ) = ( σ 2 TESS ,i + σ 2 0 ) δ ij + n ' ∑ k =1 σ 2 k e -| ∆ t ij | τ k cos 2 πf ' k ∆ t ij , ∆ t ij = t i -t j , (2)  \nwhere σ TESS ,i is TESS-provided photometric uncertainty at t i , δ ij is Kronecker delta, while σ k , τ k , and f ' k = 1 /P ' k are free parameters. This model involves the WN portion ( σ 0 ) and n ' independent quasiperiodic components (purely NWN part). By a visual look at the periodograms (e.g. see Fig. 3 discussed below) we noticed three possible spectrum bands that are concentrated near the frequencies f ' 1 /similarequal 0, f ' 2 /similarequal 10 d -1 , and f ' 3 /similarequal 20 d -1 . We set these frequencies as starting values for the GP fitting, thus assuming n ' = 3.  \nThen we run a sequential periodogram algorithm according to the following scheme. We started from k 3 as specified above, and from µ 2 involving two harmonics with fixed frequencies as per von Essen et al. (2020): f 1 = P b (orbital phase variation) and f 2 = P b / 2 (ellipsoidal variation). After  \nFigure 1. Spectral window of the TESS 18 (left) and TESS 58 (right) data, for three frequency ranges (from top to bottom). The (quasi)Nyquist alias peak at 2 f Ny = 720 d -1 is seen in the bottom panels.  \n<!-- image -->  \nthat the model (1) can be built up iteratively by adding harmonics that refer to the tallest peak in the subsequent residual periodogram. These periodograms must take the NWN into account via (2). For that we used the formalism of the likelihood-ratio periodogram (LRP) introduced by Baluev (2009, 2013b) and implemented in the PlanetPack software (Baluev 2013a, 2018). The LRP considers two rival models of the data, a base one ( µ n on n th iteration) and an alternative one involving an extra sinusoid (i.e. µ n +1 ). The likelihood ratio statistic is then computed for these models and its logarithm is plotted against the frequency of the last added sinusoid. Such a graph represents our LRP that generalises the classic Lomb (1976) - Scargle (1982) periodogram. The LRP can be computed for the classic WN case as well as for the NWN model (2), by defining the likelihood function appropriately. After each iteration all free parameters of so far detected harmonics are refined by refitting with the new model µ .  \nVery importantly, the spectral window in Fig. 2 does not demonstrate large side peaks (in the Nyquist range at least). The height of the tallest side peak is /similarequal 0 . 1. Therefore, the aliasing appears insignificant in our task, and we may safely assume that the maximum periodogram peak corresponds to the true (rather than alias) frequency. This makes all the analysis easier, for example it is unnecessary to care about the order in which the harmonics are extracted (which could be important otherwise).  \nHowever, dealing with NWN periodograms implies nu-  \nmerous NWN model refitting, which is computationally demanding. The likelihood function involves operations on very large covariance matrices ( ∼ 10 4 × 10 4 ), rendering the LRP computation too slow for being accomplished in full. But we tried to bypass this limitation using the following approach. On each iteration we computed, at first, the fast WN-only version of the LRP, and determined its 256 tallest peaks. Then for each of these peaks we run an NWN maximum-likelihood fitting, starting it from the associated WN best fit. This allowed us to compute heights of these same peaks in the NWN treatment, and then to select the tallest one among them.  \nSuch an approach is potentially vulnerable because it may, in theory, miss the actual best-significance harmonic. Even 256 peaks that appeared tallest with the WN model do not necessarily contain the one that should appear tallest in the NWN case. However, the probability of such a condition decreases when the number of peaks being inspected grows, so the validity of the results can be tested by a probe increase of that number (in the end of iterations sequence).  \nWe launched our harmonics detection iterations assuming the frequency range from f = 0 to f = 50 d -1 , which is a priori reasonable for δ Scuti stars, and also contained all apparently noticeable periodogram power. As we specified above, the NWN model (2) was initially set to have n ' = 3 terms, but somewhere in the middle of the sequence the third NWN term (with f ' /similarequal 20 d -1 ) appeared, basically, decomposed into harmonics. We started to obtain various  \ndegeneracy issues regarding this term, for example it was trying to either model just a single sinusoid ( τ → 0), or to mimic a portion of the white noise ( τ →∞ ), or to model a portion of another NWN term. When this behaviour started, we manually removed this third term from the NWN model, and continued the sequence further with k 2 (∆ t ) instead of k 3 (∆ t ). Finally, when no significant harmonics left (in the NWN treatment), we reprocessed the last iteration inspecting 1024 tallest peaks instead of just 256, and in the entire Nyquist range (max frequency of 360 d -1 instead of 50 d -1 ). This did not result in detection of additional harmonics, ensuring that we did not miss any.  \nThis algorithm was applied to TESS 18 and TESS 58 data separately, and gave similar results. Some graphs illustrating the initial and final iterations are shown in Fig. 3. Each panel demonstrates the WN periodograms (raw as well as smoothed one), the theoretic power spectrum of the current best fitting NWN model based on (2), and the set of NWN peaks that were dealt with on this iteration. On the last iteration, when we inspected 1024 peaks, those peaks are covering the entire Nyquist range without significant holes or blind segments, so it is very unlikely that some statistically significant harmonics escaped our analysis. Please also notice remarkable similarity between our fitted NWN power spectrum and the smoothed WN periodogram.  \nSome explanation should also be given regarding how we computed the statistical significance of periodogram peaks in the NWN treatment. As explained in (Baluev 2009), for the WN type of the LRP we may use the following approximation of the False Alarm Probability (FAP) from (Baluev 2008):  \nFAP( z ) z /greatermuch 1 /similarequal W e -z √ z, W = ( f max -f min ) T eff , (3)  \nwhere z is the height of the tallest peak, and T eff is an effective time span proportional to the sample variance of the timings t i . As argued (without proofs) in (Baluev 2015), the approximate formula (3) remains applicable for NWN models as well, but with T eff somewhat differently determined. Delisle et al. (2020) developed this study further, partly confirming this conclusion and providing a more general definition for T eff . However, in our task it appeared that T eff always remains practically identical, regardless of the noise model or which approximation to adopt. We obtained T eff about 24 d for TESS 18 and about 28 d for TESS 58 .  \nIn this study we aimed to obtain a more accurate final model of the WASP-33 variability, possibly even by the cost of larger probability of false positive harmonics. So we adopted a relatively mild FAP threshold of 5% (the iterations continued until the right hand side of (3) exceeded 0 . 05).  \nTable 1 lists all the harmonics that we detected in TESS 18 and in TESS 58 . 1 It is highly important that both these sectors revealed very similar sets of the harmonics. Nearly all harmonics can be cross-identified between TESS 18 and TESS 58 by comparing their frequencies (honouring uncertainties). Usually, the amplitudes of the cross-identified harmonics appeared similar as well, although some of A k  \ndemonstrated big changes between TESS 18 and TESS 58 . It is impossible to conclusively compare the phases of these harmonics, because of too large uncertainty accumulated over three years passed between these TESS sectors. In general, one may conclude that WASP-33 variability pattern demonstrates remarkable, though not entirely perfect, stability over the timescale of a few years at least.  \nAll harmonics in Table 1 are indexed in the periodincrease order as H k . Some of them, that are also indexed as F k , were cross-identified with those given by von Essen et al. (2020).  \nIn Fig. 3, top panels, we plot all the detected harmonics in the frequency - amplitude plane (with uncertainties). Harmonics that were cross-identified by their frequencies are additionally labelled with circles. We can see that in a given TESS sector just a few harmonics miss a counterpart in the another one. Such cases might indicate that the given harmonic is a false positive or vice versa, that the counterpart harmonic failed to overcome the FAP = 5% threshold level.  \nIn the bottom panels of Fig. 3 we plot the cumulative function of the harmonics power A 2 k / 2 against the amplitude A k . From these graphs we can see that ∼ 70% of the total power (variance) of the compound harmonics signal is contained in just 4 or 5 components of largest amplitude ( A k /greaterorsimilar 600 ppm), while the remaining ∼ 30% of the total power is distributed among the rest ( A k /lessorsimilar 400 ppm). The total power of the compound signal corresponds to ∼ 1200 ppm for TESS 18 and TESS 58 both.  \nTable 2 contains the best fitting NWN parameters of the model (2). The total noise variance corresponds to ∼ 700 ppm and appears nearly the same for TESS 18 and TESS 58 . Within it, the second NWN term appears nearly the same, although the first one reveals a significant change between the sectors, as well as the WN term does. But because of the stable total variance, this difference looks like a power exchange between the WN and the first NWN term over time. The total model variance becomes therefore ∼ 1400 ppm for the both TESS sectors, and the variance of these data themselves matches this level closely, confirming that our models are consistent and close to being statistically exhaustive.  \nFinally, in Fig. 3 we show the power spectrum of the compound model 'harmonics + NWN', as well as its correlation function. The power spectrum is normalized so that it can be directly compared with the periodograms in Fig. 3. The NWN has a continuous power spectrum, but harmonics are expressed by delta functions there, so we plot them as vertical arrows starting from the NWN graph. In the correlation function plots we show separate graphs for the WN and NWN portions, and for the harmonics. Notice that correlation function of a single sinusoid in (1) is A 2 k 2 cos 2 πf k ∆ t .", '4 IMPROVING TRANSIT TIMING ACCURACY': "Now let us consider how our out-of-transit photometric model of WASP-33 can be used to detrend the transit lightcurves. The most direct and robust way to achieve this is to fit the entire TESS dataset (OTD and ITD) with this model, augmented by the models of all the transits. However, this way it would require too slow nonlinear fits, and besides we have a goal in mind to employ our model with  \nFigure 2. Periodogram analysis and extraction of the harmonics.  \n<!-- image -->  \nTable 1. Estimated harmonics parameters.Table 2. Estimated noise parameters.  \nany other transit observations, so it is necessary to develop  \nan approach that avoids using full TESS data directly. Our model contains two components: the deterministic  \none µ ( t ) (harmonics) and the GP part (NWN) expressed by the correlation function k (∆ t ). Regarding the deterministic part, its uncertainty can be computed based on the covariance matrix of fitted parameters. This uncertainty remains in the range 100 ppm to 250 ppm (for TESS 18 and TESS 58 , resp.) Relatively to the average TESS-reported uncertainty for these data, ∼ 30000 ppm, this appears below 1%. Therefore, possible uncertainty in µ ( t ) is surely negligible in comparison with the best expected photometric accuracy. So we can simply subtract our best fitting µ ( t ) from the ITD, as if it was a predefined function. We will refer this part of the ITD processing as detrending for shortness.  \nRegarding the NWN part of the model, its smallest τ k from Table 2 are about 3-4 h, which is similar to the exoplanet transit duration. This makes it tentative to construct  \nFigure 3. All the detected harmonics in the frequency - amplitude plane and their contributions to the total cumulative power.  \n<!-- image -->  \nsome kind of a GP prediction model (based on the OTD subset) and then to subtract it from the ITD, just like µ ( t ). We undertook an attempt to compute an OTD-conditioned GP predictive model with parameters taken from Table 2. However, the practical predictability timescale appeared significantly smaller than τ k , namely below 1 h. Along the transit duration the variance of this prediction remained mostly close to the total variance of the NWN (i.e. to the unconditional variance, as if OTD were not used). Therefore, such a predictive model would not be very helpful, but simultaneously this means that we would not loose significant information if we do not use the OTD directly when processing the ITD. This enabled us to apply a simplified approach with maximum-likelihood fitting of the ITD adopting a fixed covariance matrix (which is defined by the NWN parameters from Table 2).  \nImportantly, this approach is equivalent to the whitening in which the data are passed through a linear transform that renders them white. 2 With our GP-fitting algorithm  \n2 If V = Var x is the covariance matrix of the data x , and r is the vector of their residuals then the whitening transform can be r ' = V -1 2 r . This implies Var r ' = I being the unit matrix, and | r ' | 2 = r T V -1 r , which forms the only variable part of the likelihood function when V is fixed.  \nit is unnecessary to apply whitening transform explicitly, but nonetheless we will refer this part of ITD processing as whitening.  \nSummarising, our ITD processing consists of two cleaning steps, detrending and whitening. In Fig. 4 we show how they work for a few sample transits. In these graphs we plot the TESS data before processing (green points), the sumof-harmonics µ (orange curve), the detrended ITD (black points, shifted down) together with transit model (orange curve), and detrended and whitened ITD, also together with their model (which is also whitened). One can see that detrending can remove a large portion of the host star variability, however some residual nuisance variation still remains and should be removed by whitening.  \nThe effect of this cleaning on transit timings is demonstrated in Fig. 4, where we plot the TTV residuals derived from the initial (raw) ITD, detrended ITD and detrended and whitened ITD. One can see that the TTV scatter is remarkably reduced after each step. The initial TTV residuals r.m.s. was 63 s. After detrending it became 36 s, while after the whitening it finally dropped to 18 s. In total we have a rather dramatic r.m.s. improvement by a factor of 3 . 5. The last best-accuracy set of transit timings is given in Table 3.  \nAfter that we investigated our improved timing data for possible variations, but with a null result: we did not find  \nFigure 4. Power spectrum and correlation function of the compound model 'harmonics+NWN'.  \n<!-- image -->  \nTable 3. WASP-33 b transit timings derived from the TESS ITD.  \nany significant quadratic term in the TTV trend and we did not find any significant periodic TTV.  \nFinally, fitting the detrended and whitened ITD allowed us to obtain an improved set of WASP-33 system parameters. They are given in Table 4. This fit used the stellar mass estimate by Collier Cameron et al. (2010), but ignored their stellar radius since we need to know the correlation between M /star and R /star to correctly include them both. Our estimate for  \nthe stellar density appeared somewhat smaller though (by ∼ 2 . 5 σ ).  \nWe treated the limb darkening coefficients as fittable parameters, but for a comparison in Table 4 we also give their reference values (i) computed from the tables by Claret (2017) using T eff , log g , and [ Fe/H ] by Collier Cameron et al. (2010) and (ii) computed by von Essen et al. (2020). These reference sets both differ from  \nFigure 5. Detrending and whitening for a few WASP-33 b transits. Notice that whitening normalizes the data to unit variance, therefore changing their scale.  \n<!-- image -->  \nTable 4. Best fitting WASP-33 system parameters.  \n∗ Epoch BJD TDB = 2459000.  \n∗∗ Taken from (Collier Cameron et al. 2010).  \nFigure 6. WASP-33 transit timing residuals before and after ITD cleaning. Notice the different ordinate range in the plots.  \n<!-- image -->  \nthe fitted one, as well as from each other. This may indicate some issues regarding their reliability for not yet clear reasons.", '5 DISCUSSION': "The approach based on the explicit modelling of WASP33 variability appears extremely successful in improving its exoplanet transit timings accuracy. In total, we obtained a dramatic TTV r.m.s. improvement by a factor of 3 . 5. However, it is important that we must take into account the deterministic as well as non-deterministic (NWN) parts of this variability. Both of them incorporate a similar amount in the total TTV error budget.  \nOur results still cannot be directly applied to observations from other sites than TESS, because the model parameters somewhat change on the timescale of a few years. Because of that we cannot use TESS-derived model with data obtained e.g. years prior to TESS. However, the frequencies of the harmonics and of the QP part of the NWN look stable, so it may be possible to apply our model to ground-based WASP-33 photometry assuming that only these frequencies are fixed and to refit harmonics amplitudes and phases, as well as NWN variances. But possible ways to reprocess old data using TESS-based models need more investigation.  \nOur results regarding the set of frequencies in the de-  \nterministic part of the WASP-33 variability are generally similar to von Essen et al. (2020). We found 43 harmonics in total, although only 30 of them are surely reliable because they were revealed in the both datasets, TESS 18 and TESS 58 . von Essen et al. (2020) detected 29 harmonics in TESS 18 , which is nearly the same number. However, there are a few individual differences. Based on our analysis we cannot confirm that their harmonics F21, F25, and F29 are statistically significant, even in a single TESS sector. On contrary, our harmonics H22, H32, H37, and H40 were revealed both in TESS 18 and TESS 58 , but left undetected by von Essen et al. (2020). Among our harmonics that were revealed in TESS 18 , but not in TESS 58 , only one (H3 = F28) appeared in von Essen et al. (2020). Among harmonics that were detected in TESS 58 , but not in TESS 18 , none appeared in von Essen et al. (2020). This may render our 13 'singlesector' harmonics less reliable, although this does not necessarily mean that all of them, or their majority, are false positives. Based on our FAP threshold of 5%, the expected number of false positives after 36 or 37 applications of the significance test should be just one or two (for each TESS sector). And these false positives are, most probably, among the 'single-sector' harmonics.  \nRegarding the secondary eclipses of WASP-33 b , detected by von Essen et al. (2020), we cut them away from our OTD, so they did not affect our analysis and we did not consider them here. But two other planet-related effects re-  \nported by (von Essen et al. 2020), namely the orbital phase variation and the ellipsoidal variation, could affect our results, so we included them in the model µ from the beginning of our analysis. For the orbital phase variation our final amplitude estimate was (110 ± 100) ppm and (32 ± 44) ppm (resp. in TESS 18 and TESS 58 ), while von Essen et al. (2020) gave (100 . 4 ± 13 . 1) ppm. We can see a large disagreement in the uncertainty, most probably related to the different NWN treatment by von Essen et al. (2020) which we believe was oversimplified. Based on our analysis we cannot confirm that there is a statistically significant evidence in favour of the orbital phase variation in TESS observations. For the ellipsoidal variation we obtained amplitude estimations of (65 ± 41) ppm and (97 ± 31) ppm, for the two TESS sectors. von Essen et al. (2020) gave 27 . 4 ppm in the text, although we think this could be a typo, because from their Fig. 9 this amplitude follows to be ∼ 40 ppm. In any case, their value was theoretically predicted rather than estimated from the data directly. The value of ∼ 40 ppm is within at most two-sigma from our estimations. As we can see, in our analysis this effect comparable to the uncertainty and cannot be reliably confirmed observationally.", 'DATA AVAILABILITY': 'The data underlying this article are available in the MAST (Mikulski Archive for Space Telescopes), at https://archive.stsci.edu/ , in the article and in its online supplementary material.', 'REFERENCES': '```\nAntoci V., et al., 2019, MNRAS, 490, 4040 Baglin A., Breger M., Chevalier C., Hauck B., Le Contel J. M., Sareyan J. P., Valtier J. C., 1973, A&A, 23, 221 Baluev R. V., 2008, MNRAS, 385, 1279 Baluev R. V., 2009, MNRAS, 393, 969 Baluev R. V., 2011, Celest. Mech. Dyn. Astron., 111, 235 Baluev R. V., 2013a, Astron. & Comput., 2, 18 Baluev R. V., 2013b, MNRAS, 429, 2052 Baluev R. V., 2015, MNRAS, 446, 1493 Baluev R. V., 2018, Astron. & Comput., 25, 221 Borucki W. J., et al., 2010, Science, 327, 977 Breger M., 1979, PASP, 91, 5 Breger M., Stockenhuber H., 1983, Hvar Observ. Bull., 7, 283 Breger M., et al., 1999, A&A, 349, 225 Campbell W. W., Wright W. H., 1900, ApJ, 12, 254 Chakrabarty A., Sengupta S., 2019, AJ, 158, 39 Claret A., 2017, A&A, 600, A30 Colacevich A., 1935, PASP, 47, 231 Collier Cameron A., et al., 2010, MNRAS, 407, 507 Cox J. P., 1963, ApJ, 138, 487 Csizmadia S., Smith A. M. S., K\'alm\'an S., Cabrera J., Klagyivik P., Chaushev A., Lam K. W. F., 2023, A&A, 675, A106 Delisle J. B., Hara N., S\'egransan D., 2020, A&A, 635, A83 Eggen O. J., 1956, PASP, 68, 238 Fath E., 1935, PASP, 47, 232 Foster G., 1995, AJ, 109, 1889 Garg A., et al., 2010, AJ, 140, 328 Herrero E., Morales J. C., Ribas I., Naves R., 2011, A&A, 526, L10 K\'alm\'an S., B\'okon A., Derekas A., Szab\'o G. M., Heged"us V., Nagy K., 2022, A&A, 660, L2\n```  \n```\nLomb N. R., 1976, Ap&SS, 39, 447 Milligan H., Carson T. R., 1992, Ap&SS, 189, 181 Poleski R., et al., 2010, Acta Astronomica, 60, 1 Scargle J. D., 1982, ApJ, 263, 835 von Essen C., et al., 2014, A&A, 561, A48 von Essen C., Mallonn M., Borre C. C., Antoci V., Stassun K. G., Khalafinejad S., Tautvaiˇsien˙e G., 2020, A&A, 639, A34\n```  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.'}