Issue 
A&A
Volume 555, July 2013



Article Number  A92  
Number of page(s)  12  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201321407  
Published online  05 July 2013 
Qatar1: indications for possible transit timing variations ^{⋆}
^{1}
Hamburger Sternwarte, University of Hamburg,
Gojenbergsweg 112,
21029
Hamburg,
Germany
email:
cessen@hs.unihamburg.de
^{2}
Dept. of Astronomy, Box 351580,
University of Washington,
Seattle,
WA
98195,
USA
Received:
5
March
2013
Accepted:
22
May
2013
Aims. Variations in the timing of transiting exoplanets provide a powerful tool for detecting additional planets in the system. Thus, the aim of this paper is to discuss the plausibility of transit timing variations (TTVs) on the Qatar1 system by means of primary transit light curves analysis. Furthermore, we provide an interpretation of the timing variation.
Methods. We observed Qatar1 between March 2011 and October 2012 using the 1.2 m OLT telescope in Germany and the 0.6 m PTST telescope in Spain. We present 26 primary transits of the hot Jupiter Qatar1b. In total, our light curves cover a baseline of 18 months.
Results. We report on indications for possible longterm TTVs. Assuming that these TTVs are true, we present two different scenarios that could explain them. Our reported ~190 days TTV signal can be reproduced by either a weak perturber in resonance with Qatar1b, or by a massive body in the brown dwarf regime. More observations and radial velocity monitoring are required to better constrain the perturber’s characteristics. We also refine the ephemeris of Qatar1b, which we find to be T_{0} = 2456157.42204 ± 0.0001 BJD_{TDB} and P = 1.4200246 ± 0.0000007 days, and improve the system orbital parameters.
Key words: methods: data analysis / methods: observational / techniques: photometric / eclipses
Tables of the transit observations are only available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/qcat?J/A+A/555/A92
© ESO, 2013
1. Introduction
The discovery of Neptune in 1846 was a milestone in astronomy. The position at which the planet was detected by Galle had been predicted – independently – by Le Verrier and Adams, who attributed the observed irregularities of Uranus’ orbit to the gravitational attraction of another perturbing body outside Uranus’ radius. Thus, Neptune became the first planet to be predicted by celestial mechanics before it was directly observed.
Similar to the case of Neptune, in the realm of exoplanets an additional planet or exomoon can reveal itself by its gravitational influence on the orbital elements of the observed planet. Most notably, for transiting planets a perturber is expected to induce shortterm transit timing variations (TTVs; Holman & Murray 2005; Agol et al. 2005); because timing measurements can be carried out quite accurately, the TTV searches are quite sensitive, and it is possible to detect additional objects in the stellar system and to derive their properties with TTVs.
The detection of TTVs in groundbased measurements requires both a sufficiently long baseline and good phase coverage. The search for planets by TTVs has been a major activity in groundbased exoplanet research. So far, TTVs have been claimed in WASP3b (Maciejewski et al. 2010), WASP10b (Maciejewski et al. 2011), WASP5b (Fukui et al. 2011), HATP13 (Pál et al. 2011; Nascimbeni et al. 2011, but see Fulton et al. 2011), and OGLE111b (Díaz et al. 2008, but see Adams et al. 2010).
Summary of our observations carried out with OLT (top) and PTST (bottom).
In recent years, the search for TTVs in extrasolar planetary systems has entered a new era, marked by the advent of the spacebased observatories CoRoT and Kepler, which provide photometry of unprecedented accuracy. On one hand, the Kepler team has already presented more than 40 TTVs in exoplanetary systems, such as Kepler9 (Holman et al. 2010), Kepler11 (Lissauer et al. 2011), Kepler19 (Ballard et al. 2011), Kepler18 (Cochran et al. 2011), and Kepler 29, 30, 31 and 32 (Fabrycky et al. 2012). On the other hand, Steffen et al. (2012) searched for planetary companions orbiting hotJupiter planet candidates and found that most of the systems show no significant TTV signal. Yet, the Kepler satellite observes only a small area of the sky and has a limited lifetime, thus TTV will continue to play a major role also in groundbased studies.
The transiting Jovian planet Qatar1b was discovered by Alsubai et al. (2011) as the first planet found within The Qatar Exoplanet Survey. Qatar1A is an old (>4 Gy) K3V star with 0.85 M_{⊙} and 0.82 R_{⊙}, orbited by a closein (0.023 AU) planet with a mass of 1.1 M_{J} and a period of 1.42 d. Qatar1b has a radius of 1.16 R_{J} and an inclination angle of 83.47°, implying a nearly grazing transit. Therefore Qatar1 offers an outstanding opportunity to search for TTVs for several reasons: Its large inclination and closein orbit yield a highimpact parameter, making the shape and timing of the transit sensitive to variations of the orbital parameters. Its large radius and short orbital period allow for the study of a large number of deep transits over hundreds of epochs.
We have therefore observed transits in Qatar1b over the past few years, which we describe here. We specifically describe the observational setup in Sect. 2 as well as the data reduction process. We subsequently present the details of our light curve analysis and transit modeling in Sect. 3 and describe our TTV analysis in Sect. 4. We interpret our results in Sect. 5, where we present different dynamical scenarios that were tested to suggest the existence of an additional body in the system. In Sect. 6 we conclude.
2. Observations and data reduction
Our observations comprise 18 transits of Qatar1 obtained using the 1.2 m OskarLühning telescope (OLT) at Hamburg Observatory, Germany, and 8 transits using the 0.6 m Planet Transit Study Telescope (PTST) at the Mallorca Observatory in Spain.
The OLT data were taken between March 2011 and October 2012 using an Apogee Alta U9000 CCD with a 9′ × 9′ field of view. Binning was usually 4 × 4, but the first transit observations were taken in a 2 × 2 configuration. The binning was increased to reduce individual exposure times. The data were obtained with typical exposure times between 45 and 300 s depending on the night quality, the binning configuration, and the star’s altitude. All exposures were obtained using a JohnsonCousins Schuler R filter. Because Qatar1 is circumpolar at Hamburg’s latitude, typical airmass values range from 1 up to 1.9, while the average seeing value is 2.5 arcsec.
The PTST data were taken between May and August 2012 using a Santa Barbara CCD with a 30′ × 30′ field of view in a 3 × 3 binning and a Baade Rband filter setup. The exposure times range from 50 to 90 s. The observations were obtained with airmass values between 1.1 and 3 and typical seeing values of 2 arcsec. In Table 1 we summarize the main characteristics of our observations obtained at both sites. Combining OLT and PTST data, the observations cover a total of 416 epochs.
Calibration images such as bias and flat fields were obtained on each observing night. We used the IRAF task ccdproc for bias subtraction and flatfielding on the individual data sets, followed by the task apphot to carry out aperture photometry on all images including individual photometric errors. We measured fluxes using different apertures centered on the target star and six more stars with a similar brightness as Qatar1, which are present in both field of views. The apparent brightness of the reference star chosen to produce the differential light curves is very close to Qatar1. Multiband photometry of the same field of view reveals no significant difference between both stars as a function of photometric color, either suggesting that they are of similar spectral type. This minimizes any atmospheric extinction residuals on the differential light curves, and in addition, the apparent proximity of the two stars (~2 arcmin) minimizes systematic effects related to vignetting, comatic aberration, or CCD temperature gradients, which all increase with increasing distance from the telescope’s optical axis (i.e., the center of the chip). We also checked the constancy of the reference star against the other five comparison stars and chose as final aperture the one that minimized the scatter in the resulting light curves. The differential light curves were then produced by dividing the flux of the target star by that of the reference star.
Typical sky brightness values per binned pixel were of about 3500 counts for OLT, and 2000 for PTST. Figure 1 shows the PTST field of view (light background) superposed on the field of view of OLT (dark background). Qatar1 lies at the center of the field of view, marked with a large red circle. The five comparison stars are indicated with green squares and the reference star, indicated with a small blue circle, is the one used to produce the differential light curves.
After the reduction process, we fitted a straight line to the outoftransit data points to correct for any residual systematic trend and to normalize the differential light curves. It is worth mentioning that we calculated the timing offsets using both raw and normalized data. We therefore confirm that the normalization process does not produce any shifts in the midtransits but only more accurate planetary parameters. This might be because the amplitude of any systematic trend present in our light curves was smaller than ~2 mmag.
Fig. 1 OskarLühning telescope (dark background) and Planet Transit Search Telescope (light background) field of views. Qatar1 is indicated with a large red circle centered on the OLT field of view, along with the comparison stars as green squares. The small blue circle upward and to right of Qatar1 indicates the reference star used to produce the differential photometry. 
3. Data analysis
3.1. Data preparation
IRAF provides heliocentric corrections, therefore our time stamps are given as heliocentric Julian dates (HJD_{UTC}) and are converted into barycentric Julian dates (BJD_{TDB}) using the web tool provided by Eastman et al. (2010)^{1}.
3.2. Fit approach
We fitted the transit data with the transit model developed by Mandel & Agol (2002), making use of their occultquad FORTRAN routine^{2}. From the transit light curve, we can directly infer the following parameters: the orbital period P, the midtransit time T_{0}, the radius ratio p = R_{p}/R_{s}, the semimajor axis (in stellar radii) a/R_{s}, and the orbital inclination i. For our fits we assumed a quadratic limbdarkening law with fixed coefficients u_{1} and u_{2}.
3.3. Limbdarkening coefficients
Alsubai et al. (2011) presented a spectroscopic characterization of Qatar1 based on comparing their observed spectrum with synthetic template spectra and suggested that Qatar1 has an effective temperature T_{eff} of 4861 ± 125 K, a surface gravity log g of 4.536 ± 0.024, and solar metallicity [Fe/H] of 0.20 ± 0.10. Because our observations with OLT and PTST were obtained using different (nonstandard) filter sets, we decided to calculate angleresolved synthetic spectra from spherical atmosphere models using PHOENIX (Hauschildt & Baron 1999; Witte et al. 2009) for a star with effective temperature T_{eff} = 4900 K, [Fe/H] = 0.20 and log g = 4.5, thereby closely matching the spectroscopic parameters of Qatar1A. We then convolved each synthetic spectrum with the OLT and PTST filter transmission functions and integrated in the wavelength domain to compute intensities as a function of μ = cosθ, where θ is the angle between the line of sight and the radius vector from the center of the star to a reference position on the stellar surface. The thus derived intensities were then fitted with a quadratic limbdarkening prescription, viz. to obtain the u_{1} and u_{2} limbdarkening coefficients. As an aside, we note that the best approach to obtain limbdarkening coefficients would be to fit a more sophisticated biparametric approximation to the stellar intensities (Claret & Hauschildt 2003) to the PHOENIX intensities, which produces the smallest deviations in the generation of the limbdarkening coefficients, which fits the stellar limb better. However, to introduce such a limbdarkening law in the production of primary transit light curves would be computationally cumbersome. Moreover, Qatar1 is a nearly grazing system, where changes in limb darkening do not strongly affect the primary transit light curves. We checked that the time that the planet spends in the smallμ regime is very short. Furthermore, at a given time the planet covers a broad range of μvalues. Neglecting these points in determining the limbdarkening coefficients will probably not affect the subsequent parameter determination. To calculate the OLT and PTST limbdarkening coefficients, we fitted a quadratic limbdarkening law to PHOENIX intensities, neglecting the data points between μ = 0 and μ = 0.1. Figure 2 shows the OLT and PTST limbdarkening normalized functions, and Table 2 the fitted limbdarkening coefficients.
Fig. 2 Qatar1 normalized intensities considering the OLT and PTST filter transmission functions. Data points on the right of the vertical dashed line were not considered in the fitting procedure. 
Bestfit limbdarkening coefficients (LDCs) for OLT and PTST, along with the 1σ errors.
3.4. Parameter errors
To determine reliable errors for the fit parameters given the mutual dependence of the model parameters, we explored the parameter space by sampling from the posteriorprobability distribution using a Markovchain MonteCarlo (MCMC) approach.
The neargrazing transit geometry of Qatar1 introduces a substantial amount of correlation between the system parameters, which can easily render the MCMC sampling process inefficient. To cope with these difficulties we used a modification to the usually employed MetropolisHastings sampling algorithm, which is able to adapt to the strong correlation structure. This modified MetropolisHastings algorithm is described in the seminal work by Haario et al. (2001) and has become known as the adaptive Metropolis (AM) algorithm. It works like the regular MetropolisHastings sampler, but relies on the idea of allowing the proposal distribution to depend on the previous values of the chain.
The regular MetropolisHastings algorithm produces chains based on a proposal distribution that may depend on the current state of the chain. The proposal distribution is often chosen to be a multivariate normal distribution with fixed covariance. The main difference between the regular Metropolis sampler and AM is that AM updates the covariance matrix during the sampling, e.g. every 1000 iterations. Because the proposal distribution is thus tuned based on the entire sampling run, AM chains in fact lose the Markov property. Nonetheless, it can be shown that the algorithm retains the correct ergodic properties under very general assumptions, i.e., the covariance matrix stabilizes during the sampling process and AM chains properly simulate the target distribution (Haario et al. 2001; Vihola 2011). Although adaptive MCMC algorithms are a recent development, they have successfully been used by other authors in the astronomical community, for instance by Balan & Lahav (2009), Irwin et al. (2010) and Irwin et al. (2011).
Our MCMC calculations make extensive use of routines of PyAstronomy^{3}, a collection of Python routines providing a convenient interface for fitting and sampling algorithms implemented in the PyMC (Patil et al. 2010) and SciPy (Jones et al. 2001) packages. The AM sampler is implemented in the PyMC package, which is publicly available for download. We refer to the detailed online documentation^{4}.
We checked that AM does yield correct results for simulated data sets with parameters close to Qatar1, and found that this approach showed fast convergence and was efficient. To express our lack of more a priori knowledge regarding the Qatar1 system parameters, we assumed uninformative uniform prior probability distributions for all parameters, but we found that the parameters are well determined, so that the actual choice of the prior is unimportant for our results.
Fig. 3 OLT and PTST light curves (top panels). Superposed is the Mandel & Agol (2002) primary transit lightcurve model in a continuous line, considering the parameters obtained in Sect. 3.2. As indicated by the dates, the light curves evolve in time from bottom to top and from left to right. The residuals (bottom panels) have been calculated after subtracting the best transit fit parameters. 
Fig. 3 continued. 
Fig. 3 continued. 
In Fig. 3 we show our 26 obtained light curves, along with the residuals after removing the primary transit feature. Our final relative photometry is available in its entirety in machinereadable form in the electronic version of this paper. First columns contain BJD_{TDB}, second columns normalized flux, and third columns individual errors.
3.5. Correlated noise
Carter & Winn (2009, and references therein) studied how timecorrelated noise affects the estimation of the parameter of transiting systems. To quantify whether and to what extent our light curves are affected by red noise, we reproduced part of their analysis as follows: First, we produced residuals from our final fits by subtracting the primary transit model from each light curve. We individually divided each light curve into M bins of equal duration. Since our data are not always equally spaced, we calculated a mean value N of data points per bin. If the data are not affected by red noise, they should follow the expectation of independent random numbers, where σ_{1} is the sample variance of the unbinned data and σ_{N} is the sample variance (or rms) of the binned data, with the following expression: where is the mean value of the residuals per bin, and is the mean value of the means.
If correlated noise is present, then each value σ_{N} will differ by a factor β_{N} from their expectation. The parameter β, an estimation of the strength of correlated noise in the data, is found by averaging β_{N} over a range Δn corresponding to time scales that are judged to be most important. For data sets free of correlated noise, we expect β = 1. For primary transit observations, Δn is the duration of ingress or egress. For Qatar1, the time between first and second contact (or equivalently, the time between third and fourth contact) is ~15 min.
In Fig. 4 we show the results of our correlated noise analysis for nine of the longest light curves. Black lines represent the expected behavior in the absence of red noise, and red and green lines represent the variance of the binned data for OLT and PTST, respectively, as a function of bin size. As expected, the larger the bin size, the smaller the rms. For each light curve we calculated β considering Δn = 15 minutes. For the OLT and PTST primary transits, β lies between β = 0.78 and β = 1.33, with ⟨ β ⟩ = 1.038. Thus, there is no evidence for significant correlated noise in our light curves.
Finally, Pont et al. (2006) suggested to enlarge individual photometric errors by a factor β to account for systematic effects on the light curves. This would increase the parameter errors without changing the parameter estimates. Since our light curves do not present any strong evidence of correlated noise, we did not modify the individually derived errors.
Fig. 4 Qatar1 rms in parts per million (ppm) of the timebinned residuals as a function of bin size in logarithmic scale. Red and green lines correspond to OLT and PTST data respectively, and black lines show the expected behavior under the presence of uncorrelated noise. 
3.6. Effects of exposure times and transit coverage
The accuracy in determining of the midtransit time is affected by the number of data points during transit and the number of offtransit data points, even more so when a normalization is involved.
To study whether the midtransits are affected by the transit observation duration, we computed the timing residual magnitudes as a function of the number of data points per transit. If any systematic effect dominates the light curves, a larger timing offset for those transits that are sampled the least is expected. We show our results in Fig. 5. Since the light curves are normalized and the midtransits might be sensitive to normalizations, we also computed the timingresidual magnitudes as a function of the number of offtransit data points. We used the Pearson correlation coefficient r to quantify the correlation of the midtransit offsets. In both cases this was r ~ −0.11 for all data points, and r ~ 0.05, which rules out the 0th epoch, which was expected to be zero by construction. We found no significant correlation.
Fig. 5 Timingresidual magnitudes in minutes as a function of in and offtransit data point number (top) and offtransit data point number only (bottom). The outermost point close to zero is the primary transit, which was selected to be the 0th epoch. Midtransit errors are not plotted to avoid visual contamination. 
Kipping (2010) studied the effects of finite integration times on the determination of the orbital parameters. He showed that the time difference between the midtransit moment and the nearest light curve data point might cause a shift in the TTV signal, which is expected to be one half of the rate sampling. To test how significant this effect is over our light curves, we calculated the mean exposure times per observing night, which was about 80 s. Figure 6 shows the mean exposure times per night, versus the magnitude of the timing residuals. Most of the data points lying around the halfmean exposure time (~35 s) were identified to have the most accurate midtransit timings, with exposure times of about one minute, which makes it unlikely that they are affected by a sampling effect. The Pearson correlation coefficient of r = −0.15 again reveals no significant correlation.
Fig. 6 Timingresidual magnitudes in minutes for OLT (red circles) and PTST (green squares), as a function of mean exposure times. The horizontal dashed line indicates half of the mean exposure time. 
3.7. Results
In general terms, the parameters P and T_{0} are usually correlated with each other, but are uncorrelated with the remaining transit parameters. Therefore, these parameters can be determined quite accurately without interference of the remaining ones. We used a NelderMead simplex to approach the bestfit solution, which is provided as the starting values of the MCMC sampler. The AM algorithm then samples from the posterior distribution for the parameters to obtain error estimates. After 10^{8} iterations we discarded a suitable burnin (typically 10^{6} samples) and determined the combination of parameters resulting in the lowest deviance. We consider the lowest deviance as our global bestfit solution. The errors were derived from the 68% highest probability density or credibility intervals (1σ). Our results are summarized in Table 3. None of our fit parameters are consistent with those reported by Alsubai et al. (2011); a possible explanation for this inconsistency might be that these authors determined the system parameters using only four light curves, two of which were incomplete, obtained in less than two weeks. During the revision of this paper, Covino et al. (2013) presented highprecision radial velocity measurements from which the RossiterMcLauglin effect was observed. The authors also obtained five new photometric transit light curves, from which the orbital parameters of the system were improved. With respect to the orbital period, one of the most important parameters for determining TTVs, our bestfitted orbital period seems to be consistent within errors with the one found by Covino et al. (2013).
Bestfit (lowest deviance) parameters for our 26 transits of Qatar1 after 10^{8} MCMC samples.
OLT (top) and PTST (bottom) fitted midtransits with 1σ errors and the OC data points.
Fig. 7 OC diagram for Qatar1 in minutes (left axis) and days (right axis) considering the OLT (red circles) and PTST (green squares) data points, along with the timing residuals. Our initial bestfitting model is overplotted with a continuous black line, and the second one with a dashed black line. 
To obtain individual midtransit times T_{mid,i}, we considered the bestfit values of p, a, i, u_{1}, u_{2}, and P. We specified Gaussian priors on these parameters and refitted each one of the 26 individual transit light curves for the midtransit time T_{mid,i}. Our results are listed in Table 4 together with the derived 1σ errors on these transit times. To compute the timing deviations compared with a constant period we fitted the observed midtransit times T_{o,i} to the expression finding the ephemeris as bestfitting values. All errors are obtained from the 68.27% confidence level of the marginalized posterior distribution for the parameters. With these parameters we computed the OCvalues also listed in Table 4. We used the available simultaneous observations separated two months from each other as a diagnostics of our fitting procedure. Both midtransits (epochs –43 and 0) are consistent with each other within the errors.
4. Transit timing variation analysis
4.1. OC diagram
In the absence of any timing variations we expect no significant deviations of the derived OCvalues from zero. Testing the null hypothesis OC = 0 with a χ^{2}test, we found = 2.56 with 24 degrees of freedom. This high value led us to reject the null hypothesis that OC vanishes.
We then applied the LombScargle periodogram (Lomb 1976; Scargle 1982; Zechmeister & Kürster 2009) to search for any significant periodicity contained in the OC diagram. Figure 7 shows the data points used to perform the periodogram, along with a list in Table 4. The resulting periodogram reveals a first peak at ν_{TTV,1} = 0.00759 ± 0.00075 cycl P^{1} (corresponding to a period of 187 ± 17 days) and a second one at about half the frequency ν_{TTV,2} = 0.00367 ± 0.00059 cycl P^{1} (corresponding to 386 ± 54 days). For both periodic signals, new ephemeris were refitted using the linear trend and a sinusoidal variation in the form For the most significant frequency, the fitted amplitude and phase are A_{TTV,1} = 0.00052 ± 0.0002 days (equivalently, 0.75 ± 0.28 min) and φ_{TTV,1} = 0.04 ± 0.05. We then recalculated for the first frequency. Using an Ftest to check the significance of the fit improvement we obtained a pvalue of 0.02, indicating that the sinusoidal trend does indeed provide a better description than the constant at 2.8σ level. We can also compare the models using the Bayesian information criterion, BIC = χ^{2} + klnN, which penalizes the number k of model parameters given N = 26 data points. We obtained BICs of 68.1 and 57.2 for the constant and sinusoidal trend, respectively. Since the BIC is no more than a criterion for model selection among a set of models, the method one more time favors the sinusoidal trend.
ETDfitted instant of minima with 1σ errors and the OC data points.
We finally estimated the falsealarm probability (FAP) of the TTV signal, using a bootstrap resampling method by randomly permuting the midtransit values (5 × 10^{5} times) along with their individual errors, fixing the observing epochs and calculating the LombScargle periodogram afterward. We estimated the FAP as the frequency with which the highest power in the scrambled periodogram exceeds the maximal power in the original periodogram. In this fashion we estimated an FAP of 0.05% for our observed TTV signal, consistent with our previous estimates.
To further check whether the addition of the sinusoidal variation provides an explanation for the timing offsets, we made use of the midtransits of Qatar1 available in the Exoplanet Transit Database^{5} (ETD; Poddaný et al. 2010) to reinforce our results. The ETD data are very heterogeneous, we therefore converted the best 26 midtransits from HJD_{UTC} to BJD_{TDB}. Each ETD primary transits has a data quality indicator ranging from 1 (small scatter and good timesampling) to 4 (the opposite). For our analysis we only selected the transits with quality flags 1 (5 primary transits) and 2 (21 primary transits). These transits span a total of ~500 epochs. We produced the OC diagram (Table 5) using the following fitted ephemeris: The primary transit selected for the 0th epoch is the most accurate one of all available light curves. We produced a periodogram using the resulting timing offsets and found one peak at ν_{TTV,ETD} = 0.00784 ± 0.0011 cycl P^{1} (corresponding to a period of 181 ± 22 days) with an amplitude of A_{TTV,ETD} = 0.0009 ± 0.0005 days (equivalently, 1.29 ± 0.72 min) and a phase of φ_{TTV,ETD} = 0.02 ± 0.03. Within the errors, the frequency, amplitude, and phase are consistent with the ones found using OLT and PTST data. Figure 8 shows two periodograms, the one on top produced using our data, and the one on bottom produced with ETD data.
Fig. 8 Zechmeister & Kürster (2009) periodograms generated from the timing residuals of Qatar1, showing a peak at ν_{TTV,1} = 0.00759 ± 0.00075 cycl P^{1} (our data, top panel) and ν_{TTV,ETD} = 0.0078 ± 0.0011 cycl P^{1} (ETD data, bottom panel). Vertical lines indicate the 1σ error on ν_{TTV,1}. The parameters derived from our periodogram are maximum power =6.1, and FAP for the maximumpower peak =0.19%. 
4.2. Error analysis
To test the reliability of our error estimates for each midtransit time and to check the influence of the error magnitude on the estimated peak frequencies, we iteratively constructed new LombScargle periodograms randomly increasing the individual midtransit errors mostly by a factor of two.
At each iteration, we first added the absolute value of a random number that was drawn from a normal distribution (μ = 0, σ = 0.0004 days) to each midtransit error, and calculated the leading frequency afterward. After 5 × 10^{5} of such iterations, we found that the only effect on the error increment is the permutation of the leading frequency from ν_{TTV,1} to ν_{TTV,2}. This permutation occurred only ~9% of the times.
In Fig. 9 we show the resulting histogram of the calculated leading frequencies. They all fall inside the 1σ errors of ν_{TTV,1} and ν_{TTV,2}, estimated from our original fits. Thus, the derived midtransit errors do not significantly affect the outcome in terms of dominating frequencies.
Fig. 9 Histogram of the leading frequencies obtained after incrementing the individual midtransit errors, in units of 10^{3}. Vertical continuous lines indicate our best two leading frequencies ν_{TTV,1} and ν_{TTV,2}, along with 1σ errors (dashed vertical lines). 
Fig. 10 Simulated TTV standard deviation in minutes of Qatar1b as a function of the outer perturber period, for two values of eccentricity and three values of perturber masses. The gray area on the left correspond to perturber orbits that would make the Qatar1b orbit unstable. The mean motion resonances 2:1, 3:1, and 5:2 are indicated with vertical dashed lines. 
Fig. 11 OC diagram for OLT (red circles) and PTST (green squares) data points plus Nbody solutions for two different dynamical scenarios (black lines) and our best initial sinusoidal fit (blue dashed lines) artificially shifted in phase for a better comparison. 
Three possible solutions for our two main TTV signals for three different dynamical scenarios.
5. Interpretation of the observed transit timing variations
For the purposes of this section we assumed that the reported TTVs are real and explored possible physical scenarios that would explain the observed variations. We considered both the ~190 and ~380 day periods, although the ETD data set is consistent only with the former period. As is well known, a suitably placed third body in the Qatar1 system can lead to perturbations in the orbit of Qatar1b, which manifest themselves as TTVs. Since there are no complete analytical solutions to the threebody problem, we made use of a numerical integration scheme. Because a planetary orbit is determined by a large number of parameters (i.e., the eccentricity, the longitude of the nodes, the orbital inclination, among many others), it is a real challenge to estimate true orbital parameters by fitting TTV models to groundbased midtransit shifts only. We therefore started with the simplest approach and assumed coplanar orbits using the inclination derived from the light curve fitting, along with the bestfit orbital period, and the transiting planet mass and eccentricity obtained by Alsubai et al. (2011) as fixed. Given trial mass, orbital period, eccentricity, longitude of periastron and ascending node, and the time of periastron passage, the Nbody code calculates the time of occurrence of the primary transits taking into account the interaction between the planets. We specifically considered two scenarios, a putative lowmass planet in a resonant orbit relative close to Qatar1b, and a putative highmass planet or brown dwarf in an 190day orbit around Qatar1.
5.1. Weak perturber in resonance with Qatar1b
Two planetary bodies in orbital resonance with each other will experience longterm changes in their orbital parameters. On one hand, since only a relative short stretch of data (less than two years) is available, expecting a unique solution is too ambitious. On the other hand, most of the perturber’s dynamical setups will translate into negligible TTVs.
Taking into account these difficulties, we calculated the modeled TTV scatter for a twoplanetary system, to study in advance which would be the parameter space giving rise to TTVs similar to the scatter of our data. We first considered both bodies in circular orbits and then the perturber in an eccentric orbit. The mass of the perturber was systematically changed between 1, 8, and 15 M_{⊕}. Figure 10 shows our results. The vertical dashed lines indicate the 2:1, 3:1, and 5:2 resonances. In these regions, perturbers of the order of several Earth masses would produce TTVs similar to the observed ones.
We thus investigated some possible solutions that would satisfy the first two TTV periods assuming the specific resonant orbits indicated in Fig. 10. As an example we show in Fig. 11 and Table 6 three arbitrary cases; similar solutions can be derived with other choices of resonant orbits. To produce Fig. 11 we ran our Nbody simulation with specified initial conditions (cf., Table 6) and compared the results with the observed OCvalues and our (original) sinusoidal fit (cf., Fig. 11). As is clear from Fig. 11, all these solutions have χ^{2} values similar to our first fitting attempt. Thus, it is futile to produce a proper fitting procedure.
Fig. 12 TTV amplitude as a function of the mass of the perturber for different eccentricity values. The black continuous line shows the bestfitted amplitude to our OC diagram, along with 1σ errors (dashed lines). 
5.2. Massive perturber in a 190day orbit
As an alternative to a lowermass perturber in a resonant orbit, we also considered a more massive perturber in a nonresonant 190day orbit. In this case the TTV amplitude strongly depends on the mass and the eccentricity of the assumed perturber. Again considering the simplest case of coplanar orbits and both longitude of periastron and longitude of the ascending node fixed to zero, we computed the amplitude of the TTV signal as a function of the perturber mass and its orbital eccentricity using our Nbody code, and we compared the predicted TTV amplitudes with those observed. Our results are shown in Fig. 12, where we plot the expected TTV amplitude as a function of the perturber mass for various eccentricities in the range between zero and 0.8. For low eccentricity values we require perturber masses of more than 80 Jupiter masses, i.e., we would require a lowmass star. For eccentricities above 0.6 a brown dwarf could explain the observed TTV amplitudes. It is obvious that such an object would lead to RV variations in the host of Qatar1, which could in principle be observed. Since for a very eccentric orbit the
RV variations are concentrated around the periastron passage and since only ten days of RV measurements are available for Qatar1 so far, any longterm RV variations of Qatar1 are unknown so far and could confirm or reject the presence of a massive perturber.
6. Conclusions
Our analysis of the midtiming residuals of Qatar1 taken during almost two years indicate that the orbital period of the exoplanet is not constant. The observed longterm timing variations are highly significant from a statistical point of view and can be explained by very different physical scenarios. RV monitoring of Qatar1 will provide an upper limit to the mass of a possible perturber and continued timing observations of Qatar1 are required to better delineate the solution space for the possible perturber geometries.
http://www.hs.unihamburg.de/DE/Ins/Per/Czesla/ PyA/PyA/index.html
Acknowledgments
C. von Essen acknowledges funding by the DFG in the framework of RTG 1351, and P. Hauschildt, S. Witte and H. M. Müller for discussions about the PHOENIX code and limbdarkening effects.
References
 Adams, E. R., LópezMorales, M., Elliot, J. L., Seager, S., & Osip, D. J. 2010, ApJ, 714, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Agol, E., Steffen, J., Sari, R., & Clarkson, W. 2005, MNRAS, 359, 567 [NASA ADS] [CrossRef] [Google Scholar]
 Alsubai, K. A., Parley, N. R., Bramich, D. M., et al. 2011, MNRAS, 417, 709 [NASA ADS] [CrossRef] [Google Scholar]
 Balan, S. T., & Lahav, O. 2009, MNRAS, 394, 1936 [NASA ADS] [CrossRef] [Google Scholar]
 Ballard, S., Fabrycky, D., Fressin, F., et al. 2011, ApJ, 743, 200 [NASA ADS] [CrossRef] [Google Scholar]
 Carter, J. A., & Winn, J. N. 2009, ApJ, 704, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Claret, A., & Hauschildt, P. H. 2003, A&A, 412, 241 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cochran, W. D., Fabrycky, D. C., Torres, G., et al. 2011, ApJS, 197, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Covino, E., Esposito, M., Barbieri, M., et al. 2013, A&A, 554, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Díaz, R. F., Rojo, P., Melita, M., et al. 2008, ApJ, 682, L49 [CrossRef] [Google Scholar]
 Eastman, J., Siverd, R., & Gaudi, B. S. 2010, PASP, 122, 935 [NASA ADS] [CrossRef] [Google Scholar]
 Fabrycky, D. C., Ford, E. B., Steffen, J. H., et al. 2012, ApJ, 750, 114 [NASA ADS] [CrossRef] [Google Scholar]
 Fukui, A., Narita, N., Tristram, P. J., et al. 2011, PASJ, 63, 287 [NASA ADS] [Google Scholar]
 Fulton, B. J., Shporer, A., Winn, J. N., et al. 2011, AJ, 142, 84 [NASA ADS] [CrossRef] [Google Scholar]
 Haario, H., Saksman, E., & Tamminen, J. 2001, Bernoulli, 7, 223 [CrossRef] [MathSciNet] [Google Scholar]
 Hauschildt, P. H., & Baron, E. 1999, J. Comput. Appl. Math., 109, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Holman, M. J., & Murray, N. W. 2005, Science, 307, 1288 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Holman, M. J., Fabrycky, D. C., Ragozzine, D., et al. 2010, Science, 330, 51 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Irwin, J., Buchhave, L., Berta, Z. K., et al. 2010, ApJ, 718, 1353 [NASA ADS] [CrossRef] [Google Scholar]
 Irwin, J. M., Quinn, S. N., Berta, Z. K., et al. 2011, ApJ, 742, 123 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientific tools for Python, http://www.scipy.org [Google Scholar]
 Kipping, D. M. 2010, MNRAS, 408, 1758 [NASA ADS] [CrossRef] [Google Scholar]
 Lissauer, J. J., Fabrycky, D. C., Ford, E. B., et al. 2011, Nature, 470, 53 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Lomb, N. R. 1976, Ap&SS, 39, 447 [NASA ADS] [CrossRef] [Google Scholar]
 Maciejewski, G., Dimitrov, D., Neuhäuser, R., et al. 2010, MNRAS, 407, 2625 [NASA ADS] [CrossRef] [Google Scholar]
 Maciejewski, G., Dimitrov, D., Neuhäuser, R., et al. 2011, MNRAS, 411, 1204 [NASA ADS] [CrossRef] [Google Scholar]
 Mandel, K., & Agol, E. 2002, ApJ, 580, L171 [NASA ADS] [CrossRef] [Google Scholar]
 Nascimbeni, V., Piotto, G., Bedin, L. R., et al. 2011, A&A, 532, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pál, A., Sárneczky, K., Szabó, G. M., et al. 2011, MNRAS, 413, L43 [NASA ADS] [Google Scholar]
 Patil, A., Huard, D., & Fonnesbeck, C. J. 2010, J. Stat. Software, 35, 1 [Google Scholar]
 Poddaný, S., Brát, L., & Pejcha, O. 2010, New A, 15, 297 [Google Scholar]
 Pont, F., Zucker, S., & Queloz, D. 2006, MNRAS, 373, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Scargle, J. D. 1982, ApJ, 263, 835 [NASA ADS] [CrossRef] [Google Scholar]
 Steffen, J. H., Ragozzine, D., Fabrycky, D. C., et al. 2012, Proc. Nat. Acad. Sci., 109, 7982 [NASA ADS] [CrossRef] [Google Scholar]
 Vihola, M. 2011, Stochastic Processes and their Applications, in press [arXiv:0903.4061] [Google Scholar]
 Witte, S., Helling, C., & Hauschildt, P. H. 2009, A&A, 506, 1367 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
Bestfit limbdarkening coefficients (LDCs) for OLT and PTST, along with the 1σ errors.
Bestfit (lowest deviance) parameters for our 26 transits of Qatar1 after 10^{8} MCMC samples.
OLT (top) and PTST (bottom) fitted midtransits with 1σ errors and the OC data points.
Three possible solutions for our two main TTV signals for three different dynamical scenarios.
All Figures
Fig. 1 OskarLühning telescope (dark background) and Planet Transit Search Telescope (light background) field of views. Qatar1 is indicated with a large red circle centered on the OLT field of view, along with the comparison stars as green squares. The small blue circle upward and to right of Qatar1 indicates the reference star used to produce the differential photometry. 

In the text 
Fig. 2 Qatar1 normalized intensities considering the OLT and PTST filter transmission functions. Data points on the right of the vertical dashed line were not considered in the fitting procedure. 

In the text 
Fig. 3 OLT and PTST light curves (top panels). Superposed is the Mandel & Agol (2002) primary transit lightcurve model in a continuous line, considering the parameters obtained in Sect. 3.2. As indicated by the dates, the light curves evolve in time from bottom to top and from left to right. The residuals (bottom panels) have been calculated after subtracting the best transit fit parameters. 

In the text 
Fig. 3 continued. 

In the text 
Fig. 3 continued. 

In the text 
Fig. 4 Qatar1 rms in parts per million (ppm) of the timebinned residuals as a function of bin size in logarithmic scale. Red and green lines correspond to OLT and PTST data respectively, and black lines show the expected behavior under the presence of uncorrelated noise. 

In the text 
Fig. 5 Timingresidual magnitudes in minutes as a function of in and offtransit data point number (top) and offtransit data point number only (bottom). The outermost point close to zero is the primary transit, which was selected to be the 0th epoch. Midtransit errors are not plotted to avoid visual contamination. 

In the text 
Fig. 6 Timingresidual magnitudes in minutes for OLT (red circles) and PTST (green squares), as a function of mean exposure times. The horizontal dashed line indicates half of the mean exposure time. 

In the text 
Fig. 7 OC diagram for Qatar1 in minutes (left axis) and days (right axis) considering the OLT (red circles) and PTST (green squares) data points, along with the timing residuals. Our initial bestfitting model is overplotted with a continuous black line, and the second one with a dashed black line. 

In the text 
Fig. 8 Zechmeister & Kürster (2009) periodograms generated from the timing residuals of Qatar1, showing a peak at ν_{TTV,1} = 0.00759 ± 0.00075 cycl P^{1} (our data, top panel) and ν_{TTV,ETD} = 0.0078 ± 0.0011 cycl P^{1} (ETD data, bottom panel). Vertical lines indicate the 1σ error on ν_{TTV,1}. The parameters derived from our periodogram are maximum power =6.1, and FAP for the maximumpower peak =0.19%. 

In the text 
Fig. 9 Histogram of the leading frequencies obtained after incrementing the individual midtransit errors, in units of 10^{3}. Vertical continuous lines indicate our best two leading frequencies ν_{TTV,1} and ν_{TTV,2}, along with 1σ errors (dashed vertical lines). 

In the text 
Fig. 10 Simulated TTV standard deviation in minutes of Qatar1b as a function of the outer perturber period, for two values of eccentricity and three values of perturber masses. The gray area on the left correspond to perturber orbits that would make the Qatar1b orbit unstable. The mean motion resonances 2:1, 3:1, and 5:2 are indicated with vertical dashed lines. 

In the text 
Fig. 11 OC diagram for OLT (red circles) and PTST (green squares) data points plus Nbody solutions for two different dynamical scenarios (black lines) and our best initial sinusoidal fit (blue dashed lines) artificially shifted in phase for a better comparison. 

In the text 
Fig. 12 TTV amplitude as a function of the mass of the perturber for different eccentricity values. The black continuous line shows the bestfitted amplitude to our OC diagram, along with 1σ errors (dashed lines). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.