WASP-35 and HAT-P-30/WASP-51: re-analysis using TESS and ground-based transit photometry

WASP-35b was discovered by Enoch et al. (2011) in the Wide Angle Search for Planet (WASP) project. This transiting exoplanet system has a inflated ($M_{p}=0.72\pm 0.06M_{Jup}$ and $R_{p}=1.32\pm 0.05R_{Jup}$) hot Jupiter which orbits a metal-poor ($[Fe/H]=-0.15\pm 0.09$) host star with a period of 3.16 days. Through analyzing their photometric and spectroscopic observations at a series of ground-based telescopes, Enoch et al. (2011) confirmed the planetary nature of WASP-35b and suggested that the host star of WASP-35 was lack of stellar activity. Later, Mortier et al. (2013, 2014) updated the physical parameters of the host star using a spectroscopic analysis, Kokori et al. (2021) re-analyzed and updated the ephemeris based on the light curves of Enoch et al. (2011), and Shan et al. (2021) refined and updated the ephemeris using the data of TESS photometry.

HAT-P-30b (also known as WASP-51b) was discovered independently by RN518 in the Hungarian-made Automatic Telescope Network (HATNet) and Enoch et al. (2011) in the WASP project, respectively. HAT-P-30b is a transiting hot Jupiter ($M_{p}=0.711\pm 0.028M_{Jup}$ and $R_{p}=1.340\pm 0.065R_{Jup}$) orbiting a late F-type host star with a period of 2.81 days and the orbit is circular. HAT-P-30 has a relatively bright (V=10.4) host star, therefore it is very suitable for using transmission spectra to study the planetary atmosphere. RN518 studied the Rossiter-McLaughlin effect of HAT-P-30 system, and found the orbit of HAT-P-30b was highly tilted with a sky-projected angle between the star’s spin axis and the planet’s orbit normal, $\lambda=73.5\pm 9.0\arcdeg$. Maciejewski et al. (2016), Saha et al. (2021), Wang et al. (2021) and Edwards et al. (2021) acquired new transit light curves of HAT-P-30b based on several ground-based telescopes and measured the physical parameters. Davoudi et al. (2020) analyzed several relative high quality transit light curves of HAT-P-30b collected from the Exoplanet Transit Database (ETD) ¹¹1http://var2.astro.cz/ETD/. Kokori et al. (2021) re-analyzed and updated the ephemeris based on the light curves of Maciejewski et al. (2016), RN518 and Enoch et al. (2011), and found no TTVs for HAT-P-30b.

TESS is an all-sky space survey, which is designed to search for transiting exoplanets orbiting the bright and nearby stars, these planet systems are enable us to carry out the follow-up observations and study their physical properties. TESS performs time-series photometry to monitor at least 200,000 main-sequence stars using four 100 mm telescopes with wide-field optical CCD cameras. The broad passband filter of TESS covers from 600 to 1000 nm and the combined field of view of each sector is 24 $\times$ 96 square degrees. Since its launch in 2018, TESS has completed 44 sectors of observations and collected a lot of high-quality transiting light curves.

WASP-35 was observed by TESS in Sectors 5 (2018 November 15 - 2018 December 11) and Sectors 32 (2020 November 19 - 2020 December 17), and HAT-P-30 was observed by TESS in Sectors 7 (2019 January 7 - 2019 February 2) and Sectors 34 (2021 January 13 - 2021 February 9). The data were collected in a 2 minute cadence and reduced by the pipeline developed by the Science Processing Operations Center (SPOC; Jenkins et al., 2016). We downloaded the light curve files from the archives at Mikulski Archive for Space Telescopes (MAST)²²2https://archive.stsci.edu/mission/tess/ and accessed the light curves using the Lightkurve python package (Lightkurve Collaboration et al., 2018). These data had been corrected for the instrumental systematic variations by using the SPOC pipeline, but some long-term trends still remained in the data. In order to alleviate the influences from the remaining trends, we employed Gaussian Process (GP; Rasmussen & Williams, 2006) to model these data (see below for further details). At last, we obtained 15 complete transit light curves of WASP-35 and 16 complete transit light curves of HAT-P-30. The time system of these light curves is barycentric Julian date (BJD) which can be used directly to determine the system parameters and ephemerides. The examples of raw TESS light curves are displayed in the top panels of Figure 1 and Figure 2, respectively.

We observed the transit events of WASP-35 and HAT-P-30 using the Andor 4K$\times$4K CCD camera attached to the YO-1m on 2018 November 20 and 18, respectively. The Johnson-Cousins R filter was utilized in both observations and the field of view is 15$\times$15 arcmin². During two observations, both weather conditions and instrument statuses were good, and the exposure times were 120s. The transit event of HAT-P-30b occurred at sunset, which resulted in a relatively larger dispersion at the beginning of the observation (see Figure 3). In addition, we collected two sets of high quality light curves of HAT-P-30 from the CDS database, these light curves were observed by RN518 using the KeplerCam imager of the 1.2 m telescope at the Fred Lawrence Whipple Observatory (FLWO-1.2m).

The data reduction follows the standard procedures which were described in Wang et al. (2013). Our data reduction pipeline is based on the IRAF package, written in the python language by means of the pyRAF interface, including image checking, image trimming, bias subtraction, flat-field correction, cosmic ray removal, establish astrometric solution, aperture photometry and systematic error correction. To reach a higher precision, the pipeline attempts a series of photometric apertures for the target and reference stars and we pick out an optimal one to minimize the dispersion of each light curve. The transit light curve is obtained by using the optimal aperture.

For the ground-based observations, the shallow transit signals are easy to be drowned out by the systematic errors due to the effect of the telluric atmosphere and imperfect observation conditions, so it is necessary to remove the systematic errors hiding in above transit light curves. First, we employed the coarse de-correlation (Collier Cameron et al., 2006) and SYSREM algorithms (Tamuz et al., 2005) to diagnose and correct them by using the pipeline mentioned above. At this point, there are still some long-term systematic trends in the transit light curves, we employed GP to diagnose and correct them. Furthermore, in order to obtain accurate ephemerides of targets, we converted the observing time into barycentric Julian date by using the method proposed by Eastman et al. (2010).

Gaussian Process (Rasmussen & Williams, 2006) is a popular method which has been widely used to model time series data in the exoplanet community for about a decade, such as modeling stellar activity signals in RV data and correcting instrumentally induced systematic errors in transit light curves (Gibson et al., 2012; Evans et al., 2015; Murgas et al., 2020; Hurt et al., 2021; Langellier et al., 2021). Instead of modeling instrumental systematic errors as a deterministic function with auxiliary measurement parameters, GP provides a non-parametric method to model systematic errors from the observed data set. Here, we used the code of juliet (Espinoza et al., 2019) to corrected long-term trends and instrumental systematic errors in TESS photometry data and the residual systematic errors in the ground-based photometry data.

We modeled the photometric light curves in flux space separately depending on the instrument and the time of the observation by using juliet code. For the photometry data sets, juliet establishes a common model in which light curve is modeled as a linear combination of the transit model and the noise model:

where $\mathcal{M}(t)$ is the transit model with the dilution factor for the given instrument and the mean offset of out-of-transit flux, $\epsilon(t)$ is the noise model which is being modeled by GP in the photometric light curve. Juliet uses the analytic transit light curve model of Mandel & Agol (2002) with a quadratic limb-darkening law and the limb-darkening coefficients ($q_{1}$ and $q_{2}$) proposed by Kipping (2013) to model the transit signals and employs batman package (Kreidberg, 2015) to implement the transit light curve model.

We checked the nearby stars around WASP-35 and HAT-P-30 using Gaia DR3 database (Gaia Collaboration et al., 2021), both of the fields do not seem to be especially crowded, and there was no other source bright enough to make a perceptible impact, so we fixed the dilution factors of both of the targets to 1. For TESS photometry data, the state of the instrument may be different in different observation sectors, so we modeled the TESS data sector by sector, respectively. For the ground-based observations, the state of ground-based instrument and the weather is changing from day to day, so we modeled the ground-based data day by day, respectively. In all of the model processes, we chose a celerite (approximate) Matern multiplied exponential kernel within Juliet implemented by celerite package (Foreman-Mackey et al., 2017). For the analytic transit model part, the parameters are the orbital period $P$, the mid-transit time $T_{0}$, scaled semi-major axis $a/R_{A}$, the eccentricity $e$ (fixed to 0), the argument of periastron $\omega$ (fixed to 90 degrees), the impact parameter $b$, the planet-to-star radius ratio $R_{b}/R_{A}$, the limb-darkening coefficients $q_{1}$ and $q_{2}$, the mean out-of-transit flux, the dilution factors (fixed to 1), and the jitter. The input values we used were obtained by Enoch et al. (2011) for WASP-35 and RN518 for HAT-P-30 respectively, and a wide normal prior was used in the GP. For the GP component, the hyperparameters were the amplitude of the GP, two length scales corresponded to the Matern and the exponential part. We used a log-uniform prior for these hyperparameters, the amplitude of the GP varied from $10^{-6}$ to $10^{6}$, and both of the two length scales varied from $10^{-3}$ to $10^{3}$. Juliet uses PyMultiNest package (Buchner et al., 2014; Feroz et al., 2009) with 500 live points to explore the parameter space. Finally, we obtained the raw light curves and the full median posterior model, namely the transit model plus the median GP process. So we could subtract the GP part of the model from the raw light curve to correct the residual systematic errors in photometry data and obtained the final light curves. The transit light curves before and after the GP correction with the best-fitting models and the residuals are shown in Figure 1, Figure 2, Figure 3 and Figure 4.

Transit timing variations are generally parameterized with the deviations between the observed transit times and their expected values assuming a Keplerian motion for the planet. The TTVs for these two targets were computed based on our new observations, TESS photometry, transit observations collected from ETD as well as transit observations in the literature. In tables 1 and 2, there exist discrepancies between the transit times derived from simultaneous TESS and ground-based observations for some transit events, which may mean that the timing uncertainties were slightly underestimated, especially for the ground-based transit light curves which have lower quality than the TESS ones. Thus, we added a timing jitter to cope with this issue in the following TTV modeling. Refined orbital periods were obtained by linearly fitting all available transit times, in terms of a total of 55 transit times of HAT-P-30b, and 22 of WASP-35b. Figure 7 presents all the TTV measurements for transiting exoplanet WASP-35b, and Figure 8 for HAT-P-30b.

As the degeneracy in the explanation for hot Jupiter’s TTVs mentioned in Section 1, we used four different models, namely linear ephemeris model, orbital decay model, the elliptical orbit precession model, and planetary interaction-induced TTV model, to fit the timing data of HAT-P-30b and WASP-35b. In this section we describle our analyses on the first three models, while the last one will be presented in Section 5.2. Our timing analysis was partly similar to Patra et al. (2017), Yee et al. (2020) and Turner et al. (2021).

The first model assumes a constant period:

where $T_{0}$ is the reference mid-transit time and $E$ is the cycle number.

The second model assumes the planet has a constant orbital period decay rate:

where $dP/dE$ is the decay rate.

The third model assumes the planet has a nonzero orbital eccentricity $e$ and the argument of pericenter $\omega$ is precessing uniformly over time:

where $P_{s}$ is the sidereal period, $P_{a}$ is the anomalistic period and $d\omega/dE$ is the precession rate (Giménez & Bastero, 1995).

For the these models, we preformed the MCMC calculations to find the best-fitting parameters using the emcee package (Foreman-Mackey et al., 2013) and ran 50,000 MCMC steps with 1,000 burn-in steps to ensure convergence. The results of three timing model fitting are listed in Table5 and Table 6, and Figure 7 and Figure 8 show the transit timing data with the best-fitting orbital decay and apsidal precession models.

We used the Bayesian Information Criterion (BIC; Schwarz, 1978) as the penalty to compare two different best-fit models, and the BIC is defined as:

where $k$ is the number of free parameters and $n$ is the number of data points.

The results from our TTV measurements (see Section 5.1) allow to infer the upper mass limit for an additional perturbing planet in each of these two systems, which assuming that all the TTV singals origin from the gravitational interactions among planets involved. Although inverting the high-cadence and high signal-to-nosie TTV signals could be used to measure the mass and eccentricity of the perturbing exoplanets, obviously that is not our cases for HAT-P-30b and WASP-35b. However, the perturbation from the additional planet can be approximately quantified by the RMS of TTVs for our sparse measurements. The TTV effects are strongly amplified for orbital configurations in (or near) Mean-Motion Resonances (MMR) (Agol et al., 2005; Holman & Murray, 2005), in which the detection of a low-mass planetary perturbing body will in principle be permitted.

An upper mass limit can be obtained by employing a N-body code to perform direct orbit integrations, which is widely used in the literature (e.g., Wang et al., 2017, 2018a, 2018b, 2021; Cortés-Zuleta et al., 2020). Within the framework of the three-body orbital configuration, we numerically integrated the orbits of each of these two hot Jupiters and a hypothetic perturbing planet around their host star. To meet our need, we modified the TTV inversion code of Sun et al. (2019), which employs TTVFast (Deck et al., 2014) to perform direct orbit integrations. Summarily, our new code caculates the RMS of synthetic TTVs, whose epochs are selected to be identical to those of measured TTVs, considering a transiting exoplanet perturbed by an $1M_{\oplus}$ planet on an arbitrary orbital architecture. Such a configuration will help simplify the computing of upper mass limit, as the amplitude of a TTV pattern is linearly proportion to the mass of perturbing planet. Hence, the upper mass limit of perturbing planet could be well estimated by the product of the RMS of measured TTVs divided by that of synthetic TTVs.

This code was then applied to a series of orbital periods of the perturbing planet while fixing all the other orbital parameters. As the TTV patterns strongly depended on the perturber’s mass, the orbital period, the eccentricity, and the mutual inclination of the orbit (Agol et al., 2005; Holman & Murray, 2005; Nesvorný & Morbidelli, 2008; Lithwick et al., 2012; Xie et al., 2014), we performed orbital integrations in terms of three different orbital architectures for hypothetical perturbers in this study: (a) initially on coplanar and circular orbit; (b) initially on coplanar and slightly eccentric orbit (i.e., $e_{c}=0.2$); (c) initially on inclined and slightly eccentric orbit (i.e., $i_{c}=i_{b}-30\arcdeg$, $e_{c}=0.2$, where $i_{b}$ and $i_{c}$ denote the inclinations of the known transiting planet and the perturber, respectively). For clarity, hereafter we label these three different orbital architectures as Case a, Case b, and Case c, respectively. Except the orbital period of the perturbing planet, the remaining orbital elements (i.e., the longitude of the ascending node $\Omega$, the argument of periastron $\omega$, and the mean anomaly $M$ at a reference time) were fixed to selected values, $\Omega_{c}=\Omega_{b}=0$, $\omega_{c}=\omega_{b}=90\arcdeg$, $M_{c}=M_{b}+180\arcdeg$; the orbital period of the hypothetical perturber $P_{c}$ was searched from 1 day to 10$P_{b}$ with a step of 0.007 days. We expect that the first setting would provide a most conservative estimate of the upper mass limit of a possible perturber (Bean, 2009; Fukui et al., 2011; Hoyer et al., 2011, 2012; Cortés-Zuleta et al., 2020).

In addition, the RVs of the host star induced by a hypothetical perturber could place constraints on the perturber’s mass. If the planets are on non-interacting Keplerian orbits, the RVs of the host star are the sum of the RVs arosed by each planetary component’s Keplerian motion (Ford, 2006). Therefore, the residuals of both HAT-P-30 and WASP-35’s RVs after removing the contributions from both known planets could well present the RVs induced by additional bodies. Placing good constraints on the mass of the additional body through fitting RV curves requires a good coverage in its orbit phase and high signal-to-noise RV measurements, however, which is not the cases for the hypothetical perturbers in both HAT-P-30 and WASP-35. Therefore, we used the RMS of the residuals of both HAT-P-30 and WASP-35’s RVs instead of the residuals to statistically constrain the mass of the hypothetical perturber. The amplitude of RV curves $K$ arosed by a planet on its host star is scaled with the following formula:

where $M_{p}$, $M_{*}$, and $P_{orb}$ are the mass of the planet, the mass of the host star, and the orbital period of the planet, respectively. For the orbital eccentricties (i.e., $e_{c}\leq 0.2$) of the perturbers we assumed, the RMS of RV curves statistically equals about $\sqrt{2}/2$ times of the RV amplitude. We obtained the mass limits, that is, the yellow curves in Figure 10 and Figure 11 of the hypethetical perturbers based on the RMS of both HAT-P-30 and WASP-35’s RV residuals derived in Section 4.

In the following parts, we present the results of each system for which we have calculated the RMS scatter of TTVs ($\rm TTV_{\rm RMS}$) on a series of orbital periods of a perturbing planet. In each of the three cases, we found the common instability regions located in the proximity of the transiting planet with MEGNO color-coded as yellow (corresponding to $\langle Y\rangle>3.5$).

Transiting exoplanets offer the possibility of characterising the planetary atmosphere using observations at different orbital phases. These observations include the transmission spectrum during the transit (Charbonneau et al., 2002; Deming et al., 2013), the emission spectrum at secondary eclipse (Showman & Guillot, 2002; Cho et al., 2003) and the phase curve throughout the orbit(Cowan & Agol, 2008; Mallama, 2009). The transmission spectrum probes the wavelength-dependent extinction due to the planetary atmosphere at its day–night terminator region by using photometric or spectroscopic observations (Seager & Sasselov, 2000; Brown, 2001; Ehrenreich et al., 2006). For hot Jupiters, they have relatively low densities, large radius and hot equilibrium temperatures, which will lead to large atmospheric scale heights, namely, the planets have significant atmospheric structures and thus are ideal targets for atmospheric studies. Both WASP-35b and HAT-P-30b belong to hot Jupiters. Therefore, based on the new system parameters, we attempt to predict the atmospheric properties of WASP-35b and HAT-P-30b.

For the transmission spectra, the amplitude of the transmission signals in a cloud-free atmosphere is proportional to the atmospheric scale height $H$:

where $k_{B}$ is Boltzmann’s constant, $T_{eq}$ is the temperature of the planet, $\mu_{m}$ is the mean molecular mass and $g_{b}$ is the surface gravity of the planet (Seager & Sasselov, 2000; Brown, 2001; Ehrenreich et al., 2006). According to de Wit & Seager (2013), hot Jupiters should be H/He-dominated atmosphere, and the mean molecular mass is approximately 2.3 amu, so the atmospheric scale height for WASP-35b is about to 550 km and 750 km for HAT-P-30b. So, both WASP-35b and HAT-P-30b should have distinct atmospheres, the spectral features in transmission spectra could have amplitude of 5-10$H$ which are suitable for studying the properties of their atmospheres (Sing et al., 2013).

For hot Jupiters, the spectral features in the optical transmission spectrum are expected to have Na, K, TiO and VO. According to the dichotomy of highly irradiated, close-in giant planets proposed by Fortney et al. (2008) and the newly derived system parameters, both WASP-35b and HAT-P-30b are at the boundary between the pL and pM classes. Therefore we could expect that alkali metal absorption lines as well as TIO and VO likely appear in their high-precision transmission spectra. Madhusudhan et al. (2011) found the abundances of TiO and VO were strongly affected by the C/O ratio. In general, the upper atmospheric absorptions in the optical may lead to thermal inversions, such as the strong optical absorption of TiO and VO (Fortney et al., 2008) as well as the absorption of alkali metal at high C/O ratios (Mollière et al., 2015). H₂O, CO, CH₄, CO₂ and H₂ are the major chemical and spectroscopic species which dominant the C/O ratio in the hot Jupiter atmosphere and the series of reactions can be summarized as

WASP-35b and HAT-P-30b have relative high temperatures ($T_{eq}>1400$K), a simple way to estimate the C/O ratio of the atmosphere is observing the transmission spectra of H₂O and CH₄ in the near- to mid-infrared. Both of the systems have zero eccentricities, this means the planet remains at a constant distance from its host star, which avoids orbital-induced thermal response of the planetary atmospheric temperature and hence chemistry (Visscher, 2012). Clouds and haze are common in the planetary atmospheres and strongly impact on observations of transmission spectra (Sing et al., 2013; Spyratos et al., 2021). Uniform global clouds can obscure the absorption features of prominent chemical species, patchy clouds can mimic high mean molecular weight atmospheres (Fortney, 2005; Line & Parmentier, 2016) and haze can produce a significant slope at optical passband in transmission spectra (Lecavelier Des Etangs et al., 2008). According to the theory of Lecavelier Des Etangs et al. (2008), assuming a power law for the cross section with wavelength in the form of $\sigma=\sigma_{0}(\lambda/\lambda_{0})^{\alpha}$, the slope induced by Rayleigh scattering is given by

For the cloud-free atmosphere, the Rayleigh scattering is caused by H₂ and He, the expected value of $\alpha$ is -4. For the cloud atmosphere, the observable atmosphere would be smaller than pure gaseous atmosphere, so that the amplitude of absorption features would be small or the transmission spectra would be flatten and featureless.