1. INTRODUCTION

Since August 2018, the Transiting Exoplanet Survey Satellite (TESS, ^{Ricker et al. 2015}) has been delivering high-precision photometric observations in a broad optical band (0.6 - 0.95 µm) for a large sample of bright stars from the southern and northern hemispheres. The wavelength coverage of TESS allows measurements of the combined reflected and thermally emitted planetary light as a function of longitude.

The exoplanet WASP-121b was discovered by ^{Delrez et al. (2016)}, with a period of ≈ 1.275 days, and it is one of the hottest transiting planets known to date. Its bright host F6-type star (V = 10.4), short orbital period, and inflated radius (α/R _s = 3.674, R _p = 1.865R _J ) makes it one of the best targets for investigating its atmosphere with various techniques. Moreover, due to its short orbital period, it is likely that WASP-121b is tidally locked to its host star (^{Daylan et al. 2021}), which makes it probable to have atmospheric features detectable in the averaged planetary flux (^{Showman & Guillot 2002}). Several studies have measured WASP-121b’s primary transit (when an exoplanet passes in front of its host star) (e.g., ^{Delrez et al. 2016}; ^{Evans et al. 2016}, ²⁰¹⁸). By using optical and near-infrared photometry, the depth of its secondary eclipse (i.e., when an exoplanet is occulted by its host star) was determined by ^{Delrez et al. (2016)}; ^{Kovacs & Kovacs (2019)}; ^{Garhart (2019)}. The dayside and nightside temperatures of WASP-121b were measured as 2870K and < 2200K, respectively, according to an analysis of the thermal emission (^{Bourrier et al. 2020}), which is to be expected given the planet’s proximity to its host F-type star. The reflection component was not included in ^{Bourrier et al. (2020)} phase curve, but we take it into account in our comprehensive full phase curve model. WASP-121b’s geometric albedo was estimated as 0.70-0.040+0.037 based on the optical phase curve analysis by ^{Daylan et al. (2021)}. In our analysis, the reflection component is also calculated simultaneously with other parameters to highlight the correlations between all of the constrained parameters.

The main objective of the current study is to learn more about the thermal emission and atmospheric structure of WASP-121b by performing our comprehensive joint model and by comparison of our results with previous studies like ^{Daylan et al. (2021)}; ^{Bourrier et al. (2020)}. To achieve this, we analyze the full-orbit optical TESS phase curve and model the primary transit, secondary eclipse, and four main components of the phase curve, which include tidal ellipsoidal distortion, thermal, and reflected emission of the planet. We also calculated rotational Doppler beaming and discovered that it is not significant given the precision of the light curves. We can determine the uncertainty and correlations among all constrained parameters using our comprehensive joint model, which allows us to extract information from all parameters at the same time.

Here, we describe our WASP-121b analysis by presenting our comprehensive phase curve model and comparing our findings to previous measurements. The paper is organized as follows; in § 2 we describe the observations and data reduction methods that were used. In § 3 we describe in detail the four different components that were used to characterise the phase curve of WASP-121b. In § 4 we present our joint model as well as the fitting procedure we employed to acquire our results. We provide our physical parameters derived from the TESS observations in § 5, and discuss our results providing a brief summary in § 6.

2. OBSERVATIONS AND DATA REDUCTION

Between the 8th of January and 1st of February 2019, the TESS camera 3 monitored WASP-121 (also known as TIC 22529346) throughout its sector number 7. The observation span was 24.5 days and included 18 primary transits of WASP-121b.

Photometric data were processed through the Science Processing Operations Center (SPOC) pipeline (^{Jenkins 2017}). In this study, we decided to use PDCSAP (Pre-Search Data Conditioning) light curves because they are corrected for instrumental systematic noise which is present in the Simple Aperture Photometry (SAP) light curves; thus PDCSAP light curves show considerably less scatter and short timescale flux variation (^{Smith et al. 2012}; ^{Stumpe et al. 2014}). The PDCSAP light curve of WASP-121 was also used in other studies investigating the phase curve of WASP-121b such as ^{Bourrier et al. (2020)}.

The PDCSAP photometry is presented in the upper panel of Figure 1, which shows the remaining systematics in the data at short time scales, particularly in Sector 7’s second orbit. Instrumental effects including changes in the thermal state of TESS and pointing instabilities cause these remaining systematics. The PDCSAP light curve’s median was used to normalise the data. To have a fair comparison with ^{Bourrier et al. (2020)}, we did exactly the same steps in preprocessing of data. Although the dominant systematics were corrected by default in the PDCSAP light curve, we corrected it further for the remaining systematics. To do this, we used the median detrending algorithm with a window length of one orbital period to smooth the PDCSAP light curve, keeping variability at the planetary period and minimising the effect of normalisation on the phase curve of WASP-121b. If we choose a smaller window length of one orbital period, then it is very likely that the signal will be absorbed and removed from the atmosphere. We followed the same processes using ^{Bourrier et al. (2020)}; the regression is shown in Figure 1 and was implemented using the Python package wotan (^{Hippke et al. 2019}). We also performed phase folding at the orbital period of WASP-12b and binned every 50 datapoints after detrending; these reprocessed data were used in our further analysis. The reprocessed data are shown in Figure 2.

Fig. 1 (Top) TESS light curve (PDCSAP flux) of WASP-121. The PDCSAP photometry is indicated with black dots, and the solid blue line shows the trend obtained by applying a detrending filter determined by wotan. (Bottom) PDCSAP light curve after normalization by its median detrending. The color figure can be viewed online. 

Fig. 2 (Top) Our reprocessed data of WASP-121 (blue dots) compared to our best fitted full phase curve model (black curve). (Middle) zoom of the secondary eclipse and phase curve variations with the reflection modulation (dashed orange curve), ellipsoidal distortion (dashed red curve), and baseline (dotted blue line). (Bottom) The best fitted model’s corresponding residuals. The color figure can be viewed online. 

3. PHASE CURVE

In addition to the primary transit light curve and secondary eclipse, photometric observations reveal additional variation induced by the orbiting exoplanet over the full planetary orbit. This variation can be decomposed into several components, namely thermal emission, reflected light, Doppler beaming, and ellipsoidal variation. In this study, we assume that the phase curve variation is a combination of thermal emission, reflection, and ellipsoidal variation, and we ignore the Doppler beaming; the reason will be explained in § 3.3.

4. MODEL AND FITTING PROCEDURE

Our joint model consists of primary transit, secondary eclipse, and phase curves that incorporate the thermal emission, reflection, and ellipsoidal variations. We also included a constant baseline to compensate for any normalization bias. For the primary transit and secondary eclipse we used the Python packages batman (^{Kreidberg 2015}) and for the thermal emission, we used spiderman (^{Louden & Kreidberg 2018}). Our thermal model is based on a semiphysical model of Zhang & Showman (see § 3.1) implemented by spiderman. The reflection is modeled as equation 1 (see § 3.2) and ellipsoidal variations are modeled as equation 7 (see § 3.4). Performing a joint model analysis allows us to extract information about all parameters simultaneously from the data sets. It also gives us the ability to assess the uncertainty and correlations between all of the constrained parameters.

To determine the parameters of the full phase curve, we fitted our joint model to the reprocessed data (see § 2). The best fit parameters and their associated uncertainties are determined using a Markov Chain Monte Carlo (MCMC) approach using the affine invariant ensemble sampler emcee package (^{Foreman-Mackey et al. 2013}).

We fit for R _p /R _s , α/R _s , i, u ₁, u ₂, ξ, T _N , ∆T _DN , additive baseline, secondary eclipse depth, A _ref , ∆ _P and A _ellip . The priors of u ₁, u ₂, ξ, T _N , and ∆T _DN are equal to those of ^{Bourrier et al. (2020)}. The priors of R _p /R _s , a/R _s , and i had normal priors in ^{Bourrier et al. (2020)}, and we chose an uninformative uniform prior for them. The additional parameters in our model have a wide uninformative prior, allowing us to obtain their best estimation.

Table 1 provides information on individual prior distributions that were chosen. We fix the transit epoch, T ₀ and orbital period, P because we use one sector of TESS data which covers about 24 days, whereas the period from ^{Delrez et al. (2016)} takes into account years of WASP data, which provides more information on the period. Considering the Lucy-Sweeney bias (^{Lucy & Sweeney 1971}), we adopt a circular orbit by fixing the eccentricity e, to 0, and the argument of periastron, ω to values obtained by ^{Bourrier et al. (2020)}. To generate the posterior distributions, we ran 700 walkers over 1000 steps with a burn-in phase of the 20% sample. The walkers are plotted and visually inspected for convergence. We estimated the median and standard deviation from the posterior distributions at 1σ, which contains 68% of the posterior distribution, for our best fitted values and uncertainties.

5. RESULT

The results of our joint model fitting to the reprocessed data are shown in Table 1. The best fitted model’s reduced chi-squared χ2 (i.e., RMS of the residuals per degree of freedom) is 1.29, indicating a good fit to the TESS photometry. Figure 2 shows the reprocessed data, as well as the best fitted model of the full phase curve. Our best-fitted model’s residuals still exhibit some correlated noise, which could be uncorrected TESS systematic noise. Identical correlated noise was also present in the residuals of the best-fitted ^{Bourrier et al. (2020)}. The corner plot for our retrieved posterior distributions from the joint model fit is shown in Figure 3. In addition, when we fit our joint model to unbinned data, we found that our results are generally consistent. Furthermore, we binned every 80 datapoints and fitted our joint model, and the results were consistent with the values presented in Table 1. This test demonstrates that the results provided in this work are robust to the binning effect.

Fig. 3 Retrieved posterior distributions by fitting our joint model to the phase curve of the WASP-121b. the black points indicate the best-fit values, and the colours of the contours highlight the 1 and 2σ simultaneous 2D confidence regions which contain respectively 39.3% and 86.5% of the samples. The solid black line corresponds to the median values, while dashed black lines show the 1σ highest density intervals. 

We calculated a planetary radius (in stellar radii), (R _p /R _s ), of 0.1234-0.0005+0.0005 and a reasonably large secondary eclipse depth with amplitudes of 489-10+16 ppm. Our measured secondary eclipse depth value is within 1σ of the value reported in the ^{Daylan et al. (2021)}. However, our estimated value is larger (1.6σ) than ^{Bourrier et al. (2020)} measured value. Other orbital and transit parameters agree well with previously published values in the literature (^{Wong et al. 2020}; ^{Daylan et al. 2021}; ^{Bourrier et al. 2020}).

Our estimation of the geometric albedo using equation 2 is 0.069-0.02+0.06 which is consistent with the estimate of ^{Daylan et al. (2021)} of 0.070-0.040+0.037 ppm. ^{Mallonn et al. (2019)} estimated a geometric albedo of A _g = 0.16 ± 0.11 in the z’ band.

We measured the ratio of radiative versus advection time scale of atmospheric height as ξ=-0.022-0.141+0.159 , which is consistent with zero. This implies that there is no thermal redistribution between WASP-121b’s day and night sides, resulting in a larger day-night temperature difference. We measured the temperature of the night and day sides, 2236-38+97 and 2941-61+147 , respectively, which are in agreement with the values published in ^{Bourrier et al. (2020)} and ^{Daylan et al. (2021)}. ^{Parmentier et al. (2018)} and ^{Evans et al. (2017)} by fitting the blackbody model to Spitzer and Hubble Space Telescope WFC3 observations could measure a dayside temperature of 2650±10K and 2700±10K, respectively.

Based on our best fitted model, we estimated the ellipsoidal variation amplitude to be 20-3+2 ppm, which is more in line with the theoretical estimate of 20 ppm based on equation 5 and slightly larger than ^{Daylan et al. (2021)} estimation which was 8-6+12 ppm.

The most remarkable result of our study is the simultaneous measurement of the primary transit, the secondary eclipse, and the robust detection of the total phase curve component corresponding to thermal emission, reflected light, and ellipsoidal variation (see Table 1). The three components and full-phase curve are plotted in the middle panel of Figure 2.

In our analysis, we also experimented what would happen if we let the eccentric e and ω free in our joint fit. In this case we obtained that these results are consistent with the values reported in Table 1 at about ≈ 1σ. We obtain eccentricity constraints: e=0.0024-0.0024+0.0041 And ω=9.05-1.06+2.32 deg which are consistent with those published in ^{Bourrier et al. (2020)}. According to Lucy-Sweeny bias (^{Lucy & Sweeney 1971}), in order to measure a non-zero eccentricity with 95% confidence, a result of e > 2.45σ _e is required, where σ _e is the standard deviation of the eccentricities (^{Eastman et al. 2013}). As a result, we can confidently rule out WASP-121b non-zero eccentricity.

In addition to our total phase curve model, we investigated a scenario in which the planetary flux is purely reflective. To approximate the planetary reflection, we used the Lambertian reflection model implemented in spiderman and characterized by a geometric albedo A _g . Using this scenario we estimated the geometric albedo to be 0.46-0.035+0.036 which is significantly (3σ) larger than the estimate of the geometric albedo reported by ^{Bourrier et al. (2020)}. It is quite close to Ag=δ(αp/Rp)2=0.47-0.03+0.03 which is the value derived from the TESS secondary eclipse depth δ=489-10+16ppm. The reduced χ ² of this purely reflective scenario is 1.8.

6. SUMMARY AND CONCLUSIONS

In this work, we presented our full phase curve model for analyzing the transiting ultra-hot Jupiter WASP-121b utilizing one sector of TESS observations. There were two reasons for using only one sector of TESS. The first is that different TESS sectors have different systematic noises, and combining several sectors may introduce additional complications and difficulties in our joint modelling. The second and most important reason is that we wanted to use the same data set as in ^{Bourrier et al. (2020)} and ^{Daylan et al. (2021)} so that we could assess how much improvement we could get from using a more complete model. We first used the median detrending technique with a window length of one orbital period of WASP-121b to conduct a smooth detrending of the TESS data, in order to have comparable data with ^{Bourrier et al. (2020)}, who also did the exact same steps. We binned every 50 data points after phase folding at the orbital period, as ^{Bourrier et al. (2020)} performed previously. In our subsequent analysis, we used these reprocessed data. Then we fitted our joint model to the reprocessed data. Our joint model consists of primary transit, secondary eclipse, and phase curves that incorporate the thermal emission, reflection, and ellipsoidal variations.

We reliably measured the secondary eclipse with a depth of 489-10+16 ppm after eliminating systematic noise. The combination of thermal emission and reflection in the TESS bandpass results in a relatively significant secondary eclipse depth of WASP-121b. Due to the strong stellar irradiation and low geometric albedo, the secondary eclipse is expected to be mostly dominated by the planet thermal emission.

Our measurement of the ξ=-0.022-0.141+0.159 is statistically consistent with zero. This value indicates that the atmosphere of WASP-121b has inefficient thermal redistribution from dayside to nightside, which is consistent with results in the literature (^{Bourrier et al. 2020}; ^{Daylan et al. 2021}) and with theoretical models (^{Komacek et al. 2017}; ^{Perez-Becker & Showman 2013}). WASP-121b’s maximum temperature region is located near the sub-stellar point due to inefficient thermal redistribution, as advection does not redistribute heat across longitudes (^{Zhang & Showman 2017}). The inefficient thermal redistribution also results in substantial differences in the night and dayside temperatures (734-55+28K) of WASP-121b. Our measured dayside temperature of 2941-150+61K for WASP-121b places it in the ultra-hot Jupiter class (^{Parmentier et al. 2018}; ^{Bell & Cowan 2018}).

In this study, we did not assume that the flux of WASP-121b measured by TESS was exclusively thermal emission, and we took into account reflected light. Our best fitted joint model yielded a low geometric albedo of 0.069-0.02+0.06 , indicating that reflection in the TESS passband of WASP-121b is not negligible, which was ignored by ^{Bourrier et al. (2020)}. Our estimated low geometric albedo value is consistent with ^{Daylan et al. (2021)} and other hot Jupiters, in particular, irradiated hot Jupiters at the same wavelength as ^{Schwartz & Cowan (2015)}. It is also consistent with other short-period hot Jupiter planets, such as WASP-18b (A _g < 0.048 at 2σ; ^{Shporer et al. (2019)}), Qatar-2b ( A _g < 0.06 at 2σ; ^{Dai et al. (2017)}), and WASP-12b (A _g < 0.064 at 97.5% confidence; ^{Bell et al. (2017)}). Considering the fact that the bandpass of TESS is close to the wavelength region where the host star is brightest, the bond albedo is small when the geometric albedo is small (^{Shporer et al. 2019}).

The amplitudes of the ellipsoidal variation and Doppler beaming are significantly smaller than those of reflected light and thermal emission, according to theoretical estimates (see Figure 2). We did not incorporate Doppler beaming in our phase curve model because our theoretical estimation of the amplitude of Doppler beaming yields a value of ≈ 2 ppm , which is not significant given the precision of the photometric data. Finally, our best fitted joint model also provided us with an estimate of the amplitude of the ellipsoidal variation, which is consistent with theoretical expectations.

The hot host star and the short orbital period of WASP-121b cause it to be highly irradiated. Furthermore, the lack of a statistically significant phase shift, poor heat distribution, and low albedo are all compatible with other highly irradiated giant gas planets. This study demonstrated that our model may be used to explore the full phase curves of transiting systems. The fact that the WASP-121b phase curve modulations were clearly detected shows that TESS data are sensitive to photometric variations in systems with short periods and massive planets.

More TESS data from extended missions or from other existing facilities like the CHaracterising ExOPlanet Satellite (CHEOPS) (^{Benz et al. 2021}) will also enable a more in-depth study of exoplanets’ full phase curve. Our WASP-121b retrieval analysis provides a glimpse into the comprehensive analysis of the full orbital phase curve which can be performed by combining optical and thermal infrared observations, near-infrared emission using existing facilities like the ARIEL (^{Tinetti et al. 2018}), and upcoming facilities with higher resolution, such as the James Webb Space Telescope (JWST) (^{Gardner et al. 2006}).