APOSTLE: LONGTERM TRANSIT MONITORING AND STABILITY ANALYSIS OF XO-2b

Observational efforts on two fronts (ground-based and space-based) have revealed a great diversity in the properties of planets (and planet candidates), and posed several new questions about the planet formation process. The search for planets and the careful measurement of planetary properties are the observational foundations upon which the physics of planet formation may be understood. The exoplanet community has co-opted the term "architecture" to embody several properties of planetary systems such as multiplicity and orbital parameters like eccentricities (e), inclinations (i_orb), semi-major axes (a), and orbital periods. The distributions of individual planetary properties like masses and radii (M_p and R_p, respectively) may also be included in this term. The primary goal of the Apache Point Observatory Survey of Transit Lightcurves of Exoplanets (APOSTLE) is to catalog transit light curves at high precision in order to (1) measure system parameters and (2) look for transit timing variations (TTVs) that may indicate the presence of additional planets.

The transit technique applies to those systems where the orbital inclination of an exoplanet is close to 90° (i.e., edge-on) with respect to the observer's sky-plane (see the discovery paper Charbonneau et al. 2000). In this case, the observer sees a u-shaped dip in the starlight caused when the planet eclipses the star (Winn 2011). The objective of several ground-based and space-based efforts focused on transit observations is to catalog and improve measurements of system parameters, which in turn gives us an improved picture of the architecture of the planetary system. This process is key toward developing theories of planet formation that can adequately explain the origin and evolution of all planetary systems (including our own).

The APOSTLE target XO-2b is a Hot-Jupiter on a 2.6 day orbit around an early type K dwarf (V = 11.2; Burke et al. 2007). The planet is known to have a mass of 0.555 M_Jup and a radius of 0.992 R_Jup (Southworth 2010). The system is not known to have any other planet-mass objects. However, Narita et al. (2011) claim to detect long-term radial velocity variation that deserves further follow-up. Transit timing measurements from the ground are consistent with a linear ephemeris (Fernandez et al. 2009; Sing et al. 2011; Crouzet et al. 2012), although the cited studies note that there are statistically significant deviations in the measurement of the orbital periods. A search for additional eclipses using the EPOXI mission by Ballard et al. (2011) had inconclusive results since the eclipse was not fully sampled. Spitzer observations of the secondary eclipse of XO-2b show IRAC 3.6, 4.5, and 5.8 μm fluxes which are consistent with the presence of a temperature inversion (Machalek et al. 2009). There have also been detections of optical absorbers in the planetary atmosphere such as potassium from narrow band optical transmission spectrophotometry (Sing et al. 2011). Early theoretical studies indicated that stellar insolation levels directly influenced the presence or absence of a thermal inversion layer (depending on the survival of atmospheric absorbers; Hubeny et al. 2003; Burrows et al. 2007; Fortney et al. 2008). The planetary atmosphere classification system developed by Fortney et al. (2008) places XO-2b in an transition zone between planets with (pM) and without (pL) thermal inversions. XO-2b is one of a handful of Hot-Jupiters in this region. However, observational evidence suggests (Machalek et al. 2008) that more complicated models need to be considered. Irradiation from stellar activity may need to be included (Knutson et al. 2010) and atmospheric chemistry may also need to be considered to provide a more complete picture (Madhusudhan 2012). The host star of the XO-2 system is known to be a high metallicity ([Fe/H] = 0.45 ± 0.02), high proper motion star (μ_Tot = 157 mas yr⁻¹) in a visual binary (Burke et al. 2007). Spectral activity indices show that XO-2 is fairly inactive compared to other stars of similar spectral type (Knutson et al. 2010).

A primary goal of APOSTLE and other campaigns that monitor transiting exoplanets is to search for TTVs that reveal the presence of unseen companions (Agol et al. 2005; Holman & Murray 2005). In principle, Earth-mass planets in or near mean motion resonances could perturb the orbit of the transiting planet enough to produce a sinusoidal oscillation in the mid-points of transits. However, the full range of stable orbits that can produce a detectable TTV signal has never been explored. The detection of TTVs by Kepler (e.g., Holman et al. 2010) has demonstrated that stable systems are capable of producing TTVs, and other studies have explored a limited range of architectures (Haghighipour & Kirste 2011), but the systematic exploration of parameter space of an Earth-mass planet in orbit near a hot Jupiter has not been undertaken. Here we examine 3.6 million possible masses and orbits of an approximately Earth-mass planet orbiting in or near the 2:1 outer mean motion resonance through N-body simulations. We find that stable orbits that also produce detectable TTVs do exist in the XO-2 system, and hence can exist in similar Hot-Jupiter systems.

In this paper we report observations of 10 transits of XO-2b, taken as part of APOSTLE. In Section 2 we outline our observations. In Section 3 we briefly outline (1) the data reduction, photometry (Section 3.1), (2) the transit model (Section 3.2) and light curve fitting (Section 3.3); both processes have also been described in previous work (Kundurthy et al. 2013). In Section 4 we present our estimates of the system parameters for XO-2b and in Sections 4.1 and 4.2 we present results from our study of transit depth variations (TDVs) and TTVs. In Section 5 we present results form N-body simulations used to study the stability of hypothetical planetary configurations at the 2:1 mean-motions resonance. Finally, in Section 6 we summarize our findings.

XO-2b was observed by members of the APOSTLE team on 10 occasions over a timespan of 3 yr from early 2008 until the spring of 2011. All observations of XO-2 were carried out using AGILE, a high-speed frame-transfer CCD (Mukadam et al. 2011), on the ARC⁴ 3.5 m telescope at Apache Point, New Mexico. The summary of observations is given in Table 1. XO-2 was an early APOSTLE target and was observed using a variety of instrumental settings, as the team had not converged on an optimal observing strategy prior to 2010. Early observations were made in the I-band (λ₀ = 805 nm; Cousins 1976; Bessell 1990) with several data sets taken with short read-out (exposure times, Column 5 in Table 1). The short read-out mode allowed for fine sampling of the light curve, but, due to the lower signal-to-noise and the unsuitability for characterizing systematics, this observing mode was abandoned. In addition, the I-band images also contained a strong contribution from a fringe pattern due interference from the backscattering of atmospheric lines within the CCD's pixels. The removal of this fringe pattern is discussed in Section 3.1. Early in 2010, the observing strategy changed to longer readouts (typically 45–75 s) to reduce the level of correlated noise (red-noise), and were made using the r'-band  which is similar to the SDSS⁵ r filter (λ₀ = 626 nm; Fukugita et al. 1996), to reduce the influence of fringes. In the long read-out mode the telescope was defocused to spread the stellar point-spread function (PSF) across multiple pixels, which minimized the systematics caused by pixel-to-pixel wandering of the PSF over the imperfect flatfield. The longer exposures also allowed for a greater count rate which maximized the signal-to-noise per image. The count rate was kept below AGILE's non-linearity limit of ∼52k ADU and well below its saturation level of 61k ADU by small adjustments to the telescope's secondary focus during observations.

XO-2 (TYC 3413-5-1) and its visual binary companion TYC 3413-210-1 (separated by ∼30'') were the only bright stars that fit in AGILE's field of view. Both stars are of identical spectral type (K0V) and nearly identical brightnesses in the filters used by APOSTLE, with their Johnson R and I magnitudes at 10.8 and 10.5, respectively (Monet et al. 2003). Our uncalibrated differential photometry also showed good agreement in their brightnesses, with the out-of-eclipse, un-normalized flux ratios being different by only 4% and 3% in the I-band and r'-band, respectively (Column "Flux Norm." in Table 1). The observations were made over a variety of observing conditions (Column "Obs. Conditions" in Table 1). The observing conditions are classified as "Clear" or "Poor Weather" with the former implying good data with few or no interruptions in data collection, and the latter indicating that the we experienced cloud cover or poor seeing conditions resulting in lower quality data. The tabulated transits are those for which we were able to capture the whole transit or at least a partial transit. Partial transits are those where portions of the in-eclipse light curve were lost due to bad weather or instrumental failure. Several data points were lost for the nights of UTD 2008 September 22 (transit 4), 2009 February 6 (transit 5), 2009 March 12 (transit 6), and 2010 December 27 (transit 8). However, we do include these nights since we did manage to obtain reasonable portions of the in-eclipse and out-of-eclipse data, which make it possible to determine transit properties (albeit at a loss of accuracy and precision).

This section outlines various stages in the analysis of light curves, starting with (1) the image reduction and photometry, and (2) the transit model and Markov Chain Monte Carlo (MCMC) analyzer (see also Kundurthy et al. 2013).

3.1. Reduction and Photometry

APOSTLE data were reduced using a pipeline developed specifically for AGILE images. The pipeline (written in IDL⁶) performs pixel-by-pixel error propagation, and image processing specific to the AGILE CCD. In addition to the standard photometric reduction steps like dark subtraction and flat-fielding, the pipeline performs non-linearity and fringe corrections specific to AGILE. The details on AGILE's non-linearity correction are described in Section 3 of Kundurthy et al. (2013). Some of the initial observations by APOSTLE were carried out using short exposures (see Column 5 "Exp.") in Table 1. We binned these data by averaging the flux ratios in 45 s bins.

For several of the initial XO-2 data taken using AGILE's I-band, a fringe pattern had to be subtracted to create science images. Photons from strong atmospheric lines (in the I-band bandpass) backscatter within the CCD pixels, and owing to the variable thickness of the pixel array, interference between these photons creates the fringe pattern. Since the fringe pattern is convolved with the illumination pattern of the CCD, the fringe correction has to be applied after dark-subtraction and flat-fielding. During the initial characterization of the AGILE CCD, an empirical fringe frame was produced by median combining dithered frames on the dark sky, where the flux contributions from stars were removed by outlier rejection. The resulting combined frame served as a "Model" fringe pattern (F), normalized to have a median of zero, and an amplitude of one, such that it could be scaled to match the fringe patterns on science frames. The science frame affected by the fringe pattern ($T^{^{\prime }}$) is assumed to be a linear combination of the fringe-less science image (T) and a fringe pattern:

$$\begin{equation} T^{^{\prime }} = T + a_{0} F, \end{equation} \tag{ 1 }$$

where F is the model fringe frame and a₀ is a scaling factor describing the amplitude of the fringing on a given frame. The fringe amplitude is estimated by minimizing χ² and fitting for a₀ in the above model. The corrected images ($T = T^{^{\prime }} - a_{0}F$), were found to be sufficiently corrected of fringes after visual examination. We found the fringe amplitudes (a₀) on XO-2 science frames to always be smaller than the standard deviation (i.e., scatter) of the global sky background on each frame; a₀ ranged between 6%–13% of the scatter in the sky for XO-2 data.

We extracted photometry from an optimal circular aperture centered around the target and comparison stars. In addition we extracted the counts on image products like the master-dark and master-flat, to serve as nuisance parameters for detrending, using the same aperture and centroids from photometry. The centroids were derived using SExtractor (Bertin & Arnouts 1996), which allowed for the use of customized donut-shaped convolution kernels for defocused PSFs. Science frames where pixels inside a photometric aperture exceeded AGILE's saturation limit of 61k were rejected. Images at the other extreme, where the stars were obscured by clouds, and resulted in low signal-to-noise measurements were also rejected (i.e., where individual photometric errors were >5000 ppm). The fraction of rejected frames per night is listed in Column 9 "%Rej." in Table 1. The optimal aperture was selected after extracting photometry on a list of circular apertures with radii between 5–50 pixels at an interval of 1 pixel. The optimal aperture was selected where the scatter in the residuals of the detrended light curve minus a trial transit model (based on values from the literature) was minimized. The correction function F_cor (or detrending function) is modeled as a linear sum of nuisance parameters as described by the following equation:

$$\begin{equation} F_{{\rm cor},i} = \sum _{k=1}^{N_{{\rm nus}}} c_{k} X_{k,i}, \end{equation} \tag{ 2 }$$

where X_{k, i} are the nuisance parameters, c_k are the corresponding coefficients. The index k counts over the number of nuisance parameters N_nus. The detrending coefficients are chosen by minimizing the χ² between the observed data (O), a model function (M), and the correction function,

$$\begin{equation} \chi ^2 = \sum _{j}^{N_{{\rm all}}} \frac{(O_{j} - M_{j} - F_{{\rm cor},j})^2}{\sigma _{j}^2}, \end{equation} \tag{ 3 }$$

here j is the index which counts over the total number of data points (N_all). The observed, model and correction function terms are all in normalized flux ratio units. The list of parameters used for detrending each light curve in the XO-2 dataset are presented in Table 2, with variable definitions given in the footnotes. A suitable set of detrending parameters were selected for a given night by running several manual trials of linear least squares minimization on single transit light curves and its corresponding model light curve and correction function (Equation (2)). The model parameters were fixed at reasonable values and only the coefficients (c_k) were fit. Those nuisance parameters which returned large uncertainties on the coefficients were excluded since this indicated a poor match to the noise trend in a light curve. This method of selecting nuisance parameters is admittedly ad-hoc; automated parameter selection techniques could be applied in the future. The detrending is performed using the final set of nuisance parameters in conjunction with fitting the transit parameters in order to ensure that the correction process is not biased by the trial model used in the selection of nuisance parameters.

Table 2. Nuisance Parameters Used for XO-2 Nights

T#	Nuisance Parameters Used (FLDC, OLDC, and MDFLDC)
T1	airmass, msky1, msky2, gsky, x1, y1, x2, y2, sD1, sD2, sF1, sF2, a0
T2	airmass, msky1, msky2, gsky, x1, y1, x2, y2, sD1, sD2, sF1, sF2, a0
T3	airmass, msky1, msky2, gsky, x1, y1, x2, y2, sD1, sD2, sF1, sF2, a0
T4	airmass, msky1, msky2, gsky, x1, y1, x2, y2, sD1, sD2, sF1, sF2, a0
T5	airmass, x1, y1, sD1, sF1, a0
T6	airmass, x1, y1, sD1, sF1, a0
T7	airmass, gsky, x1, y1, sD1, sF1
T8	airmass, gsky, x1, y1, sD1, sF1
T9	airmass, gsky, x1, y1, sD1, sF1
T10	airmass, gsky, x1, y1, sD1, sF1

Notes. airmass: atmospheric column; (x₁, y₁), (x₂, y₂): centriods of Target (1) and Comparison (2); msky₁, msky₂: median sky around Target (1) and Comparison (2); gsky: global median sky; sD₁, sD₂: sum of counts in aperture on master-dark; sF₁, sF₂: sum of counts in aperture on master-flat; a₀: fringe scaling (see Equation (1)).

Download table as:  ASCIITypeset image

The 10 transits of XO-2b are shown in Figure 1 in normalized flux ratios (with offsets for clarity). Shown are the detrended lightcurves, the corresponding errors, and the best fit transit model; the data shown in Figure 1 are listed in Table 3. The plotted data result from the data reduction and model fitting processes described in Sections 3.1–3.3.

Zoom In Zoom Out Reset image size

Table 3. APOSTLE Light Curve Data for XO-2

T#	T − T0	Norm. Fl. Ratio	Err. Norm. Fl. Ratio	Model Data
(1)	(2)	(3)	(4)	(5)
1	−0.1341616	1.000004100	0.001207284	1.000000000
1	−0.1315574	0.999797978	0.001102731	1.000000000
1	−0.1305157	0.997707769	0.001047412	1.000000000
1	−0.1299949	0.999906661	0.001190794	1.000000000
1	−0.1289532	0.998693635	0.001118320	1.000000000
1	−0.1279116	0.998939727	0.001098241	1.000000000
1	−0.1273907	0.999021973	0.000987372	1.000000000
1	−0.1268699	1.001381660	0.001209313	1.000000000
1	−0.1263491	0.998899046	0.001076096	1.000000000
1	−0.1258282	1.001212353	0.001160014	1.000000000
.	.	.	.	.
.	.	.	.	.

Notes. Column 1: transit number; Column 2: time stamps − mid transit times (BJD); Column 3: normalized flux ratio; Column 4: error on normalized flux ratio; Column 5: model data.

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

3.2. MultiTransitQuick

3.3. Transit MCMC

This section describes results from our execution of the two chains for the parameter set $\boldsymbol{\theta }_{{\rm Multi\hbox{-}Filter}}$, and one for the $\boldsymbol{\theta }_{{\rm Multi\hbox{-}Depth}}$ described in Section 3.3. Post processing statistics and other data for these chains are listed in Table 4. The columns "N_free," "Chain Length," "Corr. Length" and "Eff. Length" list the number of free parameters, the length of the cropped chain, the correlation and effective lengths, respectively. All chains were run for approximately 2 million steps, but about 100,000 of the initial steps were removed to account for "burn-in" and selection rate stabilization. The XO-2b chain with OLDC has a low effective length indicating poor Markov chain statistics. The FLDC chains (both $\boldsymbol{\theta }_{{\rm Multi\hbox{-}Filter}}$ and $\boldsymbol{\theta }_{{\rm Multi\hbox{-}Depth}}$) model satisfy the condition of a well sampled posterior distribution (effective length is >1000). The final two columns list the goodness of fit (i.e., lowest χ² in the MCMC ensemble) and degrees-of-freedom (DOF) from the respective chain. Parameters from all chains had Gelman–Rubin $\hat{\rm R}$-statistics close to 1 indicating that the parameter space was covered evenly (though the OLDC chain was not sampled finely enough, based on the auto-correlation data).

Table 4. TMCMC Chains for XO-2

Chain	Model Vector	Nfree	Chain Length	Corr. Length	Eff Length	χ2	DOF
FLDC	$\boldsymbol{\theta }_{{\rm Multi\hbox{-}Filter}}$	14	1,900,001	100	19,000	3272.94	3312
OLDC	$\boldsymbol{\theta }_{{\rm Multi\hbox{-}Filter}}$	18	1,900,001	4,949	383	3268.48	3308
MDFLDC	$\boldsymbol{\theta }_{{\rm Multi\hbox{-}Depth}}$	22	1,900,001	588	3,231	3197.53	3304

Download table as:  ASCIITypeset image

The resulting best-fit parameter estimates are listed in Table 5 for the Multi-Filter models, and in Table 6 for the Multi-Depth models. These tables also list the derived system parameters. The transformation between the MTQ parameters to the derived system parameters are described in Carter et al. (2008) and Kundurthy et al. (2011). Contour plots showing the joint probability distributions (JPDs) for the fit and derived parameters are shown in Figures 2 and 3, respectively. There are no strong correlations between the fit parameters t_T, t_G and the transit depths, D_I and $D_{r^{\prime }}$, as seen in Figure 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 5. XO-2 Parameters for $\boldsymbol{\theta }_{{\rm Multi\hbox{-}Filter}}$

Parameter	FLDC	OLDC	Unit
MTQ Parameters
tG	0.0107 ± 0.0002	0.0108 ± 0.0002	days
tT	0.1008 ± 0.0001	0.1004 ± 0.0002	days
D(I)	0.0126 ± 0.0001	0.0126 ± 0.0001	...
$D_{{(r^{\prime })}}$	0.0130 ± 0.0001	0.0131 ± 0.0001	...
v1(I)	(0.5944)	0.5440 ± 0.0313	...
v1(r')	(0.6980)	0.6411$^{+0.0332}_{-0.0249}$	...
v2(I)	(0.1452)	0.2806 ± 0.0904	...
v2(r')	(0.3524)	0.4786$^{+0.0734}_{-0.0992}$	...
Derived Parameters
(Rp/R⋆)(I)	0.1030 ± 0.0003	0.1033 ± 0.0004	...
$(R_{p}/R_{\star })_{{(r^{\prime })}}$	0.1024 ± 0.0003	0.1029 ± 0.0004	...
b	0.17$^{+0.04}_{-0.02}$	0.21 ± 0.03	...
a/R⋆	8.14 ± 0.06	8.11 ± 0.07	...
iorb	88.79 ± 0.15	88.53$^{+0.02}_{-0.10}$	(deg)
ν/R⋆	19.55 ± 0.15	19.48 ± 0.17	days−1
ρ⋆	1.49 ± 0.03	1.48 ± 0.04	g cm−3
P (2.6159 days +)	−3467 ± 22	−3466 ± 22	ms

Download table as:  ASCIITypeset image

Table 6. XO-2 Parameters for $\boldsymbol{\theta }_{{\rm Multi\hbox{-}Depth}}$

Transit Depths	Value	Units	Rp/R⋆	Value	Units
		I-band
(D)1	0.01338 ± 0.00013	...	(Rp/R⋆)1	0.1061 ± 0.0005	...
(D)2	0.01300 ± 0.00017	...	(Rp/R⋆)2	0.1047 ± 0.0007	...
(D)3	0.01226 ± 0.00013	...	(Rp/R⋆)3	0.1019 ± 0.0005	...
(D)4	0.01199 ± 0.00013	...	(Rp/R⋆)4	0.1009 ± 0.0005	...
(D)5	0.01209 ± 0.00039	...	(Rp/R⋆)5	0.1013 ± 0.0015	...
(D)6	0.01267 ± 0.00015	...	(Rp/R⋆)6	0.1035 ± 0.0006	...
		r'-band
(D)7	0.01326 ± 0.00009	...	(Rp/R⋆)7	0.1032 ± 0.0004	...
(D)8	0.01320 ± 0.00015	...	(Rp/R⋆)8	0.1030 ± 0.0006	...
(D)9	0.01307 ± 0.00009	...	(Rp/R⋆)9	0.1026 ± 0.0003	...
(D)10	0.01300 ± 0.00013	...	(Rp/R⋆)10	0.1024 ± 0.0005	...
Other MTQ Parameters
Parameter	Value	Units	Parameter	Value	Units
tG	0.0108 ± 0.0001	days	tT	0.1008 ± 0.0001	days
v1(I)	(0.5944)	...	v1(r')	(0.6980)	...
v2(I)	(0.1452)	...	v2(r')	(0.3524)	...
Derived Parameters
b	0.24 ± 0.02	...	a/R⋆	8.02$^{+0.03}_{-0.04}$	...
iorb	88.29 ± 0.15	o(deg)	ν/R⋆	19.27$^{+0.08}_{-0.10}$	days−1
ρ⋆	1.43 ± 0.02	g cm−3	...	...	...

Download table as:  ASCIITypeset image

System parameters agree with previously published values in the literature, as seen by the overlap of the uncertainties in the JPD plot (Figure 3). In addition APOSTLE's measurements give tighter constraints on several of the system parameters. The errors on a/R_⋆, i_orb, the impact parameter (b) and the stellar density (ρ_⋆) are more precise than previous measurements by a factor of ∼3 (Burke et al. 2007; Fernandez et al. 2009; Sing et al. 2011). However, previous studies using the combination of the TMCMC Markov chain analyzer and MTQ transit model have resulted in underestimated errors since we do not include rednoise analysis (Carter & Winn 2009). Hence, the errors presented in Table 6 and Figure 3 do not truly reflect improvements in the measurements of parameters. More conservative constraints were placed on a subset of these system parameters using the TAP package. Comparisons of some parameters and their uncertainties are presented in Table 7.

Table 7. Comparison of Estimates of System Parameters for XO-2b

Parameter	TMCMC	TAP	B07	F09	S11	C12	Units
a/R⋆	8.14 ± 0.06	8.13 ± 0.10	8.20 ± 0.20	8.13 ± 0.20	7.83 ± 0.17	7.99 ± 0.07	...
iorb	88.79 ± 0.15	88.80 ± 0.61	88.90 ± 0.75	...	87.62 ± 0.51	88.01 ± 0.33	(deg)
b	0.172$^{+0.040}_{-0.021}$	0.171 ± 0.085	0.158 ± 0.110	0.160 ± 0.110	0.324 ± 0.070	0.280 ± 0.044	...
ρ⋆	1.49 ± 0.03	1.49 ± 0.05	1.52 ± 0.11	1.48 ± 0.10	1.33 ± 0.09	1.41 ± 0.04	g cm−3

Notes. TMCMC and TAP values are from independent analysis of APOSTLE light curves; B07: Burke et al. 2007; F09: Fernandez et al. 2009; S11: Sing et al. 2011; C12: Crouzet et al. 2012.

Download table as:  ASCIITypeset image

It is clear that using TAP on the APOSTLE dataset and accounting for rednoise provides more conservative estimates of the system parameters when compared to TMCMC values. TAP errors for a/R_⋆ and ρ_⋆ are better than those reported by Sing et al. (2011) by up to ∼40%, whose observations were made using the 10.4 m Gran Telescopio Canarias (GTC). One must note the Sing et al. (2011) data result from narrow-band photometry and hence would have much lower photometric precision when compared to broadband observations from the same telescope. The resulting TAP errors for the impact parameter (b) and the orbital inclination are larger than those from the GTC study by ∼20%. Thus we report only some improvements in the measurements of system parameters. Our system parameters agree with the current best estimates reported in Crouzet et al. (2012), but our error bars are larger by factors ⩽2.

4.1. Transit Depth Analysis

Figure 4 shows transit depth versus transit epoch for the I-band (top panel) and r'-band (bottom panel). The overall variations in the I-band depth are ∼0.05% compared to the 0.01% uncertainty in D_(I) from the joint fit to depths in the FLDC chain. Depth variations can be caused by spots. Even though the variations we present are significant we refrain from claiming the detection of spot-modulation. These deviations are likely due to the incomplete sampling of several transits. Spots influence stellar brightness to a greater extent at shorter wavelengths, so the r'-band would be more conducive to showing depth variations. However, the overall depth variation in the r'-band light curves is of the order of 0.01% and is consistent with the errors on $D_{(r^{\prime })}$ from the fit reported by the FLDC chain (Table 5). The variations seen in the I-band depths are difficult to explain. They can either be due to (1) real brightness variations caused by spot-modulations, or (2) variations arising from the transit model's inability to accurately constrain transit depths and errors given the incomplete sampling of light curves (like transits 4–6). We note that errors on the transit depth from our Multi-Depth LDC chains are more sensitive to the incompleteness of in-eclipse data rather than the photometric precision of a given light curve. For example, transit numbers 4 and 6 have 86% and 72% of the in-eclipse transit data sampled, respectively; they also have light curve residuals of 939 ppm and 405 ppm, respectively. Thus even though transit 6 has significantly better photometric precision, the depth for transit 4 is constrained slightly better due to the fact that more in-eclipse points are available. A more dramatic difference can be seen with transit 5, which has only 47% of the in-eclipse light curve sampled and the worst depth constraint. Other than transit 5 all transits have eclipses sampled to better than 70%, hence the depth uncertainties are all comparable. Even though the Multi-Depth LDC chain satisfied the convergence criteria described in Section 3.3, the effective length of the chain is far lesser than the value derived for the Multi-Filter LDC chain (see Table 4), indicating the lower significance of these MCMC results. In addition one must note the lack of rednoise analysis in TMCMC, which implies that the errors presented Figure 4 are probably underestimates.

Zoom In Zoom Out Reset image size

4.2. Transit Timing Analysis

Several planetary systems have been observed to have TTVs (Holman et al. 2010; Lissauer et al. 2011a; Ballard et al. 2011; Nesvorný et al. 2012). In certain configurations these variations can be on the order of minutes and can be easily seen in an diagram showing observed minus computed transit times (i.e., the difference between the measured times and times expected from purely Keplerian orbital periods). Using TTVs to look for additional planets was first proposed by Agol et al. (2005) and Holman & Murray (2005), who showed that unseen planetary siblings can gravitationally influence the orbits of a known transiting planet. The TTV signals are known to be especially strong if the unseen companion lies close to mean motion resonance with the transiting planet.

Available transit times from the literature for XO-2b include those of the discovery paper by Burke et al. (2007), and the follow-up observations by Fernandez et al. (2009), Sing et al. (2011), and Crouzet et al. (2012). We excluded the 11 transit times from Burke et al. (2007) since the reported timing errors were on the order of ∼3 minutes, and are too large to provide meaningful constraints on the orbital ephemeris. We do include transit times from the other three studies. The time coordinate BJD (TDB) has become the standard system used for transit timing studies (Eastman et al. 2010) and all transit times used in this study were brought to this system. The timestamps of all APOSTLE data were converted to BJD (TDB) in the customized reduction pipeline (Kundurthy et al. 2011). The APOSTLE pipeline's time conversions have been verified by comparison to the commonly used time conversion routines made available by Eastman et al. (2010).

The observed minus computed mid-transit times (O−C) are plotted in Figure 5. There are two versions of APOSTLE transit times presented, one from the FLDC chain (TMCMC + MTQ), and the other from the TAP fit to APOSTLE light curves. These transit times are tabulated in Table 8. A linear ephemeris was fit to all the data (including literature values) using the equation,

$$\begin{equation} T_i = T0 + {\rm Epoch}_i \times P, \end{equation} \tag{ 4 }$$

resulting in a best fit ephemeris of

$$\begin{eqnarray*} P = 2.61585988 \pm 0.00000016 & \; {\rm days} \\ T0 = 2454406.720516 \pm 0.000046 & \; {\rm BJD} \end{eqnarray*}$$

with a goodness of fit χ² = 105.14 for 31 DOF. The reduced chi-squared χ²/ν = 3.39 indicates that the linear fit is not robust, and is either due to the large timing deviations in the data or underestimated errors. For example, several of the Sing et al. (2011) data points lie far from the zero O−C line, with the largest deviation being ∼106 s The scatter in the O−C values as a whole is ∼39 s which makes the 106 s deviation fall within the 3σ confidence interval of the collective data set. Though the linear ephemeris does not precisely fit the timing data, the level of variation is not significant enough for us to claim unseen planets as the cause.

Zoom In Zoom Out Reset image size

Table 8. APOSTLE Transit Times for XO2

Epoch	T0 (TMCMC)	σT0	T0 (TAP)	σT0
	2,400,000+ (BJD)	(BJD)	2,400,000+ (BJD)	(BJD)
26	54474.73242	0.00011	54474.73252	0.00021
39	54508.73875	0.00013	54508.73877	0.00020
47	54529.66596	0.00009	54529.66604	0.00016
148	54793.86845	0.00024	54793.86836	0.00041
177	54869.72753	0.00023	54869.72735	0.00050
190	54903.73397	0.00014	54903.73403	0.00020
416	55494.91796	0.00015	55494.91794	0.00021
440	55557.69891	0.00020	55557.69892	0.00025
453	55591.70507	0.00008	55591.70511	0.00015
466	55625.71093	0.00014	55625.71112	0.00019
Fit	Period (days)	σP	T0 (BJD)	σT0
TMCMC	2.615860095	±0.000000209	2454474.7327333	±0.0000599
TAP	2.615860014	±0.000000346	2454474.7327964	±0.0001028

Download table as:  ASCIITypeset image

In order to compare the ephemerides derived with and without rednoise analysis, we fit for a linear ephemeris to the APOSTLE transit times from TMCMC and TAP, respectively (presented in Table 8). The difference between the periods derived for these subsets and the period derived from all available transit times was <7 ms. The results from the fits to the TMCMC and TAP transit times are presented at the bottom of Table 8. The reduced χ²s for linear ephemeris fits to the TMCMC and TAP subsets were 3.96 and 1.26, respectively, confirming that TAP gives more conservative errors for the transit times due to the red-noise analysis.

We were also able to rule out sinusoidal trends in the data by running a generalized Lomb–Ccargle analysis on the O − C data (Zechmeister & Kürster 2009). We fit for a period, amplitude, and phase offset on two sets of O − C data. The first set included the literature dataset and all APOSTLE measurements, and the second included only the TAP measurements of APOSTLE. We found that a sinusoid of period ∼19 days and amplitude of ∼32 s improved the fit when the entire timing data set was used. The sinusoid fit yielded a χ² =81.8, which was an improvement of Δχ² = 23.3 compared to the linear ephemeris (see above). However, the reduced χ² remains greater than 1, making it a non-robust fit. The periodogram analysis on only the TAP data yielded a period and amplitude of 8.5 days and 27.4 s, respectively. The Δχ² showed improvement compared to the linear ephemeris model by ∼7.6. In a manner similar to Becker et al. (2013) we tested the significance of this sinusoidal fit by repeating the analysis on 10⁵ randomly cycled permutations of the amplitudes of the O − C measurements, keeping the epochs fixed. We found that 76.7% of the fits showed improved Δχ² at amplitudes greater than or equal to 27.4 s, indicating that the periodicity is not likely from a real TTV.

For the case when planets are in mean-motion resonance (MMR), Agol et al. (2005) showed that the analytic expression, δt_max ∼ (P/4.5j)(m_pert/(m_pert + m_trans)), can roughly estimate the amplitude of the timing deviation (δt_max). The quantities m_pert, m_trans, P and j are the mass of the unseen perturber, the mass of the transiting planet, the orbital period of the transiting planet and the order of the resonance, respectively. For the XO-2 system, we can rule out possible system configurations given the amplitude of the weak sinusoidal fit to the TTV data (∼27.4 s). Using the orbital period from Table 8 and the mass of XO-2b, m_trans = 0.565 M_Jup (Fernandez et al. 2009), we compute the maximum mass perturber that could exist in the XO-2 system in the 2:1 MMR to be ∼0.2 M_⊕, i.e., additional planets with M_p < 0.2 M_⊕ may exist near the 2:1 MMR given these data. At higher order resonances, this maximum mass (for a possible perturber) is larger. A more detailed analysis of TTV needs to account for a variety of orbital configurations of hypothetical companions, which is addressed in the following section.

In order to evaluate the stability of systems that produce a detectable TTV, we integrated 3.6 × 10⁶ orbital configurations designed to mimic the XO-2 system. Each trial consisted of the known planet XO-2b and a hypothetical terrestrial-like companion.

The parameter space we cover is presented in Table 9. In this table Δ is the interval between values, all of which were varied uniformly in the range between "Min" and "Max," producing N bins. In this table, m is mass of the hypothetical terrestrial exoplanet, e is eccentricity, i_orb is the orbital inclination measured from the plane perpendicular to the sky, Ω is the longitude of ascending node, ω is the argument of pericenter, and M is the mean anomaly. Since our timing precisions are at the level of a signal from a terrestrial planet, we limited the mass coverage to the "super-Earth" range of 1–10 M_⊕ planets. Our semi-major axis coverage spans the 2:1 resonance, with one choice directly at the commensurability. All other parameters were varied so as to fully cover all possible architectures, but still to keep the entire program tractable. Each simulation required ≳ 1 hour of CPU time, or ∼7 million CPU hours total.

For the N-body simulations, we used the symplectic integrator code in the HNBody package (Rauch & Hamilton 2002),⁷ which includes general relativity. We integrated each trial for 10⁶ orbits of XO-2b (∼7200 yr), which is enough to identify >95% of unstable orbits for other exoplanet system (Barnes & Quinn 2004). Each simulation could produce one of three outcomes: stable, unstable, or fail. Stable configurations lost no planets due to gravitational perturbations, while unstable ones did. Failed systems did not conserve energy to better than 1 part in 10⁴, which is required for symplectic integrators (Barnes & Quinn 2004).

In Figure 6, a sample slice through the data set is shown. We limited the visualization to cases with m ∼ 1 M_⊕ planets with i_orb = Ω = ω = 0, and M = 144°. Green squares designate stable configurations, red unstable, and yellow failed. Most simulations from the yellow bins are probably unstable, as close approaches between planets can violate the algorithm's underlying assumption that the gravitational force from the star is much larger than that from any other planet. We note that the average fractional change in energy (dE/E) for stable cases was 7 × 10⁻⁹.

Zoom In Zoom Out Reset image size

The trials shown in Figure 6 show the that, for this value of M, the mean motion resonance at 0.0587 AU stabilizes the system. At all values of a and low e, the system is stable, as predicted by Hill stability theory (see Marchal & Bozis 1982; Gladman 1993; Barnes & Greenberg 2006; Kopparapu & Barnes 2010). However at the resonance, stability is likely for e ⩽ 0.7. As expected, we find that only certain values of M predict this "tongue" of stability.

The contour lines in Figure 6 show the values of the TTV signal in seconds. For the broad stable region at e ≲ 0.1, the signal is at or below the detection limit. For larger e, especially in the resonance, the value can be much larger. For the bin at a = 0.0587 AU, and e = 0.7, the signal magnitude is close to 2377 s, or nearly 40 minutes. Clearly, the APOSTLE project could have detected an Earth-mass companion if it were in a favorable orbit. We do note, however, that an Earth-like planet with large eccentricity is likely to be rapidly tidally circularized due to the close proximity (Rasio et al. 1996; Jackson et al. 2008; Barnes et al. 2010), although some eccentricity would be maintained by the resonance.

As noted in Table 8 (Section 4.2) the median TTV errors obtained from APOSTLE light curves of XO-2b were ∼20 s. Of the 3.6 million parameter configurations, ∼1.14 million were found to be stable, and approximately 51% of these stable configurations produce TTVs on the order of 20 s or less and thus would not be easily detected by APOSTLE or other surveys with similar capability. The vast amount of data generated by the stability simulations cannot be adequately discussed in this text, hence we provide the results from all 3.6 million simulations as part of an online-only table, a segment of which is displayed in Table 10. The first column "Sim. ID" denotes the simulation number and columns numbered 2–8 are the input configurations for a given simulation. The range of these simulation configurations were summarized earlier in this section and in Table 9. Columns 9–13 are the simulation statistics. Column 9 indicates the stability outcome of a given simulation, with 0 or 1 denoting unstable and stable simulations respectively. Column 10 lists the energy conservation of a given simulations as noted by the maximum fractional change in energy over a given step (dE/E). Columns 11–13 list the standard deviations in the transit timing variations, transit duration variations and the variations in the impact parameter, (b) respectively.

• 1.  
  Photometric precision. The XO-2 system was observed over a variety of observing conditions over a period of 3 yr between 2008–2011. The best photometric precisions we obtained with our I-band and r'-band observations were 405 and 354 ppm, respectively (see Table 1).
• 2.  
  XO-2b system parameters. Our analysis of the 10 transit light curves yielded estimates of system parameters that agree with measurements presented by other studies (see Figure 3). We were able to improve the constraints on a/R_⋆ and ρ_⋆ by ∼40% compared to the previous best measurements from the ground (Sing et al. 2011). The measurements are presented in Table 7; see the TAP values. We could not get the Markov chains to converge while fitting for stellar limb-darkening parameters, echoing the results from previous studies.
• 3.  
  Search for transit depth variations. Our Multi-Depth fits show some variations in the transit depth over transit epoch (see Figure 4 and Table 6) for the six I-band light curves. Since three out of the six light curves were not fully sampled we cannot confidently assert real variability in the data. The seen variations could be shortcomings of the transit model's ability to fit a set of light curves with both complete and incomplete data. We do not see similar variations in the r'-band, where one may expect spot-modulated variability to appear, due to the greater spot-to-star contrast in the r'-band. There are no known reports of stellar activity on XO-2, hence the r'-band results are consistent with this fact.
• 4.  
  Search for transit timing variations. The XO-2b dataset contains light curves with some of the best photometric precisions achieved with ground-based observations. Since photometric precision directly translates to transit timing precision we are able to report timing precisions as low as 12 s (after red-noise analysis, see Table 8). We were unable to detect significant timing deviations for XO-2b in our data. The linear fit was not robust, with χ²/DOF being significantly greater than 1, indicating large scatter around the linear ephemeris fit. The transit times derived from the APOSTLE light curves using the red-noise analysis of Carter & Winn (2009) resulted in more conservative errors than those derived using TMCMC. A linear ephemeris is consistent with the transit time measurements reported from the TAP analysis. The overall variation in the O − C values from our rednoise analysis was on the order of ∼39 s. From a sinusoidal fit to APOSTLE's O − C data, we obtained a TTV amplitude of ∼27.4 s. However, checking the goodness-of-fit of sinusoids to random rearrangements of the data show that the detected periodicity unlikely to be real. The resulting amplitude rules out planetary companions more massive than 0.2 M_⊕ near the 2:1 MMR, and larger companions near higher order resonances.We conclude that the set of transit times published in the literature for XO-2b and other transiting systems in general are not suited for transit timing analysis. Lacking red-noise analysis leads to underestimated timing errors and may lead to premature reporting of timing variations. A proper analysis of transit times would need a simultaneous analysis of transit light curves using a transit model which is (1) suited for Bayesian inference (i.e., with a fairly uncorrelated parameter set; Carter et al. 2008) and (2) a transit model which can adequately account for red-noise in the data (like TAP; Gazak et al. 2012). Using large data sets of transit light curves may be inefficient due to the slowness of Markov chains with the addition of model parameters. In addition to the development of more detailed models, the utilization of fast Markov chain algorithms (for e.g., Foreman-Mackey et al. 2013) is also recommended.Our lack of detection of TTVs in the data for the Hot-Jupiter XO-2b is also consistent with Kepler's findings that (1) Hot-Jupiters tend to lack other planetary siblings (Latham et al. 2011; Steffen et al. 2012) and (2) members of multi-planet systems with short period planets (Period <10 days) are more likely to be Hot-Neptunes (Latham et al. 2011; Lissauer et al. 2011b).
• 5.  
  Dynamical study. We ran 3.6 million N-body simulations of possible multi-planet configurations near the 2:1 resonance in order to test for (1) orbital stability and (2) detectability of TTV signals. We varied several properties of the hypothetical companion to XO-2b (see Table 9), to look for patterns in the resulting stability and transit times of the simulations over parameter space. Of the several stable configurations, we find that ∼51% of the simulations would display TTVs weaker than the precision limits of our survey and hence remain undetectable. The entire table of simulation statistics are presented online (see Table 10).