HAT-P-11b: A SUPER-NEPTUNE PLANET TRANSITING A BRIGHT K STAR IN THE KEPLER FIELD^*

The 84 × 84 region around GSC 03561–02092, a field internally labeled as "G155," was observed on a nightly basis in two seasons, whenever weather conditions permitted. First, during the Fall of 2004, we acquired 1213 frames with the HAT-6 instrument located at FLWO, and 4091 frames with HAT-9 located at SMA, Mauna Kea. We revisited the field in 2005, and acquired an additional 6166 frames with HAT-9. Altogether we gathered 11470 5 minute exposures at a 5.5 minute cadence. This unusually rich data set was motivated by the overlap of G155 with the field of view (FOV) of the future Kepler mission.

The calibration of the HATNet frames was done utilizing standard procedures based on IRAF.¹⁴ The calibrated frames were then subjected to star detection and astrometry, as described by Pál & Bakos (2006). Aperture photometry using three apertures was performed on each image at the fixed positions of the stellar centroids, as derived from the Two Micron All Sky Survey (2MASS) catalog (Skrutskie et al. 2006) and the individual astrometric solutions relating the pixel coordinates to the world coordinate (ICRS) system of 2MASS. We extracted photometry for all 125,000 stars down to I ≲ 14 that fell in the field. The raw instrumental magnitudes m_r of each individual frame were transformed to a reference frame by fourth-order polynomials in X and Y, and first order in color. The fitted magnitudes m_f yielded by the above smooth fit are used for generating the time versus m_f light curves.

These light curves have a noise characteristic that is sometimes referred to as pink noise (Pont et al. 2006), because it is a combination of Gaussian white noise (due to e.g., photon noise) and a correlated red noise (due to e.g., spatial drift of stars and uncorrected flatfield effects). The red noise is often referred to as trends or systematic variations. Trends in the light curves have an adverse effect on detecting shallow transiting signals: they mask the real signal, and when phase-folded with a trial period they can also mimic transit signals. Furthermore, for photometric follow-up data, the trends can lead to large systematic errors in the transit parameter determination, such as in the impact parameter or the depth of the transit. Thus, proper treatment of systematic variations is crucial for both discovery and accurate follow-up characterization, especially for the shallow (4.2 mmag) transit of a hot Neptune presented in this paper. The two basic methods we have employed are the external parameter decorrelation (EPD) technique, briefly described earlier in Bakos et al. (2007), and the trend filtering algorithm (TFA; Kovács et al. 2005). The technical details and some new definitions are given in the Appendix.

The HATNet light curves were decorrelated against trends using the EPD technique (in constant EPD mode), and subsequently by a simple global TFA (without reconstruction, and separately from the EPD). For the brightest stars in the field we achieved a photometric precision of 2.9 mmag at 5.5 minute cadence. The light curves were searched for periodic box-like signals using the box least-squares transit search method (BLS; see Kovács et al. 2002). The BLS frequency spectrum of GSC 03561–02092 (also known as 2MASS 19505021+4804508; α = 19^h50^m5021, δ = +48^d04^m508; J2000) showed a number of significant peaks, the highest one at 0.03148 [c/d] (31.76 day period), and the second one at ∼0.2046 [c/d] (4.8878 days). Fortunately, we inspected the second peak, and found the corresponding transit (Figure 1) with a depth of 4.2 mmag, worthy of follow-up. The dip had a relative duration (first to last contact) of q ≈ 0.0206, equivalent to a total duration of Pq ≈ 2.41 hr. The signal was confirmed with subsequent reconstructive TFA using a trapeze-shaped model function, being the most significant periodicity in the data (the 31.76 day signal was probably due to a left-over systematic). The rms of the TFA processed light curve was ∼3.2 mmag, making the detection of the 4.2 mmag transit dip rather challenging, but still robust with S/N = 11.4. To our best knowledge, this is the shallowest transit signal found by ground-based transit searches that belongs to a confirmed planet. HD 149026, with an even shallower transit signal (Winn et al. 2008) was first detected via spectroscopy (Sato et al. 2005). The second shallowest planetary transit found by the transit search method is HAT-P-2b (Bakos et al. 2007). We note that on occasion HATNet has found even shallower transiting signals (down to 2 mmag depth) where the transit was real, but was due to a blended system.

Zoom In Zoom Out Reset image size

Following the procedures described in Latham et al. (2009), we used the CfA Digital Speedometers (DS; Latham 1992), mounted on the 1.5 m Wyeth Reflector at the Oak Ridge Observatory in Harvard, Massachusetts, and on the 1.5 m Tillinghast Reflector at FLWO, to obtain low S/N ratio high-resolution spectra. Altogether five spectra were obtained between 2001 December 17 and 2007 June 5 (the first observations were obtained for a different project to survey dwarfs in the solar neighborhood that might be suitable targets for SETI). The signal-to-noise ratios ranged from 10 to 18 per spectral resolution element of 8.5 km s⁻¹.

Reconnaissance spectroscopy is an important step in weeding out astrophysical false positive systems that mimic planetary transits. Stellar parameter determination can distinguish between dwarf and giant stars via the measurement of surface gravity log g_⋆, and thus eliminate systems where the light curve of an eclipsing binary is blended with a giant star, or where the transit signal seen in our data may not be real, as it is physically not feasible to orbit around (and outside) the giant star with such a short period. In addition, large RV variations (of the order of several km s⁻¹) are indicative of orbital motion due to stellar companions rather than planets orbiting a star. Fine analysis of the spectra can also reveal stellar triple systems. Finally, the rapid rotation of a host star, as may be indicated by the rotational broadening of the spectra, is often correlated with a stellar companion that is massive enough to synchronize the rotation with the orbital motion.

HAT-P-11 survived all these steps, and the RV measurements showed an rms residual of 0.29 km s⁻¹, consistent with no detectable RV variation. Initial atmospheric parameters for the star, including the effective temperature T_eff⋆ = 4750  ±   125 K, surface gravity log g_⋆ = 4.5  ±   0.25 (cgs), and projected rotational velocity vsin i = 0.0 km s⁻¹, were derived as described by Latham et al. (2009). The mean line-of-sight velocity of the star was Γ = −63.56  ±   0.29 km s⁻¹ on an absolute scale.

Photometric follow-up observations are important (1) to perform independent confirmation of the initial detection, (2) to allow for accurate characterization of the system (see also Section 5), (3) to help eliminate blend scenarios (Section 4), and, if multiple events are observed, (4) to possibly average out extra variations of the light curve due to spots. Since the transit of HAT-P-11 was originally detected with the 11 cm diameter HATNet telescopes (although when phase-folding many events), it is certainly feasible to confirm a single transit with a 1 m class telescope. However, acquiring a high-quality observation of the 4.2 mmag transit dip in the light of a star as bright as V = 9.59 is quite challenging. An accurate and highly precise light curve is essential for eliminating blend scenarios, and then for determining the physical parameters of the transiting planet–star system. In general, the shallower the transit, the wider the range of possible blends (i.e., more than two body systems) that can mimic the observed transit. Some of these can only be distinguished by subtle effects, such as the duration of ingress/egress, and accurate depth and shape of the transit.

For these reasons, we launched an extensive photometry follow-up campaign, and attempted transit observations of HAT-P-11b altogether 12 times in 2007 and 2008, leading to ∼10 successful observations of partial or full transits (see Table 1). We primarily used the FLWO 1.2 m telescope and the KeplerCam CCD in Sloan z band, with exposure times of ∼10 s and read-out time of 12 s. The last observation on 2008 November 20/21 MST was taken through Sloan r band to complement our blend analysis (Section 4). We also observed transits 3 times using the 0.6 m Schmidt telescope of Konkoly Observatory at the Piszkéstető Mountain Station. These observations were taken through I_C band, and the only useful data set proved to be the one from 2008 September 3/4 CET.

Data were reduced in a similar manner as for the HATNet data (Section 2.1), and as described in Bakos et al. (2009). Following bias and flat calibration, we derived an initial first-order astrometrical transformation between the ∼750 brightest stars and the 2MASS catalog. The position of HAT-P-11 in the 2MASS catalog was corrected for each observed epoch before deriving the astrometry due to its moderately high proper motion (263.32  ±   1.31 mas yr⁻¹; Perryman et al. 1997). Photometry was carried out for all stars in the field using three apertures that were adjusted each night to match the observing conditions (sky background, profile width). For the nights with good seeing, one aperture was kept small to avoid the faint neighbor 2MASS 19505049+4805017 currently at approx 6'' distance (see Section 4). Instrumental magnitudes were transformed to a photometric reference frame (selected to be at low air mass, low sky background, etc.) using a first-order polynomial in the X, Y pixel coordinates and the 2MASS J − K color of ∼600 stars. The transformation was iteratively determined by weighting with individual Poisson noise errors of the stars, and rejecting 3σ outliers, plus eliminating the main target and its faint neighbor from the fit. The smooth fit was repeated based on the ∼180 best stars, using the median of the individual light-curve magnitude values as a new reference system, and weighting the fit with the rms of the light curves (i.e., substituting the former Poisson errors).

The FLWO 1.2 m KeplerCam observations usually result in high-quality photometry of ∼1% deep planetary transits (see e.g., Latham et al. 2009; Bakos et al. 2009), because of the large FOV (= 23' × 23') with many potential comparison stars, fine pixel resolution (0336 pixel⁻¹), good quality and high quantum efficiency CCD (monolithic, 4k × 4k Fairchild 486 chip), good sky conditions from FLWO, the fast readout through four channels, and careful data processing. The performance on this very shallow transit of a fairly bright star, however, turned out to be slightly sub-optimal with some residual trends, necessitating diverse application of the EPD and TFA methods. Since these were part of our global modeling of the data, including the HATNet and RV data, they are detailed later in Section 5 and in the Appendix. The follow-up light curves after applying simultaneous EPD–TFA with per night coefficients for EPD and global coefficients for TFA are displayed in Figure 2; the photometric data is provided in Table 2. The model function (denoted as $m_0(\vec{p},t_i)$ in Equation (A1) in the Appendix) in the simultaneous EPD–TFA fit is an analytical transit model from Mandel & Agol (2002).

Zoom In Zoom Out Reset image size

We started observations of HAT-P-11 on 2007 August 22 with the HIRES instrument (Vogt et al. 1994) on the Keck-I telescope located on Mauna Kea, Hawaii. It was soon realized that the target is unusually complicated because of the small RV amplitude with respect to the moderate velocity jitter due to the active K dwarf star, and also because of long-term trend(s) present. Thus, HAT-P-11 has been extensively observed over the past two years, and we have gathered altogether 50 spectra and three template observations (see Table 3). This is 5–10 times the number of spectra collected for a typical HATNet transit candidate.

The width of the spectrometer slit used on HIRES was 086, resulting in a resolving power of λ/Δλ ≈ 55, 000, with a wavelength coverage of ∼3800–8000 Å. The iodine gas absorption cell was used to superimpose a dense forest of I₂ lines on the stellar spectrum and establish an accurate wavelength fiducial (see Marcy & Butler 1992). Relative RVs in the solar system barycentric frame were derived as described by Butler et al. (1996), incorporating full modeling of the spatial and temporal variations of the instrumental profile. The final RV data and their errors as a function of barycentric Julian date are listed in Table 3. It is reassuring that a simple Fourier analysis of the RV data without prior assumption on any periodic event yields a primary period of 4.888  ±   0.019 days, confirming the same periodic phenomenon that is present in the independent data set consisting of the discovery and follow-up photometry.¹⁵

Based on the numerous observations of a transit-like event in the photometry, our initial physical model was thus a single body orbiting and transiting the star, causing both the RV variations and the transits. The RV variations induced on the parent star can be characterized by six parameters: the period P, the center of transit T_c (i.e., the phase of the RV curve), the semi-amplitude K, the RV zero-point γ, and the Lagrangian orbital elements (k, h) = e × (cos ω, sin ω). An orbital fit to the RV data without any constraints from the photometry yielded the following values: P = 4.8896  ±   0.0017 days, T_c = 2454552.04  ±   0.16 (BJD), K = 13.1  ±   2.1 m s^-1, γ = 0.9  ±   0.8 m s^-1 (arbitrary scale), k = 0.28  ±   0.12, and h = 0.09  ±   0.09.

Later (Section 5) we present a global modeling of the data, where joint analysis of the photometry and RV data is performed. This yields precise P and T_c parameters that are almost entirely constrained by the photometric data (notably the sharp ingress/egress features), with virtually no coupling from the RV data. The knowledge of period and phase from photometry is a very tight constraint in the analysis of the RV data, yielding significant improvement on determination of the other orbital parameters. It is thus well justified to analyze the RV data separately by fixing the ephemeris to that determined by the photometry, and assuming that there is an eccentric orbital motion present. In addition, we can also investigate whether our initial model correctly described the physics of the system, and whether there are other signals present in the RV data.

Along these lines, our refined model for the RV data was an eccentric orbital motion with P and T_c fixed to those values found from the global modeling (and primarily constrained by the photometry), plus a sinusoidal motion of the form A₂sin(f₂t) + B₂sin(f₂t) (equivalent to A'₂sin(f₂t +  ϕ), but linear in the fitted A₂ and B₂ parameters). We searched the f₂ domain by fitting γ, K, k, h (eccentric orbit parameters) and A₂, B₂ at each f₂, and noting the χ². This "frequency scan" located a number of significant peaks with small inverse χ² (Figure 3, upper panel), the most significant being f₂ = 0.9602  ±   0.0009. We suspected that the peaks around 1 [c/d] might be due to a long-term trend in the RV, sampled with a daily and lunar-cycle periodicity. To test this effect, we generated a mock signal by co-adding an analytic eccentric orbital motion, a long-term RV drift, and Gaussian noise, and sampled these at the exact times of the 50 Keck observations. We ran the same fitting procedure (Keplerian plus sinusoidal motion) in the above-described frequency scan mode on the mock data. Indeed, the drift appeared as a number of significant frequencies around 1 [c/d] (Figure 3(b)), their exact location depending on the amplitude of the drift and the random seed used to generate the noise. Then, we repeated the above analysis including a long-term drift in the fitted model function while performing the frequency scanning procedure (i.e., fitting at each frequency value that is being stepped on a grid).

Zoom In Zoom Out Reset image size

The results are plotted in panel(c) of Figure 3, on the same scale as used in panels (a) and (b). It can be clearly seen that the magnitude of the residuals is definitely smaller and the strong structures disappeared. Using the mock data set, the results are almost the same (Figure 3(d)), confirming our assumptions that a long-term drift in the data causes aliases around 1 [c/d]. This is in line with the experience based on long-term monitoring of stars with Keck/HIRES: in nearly every case where there is either a planet or a long-term trend, a spurious spike appears in the periodogram at close to one day.

Note, however, that in spite of the clearing of peaks around 1 [c/d], the 0.9602 [c/d] peak remained present even after the simultaneous fit of a Keplerian orbit and a drift to the Keck data.

It is also interesting that the χ² of the fit with Keplerian plus f₂ sinusoid component is somewhat smaller than the fit with Keplerian plus long-term trend. The f₂ = 0.9602 [c/d] periodicity may be a real physical effect, or can be an alias of a real physical effect with a different period (rotation period, activity of star), or an alias of a systematic effect (lunar cycle, instrumental effect). Altogether we have three components to consider in the final model of the RV variations: (1) Keplerian orbit, (2) long-term drift, (3) f₂ = 0.9602 [c/d] periodicity.

The model of a Keplerian orbit and a long-term drift ((1)+(2)) has simple underlying physics by assuming an inner and an outer planets. The results of the orbital fit for this basic model are exhibited in Figure 4. The data folded with the P = 4.8878162 day period after removal of the best-fit drift is plotted in the top panel. Superimposed is the Keplerian orbit that is clearly eccentric. The orbital elements for this fit were: K = 11.6  ±   1.2 m s⁻¹, k = ecos ω = 0.201  ±   0.049, h = esin ω = 0.051  ±   0.092, γ = −0.4  ±   0.8 m s⁻¹, and the best-fit value for the linear drift is $G_1 = 0.0297 \,\,{\pm}\,\, 0.0050\,\rm m\,s^{-1}\,day^{-1}$. Note that the K amplitude of the orbit is only 11.6  ±   1.2 m s⁻¹. In order to have a reduced χ² value of 1.0, a jitter of 5.01 m s⁻¹ has to be added to the formal errors.

Zoom In Zoom Out Reset image size

We also derived an orbital fit with a full model of all three components (Keplerian, drift, short-period sinusoidal) so as to check the effect of f₂ on the orbital parameters. We found that the change in the orbital elements is insignificant, with a small decrease in their respective error bars: K = 10.91  ±   0.96 m s⁻¹, k = ecos ω = 0.225  ±   0.038, h = esin ω = 0.069  ±   0.072, drift $G_1 = 0.0176 \,\,{\pm}\,\, 0.0056\,\rm m\,s^{-1}\,day^{-1}$, and the amplitude of the f₂ sinusoid is 5.7  ±   1.3 m s⁻¹. The jitter value changed to 3.89 m s⁻¹.

We also investigated models with nonlinear drift, as characterized by G₂ quadratic and G₃ cubic terms, and found that these are insignificant based on the present data. We also checked for correlations between the RV residuals from the best fit, and the spectral bisector spans (BSs), and the S activity index, but found that these correlations were insignificant.

Given the fact that the origin of the f₂ periodicity (component 3) is unknown, and it is suspected to be an alias that may diminish by taking more data, plus the key orbital elements do not change significantly by taking it into account; in the rest of this paper, we adopted the simpler model of a Keplerian orbit (with $K=11.6\,\,{\pm}\,\, 1.2\,\rm m\,s^{-1}$) and a long-term drift without the short-term periodicity.

We employed the Spectroscopy Made Easy (SME) package of Valenti & Piskunov (1996) along with the atomic-line database of Valenti & Fischer (2005) to derive an initial value for the stellar atmospheric parameters. In this analysis, we used the iodine-free template spectrum obtained by the HIRES instrument on Keck I. The SME analysis of stellar spectra resulted in a stellar surface gravity log g_⋆ = 4.7  ±   0.1 (CGS), metallicity [Fe/H] = 0.32  ±   0.06 dex, effective temperature T_eff⋆ = 4850  ±   50 K, and the projected rotational velocity vsin i = 0.5  ±   0.5 km s⁻¹. However, the vsin i value depends inversely on the assumed value of macro-turbulence v_mac in the star. The value of v_mac assumed in the SME derivation was 2.57 km s⁻¹, but for stars of type K0 or later, v_mac may be as small as 1.3 km s⁻¹ (e.g., Gray 1988). With that value of v_mac, vsin i could be as high as 2.7 km s⁻¹. Because of the indeterminacy of v_mac for HAT-P-11, we can conclude solely that HAT-P-11 is a relatively slow rotator and assign it a projected rotational velocity vsin i = 1.5  ±   1.5 km s⁻¹, which encompasses both extremes of vsin i noted here.

At this stage, we could use the effective temperature as a color indicator and the surface gravity as a luminosity indicator, and determine the stellar parameters based on these two constraints using a set of isochrones. However, it has been shown (Sozzetti et al. 2007) that log g_⋆ has a subtle effect on the spectral line shapes, and is usually not the best luminosity indicator. For planetary transits, the a/R_⋆ normalized semimajor axis and related ρ_⋆ mean stellar density typically impose a stronger constraint on possible stellar models (Sozzetti et al. 2007). (The use of a/R_⋆ is discussed later in Section 3.2.) However, for HAT-P-11 there is an even better luminosity indicator, since it is a bright and nearby star with known parallax of small relative error. HAT-P-11 enters the Hipparcos catalog as HIP 97657 with reported parallax of 27.50  ±   0.96 mas (Perryman et al. 1997), equivalent to a distance of 36.4  ±   1.3 pc and distance modulus of Δ = 2.80  ±   0.08. Combination of the distance information with the apparent brightnesses in various photometric bands yields the absolute magnitude of the star, which is a tighter constraint on the luminosity than the log g_⋆ or a/R_⋆ constraints.

As regards apparent magnitude, the TASS (Droege et al. 2006) photometry for this star is V_TASS = 9.587  ±   0.071 and I_TASS = 8.357  ±   0.050, yielding (V − I)_TASS = 1.230  ±   0.087. On the other hand, the Yonsei–Yale (YY) stellar evolution models (Yi et al. 2001; Demarque et al. 2004) indicate (V − I)_YY,expect = 0.98  ±   0.03 for a star with T_eff⋆ ∼ 4850 K (based on the SME results). The 3σ inconsistency cannot be explained by interstellar reddening, since this is a close-by star. Thus, we opted not to use the TASS photometry in this analysis.

The 2MASS catalog (Skrutskie et al. 2006) provided much better agreement. The magnitudes reported in the 2MASS catalog have to be converted to the standard ESO system, in which the stellar evolution models (specifically the YY models) specify the colors. The reported magnitudes for this star are J_2MASS = 7.608  ±   0.029, H_2MASS = 7.131  ±   0.021, and K_2MASS = 7.009  ±   0.020, which is equivalent to J = 7.686  ±   0.033, H = 7.143  ±   0.028, and K = 7.051  ±   0.021 in the ESO photometric system (see Carpenter 2001). Thus, the converted 2MASS magnitudes yield a color of (J − K) = 0.635  ±   0.043 that is within 1σ of the expected (J − K)_YY,expect = 0.59  ±   0.02. We thus relied on the 2MASS K apparent magnitude and the parallax to derive an absolute magnitude of M_K = 4.18  ±   0.07. The choice of K band was motivated by the longest wavelength with smallest expected discrepancies due to molecular lines in the spectrum of this K4 dwarf.

In practice, the isochrone search for the best-fit stellar parameters was done in a Monte Carlo way, by assuming Gaussian uncertainties for the Hipparcos parallax, T_eff⋆, [Fe/H], the apparent 2MASS magnitudes, and the conversion coefficients by Carpenter (2001) that transform the 2MASS magnitudes to the standard K band. Lacking information, possible correlations between some of these parameters (e.g., T_eff⋆ and [Fe/H]) were ignored. A large set (∼5000) of random Δ (distance modulus), T_eff⋆, [Fe/H], and K values were generated, and for each combination we searched the stellar evolutionary tracks of the YY models for the best-fit stellar model parameters (such as M_⋆, R_⋆, etc.). For an unevolved K star there was no ambiguity in the solution, i.e., we did not enter a regime where isochrones cross each other. Certain parameter combinations in the Monte Carlo search did not match any isochrone. In such cases (∼14% of all trials) we skipped to the next randomly drawn parameter set. At the end. we derived the mean values and uncertainties of the physical parameters based on their a posteriori distribution. We also refined the stellar surface gravity. The new value log g_⋆ = 4.59  ±   0.03 agrees well with the earlier SME value, and has much smaller uncertainty confirming our previous assumptions that the absolute magnitude for this star is a better luminosity indicator than the surface gravity.

We then repeated the SME analysis by fixing log g_⋆ to the new value, and only adjusting T_eff⋆, [Fe/H] and vsin i. This second iteration yielded T_eff⋆ = 4780  ±   50 K and [Fe/H] = +0.31  ±   0.05. It also yielded vsin i = 0.3  ±   0.5 km s⁻¹, but for the same reason as given earlier we adopt the relaxed range vsin i = 1.5  ±   1.5 km s⁻¹. We accepted the above values as the final atmospheric parameters for this star. We then also repeated the isochrone search for stellar parameters, yielding M_⋆ = 0.809^+0.020_−0.027 M_☉, R_⋆ = 0.752  ±   0.021 R_☉, and L_⋆ = 0.26  ±   0.02 L_☉. Along with other stellar parameters, these are summarized in Table 4. The stellar evolutionary isochrones for metallicity [Fe/H] = +0.31 are plotted in the right panel of Figure 5, with the final choice of effective temperature T_eff⋆ and the absolute magnitude M_K marked, and encircled by the 1σ and 2σ confidence ellipsoids.

Zoom In Zoom Out Reset image size

The justification for our choice of using the parallax as the luminosity indicator is demonstrated by Figure 5. In the left panel, we plot a set of [Fe/H] = +0.31 Yale isochrones as a function of effective temperature, with the vertical axis being log g_⋆, and the observed values with their respective error ellipsoids overlaid. The middle panel shows the same set of isochrones with the a/R_⋆ luminosity indicator on the vertical scale. Here the geometric semimajor axis a/R_⋆ is determined from the photometric transit and RV data (as shown later in Section 3.2), and its error is significantly increased by the uncertainties in the eccentricity of the RV data. The right panel shows the same set of isochrones, with vertical axis being the absolute K magnitude. Overlaid is the observational constraint based on the apparent magnitude and Hipparcos parallax, along with the respective 1σ and 2σ error ellipsoids. It is clearly seen that the right panel imposes the tightest constraint on the stellar parameters. We note here that an isochrone search using the reported J − K color (instead of T_eff⋆ as color indicator) also agrees with the evolutionary models; however, the relative volume of its confidence ellipsoid is somewhat larger.

The effective temperature from the SME analysis and the surface gravity derived above correspond to a K4V star using Gray (1992). The B − V color index from the YY isochrones is 1.063  ±   0.024, also consistent with the K4 spectral type of Gray (1992). We also ran the Ramírez & Meléndez (2005) and Casagrande et al. (2006) temperature calibrations in reverse to back out the B − V required to produce the SME temperature. We got B − V = 1.025  ±   0.023 and 1.067  ±   0.025, respectively. For the final value, we accepted their average: B − V = 1.046  ±   0.024.

3.1.1. Baraffe Isochrones

We also investigated the dependence of stellar parameters on the choice of isochrones. Since HAT-P-11 is a K dwarf, we used the Baraffe et al. (1998) isochrones that are usually a better choice for late-type dwarfs. Baraffe et al. (1998) presented three sets of isochrones with different mixing length parameters. The one with α = 1.0 is a better match for low-mass stars (such as late K or M dwarfs), the one with α = 1.9 matches the Sun, and the third set with α = 1.5 is in between. In Figure 6, we plot these isochrones as a function of T_eff⋆, with the vertical axis being absolute K magnitude in the Bessel Brett system. Both the 2MASS K_s band and the Baraffe CIT systems were transformed to K using the relations in Carpenter (2001). The Baraffe isochrones for α>1.0 are given only for solar metallicity, while those for α = 1.0 are available for solar and sub-solar (Baraffe et al. 1997) metallicity only, but not for the metal-rich composition of HAT-P-11 ([Fe/H] = +0.31  ±   0.05). To investigate the effect of metallicity, we assumed the same qualitative behavior (i.e., opposite shift in the absolute magnitude–temperature plane) for isochrones with metallicity increased by +0.3, independent of their mixing lengths. In Figure 6, we overplot a metal-poor ([Fe/H] = −0.3, α = 1.0) isochrone, and conclude that an opposite change in metallicity to match that of HAT-P-11 would move the α = 1.0, [Fe/H] = +0.3 isochrone away from the observational values of K and T_eff⋆. Because of the lack of metal-rich models from Baraffe, and the lack of the knowledge of the proper mixing length parameter, we omit any quantitative conclusions from these models. Altogether, we found a better match with the YY isochrones, and thus accepted YY-based stellar parameters as final values (Table 4).

Zoom In Zoom Out Reset image size

As noted by Sozzetti et al. (2007) and Torres et al. (2008), a possible luminosity indicator for host stars of transiting planets is the a/R_⋆ quantity, where a is the relative semimajor axis, and R_⋆ is the radius of the host star. Here a/R_⋆ is in simple relation with the mean stellar density, if we assume that the planetary mass is much smaller than the stellar mass:

$$\begin{equation} \rho _\star = 0.01892 \frac{(a/R_\star)^3}{(P/{\rm day})^2}\,\rm g\,cm^{-3}. \end{equation} \tag{ 1 }$$

Analysis of the transit light curve yields—among other parameters—the quantity ζ/R_⋆  that is related to the time spent in between the planetary center crossing the limb of the star as T_dur = (2ζ/R_⋆)⁻¹. For circular orbits and equatorial transits a/R_⋆ is a function of ζ/R_⋆ and the orbital period, and thus can be determined in a straightforward way. For eccentric orbits, the relation between ζ/R_⋆ and a/R_⋆, as based on Tingley & Sackett (2005):

$$\begin{equation} \frac{\zeta }{R_\star }=\frac{a}{R_\star }\frac{2\pi }{P}\frac{1}{\sqrt{1-b^2}} \frac{1+h}{\sqrt{1-e^2}}. \end{equation} \tag{ 2 }$$

Thus, in addition to the light-curve parameters, a/R_⋆ also depends on the orbital eccentricity and argument of pericenter, and uncertainties in these parameters propagate into the error of a/R_⋆.

In the case of HAT-P-11, the RV data show a significant eccentricity of e = 0.198  ±   0.046 (Section 5). Although ζ/R_⋆ and the $\sqrt{1-b^2}$ terms have small uncertainties (1.0% and 4.8%, respectively), the uncertainty in the $(1+h)/\sqrt{1-e^2}$ term is higher, yielding a final value of a/R_⋆ = 14.6^+1.7_−1.4. The significant error can be credited to the uncertainties in the orbital parameters, as caused by the small RV amplitude and the stellar jitter. For comparison, if a/R_⋆ is calculated backward from the Hipparcos-parallax-based stellar evolutionary modeling, we get a/R_⋆ = 15.58^+0.17_−0.82. It is reassuring that this is consistent with that derived from the global modeling of the data, but the uncertainties are ∼3–4 times smaller. This justifies the choice of the Hipparcos-parallax-based luminosity indicator over using the a/R_⋆ constraint from global modeling of the data. This is also confirmed by comparing panels (b) and (c) of Figure 5, where the confidence ellipsoids from the parallax constraint (panel c) impose a tighter constraint on the isochrones.

The EPD (pre-TFA) HATNet light curve of HAT-P-11 shows significant periodic variations with P ≈ 29.2 days and a peak-to-peak amplitude of 6.4 mmag in I-band (Figure 7). This period is detected both in the autocorrelation function and the Fourier power spectrum of the light curve. The period is detected in both seasons covered by the HATNet light curve, and the signal remains in phase across both seasons as well. While this period is suspiciously close to the lunar cycle, we note that the light curve is relatively unchanged by the EPD procedure, i.e., these variations do not appear to correlate with any external parameters, including the sky background. The TFA procedure does suppress this signal when several hundred template stars are used. However, this is often the case when applying non-reconstructive TFA to other long-period variable stars such as Cepheids or Miras. To test whether or not this signal is due to some non-astrophysical systematic variation, we searched for other light curves exhibiting a strong peak in the Fourier power spectrum near f = 0.034 [c/d]. We found that nine out of the 5000 brightest stars with I ≲ 10.0 showed such a peak; however, none of these stars are in phase with HAT-P-11. Moreover, the distribution of peak frequencies does not show a pile-up at f = 1/P = 0.034 day⁻¹ relative to other frequencies, which one might expect to see if this variation were a systematic trend. As a final test, we attempted to recover the signal after applying different TFA template sets. We tried 100 disjoint template sets of 140 light curves. In all cases, the P ≈ 29.2 day signal was recovered in the autocorrelation function. We conclude that the signal is not a systematic variation that is present in the light curves of many other stars, and if it is not of astrophysical origin then it is due to a phenomenon that affects the light curve of HAT-P-11 in an apparently unique manner.

Zoom In Zoom Out Reset image size

A likely interpretation of the variation seen in Figure 7 is that it is due to the rotational modulation of starspots on the surface of HAT-P-11. The 6.4 mmag amplitude of the variation is consistent with other observations of spotted K dwarf stars. Note that there is no significant f₂ = 0.034 [c/d] periodicity in the RV data, when it is modeled as the combination of a Keplerian orbit plus a sinusoid component (see Figure 3). The secondary peaks in the autocorrelation function in the figure, plus the co-phasing of the 29 day variation over two observing seasons (left curve of Figure 7) indicate that individual starspots or starspot groups persist for at least several rotations.

If the variation is indeed due to rotational modulation by starspots, then the rotation period of HAT-P-11 is P ≈ 29.2 days. This may be compared with the rotation period predicted from the B − V color of the star and the Ca ii emission index S, using relations of Noyes et al. (1984). From Table 3, the median value of S observed from HIRES is 〈S〉 = 0.61. This value and B − V = 1.046  ±   0.024 as derived earlier yield a rotation period of P_calc = 24.2 days. The uncertainty in this calculation is difficult to quantify. The relations given in Noyes et al. (1984) were based on a theoretically motivated fit to the rotation period determined for a number of lower main-sequence stars from the rotational modulation of their chromospheric emission. For the 18 K stars in their sample with measured rotation periods P_obs, the rms difference between P_obs and the calculated value P_calc is 3.5 days. However, Noyes et al. (1984) pointed out that the empirical relation for P_calc as a function of S and B − V becomes rather unreliable for B − V ≳ 1.0, because of a paucity of stars with observed rotation periods in this color range. Hence, the true uncertainty is doubtless somewhat larger than 3.5 days. We conclude that a rotation period of 29.2 days is consistent with expectations from the star's B − V color and 〈S〉 value.

A 29 day rotation period of HAT-P-11, coupled with its radius, implies an equatorial velocity v_eq = 1.3 km s⁻¹. Let us assume that the stellar rotational axis is inclined at i ∼ 90° to the line of sight. This is supported by the geometry of our solar system and the generally similar geometries of the stellar systems with transiting planets whose projected rotational axis inclinations have been measured through the Rossiter–McLaughlin effect. We would then expect vsin i ∼ 1.3 km s⁻¹, consistent with the value of vsin i reported in Table 4.

Figure 7 shows the time history of the Ca ii S index. Figure 8 shows the very prominent emission cores of the Ca ii H and K lines observed at three levels of activity during the time of the HIRES observations. A long-term variation of S is apparent, with timescale close to the length of the current data set. As discussed later, this may be due to long-term variations of stellar activity analogous to the solar activity cycle. Phasing of the S data at a period of 29.2 days shows no evidence for periodic behavior. If the star is really rotating at a period of 29 days, this suggests that the chromospheric emission in the H and K lines is spread nearly uniformly in longitude over the star.

Zoom In Zoom Out Reset image size

The color B − V = 1.046 of HAT-P-11 and its median level of chromospheric emission 〈S〉 = 0.61 imply a chromospheric emission ratio R'_HK given by log R'_HK = −4.584 (Noyes et al. 1984). From this, we may crudely estimate an age for the star using the inverse square root relation between activity and age as originally posited by Skumanich (1972). Using a fit of this relation to the Sun, the Hyades, the Ursa Major group, and 412 individual lower main-sequence stars as derived by Soderblom et al. (1991), we determine a "chromospheric age" T_cr∼ 1.25 Gyr. The uncertainty in this estimate is difficult to quantify, especially because most of the stars in the sample had B − V < 1.0, but it does suggest that based on its chromospheric emission level the star is likely to be at the low end of the age range given in the isochrone fit discussed in Section 3.1.

As always in determining whether the signature from combined photometric and RV variations in a star is due to a transiting planet, it is necessary to exclude the possibility that the entire combined information is due to a set of circumstances that falsely give rise to the characteristic signature of a transiting planet. Following Torres et al. (2007), we first explored the possibility that the measured RV variations are caused by distortions in the spectral line profiles due to a nearby unresolved faint eclipsing binary, whose relatively large RV variations mix with the non-varying spectrum of the primary and give rise to an apparent RV signal with the observed small amplitude. In this case it has been shown (e.g., Queloz et al. 2001) that the BS of spectral lines in the blended spectrum varies in phase with the RV signal itself, with similar amplitude. We have carried out an analysis of the BS based on the Keck spectra as described in earlier HATNet detection papers (see Bakos et al. 2007), and do not see any statistically significant correlation between the line BSs and the measured radial velocities (confer Figure 4 and Table 3), thus providing no support for the hypothesis that the signal is caused by a blended eclipsing stellar system.

However, because of the small amplitude of the RV signal (and hence small amplitude of expected BS variations if the apparent RV signal is indeed due to a blend), coupled with the large jitter as seen in the RV residuals in Figure 4, the above test cannot completely rule out contamination by an unresolved binary system. Thus, we sought other independent ways of excluding blend scenarios. In the following subsections, we attempt to model the system as a hierarchical triple, where the smallest component is either a star or a Jupiter-sized planet, and also as a background eclipsing binary blended with the light of a foreground K dwarf. We show that none of these models are consistent with all of the available observations.

We consider the possibility that HAT-P-11 is a hierarchical triple system, having a deeper intrinsic eclipse of two bodies diluted by the light of the bright K dwarf. Note that the bright K dwarf with well known properties (parallax, proper motion, etc.) in this putative triple system is referred to as HAT-P-11 in the following discussion. To rule out the hierarchical triple scenario, we attempt to fit a blend model to the observations. There are two scenarios for the eclipsing system in the triple system. The first one is a stellar eclipsing binary. The second case is a transiting hot Jupiter orbiting a low-mass star, with the few percent deep transit diluted by HAT-P-11 to become only 4.2 mmag. We note that the shallower the transit event, the larger the variety of configurations that can match the observations. HAT-P-11 represents a transition toward millimag transit events, the blend analysis of which will certainly be challenging.

In each case (stellar binary, versus transiting hot Jupiter blend), we fit the follow-up z, r, and I_C-band light curves  together with the HATNet I-band light curve following a procedure similar to that described by Torres et al. (2005). We assume that the bright star (HAT-P-11), which is not eclipsed, has the mass, metallicity, age, and distance determined in Section 3. We also assume the components of the eclipsing system have the same metallicity, age, and distance as the bright star, i.e., they form a hierarchical triple. We fix the ephemeris, orbital eccentricity, and argument of periastron to the values determined in Section 5.4, and vary the masses of the two components and the inclination of the orbit. The downhill simplex algorithm is used to optimize the free parameters. For the case where one of the components is a planet, we vary the radius of the planet rather than its mass. For the eclipsing binary case, the magnitudes and radii of the stars are taken from the Padova isochrones (Girardi et al. 2000) applying the radius correction described by Torres et al. (2005). We use these isochrones to allow for stars smaller than 0.4 M_☉. For the planet case, we use the YY isochrones to be consistent with the single-star planet modeling.

Figure 9 shows the best-fit model for the eclipsing binary case and two illustrative models for the blended hot-Jupiter case overplotted on the z-band light curve. We are unable to fit the light curve using a combination of three, physically associated stars, but are able to fit it as a planet transiting one component of a binary star system, so long as the star hosting the planet has M > 0.72 M_☉. Interestingly, 0.72 M_☉ emerges as a critical mass (and corresponding critical radius, luminosity) for the eclipsed star, and splits the parameter space into two disjoint domains. As the mass of the host star is reduced below 0.72 M_☉ (and its radius is decreased), the only way to maintain the depth and duration of the already central (b = 0, i.e., longest possible) transit is to increase the radius of the planet. This, however results in a longer ingress and egress time than is allowed by the observations. Above 0.72 M_☉ the radius of the planet needed to fit the transit depth changes slowly, while raising the impact parameter of the planet can compensate to fit the transit duration (see Figure 9). While a model of this form fits the light curve, the V-band light ratio of the planet-hosting star to the brighter 0.81 M_☉ star would be >0.46. The spectrum of a second star this bright should have been easily detected in the Keck and DS observations, and therefore we can rule out both of these blend scenarios.

Zoom In Zoom Out Reset image size

An alternate model of an astrophysical false positive that would make the contaminating star much fainter relative to the main object is one involving a background eclipsing binary (chance alignment). We modeled this case using the same methodology described above, fitting the combined z-band follow-up light curves as the sum of the light from three stars. The properties of the main star were held fixed as before, and those of the eclipsing binary components were constrained to lie on the same isochrone, which for simplicity was taken to be the same as the main star. Extensive tests indicated that the detailed shape of the resulting synthetic light curve is fairly insensitive to the distance at which we place the eclipsing binary behind the target, and thus acceptably good fits to our z-band photometry can be achieved for a wide range of distances, given the measurement precision. Figure 10 shows an example of such a fit, in which the eclipsing binary is placed about 60 pc behind the main star and is composed of a 0.64 M_☉ star (spectral type K7V–M0V) orbited by a 0.13 M_☉ stellar companion. This model is not only consistent with the photometry, but it predicts an optical brightness for the primary of the eclipsing binary of only 3% relative to the main star. Detecting such a faint set of spectral lines in our Keck spectra would be challenging. Furthermore, if we place the eclipsing pair 110 pc behind the main star, the light curve fit is still about the same, but the relative brightness decreases by a factor of 2, making the eclipsing binary spectroscopically undetectable. Spectral line bisector variations predicted by this second model (see, e.g., Torres et al. 2005) would be at the level of the scatter in our actual measurements, and thus could not be entirely ruled out either.

Zoom In Zoom Out Reset image size

While this blend scenario appears to satisfy all observations, it implies the presence of an eclipsing binary ∼4.5 mag fainter in the optical very near the present location of our target. It is fortunate that HAT-P-11 has a large proper motion (0264 yr⁻¹; Perryman et al. 1997). Using Palomar Observatory Sky Survey plates from 1951 (POSS-I, red and blue plates), we can view the sky at the current position of HAT-P-11 unobstructed, since the target was 15'' away. Figure 11 shows a stamp from the POSS-I, POSS-II plates, and also a current observation with the FLWO 1.2 m telescope. The proper motion of HAT-P-11 is apparent. A number of faint stars are marked on the figure. Note that the faintest possible star that can cause the observed 4.2 mmag variation in the combined light curve would be 6 magnitudes fainter than HAT-P-11, and in this extreme situation the blended faint star would need to completely disappear during its transit. There is no faint star down to ∼19 mag within ∼5'' of the current position of HAT-P-11. The closest star marked "2" in Figure 11 is 2MASS 19505049+4805017, with z ≈ 14.4, or about z ≈ 5.6 mag fainter than HAT-P-11. This is well resolved in some of our follow-up photometry observations, such as on the 2007 September 2 night, and its brightness was constant, with rms much smaller than the required amplitude of its variation should be (∼1.3 mag) to cause the observed 4.2 mmag combined dip.

Zoom In Zoom Out Reset image size

In summary, we (1) have not seen bisector variations correlated with the RV, (2) investigated unresolvable hierarchical triple systems (blended eclipsing binary or transiting hot Jupiter), and found all to be incompatible with the photometric and spectroscopic data together, and (3) were able to exclude chance alignment of a background eclipsing binary based on the high proper motion of the star. We therefore conclude that the transit signal, and the synchronized RV signal, indeed are both due to a sub-stellar companion transiting the K4 dwarf HAT-P-11.

Our model for the follow-up light curves used analytic formulae based on Mandel & Agol (2002) for the eclipse of a star by a planet, where the stellar flux is described by quadratic limb darkening. We did not attempt to fit for the limb-darkening coefficients, as they require very high quality data to resolve their degeneracy with other parameters. Instead, the limb-darkening coefficients were derived from the second iteration SME results (Section 3), using the tables provided by Claret (2004) for z, r (FLWO 1.2 m), and I_C bands (Konkoly 0.6 m Schmidt). The transit shape was parameterized by the normalized planetary radius p ≡ R_p/R_⋆, the square of the impact parameter b², and the reciprocal of the half duration of the transit¹⁶ζ/R_⋆. We chose these parameters because of their simple geometric meanings and the fact that these show negligible correlations (see Carter et al. 2008; Bakos et al. 2007; Pál 2008b).

Our model for the HATNet data was a simplified version of the Mandel & Agol (2002) analytic functions, because (1) the number of in-transit data points in the HATNet light curve is significantly less than the same number for the follow-up light curves, and (2) the individual errors on the points are much worse for HATNet than for the follow-up. Thus, while the HATNet data are efficient in constraining the ephemeris, they are not comparable to the follow-up data in determining other transit parameters that depend on the exact shape of the transit (such as depth, or duration of ingress). For these reasons, we employed a simple transit model that neglects limb darkening. This model function, which we label as "P1P3," originates from an expansion of the exact Mandel & Agol (2002) model function by Legendre polynomials, resulting in a functional form (omitting scaling factors) of the ingress (and symmetrically, the egress) of x → (21x − 5x³)/16, where the independent variable x is scaled to be 1 at first (third) contact, and −1 at second (fourth) contact. For reference, the trapeze model has a simple linear function describing the ingress: x → x. The "P1P3" model approximates the Mandel & Agol (2002) function to better than 1%, and our experience shows that such an approximation yields the same timing precision as more sophisticated model functions. In addition, using the simplified model yields much faster computation. A similar transit model approximation using hyperbolic tangent functions was described earlier by Protopapas et al. (2005). The depth of the HATNet transits was also adjusted in the fit (independent of the follow-up data), and we found that the depth derived from the HATNet light curve is comparable (Section 5.4) to the depth of the (EPD and TFA reconstructed) follow-up transit measurements, namely, the ratio of the two was B_inst = 0.95  ±   0.08. In general, this "instrumental blend" parameter has to be adjusted independently since the possible contamination by nearby stars in the wide-field images may yield shallower transits.

The adopted model for the RV variations was described earlier in Section 2.4. The RV curve was parameterized by an eccentric Keplerian orbit with semi-amplitude K, RV zero-point γ, Lagrangian orbital elements (k, h) = e × (cos ω, sin ω), plus a linear trend G₁.

We assumed that there is a strict periodicity in the individual transit times. In practice, we introduced the first transit center as T_A = T_c,−231 and the last transit center as T_B = T_c,+91, covering all of our measurements with the HATNet telescopes, the FLWO 1.2 m telescope and the Schmidt telescope. The transit center times for the intermediate transits were interpolated using these two epochs and the N_tr transit number of the actual event under the assumption of no transit timing variations (TTVs). The model for the RV data contained the ephemeris information through the T_c,−231 and T_c,+91 variables; for instance, the period assumed during the fit was P = (T_c,+91 − T_c,−231)/(91 + 231). The other coupling between the transit photometry and the RV data is via the k and h Lagrangian orbital elements that determine the relation between the photometric and orbital ephemeris, plus have a minor effect on the transit shape. Altogether, the 12 parameters describing the physical model were T_c,−231, T_c,+91, R_p/R_⋆, b², ζ/R_⋆, K, γ, k = ecos ω, h = esin ω, and G₁. As mentioned above, the parameters were extended with the instrumental blend factor B_inst, and the HATNet out-of-transit (OOT) magnitude, M_0,HATNet.

A "joint" fit of the physical model to the global data set has been performed routinely for recent discoveries (Bakos et al. 2009; Latham et al. 2009). Here we extended our physical model that describes the transit events and the radial velocities with an instrumental model that describes the systematic variations of the data. Since the HATNet photometry has been already EPD and TFA corrected (Section 2.1), and the Keck RVs have been investigated for trends (Section 2.4), but deemed to have insufficient number to decipher any significant systematic variation (50 points compared to $\mathcal {O}(10^4)$ points for photometry), we only modeled systematic variations of the follow-up photometry (Section 2.3). As mentioned earlier, the effect of systematics on the shallow transit in the photometric follow-up data is enhanced, and can impact the basic planetary parameters, such as the planetary radius. We experimented with several models, since the treatment of systematics may improve precision, but not necessarily the accuracy of the results. For reference, we call the simplest model, with no EPD and no TFA, as model "E0T0" (for technical details, see the Appendix).

Model "E0TRG" was independent (i.e., not simultaneous) reconstructive global TFA of the FLWO 1.2 m z-band data (eight nights, see Table 1), assuming a trapeze-shaped $m_0(\vec{p},t)$ model function in Equation (A2) with three parameters (total duration, ingress duration, and depth). The period was iteratively refined using both the FLWO 1.2 m and the HATNet data, but reconstructive TFA was only performed on the FLWO photometry. We used the magnitude values transformed to the selected reference frame for each night (see Section 2.3), and the zero-point shift for each night was determined separately. We increased the number of template stars till the S/N of the trapeze reached its maximum (S/N = 39), requiring 246 stars with z ≲ 14. Note that this was a global TFA, in the sense that there was one c_k TFA coefficient per star in Equation (A3) for all nights combined (see Appendix A.2). The reconstructed light curve was then fed into the joint fit.

Model "ELT0" had EPD parameters in the simultaneous fit, together with the parameters of the physical model, but no TFA-based detrending, i.e., all the systematic variations in the follow-up photometry were modeled as due to external parameters. We chose five such parameters, namely the hour angle (characterizing a monotonic trend that linearly changes over time), the square of the hour angle, the stellar profile sharpness parameter, S = (2.35/FWHM)², the X and Y pixel position on the chip (each with two coefficients). The exact functional form of the above parameters contained eight coefficients, including the auxiliary OOT magnitude of the individual events. The EPD parameters were independent for all 10 nights, implying 80 additional coefficients in the global fit.

Model "E0TL" had local TFA parameters in the instrumental model for each night, and no EPD was used. We selected 23 TFA template stars that were present in all observations, representing faint and bright stars spread across the chip. The total number of c_k coefficients in the simultaneous fit was thus 230. Model "E0TG" had global TFA parameters for all nights, using 211 stars that were present on all frames. Note that without simultaneous fit using an underlying physical model (m_j(t_i)) in Equation (A2)) that takes into account the different shape of the light curves acquired through different photometric bands, such a global TFA fit would not be possible.

We tested variants of the above methods. Model "ELTL" was a simultaneous fit of EPD and TFA parameters in local EPD mode and local TFA mode with altogether 80 + 230 parameters in addition to the 12 physical parameters. Model "ELTG" was a simultaneous fit of local EPD and global TFA parameters with an additional 80 + 211 instrumental parameters. The parameter sets were extended in all cases with the OOT magnitudes for each follow-up light curve and the HATNet photometry. The number of fitted parameters was much smaller than the number of data points (∼5000).

It is computationally rather intensive to perform a fit with ∼300 parameters and to determine their respective error distributions.

Since the majority of the fitted parameters appear as linear terms in the final form of the model functions, we minimized χ² in the parameter space by using a hybrid algorithm, namely, by combining the downhill simplex method (aka AMOEBA; see Press et al. 1992) with the classical linear least-squares (CLLS) algorithm. The simplex itself was propagating in the "nonlinear" hyperplane of the parameter space, while in each point of the hyperplane the value of χ² was computed using classical linear CLLS minimization. This method yielded a very good convergence even with the large number of free parameters.

Uncertainties of the parameters were derived using the Markov Chain Monte Carlo method (MCMC; see Ford 2006) by starting two types of chains from the best-fit value. The first one was the classical MCMC chain, letting all parameters vary. The other type of chain was generated by only allowing the variation in a hyperplane of the parameter space consisting of the nonlinear parameters and of some additional linear parameters where error calculation was desired (examples are K, γ or G₁). In this "Hyperplane-CLLS" chain, the remaining linear parameters were determined at each step via CLLS. The H-CLLS chain has the advantage that transition probabilities are much higher in the separated hyperplane, and thus the convergence and computing time requirements are also smaller by an order of magnitude. Furthermore, the error distribution of parameters for the hyperplane variables was identical to those determined by classical MCMC. The distribution for the parameters outside the hyperplane was distorted, but these were considered as auxiliary, where such distribution is irrelevant (e.g., the EPD and TFA coefficients). Therefore, the final distribution of the physical parameters was derived by using H-CLLS chains.

For both types of the chains the a priori distributions of the parameters were chosen from a generic Gaussian distribution, with eigenvalues and eigenvectors derived from the Fisher covariance matrix for the best-fit value. All functions that appear in the physical and instrumental model are analytical, and their partial derivatives are known. The parametric derivatives of the transit light-curve model of Mandel & Agol (2002) are found in Pál (2008b), while the parametric derivatives of the eccentric RV curves can be obtained by involving the implicit function theorem, and can be expressed as functions of the eccentric anomaly. Thanks to these analytic properties, the derivation of the Fisher matrix was straightforward.

We performed a joint fit on the HATNet light curve, the follow-up photometry, and the Keck radial velocities (the "data"), using the physical model described in Section 5.1, as extended by the various models for systematics (Section 5.2), and using the algorithms described in Section 5.3.

We carried out the fit considering each model describing the systematics in the follow-up photometry separately in order to investigate their effect on the final results. As mentioned earlier, this is important due to the very shallow transit signal and the large relative amplitude of systematic variations. We used the m_f-fitted magnitude values (Section 2.1) for the photometry follow-up light curves (without prior EPD or TFA). Since only the follow-up photometric data component was described by the different systematic models (see Section 5.2), we focused on those light-curve parameters that are affected by the photometric follow-up, namely, the transit center time of the last event T_c,+91, the radius ratio of the planet to star, the (square of the) impact parameter and the transit duration parameter ζ/R_⋆. The other adjusted parameters were not relevant in this comparison, since they are not affected by the follow-up light curves. Note that with the exception of the E0TRG model, all EPD and TFA parameters were fit simultaneously with the physical model parameters.

The best-fit values and the respective uncertainties are summarized in Table 5. It is clear that EPD and TFA both significantly reduce the uncertainties in all the parameters. Note that the implementation of E0TRG was somewhat sub-optimal. Because the reconstructive trapeze shape was not bandpass dependent (as opposed to the full Mandel & Agol (2002) model used in the other cases), only a fraction of the follow-up were used, namely, all photometry in z band. When using EPD and TFA together, the global TFA (ELTG) performs better than local TFA (ELTL) for each night. The improvements by EPD (ELT0) and TFA (E0TG) separately are roughly the same. The best model, based on the formal error bars, is the EPD and global TFA together (ELTG); here the unbiased error in the light-curve parameters decreased by a factor of 1.5.

The final value for the parameters also changed by 3σ in the sense that both R_p/R_⋆ and b² decreased. This means that the treatment of the systematics recovers a sharper transit with smaller impact parameter and smaller planetary radius. If the stellar parameters were determined using the a/R_⋆ luminosity constraint (Section 3.2), then ignoring the correction for systematics would lead to higher impact parameter, larger R_p/R_⋆, and smaller a/R_⋆, corresponding to smaller density (Equation (1)). One may expect that the effect of systematics is slowly diminished by accumulating more data taken with different instruments. This is consistent with our experience, that is the initial follow-up light curves without proper treatment of systematics indicated larger impact parameters, and a lower density dwarf star that was inconsistent with the Hipparcos-based stellar parameters (Section 3.1). By accumulating more data this inconsistency diminished (Section 3.2). Because both EPD and TFA point toward this direction, and also deliver the smallest formal uncertainty, we expect that the most accurate model will be the ELTG in Table 5. This conjecture will be decided with accurate ground-based follow-up photometry (such as presented recently in Johnson et al. 2009), and ultimately when the Kepler space mission returns the light curve for HAT-P-11. This will have important implications on optimal trend removal for other transiting systems, where such a space-born "reference" will not be available.

The results for the simultaneous fit using the ELTG systematic model for the follow-up light curves were the following: T_c,−231 = 2453217.75466  ±   0.00187 (BJD), T_c,+91 = 2454605.89132  ±   0.00032 (BJD), K = 11.6  ±   1.2 m s⁻¹, k ≡ ecos ω = 0.201  ±   0.049, h ≡ esin ω = 0.051  ±   0.092, R_p/R_⋆ = 0.0576  ±   0.0009, b² = 0.120  ±   0.087, ζ/R_⋆ = 22.15  ±   0.15 day^-1, and γ = −0.4  ±   0.8 m s⁻¹, B_instr = 0.95  ±   0.08, G₁ = 0.0297  ±   0.0050 m s⁻¹ day⁻¹, M_0,HATNet = 8.35892  ±   0.00003 (instrumental HATNet OOT magnitude). The combined and phase-binned light curve of all z-band FLWO 1.2 m observations is shown in Figure 12. The planetary parameters and their uncertainties can be derived by the direct combination of the a posteriori distributions of the light curve, RV, and stellar parameters (see also Pál et al. 2008a). We found that the mass of the planet is M_p = 0.081  ±   0.009 M_J = 25.8  ±   2.9 M_⊕, the radius is R_p = 0.422  ±   0.014 R_J = 4.73  ±   0.16 R_⊕, and its density is ρ_p = 1.33  ±   0.20 g cm⁻³. The final planetary parameters are summarized at the bottom of Table 6.

Zoom In Zoom Out Reset image size

Table 6. Orbital and Planetary Parameters

Parameter	Value
Light Curve Parameters
P (days)	4.8878162  ±   0.0000071
Tc (BJD)	2454605.89132  ±   0.00032
T14 (days) a	0.0957  ±   0.0012
T12 = T34 (days) a	0.0051  ±   0.0013
a/R⋆	15.58+0.17−0.82
ζ/R⋆	22.15  ±   0.15
Rp/R⋆	0.0576  ±   0.0009
b ≡ acos i/R⋆	0.347+0.130−0.139
i (deg)	88.5  ±   0.6
Spectroscopic Parameters
K (m s−1)	11.6  ±   1.2
γ (km s−1)	−0.4  ±   0.8
G1 (m s−1/day)	0.0297  ±   0.0050
kRVb	0.201  ±   0.049
hRVb	0.051  ±   0.092
kCc	0.190  ±   0.046
hCc	−0.016  ±   0.056
e	0.198  ±   0.046
ω	3552  ±   173
Secondary Eclipse Parametersd
Ts (BJD)	2454608.96  ±   0.15
Ts,14	0.1006  ±   0.0130
Ts,12	0.0054  ±   0.0013
Planetary Parameters
Mp (MJ)	0.081  ±   0.009
Rp (RJ)	0.422  ±   0.014
C(Mp, Rp)d	0.025
ρp (g cm−3)	1.33  ±   0.20
a (AU)	0.0530+0.0002−0.0008
log gp (cgs)	3.05  ±   0.06
Teq (K)	878  ±   15
Θ	0.025  ±   0.003
Fper (erg s−1 cm−2) e	2.04 × 108  ±   2.78 × 107
Fap (erg s−1 cm−2) e	9.11 × 107  ±   9.53 × 106
〈F〉 (erg s−1 cm−2) e	1.34 × 108  ±   9.38 × 106

Notes. ^aT₁₄: total transit duration, time between first to last contact; T₁₂ = T₃₄: ingress/egress time, time between first and second, or third and fourth contact. ^bLagrangian orbital elements, based purely on the RV data. ^cRefined values of h and k, derived from RV analysis and the constraint given by the $C_{\rm kh}\equiv \sqrt{1-e^2}/(1+h)$ value, resulting from the light-curve modeling and stellar evolution analysis. ^dCorrelation coefficient between the planetary mass M_p and radius R_p. ^eOccultation parameters and incoming flux per unit surface area in periastron, apastron, and average for the orbit were calculated using the refined values k_C and h_C.

Download table as:  ASCIITypeset image

An interesting aspect of the HAT-P-11 system is that the light-curve analysis, the Hipparcos parallax, and the theoretical isochrones together provide an extra constraint on the orbital eccentricity. In our analysis, the orbital eccentricity and argument of pericenter are characterized by the Lagrangian orbital elements k = ecos ω = 0.201  ±   0.049 and h = esin ω = 0.051  ±   0.092. If the analysis relied purely on the RV data, then the errors in k and h would be similar, the error of h being somewhat smaller (see Section 2.4). However, when the data are complemented with photometry of transit events (as in our case), the period and phase of the orbit become tightly constrained. As a result, the k and h orbital elements will have different uncertainties, since the mean longitude at the transit is constrained more by cos ω than sin ω (this is due to the same underlying reason that the phase lag between the secondary and primary transits is also proportional to ecos ω and does not strongly depend on esin ω). Indeed, the a posteriori distribution of h and k as yielded by the joint analysis (Section 5.4) is asymmetric, and the uncertainty in k is half of the uncertainty in h (Figure 13, large oval around the filled diamond).

Zoom In Zoom Out Reset image size

By re-arranging Equation (2), we can define C_kh as

$$\begin{equation} C_{\rm kh}\equiv \frac{\sqrt{1-h^2-k^2}}{1+h}= \frac{a}{R_\star }\left(\frac{\zeta }{R_\star }\right)^{-1}\frac{n}{\sqrt{1-b^2}}, \end{equation} \tag{ 3 }$$

where n = 2π/P is the mean motion (see also Ford et al. 2008). The period P, impact parameter b, and ζ/R_⋆ are well constrained by the photometry of multiple transit events. The normalized semimajor axis a/R_⋆ is determined from the apparent magnitudes, Hipparcos parallax and isochrones, as described in Section 3.2. These together yield an estimate of C_kh = 0.931  ±   0.063, independent of the RV data, and thus provide an additional constraint on h and k in the form of a hyperbola on the k − h plane. Taking into account the uncertainties in C_kh, this appears as a stripe in the k − h plane (Figure 13). This constraint can be used to decrease the uncertainty of h by a factor of 2, so that it becomes similar to that of k.

In practice, we have a Monte Carlo distribution for k, h, and C_kh. For each k and h value in MC distribution, we determine the closest k', h' point of the hyperbola (there is a unique solution for this). Since the value for C_kh constrains only one degree of freedom, we can treat (k, h) and (k', h') as independent. Therefore, we define k_C and h_C as the mean of k, k' and h, h', respectively, and consider these as an improved estimate of the orbital eccentricity parameters. The result is an MC distribution for k_C and h_C, where their mean values and errors can be derived: k_C = 0.190  ±   0.046 and h_C = −0.016  ±   0.056. Indeed, the errors in k_C and h_C are now comparable.

Based on k_C and h_C, the refined orbital eccentricity is e = 0.198  ±   0.046, and the argument of pericenter is ω = 3552  ±   173. Note that the above method is unique, and can be applied for eccentric transiting systems with well known parallax but poorly determined RV orbit. It is likely that the Kepler mission will find similar cases. A counter example is HAT-P-2b (Bakos et al. 2007), where the RV amplitude is large compared to the uncertainties (1 km s⁻¹), and the error in the parallax is larger, thus h cannot be refined. The current analysis points toward an even more general way of simultaneously fitting to all data, including the stellar isochrone search with various constraints (e.g., parallax information with error bars).

The stellar parameter determination (Section 3), the blend analysis (Section 4), and the global modeling of the data presented in Section 5 together show the detection of a 25.8  ±   2.9 M_⊕ and 4.73  ±   0.16 R_⊕ planet that orbits a bright (V = 9.59) K4 dwarf star on an eccentric (e = 0.198  ±   0.046) orbit with P = 4.8878 day period, causing a 11.6  ±   1.2 m s⁻¹ RV variation of the host star and a 4.2 mmag transit as it passes in front of the star. The planet is somewhat more massive than our own Uranus (14.3 M_⊕) and Neptune (17.1 M_⊕), and slightly more massive than GJ 436b (22.6 M_⊕). Based on its mass and radius, it is justified to use the metonym super-Neptune for classification of HAT-P-11b, just like for GJ 436b, the only similar object known so far. For the location of ∼Neptune-mass objects, see Figure 14.

Zoom In Zoom Out Reset image size

When compared to models of Fortney et al. (2007), HAT-P-11b is much smaller in radius than similar mass planets with 50% rock/ice core and 50% H/He envelope, and HAT-P-11b is much bigger than pure rock/ice planets without a gas envelope. The mismatch with such pure-rock/ice "super-Earths" is confirmed by models of Valencia et al. (2007). Overlaid in Figure 14 are super-Earth planet models with Earth-like composition (lower dotted line) and 50% by mass $\rm H_2 O$ content and Earth-like Si/Fe ratio (upper dotted line). Both models fall clearly below the observational values for HAT-P-11b. As the expected limit for gravitational capture of H/He is ∼10 M_⊕ (Rafikov 2006), it is hardly surprising that HAT-P-11b has some H/He, even if evaporation played an important role during its evolution. HAT-P-11b receives only modest irradiation when compared to hot Jupiters. The time-averaged flux over its eccentric orbit is 1.34 × 10⁸  ±   9.38 × 10⁶ $\rm erg\,s^{-1}\,cm^{-2}$, and the equivalent semimajor axis where similar irradiation would be received from a Sun-like star is a_equiv ∼ 0.1 AU. Using the models from Fortney et al. (2007), the radius of HAT-P-11b is indeed consistent with a 90% heavy element planet with an irradiation corresponding to a_equiv = 0.1 AU. In these models, the heavy elements are located in the core, and the envelope is metal-free.

We also compared the observational values to the thorough theoretical work of Baraffe et al. (2008). Planetary isochrones for various metal content, age, and irradiation are also plotted in Figure 14. As Fortney et al. (2007) note, planetary radii and their evolution are scarcely affected when irradiation is small, and a_equiv ≳ 0.1 AU; thus, it is expected that HAT-P-11b will match non-irradiated models better, and indeed, the best match is with extreme metal rich non-irradiated Baraffe models with Z = 0.9 fraction of heavy elements (rock and ice), and a 10% H/He envelope. The non-irradiated metal-rich models of various ages between 3 and 7 Gyr lie very close to each other and are all within the observational error bars. However, the non-irradiated models with significant (50%) H/He envelope (top of figure) are far from matching the observations. Baraffe et al. (2008) note that, in general, the distribution of the heavy elements (core versus envelope) can have considerable impact on the planetary radius, but in this extreme metal-rich scenario, the freedom to vary their distribution is limited. Altogether, HAT-P-11b appears to be a super-Neptune planet with Z = 0.9 and with a 10% H/He envelope.

It is interesting to compare HAT-P-11b to GJ 436b—the only other super-Neptune with known mass and radius. These two planets are very similar in some physical properties in spite of the very different environment. There are numerous system parameter determinations for GJ 436b available in the literature (see compilation and Table A11 of Southworth 2009). Starting from the discovery of the transits by Gillon et al. (2007a) (22.6 M_⊕, 3.95 R_⊕), planetary masses range from M_p = 22.26 M_⊕ to 24.8 M_⊕, and planetary radii range from 3.95 R_⊕ to 4.9 R_⊕ (Bean et al. 2008). Some of these earlier estimates, including that of Bean et al. (2008), are plotted in Figure 14. The RV amplitude of GJ 436 (K = 18.34  ±   0.52 m s⁻¹; T08) is larger and better determined than HAT-P-11b (K = 11.6  ±   1.2 m s⁻¹), and the transits are deeper (6.5 mmag versus 4.2 mmag), and have even been observed from space by Spitzer (Gillon et al. 2007b; Deming et al. 2007) and Hubble Space Telescope (Bean et al. 2008). The reasons for the larger K value and transit depth of GJ 436, however, are due to the host star being a small mass and radius M dwarf (M_⋆ ≈ 0.45 M_☉, R_⋆ ≈ 0.46 R_☉, T08 and Southworth 2009). Uncertainties in both the planetary mass and radius remain considerable because of the uncertainty in the stellar parameter determination for M dwarf stars (T08). Altogether, the masses and radii of GJ 436b and HAT-P-11b are very close (to within 2σ). GJ 436b has a slightly smaller and better determined mass. HAT-P-11b has a comparable and somewhat better determined radius.

There are other similarities between the two planets. Both orbit on short period mildly eccentric orbits, GJ 436b: e = 0.16  ±   0.02 (Gillon et al. 2007a) and HAT-P-11b: e = 0.198  ±   0.046. The fact that they have not circularized yet is interesting, as the circularization timescale appears shorter than the age of the host stars (e.g., Matsumura et al. 2008, and references therein), although the stellar ages are ill determined due to the unevolved M and K dwarf host stars. As Matsumura et al. (2008) noted on the origins of eccentric close-in planets, some of these planets may have larger than expected Q specific dissipation functions (depending on the planetary structure, and characteristics of the tides). The fact that neither of these hot Neptunes has circularized yet may tell us something about the formation and planetary structure of Neptunes. We note that the putative HAT-P-11c planet (Section 6.2) causing the long-term drift in the data appears to have a large semimajor axis, because the RV drift does not show any significant nonlinear coefficient over the one year timescale of the observations (Section 2.4). Thus, it is unlikely that it is causing eccentricity pumping. Whether the Kozai mechanism (Fabrycky & Tremaine 2007) plays an important role will be subject to further investigations when the long-term drift in the RVs is better characterized. Both GJ 436b and HAT-P-11b have similar surface gravities, log g_p = 3.11  ±   0.04 (cgs, T08) for GJ 436b versus log g_p = 3.05  ±   0.06 for HAT-P-11b. They also appear to have similar mean densities, $\rho _{p}= 1.69\,\,{\pm}\,\, ^{0.14}_{0.12}\,\rm g\,cm^{-3}$ (T08) for GJ 436b and $1.33\,\,{\pm}\,\, 0.20\,\rm g\,cm^{-3}$ for HAT-P-11b.

It is surprising that in spite of these similarities (mass, radius, density, surface gravity, eccentricity), the two planets are in grossly different environments. Consequently, the GJ 436b–HAT-P-11b pair must yield tight constraints on theories that should reproduce their properties in spite of the different environments. GJ 436b orbits an M dwarf with half the mass and radius of HAT-P-11, one-tenth of the luminosity (L_⋆ = 0.026 L_☉ versus 0.26  ±   0.02 L_☉), and colder effective temperature of T_eff⋆ = 3350  ±   300 K (T08) versus 4780  ±   50 K for HAT-P-11b. The smaller semimajor axis of GJ 436b compensates to some extent for the lower luminosity host star, but still, the irradiation they receive is grossly different (GJ 436b: $3.2\times 10^{7}\rm erg\,s^{-1}\,cm^{-2}$, HAT-P-11b: $1.34\times 10^{8}\,\rm erg\,s^{-1}\,cm^{-2}$, i.e., 4 times more, as time integrated over the eccentric orbit). These lead to different expected equilibrium temperatures: 650  ±   60 K for GJ 436b (T08) versus 878  ±   15 K, assuming complete heat redistribution. Also, the different stellar effective temperatures lead to different spectral distributions of the infalling flux. It will be interesting to compare the atmospheric properties of the two planets to see the effect of different integrated flux and spectral flux distributions on these otherwise similar planets.

While both GJ 436b and HAT-P-11b have similar mass and radius, and inferred structure of Z = 0.9 heavy element content with 10% H/He envelope, it is worth noting that the metallicity of the host stars, and presumably the environment these planets formed in, is different; [Fe/H] = −0.03  ±   0.20 for GJ 436 (T08) and [Fe/H] = +0.31  ±   0.05 for HAT-P-11.

From the observational point of view, HAT-P-11 is a bright star (V ∼ 9.58), more than a magnitude brighter than GJ 436. The transits of HAT-P-11b are at a small impact parameter (b = 0.347^+0.130_−0.139) as compared to the near-grazing transit of GJ 436b (b = 0.848). The total duration of the transit for HAT-P-11b is about twice (0.0957  ±   0.0012 days) that of GJ 436b (0.042 days). We can thus expect very high-quality follow-up observations of HAT-P-11b in the future.

The RV of HAT-P-11 shows, in addition to the 4.8878 day period induced by HAT-P-11b,  a significant trend of $G_1 = 0.0297 \,\,{\pm}\,\, 0.0050\rm \,\rm m\,s^{-1}\,day^{-1}$. This drift could be due to a second planet in the system, corresponding to $M_2\sin i_2/a_2^2 \sim \dot{G}_1 / G = (0.061 \,\,{\pm}\,\, 0.01) M_{\rm J}\,\mathrm AU^{-2}$, where the subscript "2" refers to HAT-P-11c, and G is the gravitational constant. Such a conjecture is scarcely surprising: for example, Bouchy et al. (2009) have noted that 16 out of the 20 hot Neptune planets discovered heretofore are members of multiple planet systems. Interestingly enough, as pointed out by Bouchy et al. (2009), this fraction (80%) is significantly higher than 23%, which is the fraction of all known exoplanets which are in multiple systems. We also recall that the significant orbital eccentricity of HAT-P-11b is reminiscent of that of GJ 436b, and that the latter has motivated the speculation that there may be another planet in the system, sufficient to maintain the eccentricity of GJ436b (Ribas et al. 2008; Bean & Seifahrt 2008). Finally, as was also pointed out by Bouchy et al. (2009), of seven multiple planet systems containing both a gas-giant planet and a hot Neptune planet, all have super-solar metallicities. The high metallicity of HAT-P-11, [Fe/H] = +0.31  ±   0.05, then strengthens the suspicion that there may be multiple planets in the system.

Continued monitoring of the RV of HAT-P-11 over the next few years should indicate whether the drift persists and whether it shows evidence of curvature of the orbit. Another way to establish the presence of a second body in the system is to search for TTVs in repeated transits (Agol et al. 2005; Holman & Murray 2005).

The fact that HAT-P-11 will lie on one of the detector chips for the Kepler spacecraft presents remarkable opportunities for scientific follow-up of the discovery presented here. If all goes well, Kepler should be capable of 1 minute cadence photometry at a precision of 0.1 mmag for this bright star. Currently, Kepler is expected to be operational for a minimum of 3.5 years, and perhaps longer; in that time it should have more than 250 transits by HAT-P-11b. In addition, it should show photometric variability due to surface inhomogeneities such as starspots, spanning more than 40 rotations if indeed the rotation period of the planet is verified to be about 29.2 days. Here we list some opportunities for important scientific follow-up with Kepler:

Kepler transit data. With a depth of 4.2 mmag, the transit light curve of HAT-P-11b should be measurable by Kepler with excellent precision for a single light curve and extraordinary precision when many light curves are combined. To maximize the accuracy of parameters derived for the HAT-P-11 system, it is important that, in addition to the Kepler observations, precise ground-based RV data be obtained through the duration of the Kepler observations.

The possibility of a second planet in the system can be tested through a search for TTVs in the large number of transits to be observed by Kepler. As noted above, RV data obtained over the same time frame could provide independent evidence for such a planet, and comparison of the two data sets would be of obvious importance.

Kepler OOT photometry. The Kepler data should reveal the detailed modulation of photospheric brightness—presumably due to starspot activity—with great precision and time coverage over the lifetime of the mission. This could yield much useful information about the star, including its rotation and even differential rotation (Fröhlich 2007). Furthermore, if multiple transits of HAT-P-11b pass in front of the same starspot on the face of the slowly rotating star, then extraordinary precision can be achieved in determining the rotation rate (Silva-Valio 2008). Comparison with concurrent ground-based measurements of the chromospheric emission S index will enable investigation of how photospheric spots correlate with chromospheric emission on stars other than the sun. As noted earlier, the data on S index reported here show a long-term (∼450 day) variation reminiscent of solar activity cycles; monitoring both photospheric emission variations from Kepler broadband photometry and chromospheric emission variations from ground based telescopes could prove to be useful.

Kepler asteroseismology. Kepler will obtain extraordinary information about the internal structure and evolution of stars from asteroseismic investigations, which measure the frequencies and amplitudes of acoustic oscillation modes ("p-modes") through small changes in the integrated brightness of the star. The "large splitting" of p-mode frequencies yields direct information on the radius of the star, and their "small splitting" gives information on the H/He ratio in the core and hence the age of the star (e.g., Christensen-Dalsgaard 2004). A great deal of additional information is possible if the star is accompanied by a transiting planet (like HAT-P-11b); in this case, the transit light curves yield completely independent information on the radius and mean density of the star. An additional bonus in the case of HAT-P-11 is that it has an excellent Hipparcos parallax. Combining these independent pieces of information should lead to very detailed knowledge of its internal structure, including second-order near-surface effects not accounted for in normal interior structure models (Kjeldsen et al. 2008). Information on the age of the star, if it can be obtained from measurement of the small splitting, could provide a useful comparison with the age inferred either from the global (light curve, RV, stellar isochrone) modeling of the data or from the Ca HK index.

As noted above, there is a suggestion from the variation of its S index that HAT-P-11 is undergoing long-term variations of its surface magnetic activity somewhat reminiscent of the solar activity cycle. On the Sun, solar cycle-related variations of solar magnetic activity are known to be correlated with frequency variations of p-modes (e.g., Woodard & Noyes 1985; Woodard et al. 1991), because the sub-surface structure is modified in regions of solar magnetic activity. Monitoring the p-mode frequencies of HAT-P-11 over the course of the Kepler mission, along with continued ground-based monitoring of the S-index, may show whether similar behavior occurs in the sub-surface layers of another star with significantly different structure.

Finally, we note that two other stars with known transiting planets, namely, TrES-2 (O'Donovan et al. 2006) and HAT-P-7 (Pál et al. 2008a), also lie on the Kepler detector array, so it should be possible to carry out studies of these stars and their planets similar to those described above for the HAT-P-11 system. It would be interesting to compare the results for these two stars, with spectral type F and G, respectively, with those for the K star HAT-P-11. Almost certainly many more TEP systems will be discovered by Kepler, but those systems that are known in advance of the mission should be especially valuable since observations with the 1 minute cadence can be scheduled from the start of the mission, plus the longer baseline improves our sensitivity to detect secular changes.

Although briefly described earlier (Bakos et al. 2007), it is worth defining the EPD technique, since it is extensively used in this work. The EPD method attempts to determine the actual functional form of E in Equation (A1). The EPD effects are treated and determined as specific for each star, i.e., no information from other stars are used (this is a key difference when compared to TFA that uses a template of other stars—see below).

Constant (or simple) EPD. The simplest form of EPD is when we assume that the underlying m₀ signal is constant, and also that 〈E〉 ≳ 〈T〉 (i.e., the effects to be corrected by EPD are of the same order as those corrected by TFA). A typical application is for planetary transits in survey mode, since such transits are short events compared to the total orbital time, and have mostly been observed around stars with much smaller variation than the transit amplitude (note that there are exceptions; see Alonso et al. 2008); thus the underlying signal can be approximated with its median ∼95% of the time. The E relation between δm(t_i) ≡ m(t_i) − m₀(t_i) = m(t_i) − 〈m(t_i)〉_m and $\vec{e_i(t_i)}$ parameters is sought via a χ² minimization, usually by nonlinear least-squares method with outlier rejection (〈〉_m denotes the median value). Upon determining E, we get an EPD corrected signal δm_EPD(t_i) = δm(t_i) − E(t_i).

Our experience with HATNet survey data has shown that the effects of E and T are indeed comparable, and thus we use constant EPD. Specifically, the external parameters are the X,Y sub-pixel position, the background and its standard deviation, the stellar profile parameters characterizing the point-spread function width (called S) and its elongation (D and K), plus the hour angle and zenith distance.

Simultaneous EPD. Used primarily in the analysis of the photometric follow-up data, where most of the observations are centered on the transits; thus, the assumption of constant signal does not hold. One possibility is to use the OOT section of the light curve to determine the EPD parameters in a constant EPD mode, and then apply the correction for the in-transit section of the light curve as well. This is sub-optimal, as the OOT may be very short compared to the in-transit section, or may be missing. Thus, the key difference compared to constant EPD is that we use all data to recover the dependence on external parameters, i.e., $m_0(\vec{p},t_i)$ in Equation (A1) is not constant, but assumed to be a function of other parameters and time. The fitting procedure determines both the E correlation and the optimal set of $\vec{p}$ parameters simultaneously.

The signal after (or without) the EPD procedure still contains the general systematic variations denoted as T(t_i) in Equation (A1). These are suppressed by the Trend Filtering Algorithm (TFA; see Kovács et al. 2005), assuming that certain other stars in an M element light-curve "template" show similar variations.

To recall, for selected star j, TFA minimizes the following expression:

$$\begin{equation} D = \sum _{i=1}^{N} \big\lgroup m_j(t_i) - m_{j,0}(\vec{p},t_i) - F_M(t_i) - E_j(\vec{e}(t_i)) \big\rgroup ^2, \end{equation} \tag{ A2 }$$

where

$$\begin{equation} F_M(t_i) = \sum _{k=1,k\ne j}^{k=M}{c_k m_k(t_i)} \end{equation} \tag{ A3 }$$

is the TFA filter function of M elements. The notation is like above in Equation (A1), $m_{j,0}(\vec{p},t_i)$ is the model function, and E_j is the EPD correction.

There are many variants of TFA, and they have been used widely in this work. The signal search in the HATNet data is done in two parallel steps. First via simple (or constant) TFA where there is no assumption on the periodicity, the model function m_j,0 is a constant, E_j is zero, as EPD has been performed as an independent step before TFA, and the c_k coefficients are sought. Second, via reconstructive TFA, where the simple TFA is followed by a frequency search, the signal is phase-folded with the most significant frequency, the model function m_j,0 is fitted to the folded data, the model function is unwrapped to the original time base, and D is minimized again to iteratively determine the c_k coefficients. The third method, introduced in this work, is simultaneous TFA, whereby the c_k TFA coefficients and the functional dependence of $m_{j}(\vec{p},t_i)$ on $\vec{p}$ parameters are simultaneously determined.

EPD and TFA can be performed sequentially, or even simultaneously. In the analysis of the photometric follow-up data (see Section 2.3), we implemented such a simultaneous EPD–TFA, where the $E_j(\vec{e})$ EPD function is simultaneously fitted with the c_k TFA template coefficients and the m_j model function.

Both TFA and EPD can be performed globally, using one set of coefficients for the entire data set, and locally, when Equation (A2) is split up into smaller data blocks, such as one-night segments, and the c_k TFA coefficients and the EPD function parameters are fitted for each data block separately to allow for changing systematics. Conversely, for global TFA systematic variations of template stars have to "match" those of the main target (minus the model) for all nights.

Finally, we note that if we do not use EPD (eliminate the term for E in Equation (A1)), but apply TFA only, then naturally, TFA will take care of some of the systematics that otherwise would have been corrected by EPD.