Transiting Exoplanet Monitoring Project (TEMP). II. Refined System Parameters and Transit Timing Analysis of HAT-P-33b

As the length of the catalog mounts,¹⁵ so to does the importance of characterizing the alien worlds. A better understanding of the extrasolar planets' compositions, their formation and their evolution constitutes a grand challenge for the twenty-first century. With these larger goals as motivation, we initialized the Transit Exoplanets Monitoring Project (TEMP) to specifically study the transiting exoplanet systems with high-precision photometric follow-up observations (Wang et al. 2016).

High-precision photometric follow-ups lead to more accurate measurements of planetary radii and orbital inclinations, and combined with the RV method permits determinations of planetary masses, which in turn give densities, and hence the planetary compositions (Sato et al. 2005). With improved photometry, we can determine more precise orbital ephemerides, which streamline future research studies, including those that draw on the Rossiter–McLaughlin effect (Nutzman et al. 2011; Sanchis-Ojeda & Winn 2011; Sanchis-Ojeda et al. 2013), transmission spectra (Mancini et al. 2016) and spectroscopy at secondary eclipse (Star Cartier et al. 2016). Furthermore, we can perform transit timing variation (TTV) analysis with high-precision photometric data. These provide us the powerful tools to detect close-in companions in known hot-Jupiter systems and hence enable the zeroth-order test of competing formation scenarios for hot Jupiters (Lin et al. 1996; Bodenheimer et al. 2000; Wu & Murray 2003; Ford & Rasio 2008; Nagasawa et al. 2008; Wu & Lithwick 2011; Batygin et al. 2016). Moreover, with TTVs in hand, we can confirm the planetary nature and measure masses for planets in multi-transiting systems (Lithwick et al. 2012; Xie et al. 2014; Hadden & Lithwick 2017). Most of the multi-transiting systems that are currently known were detected by the Kepler space telescope (Fabrycky et al. 2014). The host stars, however, of the Kepler-detected systems are mostly too faint for feasible RV follow-up observations from the ground using small to medium aperture optical telescopes, incentivizing the search for TTVs among these systems.

In addition, we can confirm candidate exoplanets and refine their orbital ephemerides through photometric monitoring of planets that have been observed for a limited number of transits in the K2 data sets (Howell et al. 2014) and/or from the forthcoming TESS mission (Ricker et al. 2015). High-precision photometric follow-ups in multi-band also provide a means to determine the chemical compositions and atmospheric properties for exoplanets (Fukui et al. 2013; Lendl et al. 2013; Mancini et al. 2013; Sing et al. 2016).

In the first stage of TEMP, we have focused primarily on monitoring the hot Jupiters found by ground-based transiting surveys, concentrating on those for which only limited photometric follow-up observations have been published. In many cases, the parameters for such systems are both imprecise and incomplete. These targets, therefore, offer an optimal scientific benefit and studies similar to these provided by TEMP have been presented by other groups (Becker et al. 2015; Seeliger et al. 2015; Collins et al. 2017), demonstrating that the TTVs provide a key avenue for insights into exoplanetary formation and evolution. In this paper, we present the scientific results that emerged from a monitoring campaign on HAT-P-33b. This system was chosen because of its large RV residuals and the existence of only sparse data in the form of incomplete light curves (Hartman et al. 2011).

HAT-P-33b was discovered by Hartman et al. (2011), who found the planet to be a highly inflated hot Jupiter (${M}_{{\rm{P}}}=0.763\,\,{M}_{{\rm{J}}},{R}_{{\rm{P}}}=1.827\,\,{R}_{{\rm{J}}}$) transiting a late-F dwarf star (${M}_{* }=1.403\,{M}_{\odot },{R}_{* }=1.777\,{R}_{\odot }$) with an orbital period of $3.474474\,\mathrm{days}$. Seven light curves were presented in their work, but only one is complete. Following the discovery by Hartman et al. (2011), four more RV measurements (obtained also using Keck HIRES) were presented in Knutson et al. (2014) bringing the total number of RV observations to be 26. The extended data showed no evidence of long-period companions within the HAT-P-33b system.

In this work, we present ten new light curves of eight different transits of HAT-P-33b. The light curves are all complete, with a typical photometric precision better than $2.0\,\mathrm{mmag}$ save one, which is partial and has a precision of $3.0\,\mathrm{mmag}$. Based on our photometric data and the extended RV measurements (Knutson et al. 2014), we revisit the system parameters, refine the orbital ephemeris, and explore the possibility of existence of additional planets in the system.

This paper is organized as follows. We describe the photometric observations and data reduction in Section 2. The data analysis is presented in Section 3. In Section 4, we give the results and discussion. Finally, a brief summary of our work is presented in Section 5.

We have recorded a total of 10 light curves of eight different transits events observed by two telescopes (a $60/90\,\mathrm{cm}$ Schmidt and a $60\,\mathrm{cm}$ telescope) at Xinglong Station operated by National Astronomical Observatories of China (NAOC) between 2013 December and 2016 January. Two of the transit events were observed by the two telescopes simultaneously.

The first seven transit events were monitored by the $60/90\,\mathrm{cm}$ Schmidt telescope. It has a $4{\rm{K}}\times 4{\rm{K}}$ CCD with a $\sim 94^{\prime} \times \sim 94^{\prime}$ field of view, which gives a pixel scale of $1\buildrel{\prime\prime}\over{.} 38\,\mathrm{pixel}{}^{-1}$ and a typical readout time of $93\,{\rm{s}}$ (Zhou et al. 1999, 2001). To reduce the readout times and increase the duty cycle of the observations, the images were windowed down to 512 × 512 pixels with 1 × 1 binning, resulting in a reduced readout time of $4\,{\rm{s}}$.

The transit events that occurred on UT 2014 February 27 and UT 2014 March 6 were also simultaneously observed by the $60\,\mathrm{cm}$, and the last transit in our sequence was also monitored with this telescope. The $60\,\mathrm{cm}$ telescope is equipped with a 512 × 512 CCD and covers a field of view of $17^{\prime} \times 17^{\prime}$, resulting in a pixel scale of $1\buildrel{\prime\prime}\over{.} 95\,\mathrm{pixel}{}^{-1}$. No windowing with 1 × 1 binning was performed during these observations, giving a standard readout time of $3\,{\rm{s}}$.

It is common practice to defocus the telescope in order to optimize signal-to-noise and to keep photoelectron counts within the CCD's range of linear response (Southworth et al. 2009). Broadened stellar point-spread functions are less sensitive to focus or telescope pointing changes, which would otherwise cause systematic errors. Defocusing produces longer exposure times, which increase the duty cycle of observations and reduce Poisson or scintillation noise (Hinse et al. 2015).

In our observations of HAT-P-33b, which is a ${V}_{\mathrm{mag}}=11.19$ star, we slightly defocused our telescopes. The linear range of the CCD is maintained for target counts less than 30,000. For the sake of conservatism, we keep our target at a typical count of 20,000 which is reached within 15 s in a clear night. Fifteen seconds is too short, however, to achieve optimal reduction of Poisson and scintillation noise, which further motivates our defocusing of the CCD images. The background count is about 300 within our typical 60 s exposure time. We adjust exposure times throughout each data-taking session in order to maintain counts that fall within the linear regime of the CCD. The exposure time, however, was kept fixed during the ingress and egress phases to avoid affecting the precision of transit timing, which is a critical aspect of our work. The telescope time stamp server was synchronized on a nightly basis with the US Naval Observatory (USNO) time.¹⁶ Timings are measured accurately to within one second and recorded using the UTC time standard. A summary of our observations is listed in Table 1.

Table 1.  Overview of Observations and Data Reduction

Date	Time	Telescopea	Filter	Frames	Exposure	Read	Airmass	Moon	Comp.	Apertureb	Scatterc
(UTC)	(UTC)				(s)	(s)		illum.	Stars	(pixels)	(mmag)
2013 Dec 09	14:40:56-20:05:27	Schmidt	R	247	60, 80	4	$1.41\to 1.01\to 1.07$	0.51	5	10	2.0
2014 Feb 27	10:54:42-18:21:57	$60\,\mathrm{cm}$	R	373	60, 70	3	$1.13\to 1.01\to 2.01$	0.04	3	11	2.0
2014 Feb 27	11:01:37-17:54:05	Schmidt	R	295	60, 110	4	$1.12\to 1.01\to 1.75$	0.04	3	14	1.9
2014 Mar 06	10:48:53-18:02:13	$60\,\mathrm{cm}$	R	397	60	3	$1.09\to 1.01\to 2.10$	0.32	2	10	1.6
2014 Mar 06	11:09:38-17:27:35	Schmidt	R	251	80, 110	4	$1.06\to 1.01\to 1.76$	0.32	2	14	1.7
2015 Jan 16	14:52:35-20:47:16	Schmidt	R	343	55	4	$1.03\to 1.01\to 1.81$	0.19	4	14	1.7
2015 Jan 23	13:23:02-19:45:43	Schmidt	R	322	52	4	$1.12\to 1.01\to 1.56$	0.14	4	11	1.9
2015 Jan 30	12:23:56-18:56:37	Schmidt	R	465	40, 50	4	$1.18\to 1.01\to 1.44$	0.84	4	10	1.9
2015 Feb 13	10:40:54-16:35:31	Schmidt	R	353	45, 55	4	$1.34\to 1.01\to 1.15$	0.34	2	13	1.8
2016 Jan 09	12:16:06-15:15:30	$60\,\mathrm{cm}$	R	695	10	3	$1.55\to 1.04$	0.00	3	10	3.0

Notes.

^aThe telescopes belong to the Xinglong Station operated by National Astronomical Observatories of China (NAOC). ^bThe aperture indicates the aperture diameter around stars. ^cThe scatter indicates the rms of residuals from our best-fitting model.

Download table as:  ASCIITypeset image

All the data have been calibrated using a standard procedure, including overscan correction and flat-fielding for data from the Schmidt telescope, as well as bias correction and flat-fielding for data from the $60\,\mathrm{cm}$ telescope. Twilight sky flats were obtained by the $60\,\mathrm{cm}$ telescope, whereas dome flats were taken with the Schmidt. We used SExtractor (Bertin & Arnouts 1996) to perform aperture differential photometry. All the stars in the field with enough flux were tested for photometric non-variability, and the most favorable sources were chosen as reference stars. With the reference stars, we obtained the differential light curve, which has the smallest root-mean-square (rms) scatter for each transit, by manually varying the aperture diameter from 8 to 16 pixels. A summary of the aperture photometry is given in Table 1. We then removed trends that may be caused by the variation of airmass and intrinsic stellar variability, by performing a linear fit to the out-of-transit data. To maintain timing consistency, we converted the UTC time stamps to Barycentric Julian Date in the TDB time standard (${\mathrm{BJD}}_{\mathrm{TDB}}$) for each light curve using the online procedure.¹⁷ The final set of 10 recorded light curves are listed in Table 2; in total, these data comprise 3732 measurements.

Table 2.  Photometry of HAT-P-33

${\mathrm{BJD}}_{\mathrm{TDB}}$ a	Relative Flux	Scatter	Telescope	Filter
2456636.117240	0.9990	0.0020	Schmidt	R
2456636.117981	0.9990	0.0020	Schmidt	R
2456636.118722	1.0000	0.0020	Schmidt	R
2456636.121361	1.0000	0.0020	Schmidt	R
2456636.122877	0.9963	0.0020	Schmidt	R
2456636.123618	1.0018	0.0020	Schmidt	R
2456636.124359	1.0009	0.0020	Schmidt	R
2456636.125099	1.0018	0.0020	Schmidt	R
2456636.125840	1.0036	0.0020	Schmidt	R
2456636.126569	1.0000	0.0020	Schmidt	R
2456636.127310	0.9972	0.0020	Schmidt	R
2456636.128051	0.9990	0.0020	Schmidt	R
2456636.128803	1.0000	0.0020	Schmidt	R
2456636.129521	0.9972	0.0020	Schmidt	R
2456636.130273	0.9990	0.0020	Schmidt	R

Note.

^aAll the timing throughout the paper is based on ${\mathrm{BJD}}_{\mathrm{TDB}}$, calculated from Coordinated Universal Time (UTC) using the procedure developed by Eastman et al. (2010).

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

We applied the EXOFAST¹⁸ (a fast exoplanetary fitting package in IDL) developed by Eastman et al. (2013) to perform our data modeling. The package can simultaneously fit transit and RV data with given priors, robustly deriving the parameter values and their uncertainties using the differential evolution Markov chain Monte Carlo (DE-MC) algorithm. At each Markov chain step, EXOFAST employs the Torres relations to calculate ${M}_{* }$ and ${R}_{* }$ with given ${T}_{\mathrm{eff}},[\mathrm{Fe}/{\rm{H}}]$, and $\mathrm{log}({g}_{* })$ (Torres et al. 2010).

To revisit the system parameters of HAT-P-33b, we performed a global fit based on our seven light curves and the extended RV measurements from Knutson et al. (2014). The priors of the system parameters used in the fit were obtained from Hartman et al. (2011) and are presented in Table 3. We also obtained the priors for the limb-darkening parameters in the R band (u₁ = 0.2631, u₂ = 0.3155), following the description in Claret & Bloemen (2011). As a first step, EXOFAST fitted the RV and transit data sets independently and scaled the uncertainties to obtain a reduced ${\chi }_{\mathrm{red}}^{2}=1$ for each best-fitting model. Then it performed a global fit based on both data sets. A total of 32 simultaneous chains were constructed in our fit, each having a maximum of 100,000 steps. As described in Eastman et al. (2013), the Markov chains are considered to have converged when both the Gelman–Rubin statistic is less than 1.01 and the number of independent draws is greater than 1000 for all parameters. Only after passing this test six consecutive times can the chains be considered well-mixed and EXOFAST will stop. As a final step, we evaluated the well-mixed results to obtain best-fitting values with $1\sigma$ error bars for system parameters, which are also listed in Table 3.

Table 3.  System Parameters for HAT-P-33

Parameter	Units	This Work	Hartman et al. (2011)	Knutson et al. (2014)
Stellar Parameters:
M*	Mass (M⊙)	${1.42}_{-0.15}^{+0.16}$	1.403 ± 0.096	1.403 ± 0.096a
${R}_{* }$	Radius (R⊙)	${1.91}_{-0.20}^{+0.26}$	1.777 ± 0.280	⋯
${L}_{* }$	Luminosity (L⊙)	${5.7}_{-1.6}^{+2.3}$	${4.73}_{-1.25}^{+1.87}$	⋯
${\rho }_{* }$	Density (cgs)	${0.289}_{-0.081}^{+0.098}$	⋯	⋯
$\mathrm{log}({g}_{* })$	Surface gravity (cgs)	${4.030}_{-0.090}^{+0.079}$	4.09 ± 0.11	4.09 ± 0.11a
${T}_{\mathrm{eff}}$	Effective temperature (K)	${6460}_{-290}^{+300}$	6401 ± 88	6401 ± 88a
$[\mathrm{Fe}/{\rm{H}}]$	Metalicity	0.01 ± 0.31	0.05 ± 0.08	0.05 ± 0.08a
Planetary Parameters:
e	Eccentricity	${0.180}_{-0.096}^{+0.11}$	0.148 ± 0.081	${0.13}_{-0.1}^{+0.19}$
${\omega }_{* }$	Argument of periastron (degrees)	${88}_{-34}^{+33}$	96 ± 119	15 ± 22
P	Period (days)	3.47447472 ± 0.00000088b	3.474474 ± 0.000001	⋯
a	Semimajor axis (au)	0.0505 ± 0.0018	0.0503 ± 0.0011	⋯
MP	Mass (MJ)	${0.72}_{-0.12}^{+0.13}$	0.763 ± 0.117	0.65 ± 0.14
RP	Radius (RJ)	${1.87}_{-0.20}^{+0.26}$	1.827 ± 0.290	⋯
${\rho }_{{\rm{P}}}$	Density (cgs)	${0.134}_{-0.042}^{+0.053}$	${0.15}_{-0.05}^{+0.11}$	⋯
$\mathrm{log}({g}_{{\rm{P}}})$	Surface gravity	${2.70}_{-0.11}^{+0.10}$	2.75 ± 0.13	⋯
Teq	Equilibrium Temperature (K)	${1920}_{-120}^{+140}$	1838 ± 133	⋯
Θ	Safronov Number	${0.0271}_{-0.0050}^{+0.0056}$	${0.030}_{-0.007}^{+0.005}$	⋯
$\langle F\rangle$	Incident flux (109 erg s−1 cm−2)	${2.96}_{-0.65}^{+0.84}$	${2.58}_{-0.61}^{+0.93}$	⋯
RV Parameters:
$e\cos {\omega }_{* }$	⋯	${0.004}_{-0.086}^{+0.092}$	0.040 ± 0.078	${0.114}_{-0.097}^{+0.16}$
$e\sin {\omega }_{* }$	⋯	${0.154}_{-0.096}^{+0.11}$	0.073 ± 0.138	${0.015}_{-0.023}^{+0.13}$
TP	Time of periastron (${\mathrm{BJD}}_{\mathrm{TDB}}$)	${2457046.20}_{-0.23}^{+0.22}$	⋯	⋯
K	RV semi-amplitude (m s−1)	78 ± 12	82.8 ± 12.0	${72}_{-16}^{+19}$
${M}_{P}\sin i$	Minimum mass (MJ)	${0.72}_{-0.12}^{+0.13}$	⋯	⋯
${M}_{{\rm{P}}}/{M}_{* }$	Mass ratio	${0.000484}_{-0.000076}^{+0.000077}$	⋯	⋯
γ	Systemic velocity (m s−1)	−7 ± 11	⋯	⋯
$\dot{\gamma }$	RV slope (m s−1 day−1)	−0.024 ± 0.018	⋯	$-{0.021}_{-0.023}^{+0.02}$
Primary Transit Parameters:
TC	Time of transit (${\mathrm{BJD}}_{\mathrm{TDB}}$)	2456035.137750 ± 0.000272b	2455100.50255 ± 0.00023	⋯
${R}_{{\rm{P}}}/{R}_{* }$	Radius of planet in stellar radii	${0.10097}_{-0.00052}^{+0.00056}$	0.1057 ± 0.0011	⋯
$a/{R}_{* }$	Semimajor axis in stellar radii	${5.69}_{-0.59}^{+0.58}$	${6.08}_{-0.72}^{+0.98}$	⋯
u1	linear limb-darkening coeff	0.264 ± 0.026	⋯	⋯
u2	quadratic limb-darkening coeff	0.315 ± 0.037	⋯	⋯
i	Inclination (degrees)	${88.2}_{-1.3}^{+1.2}$	${86.7}_{-1.2}^{+0.8}$	${86.7}_{-1.2}^{+0.8}$ a
b	Impact Parameter	${0.151}_{-0.098}^{+0.10}$	0.325 ± 0.002	⋯
δ	Transit depth	0.01020 ± 0.00011	⋯	⋯
TFWHM	FWHM duration (days)	${0.16354}_{-0.00072}^{+0.00070}$	⋯	⋯
τ	Ingress/egress duration (days)	${0.01707}_{-0.00036}^{+0.00080}$	0.0194 ± 0.0002	⋯
T14	Total duration (days)	${0.18075}_{-0.00089}^{+0.00097}$	0.1836 ± 0.0007	⋯
PT	A priori non-grazing transit prob	${0.188}_{-0.035}^{+0.054}$	⋯	⋯
${P}_{T,G}$	A priori transit prob	${0.231}_{-0.043}^{+0.066}$	⋯	⋯
F0	Baseline flux	0.999837 ± 0.000061	⋯	⋯
Secondary Eclipse Parameters:
TS	Time of eclipse (${\mathrm{BJD}}_{\mathrm{TDB}}$)	${2457044.48}_{-0.19}^{+0.21}$	2455102.330 ± 0.175	⋯
bS	Impact parameter	${0.21}_{-0.13}^{+0.14}$	⋯	⋯
${T}_{S,\mathrm{FWHM}}$	FWHM duration (days)	${0.219}_{-0.037}^{+0.053}$	⋯	⋯
${\tau }_{S}$	Ingress/egress duration (days)	${0.0239}_{-0.0044}^{+0.0064}$	0.0230 ± 0.0085	⋯
${T}_{S,14}$	Total duration (da2ys)	${0.243}_{-0.041}^{+0.059}$	0.2090 ± 0.0480	⋯
PS	A priori non-grazing eclipse prob	${0.1379}_{-0.0031}^{+0.0046}$	⋯	⋯
${P}_{S,G}$	A priori eclipse prob	${0.1689}_{-0.0038}^{+0.0057}$	⋯	⋯

Notes. The published system parameters of HAT-P-33 from the literature (Hartman et al. 2011; Knutson et al. 2014) are presented for comparison.

^aIn Knutson et al. (2014), the stellar parameters (including ${M}_{* },\mathrm{log}({g}_{* }),{T}_{\mathrm{eff}},[\mathrm{Fe}/{\rm{H}}]$) and orbital inclination (i) were adopted from Hartman et al. (2011). ^bWe got P and T_C through a linear fit based on the mid-transit times, which are calculated from the new and published light curves (see Section 4.2).

Download table as:  ASCIITypeset image

To accurately measure the mid-transit times for all 17 light curves, we separately performed a fit for each light curve in conjunction with the extended RV measurements from Knutson et al. (2014). The time stamps of the published light curves were converted to ${\mathrm{BJD}}_{\mathrm{TDB}}$ for reasons of consistency. In these fits, we fixed the system parameters to the values obtained from the aforementioned global fit, excepting ${T}_{{\rm{c}}}$ and baseline flux of the light curve (F₀), which were allowed to float instead. For the published light curves from Hartman et al. (2011), the limb-darkening parameters were fixed to different values in diverse bands¹⁹ during fitting processes. After a fitting process, similar to the global fit mentioned above, we had an estimate for the mid-transit time for each transit event.

4.1. System Parameters

As the result of global fit, the final parameters for the HAT-P-33 system together with the results from previous work (Hartman et al. 2011; Knutson et al. 2014) are listed in Table 3. The resulting best-fitting models for the combined photometric and RV data are plotted in Figures 1 and 2, separately.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As expected, our RV parameters are consistent with those of Knutson et al. (2014), which resulted from the same RV data sets. These RV parameters also agree with the results from Hartman et al. (2011), though their RV data set contains four fewer points. As with Knutson et al. (2014), we did not find a long-period trend in the RV residuals, so we give the minimum mass of a potential planetary perturber following the convention defined by Wright et al. (2007).

The resulting transit parameters also agree with those from Hartman et al. (2011) except some with slight differences, including a lower impact parameter (b), a smaller value for the planet to star radius ratio (${R}_{{\rm{P}}}/{R}_{* }$), a shorter ingress/egress duration (τ), a shorter total duration (T₁₄), and a larger inclination (i). Comparing to the published work, which was based on only one full-transit light curve, our results are more robust, as a consequence of being based on the seven complete light curves.

Our stellar parameters show agreement with those of Hartman et al. (2011), which were chosen as the spectroscopic priors for the global fit in advance. Finally, the planetary parameters of HAT-P-33b calculated based on the derived RV, transit, and stellar parameters also agree well with those in Hartman et al. (2011).

4.2. Mid-transit Times

In order to revisit the orbital ephemeris and seek TTV signals for the HAT-P-33b system, we acquired accurate mid-transit times (${T}_{{\rm{c}}}$) through separately fitting each light curve. The best-fitting models are shown in Figure 3 and the resulting mid-transit times are listed in Table 4, with uncertainties obtained with the DE-MC method. We fitted obtained transit times with a linear function of transit epoch number (N),

$$\begin{eqnarray}&&{T}_{{\rm{c}}}[N]={T}_{{\rm{c}}}[0]+{NP},\end{eqnarray} \tag{ 1 }$$

where P is the planetary orbital period, ${T}_{{\rm{c}}}[0]$ represents the zero epoch. The best-fitting values are

$$\begin{eqnarray}&&{T}_{{\rm{c}}}[0]=2456035.137750\pm 0.000272\,[{\mathrm{BJD}}_{\mathrm{TDB}}],\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&P=3.47447472\pm 0.00000088[\mathrm{days}].\end{eqnarray} \tag{ 3 }$$

Our orbital ephemeris agree well with the result from Hartman et al. (2011).

Zoom In Zoom Out Reset image size

Table 4.  Mid-transit Times for HAT-P-33b

Epocha	Telescopeb	${T}_{{\rm{c}}}$	${\sigma }_{{T}_{{\rm{c}}}}$	O − C
		(${\mathrm{BJD}}_{\mathrm{TDB}}$)	(s)	(s)
−557	FLWO	2454099.85310	41.47	−192.40
−462	FLWO	2454429.93117	47.52	64.29
−427	FLWO	2454551.53949	50.11	211.56
−357	FLWO	2454794.75180	31.70	112.99
−355	FLWO	2454801.69988	63.94	56.88
−123	FLWO	2455607.77606	51.84	−112.18
−121	FLWO	2455614.72513	31.24	−101.77
173	Schmidt	2456636.22181	62.21	−5.85
196	$60\,\mathrm{cm}$	2456716.13465	59.62	−12.65
196	Schmidt	2456716.13498	56.16	15.86
198	$60\,\mathrm{cm}$	2456723.08352	40.15	−19.52
198	Schmidt	2456723.08417	50.52	36.64
289	Schmidt	2457039.26162	39.74	58.23
291	Schmidt	2457046.20981	39.74	−7.39
293	Schmidt	2457053.15793	40.61	−79.05
297	Schmidt	2457067.05749	37.15	64.46
392	$60\,\mathrm{cm}$	2457397.13183	57.02	−1.12

Notes.

^aThe first seven time points are obtained from the published light curves (Hartman et al. 2011) through separate fits, the others are from our photometric data. As mentioned above, the epochs (239, 241) were followed by two telescopes simultaneously. ^bFor more information about the FLWO telescope, see Hartman et al. (2011).

Download table as:  ASCIITypeset image

To get conservative uncertainty estimates for a more reliable future observation schedule, the uncertainties for the mid-transit times during the fitting were rescaled through a common factor to get ${\chi }^{2}/{N}_{\mathrm{dof}}=1$. However, the uncertainties of mid-transit times listed in Table 4 were not rescaled in this way, nor were the error bars plotted in Figure 4.

Zoom In Zoom Out Reset image size

Figure 4 displays the deviations of mid-transit times from the linear orbital ephemeris (Equations (1)–(3)), with an rms of $93.68\,{\rm{s}}$. This value is largely affected by the mid-transit times derived from the published light curves, which gives an rms of $144.12\,{\rm{s}}$ over a $4\,\mathrm{year}$ time span. As a contrast, the rms of mid-transit times derived from our data is only $41.87\,{\rm{s}}$ within a time span of $3\,\mathrm{year}$. The large deviation of mid-transit times from the published light curves may be mainly caused by their incomplete coverages. In total, most of the mid-transit times are in the $\pm 3\sigma$ errors range of the orbital ephemeris. Especially for the mid-transit times resulted from our photometric data, which are very consistent with the $\pm 1\sigma$ errors.

4.3. Limits On Additional Perturbers

Although neither significant TTVs nor a residual RV signal were found, we can place the upper mass limits of a potential close-in perturbing planet in the HAT-P-33 system. The results are shown in Figure 5.

Zoom In Zoom Out Reset image size

The host star (HAT-P-33) is an active late-F dwarf (Hartman et al. 2011), which has a large RV uncertainty $(\mathrm{rms}\,=47.5\,{\rm{m}}\,{{\rm{s}}}^{-1})$. The mass limits based on the RV residuals following the convention in Wright et al. (2007) is thus very loose, as indicated by the black dashed line in Figure 5, which can only exclude a perturber with a mass larger than $0.6\,\,{M}_{{\rm{J}}}$ near the 1:5 resonances (0.69 day orbit) or $2.0\,\,{M}_{{\rm{J}}}$ near the 5:1 resonances (17.37 day orbit).

Fortunately, TTV measurements are less sensitive to the stellar activity than are Doppler measurements. We made use of the MERCURY6 orbit integration package Chambers (1999) to place upper mass limits of a potential perturbing body. The TTV data exhibited an rms scatter of $93.68\,{\rm{s}}$.

In our simulations, we assumed that the orbits for both the known hot Jupiter and a potential perturber are coplanar and circular, which gives the most conservative estimate of upper mass limits of the potential perturber (Bean 2009; Fukui et al. 2011). The arguments of periastron, ω, the ascending nodes, Ω, and the initial mean anomalies, M₀, of the known hot Jupiter and a potential perturber are fixed to $\omega =88^\circ$ (from Table 3), ${\rm{\Omega }}=270^\circ$, and ${M}_{0}=0^\circ$. We explored the mass space of the potential perturber for both interior and exterior orbits with the orbital period ratio from 1/5 to 5 times (0.69–$17.37\,\mathrm{days}$) that of HAT-P-33b, which is equivalent to a semimajor axis range from 0.017 to $0.154\,\mathrm{au}$. We incremented the perturber's semimajor axis by $0.00001\,\mathrm{au}$. The resolution is enough to depict the constraints on the perturber mass in the resonant configurations, that the TTV signals are significantly sensitive to (Holman & Murray 2005; Agol & Steffen 2007). In each increment of a, we obtained the upper mass limit of the potential perturber by iterative linear interpolation with an initial mass of $1.0\,{M}_{\oplus }$ and a convergence tolerance of $1.0\,{\rm{s}}$ for the TTVs.

Comparing to the loose limits by RV data, the mass limits from our TTV measurements are much tighter near the low-order mean-motion resonances, as illustrated by the black solid line in Figure 5. We can exclude the existence of a perturber with mass larger than 0.6, 0.3, 0.5, 0.5, and $0.3\,{M}_{\oplus }$ near the 1:3, 1:2, 2:3, 3:2, and 2:1 resonances, respectively.

In Figure 5, we also present the the dynamical stability in the hypothetical three-body system through the Mean Exponential Growth of Nearby Orbits (MEGNO) Index (Goździewski et al. 2001; Cincotta et al. 2003; Hinse et al. 2010). The resulting dynamical stability map agrees well with that obtained by the analytic method described in Barnes & Greenberg (2006).

We initiated the TEMP to study the known exoplanets in great detail, with specific goals of obtaining a better grasp of planetary interior structures, formation, and evolution.

One of the initial targets for TEMP, HAT-P-33b, has been observed by two telescopes from 2013 December to 2016 January. In total, we obtained 10 light curves of eight different transit events, thereby substantially enriching the photometric database of HAT-P-33b.

To revisit the system parameters of HAT-P-33b, we have performed a global fit based on our new light curves and the expanded RV data (Knutson et al. 2014). Though most of the results agree well with those from the published work (Hartman et al. 2011; Knutson et al. 2014), some slight discrepancies still exist in the transit parameters.

We also separately conducted fits for the 17 light curves to obtain precise mid-transit times. Along with these, we revisited the orbital ephemeris for HAT-P-33b, which agrees well with that in Hartman et al. (2011).

Though no substantial TTV signal has been found from the linear orbital ephemeris of HAT-P-33b, we can constrain the upper mass limits of a potential close-in perturbing planet based on the measured TTVs with an rms scatter of $93.68\,{\rm{s}}$. The restriction is much stronger near the low-order mean-motion resonances. We can exclude the existence of a planet with mass larger than 0.6, 0.3, 0.5, 0.5, and $0.3\,{M}_{\oplus }$ near the 1:3, 1:2, 2:3, 3:2, and 2:1 resonances, respectively. However, we still cannot rule out the existence of additional close-in planets in the non-resonant area. Whether additional planets frequently exist in the nearby non-resonant area of hot Jupiters is an open question. Further work is needed to answer this question, and hence to better reveal the nature of planetary formation and evolution.

This research is supported by the Strategic Priority Research Program: The Emergence of Cosmological Structures of the Chinese Academy of Sciences (grant No. XDB09000000); the National Basic Research Program of China (Nos. 2013CB834900, 2014CB845704, 2013CB834902, and 2014CB845702); the National Natural Science Foundation of China (under grant Nos. 11333002, 11433005, 11373033, 11503009, 11003010, 11373035, 11203034, 11203031, 11303038, 11303043, 11073032, 11003021, 11603035, 11673027, 11603034, and 11173016); the Main Direction Program of Knowledge Innovation of Chinese Academy of Sciences (No. KJCX2-EW-T06); Japan Society for Promotion of Science (JSPS) KAKENHI Grant Numbers JP25247026. This work is also supported by the External Cooperation Program of Chinese Academy of Sciences (grant No. 114A11KYSB20160057). Songhu Wang gratefully acknowledges the award of a Heising-Simons 51 Pegasi Postdoctoral Fellowship. Tobias C. Hinse acknowledges KASI research grant 2016-1-832-01. Numerical computations were partly carried out using the SFI/HEA Irish Center for High-End Computing (ICHEC) and the 3rd generation Polaris High-Performance Computing cluster at KASI/South Korea. Research at the Armagh Observatory is funded by the Department of Culture, Arts and Leisure (DCAL).

Facilities: Beijing:Schmidt - , Beijing:0.6m - .

Software: SExtractor (Bertin & Arnouts 1996), EXOFAST (Eastman et al. 2013).