Three's Company: An Additional Non-transiting Super-Earth in the Bright HD 3167 System, and Masses for All Three Planets

The NASA K2 mission uses the Kepler spacecraft to observe a series of fields, called campaigns, around the ecliptic plane. Near-continuous, high-precision photometry is obtained on 10,000–20,000 targets per campaign, most targets having 30 minute integrations. Campaign 8 (C8) was observed for 80 days from 2016 January 04 to 2016 March 23. The calibrated pixels were downloaded from the Mikulski Archive for Space Telescopes (MAST) and processed in the same fashion as Crossfield et al. (2015). In brief, following the methods of Vanderburg & Johnson (2014) and Vanderburg (2014), the photometry is divided into six roughly equal segments, and each is decorrelated against the location of the photocenter of the light using a 1D Gaussian process. The major systematic in the photocenter location is the roll of the spacecraft around the telescope foresight, which is corrected approximately every six hours. By switching antennae at the start of Campaign 8, the magnitude of the roll was reduced significantly from that seen in Campaign 7, resulting in overall higher quality light curves with higher precision.³⁵ One of the targets observed in Campaign 8 was HD 3167, a bright (V = 8.9, K = 7.0), nearby (46 pc), K0 dwarf star, also designated as EPIC 220383386. The detrended photometry is shown in the top panel of Figure 1.

Zoom In Zoom Out Reset image size

Three transits of a long-period, relatively deep (∼1 mmag) planet candidate were first detected in a by-eye search of the brightest targets in C8, marked in cyan in the top panel of Figure 1. On closer inspection, shallower transits at a significantly shorter period were also detected, marked in red. Using the TERRA algorithm (Petigura et al. 2013), two signals were found with periods of 0.959609 days (shown in the bottom left panel of Figure 1) and 29.8479 days (shown in the bottom right panel), with transit depths of 294 ppm and 946 ppm, respectively. These signals were subsequently reported by Vanderburg et al. (2016b) as HD 3167 b and c respectively. After removal of those signals, no additional transiting signals were found with an SNR above $5\sigma$, corresponding to ∼0.8 ${R}_{\oplus }$.

In addition, we performed several tests of the photometry to rule out obvious false positive scenarios prior to acquiring expensive, high-precision RV measurements. These included an adaptation of the model-shift uniqueness test, originally designed for Kepler data and described in Section 3.2.3 of Coughlin et al. (2016). In brief, the test searches for other significant transit-like events in the light curve when phased to the period of the putative planet signal: false positives will often show multiple significant events across all phases due to the higher levels of correlated noise. Both planets b and c passed the model-shift uniqueness test. We also included an adaptation of the Locality Preserving Projections (LPP) test, described in Thompson et al. (2015) for Kepler data, which uses dimensionality reduction and k-nearest neighbors to measure how similar a putative signal is to a planetary transit signal. Both planets b and c also passed the LPP test.

Finally, we examined the photocenters of light during the transits of planet b: significant motion of the photocenter of light away from the location of the putative host star during transit is a powerful technique for detecting false positive events due to background eclipsing binaries. For the original Kepler mission, the spacecraft pointing stability was so high that this method could be used to identify false positives lying well within the same pixel as the target star. We adapt the difference imaging technique used in Kepler (Bryson et al. 2013) to K2. Due to the strong roll motion in K2, there is a large change in the light distribution between two consecutive cadences. even in the absence of a transit. Instead, for each in-transit cadence, we look for out-of-transit cadences at the same roll angle and separated by exactly one thruster firing event. The roll angle is measured by the same technique as Vanderburg & Johnson (2014). By requiring the out-of-transit cadences to be close in time, we minimize of the impact of motion perpendicular to the roll axis due to, e.g., differential velocity abberation. However, the K2 roll motion is not exactly repeatable, and not all in-transit cadences have out-of-transit cadences that meet our requirement both before and after the transit. For cadences that do, we fit the PRF model of Bryson et al. (2010) to the in-transit and difference image, and compute the shift in the photocenter. We average over all cadences for which a difference can be computed, and calculate the probability that the observed distribution of offsets is consistent with the hypothesis that the location of the transit is consistent with the location of the target star. For simplicity, the distribution of offsets is assumed to be Gaussian in both row and column. At ninth magnitude, HD 3167 is highly saturated, resulting in large scatter in the row direction due to the bleed of saturated pixels, and the distribution is highly non-Gaussian. Nevertheless, Figure 2 shows no strong evidence that the source of the transit for planet b is offset from the target. There are only three transits of HD 3167 c, and one of those gives a poor fit to the photocenter location, so we do not perform the photocenter analysis on this planet.

Zoom In Zoom Out Reset image size

To determine the stellar parameters of the host star, we obtained three spectra of HD 3167 using Keck-HIRES with an S/N of ∼260 at 6000 Å, without using the iodine cell as is typical of the precision RV observations; Figure 3 shows a segment of a spectrum in the region of the Mg b triplet. We derived the stellar properties using the spectral forward-modeling procedure and line list of Brewer et al. (2016). We first fit for ${T}_{\mathrm{eff}}$, $\mathrm{log}g$, [M/H], and Doppler broadening using a scaled solar abundance pattern except for the alpha elements calcium, silicon, and titanium. We then fixed the stellar parameters and solved for the abundances of 15 elements. Finally, we repeated the process using this new abundance pattern. The results from fitting the three different spectra were nearly identical for all parameters. We then apply the empirical corrections from Brewer et al. (2016) to obtain the final parameters, summarized in Table 1.

Zoom In Zoom Out Reset image size

Table 1.  HD 3167 Stellar Parameters

Parameter	Value	Units
R.A.	00:34:57.52	hh:mm:ss
Decl.	+04:22:53.3	dd:mm:ss
EPIC ID	EPIC 220383386	⋯
2MASS ID	2MASS J00345752+0422531	⋯
V	8.941 ± 0.015	mag
K	7.066 ± 0.020	mag
Spectral Type	K0 V	⋯
${T}_{\mathrm{eff}}$	5261 ± 60	K
log g	4.47 ± 0.05	log10 (cm s−2)
R${}_{\star }$	0.872 ± 0.057	${R}_{\odot }$
M${}_{\star }$	0.866 ± 0.033	${M}_{\odot }$
${\rho }_{\star }$ a	1.902 ± 0.092	g cm−3
${\rho }_{\star ,b}$ b	${1.40}_{-0.79}^{+0.52}$	g cm−3
${\rho }_{\star ,c}$ c	${1.39}_{-0.94}^{+0.65}$	g cm−3
Distance	45.8 ± 2.2d	pc
[Fe/H]	0.04 ± 0.05	⋯
v sin i	1.7 ± 1.1	km s−1
log R'HK	−5.04	⋯

Notes.

^aSpectroscopically derived. ^bDerived from transit light-curve fit to planet b. ^cDerived from transit light-curve fit to planet c. ^dvan Leeuwen (2007).

Download table as:  ASCIITypeset image

The analysis procedure has been shown to recover gravities consistent with those of asteroseismology with an rms scatter of 0.05 dex (Brewer et al. 2015) and we adopt this as the uncertainty in $\mathrm{log}g$. Brewer et al. (2016) shows that there is a 39 K offset with temperatures derived from well-measured angular diameters. We add this in quadrature to their 25 K statistical uncertainties for a total uncertainty of 46 K. The statistical uncertainty in the [Fe/H] measurement is only 0.01 dex but the empirical correction at this temperature is 0.09 dex. We adopt half of the offset, 0.05 dex, as our uncertainty in [Fe/H]. Finally, we compare the results of the analysis to those given by the Stellar Parameter Classification tool (SPC; Buchhave et al. 2012, 2014) and SpecMatch (Petigura et al. 2015) and find they agree to within 1σ. Following the procedure in Crossfield et al. (2016), we use the free and open source isochrones Python package (Morton 2015) and the Dartmouth stellar evolution models (Dotter et al. 2008) to estimate the stellar radius and mass given in Table 1. The resulting stellar density is consistent with values derived in the transit analyses. HD 3167 was not included in Gaia Data Release 1 (Gaia Collaboration et al. 2016), possibly due to the incompleteness at the bright end or the poorer coverage along the ecliptic, where the K2 mission observes by necessity. However, future Gaia releases should produce a precise distance and allow for stronger constraints on the stellar parameters. From the HIPPARCOS parallax (Perryman et al. 1997), and following the same procedure as Brewer et al. (2016), we derive an age for HD 3167 of 7.8 ± 4.3 Gyr.

The K2 data show some longer-term variability that may be caused by stellar rotation (see Figure 1 of Vanderburg et al. 2016b). Examining the auto-correlation function of the light curve reveals a broad peak from 20 to 35 days, with a maximum at 27.2 days. The rotational velocity of 1.7 ± 1.1 km s⁻¹ is fairly poorly constrained, and allows a range of rotational periods from 10 to 40 days. These values are broadly consistent with the expected value for a field K-dwarf (see, e.g., Newton et al. 2016). We examine the correlations between stellar activity indicators and the measured radial velocities in Section 3.2.1.

2.2.1. Proper Motion

The proper motion of HD 3167 is quite large (107 mas yr⁻¹ in R.A. and −173 mas yr⁻¹ in decl.; Huber et al. 2016). In the 63 years since, the 1953 Palomar Observatory Sky Survey (POSS) images, HD 3167 has moved more than 125, enabling us to utilize archival POSS data to search for background stars that are now, in 2016, hidden by HD 3167. Using the 1953 POSS data, shown in the top panel of Figure 4, we find no evidence of a background star at the current postion of HD 3167 to a differential magnitude of ∼5 mag, shown in the bottom panel of Figure 4. Because HD 3167 is saturated in the POSS images, this sensitivity was estimated by placing fake sources at the epoch 2016 position of HD 3167 in the epoch 1953 image and estimating the 5σ threshold for detection. The photometric scale of the image (and hence, the magnitudes of the injected test stars) was set using the star located 1' to the southeast of HD 3167, which has an optical magnitude of approximately B = 15.5. This analysis does not rule out the most extreme background eclipsing binaries (a 50% eclipsing binary would produce a 1 mmag transit at a differential magnitude of 6.8 mag), but was sufficient for us to instigate the high-precision RV campaign.

Zoom In Zoom Out Reset image size

2.2.2. Adaptive Optics

We obtained near-infrared adaptive optics images of HD 3167 at Keck Observatory on the night of 2016 July 14 UT. Observations were obtained with the 1024 × 1024 NIRC2 array and the natural guide star system; the target star was bright enough to be used as the guide star. The data were acquired using the narrow-band Br-γ filter using the narrow camera field of view with a pixel scale of 9.942 mas/pixel. The Br-γ filter has a narrower bandwidth (2.13–2.18 $\mu {\rm{m}}$), but a similar central wavelength (2.15 $\mu {\rm{m}}$) compared the Ks filter (1.95–2.34 $\mu {\rm{m}};$ 2.15 $\mu {\rm{m}}$) and allows us to observe HD 3167 without saturation. A three-point dither pattern was utilized to avoid the noisier lower left quadrant of the NIRC2 array. The three-point dither pattern was observed with 10 coadds and a 0.726 s integration time per coadd for a total on-source exposure time of 65 s.

HD 3167 was measured with a resolution of 0050 (FWHM). No other stars were detected within 4'' of HD 3167. In the Br-γ filter, the data are sensitive to stars that have K-band contrast of ΔK = 3.4 mag at a separation of 01 and ΔK = 8.0 mag at 05 from the central star. We estimate the sensitivities by injecting fake sources with a signal-to-noise ratio of 5 into the final combined images at distances of N × FWHM from the central source, where N is an integer. The 5σ sensitivities, as a function of radius from the star, are shown in Figure 5. Beyond 4'', there are no additional stars visible in 2MASS out to a radius of ∼20''.

Zoom In Zoom Out Reset image size

2.2.3. RV Measurements

After the identification and validation of the two transiting planet signals, a high-cadence observing campaign was rapidly launched in order to obtain mass measurements while C8 was still visible. The final data set includes observations obtained with Keck/HIRES, Automated Planet Finder (APF)/Levy, and HARPS-N, described below. The full set of RV measurements is given in Table 2.

Table 2.  Radial Velocities

HJD${}_{\mathrm{UTC}}$	RVa	Unc.b	Inst.
(−2440000)	(m s−1)	(m s−1)
17576.70031	19526.75	0.75	HARPS-N
17576.96772	−8.69	0.93	APF
17578.69042	19521.29	1.51	HARPS-N
17578.95901	−12.99	1.05	APF
17579.67410	19519.00	1.58	HARPS-N
17579.94894	−5.76	1.33	APF
17580.11952	−5.44	0.56	HIRES
17580.71980	19520.10	1.20	HARPS-N
17581.10524	−0.98	0.52	HIRES

Notes.

^aZero-point offsets between instruments have not been removed and must be fit as free parameters when analyzing this data set. ^bStellar jitter has not been incorporated into the uncertainties.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Our observational setup for both Keck/HIRES and the APF/Levy was essentially identical to those described in Fulton et al. (2016) and Burt et al. (2014). We collected a total of 60 RV measurements using Keck/HIRES (Vogt et al. 1994), and 116 measurements using the Levy Spectrograph on the APF (Radovan et al. 2014; Vogt et al. 2014) at Lick Observatory between 2016 July 7 and 2016 December 2. For all of the Keck/HIRES measurements, we collected three consecutive exposures in order to mitigate the affects of stellar oscillations (Dumusque et al. 2011). The three measurements were then binned together before a jitter term, which includes a contribution from stellar jitter, is added in quadrature during the modeling process (see Section 3); this technique was not necessary at APF due to the smaller telescope aperture and longer exposure times. Whenever possible, we observed HD 3167 two times during a single night with maximum temporal separation to improve phase coverage for HD 3167 b.

Each Doppler spectrum was taken through a cell of gaseous iodine that imprints a dense forest of molecular absorption lines onto the stellar spectrum and serves as both a wavelength and point-spread function (PSF) reference. The slits chosen provided spectral resolving power of $R\sim {\rm{70,000}}$ and $R\sim {\rm{100,000}}$ for Keck and APF respectively. A series of iodine-free spectra were also collected using a narrower slit on both instruments ($R\sim {\rm{85,000}}/{\rm{120,000}}$ for Keck/APF). These spectra were deconvolved with the instrumental PSF and used as models of the intrinsic stellar spectrum. We modeled each RV observation as the deconvolved intrinsic stellar spectrum shifted by a best-fit RV and multiplied by an ultra-high-resolution iodine transmission spectrum. This is then convolved with an instrumental PSF, which is modeled as the sum of 13/15 Gaussians for Keck/APF (Butler et al. 1996). We reject measurements with SNR < 45, mid-exposure times before or after 13° twilight, and measurements collected when the star was within 20° of the moon. The rejected observations are not included in Table 2.

We also observed HD 3167 with the HARPS-N spectrograph (Cosentino et al. 2012) located at the 3.58 m Telescopio Nazionale Galileo on the island of La Palma, Spain. HARPS-N is a stabilized spectrograph designed for precise RV measurements. We observed HD 3167 76 times between 2016 July 7 (independently beginning the same night as the HIRES/APF campaign) and 2016 December 7, obtaining high-resolution optical spectra with a spectral resolving power of R = 115,000. Most of our observations consisted of 15 minute integrations, which yielded formal photon-limited Doppler uncertainties between 0.6 and 1.6 m s⁻¹. Similarly to the Keck/HIRES measurements, we typically observed HD 3167 two times per night, separated by a couple of hours, in order to better sample the inner planet's orbit; on several occasions, we observed HD 3167 up to six times per night. We measured radial velocities by calculating a weighted cross-correlation function between the observed spectra and a binary mask (Baranne et al. 1996; Pepe et al. 2002).

One factor that extended our RV campaign was the aliasing between the ∼1 day orbital period of planet b and the ∼1 month orbital period of outer planet c. Given the restrictions on observing enforced by the diurnal cycle and the tendency for telescope time to be allocated approximately monthly around the full moon, it was difficult at any single longitudinal site to secure the required phase coverage to break the degeneracy between planets b and c. Figure 6 shows, for each of the three telescopes, the b and c phase combinations of the observations. The HARPS-N observations, shown as yellow diamonds, represent the most precise measurements in our RV sample, but have large bands of phase combinations that are unsampled. Similarly, the HIRES measurements, shown as black open circles, do not cover the full range of phase combinations. Early analyses of the radial velocities from either of these sites individually led to degeneracies in the RV semi-amplitudes, and therefore masses, of the b and c planets. The APF observations, shown as green points, for which there is the most regular access to the telescope, provide comprehensive coverage of the phase combinations of planets b and c. By combining the higher precision but limited phase coverage observations from HARPS-N and HIRES with the lower precision but broad phase coverage of APF, we break the degeneracies and constrain the orbital solution as discussed below.

Zoom In Zoom Out Reset image size

We analyzed the transit signals for planets b and c independently in our light curve, using the same modeling, fitting, and MCMC procedures as described in Crossfield et al. (2016). As in that analysis, eccentricity was held to zero; for the RV analysis described in Section 3.2, we allowed the eccentricity of planet c to float. The results are shown in Table 5 and are consistent with the parameters given by Vanderburg et al. (2016b) for planets b and c. We examine the transit times of planet b and find no evidence of variations above the level of ∼15 minutes, shown in Figure 7. Occasional outliers are present in the individually derived transit times, but we conclude that these are likely a result of the low cadence of the Kepler observations combined with a non-perfect detrending. We exclude cadences affected by spacecraft thruster firings prior to analysis. In addition, we apply the cosmic-ray detection algorithm for K2 photometry developed by Benneke et al. (2016), but do not identify any cosmic-ray events as the source for the outliers in the transit timing.

Zoom In Zoom Out Reset image size

We analyzed the RV time series using the publicly available RV fitting package RadVel(B. J. Fulton & E. A. Petigura 2017, in preparation).³⁶ RadVel is written in object-oriented Python and is designed to be highly extensible, flexible, and documented for easy adaptation to a variety of maximum-likelihood fitting and MCMC applications. The standard version of RadVel downloadable from GitHub³⁷ includes a pipeline that is capable of modeling multi-planet, multi-instrument RV time series utilizing a fast Keplerian equation solver written in C.

Our likelihood function for this analysis follows that of Sinukoff et al. (2016):

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }=-\displaystyle \sum _{i}\left[\displaystyle \frac{{({v}_{i}-{v}_{m}({t}_{i}))}^{2}}{2({\sigma }_{i}^{2}+{\sigma }_{{\rm{j}}}^{2})}+\mathrm{ln}\sqrt{2\pi ({\sigma }_{i}^{2}+{\sigma }_{{\rm{j}}}^{2})}\right],\end{eqnarray} \tag{ 1 }$$

where v_i are the gamma-subtracted velocity measurements (${v}_{i}={v}_{i,\mathrm{inst}}-{\gamma }_{\mathrm{inst}}$, where ${\gamma }_{\mathrm{inst}}$ is an instrument-dependent term) with associated uncertainties ${\sigma }_{i}$, and ${v}_{m}k({t}_{i})$ is the Keplerian model at time t_i.

We first find the maximum-likelihood model using the Powell minimization technique (Powell 1964) then perturb the best-fit parameters by 1 part in 10⁵ to start 50 parallel MCMC chains. RadVel incorporates the affine-invariant sampler of the emcee package (Foreman-Mackey et al. 2013). The Gelman-Rubin (Gelman et al. 2003) and T_z statistics (Ford 2006) are checked in real-time during the MCMC exploration. The chains are deemed well-mixed and the MCMC is halted when the Gelman-Rubin statistic is within 3% of unity and ${T}_{z}\gt 1000$ for all free parameters. We chose to parameterize the Keplerian orbits using $\sqrt{e}\sin \omega$ and $\sqrt{e}\cos \omega$ instead of e and ω in order to increase convergence speed. We assigned uniform priors to $\sqrt{e}\sin \omega$, $\sqrt{e}\cos \omega$, velocity semi-amplitudes (K), and the zero-point offsets (γ). The jitter terms for each instrument (${\sigma }_{{\rm{j}}}$) are defined in Equation (2) of Fulton et al. (2015), and serve to capture the stellar jitter and instrument systematics such that the reduced ${\chi }^{2}$ of the best-fit model is close to 1. The ${\chi }^{2}$ values in Table 3 are reported without including the jitter terms, since including them would artificially reduce the final ${\chi }^{2}$ values. Gaussian priors were assigned to the ephemerides of the two transiting planets using the values reported in Vanderburg et al. (2016b). We examine the fits for system architectures from zero to three planets and choose the three-planet solution favored by the Bayesian information criterion (see Table 3 for details). The median values and the 68% credible intervals of the three-planet solution are reported in Table 4. The best-fit three-planet Keplerian model is shown in Figure 8, and the correlations between the derived planet masses and densities is shown in Figure 9. As noted in Figure 6, the incomplete coverage of the phase combinations of planets b and c leads to a slight degeneracy in the derived masses.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3.  Model Comparison

Statistic	0 planets	1 planet	2 planets	3 planets (adopted)
${N}_{\mathrm{data}}$ (number of measurements)	252	252	252	252
${N}_{\mathrm{free}}$ (number of free parameters)	6	9	14	19
rms (rms of residuals in m s−1)	4.71	4.22	3.52	3.16
${\chi }^{2}$ (assuming no jitter)	770.54	573.94	450.02	293.6
${\chi }_{\nu }^{2}$ (assuming no jitter)	3.13	2.36	1.89	1.26
$\mathrm{ln}{ \mathcal L }$ (natural log of the likelihood)	−736.59	−701.96	−662.21	−621.8
BIC (Bayesian information criterion)	1484.71	1418.45	1343.95	1268.13

Download table as:  ASCIITypeset image

Table 4.  The MCMC Posterior Values for the Three-planet Solution

Parameter	Value	Units
Orbital Parameters
Pb	0.959641 $\pm \ 1.1e-05$	days
$T{\mathrm{conj}}_{b}$	2457394.37454 ± 0.00044	JD
eb	≡0.0	⋯
${\omega }_{b}$	≡0.0	radians
Kb	${3.58}_{-0.25}^{+0.26}$	m s−1
Pc	29.8454 ± 0.0012	days
$T{\mathrm{conj}}_{c}$	${2457394.9787}_{-0.0011}^{+0.0012}$	JD
ec	<0.267	⋯
${\omega }_{c}$	$-{3.2}_{-1.9}^{+2.0}$	radians
Kc	2.24 ± 0.28	m s−1
Pd	${8.492}_{-0.024}^{+0.023}$	days
$T{\mathrm{conj}}_{d}$	2457806.1 ± 0.5	JD
ed	<0.36	⋯
${\omega }_{d}$	−3.2 ± 1.4	radians
Kd	2.39 ± 0.24	m s−1
Modified MCMC Step Parameters
$\sqrt{e}\cos {\omega }_{b}$	≡0.0	⋯
$\sqrt{e}\sin {\omega }_{b}$	≡0.0	⋯
$\sqrt{e}\cos {\omega }_{c}$	0.001 ± 0.15	⋯
$\sqrt{e}\sin {\omega }_{c}$	0.01 ± 0.24	⋯
$\sqrt{e}\cos {\omega }_{d}$	$-{0.14}_{-0.19}^{+0.23}$	⋯
$\sqrt{e}\sin {\omega }_{d}$	0.002 ± 0.23	⋯
Other Parameters
${\gamma }_{\mathrm{HIRES}}$	$-{0.9}_{-0.47}^{+0.46}$	m s−1
${\gamma }_{\mathrm{APF}}$	$-{0.51}_{-0.37}^{+0.36}$	m s−1
${\gamma }_{\mathrm{HARPSN}}$	19528.8 ± 0.23	m s−1
$\dot{\gamma }$	≡0.0	m s−1 day−1
$\ddot{\gamma }$	≡0.0	m s−1 day−2
${\sigma }_{\mathrm{HIRES}}$	${3.42}_{-0.35}^{+0.4}$	${\rm{m}}\ {{\rm{s}}}^{-1}$
${\sigma }_{\mathrm{APF}}$	${3.45}_{-0.27}^{+0.3}$	${\rm{m}}\ {{\rm{s}}}^{-1}$
${\sigma }_{\mathrm{HARPSN}}$	${1.4}_{-0.19}^{+0.22}$	${\rm{m}}\ {{\rm{s}}}^{-1}$

Note. The measured system velocity (γ) and the derived jitter term (${\sigma }_{\mathrm{jit}}$) are quoted for each of the three instruments. The reference epoch for γ, $\dot{\gamma }$, $\ddot{\gamma }$ is 2457652.6. 507,500 links were saved.

Download table as:  ASCIITypeset image

Table 5.  HD 3167 Planet Parameters

Parameter	Value	Units
(1)	(2)	(3)
Planet b
Period	${0.959641}_{-0.000012}^{+0.000011}$	days
Transit mid-point	2457394.37454 ± 0.00043	BJDTDB
${R}_{p}/{R}_{\star }$	${0.01744}_{-0.00089}^{+0.00170}$	⋯
$a/{R}_{\star }$	${4.082}_{-0.986}^{+0.464}$	⋯
b	${0.47}_{-0.32}^{+0.31}$	⋯
i	${83.4}_{-7.7}^{+4.6}$	deg
e	0 (fixed)	⋯
Transit depth	294	ppm
t14	${1.622}_{-0.074}^{+0.060}$	hr
Rp	${1.70}_{-0.15}^{+0.18}$	R⊕
K	${3.58}_{-0.26}^{+0.25}$	m s−1
Mp	5.02 ± 0.38	M⊕
ρ	${5.60}_{-1.43}^{+2.15}$	g cm−3
a	0.01815 ± 0.00023	au
${S}_{\mathrm{inc}}$	${1625}_{-222}^{+244}$	S⊕
Planet c
Period	29.8454 ± 0.0012	days
Transit mid-point	2457394.9788 ± 0.0012	BJDTDB
${R}_{p}/{R}_{\star }$	${0.0313}_{-0.0018}^{+0.0045}$	⋯
$a/{R}_{\star }$	${40.323}_{-12.622}^{+5.549}$	⋯
b	${0.50}_{-0.33}^{+0.31}$	⋯
i	${89.3}_{-0.96}^{+0.5}$	deg
e	<0.267	⋯
Transit depth	946	ppm
t14	${5.15}_{-0.19}^{+0.26}$	hr
Rp	${3.01}_{-0.28}^{+0.42}$	R⊕
K	${2.23}_{-0.28}^{+0.29}$	m s−1
Mp	${9.80}_{-1.24}^{+1.30}$	M⊕
ρ	${1.97}_{-0.59}^{+0.94}$	g cm−3
a	0.1795 ± 0.0023	au
${S}_{\mathrm{inc}}$	${16.6}_{-2.3}^{+2.5}$	S⊕
Planet d
Period	8.509 ± 0.045	days
$T{\mathrm{conj}}_{d}$	${2457806.07}_{-0.50}^{+0.52}$	BJDTDB
e	<0.36	⋯
Mp sin i	6.90 ± 0.71	M⊕
a	0.07757 ± 0.00027	au
${S}_{\mathrm{inc}}$	88.9 ± 6.2	S⊕

Note. t₁₄ is the total transit duration from the first to fourth contact. ${S}_{\mathrm{inc}}$ is the irradiation at the surface of the planet in units of the irradiation at Earth. For Planet d, $T{\mathrm{conj}}_{{\rm{d}}}$ is the time of inferior conjunction.

Download table as:  ASCIITypeset image

3.2.1. Search for a Third Planet

We search for additional planets in the RV data using the automated planet discovery pipeline described in Fulton et al. (2016) and Howard & Fulton (2016). In brief, this pipeline utilizes a custom implementation of the two-dimensional Keplerian Lomb–Scargle periodogram (2DKLS; O'Toole et al. 2009). Periodogram power is defined as a change in ${\chi }^{2}$ relative to a baseline ${\chi }^{2}$. For this particular search, the baseline ${\chi }^{2}$ is derived from the best two-planet model fit. The periodogram, shown in the top panel of Figure 10, demonstrates the change to the fit when adding a third planet as a function of the orbital period of that planet. Offsets between data from different instruments and inhomogeneous measurement uncertainties are incorporated into ${\chi }^{2}$. In order to assess the significance of peaks in the periodogram, we determine an empirical false alarm probability (eFAP) by fitting a log-linear function to the distribution of values in a given periodogram.

Zoom In Zoom Out Reset image size

We find a significant peak with eFAP ∼ 0.3% and a period of ∼8.5 days in the 2DKLS periodogram of the combined RV data set when searching for a third Keplerian signal. When we add this additional Keplerian into the MCMC fits described in Section 3.2, we see an improvement in the Bayesian Information Criterion (BIC, Liddle 2007) of 76, which indicates that the three-planet model is highly favored over the two-planet model.

We also calculate the 2DKLS periodogram for each instrument independently. In the 2DKLS periodogram for the APF data, we find that the highest periodogram value similarly falls at a period of ∼8.5 days, with an eFAP ∼ 20%. We find that the highest peak in the 2DKLS periodogram of the HARPS-N data falls at a period of ∼11 days, which is near an alias of 8.5 days caused by the sampling being concentrated around lunar cycles (1/8.5 days–1/29.5 days = 1/11.9 days). The second highest peak in the 2DKLS periodogram of the HARPS-N data falls at a period of 8.4 days. The HIRES data also shows an insignificant peak with a period of ∼11 days. The APF data, which has a much more uniform sampling due to the semi-dedicated nature of the telescope, is critical to break the monthly alias and reveal the true period of the third planet.

In order to examine whether the 8.5 day signal could be caused by a window function effect, in the fashion of α Cen Bb (Rajpaul et al. 2016), we perform the following test: using the real observing times, we generate a simulated RV curve from the properties of the two transiting planets. For each point, we generate an uncertainty drawn from a normal distribution of the quadrature sum of the observation error and the instrument jitter for the instrument that obtained that observation. We run the simulated RV curve through 2DKLS, and after removing the two known signals, we see no significant remaining power in the 7–10 day range, implying that the observed 8.5-day signal in the real data is not caused by a window function effect of the observations. We show the results of this final search in the lower panel of Figure 10.

We also examined whether the 8.5-day signal could be caused by stellar activity, since the period is potentially near an integer alias of the stellar rotation period. The best-fit jitter value for Keck/HIRES is surprisingly large in comparison to that from the HARPS-N data set. Long-term Keck/HIRES monitoring of stars with similar spectral types and activity levels show jitter as low as 1.8 m s⁻¹. Inspection of the residuals in Figure 8 show a systematic structure that appears to be present in only the Keck/HIRES data set. These correlated residuals are the source of the inflated jitter. We collected iodine-free template observations for this star on three different occasions and recalculated the velocity time series using each of the different templates. The results were comparable in each case and the structure in the residuals did not change significantly. We also searched for correlations of the velocity residuals with environmental and pipeline parameters. The Keck/HIRES velocity residuals are weakly correlated with both barycentric correction and S value. We tried subtracting a linear trend from RV against barycentric correction and/or S value by adding a term into the likelihood in the MCMC fit, but found only very modest improvement to the final jitter value and no significant difference to the final results. Since the structure in the residuals appears to be quasi-periodic and weakly correlated with the S value, we suspect that the source of the large jitter is likely caused by rotational modulation of starspots. The iodine technique used to extract the velocities from the Keck/HIRES and APF spectra could be more sensitive to the line-shape distortions produced by these starspots compared to the cross-correlation technique used to extract the velocities from the HARPS-N spectra. As shown in Figure 10, the signal of HD 3167 d is present in the HARPS-N and APF data, which do not show systematic structure in their residuals, so we are confident that the signal of planet d is not caused by stellar activity. We investigated this further by examining the stacked periodogram of the radial velocities (Mortier & Collier Cameron 2017) and noting that the strength of the 8.5-day signal peak in the periodogram increases with the addition of more data, as distinct from the behavior of a peak caused by quasi-periodic stellar activity.

The measured mass and radius of HD 3167 b (5.02 ± 0.38 M${}_{\oplus }$, ${1.70}_{-0.15}^{+0.18}$ R${}_{\oplus }$) indicate a bulk density of ${5.60}_{-1.43}^{+2.15}$ g cm⁻³; consistent with a predominantly rocky composition, but potentially having a thin envelope of H/He or other low-density volatiles. Figure 11 shows HD 3167 b in comparison with other small exoplanets with masses measured to better than 50% precision; the lines show the composition models of Zeng et al. (2016). We randomly draw 100,000 planet masses and radii from our posterior distributions, and compare them to the mass–radius relation of Fortney et al. (2007) for pure rock, finding results that are consistent with the models of Zeng et al. (2016). Assuming that the planet is a mixture of rock and iron, we compute the iron mass fraction from each random draw using Equation (8) of Fortney et al. (2007). We conclude that the iron mass fraction is smaller than 15% at 68% confidence and smaller than Earth's iron mass fraction (33%) at 85% confidence, under the assumption that the planet is a mixture of rock and iron, with no volatiles. The radius, ${1.70}_{-0.15}^{+0.18}$ R${}_{\oplus }$, brackets the putative transition radius from likely rocky to likely volatile rich at 1.6 R${}_{\oplus }$ proposed by Rogers (2015). Planetary envelopes in such close proximity to the host star are predicted to be stripped away, either through photo-evaporation (e.g., Owen & Jackson 2012; Lopez & Fortney 2014; Chen & Rogers 2016; Lopez 2016) or Roche lobe overflow (e.g., Valsecchi et al. 2014). Our constraints are consistent with the notion that ultra-short-period planets are predominantly rocky.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

HD 3167 c has a mass and radius of ${9.80}_{-1.24}^{+1.30}$ M${}_{\oplus }$ and ${3.01}_{-0.28}^{+0.42}$ R${}_{\oplus }$ respectively, also shown in Figure 11. The resulting bulk density of HD 3167 c is ${1.97}_{-0.59}^{+0.94}$ g cm⁻³. The mass and radius can be explained by a wide range of compositions, all of which include low-density volatiles such as water and H/He (Adams et al. 2008; Rogers & Seager 2010; Valencia et al. 2013). The planet evolution models of Lopez (2016) are consistent with an Earth-composition core surrounded by a H/He envelope comprising ∼2% of the total planet mass. Alternatively, the planet might be mostly water. With a K-band magnitude of 7, HD 3167 is amenable to transmission spectroscopy observations to detect the atmospheric constituents of planet c, discussed in Section 4, which will help to break compositional degeneracies. HD 3167 c receives an incident flux that is $\approx 16$ times that of Earth, and is much less susceptible to atmospheric photo-evaporation than planet b. Planet b could be a remnant core of a planet similar to planet c.

Unlike the case of Kepler-36, with an orbit tightly constrained by transit timing variations, the RV orbital fit of HD 3167 is not sufficiently precise to deduce the behavior of resonant harmonics directly. Thus, we approach this question from an alternative viewpoint—namely, we employ numerical experiments to examine whether the conditions required to establish such a resonant lock could have occurred in the system's evolutionary history. It is well known that MMRs arise from smooth convergent migration (in this case, likely due to interactions with the protoplanetary nebula), and the probability of capture depends both on the planetary eccentricities at the time of the resonant encounter, as well as the migration rate (Henrard 1982; Borderies & Goldreich 1984). Application of adiabatic theory (Neishtadt 1975) shows that resonance capture probability diminishes with increasing eccentricity and/or increasing migration rate (Batygin 2015). Accordingly, in our simulations, we circumvent the former issue by assuming that the planets approach one another on initially circular orbits, and only retain the migration rate as an adjustable parameter.

To facilitate orbital convergence and damping, we have augmented a standard gravitational N -body code with fictitious accelerations of the form (Papaloizou & Larwood 2000)

$$\begin{eqnarray}&&\displaystyle \frac{d{\boldsymbol{v}}}{{dt}}=-\displaystyle \frac{{\boldsymbol{v}}}{{\tau }_{\mathrm{mig}}}-\displaystyle \frac{2\,{\boldsymbol{r}}}{{\tau }_{\mathrm{dmp}}}\displaystyle \frac{({\boldsymbol{v}}\cdot {\boldsymbol{r}})}{({\boldsymbol{r}}\cdot {\boldsymbol{r}})},\end{eqnarray} \tag{ 2 }$$

where ${\tau }_{\mathrm{mig}}$ and ${\tau }_{\mathrm{dmp}}$ are the migration and damping timescales, respectively. For definitiveness, migration torque was only applied to the outer planet, while damping torques were exerted upon both planets. Additionally, the gravitational potential of the central star was modified to account for the leading-order effects of general relativity (Nobili & Roxburgh 1986). The simulations employed the Bulirsch–Stoer algorithm (Press et al.1992), and initialized the orbits in the plane, with random mean anomalies, $\sim 5 \%$ outside of the exact 7/2 resonance.

We have carried out a sequence of numerical experiments with ${\tau }_{\mathrm{mig}}$ ranging from the nominal type-I migration timescale of $\sim 5000$ years (Tanaka et al. 2002) to ${\tau }_{\mathrm{mig}}=3$ Myr (i.e., a typical protoplanetary disk lifetime; Armitage 2010), and with ${\tau }_{\mathrm{dmp}}=\infty$ as well as ${\tau }_{\mathrm{dmp}}={\tau }_{\mathrm{mig}}/100$ (Lee & Peale 2002). We tested each parameter combination with 10 cloned simulations, and did not observe capture into a 7/2 MMR a single time. As a consequence, we conclude that it is unlikely that the planets are presently affected by the nearby 7/2 resonance, and the orbital proximity to this commensurability is coincidental.

With the possibility of resonant interactions disfavored, we proceed with a purely secular (i.e., orbit-averaged) treatment of the dynamics. A specific question we now seek to address concerns the mutual inclinations within the system. In other words, what extent of misalignment among the angular momentum vectors of the planetary orbits is required for planet d to elude transit, while allowing planets b and c to transit simultaneously? Although an exact answer to this question can, in principle, be attained from numerical integrations, such calculations require a more precise knowledge of the input parameters (e.g., eccentricities, longitudes of periastron, etc.) than what is presently available. Consequently, here we settle for an approximate answer, which we deduce analytically from secular perturbation theory.

A conventional approach to modeling the long-term behavior of planetary systems that reside outside of mean-motion commensurabilities, is to replace the planetary orbits with massive wires and compute the resulting exchange of angular momentum (MD99). We note that, formally, this is equivalent to averaging the governing Hamiltonian over the mean longitudes (Morbidelli 2002). In the limit of low eccentricities and mutual inclinations (specifically, to second order in either quality), the inclination and eccentricity dynamics become decoupled, meaning that the uncertainties of the RV fit do not strongly affect the following calculations.

Within the context of this so-called Lagrange–Laplace secular theory (see Brouwer & Clemence 1961, for a complete discussion), the equations of motion for the complex inclination vector $z=i\,\exp ({\imath }\,{\rm{\Omega }})$, where i is the inclination and Ω is the ascending node, simplify to a linear eigenvalue problem:

$$\begin{eqnarray}&&\displaystyle \frac{{{dz}}_{j}}{{dt}}={\imath }\displaystyle \sum _{k=1}^{N}{B}_{{jk}}{z}_{k},\end{eqnarray} \tag{ 3 }$$

where the indexes run over the planets, and N = 3. The interaction coefficients B_jk depend exclusively on the planetary masses as well as the semimajor axis ratios, and comprise a matrix ${\boldsymbol{B}}$ that fully encapsulates the dynamics:

$$\begin{eqnarray}{B}_{{jj}} & = & -\displaystyle \frac{{n}_{j}}{4}\displaystyle \sum _{k=1,k\ne j}^{N}\displaystyle \frac{{m}_{k}}{{M}_{\star }}{\alpha }_{{jk}}{\bar{\alpha }}_{{jk}}{b}_{3/2}^{(1)}({\alpha }_{{jk}})\\ {B}_{{jk}} & = & \displaystyle \frac{{n}_{j}}{4}\displaystyle \frac{{m}_{k}}{{M}_{\star }}{\alpha }_{{jk}}{\bar{\alpha }}_{{jk}}{b}_{3/2}^{(1)}({\alpha }_{{jk}}).\end{eqnarray} \tag{ 4 }$$

In the above expression, $n=\sqrt{{{GM}}_{\star }/{a}^{3}}$ is the mean orbital frequency, $\alpha \lt 1$ is the semimajor axis ratio, ${b}_{3/2}^{(1)}({\alpha }_{{jk}})$ is a Laplace coefficient of the first kind, and $\bar{\alpha }=\alpha$ if ${a}_{j}\lt {a}_{k};$ $\bar{\alpha }=1$ if ${a}_{k}\lt {a}_{j}$. With these specifications of the problem, the solution to Equation (3) can be expressed as a super-position of N linear modes:

$$\begin{eqnarray}&&{z}_{j}=\displaystyle \sum _{k=1}^{N}{\beta }_{{jk}}\exp ({\imath }{f}_{k}t+{\delta }_{k}),\end{eqnarray} \tag{ 5 }$$

where f_k and ${\beta }_{{jk}}$ denote the eigenvalues and eigenvectors of ${\boldsymbol{B}}$, respectively. The scaled amplitudes of the eigenvectors and the phases ${\delta }_{k}$ are determined entirely by the specific choice of initial conditions.

For definitiveness, here we initialize the transiting planets (b and c) in the plane (${i}_{b}={i}_{c}=0;$ ${{\rm{\Omega }}}_{b},{{\rm{\Omega }}}_{c}$ undefined), and choose our reference direction to coincide with the present-day ascending node of the inclined planet d (${{\rm{\Omega }}}_{d}=0$). Although adopting this initial condition does not lead to a general analysis of the systems possible dynamical evolution, this simplification is justified given the current observational constraints. Consequently, the only free parameter that enters our calculations is planet d's inclination. Moreover, owing to the analytic nature of our solution, the computational cost associated with any one realization of the dynamics is negligible.

To obtain an absolute lower-bound on planet d's present-day inclination, we note that given a favorable configuration of the line of nodes relative to the line of sight, any inclination greater than ${i}_{d}\gt \arctan ({R}_{\star }/{a}_{d})=1\buildrel{\circ}\over{.} 3$ will allow planet d to elude transit for some fraction of the time, potentially during the 80 day duration of the K2 observations. The greater the mutual inclination, the larger the fraction of time that planet d does not transit, rising from ∼7% for an inclination of 3°, to ∼80% for inclinations of 10°. The nodal configuration assumption therefore becomes progressively less stringent as the adopted value of i_d increases, and it is of interest to estimate the critical i_d beyond which this limitation can be alleviated altogether.³⁸ Moreover, such a calculation can further inform a maximal i_d, beyond which none of the planets co-transit.

Following Spalding & Batygin (2016), we define a mutual inclination

$$\begin{eqnarray}&&{\eta }_{{jk}}=\sqrt{{z}_{j}{z}_{j}^{* }+{z}_{k}{z}_{k}^{* }-({z}_{j}{z}_{k}^{* }+{z}_{k}{z}_{j}^{* })},\end{eqnarray} \tag{ 6 }$$

and adopt the following criterion for a pair of planets to co-transit:

$$\begin{eqnarray}&&\sin ({\eta }_{{jk}})\lt \displaystyle \frac{{R}_{\star }}{{a}_{j}}+\displaystyle \frac{{R}_{\star }}{{a}_{k}}.\end{eqnarray} \tag{ 7 }$$

Generically, as the orbits exchange angular momentum, their mutual inclinations, ${\eta }_{{jk}}$, will experience oscillatory motion. An example of this behavior, taking ${i}_{d}=20^\circ$ as an initial condition, is shown in the left panel of Figure 13. For reference, the solid lines denote the analytic solutions obtained by matrix inversion, while the dotted lines show the numerical solution computed with the N-body code described above. Although a small discrepancy exists in the oscillation frequencies computed analytically and numerically, the amplitudes of oscillation (which are the more relevant quantities for the question at hand) are well captured by secular perturbation theory.

Zoom In Zoom Out Reset image size

In the particular case shown in the left panel of Figure 13, the orbital architecture of the observed system is correctly reproduced, without assumptions about the current lines of nodes. That is in this case, given almost any nodal configuration, a viewing geometry where planets b and c co-transit, will not permit planets b and d to co-transit also. To estimate the critical inclination of planet d below which all three planets co-transit, we have computed the maximal and minimal extents of mutual inclinations between planets b and c as well as b and d, as a function of i_d. These results are shown in the right panel of Figure 13. Cumulatively, our theoretical calculations suggest that, although planet d can escape transit for inclinations as small as ∼13, for inclinations above ∼15° the allowed range of nodal alignment that would result in planet d transiting becomes so vanishingly small that, in the absence of observed transits, we conclude that the mutual inclination that reproduces the observed orbital misalignment of the HD 3167 system is most likely greater than ∼15°.

While the flavor of secular theory employed above adequately captures the dynamics of the system over the inclination range shown in Figure 13, the Lagrange–Laplace model is well known to break down at sufficiently high inclinations. Specifically, within the context of the problem at hand, it is reasonable to expect that provided sufficiently large i_d, the system will enter the Kozai–Lidov (Kozai 1962; Lidov 1962) resonance, which can facilitate large-scale oscillations of the eccentricities. A typically quoted inclination, necessary for Kozai–Lidov oscillations to ensue, is 392. Consistently, here we find numerically that when planet d's inclination exceeds ${i}_{d}\gtrsim 41^\circ$, the system enters the Kozai–Lidov regime, and planet d's eccentricity begins to experience oscillations coupled with its argument of pericenter. The small discrepancy in the critical value of the inclination can almost certainly be attributed to the apsidal precession generated by general relativistic effects and the quadrupolar field of the inner planet b (Batygin et al. 2011), as well as the non-negligible mass of planet d itself (Naoz et al.2013).

Intriguingly, the commencement of Kozai–Lidov oscillations is not synonymous with the onset of dynamical instability. Instead, the system remains stable for at least $100$ Myr for inclinations up to ${i}_{d}\sim 60^\circ$ (an example of stable evolution with ${i}_{d}\sim 55^\circ$ is shown in Figure 14). It is only above an inclination of ${i}_{d}\sim 65^\circ$, that eccentricity oscillations become sufficiently extreme, for subsequent orbit crossing to ensue. In this regard, the dynamics of the system entails an observational consequence: if follow-up RV observations sharpen the estimate of planet d's eccentricity to a value that is close to zero, that would imply that planet d's inclination lies below ${i}_{d}\lt 40^\circ$. Conversely, significant orbital eccentricity in the system would point toward ${i}_{d}\simeq 41^\circ \mbox{--}66^\circ$ as the more likely range of orbital misalignment. Constraining the inclination to 15°–60°, under the relaxed assumption that requires no special configuration of the lines of nodes, implies a true mass of 7.1–13.8 ${M}_{\oplus }$ for HD 3167 d.

Zoom In Zoom Out Reset image size

The dynamical analysis presented herein shows that the observed orbital architecture of the HD 3167 system can be naturally explained if the orbital inclination of planet d exceeds $\sim 15^\circ$, without invoking the need for the system to be observed at a given configuration and time. An intriguing question, then, concerns the origins of such a highly misaligned orbital architecture. One distinct possibility is a transient dynamical instability, that would have led to chaotic excitation orbital inclinations. Although such a scenario is not strictly impossible, the consistency of our RV fit with circular orbits renders such an evolutionary sequence unlikely. Some additional circumstantial evidence for long-term stability is the lack of a dense, hot disk around HD 3167, like that orbiting the G8V/K0V star HD 69830 (Beichman et al. 2006), which also hosts three planets (Lovis et al. 2006). Examining the WISE photometry (Cutri et al. 2014), we find no evidence for an excess, which is expected for mature stars but the presence of which may be indicative of a recent disruptive event.

An alternative, and perhaps more plausible solution is that the orbits have inherited their inclinations from a primordially misaligned star. Over the past few years, theoretical evidence has been marshaled in support of the notion that stars can become misaligned with respect to their protoplanetary disks, during the T-Tauri stage of their lifetimes (Bate et al. 2010; Lai et al. 2011; Batygin 2012; Lai 2014; Spalding & Batygin 2014; Matsakos & Königl 2016). An attractive feature of the primordial misalignment theory is that it can simultaneously account for the observed distribution of spin–orbit misalignments of hot Jupiters (Spalding & Batygin 2015) as well as the inherent inclination dispersion (often referred to as the Kepler dichotomy Mazeh et al. 2015; Ballard & Johnson 2016) of sub-Jovian planets (Spalding & Batygin 2016). Viewed in this context, HD 3167 probably represents an evolutionary outcome of a close-in planetary system that formed in a relatively quiescent environment, and was perturbed out of orbital alignment through secular exchange of angular momentum between the planets and the young star, while retaining orbital stability.