HAT-P-20b, HAT-P-21b, HAT-P-22b, and HAT-P-23b

# HAT-P-20b–HAT-P-23b: Four massive transiting extrasolar planets1

## Abstract

We report the discovery of four relatively massive (2–7 ) transiting extrasolar planets. HAT-P-20b orbits the moderately bright V=11.339 K3 dwarf star GSC 1910-00239 on a circular orbit, with a period  d, transit epoch (), and transit duration  d. The host star has a mass of , radius of , effective temperature  K, and metallicity . The planetary companion has a mass of , and radius of  yielding a mean density of , which is the second highest value among all known exoplanets. HAT-P-21b orbits the V=11.685 G3 dwarf star GSC 3013-01229 on an eccentric ( orbit, with a period  d, transit epoch , and transit duration  d. The host star has a mass of , radius of , effective temperature  K, and metallicity . The planetary companion has a mass of , and radius of  yielding a mean density of . HAT-P-21b is a border-line object between the pM and pL class planets, and the transits occur near apastron. HAT-P-22b orbits the bright V=9.732 G5 dwarf star HD 233731 on a circular orbit, with a period  d, transit epoch , and transit duration  d. The host star has a mass of , radius of , effective temperature  K, and metallicity . The planet has a mass of , and compact radius of  yielding a mean density of . The host star also harbors an M-dwarf companion at a wide separation. Finally, HAT-P-23b orbits the V=12.432 G0 dwarf star GSC 1632-01396 on a close to circular orbit, with a period  d, transit epoch , and transit duration  d. The host star has a mass of , radius of , effective temperature  K, and metallicity . The planetary companion has a mass of , and radius of  yielding a mean density of . HAT-P-23b is an inflated and massive hot Jupiter on a very short period orbit, and has one of the shortest characteristic in-fall times ( Myr) before it gets engulfed by the star.

planetary systems — stars: individual ( HAT-P-20, GSC 1910-00239, HAT-P-21, GSC 3013-01229, HAT-P-22, HD 233731, HAT-P-23, GSC 1632-01396 ) techniques: spectroscopic, photometric

## 1. Introduction

The majority of the known transiting extrasolar planets (TEPs) have been found to lie in the 0.5  to 2.0  mass range. The apparent drop in their mass distribution at 2  has been noted by, e.g., Southworth et al. (2009), and by Torres et al. (2010). In the currently known sample, 75% of the TEPs have planetary mass , and there appears to be a minor peak in their occurrence rate at , which then sharply falls off towards higher masses. Are there any biases present against discovering massive planets? Such planets tend to be less inflated, and theory dictates that their radii shrink as their mass increases towards the brown dwarf regime. According to Baraffe, Chabrier, & Barman (2010), this reversal of the relation happens around , and falls off as (see e.g. Fortney et al., 2009). The smaller radii for massive planets yield a minor bias against discovering them via the transit method, since they produce shallower transits. Very massive planets can induce stellar variability of their host stars (Shkolnik et al., 2009), somewhat decreasing the efficiency of detecting their shallow transits via simple algorithms that expect constant out-of-transit light curves. Also, the host stars of massive planets are typically more rapid rotators: the average for host stars with planets is 3.9  (with standard deviation around the mean), whereas the same values for the massive planet hosts stars, are (with standard deviation around the mean)27. The five fastest rotators all harbor planets more massive than 2 . This presents a bias against discovering them either via radial velocity (RV) searches, which are more efficient around quiet non-rotating dwarfs, or via transit searches, where the targets may be discarded during the confirmation phase. Along the same lines, the large RV amplitude of the host star, as caused by the planetary companion, may even lead to erroneous rejection during the reconnaissance phase of candidate confirmation, since such systems resemble eclipsing binaries. Finally, there is a tendency that massive planets are more likely to be eccentric28 (Southworth et al., 2009), meaning that they require more RV observations for proper mapping of their orbits, and thus leading to a slower announcement rate. On the other hand, a strong bias for detecting such planets–compensating for most of the effects above–is the fact that the large RV amplitudes of the host stars are easier to detect, since they do not require internal precisions at the  level (see HAT-P-2, where valuable data was contributed to the RV fit by modest precision instruments yielding precision; Bakos et al. 2007). Altogether, while there are minor biases for and against detecting massive transiting planets, their overall effect appears to be negligible, and the drop in frequency at seems to be real.

Massive planets are important for many reasons. They provide very strong constraints on formation and migration theories, which need to explain the observed distribution of planetary system parameters in a wide range (Baraffe et al., 2008; Baraffe, Chabrier, & Barman, 2010), from 0.01  (Corot-7b; Queloz et al., 2009) to 26.4  (Corot-3b; Deleuil et al., 2008). Heavy mass objects necessitate the inclusion of other physical mechanisms for the formation and migration, such as planet-planet scattering (Chatterjee et al., 2008; Ford & Rasio, 2008), and the Kozai-mechanism (Fabrycky & Tremaine, 2007). They are border-line objects between planets and brown-dwarfs, and help us understand how these populations differ and overlap (see Leconte et al., 2009, for a review). For example, a traditional definition of planets is that they have no Deuterium burning, where the Deuterium burning limit is thought to be around 13 . However, there are large uncertainties on this limit due to the numerous model parameters and solutions, and the fact that Deuterium may be able to burn in the H/He layers above the core (Baraffe et al., 2008). Another possible definition of planets is based on their formation scenario, i.e. they are formed by accretion in a protoplanetary disk around their young host star, as opposed to the gravitational collapse of a molecular cloud (brown dwarfs).

Perhaps related to the formation and migration mechanisms, a number of interesting correlations involving massive planets have been pointed out. Udry et al. (2002) noted that short period massive planets are predominantly found in binary stellar systems. Southworth et al. (2009) noted that only 8.6% of the low mass planets show significantly eccentric orbits, whereas 77% of the massive planets have eccentric orbits (although low-mass systems have lower S/N RV curves, rendering the detection of eccentric orbits more difficult). Curiously, there appears to be a lack of correlation between planetary mass and host star metallicity, while one would naively think that the formation of high mass planets (via core accretion) would require higher metal content. Until this work, there was a hint of a correlation between planetary and stellar mass (e.g. Deleuil et al., 2008), in the sense that the most massive planets orbited stars, and there was a (biased) tendency that lower mass planets orbit less massive stars.

All of these observations suffer from small-number statistics and heavy biases. One way of improving our knowledge is to expand the sample of well-characterized planets. In this work we report on 4 new massive transiting planets around bright stars, namely GSC 1910-00239, GSC 3013-01229, HD 233731, and GSC 1632-01396. This extends the currently known sample of bright () and massive () transiting planets by 30% (from 13 to 17). These discoveries were made by the Hungarian-made Automated Telescope Network (HATNet; Bakos et al., 2004) survey. HATNet has been one of the main contributors to the discovery of TEPs, among others such as the ground-based SuperWASP (Pollacco et al., 2006), TrES (Alonso et al., 2004) and XO projects (McCullough et al., 2005), and space-borne searches such as CoRoT (Baglin et al., 2006) and Kepler (Borucki et al., 2010). In operation since 2003, HATNet has now covered approximately 14% of the sky, searching for TEPs around bright stars (). We operate six wide-field instruments: four at the Fred Lawrence Whipple Observatory (FLWO) in Arizona, and two on the roof of the hangar servicing the Smithsonian Astrophysical Observatory’s Submillimeter Array, in Hawaii.

The layout of the paper is as follows. In Section 2 we report the detections of the photometric signals and the follow-up spectroscopic and photometric observations for each of the planets. In Section 3 we describe the analysis of the data, beginning with the determination of the stellar parameters, continuing with a discussion of the methods used to rule out nonplanetary, false positive scenarios which could mimic the photometric and spectroscopic observations, and finishing with a description of our global modeling of the photometry and radial velocities. Our findings are discussed in Section 4.

## 2. Observations

### 2.1. Photometric detection

Table 1 summarizes the HATNet discovery observations of each new planetary system. The calibration of the HATNet frames was carried out using standard procedures correcting for the CCD bias, dark-current and flatfield structure. The calibrated images were then subjected to star detection and astrometry, as described in Pál & Bakos (2006). Aperture photometry was performed on each image at the stellar centroids derived from the Two Micron All Sky Survey (2MASS; Skrutskie et al., 2006) catalog and the individual astrometric solutions. For certain datasets (HAT-P-20, HAT-P-22, HAT-P-23) we also carried out an image subtraction (Alard, 2000) based photometric reduction using discrete kernels (Bramich, 2008), as described in Pál (2009b). The resulting light curves were decorrelated (cleaned of trends) using the External Parameter Decorrelation (EPD; see Bakos et al., 2010) technique in “constant” mode and the Trend Filtering Algorithm (TFA; see Kovács et al., 2005). The light curves were searched for periodic box-shaped signals using the Box Least-Squares (BLS; see Kovács et al., 2002) method. We detected significant signals in the light curves of the stars as summarized below:

• HAT-P-20 – GSC 1910-00239 (also known as 2MASS 07273995+2420118; , ; J2000; V=11.339). A signal was detected for this star with an apparent depth of  mmag, and a period of  days (see Figure 1). Note that the depth was attenuated by the presence of the fainter neighbor star that is not resolved on the coarse resolution (14) HATNet pixels. Also, the depth by fitting a trapese instead of the correct Mandel & Agol (2002) model, is somewhat shallower than the maximum depth in the Mandel & Agol (2002) model fit (which was 19.6 mmag; see later in § 3.3). The drop in brightness had a first-to-last-contact duration, relative to the total period, of , corresponding to a total duration of  hr. HAT-P-20 has a red companion (2MASS 07273963+2420171, ) at 6.86″ separation that is fainter than HAT-P-20 by  mag.

• HAT-P-21 – GSC 3013-01229 (also known as 2MASS 11250598+4101406; , ; J2000; V=11.685). A signal was detected for this star with an apparent depth of  mmag, and a period of  days (see Figure 2). The drop in brightness had a first-to-last-contact duration, relative to the total period, of , corresponding to a total duration of  hr.

• HAT-P-22 – HD 233731 (also known as GSC 03441-00925 and 2MASS 10224361+5007420; , ; J2000; V=9.732). A signal was detected for this star with an apparent depth of  mmag, and a period of  days (see Figure 3). The drop in brightness had a first-to-last-contact duration, relative to the total period, of , corresponding to a total duration of  hr. HAT-P-22 has a close red companion star (2MASS 10224397+5007504, ) at 9.1″ separation and magnitude fainter.

• HAT-P-23 – GSC 1632-01396 (also known as 2MASS 20242972+1645437; , ; J2000; V=12.432). A signal was detected for this star with an apparent depth of  mmag, and a period of  days (see Figure 4). Similarly to HAT-P-20, the depth was attenuated by close-by faint neighbors. The drop in brightness had a first-to-last-contact duration, relative to the total period, of , corresponding to a total duration of  hr.

### 2.2. Reconnaissance Spectroscopy

As is routine in the HATNet project, all candidates are subjected to careful scrutiny before investing valuable time on large telescopes. This includes spectroscopic observations at relatively modest facilities to establish whether the transit-like feature in the light curve of a candidate might be due to astrophysical phenomena other than a planet transiting a star. Many of these false positives are associated with large radial-velocity variations in the star (tens of ) that are easily recognized. The reconnaissance spectroscopic observations and results for each system are summarized in Table 2; below we provide a brief description of the instruments used, the data reduction, and the analysis procedure.

One of the tools we have used for this purpose is the Harvard-Smithsonian Center for Astrophysics (CfA) Digital Speedometer (DS; Latham, 1992), an echelle spectrograph mounted on the FLWO 1.5 m telescope. This instrument delivers high-resolution spectra () over a single order centered on the Mg i b triplet (5187 Å), with typically low signal-to-noise (S/N) ratios that are nevertheless sufficient to derive radial velocities (RVs) with moderate precisions of 0.5–1.0  for slowly rotating stars. The same spectra can be used to estimate the effective temperature, surface gravity, and projected rotational velocity of the host star, as described by Torres et al. (2002). With this facility we are able to reject many types of false positives, such as F dwarfs orbited by M dwarfs, grazing eclipsing binaries, or triple or quadruple star systems. Additional tests are performed with other spectroscopic material described in the next section.

Another of the tools we have used for this purpose is the FIbre-fed Échelle Spectrograph (FIES) at the 2.5 m Nordic Optical Telescope (NOT) at La Palma, Spain (Djupvik & Andersen, 2010). We used the medium-resolution fiber which produces spectra at a resolution of and a wavelength coverage of  3600-7400 Å to observe HAT-P-20. The spectrum was extracted and analyzed to measure the radial velocity, effective temperature, surface gravity, and projected rotation velocity of the host star, following the procedures described by Buchhave et al. (2010).

Based on the observations summarized in Table 2 we find that HAT-P-21, HAT-P-22 and HAT-P-23 have rms residuals consistent with no detectable RV variation within the precision of the measurements. Curiously, HAT-P-20 showed significant RV variations, even at the modest () precision of the Digital Speedometer, and the reconnaissance RV variations (including the FIES spectrum; see later) phased up with the photometric ephemeris. All spectra were single-lined, i.e., there is no evidence that any of these targets consist of more than one star. The gravities for all of the stars indicate that they are dwarfs.

### 2.3. High resolution, high S/N spectroscopy

We proceeded with the follow-up of each candidate by obtaining high-resolution, high-S/N spectra to characterize the RV variations, and to refine the determination of the stellar parameters. These observations are summarized in Table 3. The RV measurements and uncertainties are given in Table 4, Table 5, Table 6 and Table 7 for HAT-P-20 through HAT-P-23, respectively. The period-folded data, along with our best fits described below in Section 3, are displayed in Figure 6 through Figure 9 for HAT-P-20 through HAT-P-23. Below we briefly describe the instruments used, the data reduction, and the analysis procedure.

Observations were made of all four planet host stars with the HIRES instrument (Vogt et al., 1994) on the Keck I telescope located on Mauna Kea, Hawaii. The width of the spectrometer slit was , resulting in a resolving power of , with a wavelength coverage of 3800–8000 Å. We typically used the B5 decker yielding a slit, and for the last few observations on each target we used the C2 decker that enables a better sky subtraction due to the longer slit . The slit height was oriented with altitude (vertical), except for rare cases, when the slit would have run through the faint companion to HAT-P-20 or HAT-P-22. A Keck/HIRES snapshot for each planet host star is shown in Figure 5. Spectra were obtained through an iodine gas absorption cell, which was used to superimpose a dense forest of lines on the stellar spectrum and establish an accurate wavelength fiducial (see Marcy & Butler, 1992). For each target an additional exposure was taken without the iodine cell, for use as a template in the reductions. 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.

In each of Figures 69 we show also the relative index, which is a measure of the chromospheric activity of the star derived from the flux in the cores of the Ca ii H and K lines. This index was computed following the prescription given by Vaughan, Preston & Wilson (1978), and as described in Hartman et al. (2009). Note that our relative index has not been calibrated to the scale of Vaughan, Preston & Wilson (1978). We do not detect any significant variation of the index correlated with orbital phase; such a correlation might have indicated that the RV variations could be due to stellar activity, casting doubt on the planetary nature of the candidate. There is no sign of emission in the cores of the Ca ii H and K lines in any of our spectra, from which we conclude that all of the targets have low chromospheric activity levels.

### 2.4. Photometric follow-up observations

In order to permit a more accurate modeling of the light curves, we conducted additional photometric observations with the KeplerCam CCD camera on the FLWO 1.2 m telescope for each star, and with the Faulkes North Telescope (FTN) of the Las Cumbres Observatory Global Network (LCOGT) at Hawaii for HAT-P-21 only. The observations for each target are summarized in Table 1.

The reduction of these images, including basic calibration, astrometry, and aperture photometry, was performed as described by Bakos et al. (2010). We found that the aperture photometry for HAT-P-20 was significantly affected by the close-by neighbor star 2MASS 07273995+2420118 with magnitude difference at 6.86″ separation (Figure 5). Thus, we performed image subtraction on the FLWO 1.2 m images with the same toolset used for the HATNet reductions, but applied a discrete kernel with half-size of 5 pixels and no spatial variations. Indeed, for this stellar configuration, the image subtraction results proved to be superior to the aperture photometry. For all of the follow-up light curves, we performed EPD and TFA to remove trends simultaneously with the light curve modeling (for more details, see Section 3, and Bakos et al. 2010). The final time series, together with our best-fit transit light curve model, are shown in the top portion of Figures 10 through 13; the individual measurements are reported in Tables 811, for HAT-P-20 through HAT-P-23, respectively.

## 3. Analysis

### 3.1. Properties of the parent stars

Fundamental parameters for each of the host stars, including the mass () and radius (), which are needed to infer the planetary properties, depend strongly on other stellar quantities that can be derived spectroscopically. For this we have relied on our template spectra obtained with the Keck/HIRES instrument, and the analysis package known as Spectroscopy Made Easy (SME; Valenti & Piskunov, 1996), along with the atomic line database of Valenti & Fischer (2005). For each star, SME yielded the following initial values and uncertainties (which we have conservatively increased to include our estimates of the systematic errors):

• HAT-P-20 – effective temperature  K, stellar surface gravity  (cgs), metallicity  dex, and projected rotational velocity .

• HAT-P-21 – effective temperature  K, stellar surface gravity  (cgs), metallicity  dex, and projected rotational velocity .

• HAT-P-22 – effective temperature  K, stellar surface gravity  (cgs), metallicity  dex, and projected rotational velocity .

• HAT-P-23 – effective temperature  K, stellar surface gravity  (cgs), metallicity  dex, and projected rotational velocity .

In principle the effective temperature and metallicity, along with the surface gravity taken as a luminosity indicator, could be used as constraints to infer the stellar mass and radius by comparison with stellar evolution models. However, the effect of  on the spectral line shapes is rather subtle, and as a result it is typically difficult to determine accurately, so that it is a rather poor luminosity indicator in practice. Unfortunately a trigonometric parallax is not available for any of the host stars, since they were not included among the targets of the Hipparcos mission (Perryman et al., 1997). For planetary transits, another possible constraint is provided by the  normalized semi-major axis, which is closely related to , the mean stellar density. The quantity  can be derived directly from the combination of the transit light curves (Sozzetti et al., 2007) and the RV data (required for eccentric cases, see Section 3.3). This, in turn, allows us to improve on the determination of the spectroscopic parameters by supplying an indirect constraint on the weakly determined spectroscopic value of  that removes degeneracies. We take this approach here, as described below. The validity of our assumption, namely that the adequate physical model describing our data is a planetary transit (as opposed to a blend), is shown later in Section 3.2.

For each system, our initial values of , , and  were used to determine auxiliary quantities needed in the global modeling of the follow-up photometry and radial velocities (specifically, the limb-darkening coefficients). This modeling, the details of which are described in Section 3.3, uses a Monte Carlo approach to deliver the numerical probability distribution of  and other fitted variables. For further details we refer the reader to Pál (2009b). When combining  (used as a proxy for luminosity) with assumed Gaussian distributions for  and  based on the SME determinations, a comparison with stellar evolution models allows the probability distributions of other stellar properties to be inferred, including . Here we use the stellar evolution calculations from the Yonsei-Yale group (YY; Yi et al., 2001) for all planets presented in this work. The comparison against the model isochrones was carried out for each of 10,000 Monte Carlo trial sets for HAT-P-21, HAT-P-22, and HAT-P-23, and for 20,000 Monte Carlo trial sets for HAT-P-20 (see Section 3.3). Parameter combinations corresponding to unphysical locations in the H-R diagram (26% of the trials for HAT-P-20, and less than 1% of the trials for the other objects) were ignored, and replaced with another randomly drawn parameter set. For each system we carried out a second SME iteration in which we adopted the value of  so determined and held it fixed in a new SME analysis (coupled with a new global modeling of the RV and light curves), adjusting only , , and . This gave:

• HAT-P-20: ,  K, , and .

• HAT-P-21: ,  K, , and .

• HAT-P-22: ,  K, , and .

• HAT-P-23: ,  K, , and .

In each case the conservative uncertainties for and have been increased by a factor of two over their formal values, as before. For each system, a further iteration did not change  significantly, so we adopted the values stated above as the final atmospheric properties of the stars. They are collected in Table 12.

With the adopted spectroscopic parameters the model isochrones yield the stellar mass and radius, and other properties. These are listed for each of the systems in Table 12. According to these models HAT-P-20 is a dwarf star with an estimated age of  Gyr, HAT-P-21 is a slightly evolved star with an estimated age of  Gyr, HAT-P-22 is a slightly evolved star with an estimated age of  Gyr, and HAT-P-23 is a slightly evolved star with an estimated age of  Gyr. The inferred location of each star in a diagram of  versus , analogous to the classical H–R diagram, is shown in Figure 14. In all cases the stellar properties and their 1 and 2 confidence ellipsoides are displayed against the backdrop of model isochrones for a range of ages, and the appropriate stellar metallicity. For comparison, the locations implied by the initial SME results are also shown (in each case with a triangle).

The stellar evolution modeling provides color indices (see Table 12) that may be compared against the measured values as a sanity check. For each star, the best available measurements are the near-infrared magnitudes from the 2MASS Catalogue (Skrutskie et al., 2006), which are given in Table 12. These are converted to the photometric system of the models (ESO system) using the transformations by Carpenter (2001). The resulting color indices are also shown in Table 12 for HAT-P-20 through HAT-P-23, respectively. Indeed, the colors from the stellar evolution models and from the observations agree for all of the host stars within 2-. The distance to each object may be computed from the absolute magnitude from the models and the 2MASS magnitudes, which has the advantage of being less affected by extinction than optical magnitudes. The results are given in Table 12, where in all cases the uncertainty excludes possible systematics in the model isochrones that are difficult to quantify.

### 3.2. Excluding blend scenarios

Our initial spectroscopic analyses discussed in Section 2.2 and Section 2.3 rule out the most obvious astrophysical false positive scenarios. However, more subtle phenomena such as blends (contamination by an unresolved eclipsing binary, whether in the background or associated with the target) can still mimic both the photometric and spectroscopic signatures we see. In the following section we investigate whether such scenarios may have caused the observed photometric and spectroscopic features.

#### Spectral line-bisector analysis

Following Torres et al. (2007), we explored the possibility that the measured radial velocities are not real, but are instead caused by distortions in the spectral line profiles due to contamination from a nearby unresolved eclipsing binary. A bisector span (BS) analysis for each system based on the Keck spectra was done as described in §5 of Bakos et al. (2007). In general, none of the Keck/HIRES spectra suffer significant sky contamination. Nevertheless, we calculated the Sky Contamination Factors (SCF) as described in Hartman et al. (2009), and corrected for the minor correlation between SCF and BS. The results are exhibited in Figure 15, where we show the SCF–BS and RV– (BS after SCF correction) plots for each planetary system. We also calculated the Spearman rank-order correlation coefficients (denoted as for the RV vs. BS quantities and the false alarm probabilities (see Table 13). There is no correlation for HAT-P-21, HAT-P-22 and HAT-P-23, and thus the interpretation of these systems as transiting planets is clear. There is an anti-correlation present for HAT-P-20, which is strengthened when the SCF correction is applied. A plausible explanation for this is that the neighboring star at 6″ separation (see Figure 5) is bleeding into the slit, even though we were careful during the observations to keep the slit centered on the main target, and adjusted the slit orientation to be perpendicular to the direction to the neighbor. We simulated this scenario, and calculated the expected BS as a function of RV due to the neighbor, assuming that the two stars have the same systemic velocity and the seeing is 1″. Indeed, we get a slight anti-correlation from this simulation, and the range of magnitude in the BS variation is consistent with the observations. It is also possible that some of the anti-correlation is due to the fact that the slit was not in vertical angle for many of the observations. The non-vertical slit mode may result in wavelength-dependent slit losses due to atmospheric dispersion, and this could bring in a correlation with the sky background, and change the shape of the spectral lines.

### 3.3. Global modeling of the data

This section describes the procedure we followed for each system to model the HATNet photometry, the follow-up photometry, and the radial velocities simultaneously. 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, with limb darkening being prescribed by a quadratic law. The limb darkening coefficients for the Sloan  band, Sloan  band, and Sloan  band were interpolated from the tables by Claret (2004) for the spectroscopic parameters of each star as determined from the SME analysis (Section 3.1). The transit shape was parametrized by the normalized planetary radius , the square of the impact parameter , and the reciprocal of the half duration of the transit . We chose these parameters because of their simple geometric meanings and the fact that these show negligible correlations (see Bakos et al., 2010). The relation between and the quantity , used in Section 3.1, is given by

 a/R⋆=P/2π(ζ/R⋆)√1−b2√1−e2/(1+esinω) (1)

(see, e.g., Tingley & Sackett, 2005). Note the subtle dependency of  on the and Lagrangian orbital parameters that are typically derived from the RV data ( is the longitude of periastron). This dependency is often ignored in the literature, and  is quoted as a “pure” light curve parameter. Of course, if high quality secondary eclipse observations are available that determine both the location and duration of the occultation, then and can be determined without RV data. Our model for the HATNet data was the simplified “P1P3” version of the Mandel & Agol (2002) analytic functions (an expansion in terms of Legendre polynomials), for the reasons described in Bakos et al. (2010). Following the formalism presented by Pál (2009), the RVs were fitted with an eccentric Keplerian model parametrized by the semi-amplitude and Lagrangian elements and . Note that we allowed for an eccentric orbit for all planets, even if the results were consistent with a circular orbit. There are several reasons for this: i) many of the close-in hot Jupiters show eccentric orbits, thus the assumption of fixing has no physical justification (while this has been customary in early discoveries relying on very few data-points) ii) the error-bars on various other derived quantities (including ) are more realistic with the inclusion of eccentricity, and iii) non-zero eccentricities can be very important in proper interpretation of these systems.

We assumed that there is a strict periodicity in the individual transit times. For each system we assigned the transit number to a complete follow-up light curve. For HAT-P-20b this was the light curve gathered on 2009 Oct 21, for HAT-P-21b: 2010 Feb 19, HAT-P-22b: 2009 Feb 28, and HAT-P-23b: 2008 Sep 13. The adjustable parameters in the fit that determine the ephemeris were chosen to be the time of the first transit center observed with HATNet (, , , and for HAT-P-20b through HAT-P-23b respectively) and that of the last transit center observed with the FLWO 1.2 m telescope (, , , and for HAT-P-20b through HAT-P-23b respectively). We used these as opposed to period and reference epoch in order to minimize correlations between parameters (see Pál et al., 2008). Times of mid-transit for intermediate events were interpolated using these two epochs and the corresponding transit number of each event, . The eight main parameters describing the physical model for each system were thus the first and last transit center times, , , , , , and . For HAT-P-20b, HAT-P-22b, and HAT-P-23b three additional parameters were included (for each system) that have to do with the instrumental configuration (blend factor, out-of-transit magnitudes, gamma velocities; see later). For HAT-P-21b seven additional parameters were included, because it was observed in 3 different HATNet fields. These are the HATNet blend factor (one for each HATNet field), which accounts for possible dilution of the transit in the HATNet light curve from background stars due to the broad PSF (20″ FWHM), the HATNet out-of-transit magnitude (one for each HATNet field), and the relative zero-point of the Keck RVs.

We extended our physical model with an instrumental model that describes brightness variations caused by systematic errors in the measurements. This was done in a similar fashion to the analysis presented by Bakos et al. (2010). The HATNet photometry has already been EPD- and TFA-corrected before the global modeling, so we only considered corrections for systematics in the follow-up light curves. We chose the “ELTG” method, i.e., EPD was performed in “local” mode with EPD coefficients defined for each night, and TFA was performed in “global” mode using the same set of stars and TFA coefficients for all nights. The five EPD parameters were the hour angle (representing a monotonic trend that changes linearly over time), the square of the hour angle (reflecting elevation), and the stellar profile parameters (equivalent to FWHM, elongation, and position angle of the image). The functional forms of the above parameters contained six coefficients, including the auxiliary out-of-transit magnitude of the individual events. For each system the EPD parameters were independent for all nights, implying 12, 18, 12, and 36 additional coefficients in the global fit for HAT-P-20b through HAT-P-23b respectively. For the global TFA analysis we chose 20, 3, 10, and 20 template stars for HAT-P-20b through HAT-P-23b that had good quality measurements for all nights and on all frames, implying an additional 20, 3, 10, and 20 parameters in the fit for each system. In all cases the total number of fitted parameters (43, 36, 33 and 67 for HAT-P-20b through HAT-P-23b) was much smaller than the number of data points (755, 1172, 892 and 953, counting only RV measurements and follow-up photometry measurements).

The joint fit was performed as described in Bakos et al. (2010). We minimized  in the space of parameters by using a hybrid algorithm, combining the downhill simplex method (AMOEBA; see Press et al., 1992) with a classical linear least squares algorithm. Uncertainties for the parameters were derived by applying the Markov Chain Monte-Carlo method (MCMC, see Ford, 2006). This provided the full a posteriori probability distributions of all adjusted variables. The a priori distributions of the parameters for these chains were chosen to be Gaussian, with eigenvalues and eigenvectors derived from the Fisher covariance matrix for the best-fit solution. The Fisher covariance matrix was calculated analytically using the partial derivatives given by Pál (2009).

Following this procedure we obtained the a posteriori distributions for all fitted variables, and other quantities of interest such as . As described in Section 3.1,  was used together with stellar evolution models to infer a value for  that is significantly more accurate than the spectroscopic value. The improved estimate was in turn applied to a second iteration of the SME analysis, as explained previously, in order to obtain better estimates of  and . The global modeling was then repeated with updated limb-darkening coefficients based on those new spectroscopic determinations. The resulting geometric parameters pertaining to the light curves and velocity curves for each system are listed in Table 14.

Included in each table is the RV “jitter”, which is a noise term that we added in quadrature to the internal errors for the RVs in order to achieve from the RV data for the global fit. The jitter is a combination of assumed astrophysical noise intrinsic to the star, plus instrumental noise rising from uncorrected intstrumental effects (such as a template spectrum taken under sub-optimal conditions).

The planetary parameters and their uncertainties were derived by combining the a posteriori distributions for the stellar, light curve, and RV parameters. In this way we find masses and radii for each planet. These and other planetary parameters are listed at the bottom of Table 14, and further discussed in § 4.

## 4. Discussion

### 4.1. HAT-P-20b

HAT-P-20b is a very massive (