FRIENDS OF HOT JUPITERS. I. A RADIAL VELOCITY SEARCH FOR MASSIVE, LONG-PERIOD COMPANIONS TO CLOSE-IN GAS GIANT PLANETS

Our sample includes transiting planets with orbital periods between 0.7–11 days and masses between 0.06–11 M_Jup (i.e., planet with masses comparable to or larger than that of Neptune). We focus our search on a sample of twenty seven systems where there is already evidence for multi-body dynamics, including planets with eccentric orbits or orbits that are significantly tilted with respect to the star's spin axis (see Figure 1). We required that the planets in this sample have projected obliquities or eccentricities that differed from zero by more than 3σ; for convenience we refer to this as the "misaligned" sample, although we note that it also contains planets with eccentric orbits and obliquities consistent with zero. We also include a control sample of twenty four planets that appear to have well-aligned and circular orbits (i.e., within 3σ of zero), where canonical disk migration models for isolated stars could plausibly explain the presence of the observed short-period planet. Because the stellar multiplicity rate increases for more massive stars, we select our control sample to have approximately the same distribution of stellar masses as our misaligned sample in order to avoid biasing our estimates of the companion frequencies (see Figure 2 for the relative distribution and Table 1 for a list of masses for individual systems). A subset of the systems in our target list are known to exhibit radial velocity accelerations; in these cases, our data allow us to confirm and refine the properties of the long-period companion responsible for the trend.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We observed our target stars using the High Resolution Echelle Spectrometer (HIRES) (Vogt et al. 1994) on the 10 m Keck I telescope over a period of two years beginning in 2011; many of our targets also had existing HIRES observations taken prior to 2011 by other programs. We used the standard HIRES setup and reduction pipeline employed by the California Planet Search (CPS) consortium (Wright et al. 2004; Howard et al. 2009; Johnson et al. 2010). Observations were typically obtained with a slit width of 086 with integration times optimized to obtain typical signal to noise ratios of 70 per pixel. An iodine cell mounted in front of the spectrometer entrance slit provided a wavelength scale and instrumental profile for the observations (Marcy & Butler 1992; Valenti et al. 1995). We obtained a total of approximately 270 new radial velocity measurements for our target sample, with a minimum of four observations per target separated by at least six months. We then combine our data with published radial velocities obtained using other telescopes to provide the strongest possible constraints on the presence of any long-term radial velocity accelerations. We provide a summary of the radial velocity data utilized in this study in Table 2, as well as individual HIRES radial velocity measurements for each system in Table 3.

In this paper we focus on images obtained for systems with detected radial velocity accelerations; we will present a complete analysis of our AO data set including companion detections in Paper II. We obtained K band adaptive optics imaging (Wizinowich 2000) for each of our target stars using the NIRC2 instrument (Instrument PI: Keith Matthews) on Keck II in the narrow camera (10 mas pixel⁻¹) setting. We used the full 1024 × 1024 pixel field of view for most of our target stars, with the exception of several of our brightest targets where we switched to a 512 × 512 pixel subarray in order to allow for shorter integration times and avoid saturation. We utilized a standard three-point dither pattern (e.g., Bechter et al. 2013) that maximizes our spatial coverage and allows for the removal of sky and instrumental backgrounds while avoiding the lower-left quadrant on the array, which has a higher read noise level. We obtain our images in position angle mode, where the orientation of the image on the detector is kept constant as the telescope tracks, rather than using the angular differential imaging technique where the image is allowed to rotate on the detector and performing point spread function (PSF) subtraction. This maximized the efficiency of our observations while still providing deep sensitivity to low-mass stellar companions (Crepp et al. 2012).

We flat-field our images and remove hot pixels by searching for 4σ outliers at a fixed pixel position, treating each nod position separately. We calculate a median sky background using the off-nod positions and subtracting this median image from each of our science images at that nod position. Each image is then interpolated by a factor of ten and the images are stacked using the point spread function of our target star in order to align the positions. We create our final science images by taking the median flux at each pixel position in our stacked images. A summary of the observations utilized in this analysis is provided in Table 4.

In order to detect and quantify the significance of long-term accelerations in the radial velocity data we performed a uniform analysis of all 51 systems with a Differential-evolution Markov Chain Monte Carlo (DE-MCMC; Ter Braak 2006) technique similar to that of Fulton et al. (2013). The DE-MCMC algorithm speeds convergence by downgrading the importance of pre-determining optimal step sizes for each parameter. DE-MCMC runs many chains in parallel (twice the number of free parameters) and uses the difference in parameter values from two random chains in order to establish the magnitude and direction of each step. This ensures that step sizes are optimized on-the-fly to achieve ideal acceptance rates (∼18% for well-constrained fits) and high convergence rates. Step sizes for correlated parameters are automatically reduced in the direction orthogonal to the correlation which leads to fewer models being calculated in regions of parameter space that are highly disfavored by the data.

Our radial velocity model for each system was described by a minimum of eight free parameters: period (P), time of mid-transit at a particular reference epoch (T_mid), eccentricity (e), argument of periastron of the star's orbit (ω_⋆), velocity semi-amplitude (K), a relative radial velocity (RV) zero point (γ), slope ($\dot{\gamma }$), and RV "jitter." When required, we expanded on this baseline model by carrying out two-planet fits for systems where the outer companion's orbit exhibited significant curvature (HAT-P-17, WASP-8, WASP-34), and a three-planet fit for the HAT-P-13 system. For some systems, data from multiple spectrographs were included and in these cases the relative RV zero-points (γ) were fit separately for each instrument. GJ 436b has HIRES radial velocities obtained prior to the CCD upgrade, and we treat data before and after this upgrade as separate data sets with a different baseline normalization.

RV "jitter" may be dominated by instrumental effects as opposed to astrophysical noise and thus should not be expected to converge to the same value for different instruments (Isaacson & Fischer 2010). In order to prevent our fitted "jitter" parameter from being driven to abnormally large values by particularly noisy data sets we first run a set of chains with a uniform jitter value for all data sets to obtain a best-fitting model. We then run the chains again, this time scaling the jitter value at each step by σ_x/σ_CPS where σ_x is the rms of the residuals to the best-fit model from the initial run for data set x and σ_CPS is the rms of the residuals of the best-fit model for the post-upgrade HIRES data. This ensures that the measurement errors are roughly equal to the rms of the residuals to the final model for each individual data set. We also reject RV measurements from the CPS HIRES data with reported measurement errors that are greater than 10 times the median absolute deviation of all of the measurement errors for that particular star. These measurements are typically derived from very low signal-to-noise spectra where the standard HIRES extraction routine does not produce optimal results, and contribute minimally to our fits. This step generally results in the rejection of less than three outliers from each RV set.

RV measurements taken during transits of the known planet were excluded from the fit. For planets with high-cadence Rossiter measurements spanning several hours around the transit we take the error-weighted mean of the out-of-transit points and include this as a single measurement in our fits. Because we add an additional jitter term, this effectively down-weights the contribution of these high-density data sets to our fit. This conservative approach ensures that our best-fit solutions are not biased by the presence of short-term stellar variability that can cause trends in the RV measurements over several hour time scales (e.g., Albrecht et al. 2012a).

We computed 2× N DE-MCMC chains (where N is the number of free parameters in the RV model), continuously checking for convergence following the prescription of Eastman et al. (2013). We considered the chains well-mixed and halted the DE-MCMC run when the number of independent draws (T_z, as defined in Ford 2006) was greater than 1000 and the Gelman–Rubin statistic (Gelman et al. 2003; Ford 2006; Holman et al. 2006) was within 1% of unity for all parameters. In order to speed convergence (however see Section 3.2), ensure that all parameter space was adequately explored, and minimize biases in parameters that physically must be finite and positive, we step in the widely used modifications and/or combinations of orbital parameters: log(P), $\sqrt{e}\cos {\omega _{\star }}$, $\sqrt{e}\sin {\omega _{\star }}$, and log(K).

We assigned Gaussian priors to P, T_mid and secondary eclipse times where available as listed in Tables 5 and 6, and we assigned uniform priors to all other parameters. The reference epoch (abscissa) for $\dot{\gamma }$ was chosen as the mid-time of the RV time-series in order to minimize the covariance between γ and $\dot{\gamma }$. In all cases we assumed that transit timing variations caused by other known or unknown companions were negligible. The median parameter values and associated 68% confidence intervals from the DE-MCMC analysis for all systems are presented in Tables 7 and 8.

Table 7. Results from Radial Velocity Fits

Planet	e	ω⋆	ecos ω⋆	esin ω⋆	K	$\dot{\gamma }$a	Jitter	Mp	Spin-Orbit λ	Sampleb	Ref.
		(degrees)			(m s−1)	(m s−1 day−1)	(m s−1)	(MJup)	(deg)
GJ436b	0.1495$^{+0.016}_{-0.0097}$	336$^{+12}_{-11}$c	0.13654$^{+0.0004}_{-0.00047}$	−0.061$^{+0.032}_{-0.033}$	17.01 ± 0.54c	−0.00137 ± 0.00061	3.78$^{+0.32}_{-0.29}$	0.0682 ± 0.0025		Misaligned	1
HAT-P-2b	0.5079$^{+0.00093}_{-0.00079}$	186.2 ± 1.1	−0.50493$^{+0.00043}_{-0.00042}$	−0.0546 ± 0.0099	929 ± 11	−0.0938$^{+0.0067}_{-0.0069}$	31.7$^{+4.3}_{-3.5}$	9.99 ± 0.23	9 ± 10	Misaligned	3, 4
HAT-P-4b	0.004$^{+0.017}_{-0.0031}$	157$^{+110}_{-67}$	−0.00055$^{+0.00084}_{-0.0011}$	0.0001$^{+0.011}_{-0.0075}$	77 ± 3d	0.0219 ± 0.0035	9.9$^{+2.1}_{-1.6}$	0.639$^{+0.043}_{-0.042}$	−4.9 ± 11.9	Control	1, 5
HAT-P-6b	0.023$^{+0.022}_{-0.02}$	94.1$^{+37.0}_{-2.3}$	−0.00159$^{+0.00066}_{-0.00067}$	0.023 ± 0.022	120.8 ± 2.4	0.0041 ± 0.0019	6.6$^{+1.9}_{-1.6}$	1.107 ± 0.041	165 ± 6	Misaligned	1, 3
HAT-P-7b	0.0055$^{+0.007}_{-0.0033}$	204$^{+53}_{-89}$	−0.003 ± 0.002	−0.0008$^{+0.0053}_{-0.0086}$	214.3$^{+2.6}_{-2.5}$	0.0646$^{+0.004}_{-0.0038}$	12.8$^{+1.7}_{-1.4}$	1.697$^{+0.027}_{-0.026}$	155 ± 37	Misaligned	3, 6
HAT-P-8b	0.0029$^{+0.016}_{-0.0024}$	116$^{+150}_{-35}$	−7e-05$^{+0.0007}_{-0.00069}$	0.0003$^{+0.014}_{-0.0035}$	162.0$^{+4.4}_{-3.9}$d	−0.0003$^{+0.0034}_{-0.004}$	9.1$^{+3.5}_{-2.3}$	1.304$^{+0.065}_{-0.064}$	$-17^{+9.2}_{-11.5}$	Control	7, 8
HAT-P-10b	0.028$^{+0.029}_{-0.02}$	146$^{+95}_{-52}$	−0.011$^{+0.015}_{-0.026}$	0.008$^{+0.031}_{-0.017}$	75.7$^{+2.7}_{-2.6}$	−0.014$^{+0.0032}_{-0.0031}$	6.1$^{+2.1}_{-1.4}$	0.509 ± 0.023		Control	9
HAT-P-11b	0.232$^{+0.054}_{-0.053}$	7$^{+24}_{-25}$	0.213$^{+0.049}_{-0.053}$	0.028$^{+0.097}_{-0.092}$	10.2$^{+1.1}_{-1.2}$	0.0094 ± 0.0016e	5.95$^{+0.58}_{-0.52}$	0.0756 ± 0.0087	$103^{+26}_{-10}$	Misaligned	10, 11
HAT-P-12b	0.026$^{+0.026}_{-0.018}$	97$^{+220}_{-64}$	0.012$^{+0.021}_{-0.014}$	0.004$^{+0.031}_{-0.019}$	35.4 ± 1.6	−0.0004$^{+0.0018}_{-0.0019}$	4.23$^{+1.1}_{-0.92}$	0.2089$^{+0.01}_{-0.0097}$		Control	12
HAT-P-13b	0.0133$^{+0.0047}_{-0.0044}$	197$^{+32}_{-37}$	−0.0107$^{+0.0039}_{-0.0041}$	−0.0032$^{+0.0073}_{-0.0077}$	105.87 ± 0.78	0.0528$^{+0.0013}_{-0.0014}$	4.53$^{+0.54}_{-0.46}$	0.899$^{+0.03}_{-0.029}$	1.9 ± 8.6	Misaligned	13
HAT-P-14b	0.115$^{+0.015}_{-0.016}$	98.8$^{+5.4}_{-5.2}$	−0.02 ± 0.01	0.113$^{+0.015}_{-0.016}$	222.4 ± 3.7	−0.0138$^{+0.0062}_{-0.0061}$	9.1$^{+3.2}_{-2.2}$	2.316 ± 0.072	−170.9 ± 5.1	Misaligned	14, 15
HAT-P-15b	0.208$^{+0.026}_{-0.025}$	261.8$^{+2.2}_{-2.4}$	−0.0296$^{+0.0076}_{-0.0077}$	−0.206$^{+0.025}_{-0.026}$	185.8$^{+5.2}_{-5.1}$	0.0131$^{+0.0056}_{-0.0059}$	17.3$^{+3.0}_{-2.3}$	2.043$^{+0.082}_{-0.08}$		Misaligned	16
HAT-P-16b	0.0423$^{+0.01}_{-0.0077}$	215$^{+14}_{-21}$	−0.0342$^{+0.0045}_{-0.0041}$	−0.023 ± 0.015	534.1$^{+6.3}_{-6.2}$	0.005$^{+0.011}_{-0.01}$	10.6$^{+5.6}_{-3.1}$	4.22 ± 0.11	−10 ± 16	Misaligned	8, 17
HAT-P-17b	0.342 ± 0.0039	199.1 ± 1.7	−0.3229$^{+0.004}_{-0.0041}$	−0.11 ± 0.01	59.98 ± 0.79	≡ 0.0 ± 0.0f	1.5 ± 0.45	0.58 ± 0.019	$19^{+14}_{-16}$	Misaligned	18, 19
HAT-P-18b	0.106$^{+0.15}_{-0.084}$	12$^{+22}_{-21}$	0.095$^{+0.13}_{-0.082}$	0.008$^{+0.099}_{-0.018}$	25.4$^{+5.8}_{-4.5}$	0.0004$^{+0.0085}_{-0.0077}$	17.5$^{+2.5}_{-2.4}$	0.183$^{+0.034}_{-0.032}$		Control	20
HAT-P-20b	0.0158$^{+0.0041}_{-0.0036}$	327$^{+19}_{-13}$	0.013$^{+0.0023}_{-0.0025}$	−0.0084$^{+0.0053}_{-0.0052}$	1245.4$^{+6.1}_{-6.3}$	−0.0141$^{+0.0073}_{-0.0078}$	14.3$^{+4.5}_{-3.2}$	7.24 ± 0.18		Misaligned	21
HAT-P-22b	0.0064$^{+0.008}_{-0.0046}$	116$^{+180}_{-57}$	0.0002$^{+0.0045}_{-0.004}$	0.0018$^{+0.0099}_{-0.0048}$	314.4 ± 3.2	−0.0147$^{+0.0043}_{-0.0045}$	9.7$^{+2.2}_{-1.6}$	2.157 ± 0.059		Control	21
HAT-P-24b	0.033$^{+0.027}_{-0.021}$	182 ± 63	−0.02$^{+0.019}_{-0.023}$	−0.0004$^{+0.027}_{-0.028}$	86.5 ± 3.6	−0.0099$^{+0.0072}_{-0.0071}$	10.8$^{+3.2}_{-2.6}$	0.715 ± 0.035	20 ± 16	Control	3, 22
HAT-P-26b	0.14$^{+0.12}_{-0.08}$	46$^{+33}_{-71}$	0.075$^{+0.062}_{-0.065}$	0.074$^{+0.15}_{-0.097}$	8.57$^{+0.99}_{-0.97}$	0.002 ± 0.002	3.0$^{+0.74}_{-0.62}$	0.0595$^{+0.0072}_{-0.0071}$		Control	23
HAT-P-29b	0.061$^{+0.044}_{-0.036}$	211$^{+39}_{-65}$	−0.04$^{+0.034}_{-0.031}$	−0.02$^{+0.038}_{-0.057}$	77.6$^{+4.5}_{-4.6}$	0.0498$^{+0.0092}_{-0.01}$	10.8$^{+4.0}_{-2.6}$	0.773$^{+0.052}_{-0.051}$		Control	24
HAT-P-30b	0.02$^{+0.022}_{-0.014}$	114$^{+200}_{-77}$	0.008$^{+0.016}_{-0.01}$	0.002$^{+0.024}_{-0.016}$	89.8$^{+2.7}_{-2.8}$	0.0112$^{+0.0068}_{-0.0071}$	7.7$^{+2.0}_{-1.4}$	0.726 ± 0.027	73.5 ± 9.0	Misaligned	25
HAT-P-31b	0.2419$^{+0.0099}_{-0.0097}$	276.2 ± 1.8	0.0262$^{+0.0076}_{-0.0075}$	−0.2404$^{+0.0097}_{-0.0099}$	231.6$^{+2.5}_{-2.6}$	0.0054$^{+0.0072}_{-0.007}$	6.6$^{+1.7}_{-1.3}$	2.227$^{+0.089}_{-0.09}$		Misaligned	26
HAT-P-32b	0.2$^{+0.19}_{-0.13}$	58$^{+28}_{-53}$	0.076$^{+0.11}_{-0.079}$	0.15$^{+0.19}_{-0.15}$	112$^{+20}_{-21}$	−0.097 ± 0.023	64$^{+11}_{-9}$	0.79 ± 0.15	85 ± 1.5	Misaligned	3, 27
HAT-P-33b	0.13$^{+0.19}_{-0.1}$	15 ± 22	0.114$^{+0.16}_{-0.097}$	0.015$^{+0.13}_{-0.023}$	72$^{+19}_{-16}$	−0.021$^{+0.02}_{-0.023}$	53.5$^{+12.0}_{-8.1}$	0.65 ± 0.14		Control	27
HAT-P-34b	0.411$^{+0.029}_{-0.028}$	17.8 ± 7.2	0.388$^{+0.027}_{-0.028}$	0.124$^{+0.054}_{-0.051}$	364$^{+24}_{-25}$	0.071$^{+0.048}_{-0.05}$	52$^{+14}_{-10}$	3.93 ± 0.28c	0 ± 14	Misaligned	3, 28
HD149026b	0.0028$^{+0.019}_{-0.0024}$	100$^{+170}_{-11}$	−5e-05$^{+0.00036}_{-0.00045}$	0.0005$^{+0.021}_{-0.0025}$	37.9$^{+1.4}_{-1.3}$d	−0.00098$^{+0.00099}_{-0.00089}$	5.13$^{+0.85}_{-0.62}$	0.324 ± 0.011d	12 ± 7	Control	2, 3
TrES-2b	0.0036$^{+0.015}_{-0.0027}$	24$^{+63}_{-110}$	0.00076$^{+0.00053}_{-0.00052}$	0.0002$^{+0.011}_{-0.0061}$	180.1$^{+5.7}_{-5.6}$	−0.0041$^{+0.006}_{-0.0059}$	17.9$^{+4.0}_{-2.8}$	1.157$^{+0.055}_{-0.056}$	−9 ± 12	Control	1, 29
TrES-3b	0.17$^{+0.032}_{-0.031}$	270.5$^{+0.38}_{-0.32}$	0.00151$^{+0.001}_{-0.00098}$	−0.17$^{+0.031}_{-0.032}$	312$^{+13}_{-12}$d	0.08$^{+0.053}_{-0.054}$	104$^{+60}_{-31}$	1.615$^{+0.079}_{-0.077}$d		Misaligned	30
TrES-4b	0.015$^{+0.076}_{-0.012}$	80.6$^{+9.5}_{-160.0}$	0.0012$^{+0.0022}_{-0.0018}$	0.005$^{+0.082}_{-0.011}$	84 ± 10	0.015 ± 0.012	16.5$^{+5.7}_{-3.9}$	0.843$^{+0.098}_{-0.089}$	6.3 ± 4.7	Control	31, 32
WASP-1b	.0082$^{+0.026}_{-0.0072}$	91.1$^{+170.0}_{-6.3}$	3e-05 ± 0.001	0.0068$^{+0.028}_{-0.0072}$	119.6$^{+3.1}_{-3.3}$	0.0029$^{+0.0057}_{-0.0056}$	8.6$^{+3.0}_{-2.3}$	0.79 ± 0.033	−59 ± 99	Control	33, 34
WASP-2b	0.0054$^{+0.009}_{-0.0044}$	267$^{+11}_{-86}$	−0.0001$^{+0.001}_{-0.0011}$	−0.0051$^{+0.0051}_{-0.0092}$	156.7$^{+1.2}_{-1.3}$	0.0062 ± 0.0092	2.75$^{+0.75}_{-0.63}$	0.918$^{+0.027}_{-0.028}$	$153^{+11}_{-15}$	Misaligned	35
WASP-3b	0.0066$^{+0.016}_{-0.0052}$	79$^{+10}_{-130}$	0.0011$^{+0.0014}_{-0.0012}$	0.0052$^{+0.017}_{-0.0059}$	284.0$^{+5.7}_{-6.0}$	−0.0126$^{+0.0091}_{-0.0085}$	15.5$^{+4.5}_{-3.6}$	1.944$^{+0.04}_{-0.042}$	$3.3^{2.5}_{-4.4}$	Control	36, 37
WASP-4b	0.0034$^{+0.0074}_{-0.0026}$	288$^{+140}_{-21}$	0.0006$^{+0.0014}_{-0.0011}$	−0.0019$^{+0.0027}_{-0.0087}$	234.6$^{+2.2}_{-2.3}$	−0.0099$^{+0.0052}_{-0.0054}$	1.91$^{+0.29}_{-0.24}$	1.159$^{+0.063}_{-0.064}$	$-1^{+14}_{-12}$	Control	35, 38

Notes. ^aSystems with accelerations that differ from zero by more than 3σ are marked in bold. ^bThe misaligned sample consists of planets with either eccentric or misaligned orbits, while the control sample contains planets that appear to have circular and/or well-aligned orbits. ^cOur preferred value for this parameter differs from the one in the published literature; this is likely related to our treatment of the stellar jitter and Rossiter data (see Section 3.1). ^dOur value for this parameter differs from the literature, but previous fits were calculated assuming a circular orbit. ^eThe radial velocity acceleration in this system appears to be correlated with the stellar activity, and we therefore conclude that this is probably not the result of an additional companion in this system; see Section 4.1 for more details. ^fBecause the acceleration in the HAT-P-17 system has some curvature, we fit it with a two-planet solution where the linear trend slope term is fixed to zero (see Fulton et al. 2013 for the full solution).

Download table as:  ASCIITypeset image

Table 8. Results from Radial Velocity Fits Continued

Planet	e	ω⋆	ecos ω⋆	esin ω⋆	K	$\dot{\gamma }$a	Jitter	Mp	Spin-Orbit λ	Sampleb	Ref.
		(degrees)			(m s−1)	(m s−1 day−1)	(m s−1)	(MJup)	(deg)
WASP-7b	0.034$^{+0.045}_{-0.024}$	109$^{+170}_{-55}$	0.003$^{+0.026}_{-0.02}$	0.011$^{+0.055}_{-0.025}$	111.3$^{+7.3}_{-7.5}$	0.03 ± 0.04	34.6$^{+4.5}_{-3.8}$	1.131$^{+0.092}_{-0.089}$	86 ± 8	Misaligned	3, 39
WASP-8b	0.3044$^{+0.0039}_{-0.004}$	274.215$^{+0.084}_{-0.082}$	0.02237$^{+0.00032}_{-0.00031}$	−0.304 ± 0.004	221.1 ± 1.2	≡0.0 ± 0.0c	2.91$^{+0.4}_{-0.34}$	2.24$^{+0.11}_{-0.12}$	$-123^{+3.4}_{-4.4}$	Misaligned	40
WASP-10b	0.0473$^{+0.0034}_{-0.0029}$	165.6$^{+9.6}_{-8.6}$	−0.0454$^{+0.0024}_{-0.0023}$	0.0118$^{+0.0076}_{-0.008}$	568.8$^{+7.0}_{-6.7}$	−0.048$^{+0.013}_{-0.012}$	5.4$^{+1.8}_{-1.3}$	3.37 ± 0.11		Misaligned	41
WASP-12b	0.037$^{+0.014}_{-0.015}$	272.7$^{+2.4}_{-1.3}$	0.00171$^{+0.00073}_{-0.00075}$	−0.037$^{+0.015}_{-0.014}$	220.2 ± 3.1	−0.0009$^{+0.0097}_{-0.0093}$	19.5$^{+2.7}_{-2.3}$	1.39 ± 0.13	$59^{+15}_{-20}$	Misaligned	3, 42
WASP-14b	0.0822$^{+0.003}_{-0.0032}$d	251.67$^{+0.64}_{-0.75}$d	−0.02591$^{+0.00049}_{-0.00046}$	−0.078$^{+0.0034}_{-0.0032}$	987.2$^{+1.7}_{-1.8}$	0.0062$^{+0.0044}_{-0.0041}$	5.69$^{+1.1}_{-0.85}$	7.8$^{+0.45}_{-0.47}$	−33.1 ± 7.4	Misaligned	43, 44
WASP-15b	0.038$^{+0.043}_{-0.026}$	240$^{+79}_{-200}$	0.015$^{+0.027}_{-0.018}$	−0.002$^{+0.038}_{-0.044}$	61.8$^{+4.6}_{-4.5}$	0.052$^{+0.047}_{-0.044}$	4.4$^{+2.5}_{-2.3}$	0.566 ± 0.045	$-139.6^{4.3}_{-4.2}$	Misaligned	35
WASP-16b	0.015$^{+0.012}_{-0.011}$	97$^{+44}_{-20}$	−0.0009$^{+0.004}_{-0.0048}$	0.014 ± 0.013	118.9 ± 1.6	0.0056$^{+0.0071}_{-0.0072}$	2.34 ± 0.59	0.846 ± 0.033	$11^{+26}_{-19}$	Control	3, 45
WASP-17b	0.039$^{+0.05}_{-0.027}$	179 ± 120	0.006$^{+0.031}_{-0.021}$	0.0001$^{+0.047}_{-0.044}$	58.8$^{+4.4}_{-4.7}$	0.0002$^{+0.024}_{-0.026}$	11.7$^{+5.0}_{-4.4}$	0.529$^{+0.047}_{-0.048}$	$-148.7^{+7.7}_{-6.7}$	Misaligned	35, 46
WASP-18b	0.0068$^{+0.0025}_{-0.0027}$	261.1$^{+5.3}_{-7.4}$	−0.00104$^{+0.00065}_{-0.00067}$	−0.0067$^{+0.0028}_{-0.0025}$	1816.6$^{+6.1}_{-6.3}$	−0.003$^{+0.0072}_{-0.0077}$	5.1$^{+2.6}_{-1.9}$	10.47$^{+0.49}_{-0.5}$	13 ± 7	Control	3, 35
WASP-19b	0.0024$^{+0.0094}_{-0.0019}$	260$^{+15}_{-170}$	−7e-05$^{+0.00062}_{-0.00066}$	−0.0007$^{+0.0019}_{-0.01}$	254.0$^{+3.4}_{-3.3}$	0.065 ± 0.034	17.8$^{+3.2}_{-2.7}$	1.123 ± 0.036	1.0 ± 1.2	Control	47, 48
WASP-22b	0.0108$^{+0.014}_{-0.0076}$	114$^{+160}_{-56}$	0.0005$^{+0.0079}_{-0.0066}$	0.003$^{+0.018}_{-0.008}$	70.9$^{+1.5}_{-1.6}$	0.0583$^{+0.0078}_{-0.0074}$	7.2$^{+1.7}_{-1.4}$	0.569$^{+0.016}_{-0.015}$	22 ± 16	Control	49, 50
WASP-24b	0.0033$^{+0.012}_{-0.0026}$	70$^{+20}_{-150}$	0.0005$^{+0.00086}_{-0.0007}$	0.0011$^{+0.014}_{-0.0027}$	152.0 ± 3.2	−0.062 ± 0.051	3.65$^{+0.89}_{-0.8}$	1.119 ± 0.029	−4.7 ± 4	Control	7, 33
WASP-34b	.0109$^{+0.015}_{-0.0078}$	215$^{+77}_{-140}$	−0.0001$^{+0.0068}_{-0.0071}$	−0.001$^{+0.011}_{-0.017}$	71.1$^{+1.6}_{-1.7}$	≡0.0 ± 0.0c	3.2$^{+0.72}_{-0.6}$	0.57 ± 0.03		Control	51
WASP-38b	0.0329$^{+0.01}_{-0.0086}$	17$^{+27}_{-33}$	0.0284$^{+0.0076}_{-0.0078}$	0.008$^{+0.018}_{-0.016}$	252.1$^{+4.4}_{-4.3}$	−0.074 ± 0.058	11.9$^{+2.4}_{-1.9}$	2.705 ± 0.076	$7.5^{+4.7}_{-6.1}$	Misaligned	52
XO-2b	0.028$^{+0.038}_{-0.022}$	261$^{+11}_{-71}$	−0.0037$^{+0.0052}_{-0.0063}$	−0.027$^{+0.027}_{-0.039}$	93.9$^{+2.1}_{-2.2}$	0.0126$^{+0.0039}_{-0.0036}$	9.3$^{+2.4}_{-1.9}$	0.629 ± 0.017	10 ± 72	Control	1, 53
XO-3b	0.2833 ± 0.0034	346.8$^{+1.6}_{-1.5}$	0.2756 ± 0.0027	−0.0649$^{+0.0081}_{-0.008}$	1480 ± 11	−0.019$^{+0.025}_{-0.024}$	43.5$^{+8.3}_{-7.0}$	12.15 ± 0.48	37.3 ± 3.0	Misaligned	54, 55
XO-4b	0.002$^{+0.012}_{-0.002}$	240$^{+39}_{-160}$	0.00016$^{+0.00062}_{-0.00051}$	−0.0001$^{+0.0039}_{-0.0088}$	163.7 ± 4.7	0.01 ± 0.01	7.3$^{+2.4}_{-1.9}$	1.559$^{+0.052}_{-0.048}$	−46.7 ± 8.1	Misaligned	56
XO-5b	0.013$^{+0.014}_{-0.009}$	184 ± 92	−0.0031$^{+0.0073}_{-0.013}$	−0.0001$^{+0.012}_{-0.013}$	144.3$^{+2.9}_{-3.0}$	0.0041 ± 0.0029	10.7$^{+2.3}_{-1.8}$	1.051 ± 0.032		Control	57

Notes. ^aSystems with accelerations that differ from zero by more than 3σ are marked in bold. ^bThe misaligned sample consists of planets with either eccentric or misaligned orbits, the control sample contains planets that appear to have circular and/or well-aligned orbits. ^cBecause the accelerations in the WASP-8 and WASP-34 systems have some curvature, we fit then with a two-planet solution where the linear trend slope term is fixed to zero (see Tables 9 and 10 for the full solution). ^dOur values differ from those of previous fits, which did not include the measured secondary eclipse times. References for orbital inclinations and spin-orbit angles. (1) Torres et al. 2008; (2) Carter et al. 2009; (3) Albrecht et al. 2012b; (4) Pál et al. 2010; (5) Winn et al. 2011; (6) Winn et al. 2009a; (7) Simpson et al. 2011; (8) Moutou et al. 2011; (9) West et al. 2009b; (10) Bakos et al. 2010; (11) Winn et al. 2010c; (12) Hartman et al. 2009; (13) Winn et al. 2010b; (14) Torres et al. 2010; (15) Winn et al. 2011; (16) Kovács et al. 2010; (17) Buchhave et al. 2010; (18) Howard et al. 2012; (19) Fulton et al. 2013; (20) Hartman et al. 2011b; (21) Bakos et al. 2011; (22) Kipping et al. 2010; (23) Hartman et al. 2011a; (24) Buchhave et al. 2011; (25) Johnson et al. 2011; (26) Kipping et al. 2011; (27) Hartman et al. 2011c; (28) Bakos et al. 2012; (29) Winn et al. 2008b; (30) Sozzetti et al. 2009; (31) Mandushev et al. 2007; (32) Narita et al. 2010a; (33) Simpson et al. 2011; (34) Albrecht et al. 2011; (35) Triaud et al. 2010; (36) Gibson et al. 2008; (37) Tripathi et al. 2010; (38) Sanchis-Ojeda et al. 2011; (39) Hellier et al. 2008; (40) Queloz et al. 2010; (41) Johnson et al. 2009b; (42) Maciejewski et al. 2011a; (43) Joshi et al. 2009; (44) Johnson et al. 2009a; (45) Lister et al. 2009; (46) Anderson et al. 2010; (47) Hellier et al. 2011; (48) Tregloan-Reed et al. 2012; (49) Maxted et al. 2010; (50) Anderson et al. 2011; (51) Smalley et al. 2011; (52) Brown et al. 2012b; (53) Narita et al. 2011; (54) Johns-Krull et al. 2008; (55) Hirano et al. 2011; (56) Narita et al. 2010a; (57) Maciejewski et al. 2011b.

Download table as:  ASCIITypeset image

While vetting the fits for all planets we noticed that in some cases the two dimensional marginalized distributions when plotted in $\sqrt{e}\sin {\omega _{\star }}$ versus $\sqrt{e}\cos {\omega _{\star }}$ took on non-Gaussian shapes in systems for which we had many secondary eclipse times to constrain the orbital phase of the secondary eclipse (and thus ecos ω_⋆). Systems for which we had good secondary eclipse priors but the eccentricity was low took on a star-shaped appearance while a strong correlation between $\sqrt{e}\sin {\omega _{\star }}$ and $\sqrt{e}\cos {\omega _{\star }}$ emerged in systems with significant eccentricity. This correlation was much less pronounced when we plotted ecos ω_⋆ versus esin ω_⋆ (Figure 3). In order to check that our DE-MCMC algorithm was behaving as expected we created hypothetical distributions of e and ω for two cases.

Zoom In Zoom Out Reset image size

First, we created a hypothetical distribution of e as the absolute value of a Gaussian centered around zero and a uniform distribution of ω between 0 and 2π and to simulate the posterior distributions for a planet with no significant eccentricity. We then extracted the points within the distributions for which the orbital phase of the secondary eclipse calculated from e and ω was very close the median value of all secondary eclipse times calculated from the e and ω distributions (near phase = 0.5 in this case). When the entire distribution is plotted in $\sqrt{e}\sin {\omega _{\star }}$ versus $\sqrt{e}\cos {\omega _{\star }}$ we see a smooth distribution with circular contours, but when the points extracted based on the secondary eclipse times are plotted we see the star shape that closely resembles the distributions we obtain in our fits to the data.

Second, we created a hypothetical distribution of e as a normal distribution centered around a value of 0.16 with a width of 0.02 and a normal distribution of ω centered around 328 degrees with a width of 10°. These distributions are meant to mimic a planet with significant eccentricity. When points are extracted based on the secondary eclipse times in the same way as the first case we see that the points fall along a locus that matches the distributions obtained in our fits.

This test shows that our DE-MCMC algorithm is working as expected, but that the choice of parameterization to use $\sqrt{e}\sin {\omega _{\star }}$ and $\sqrt{e}\cos {\omega _{\star }}$ may not be optimal for systems with secondary eclipse measurements. However, since our DE-MCMC code continuously checks for convergence we know that the chains are converged and well-mixed. The correlated parameters will slow the convergence, but we decided that the factor of five increase in runtime was acceptable and did not change parameterization.

We use our K band NIRC2 imaging data to place upper limits on the allowed masses and orbital semi-major axes of the perturbers responsible for the measured radial velocity accelerations. Contrast curves are generated for each target as follows. First, we calculate the FWHM of the central star's point spread function in the interpolated and combined image. The maximum radius for our contrast curves is defined as the largest radial separation for which data is available at all position angles (i.e., we do not count the corners of the array). We then create a box with dimensions equal to the FWHM and step it across the array, calculating the total flux from the pixels within the box at a given position. We exclude boxes containing masked pixels¹¹ and boxes whose radial distance from the star is greater than our maximum radius limit. The 5σ contrast limit is calculated as a function of radial separation from the star by taking the standard deviation of the total flux values for boxes within a given annulus with width equal to the full width at half max of the stellar probability density function (i.e., one box width) and multiplying by five. We convert our absolute flux limits to relative delta magnitude units by taking the maximum flux value in the interpolated stellar point spread function as an estimate of the flux of the central star and calculating the corresponding relative magnitude limits for each radial distance. We show the resulting contrast curves for all of our targets in Figure 4.

Zoom In Zoom Out Reset image size

With the exception of GJ 436 and HAT-P-2, none of our target stars have directly measured parallax estimates. In most cases the discovery paper provides an estimate of the stellar properties (mass, radius, and age) from fitting stellar evolution models using constraints on the surface gravity, effective temperature, and metallicity from high-resolution optical spectroscopy and (in some cases) constraints on the stellar density from fits to the transit light curve. The distance can then be estimated using the known stellar properties and the measured apparent magnitudes in V, J, H, and K bands. We take these estimated distances and use them to convert the units of our contrast curves from separations in arc seconds to projected physical distances in AU

We convert our contrast curves from delta magnitudes in either K_s or K_p bands to mass limits for stellar companions using the latest version of the PHOENIX stellar atmosphere models (Husser et al. 2013). We assume solar metallicities for both the primary and secondary, and interpolate in the available grid of models to produce a model that exactly matches the effective temperatures and surface gravities of each star. We utilize the published temperatures and surface gravities for our primary stars, taking the best available constraints in each case. We then systematically step through the table of radius and effective temperature as a function of secondary mass for a low-mass main-sequence companion from Baraffe et al. (1998) and create matching PHOENIX models for a corresponding secondary stellar companion with those properties. The corresponding contrast ratio between the primary and secondary as a function of mass is calculated by integrating over the appropriate bandpass (either K_p or K_s). Finally, we convert our contrast curves from units of delta magnitude to secondary mass using the mass versus delta magnitude relations derived for that system.

Our approach differs from the standard approach for AO imaging searches for stellar companions (e.g., Bechter et al. 2013), which typically utilize relative K magnitude estimates from 2MASS and parallax measurements to calculate an absolute K magnitude for the primary and then interpolate in a grid of absolute K magnitudes as a function of secondary mass calculated from standard stellar evolution models at a given age (e.g., Girardi et al. 2002). Our method offers two advantages over this approach: first, we do not need a distance estimate to calculate the contrast ratio between the primary and secondary, and second, we can calculate contrast ratios in arbitrary bandpasses as needed. We validate our method by converting the K band contrast curves for HAT-P-8 and WASP-12b from Bechter et al. (2013) using our new method, and find results that are consistent to within 0.02 solar masses in both cases.

We find linear or curved trends in the measured radial velocities with slopes at least 3σ away from zero for fifteen systems listed in Tables 7 and 8. We next checked these systems to determine if any of the radial velocity trends were well-correlated with the stellar Ca ii H & K emission index S_HK. We find that one system, HAT-P-11, does exhibit a correlation with the measured S_HK and therefore conclude that this signal is likely the result of stellar activity rather than a real companion (see Figure 5). This is not surprising, as this is one of the most active stars in our sample with a log (R'_HK) value of −4.57 (Knutson et al. 2010). Of the remaining fourteen systems with evidence for an outer companion, eight have previously been reported in the published literature including: HAT-P-2 (Lewis et al. 2013), HAT-P-4 (Winn et al. 2011), HAT-P-7 (Winn et al. 2009a), HAT-P-13 (Bakos et al. 2009a; Winn et al. 2010b), HAT-P-17 (Howard et al. 2012; Fulton et al. 2013), WASP-8 (Queloz et al. 2010), WASP-22 (Maxted et al. 2010; Anderson et al. 2011), and WASP-34 (Smalley et al. 2011). We present a composite plot showing all of the detected radial velocity accelerations in Figure 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We also report new trend detections for six systems including: HAT-P-10, HAT-P-22, HAT-P-29, HAT-P-32, WASP-10, and XO-2. Finally, we do not find statistically significant accelerations in the following systems with previously reported trend detections: GJ 436 (Maness et al. 2007), HAT-P-31 (Kipping et al. 2011), and HAT-P-34 (Bakos et al. 2012). We discuss the differences between our results and those of previous studies for individual systems in Section 4.1.1 below.

4.1.1. Comparison to Previously Studies

4.1.2. Fits to Systems with Curved Radial Velocity Trends or Multiple Planets

We find three systems with evidence for curvature in the radial velocity acceleration, including: HAT-P-17, WASP-8, and WASP-34. For HAT-P-13, the second planet has a fully resolved orbit and there is an additional linear acceleration present. Our fits for HAT-P-17 follow the methods described in Fulton et al. (2013) and we obtain values that are consistent with that study; this is not surprising, as we have added just two new radial velocity measurements in our fits. We discuss our results for the other three systems individually below.

Queloz et al. (2010) reported a linear trend with a slope of 58.1 ± 1.3 m s⁻¹ yr⁻¹ for the WASP-8 system. We find that this trend has turned over in our new observations, allowing us to fit for the orbital properties of the outer companion rather than assuming a linear trend. We show our results from these new fits in Table 9 and Figure 8. We find that when we assume a circular orbit for the outer companion we obtain a reasonably well-constrained orbital solution with a period of $4339^{+850}_{-390}$ days, a radial velocity semi-amplitude K of $115.7^{+21.0}_{-9.4}$, and Msin i equal to 9.5$^{+2.7}_{-1.1}\,M_{\rm J}$. When we allow the eccentricity to vary freely in the fits our chains do not converge on a well-defined solution, and both the period and radial velocity semi-amplitude span a larger range in values ($3971^{+1100}_{-900}$ days and $110^{+24}_{-27}$ m s⁻¹, respectively). We present the results for the better-constrained circular fit in Table 9 and consider the non-circular case in more detail in Section 4.3; we note that the data are equally consistent with both circular and eccentric fits.

Table 9. Fit Parameters for WASP-8 System

Parameter	Value	Units
RV step parameters
log(Pb)	0.91162223$^{+7.6e-07}_{-7.8e-07}$	log(days)
Tc, b	2454679.33392 ± 0.00047	BJDTDB
$\sqrt{e_{b}}\cos {\omega _{b}}$	0.04055 ± 0.00064
$\sqrt{e_{b}}\sin {\omega _{b}}$	−0.5502$^{+0.0037}_{-0.0036}$
log(Kb)	2.3446 ± 0.0023	m s−1
log(Pc)	3.636$^{+0.068}_{-0.039}$	log(days)
Tc, c	2452613$^{+330}_{-610}$	BJDTDB
$\sqrt{e_{c}}\cos {\omega _{c}}$	≡0.0 ± 0.0
$\sqrt{e_{c}}\sin {\omega _{c}}$	≡0.0 ± 0.0
log(Kc)	2.061$^{+0.062}_{-0.036}$	m s−1
γ1	−57.9$^{+9.6}_{-17.0}$	m s−1
γ2	−20$^{+24}_{-37}$	m s−1
γ3	0$^{+24}_{-38}$	m s−1
$\dot{\gamma }$	≡0.0 ± 0.0	m s−1 day−1
Jitter	2.91$^{+0.4}_{-0.34}$	m s−1
RV model parameters
Pb	8.158724$^{+1.4e-05}_{-1.5e-05}$	days
Tc, b	2454679.33392 ± 0.00047	BJDTDB
eb	0.3044$^{+0.0039}_{-0.004}$
ωb	274.215$^{+0.084}_{-0.082}$	deg
Kb	221.1 ± 1.2	m s−1
Pc	4323$^{+740}_{-380}$	days
Tc, c	2452613$^{+330}_{-610}$	BJDTDB
ec	≡0.0 ± 0.0
ωc	≡90.0 ± 0.0	deg
Kc	115.0$^{+18.0}_{-9.2}$	m s−1
γ1	−57.9$^{+9.6}_{-17.0}$	m s−1
γ2	−20$^{+24}_{-37}$	m s−1
γ3	0$^{+24}_{-38}$	m s−1
$\dot{\gamma }$	≡ 0.0 ± 0.0	m s−1 day−1
Jitter	2.91$^{+0.4}_{-0.34}$	m s−1
RV derived parameters
ecos ω	0.02237$^{+0.00032}_{-0.00031}$
esin ω	−0.304 ± 0.004
Mcsin ic	9.45$^{+2.26}_{-1.04}$	MJ
ac	5.28$^{+0.63}_{-0.34}$	AU

Download table as:  ASCIITypeset image

Our new measurements also indicate that the linear trend reported in Smalley et al. (2011) for WASP-34 has turned over, allowing us to place weak constrains on the orbital period and mass of the companion for the case of a circular orbit (see Table 10 and Figure 9). We find that in these fits the companion has an orbital period of $4133^{+2500}_{-1400}$ days, a radial velocity semi-amplitude K of $196^{+280}_{-98}$, and Msin i equal to 15.5$^{+28.0}_{-8.8}\,M_{\rm J}$. We evaluate the constraints for the non-circular case separately in Section 4.3.

Table 10. Fit Parameters for WASP-34 System

Parameter	Value	Units
RV step parameters
log(Pb)	0.63525024 ± 4.5e−07	log(days)
Tc, b	2454647.55434$^{+0.00063}_{-0.00064}$	BJDTDB
$\sqrt{e_{b}}\cos {\omega _{b}}$	−0.002$^{+0.062}_{-0.064}$
$\sqrt{e_{b}}\sin {\omega _{b}}$	−0.02 ± 0.11
log(Kb)	1.8517$^{+0.0097}_{-0.01}$	m s−1
log(Pc)	3.612$^{+0.073}_{-0.059}$	log(days)
Tc, c	2454586$^{+140}_{-190}$	BJDTDB
$\sqrt{e_{c}}\cos {\omega _{c}}$	≡0.0 ± 0.0
$\sqrt{e_{c}}\sin {\omega _{c}}$	≡0.0 ± 0.0
log(Kc)	2.28$^{+0.12}_{-0.09}$	m s−1
γ1	108$^{+62}_{-37}$	m s−1
γ2	141$^{+62}_{-37}$	m s−1
$\dot{\gamma }$	≡0.0 ± 0.0	m s−1 day−1
Jitter	3.2$^{+0.72}_{-0.6}$	m s−1
RV model parameters
Pb	4.3176779 ± 4.5e−06	days
Tc, b	2454647.55434$^{+0.00063}_{-0.00064}$	BJDTDB
eb	0.0109$^{+0.015}_{-0.0078}$
ωb	215$^{+77}_{-140}$	deg
Kb	71.1$^{+1.6}_{-1.7}$	m s−1
Pc	4093$^{+750}_{-520}$	days
Tc, c	2454586$^{+140}_{-190}$	BJDTDB
ec	≡0.0 ± 0.0
ωc	≡90.0 ± 0.0	deg
Kc	189$^{+60}_{-35}$	m s−1
γ1	108$^{+62}_{-37}$	m s−1
γ2	141$^{+62}_{-37}$	m s−1
$\dot{\gamma }$	≡0.0 ± 0.0	m s−1 day−1
Jitter	3.2$^{+0.72}_{-0.6}$	m s−1
RV derived parameters
ecos ω	−0.0001$^{+0.0068}_{-0.0071}$
esin ω	−0.001$^{+0.011}_{-0.017}$
Mcsin ic	14.96$^{+6.29}_{-3.39}$	MJ
ac	5.05$^{+0.65}_{-0.46}$	AU

Download table as:  ASCIITypeset image

HAT-P-13 presents a particularly interesting case, as our fits indicate evidence for two outer companions in the system. Previous studies (Bakos et al. 2009a; Winn et al. 2010b) reported the presence of one companion with a fully resolved orbit ("c") and an additional radial velocity trend ("d"). We provide an updated estimate for the properties of companion "c," which has a period of 445.87 ± 0.12 days and an orbital eccentricity of 0.6573 ± 0.0034. We estimate a Msin i of 14.70$^{+0.48}_{-0.47}\,M_{\rm J}$ for this companion, and list the full set of fit parameters in Table 11. Dynamical studies of this system (Batygin et al. 2009; Becker & Batygin 2013) predict that this planet will perturb on the orbit of the inner hot Jupiter ("b"), resulting in the alignment of the apses of the two planetary orbits. In this scenario, the eccentricity of the inner planet can be used to constrain the inner planet's tidal Love number (a measure of its degree of central concentration). Our new fits indicate that the arguments of periapse ω for the orbits of the two inner planets are consistent but with large uncertainties on ω_b, which are primarily due to this planet's small orbital eccentricity (also see Winn et al. 2010b). We show the phased radial velocity curves in Figure 7, and provide additional constraints on the properties of companion "d" in Section 4.3. When we include a second planet in our fit to HAT-P-17 we obtain the same orbital parameters as those reported in Fulton et al. (2013), with no evidence for any additional accelerations in this system.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 11. Fit Parameters for HAT-P-13 System

Parameter	Value	Units
RV step parameters
log(Pb)	0.46482299$^{+3.2e-07}_{-3.3e-07}$	log(days)
Tc, b	2455176.53877 ± 0.00027	BJDTDB
$\sqrt{e_{b}}\cos {\omega _{b}}$	−0.096$^{+0.027}_{-0.021}$
$\sqrt{e_{b}}\sin {\omega _{b}}$	−0.031$^{+0.069}_{-0.059}$
log(Kb)	2.0248 ± 0.0032	m s−1
log(Pc)	2.64916$^{+0.00011}_{-0.0001}$	log(days)
Tc, c	2455311.82 ± 0.19	BJDTDB
$\sqrt{e_{c}}\cos {\omega _{c}}$	−0.8068 ± 0.0013
$\sqrt{e_{c}}\sin {\omega _{c}}$	0.0649 ± 0.0031
log(Kc)	2.631 ± 0.0021	m s−1
γ	−23.04$^{+0.84}_{-0.86}$	m s−1
$\dot{\gamma }$	0.0528$^{+0.0013}_{-0.0014}$	m s−1 day−1
Jitter	4.53$^{+0.54}_{-0.46}$	m s−1
RV model parameters
Pb	2.9162381$^{+2.1e-06}_{-2.2e-06}$	days
Tc, b	2455176.53877 ± 0.00027	BJDTDB
eb	0.0133$^{+0.0047}_{-0.0044}$
ωb	197$^{+32}_{-37}$	deg
Kb	105.87 ± 0.78	m s−1
Pc	445.82 ± 0.11	days
Tc, c	2455311.82 ± 0.19	BJDTDB
ec	0.6551 ± 0.0021
ωc	175.40 ± 0.22	deg
Kc	427.6 ± 2.1	m s−1
γ	−23.04$^{+0.84}_{-0.86}$	m s−1
$\dot{\gamma }$	0.0528$^{+0.0013}_{-0.0014}$	m s−1 day−1
Jitter	4.53$^{+0.54}_{-0.46}$	m s−1
RV derived parameters
ecos ω	−0.0107$^{+0.0039}_{-0.0041}$
esin ω	−0.0032$^{+0.0073}_{-0.0077}$
Mcsin ic	14.61$^{+0.46}_{-0.48}$	MJ
ac	1.258 ± 0.020	AU

Download table as:  ASCIITypeset image

Four of our targets with detected radial velocity accelerations also host candidate AO-detected companions; in this section we consider whether such companions could explain the observed radial velocity trend. For the cases where estimated spectral types are not available, we estimate the masses of the companions based on their brightness in K band relative to the primary using the same methods described in Section 3.3. We convert their projected separations on the sky to a minimum semi-major axis using the estimated distances of these systems, and then compare the estimated masses and minimum separations to the lower limit on the companion mass at that separation from the measured radial velocity trend. We calculate this lower limit using the following expression (Torres 1999; Liu et al. 2002):

$$\begin{eqnarray} M_{\rm comp} &=& 5.34\times 10^{-6} M_{\odot } \left(\frac{d}{\rm pc}\frac{\rho }{\rm arcsec}\right)^2 \nonumber \\ & & \times \left|\frac{\dot{v}}{\rm m\phantom{0}s^{-1}\phantom{0}yr^{-1}}\right|F(i,e,\omega,\phi) \end{eqnarray} \tag{ 1 }$$

where d is the distance to the system, ρ is the projected separation of the companion on the sky, $\dot{v}$ is the best-fit radial velocity trend, and F(i, e, ω, ϕ) is an equation that depends on the unknown orbital parameters of the companion. Liu et al. find that this equation has a minimum value of $\sqrt{27}/2$, which we use in our calculations here.

HAT-P-7 is known to have a common proper motion companion with a projected separation of 39 and spectral type of M5.5 (Narita et al. 2012), corresponding to a mass of approximately 0.2 M_☉. We convert this angular separation to a physical separation of 1240 AU using the estimated distance of $320^{+50}_{-40}$ pc from Pál et al. (2008). At this distance the minimum mass of the companion required to produce the measured radial velocity trend of 25.4  ±  1.3 m s⁻¹ yr⁻¹ is 540 M_☉. We therefore conclude that the observed companion cannot be responsible for the trend in this system, in agreement with the conclusion of Narita et al. (2012).

Adams et al. (2013) reported the discovery of a candidate companion to HAT-P-32 with a projected separation of 29 and a delta magnitude of 3.4 in the K_s band. If we assume that this is a bound companion at the same distance as the primary, we find an estimated mass of 0.4 M_☉. We take our distance estimate of 285 ± 5 pc from Hartman et al. (2011c), and calculate a corresponding physical separation of 830 AU for the companion. We compare this value this to the minimum mass of 318 M_☉ required to explain the radial velocity trend slope of −33 ± 10 m s⁻¹ yr⁻¹ at this separation and conclude that the companion cannot be responsible for the measured trend in this system.

Queloz et al. (2010) report a common proper motion companion to WASP-8 with a sky-projected separation of 483 ± 001 and a relative K magnitude of 2.1 from 2MASS photometry. Assuming a distance of 87 ± 7 pc (Queloz et al. 2010), we find that this companion has a physical separation of 390 AU and an estimated mass of 0.5 M_☉. This is much less than the minimum mass of 125 M_☉ needed to explain the radial velocity trend of $58.1^{+1.2}_{-1.3}$ m s⁻¹ yr⁻¹ from Queloz et al. We note that our new radial velocity data show a downward slope, indicating that the radial velocity trend reached a maximum in the past few years; this provides additional support for the hypothesis that the radial velocity trend is caused by a third, close-in body in the system.

Our preliminary K band AO imaging also resulted in the detection of a previously unknown companion to HAT-P-10 with a sky-projected separation of 034 and a relative magnitude of 2.4 in the K_p band (H. Ngo et al. in preparation). Assuming a distance of 122 ± 4 pc from Bakos et al. (2009b), we find that this corresponds to a sky-projected separation of 42 AU and a companion mass of 0.36 M_☉. This is easily consistent with the minimum mass of 0.12 M_☉ needed to explain the measured radial velocity trend of −5.1 ± 1.4 m s⁻¹ yr⁻¹; we therefore conclude that HAT-P-10 is the only system where an AO companion might explain the presence of the radial velocity trend, and we exclude this system from our subsequent analysis of the frequency of substellar companions in Section 4.4.

We next simulate RV observations of each of our stars to determine what constraints we may place on the properties of the companions responsible for the observed radial velocity accelerations. We refer to these plots as "Wright diagrams" Wright et al. (2007). For each star, we develop a logarithmically spaced 50 × 50 grid of possible companion masses and semi-major axes spanning the range 0.2 M_J < msin i < 500 M_J and 1 AU < a < 75 AU. At each mass and semi-major axis, we inject a simulated planet with a fixed eccentricity and determine the orbital parameters which allow for the best fit to the RV observations. We calculate a χ² value at each grid point for each eccentricity simulated, assuming our RV uncertainties (calculated as the quadrature sum of the reported errors and the best-fit jitter value) are random, uncorrelated, and Gaussian. We convert these likelihood values to a normalized probability, then marginalize over eccentricity. Here, we assume the long-period giant planet eccentricity distribution is well-replicated by the beta distribution (Kipping 2013):

$$\begin{equation} P_{\beta }(e;a,b)= \frac{\Gamma (a+b)}{\Gamma (a) \Gamma (b)} e^{a-1}(1-e)^{b-1} \end{equation} \tag{ 2 }$$

where P_β is the probability of a given eccentricity, Γ is a Gamma function, and a and b are constants that are fitted to the known population of long-period giant planets (a = 1.12 and b = 3.09 here).

If the trend is truly linear, we would not expect to break the degeneracy between the companion mass and semi-major axis, as the two are degenerate. For a given trend, the mass of a companion is proportional to the square of the companion's semi-major axis. In these cases, from the RV observations alone we can only place a lower limit on the mass and separation of a companion, as the orbit must be significantly longer than the RV baseline in order to produce a strictly linear trend. We place an outer limit on the companion's orbit using the contrast curves derived from our K-band AO imaging described in Section 3.3, by injecting and "imaging" artificial companions in the same manner as described by Montet et al. (2013).

A non-detection does not necessarily imply that the star does not host a giant companion; the companion may simply be too small or distant to induce a detectable RV acceleration. We can determine the likelihood of detecting a companion as a function of mass and semi-major axis. To accomplish this, we simulate planets over the grid described previously. For all stars, we inject 1000 planets at each mass and semi-major axis included in our grid. We randomly assign all other parameters following the distributions described above. We then "observe" each star by integrating the companion's orbit and calculating the magnitude of the RV signal at the times of our RV observations. Each velocity is then perturbed from the expected value by a normal variate with zero mean and standard deviation σ equal to the RV uncertainty. If the best-fit to the RV acceleration is 3σ different from zero, we consider this companion to be detected.

The end result of this analysis is a map (the Wright diagram) showing either the range of companions in mass versus semi-major axis space that are consistent with the measured radial velocity acceleration, or a map showing the region of this space where companions are excluded by the current measurements (Figure 10). Although we also presented separate two-planet fits for the HAT-P-13, HAT-P-17 (see Fulton et al. 2013), WASP-8, and WASP-34 systems, in the case of the latter two systems we assumed that the outer companion had a circular orbit in order to avoid degeneracies between the companion's eccentricity, mass, and semi-major axis in our fits. Our new maps (Figure 10) allow these partially resolved companions to have a non-zero orbital eccentricity, allowing for a broader range of orbital solutions. We do not calculate a map for the middle ("c") companion in the HAT-P-13 system, as this planet has a fully resolved orbit with well-constrained properties listed in Table 11. We list the constraints on each companion from this analysis in Table 12.

Zoom In Zoom Out Reset image size

We have determined the parameter space in mass and semi-major axis where a companion could reside for each of our fifteen systems exhibiting RV accelerations corresponding to companions whose orbits are not yet fully resolved. We can combine this information to determine the most likely underlying giant companion distribution for systems hosting short-period gas giant planets.

We assume giant planets are distributed in planet mass and semi-major axis subject to the double power law

$$\begin{equation} f(m,a) = C m^\alpha a^\beta d \ln m \, d \ln a. \end{equation} \tag{ 3 }$$

The likelihood of the data for a star with an RV trend detection is

$$\begin{equation} L_i = \int d \ln m \int d \, \ln a \, f(m,a) \, p_i(m,a), \end{equation} \tag{ 4 }$$

where p_i(m, a) is the probability of a planet at mass m and orbital semi-major axis a, as calculated using the technique of the previous section. For HAT-P-13c we have a well constrained orbit, and we approximate p_i(m, a) ≈ δ(m_i, a_i). In this case, the likelihood of the data is simply L_i = f(m_i, a_i).

For a non-detection, we are able to rule out high-mass, close-in companions, subject to our detectability simulations (D(m, a)) of the previous section. Since it remains possible that the star has a companion below our detectability limit, the likelihood of the data given a non-detection is

$$\begin{equation} L_i = 1 - \int d \ln m \, d \ln a \, D_i(m,a) \, f(m,a). \end{equation} \tag{ 5 }$$

Therefore, the total likelihood for a system of N_d detections around N stars is

$$\begin{eqnarray} \mathcal {L} &=& \prod _{i=1}^{N_d} \left (\int d \ln m \, \int d \ln a \, f(m,a) \, p_i(m,a)\right) \nonumber \\ & & \times \prod _{j=1}^{N-N_d} \left (1 - \int d \ln m \, d \ln a \, D_j(m,a) \, f(m,a)\right). \end{eqnarray} \tag{ 6 }$$

We then vary α, β, and our normalization factor C to maximize $\mathcal {L}$. This is similar to the approach taken by Cumming et al. (2008), although with their well-characterized planets, p_i was approximated by these authors as a δ function. Here, we can only assume a δ function in mass and semi-major axis for the few systems with well-characterized orbits. This approach is also functionally identical to injecting artificial planets following some distribution and matching the observed distribution to the simulated planets, in the limit as the number of injected planets approaches infinity. We calculate constraints on C, α, and β using our likelihood function described above and the emcee package developed by Foreman-Mackey (2013). The distribution of acceptable values of C, α, and β are shown in Figure 11 and listed in Table 13.

Zoom In Zoom Out Reset image size

Table 13. Maximally Likely Values for Distribution Parameters^a

Parameter	Max $\mathcal {L}$	1σ	2σ	3σ
All targets
C	7.6 × 10−4	[1.2 × 10−3, 0.013]	[2.3 × 10−4, 0.031]	[3.8 × 10−5, 0.062]
α	1.7	[0.5, 1.6]	[0.1, 2.3]	[−0.2, 3.2]
β	0.9	[0.3, 1.0]	[−0.1, 1.4]	[−0.4, 1.8]
"Misaligned" b
C	1.2 × 10−3	[0.0016, 0.018]	[3.4 × 10−4, 0.044]	[7.0 × 10−5, 0.09]
α	1.8	[0.6, 1.7]	[0.1, 2.4]	[−0.2, 3.1]
β	0.6	[−0.2, 0.6]	[−0.7, 1.1]	[−1.0, 1.5]
"Control" b
C	6.2 × 10−4	[0.0036, 0.051]	[5.1 × 10−4, 0.12]	[7.3 × 10−5, 0.24]
α	1.0	[−0.5, 0.6]	[−1.0, 1.5]	[−1.4, 2.5]
β	1.3	[0.1, 1.3]	[−0.5, 1.9]	[−1.1, 2.6]

Notes. ^aC is a normalization factor which depends on the total giant planet occurrence rate, α is the power law exponent for frequency as a function of mass, and β is the power law exponent for frequency as a function of semi-major axis. Because some distributions are significantly skewed, the maximum likelihood value is not included inside the 1σ confidence interval for all parameters. ^bSee Table 1 for a list of the systems included in the "Misaligned" and "Control" samples, and Section 4.4 for a description of how the two samples were devised.

Download table as:  ASCIITypeset image

An estimated value of the planet occurrence rate can be found by integrating f over a domain of interest. If we take a range of 1–13 M_Jup in mass and 1–20 AU in orbital semi-major axes, we find a total occurrence rate of 51% ± 10% for our sample. This suggests that there may be an additional 13 ± 5 companions in this range that were missed by our observations. If we consider a smaller range of 1–10 AU with the same mass range, the total occurrence rate is 27% ± 6%, which is lower than in our previous example but also more tightly constrained. Finally, if we consider a broader range of 0.2–13 M_Jup and 1–20 AU we estimate a higher occurrence rate of $55^{+12}_{-10}\%$ with modestly increased uncertainties. One caveat to these integrated occurrence rates is that our choice of a power law distribution may not remain accurate from Saturn mass objects all the way up to the deuterium burning limit. We chose this formalism and domain in order to facilitate comparison to previous analyses of planet occurrence rates, which have often made the same assumptions.

We compare our "misaligned" sub-sample, which contains planets with misaligned and/or eccentric orbits, to our "control" sub-sample of well-aligned planets on circular orbits. We find that the companion occurrence rates in the two sub-samples are consistent within the uncertainties, suggesting that the spin-orbit alignments and eccentricities of the short-period planets are not affected by the presence of these massive companions. All well-characterized outer planets exist in our misaligned sub-sample, so it is perhaps not surprising that the companion distribution parameters are better constrained in this sub-sample than in the control sub-sample.

We investigate the consistency of the two samples by calculating the number of RV trends we would expect to observe if the giant planets in both samples were represented by the underlying planet distribution of the "misaligned" sub-sample. For each value of C, α, and β, we calculate the number of trends we would expect to observe, given the assumed underlying planet population and our calculated ability to detect companions around each star as a function of companion mass and semimajor axis. We then weight each value according to the relative likelihood of that particular choice of parameters. In the full sample, we would expect to detect 14.4 ± 3.2 companions based on our planet distribution function (Figure 12). If the full sample of stars is described by the parameterizations of the "misaligned" (or "control") sub-sample, then we would expect to detect $14.9^{+4.6}_{-4.0}$ ($14.2^{+4.9}_{-4.3}$) companions, consistent with each other and with the main sample. We also quantify the degree of similarity between the two samples by calculating the integrated frequency for companions with masses between 1–13 M_Jup and orbital semi-major axes between 1–20 AU. We find frequencies of $46^{+12}_{-10}\%$ for the "misaligned" sample and $59^{+21}_{-15}\%$ for the "control" sample, which agree at the 1σ level. Thus, there is no evidence to suggest these two sub-samples are drawn from different populations.

Zoom In Zoom Out Reset image size

We also plotted histograms of the distribution of stellar masses and metallicities for systems with and without long-period companions (see Figure 13). We tested the significance of potential correlations between companion occurrence and either stellar mass or metallicity using the Kolmogorov–Smirnov statistic, which evaluates the probability that two samples are drawn from the same parent distribution. We find that in both cases this probability is high (21% for the two mass distributions and 67% for the two metallicity distributions), indicating that there are no statistically significant correlations between companion frequency and either of these two parameters. We note that this calculation does not take into account the variations in our sensitivity to companions for individual systems, but it is unlikely that a more complete analysis would change our conclusions on this point.

Zoom In Zoom Out Reset image size

We next consider how the multiplicity rate C for our sample compares to those of planetary systems detected by the Kepler transit survey and by radial velocity surveys, where we define multiplicity as the fraction of planetary systems containing more than one planet. Batalha et al. (2013) report the detection of 369 systems with multiple transiting planet candidates out of a total of 1797 stars with at least one transiting planet candidate. This corresponds to a multiplicity rate of approximately 21% for the Kepler sample. Tremaine & Dong (2012) perform a more detailed statistical analysis of the Kepler planet candidate sample and find a multiplicity rate ranging between 20%–50% depending on the distribution of mutual inclinations assumed in the calculation. Although both of these numbers are broadly consistent with our multiplicity rate, we note that the characteristics of the Kepler candidate multi-planet systems are dramatically different than in our sample. Our systems consist of a short-period, gas giant planet with another massive companion on a very long period (typically several years or more) orbit. In contrast, the Kepler candidate multi-planet systems typically consist of tightly packed sets of low-mass (smaller than Neptune) planets on orbits less than 100 days (e.g., Latham et al. 2011; Fabrycky et al. 2012; Steffen et al. 2013). Steffen et al. (2012) used the Kepler data set to demonstrate that hot Jupiter candidates are notably lacking in nearby, low-mass companions. For the majority of the Kepler systems, the frequency of massive long-period companions is unknown; these planets are unlikely to transit and radial velocity follow-up is challenging for most Kepler targets.

The sample of radial velocity planets provides a better basis for comparison, as it includes many systems with long-term radial velocity monitoring capable of detecting massive companions on long-period orbits (e.g., Fischer et al. 2001; Wright et al. 2007, 2009). Tremaine & Dong (2012) find a total of 162 single-planet systems and 33 multi-planet systems detected orbiting FGK dwarf stars as of 2010 August. This corresponds to a multiplicity fraction of 17% ± 3%, but this number only includes planets with fully resolved orbits. Wright et al. (2009) carry out a similar analysis including radial velocity accelerations, where they exclude accelerations in systems where there is a known stellar companion. They find that 14% of the 205 planetary systems known at the time have multiple confirmed planets, with another 14% showing evidence for an additional distant companion. They also note that this number is most likely an underestimate, as the occurrence rates for planets increase sharply toward smaller masses and radii (e.g., Howard et al. 2010; Petigura et al. 2013) where their survey has a high level of incompleteness. Our survey has a similar lack of sensitivity to very small, distant planets, and we find evidence for at least one companion around 27% of our target stars. Even after accounting for planets we might have missed, our multiplicity rate of 51% ± 10% for large companions to the transiting planets in our sample is still in reasonably good agreement with their result. Wright et al. (2009) also note that the observed excess of giant planets at short orbital periods (the "three-day pile-up") disappears for multi-planet systems, indicating that nearby gas giant companions are rare in systems with hot Jupiters. Our results support this finding, as there is no evidence for any gas giant companions interior to 1 AU for the systems in our survey.