{'fraction_non_alphanumeric': 0.048892459543211295, 'fraction_numerical': 0.0381498909675198, 'mean_word_length': 4.286494357480888, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
HAT-P-20b--HAT-P-23b: Four Massive Transiting Extrasolar Planets

We report the discovery of four relatively massive (2-7MJ) transiting extrasolar planets. HAT-P-20b orbits a V=11.339 K3 dwarf star with a period P=2.875317+/-0.000004d. The host star has a mass of 0.760+/-0.03 Msun, radius of 0.690+/-0.02 Rsun, Teff=4595+/-80 K, and metallicity [Fe/H]=+0.35+/-0.08. HAT-P-20b has a mass of 7.246+/-0.187 MJ, and radius of 0.867+/-0.033 RJ yielding a mean density of 13.78+/-1.50 gcm^-3 , which is the second highest value among all known exoplanets. HAT-P-21b orbits a V=11.685 G3 dwarf on an eccentric (e=0.2280+/-0.016) orbit, with a period of P=4.1244810+/-000007d. The host star has a mass of 0.95+/-0.04Msun, radius of 1.10+/-0.08Rsun, Teff=5588+/-80K, and [Fe/H]=+0.01+/-0.08. HAT-P-21b has a mass of 4.063+/-0.161MJ, and radius of 1.024+/-0.092RJ. HAT-P-22b orbits the V=9.732 G5 dwarf HD233731, with P=3.2122200+/-0.000009d. The host star has a mass of 0.92+/-0.03Msun, radius of 1.04+/-0.04Rsun, Teff=5302+/-80K, and metallicity of +0.24+/-0.08. The planet has a mass of 2.147+/-0.061 MJ, and compact radius of 1.080+/-0.058RJ. The host star also harbors an M-dwarf companion at a wide separation. Finally, HAT-P-23b orbits a V=12.432 G0 dwarf star, with a period P=1.212884+/-0.000002d. The host star has a mass of 1.13+/-0.04sun, radius of 1.20+/-0.07Rsun, Teff=5905+/-80K, and [Fe/H]=+0.15+/-0.04. The planetary companion has a mass of 2.090+/-0.111MJ, and radius of 1.368+/-0.090RJ (abridged).

The majority of the ∼ 90 known transiting extrasolar planets (TEPs) have been found to lie in the 0.5 M J to 2.0 M J mass range. The apparent drop in their mass distribution at ∼2 M J has been noted by, e.g., Southworth et al. (2009), and by . In the currently known sample, 75% of the TEPs have planetary mass M p < 2.0 M J , and there appears to be a minor peak in their occurrence rate at M p ≈ 2 M J , 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 M p − R p relation happens around M p ≈ 2 − 3 M J , and falls off as R p ∝ M −1/8 p (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 ef-ficiency 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 v sin i for host stars with planets M p < 2M J is 3.9 m s −1 (with 3.8 m s −1 standard deviation around the mean), whereas the same values for the massive planet hosts stars, are 10.9 m s −1 (with 12.7 m s −1 standard deviation around the mean) 13 . The five fastest rotators all harbor planets more massive than 2 M J . 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 eccentric 14 (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 m s −1 level (see HAT-P-2, where valuable data was contributed to the RV fit by modest precision instruments yielding ∼ 1 km s −1 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 2 M J 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 M J (Corot-7b; Queloz et al. 2009) to 26.4 M J (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;, and the Kozai-mechanism (Fabrycky & Tremaine 2007). They are border-line objects between planets and browndwarfs, 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 M J . 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. 13 This includes those 4 planets announced in this paper. 14 See e.g. exoplanets.org for statistics. 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 M ⋆ 1.2 M ⊙ 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, andGSC 1632-01396. This extends the currently known sample of bright (V < 13.5) and massive (M p > 2 M J ) 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 Super-WASP (Pollacco et al. 2006), TrES (Alonso et al. 2004) andXO projects (McCullough et al. 2005), and spaceborne 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 (8  I  14). 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. 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:

Photometric detection
• HAT-P-20 -GSC 1910-00239 (also known as 2MASS 07273995+2420118; α = 07 h 27 m 39.96 s , δ = +24 • 20 ′ 11.9 ′′ ; J2000; V=11.339). A signal was detected for this star with an apparent depth of ∼10.8 mmag, and a period of P =2.8753 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 ′′ pixel −1 ) 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 q =0.0268 ± 0.0003, corresponding to a total duration of P q =1.848 ± 0.019 hr. HAT-P-20 has a red companion (2MASS 07273963+2420171, J − K = 0.92) at 6.86 ′′ separation that is fainter than HAT-P-20 by ∆R = 1.36 mag.
• HAT-P-23 -GSC 1632-01396 (also known as 2MASS 20242972+1645437; α = 20 h 24 m 29.88 s , δ = +16 • 45 ′ 43.7 ′′ ; J2000; V=12.432). A signal was detected for this star with an apparent depth of ∼ 11.5 mmag, and a period of P =1.2129 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 q =0.0748 ± 0.0006, corresponding to a total duration of P q =2.178 ± 0.017 hr.  Table 1 for details), and folded with the period P = 2.8753172 days resulting from the global fit described in Section 3. The solid line shows the "P1P3" transit model fit to the light curve (Section 3.3). The bold points in the lower panel show the light curve binned in phase with a bin-size of 0.002. -Unbinned light curve of HAT-P-21 including all 28,000 instrumental I band and R band 5.5 minute cadence measurements obtained with the HAT-5, HAT-6, HAT-8 and HAT-9 telescopes of HATNet (see Table 1 for details), and folded with the period P = 4.1244808 days resulting from the global fit described in Section 3. The solid line shows the "P1P3" transit model fit to the light curve (Section 3.3). The bold points in the lower panel show the light curve binned in phase with a bin-size of 0.002.

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 km s −1 ) 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 (λ/∆λ ≈ 35,000) over a single order centered on the Mg I b triplet (∼5187Å), with typically low signal-to-noise (S/N) ra-  Table 1 for details), and folded with the period P = 1.2128841 days resulting from the global fit described in Section 3. The solid line shows the "P1P3" transit model fit to the light curve (Section 3.3). The bold points in the lower panel show the light curve binned in phase with a bin-size of 0.002. tios that are nevertheless sufficient to derive radial velocities (RVs) with moderate precisions of 0.5-1.0 km s −1 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 mediumresolution fiber which produces spectra at a resolution of λ/∆λ ≈ 46,000 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 (∼ 0.5 km s −1 ) 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.

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 periodfolded 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 0. ′′ 86, resulting in a resolving power of λ/∆λ ≈ 55,000, with a wavelength cov-erage of ∼3800-8000Å. We typically used the B5 decker yielding a 3.5 ′′ (H) × 0.861 ′′ (W ) 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 14.0 ′′ (H) × 0.861 ′′ (W ). 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 I 2 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.  Table 14). Zero phase corresponds to the time of mid-transit. The center-of-mass velocity has been subtracted. The rms around the best orbital fit is 16.2 m s −1 . Second panel: Velocity O − C residuals from the best fit. The error bars for both the top and second panel include a component from the jitter (15.7 m s −1 ) added in quadrature to the formal errors (see Section 3.3). Third panel: Bisector spans (BS), with the mean value subtracted. The measurement from the template spectrum is included (see Section 3.2). Bottom panel: Relative chromospheric activity index S measured from the Keck spectra.
In each of Figures 6-9 we show also the relative S 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 pre-  scription given by Vaughan, Preston & Wilson (1978), and as described in Hartman et al. (2009). Note that our relative S 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.

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 ∆i = 1.1 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  TABLE 4 Relative radial velocities, bisector spans, and activity index measurements of HAT-P-20.

ANALYSIS
3.1. Properties of the parent stars Fundamental parameters for each of the host stars, including the mass (M ⋆ ) and radius (R ⋆ ), 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 T eff⋆ = 4626 ± 104 K, stellar surface gravity log g ⋆ = 4.80 ± 0.06 (cgs), metallicity [Fe/H] = +0.31 ± 0.08 dex, and projected rotational velocity v sin i = 3.1 ± 0.5 km s −1 .
In principle the effective temperature and metallicity, along with the surface gravity taken as a luminosity in-  Figure 6. The parameters used in the best-fit model are given in Table 14, the RV jitter was 9.9 m s −1 , and the fit rms was  dicator, could be used as constraints to infer the stellar mass and radius by comparison with stellar evolution models. However, the effect of log g ⋆ 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 a/R ⋆ normalized semi-major axis, which is closely related to ρ ⋆ , the mean stellar density. The quantity a/R ⋆ can be derived directly from the combination of    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 log g ⋆ 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 T eff⋆ , log g ⋆ , and [Fe/H] 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 a/R ⋆ and other fitted variables. For further details we refer the reader to Pál (2009b). When combining a/R ⋆ (used as a proxy for luminosity) with assumed Gaussian distributions for T eff⋆ and [Fe/H] based on the SME deter-minations, a comparison with stellar evolution models allows the probability distributions of other stellar properties to be inferred, including log g ⋆ . 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 log g ⋆ 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 T eff⋆ , [Fe/H], and v sin i. This gave: • HAT-P-20: log g ⋆ = 4.64 ± 0.06, T eff⋆ = 4595 ± 80 K, [Fe/H] = +0.35 ± 0.08, and v sin i = 2.1 ± 0.5 km s −1 .
In each case the conservative uncertainties for T eff⋆ and [Fe/H] have been increased by a factor of two over their formal values, as before. For each system, a further iteration did not change log g ⋆ 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 6.7 +5.7 −3.8 Gyr, HAT-P-21 is a slightly evolved star with an estimated age of 10.2 ± 2.5 Gyr, HAT-P-22 is a slightly evolved star with an estimated age of 12.4 ± 2.6 Gyr, and HAT-P-23 is a slightly evolved star with an estimated age of 4.0 ± 1.0 Gyr. The inferred location of each star in a diagram of a/R ⋆ versus T eff⋆ , 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 K magnitude from the models and the 2MASS K s 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.

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 SCF (BS after SCF correction) plots for each planetary system. We also calculated the Spearman rank-order correlation coefficients (denoted as R s 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 g band, Sloan i band, and Sloan z 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 p ≡ R p /R ⋆ , the square of the impact parameter b 2 , 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 Bakos et al. 2010). The relation between ζ/R ⋆ and the quantity a/R ⋆ , used in Section 3.1, is given by (1) (see, e.g., Tingley & Sackett 2005). Note the subtle dependency of a/R ⋆ on the k ≡ e cos ω and h ≡ e sin ω 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 a/R ⋆ 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 k and h 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 K and Lagrangian elements k and h. 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 e = 0 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 a/R ⋆ ) 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 YY+a/R⋆+SME a SME = "Spectroscopy Made Easy" package for the analysis of high-resolution spectra (Valenti & Piskunov 1996). These parameters rely primarily on SME, but have a small dependence also on the iterative analysis incorporating the isochrone search and global modeling of the data, as described in the text. b YY+a/R⋆+SME = Based on the YY isochrones (Yi et al. 2001), a/R⋆ as a luminosity indicator, and the SME results.  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, N tr . The eight main parameters describing the physical model for each system were thus the first and last transit center times, R p /R ⋆ , b 2 , ζ/R ⋆ , K, k ≡ e cos ω, and h ≡ e sin ω.
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 B inst (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 M 0,HATNet (one for each HAT-Net field), and the relative zero-point γ rel 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 EPDand TFA-corrected before the global modeling, so we only considered corrections for systematics in the followup 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 (equiv- alent to FWHM, elongation, and position angle of the image). The functional forms of the above parameters contained six coefficients, including the auxiliary out-oftransit magnitude of the individual events. For each sys-tem 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 χ 2 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 a/R ⋆ . As described in Section 3.1, a/R ⋆ was used together with stellar evolution models to infer a value for log g ⋆ 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 T eff⋆ and [Fe/H]. 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 χ 2 /dof = 1 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 suboptimal 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.

HAT-P-20b
HAT-P-20b is a very massive (M p = 7.246 ± 0.187 M J = 2302.9 ± 59.5 M ⊕ ) and very compact (R p = 0.867 ± 0.033 R J ) hot Jupiter orbiting a K3 (Skiff 2009) star. HAT-P-20b is the sixth most massive, and second most dense transiting planet with ρ p =13.78 ± 1.50 g cm −3 (see Figure 16). The only planet (or brown dwarf) denser than HAT-P-20b is CoRoT-3b (Deleuil et al. 2008) with ρ p ≈ 27g cm −3 . Modeling HAT-P-20b may be a challenge, as the oldest (4 Gyr, i.e. yielding the most compact planets) Fortney et al. (2007) models with M p = 2154 M ⊕ total mass and 100 M ⊕ core-mass predict a much bigger radius (1.04 R J ). The observed radius of 0.87 M J would require a very high metal content. We note that the host star is one of the most metal rich stars that have a transiting planet ([Fe/H] = +0.35 ± 0.08). Curiously, HAT-P-20b orbits a fairly late type star (K3), as compared to most of the massive hot Jupiters that orbit ∼F5 dwarfs. It is also different from the rest of the population in that the orbit is consistent with circular at the 3σ level. The irradiation HAT-P-20b receives is one of the smallest, clearly making it a pL class exoplanet (Fortney et al. 2008): F = (2.00 ± 0.19) · 10 8 erg s −1 cm −2 , comparable to the mean flux per orbit for another "heavy" planet HD 17156b on a 21 d period orbit. The only other massive planet that receives less average flux (integrated over an orbit) is HD 80606b. HAT-P-20 is an extreme outlier in the M p -M ⋆ plane; it is a relatively small mass star harboring a very massive planet. Another outlier (albeit to a much lesser extent) with similar planetary radius and stellar mass is WASP-10b (Johnson et al. 2009;Christian et al. 2009), but this planet has less than half of the mass of HAT-P-20b (3.09 M J ). We also calculated the maximum mass of a stable moon for both the prograde and retrograde cases, and derived 0.128 M ⊕ and 8.31 M ⊕ , respectively, i.e. HAT-P-20b can harbor a fairly massive moon. An 8.31 M ⊕ retrograde moon would cause ∼ 10 s variations in the transit times, which is marginally detectable from the ground.
HAT-P-20 has a close-by faint and red companion star at ∼ 6.86 ′′ separation. Based on the Palomar sky survey archival plates, we confirm that they form a close common-proper motion pair, thus it is very likely that the two stars are physically associated. The binary has appeared in the Washington Double Star compilation (WDS) as POU2795, and was discovered by Pourteau (1933). Furthermore, based on the summary of observations in the WDS, there is already a hint of orbital motion of the companion to HAT-P-20 over the last century. The position angle of the companion changed from PA=323 • to PA=320 • over the course of 89 years (between 1909 and 1998), and it seems to be retrograde on the sky (clockwise). Thus, HAT-P-20 is yet another example of a massive planet in a binary system (Udry et al. 2002). The binary companion makes this system ideal for high precision ground or space-based studies, as it provides a natural comparison source, even though it has a later spectral type.

HAT-P-21b
With a mass of M p = 4.063 ± 0.161 M J , HAT-P-21b is the 11th most massive transiting planet. HAT-P-21b has a radius of R p = 1.024 ± 0.092 R J , mean density ρ p = 4.68 +1.59 −0.99 g cm −3 , and orbits on a moderately eccentric orbit with e = 0.228 ± 0.016, ω = 309 ± 3 • . The transits occur near apastron. As noted by Buchhave et al. (2010), 4 M J mass planets are very rare in the sample of currently known transiting exoplanets, and the only sib -TABLE 14 Orbital and planetary parameters for HAT-P-20b-HAT-P-23b is also far from grazing, where the transit is short, and other system parameters have lower accuracy. The effective temperature of the star (T eff⋆ = 5905±80 K) is close to the critical temperature of 6250 K noted recently by Winn et al. (2010), which may be a border-line between systems where the stellar spin axes and planetary orbital normals are preferentially aligned (T eff⋆ < 6250 K) and those that are misaligned (T eff⋆ > 6250 K). An alternative hypothesis has been brought up by Schlaufman (2010), where misaligned stellar spin axes and orbital normals are related to the mass of the host star. The mass of HAT-P-23 (1.13 ± 0.04) is sufficiently close to the suggested dividing line of M ⋆ = 1.2 M ⊙ , thus it will provide an excellent additional test for these ideas.

Summary
We presented the discovery of four new massive transiting planets, and provided accurate characterization of the host star and planetary parameters. These 4 new systems are very diverse, and significantly expand the sample of ∼ 13 other massive (M p 2 M J ) planets. Two of the new discoveries orbit stars that have fainter, most probably physically associated companions. The new discoveries do not tend to enhance the masseccentricity correlation, since only one (HAT-P-21b) is significantly eccentric. Also, the tentative mass-v sin i correlation noted in the Introduction is weakened by the new discoveries. The heavier mass planets (HAT-P-20b and HAT-P-21b) seem to be inconsistent with current theoretical models in that they are too dense, and would require a huge core (or metal content) to have such small radii. One planet (HAT-P-23b) is also inconsistent with the models (unless we assume that the planet is very young), but in the opposite sense; it has an inflated radius. It has been noted by Winn et al. (2010) and Schlaufman (2010) that systems exhibiting stellar spin axis-planetary orbital normal misalignment are preferentially eccentric and heavy mass planets (in addition to the key parameter being the effective temperature or mass of the host star, respectively). The four new planets presented in this work will provide additional important tests for checking these conjectures. The host stars are all bright (9.7 < V < 12.4), and thus enable in-depth future characterization of these systems.