The Transit Light Curve Project. XI. Submillimagnitude Photometry of Two Transits of the Bloated Planet WASP-4b1 Joshua N. Winn2, Matthew J. Holman3, Joshua A. Carter2, 9 0 Guillermo Torres3, David J. Osip4, Thomas Beatty2 0 2 n ABSTRACT a J 7 2 We present photometry of two transits of the giant planet WASP-4b with ] a photometric precision of 400–800 parts per million and a time sampling of P 25–40 s. The two midtransit times are determined to within 6 s. Together with E . previouslypublished times, thedataareconsistent withaconstant orbitalperiod, h p giving no compelling evidence for period variations that would be produced by a - o satellite or additional planets. Analysis of the new photometry, in combination r t with stellar-evolutionary modeling, gives a planetary mass and radius of 1.237 s ± a 0.064 M and 1.365 0.021 R . The planet is 15% larger than expected based Jup Jup [ ± on previously published models of solar-composition giant planets. With data of 1 v the quality presented here, the detection of transits of a “super-Earth” of radius 6 1.75 R would have been possible. 4 ⊕ 3 4 Subjectheadings: planetarysystems—stars:individual(WASP-4=USNO-B1.00479- . 1 0948995) 0 9 0 : v i 1. Introduction X r a Wilson et al. (2008) recently reported the discovery of WASP-4b, a giant planet that orbits and transits a G7V star with a period of 1.34 days. This discovery is notable because the planet has an unusually large radius and short orbital period, and the star is one of the 1Basedonobservationswith the 6.5mMagellanTelescopeslocatedatLasCampanasObservatory,Chile. 2DepartmentofPhysics,andKavliInstituteforAstrophysicsandSpaceResearch,MassachusettsInstitute of Technology, Cambridge, MA 02139,USA 3Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138,USA 4Las Campanas Observatory,Carnegie Observatories, Casilla 601, La Serena, Chile – 2 – brightest transit hosts that is known in the Southern sky (V = 12.5). The large radius seems to place the planet among the “bloated” planets for which there is no clear explanation (see, e.g., Burrows et al. 2007, Guillot 2008). The short period raises the possibility of observing tidal decay (Rasio et al. 1996, Sasselov 2003) and makes WASP-4b an attractive target for observationsofoccultations(secondaryeclipses) thatwouldleadtodetectionsofthereflected light and thermal emission from the planet’s atmosphere. The host star’s brightness and southern declination are important for a practical rea- son: they allow the large-aperture telescopes of the southern hemisphere to be used ad- vantageously. In this paper, we report on observations of a transit of WASP-4b with the Magellan/Baade 6.5m telescope, with the goal of deriving independent and refined param- eters for this interesting system. Previous papers in this series, the Transit Light Curve project, have achieved this goal by combining the information from many independent tran- sit observations with smaller telescopes (see, e.g., Holman et al. 2006, Winn et al. 2007). In principle, with a larger telescope, it should be possible to achieve this goal with fewer observations and also to measure precise midtransit times, which can be used to search for additional planets via the method of Holman & Murray (2005) and Agol et al. (2005). Gillon et al. (2008) recently presented photometry of WASP-4 with one of the 8.2m Very Large Telescopes, with the same motivation. This paper is organizedas follows. 2 describes theobservations anddata reduction, 3 § § describes the photometric analysis, 4 describes the results of stellar-evolutionary modeling, § and 5 discusses the newly-measured midtransit times and a refined ephemeris. Finally, 6 § § discusses the refined measurement of the planetary radius, and considers how small a planet we could have detected, given data of the quality presented here. 2. Observations and Data Reduction We observed the transits of UT 2008 Aug 19 and 2008 Oct 09 with the Baade 6.5m telescope, one of the two Magellan telescopes at Las Campanas Observatory in Chile. We used the Raymond and Beverly Sackler Magellan Instant Camera (MagIC) and its SITe 2048 2048 pixel CCD detector, with a scale of 0′.′069 pixel−1. Ordinarily this detector uses × four amplifiers, each of which reads a quadrant of the array, giving a total readout time of 23 s. We used a 2048 256 pixel subarray and a single amplifier, giving a readout time of × 10 s. We rotated the field of view to align WASP-4 and a nearby comparison star along the long axis of the subarray (parallel to the read register). The comparison star is 36′′ east and 71′′ south of WASP-4. At the start of each night, we verified that the time stamps recorded by MagIC were in agreement with GPS-based times to within one second. On each night – 3 – we obtained repeated z-band exposures of WASP-4 and the comparison star for about 5 hr bracketing the predicted transit time. Autoguiding kept the image registration constant to within 10 pixels over the course of the night. During the 2008 Aug 19 observations, WASP-4 rose from an airmass of 1.09 to 1.03, and then set to an airmass of 1.23. At first, we used an exposure time of 30 s. Shortly after midtransit the seeing improved abruptly, from a full-width at half-maximum (FWHM) of 11 pixels to 7 pixels. As a result of the higher rate of detected photons per pixel, some images were spoiled due to nonlinearity and saturation. For the rest of the night we used an exposure time of t = 15 s. exp During the 2008 Oct 09 observations, WASP-4 rose from an airmass of 1.10 to 1.03, and then set to an airmass of 1.21. Having learned our lesson in August, we deliberately defocused thetelescope so thatthe imagewidthwas dominated bythe effects ofthetelescope aberration rather than natural seeing. This was done by moving the secondary mirror by a constant displacement relative to the in-focus position determined by the guider probe. The exposure time was t = 30 s. The stellar images were “donuts” with a diameter of exp approximately 25 pixels. We used standard IRAF procedures for overscan correction, trimming, bias subtraction, and flat-field division. The bias frame was calculated from the median of 60 zero-second exposures, and the flat field for each night was calculated from the median of 60 z-band exposures of a dome flat screen. We performed aperture photometry of WASP-4 and the comparison star and divided the flux of WASP-4 by the flux of the comparison star. We experimented withdifferent aperturesizes andsky regions, aiming to minimize the variations in the out-of-transit (OOT) portion of the differential light curve. Best results were obtained with an aperture radius of 38 pixels for the 2008 Aug 19 observations and 35 pixels for the 2008 Oct 09 observations. A few images were also obtained with Johnson-Cousins VRI filters to measure the difference in color between WASP-4 and the comparison star. The instrumental magnitude differences (target minus WASP-4) were ∆V = 0.434, ∆R = 0.469, ∆I = 0.478, and ∆z = 0.472. Evidently the two stars are similar in color, with the comparison star being slightly bluer [∆(V I) = 0.044]. As described in 3, the z-band time series was corrected for − − § differential extinction between the target and comparison star by fitting a linear function of airmass to the magnitude difference. The results were consistent with the expectation that the bluer comparison star suffers from greater extinction per unit airmass. The extinction- corrected data are given in Table 1, and plotted in Fig. 1, along with the best-fitting model. Due to the abrupt seeing change on 2008 Aug 19 and the associated change in exposure – 4 – Fig. 1.— Top.—Relative z-band photometry of WASP-4 based on observations with the Magellan (Baade) 6.5m telescope. Middle.—Composite light curve. The solid line shows the best-fitting model. Bottom.—Residuals between the data and the best-fitting model. – 5 – level, we consider the 15 s and 30 s exposures as two separate time series (TS1 and TS2). Thus, together with the 2008 Oct 09 time series (TS3), there were 3 time series to be analyzed. The TS1 pre-ingress data has a median time sampling of 41 s and a standard deviation of 478 parts per million (ppm). The TS2 post-egress data has a median time sampling of 26 s and a standard deviation of 691 ppm. These noise levels are about 13% larger than the calculated noise due to photon-counting statistics, read noise, sky noise, and scintillation noise (using the approximate formulas of Reiger 1963 and Young 1967). The ratio of the observed TS1 noise to the observed TS2 noise is 0.69, which is nearly equal to (t /t )−1/2 = 0.71. The near-equality is evidence that the dominant noise source is a exp,2 exp,1 −1/2 combination of photon noise and scintillation noise, both of which vary as t . In the final exp version of TS3, which has a median time sampling of 41 s, the pre-ingress data has an rms of 475 ppm and the post-egress data has an rms of 488 ppm. These values are about 17% above the calculated noise level. In neither light curve did we detect signficant correlation between the noise and the pixel coordinates, FWHM, or shape parameters of the stellar images. 3. Photometric Analysis We fitted a transit light curve model to the data, based on the analytic formulae of Mandel & Agol (2002). The set of model parameters included the planet-to-star radius ratio (R /R ), the ratio of the stellar radius to the orbital semimajor axis (R /a), the orbital p ⋆ ⋆ inclination (i), and the midtransit time (T ). We also fitted for three parameters (∆m , k c 0 z and k ) specifying a correction for systematic errors, t ∆m = ∆m +∆m +k z +k t, (1) cor obs 0 z t where z is the airmass, t is the time, ∆m is the observed magnitude difference between obs the target and comparison star, and ∆m is the corrected magnitude difference that is cor compared to the idealized transit model. The k term specifies the differential airmass z extinction correction that was mentioned in 2. The k term was not strictly necessary t § to fit the data (it was found to be consistent with zero), but we included it in order to derive conservative error estimates for the transit times, as the parameter k is covariant t with the transit time. We allowed ∆m , k , and k to be specific to each time series, with 0 t z the exception that TS1 and TS2 had a common value of k because those data were obtained z on the same night. The final two model parameters were the coefficients (u and u ) of a quadratic limb- 1 2 darkening law, I µ = 1 u (1 µ) u (1 µ)2, (2) 1 2 I − − − − 1 – 6 – where µ is the cosine of the angle between the line of sight and the normal to the stellar photosphere, and I is the specific intensity. We allowed u and u to vary freely subject µ 1 2 only to the conditions u + u < 1, u + u > 0, and u > 0. It proved advantageous to 1 2 1 2 1 perform the fit using the linear combinations v = u +2.33 u , v = u 2.33u , (3) 1 1 2 2 1 2 − because v and v have nearly uncorrelated errors (for further discussion, see P´al 2008). 1 2 We assumed the orbit to be circular because the RV data of Wilson et al. (2008) are consistent withacircular orbit, andbecause theexpected timescale fororbitalcircularization at present, 5 4 P M a p τ = Q , (4) c p 63 (cid:18)2π(cid:19)(cid:18)M (cid:19)(cid:18)R (cid:19) ⋆ p (Rasio et al. 1996) is 0.3 Myr (Q /105) for WASP-4b, which is much shorter than the p estimated main-sequence age of the star. In this expression, Q is the the tidal dissipation p parameter (see, e.g., Goldreich & Soter 1966). This order-of-magnitude argument suggests that assuming a circular orbit is reasonable, although the value of Q for irradiated giant p planets is highly uncertain, and the expected timescale is highly approximate because it ignores the coupled evolution of the orbital distance and eccentricity (Jackson, Greenberg, & Barnes 2008). The fitting statistic was N 2 f (obs) f (calc) χ2 = i − i , (5) F (cid:20) σ (cid:21) Xi=1 f,i where N is the number of flux measurements (photometric data points), f (obs) is the ith i measurement, f (calc) is the calculated flux given a particular choice of model parameters, i and σ is the uncertainty in the ith measured flux. We determined appropriate values of f,i σ as follows. First, we multiplied the calculated errors in each time series by a constant f,i chosen to give a minimum value of χ2/N = 1. The constants were 1.20, 1.12, and 1.21 F dof for TS1, TS2, and TS3 respectively. Next, we assessed the time-correlated noise (also called “red noise”) by examining the autocorrelation function, the power spectrum, and a plot of the Allan (1966) deviation of the residuals. We also used the method described by Winn et al. (2008), in which the ratio β is computed between the standard deviation of time- averaged residuals, and the standard deviation one would expect assuming white noise. For TS1 and TS2 we found no evidence for significant correlations. For TS3 we found structure in the autocorrelation function on a time scale of 15-20 min, the approximate ingress or egress duration, giving β = 1.52. One naturally suspects that the correlated noise – 7 – represents measurement error, although it is also possible that the noise is astrophysical, arising from starspots or other stellar inhomogeneities. In support of an astrophysical origin, we find no evidence for correlated noise (β = 1) when considering only the out-of-transit data. Nevertheless we cannot draw a firm conclusion; instead we attempt to account for the correlations by multiplying the error bars of TS 3 by an additional factor of β = 1.52. Thus, for TS1, TS2, and TS3, the final values of σ were equal to the calculated error bars f,i multiplied by 1.20, 1.12, and 1.21 1.52 = 1.84, respectively. The error bars given in Table 1 × are the final values of σ that were used in the fitting process. f,i To determine the allowed ranges of each parameter, we used a Markov Chain Monte Carlo (MCMC) technique, with the Gibbs sampler and the Metropolis-Hastings algorithm, to estimate the a posteriori joint probability distribution of all the model parameters. This algorithm creates a sequence of points (a “chain”) in parameter space by iterating a jump function, which in our case was the addition of a Gaussian random deviate to a randomly- selected single parameter. After this operation, if the new point has a lower χ2 than the F previous point, the “jump” is executed: the new point is added to the chain. If not, then the jump is executed with probability exp( ∆χ2/2). When the jump is not executed, the − F current point is repeated in the chain. The sizes of the random deviates are adjusted so that 40% of jumps are executed. After creating 10 chains of 500,000 links to check for ∼ mutual convergence, and trimming the first 20% of the links to eliminate artifacts of the initial conditions, the density of the chains’ points in parameter space was taken to be the joint a posteriori probability distribution of the parameter values. Probability distributions for individual parameters are created by marginalizing over all other parameters. The results are given in Table 2, which gives the median of each distribution, along with the 68.3% lower and upper confidence limits (defined by the 15.85% and 84.15% levels of the cumulative distribution). The entries designated A are those that follow directly from the photometric analysis. The entries designated B are those that are drawn from Gillon et al. (2008) and are repeated here for convenience. The entries designated C are based on a synthesis of our modeling results and theoretical models of stellar evolution, as discussed in 4. The entries designated D are the parameters of a refined transit ephemeris based on the § two newly-measured transit times as well as some other available timing data (see 5). § As a consistency check we tried fitting the 2008 Aug 19 data and the 2008 Oct 09 data separately. We found that the results for the parameters R /R , R /a, i, u , and u were all p ⋆ ⋆ 1 2 in agreement within 1σ, suggesting that our error estimates are reasonable. – 8 – 4. Theoretical isochrone fitting The combination of transit photometry and the spectroscopic orbit (radial-velocity vari- ation) of the host star do not uniquely determine the masses and radii of the planet and the 2/3 1/3 star. There remain fitting degeneracies M M and R R M (see, e.g., Winn p ⋆ p ⋆ ⋆ ∝ ∝ ∝ 2008). Webroke these degeneracies byrequiring consistency between theobserved properties of the star, the stellar mean density ρ that can be derived from the photometric parameter ⋆ a/R (Seager & Mallen-Ornelas 2003, Sozzetti et al. 2007), and theoretical models of stellar ⋆ evolution. The inputs were T = 5500 100 K and [Fe/H] = 0.03 0.09 from Gillon et eff ± − ± al. (2008), the stellar mean density ρ = 1.694+0.017 g cm−3 derived from the results for the ⋆ −0.037 a/R parameter, and the Yonsei-Yale (Y2) stellar evolution models by Yi et al. (2001) and ⋆ Demarque et al. (2004). We computed isochrones for the allowed range of metallicities, and for stellar ages ranging from 0.1 to 14 Gyr. For each stellar property (mass, radius, and age), we took a weighted average of the points on each isochrone, in which the weights were proportional to exp( χ2/2) with − ⋆ 2 2 2 ∆[Fe/H] ∆T ∆ρ χ2 = + eff + ⋆ . (6) ⋆ (cid:20) σ (cid:21) (cid:20) σ (cid:21) (cid:20) σ (cid:21) [Fe/H] Teff ρ⋆ Here, the ∆ quantities denote the deviations between the observed and calculated values at each point. The asymmetric error bar in ρ was taken into account by using different ⋆ values of σ depending on the sign of the deviation. The weights were further multiplied by ρ⋆ a factor taking into account the number density of stars along each isochrone, assuming a Salpeter mass function. This procedure is essentially the same as that employed by Torres et al. (2008). The only difference is that we calculated the 68.3% uncertainties by assuming that the errors in T , ρ , and [Fe/H] obey a Gaussian distribution, while Torres et al. (2008) eff ⋆ took the distribution to be uniform within the quoted 1σ limits. Throughthisanalysis, wefoundM = 0.925 0.040M andR = 0.912 0.013R . The ⋆ ⊙ ⋆ ⊙ ± ± stellar agewaspoorlyconstrained, withaformallyallowed rangeof6.5 2.3Gyr andanearly ± uniform distribution. The corresponding planetary mass and radius were obtained by merg- ing the results for the stellar properties with the parameters determined in our photometric analysis, andwiththestellarradial-velocitysemiamplitudeK = 0.24 0.01kms−1 measured ⋆ ± byWilsonetal.(2008). TheresultsareM = 1.237 0.064M andR = 1.365 0.021R . p Jup p Jup ± ± Table 2 gives these results, along with the values for some other interesting parameters that can be derived from the preceding results. As a consistency check, we computed the implied stellar surface gravity and its uncertainty based on our analysis, finding logg = 4.481 0.008 ± where g is in cm s−2. This agrees with the spectroscopic determination of surface gravity, logg = 4.3 0.2, made by Wilson et al. (2008) based on an analysis of the widths of ± pressure-sensitive lines in the optical spectrum. – 9 – It is important to keep in mind that the quoted error bars for the parameters designated C in Table 2 are based on the measurement errors only, and assume that the any systematic errors in the Y2 isochrones are negligible. As a limited test for the presence of such errors, Torres et al. (2008) tried analyzing transit data for 9 systems using isochrones computed by three different groups: the Y2 isochrones used here as well as those of Girardi et al. (2000) and Baraffe et al. (1998). They found the differences to be smaller than the error bars, especially for stars similar to the Sun such as WASP-4. Similar results were found by Southworth (2009). Nevertheless, the true systematic errors in the Y2 isochrones are not known, and we have not attempted to quantify them here, although it seems plausible that the masses and radii are subject to an additional error of a few percent. 5. Transit Times and a Refined Ephemeris Precise measurements oftransit times areimportant because thegravitationalperturba- tions from other bodies in the system (such as a satellite or additional planet) could produce detectable variations in the orbital period. Based on the MCMC analysis described in the previous section, the uncertainties in the two transit times are 4.6 s and 4.9 s, making them among the most precise such measurements that have been achieved. This is a consequence of good photometric precision and fine time sampling, along with the relative insignificance of time-correlated noise and the large transit depth. Asacheckontheerrorbars,wealsoestimatedthemidtransittimeanditserrorusingthe residual permutation (RP) method, atype ofbootstrapanalysis that attempts toaccount for time-correlated errors. Fake data sets are created by subtracting the best-fitting model from the data, then adding the residuals back to the model after performing a cyclic permutation of their time indices. For each fake data set, χ2 is minimized as a function of T , ∆m , k F c 0 z and k (the other parameters being held constant, as they are nearly uncorrelated with T ). t c The distribution of the results for T is taken to be the a posteriori probability distribution c for T . For the 2008 Aug 19 data, the RP-based error bar was only 4% larger than the c MCMC-based error bar; the results were nearly indistinguishable. For the 2008 Oct 09 data, the RP-based error was 25% larger than the MCMC-based error. To be conservative, we report the larger RP-based errors in Table 1 and used those larger error bars in recomputing the ephemeris (see below). As a further check on the analysis, we allowed TS1 and TS2 to have independent values of T , ∆m , k , and k during the fitting process. In other words we fitted all the data c 0 z t but did not require that TS1 and TS2 agree on the transit time. The result was that the difference between T (TS1) and T (TS2) was 2.3 13.4 s. We also tried a similar experiment c c ± – 10 – with the 2008 Aug 19 data, splitting it into two equal parts that were fitted jointly with TS1 and TS2 (which in this case were required to agree on the midtransit time so as to provide a constraint on the transit duration). The result in this case was ∆T = 4.6 11.9 s. The c ± mutual consistency of the results suggests that our error bars are reasonable. The transits of 2008 Aug 19 and 2008 Oct 09 were separated by 38 orbital periods. By calculating ∆T /38we derive anindependent estimate oftheorbitalperiod, P = 1.3382369 c ± 0.0000024 days. The most precise determination previously reported was P = 1.3382324 ± 0.0000029 days, by Gillon et al. (2008), based on 2 years of data. Our period has comparable precision, though it is based on only 2 transits separated by 50 days. The difference between the two independent period determinations is 0.39 0.32 s. ± We fitted a linear function of epoch to the two newly-measured midtransit times along with the 5 transit times given in Table 2 of Gillon et al. (2008). The fit had χ2=7.8 with 5 degrees of freedom. The chance of finding a value of χ2 this large by chance is 17%, using the quoted 1σ error bars and assuming the errors obey a Gaussian distribution. We deem this an acceptable fit, and conclude that the available data do not provide compelling evidence for any departures froma constant period. The refined transit ephemeris is T (E) = T (0)+EP, c c with T (0) = 2,454,697.797562 0.000043[BJD], c ± P = 1.33823214 0.00000071. (7) ± Fig. 2 shows a plot of the differences between the observed and calculated transit times. 6. Discussion Our results for the planetary and stellar properties are in accord with the previous analyses by Wilson et al. (2008) and Gillon et al. (2008). In general our error bars are comparable in size or smaller than those of Gillon et al. (2008), who observed a transit with one of the Very Large Telescopes. This consistency is reassuring, especially since our error estimates are more conservative in some respects. We have accounted for uncertainty in the limb-darkening coefficients, as well as the slopes of systematic trends with time and airmass. Previous investigators assumed that these parameters were known exactly, leading to underestimated errors in any covariant parameters. Southworth (2008) demonstrated this effect for the limb-darkening coefficients in particular.1 In addition, Gillon et al. (2008) 1Our results for the limb-darkeningcoefficients areu1 =0.311 0.041andu2 =0.227 0.089. These are ± ± not too far from the tabulated values 0.265 and 0.303 given by Claret (2004) for a star with the observed