A Low-mass Exoplanet Candidate Detected by K2 Transiting the Praesepe M Dwarf JS 183

A number of focused studies and all-sky surveys have determined cluster membership probabilities for JS 183 based on proper motion and/or RVs. Kraus & Hillenbrand (2007) included JS 183 in their focused proper-motion study of the Praesepe cluster, finding a membership probability of 99.4%. Boudreault et al. (2012) similarly combined UKIDSS proper motions and photometry to obtain a cluster membership probability of 92%.

JS 183 has also been included in all-sky proper-motion surveys, with the following proper motions (μ_α, μ_δ) [mas yr⁻¹] from USNO B1, PPMX, and URAT, respectively: (−34, −14), (−39, −17), (−34, −17). These agree well with the mean Praesepe cluster motion, which has been reported as (−29, −7), (−34, −7.5), and (−36, −13) from Adams et al. (2002), Boudreault et al. (2012), and Kraus & Hillenbrand (2007), respectively. In summary, there is strong consensus in the literature that JS 183 is a high-probability member of the Praesepe cluster based on its proper motion and photometry.

The accepted age of the Praesepe cluster is ≈600 Myr (e.g., Adams et al. 2002), although an older age of ≈800 Myr was recently advocated (Brandt & Huang 2015). Praesepe has been found to have a supersolar metallicity, with most values around [Fe/H] ≈ 0.1 (for a good summary see Boesgaard et al. 2013), but with some values as high as [Fe/H] = 0.27 ± 0.10 (e.g., Pace et al. 2008). In our analysis we make use of the metallicity to estimate the stellar radius expected from empirical M dwarf radius relations, and we use the cluster age to place the JS 183 planet system in the context of other systems of known age.

Barrado y Navascués et al. (1998) conducted a comparison of Praesepe and Hyades chromospheric/coronal activity. Their reported data for JS 183 (identified as HSHJ 163) from a 1997 December Keck/HIRES observation indicated a Balmer line equivalent width EW(Hα) = 0.00 Å and X-ray luminosity log L_X < 28.07 erg s⁻¹. Adams et al. (2002) used low-resolution spectroscopy to obtain a spectral type of M3.5e, indicating activity, with EW(Hα) = −0.3 Å. Kafka & Honeycutt (2006) also report an EW(Hα) = −0.16 Å. In summary, all of these studies show Hα just very weakly in emission or filled in, suggesting modest chromospheric activity.

The photometric observations from K2 require considerable attention to the construction of proper apertures and removal of systematic noise (see, e.g., Vanderburg & Johnson 2014; Foreman-Mackey et al. 2015; Aigrain et al. 2016). Here we describe the extraction of the light curve of JS 183, the removal of systematics, the identification of the transit signal, and the characterization of out-of-transit variability.

3.1.1. Light Curve Preparation

JS 183 was observed by K2 in its Campaign 5 field, from MJD 57,139.13 to 57,213.93, acquiring 3448 photometric observations at the standard 30-minute cadence. We downloaded the light curve file for JS 183 from the Mikulski Archive for Space Telescopes (MAST), which includes both Simple Aperture Photometry (SAP) fluxes and a version of the light curve after application of the Kepler Pre-search Data Conditioning (PDC) pipeline (Twicken et al. 2010; Smith et al. 2012; Stumpe et al. 2012). This algorithm effectively models trends that are common to the ensemble of light curves, which was sufficient for the nominal Kepler mission but does not account for the diversity of pointing-related systematics present in K2 (see Vanderburg & Johnson 2014; Van Cleve et al. 2016, for details). We therefore additionally used the K2SC Gaussian process (GP) based systematics correction algorithm of Aigrain et al. (2016) to detrend both SAP and PDC light curves with respect to pointing variations. This algorithm models the stellar variability simultaneously with the pointing-dependent systematics and outputs a model for each component, so that either or both can be subtracted at will.

While the PDC light curve is usually the version of choice to look for transit signals, the PDC pipeline sometimes removes true astrophysical signals on timescales >20 days. Therefore, it is usually considered advisable to use the SAP light curve as a starting point for any detailed analysis of the transits and out-of-transit variability. However, after a visual comparison of the SAP and PDC light curves (both before and after modeling with K2SC), we opted to work with the PDC light curve throughout this paper. This decision was made because the SAP light curve is affected by long-term trends that are probably systematic, and because the K2SC modeling of the (shorter-timescale) pointing-related systematics is less successful in the SAP than in the PDC light curve for this particular object. Figure 1 shows the K2 PDC light curve for JS 183 along with the variability and systematics models from K2SC, as well as the light curve after correcting systematics only, and after also removing the variability model. Figure 2 shows the individual transits in more detail, using the detrended version of the light curve shown as version (e) in Figure 1. We estimated the 6 hr Combined Differential Photometric Precision (CDPP) of the fully detrended light curve as 274 ppm (the CDPP is the standard measure of photometric precision on transit timescales for Kepler light curves; see Christiansen et al. 2012).

Zoom In Zoom Out Reset image size

3.1.2. Transit Detection

The transits of JS 183 were identified independently by two methods. The first was visual inspection of the SAP and PDC light curves for all likely members of Praesepe observed by K2 (the transits are clearly visible in Figure 1). The second was a systematic search for transits, which we performed for all the K2SC-detrended PDC light curves for Campaign 5 (after subtraction of both systematics and variability models), in which JS 183 was identified. The details of this search are given in Pope et al. (2016), who also report the full list of transit candidates. The initial search found transits with a period of ∼10.1 days, a depth of ∼7 parts per thousand, and a duration of ∼3 hr, consistent with a transiting companion somewhat smaller than Jupiter. A full analysis of the transits is presented in Section 4.1.

Zoom In Zoom Out Reset image size

3.1.3. Out-of-transit Variability and Rotation

In order to confirm the properties of the transit, we observed JS 183 on the night of UT 2016 February 26, using a thinned, back-side-illuminated Tektronix 2048 CCD ("new2K") on the USNO, Flagstaff Station 40-inch telescope. A standard Kron−Cousins I-band filter was used with an exposure time of 120 s (based on test exposures from the previous night) and readout overhead of 51 s, resulting in a net observing cadence of 171 s. A total of 160 frames of JS 183 were obtained between UT 02:17 and UT 09:51. The frames were processed in real time with bias and flat-field frames obtained at the beginning of the night. Sky conditions were clear, although the seeing gradually worsened from about 12 FWHM at the beginning of the night to about 22 at the end of the exposure series.

As a result of these observations, we obtain a clear detection of the transit, albeit at a lower photometric precision than from K2. The detection took place about 221 days after the last of the K2 observations, an offset of 22 orbital periods. The USNO light curve is shown in conjunction with the K2 light curve in Figure 3. In Section 4.1 we show that incorporating the USNO light curve into the analysis gives a consistent solution with modeling the K2 light curve alone.

Zoom In Zoom Out Reset image size

We observed JS 183 with the Keck/HIRES spectrograph (Vogt et al. 1994) at six epochs: one in 1997 as part of another program, two in 2015, and three in 2016 as follow-up to the K2 light curve (see Table 2). All observations achieved spectral resolution R > 48,000. The 1997 observation used a wavelength range of 6300–8725 Å, while the 2015 and 2016 observations used a range of 4800–9200 Å, except for the final one. The images were processed and spectra extracted and calibrated using the makee software¹⁴ written by Tom Barlow. The spectral type estimate for the star from these high-dispersion data is M3–M4. We infer the presence of a small amount of stellar chromospheric activity based on the fact that the Hα line is completely filled in, with no apparent absorption or emission, and thus no equivalent width measurement is possible. The Hβ line is weakly in emission with EW(Hβ) = −0.45 Å.

From these spectra we obtained single-lined RVs of the host star. Heliocentric velocities were measured in the four Keck/HIRES spectra by cross-correlating with RV standards obtained at the same time as the program spectra. The mean and standard deviations of the multiorder RVs are summarized in Table 2. The systemic RV inferred from our four measurements is 35.34 ± 0.17 km s⁻¹.

The mean RV of the Praesepe cluster based on high-resolution spectra of FGK members has been determined as 37.7 km s⁻¹ (Barrado y Navascués et al. 1998) and as 34.6 ± 0.7 km s⁻¹ (Mermilliod & Mayor 1999). The mean RV of JS 183 from the observations reported here is in between these values; thus, we conclude that JS 183 is dynamically consistent with Praesepe cluster membership.

In addition, Barrado y Navascués et al. (1998) obtained medium-resolution spectra of a number of M dwarfs in Praesepe, including JS 183 (using the designation HSHJ 163). For JS 183, Barrado y Navascués et al. (1998) measure an RV of 32.4 km s⁻¹, which was consistent with Praesepe membership given the low RV accuracy from that paper, which was reported to be ∼7 km s⁻¹. In Section 4.3 we use the six RV measurements reported here to set an upper limit on the mass of the transiting planet.

Since we do not have dynamical confirmation of this planet from RV detection of the signal, we have taken two separate approaches to the analysis of the photometric observations, presented in the next two subsections. The consistent results of the two approaches provide additional evidence for the reliability of the overall system properties.

4.1.1. Primary Analysis: GP-transit Model

The K2 light curve of JS 183 is shown in Figure 1 and displays evolving starspot modulation with an amplitude of ∼2% and a period of ∼24 days. To account for the effect of this stellar variability across each transit, we simultaneously model the out-of-transit variations at the same time as fitting for the transit parameters. We therefore model the full light curve as the sum of a GP plus transit model (hereafter GP-transit model). The transit component is based on the model presented in Irwin et al. (2011), which is a modified version of the (JKT)EBOP family of models, but which uses the analytic method of Mandel & Agol (2002) for the quadratic limb-darkening law. The GP component, which is used to describe the out-of-transit variability, is based on the george package (Ambikasaran et al. 2014; Foreman-Mackey et al. 2014). For the K2 light curve, the GP component is a quasi-periodic exponential-sine-squared kernel, which is essentially a periodic kernel that is allowed to evolve over time, i.e., mimicking evolving starspot modulation.

In addition to the K2 discovery light curve, we use the USNO light curve described in Section 3.2. The single-night observation is not long enough to clearly see the evolving starspot modulation displayed in the K2 light curve. Instead, the dominant variations display a relatively rough behavior over short (≲1 hr) timescales. Accordingly, we opt for a Matern-3/2 kernel in our GP model, as its covariance properties are more suited to the observed behavior.

We simultaneously modeled the K2 and USNO light curves using our GP-transit model, with the parameter space explored through the affine-invariant Markov chain Monte Carlo (MCMC) method, as implemented in emcee (Foreman-Mackey et al. 2013). This model will be presented in detail in E. Gillen et al. (2017, in preparation). We ran the MCMC for 50,000 steps with each of 144 "walkers" and derived parameter distributions from the last 25,000 steps of each chain (after thinning each chain from inspection of the autocorrelation lengths for each parameter). In addition, the K2 model was supersampled to 1-minute cadence.

Since we do not have strong prior information about the eccentricity of this orbit, we opted to apply four different models to this system. These models are presented in Table 3 and differ only in their eccentricity, effective temperature, and surface gravity priors, the latter two giving different priors on the stellar density: (A) an eccentric orbit with a stellar density prior from the empirical relations of Mann et al. (2016a) (see Section 4.2), (B) an eccentric orbit with a stellar density prior from the SED fitting (again see Section 4.2), (C) an eccentric orbit with an uninformative stellar density prior, and (D) a circular orbit with an uninformative stellar density prior. Model A is our main model. We then test this stellar density prior first by relaxing it slightly in model B and then by relaxing it fully with an uninformative prior in model C. Finally, in model D, we keep the uninformative stellar density prior but force the orbit to circular to investigate the effect of eccentricity.

Table 3. Priors on System Parameters for Analysis in Section 4.1.1

Case	Eccentricity	Stellar Density (g cm−3)	Effective Temperature (K)	$\mathrm{log}g$
A (eccentric, empirically constrained SDa)	${ \mathcal U }(0,1.0)$	${7.266}_{-1.431}^{+1.877}$	${ \mathcal N }(\mu \,=\,3350,\sigma \,=\,50)$	${ \mathcal N }(\mu \,=\,4.79,\sigma \,=\,0.08)$
B (eccentric, SED-constrained SD)	${ \mathcal U }(0,1.0)$	${3.703}_{-2.534}^{+8.035}$	${ \mathcal N }(\mu \,=\,3350,\sigma \,=\,50)$	${ \mathcal N }(\mu \,=\,4.5,\sigma \,=\,0.5)$
C (eccentric, unconstrained SD)	${ \mathcal U }(0,1.0)$	${ \mathcal U }(0,160)$	...	...
D (circular, unconstrained SD)	δ	${ \mathcal U }(0,160)$	...	...

Notes. ${ \mathcal U }(a,b)$—uniform density with a minimum a and maximum b. ${ \mathcal N }(\mu ,\sigma )$—normal density with mean μ and standard deviation σ. δ—Dirac's delta function.

^aSD = stellar density.

Download table as:  ASCIITypeset image

The effective temperature and surface gravity priors in models A and B naturally give rise to priors on the limb-darkening coefficients. These were computed using the LDTK toolkit (Parviainen 2015), which allows us to propagate the uncertainties in these stellar parameters (and the metallicity Z = 0.1 ± 0.1) through the PHOENIX stellar atmosphere models (Husser et al. 2013) into our coefficients. The uncertainties on the limb-darkening coefficients were then inflated by a factor of 50 to account for differences between model grids. In the light curve modeling, the limb-darkening parameterization follows the triangular sampling method of Kipping (2013; see that paper for the relationship between the "q" parameters plotted in Figure 6 and standard quadratic limb-darkening coefficients). Uninformative priors were used on all other transit parameters.

Using models A–D, we can explore the information content of the data, how our priors affect our results, and the confidence with which we can characterize this system. The individual models are presented and compared below.

4.1.1.1. Individual Models: Investigating the Effect of Eccentricity and Assumed Stellar Density

Given that the planet orbital period lies close to an integer multiple of the K2 cadence and that the USNO light curve is of modest precision, the eccentricity of the orbit is unlikely to be well constrained from these data alone. As we believe the empirically determined stellar density constraint to be well-founded, we opt to use model A (eccentric orbit with a stellar density constraint from empirical relations) as the main model discussed here, but we also present results from the other three for comparison and to investigate the effect of our assumptions.

Figure 3 shows the full K2 and USNO light curves along with the GP-transit model A. The model is an acceptable fit to the large-scale structure of both light curves. In the K2 case, the GP model predicts a stellar rotation period of 23.8 days and is able to reproduce the evolving shape of the modulation. For the USNO light curve, the GP favors a smooth model for the out-of-transit variations due to the observational uncertainties, even though apparent correlations can be seen in the residuals. The significance of these variations (as well as lower-level ones present in the K2 light curve) are investigated at the end of this section. The top row of Figure 4 shows the model A transit fits to both the K2 and USNO light curves (left and right, respectively). In the K2 panel, the effect of the orbital period being close to an integer multiple of the cadence can be seen. Table 4 presents the parameters of model A (leftmost results column), and Figure 5 presents the 2D contours and 1D histograms of the transit parameter MCMC chains from which the posteriors reported in Table 4 are derived. There are significant correlations between the parameters of interest, as expected given that we only have photometric data of the transit (i.e., no secondary eclipse or RVs). Allowing an eccentric orbit drives most of the observed correlations: fixing the orbit to circular removes the correlations of interest, barring those between the cosine of the inclination (cosi) and the radius sum ((R_p+R_*)/a). The correlations arise from the fact that small changes in a given parameter can be effectively masked by corresponding changes in one or more others without significant modification to the transit shape.

Zoom In Zoom Out Reset image size

Table 4. Parameters of the Different Light Curve Models Applied to JS 183

Parameter	Symbol	Unit	Value
			Model A	Model B	Model C	Model D
Transit Parameters
Sum of radii	(R*+Rp)/a		${0.0339}_{-0.0022}^{+0.0061}$	${0.050}_{-0.013}^{+0.020}$	${0.069}_{-0.020}^{+0.026}$	${0.0403}_{-0.0062}^{+0.0175}$
Radius ratio	${R}_{{\rm{p}}}/{R}_{* }$		${0.0786}_{-0.0026}^{+0.0086}$	${0.0852}_{-0.0084}^{+0.0097}$	${0.0887}_{-0.0113}^{+0.0097}$	${0.0794}_{-0.0040}^{+0.0084}$
Orbital inclination	i	°	${89.25}_{-0.61}^{+0.53}$	88.0 ± 1.6	86.4 ± 4.2	${88.82}_{-1.30}^{+0.82}$
Orbital period	P	days	10.134589 ± 0.000084	10.134590 ± 0.000098	10.134603 ± 0.000102	10.134588 ± 0.000094
Time of transit center	Tcenter	BJD	2,457,282.6279 ± 0.0011	2,457,282.6279 ± 0.0015	2,457,282.6281 ± 0.0020	2,457,282.6278 ± 0.0012
	$\sqrt{e}\cos \omega$		0.02 ± 0.43	0.01 ± 0.52	−0.03 ± 0.6
	$\sqrt{e}\sin \omega$		$-{0.29}_{-0.28}^{+0.20}$	−0.04 ± 0.42	${0.28}_{-0.50}^{+0.38}$	⋯
Eccentricity	e		${0.24}_{-0.18}^{+0.27}$	${0.27}_{-0.19}^{+0.27}$	${0.46}_{-0.30}^{+0.21}$	⋯
Longitude of periastron	ω	°	$-{76.8}_{-69.7}^{+53.0}$	−10.0 ± 130.0	58.0 ± 75.0	⋯
K2 linear LDa coefficient	uK2		0.456 ± 0.082	0.67 ± 0.26	0.88 ± 0.38	0.60 ± 0.27
K2 nonlinear LD coefficient	${u}_{K2}^{\prime }$		0.24 ± 0.20	0.09 ± 0.29	−0.11 ± 0.39	0.07 ± 0.34
USNO linear LD coefficient	uUSNO		0.304 ± 0.079	0.27 ± 0.22	${0.28}_{-0.20}^{+0.31}$	0.60 ± 0.27
USNO nonlinear LD coefficient	${u}_{\mathrm{USNO}}^{\prime }$		0.21 ± 0.22	${0.16}_{-0.22}^{+0.31}$	${0.13}_{-0.23}^{+0.32}$	0.07 ± 0.34
Out-of-transit Variability Parameters
K2 amplitude	AK2	%	${0.230}_{-0.042}^{+0.058}$	${0.226}_{-0.042}^{+0.056}$	${0.222}_{-0.041}^{+0.057}$	${0.226}_{-0.041}^{+0.054}$
K2 timescale of SqExp term	${l}_{\mathrm{SE}K2}$	%	156.2 ± 3.4	156.2 ± 3.3	155.6 ± 3.3	156.3 ± 3.3
K2 scale factor of ExpSine2 term	${{\rm{\Gamma }}}_{\mathrm{ESS}K2}$	days	0.484 ± 0.090	${0.491}_{-0.086}^{+0.092}$	0.494 ± 0.092	0.490 ± 0.088
K2 period of ExpSine2 term	${P}_{\mathrm{ESS}K2}$	days	23.83 ± 0.29	23.84 ± 0.28	23.85 ± 0.31	23.83 ± 0.27
K2 white-noise term	σK2	%	1.371 ± 0.017	1.371 ± 0.018	1.371 ± 0.017	1.371 ± 0.017
USNO amplitude	AUSNO	%	0.36 ± 0.18	0.32 ± 0.18	0.31 ± 0.20	0.37 ± 0.17
USNO timescale	lUSNO	days	${5.0}_{-2.0}^{+2.6}$	5.9 ± 2.7	6.4 ± 2.6	5.0 ± 2.4
USNO white-noise term	σUSNO	%	1.167 ± 0.071	${1.175}_{-0.069}^{+0.079}$	${1.175}_{-0.068}^{+0.075}$	${1.173}_{-0.065}^{+0.071}$

Note.

^aLD = limb darkening.

Download table as:  ASCIITypeset image

We now compare the four different models to investigate the effect of our assumptions about the eccentricity and stellar density. The four rows of Figure 4 show (in descending order) the fits of models A, B, C, and D. All models are able to reproduce the observed transit shapes within the observational uncertainties, even though they have different shapes and durations. These transit models are integrated to the K2 cadence. To highlight the effect of this observational cadence on the depicted transit models, Figure 6 shows the unintegrated transit of model A (red line and shaded regions) and the integrated model predictions for the times of observation (blue points).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 7 presents the posterior parameter distributions of the four GP-transit models. As expected, the period and ephemeris are essentially insensitive to the choice of model and underlying assumptions. As previously mentioned, the limb-darkening coefficients of models A and B are theoretically constrained given the assumed stellar temperature, surface gravity, and cluster metallicity, whereas they are unconstrained in models C and D; the differences can be seen in the limb-darkening panels of Figure 7. In the distributions for the radius sum ((R_p+R_*)/a), radius ratio (R_p/R_*), inclination (i), and the combination terms $\sqrt{e}\cos \omega$ and $\sqrt{e}\sin \omega$, we see the effect of our assumptions about the eccentricity and stellar density. First, the stellar density: this essentially places a prior on the transit duration, which in turn constrains the radius ratio, radius sum, and $\sqrt{e}\sin \omega$ (and to lesser extents $\sqrt{e}\cos \omega$ and the orbital period P). Most importantly, model A (black), which places a tight constraint on the stellar density, constrains the transit duration such that the model favors a single-peaked radius ratio distribution around R_p/R_* ∼ 0.075–0.08. Relaxing the constraint on the stellar density (models B and C; orange and blue, respectively) allows the model to explore an additional family of solutions composing a larger planet-to-star radius ratio (and correspondingly these models explore a larger range of inclinations, radius sums, and eccentricities to remain consistent with the observed fluxes). The effect of eccentricity can be investigated by comparing models C and D (blue and green, respectively). Forcing the orbit to be circular has a similar effect to placing a tight constraint on the stellar density: model D favors a single-peaked radius ratio distribution and has tighter posterior distributions for the radius sum and inclination than models B and C.

Zoom In Zoom Out Reset image size

The GP hyperparameters are mainly insensitive to the chosen model. For the K2 light curve, the amplitude (A_K2), period (P_{ESS K2}), and white-noise term (σ_K2) are well constrained by the data, but the scale factor (Γ_{ESS K2}) and evolutionary timescale (l_{SE K2}) have Gaussian priors that essentially define their posterior distributions. The evolutionary timescale is related to the half-life of the active regions that are responsible for the observed modulation pattern; however, the exact nature of this relationship is not firmly established. It is worth noting that data covering two evolutionary periods are needed to robustly constrain this parameter, so it is not surprising that it is not well constrained in this case. However, while the evolutionary timescale and scale factor affect the other GP hyperparameters, they do not affect the final transit parameters significantly. For the USNO light curve, the GP amplitude (A_USNO) and timescale (l_USNO) hyperparameters are not well determined; the differences in the timescale posteriors arise from the tightness of the different transit model fits and are therefore driven by the stellar density and eccentricity constraints.

4.1.1.2. Investigating the Residual Variations in the K2 and USNO Light Curves

K2: The model presented above for the K2 light curve seeks to explain the large-scale structure. While it is an acceptable fit to the data, zooming in on small sections of the light curve reveals small-amplitude, short-timescale variations, which could arise from low-level stellar variations or residual systematics. Arguably, adding another GP into the model might account for these, but there would be a lot of flexibility in such a composite GP model, and it is not clear that the residual variations are significant throughout the majority of the light curve. Therefore, to test the significance of these residual variations, we first detrended the light curve with respect to the smoothly varying starspot modulations and then modeled the residual light curve using a Matern-3/2 kernel in our GP-transit model (as previously discussed, this kernel displays a relatively rough behavior, which is suited to the observed residual systematics). It appears as though the information content of the residual light curve is not sufficient to support the additional treatment of the Matern-3/2 kernel. The GP model favors a very small amplitude and either a long timescale over which to vary or a very short one, depending on the initial guess for the input scale. Neither of these is suited to the observed residuals, which have timescales of a few hours, and which are at times clearly observed above the noise. We therefore did not use a more sophisticated GP kernel for the K2 light curve.

USNO: The residuals of the USNO model are clearly seen in the bottom panel of Figure 3. The mean GP model smoothly varies over a relatively long timescale, in contrast to parts of the light curve that show variations over ≲1 hr timescales. To assess the significance of these variations, we reran our model and used the GP itself to estimate the full observational uncertainties (rather than simply allowing it to inflate the uncertainties if required). This second model converges on the same variability and transit parameter values as presented above, which suggests that the information content of the light curve is not sufficient to warrant smaller uncertainties; accordingly, we did not investigate the observed variations further.

It is worth noting, however, that the validity of the observational uncertainties, and hence the analysis above, is subject to the choice of model. Here, the GP model attempts to find the simplest way of explaining the data, given our choice of covariance kernel, and in doing so may opt to increase the uncertainties on each data point rather than the complexity of the chosen family of models if the data set as a whole is unable to support the additional complexity. Accordingly, in model A (presented in Section 4.1.1), the GP opted to slightly inflate the uncertainties by factors of ∼1.37 and 1.17 for the K2 and USNO light curves, respectively (see Table 4).

4.1.2. Sanity Check: Independent K2 Light Curve Analysis

We carry out a second independent light curve analysis to ensure the robustness of our primary analysis. This analysis uses the K2SC-detrended PDC-MAP light curve and simplifies the approach by assuming white noise. The transits are modeled using PyTransit (Parviainen 2015) with a quadratic Mandel–Agol model, and emcee is used for posterior sampling, as in the main analysis (Section 4.1.1). The files associated with this analysis can be found from GitHub at https://github.com/hpparvi/prae1b.

This analysis consider four cases for orbital eccentricity, stellar density, and $\mathrm{log}g$ (assuming a stellar radius of 0.43 R_⊙), listed in Table 5. The sampling space is 10-dimensional, composed of the zero epoch, orbital period, planet–star area ratio, impact parameter, stellar density, $\sqrt{e}\cos \omega$, $\sqrt{e}\sin \omega$, two quadratic limb-darkening parameters, and average white-noise level. All the parameters except the eccentricity and $\mathrm{log}g$ have uninformative priors, and the limb darkening uses the parameterization described by Kipping (2013).

Table 5. Priors on System Parameters for Analysis in Section 4.1.2

Case	Eccentricity	Stellar Density (g cm−3)	$\mathrm{log}g$
A (constrained $\mathrm{log}g$ 2)	${ \mathcal U }(0,0.9)$	${ \mathcal U }(0.5,15)$	${ \mathcal N }(\mu \,=\,4.81,\sigma \,=\,0.08)$
B (constrained $\mathrm{log}g$ 1)	${ \mathcal U }(0,0.9)$	${ \mathcal U }(0.5,15)$	${ \mathcal N }(\mu \,=\,4.5,\sigma \,=\,0.5)$
C (eccentric, unconstrained)	${ \mathcal U }(0,0.9)$	${ \mathcal U }(0.5,15)$	${ \mathcal U }(3,8)$
D (circular)	δ	${ \mathcal U }(0.5,15)$	${ \mathcal U }(3,8)$

Note. ${ \mathcal U }(a,b)$—uniform density with a minimum a and maximum b.${ \mathcal N }(\mu ,\sigma )$—normal density with mean μ and standard deviation σ. δ—Dirac's delta function.

Download table as:  ASCIITypeset image

The posterior sampling starts by creating a parameter vector population of 100 walkers distributed uniformly inside the prior boundaries. We clump the population close to the global posterior maximum using the Differential Evolution global optimization algorithm implemented in PyDE ¹⁵ (Parviainen 2016) and initialize the MCMC sampler with the clumped parameter vector population. The sampler is run over 15,000 iterations and then samples from every 100th iteration, after a burn-in set of 2000 iterations is selected to represent the posterior distribution. The thinning factor of 100 iterations was chosen by inspecting the autocorrelation lengths for individual parameters.

The results from the sanity check analyses are shown in Figure 8, and they agree with the primary analysis. One small difference between the two analyses is that the main analysis finds a broader and near-bimodal distribution for the planet/star radius ratio. We believe that this is because the main analysis uses GPs, whereas the analysis in this section assumes white noise. When using GPs, the transit shape does not constrain the parameter space as strongly, and that additional freedom allows for a larger variety of possible solutions. In the end, we do not believe that this difference represents a challenge to our final system parameters.

Zoom In Zoom Out Reset image size

We fit a grid of BT-SETTL stellar atmosphere models (Allard et al. 2012) to the available ugrizJHK_S + WISE 1–3 catalog photometry, spanning the wavelength range 0.35–12 μm (summarized in Table 6). The reported uncertainties on the Sloan Digital Sky Survey (SDSS) griz photometry are very small (∼0.005 mag), so we rounded up these photometric uncertainties to 0.05 mag to account for photometric variability due to modest activity in this source, as well as zero-point uncertainties of ≈0.02 mag in the SDSS photometry. The free parameters of the fit are T_eff, $\mathrm{log}\,g$, [Fe/H], and A_V; for the latter we adopted a maximum possible value of 0.085 mag based on the full line-of-sight extinction from the Schlegel et al. (1998) dust maps. We interpolate the model grid by 10 K T_eff, from the native gridding of 100 K. We do not interpolate more finely than the native gridding of 0.5 dex in $\mathrm{log}g$ and [Fe/H] because broadband SED fitting is much less sensitive to these parameters compared with T_eff.

The resulting best-fit parameters are T_eff = 3350 ± 50 K, ${A}_{V}={0.000}_{-0.000}^{+0.085}$, $\mathrm{log}g=4.5\pm 0.5$, [Fe/H] = 0.0 ± 0.5. The fit is shown in Figure 9, with a reduced χ² of 6.1. The uncertainties in the fit parameters are estimated according to the usual criterion of Δχ² = 4.72, relative to the best fit, for four fit parameters (e.g., Press et al. 1992), where we first rescaled the χ² so as to make the reduced χ² of the best fit equal to 1. This is equivalent to inflating the measurement errors by a constant factor and has the effect of increasing the parameter uncertainties accordingly. The T_eff from our SED fit is consistent with the previously reported M3.5 spectral type (Adams et al. 2002) to within ∼0.5 spectral subtype.

Zoom In Zoom Out Reset image size

We can use the fitted SED models to compute a bolometric flux at Earth, again according to the Δχ² criterion above. This gives F_bol =  2.10 ± 0.09 × 10⁻¹¹ erg s⁻¹ cm⁻². Note that this F_bol and its uncertainty are robust because the model fit is essentially an interpolation of the observed fluxes that span the majority of the SED. This F_bol, together with the nominal cluster distance of 182 ± 6 pc (van Leeuwen 2009) and the best-fit T_eff, then permits a calculation of the stellar radius via the Stefan–Boltzmann relation, giving R = 0.44 ± 0.03 R_⊙. The radius uncertainty includes the uncertainties on the bolometric flux, the parallax distance, and T_eff, in quadrature.

This is somewhat larger than the radius predicted from the empirical T_eff–R relation of Mann et al. (2016a), which gives R = 0.31 ± 0.04 R_⊙. However, the Praesepe cluster has been reported to be metal-rich, with [Fe/H] ≈ 0.1, but in some cases reported to be as high as [Fe/H] = 0.27 ± 0.10 (Pace et al. 2008). Using this range of [Fe/H] with the T_eff–[Fe/H]–R relation of Mann et al. (2016a) predicts R = 0.37 ± 0.03 R_⊙. Finally, adopting instead the ${{M}_{K}}_{S}$–[Fe/H]–R relation of Mann et al. (2016a), which those authors find has the tightest of all the empirical relations, gives R = 0.42 ± 0.01 R_⊙, in good agreement with our measured value.

Finally, the star is reported both here and in Barrado y Navascués et al. (1998) (Section 2.3) to have a filled-in Hα line, with EW(Hα) ≈ 0 Å indicating the presence of low-level chromospheric activity. Stassun et al. (2012) found that the Hα EW is directly related to the degree of radius inflation in chromospherically active low-mass stars, and their empirical relations give a radius inflation of ∼6% for EW(Hα) = 0 Å. Adjusting the Mann et al. (2016a) ${{M}_{K}}_{S}$–[Fe/H] predicted radius for this amount of inflation and including the uncertainty in the Stassun et al. (2012) relation then gives a final predicted radius of R = 0.44 ± 0.02 R_⊙, in excellent agreement with the bolometric radius value we measure above.

Next, we seek to determine the physical properties of the planet. The planet radius is tied to the stellar radius via the transit radius ratio, and the planet mass is tied to the stellar mass via the RV mass ratio. The host star's radius we have determined above (Section 4.2) to be R = 0.44 ± 0.03 R_⊙, consistent with the (activity-corrected) predicted radius from the empirical relations of Mann et al. (2016a). Similarly, we estimate the host star mass from the same empirical relations. Using the Mann et al. (2016a) relation for M_⋆ versus ${{M}_{K}}_{S}$, we obtain an estimated mass of 0.44 ± 0.04 M_⊙, where the uncertainty includes both the photometric and distance uncertainties on ${{M}_{K}}_{S}$ and the scatter in the empirical relation. This mass, together with the radius from above, then gives a stellar surface gravity of log g = 4.82 ± 0.06.

The ratio of the planet radius to the stellar radius is found to be 0.08, and with the measured stellar radius of 0.44 R_⊙, that yields a planet radius of 0.32 ± 0.02 R_J. We can also calculate an estimated equilibrium temperature for the planet based on the ratio of the stellar radius to planet semimajor axis, stellar temperature, taking the orbit to be circular, and assuming a Bond albedo of 0.3 (similar to Neptune). With those assumptions, we find T_eq = 480 ± 23 K. The full list of derived and calculated planet properties is listed in Table 7.

We have performed an analysis of our RV measurements in order to estimate an upper limit for the mass of the transiting planet. We have six RV measurements, and the standard deviation of these measurements is 1.08 km s⁻¹, which corresponds to a maximum mass for the planet of 6.7 M_Jup (assuming a zero eccentricity). However, in order to take the shape of the RV model into account, we have tested a circular-orbit model with seven free parameters, namely, the systemic RV (uniform prior between 30 and 40 km s⁻¹), orbital period, time of mid-transit and inclination (all three with Gaussian priors based on the light curve analysis), stellar mass (Gaussian prior according to the above estimation), and planet mass (with a uniform prior between 0 and 30 Jupiter masses). We have also included RV jitter (with uniform prior between 0 and 5 km s⁻¹) to account for possible additional sources of white noise. We used the emcee code (Foreman-Mackey et al. 2013) to explore the posterior distribution of the planetary mass based on the current data. We used 50 walkers and 5000 steps per walker. We then used a simple burn-in of half of the chains and merged all of them to get the final posterior distribution. This method is similar to that used in Grunblatt et al. (2016). The result of this analysis yields no detection of the planet (as expected from the RV uncertainties). However, it allows us to set an upper limit for the mass of the transiting object of 0.83 M_Jup at 99.7% confidence given the current data (assuming e = 0). The results of these tests are shown in Figure 10. Note that in this RV phase plot, the maximum expected RV displacement based on the transit model would occur at a phase of 0.75, opposite the offset of the RV observation at that phase, which therefore strongly constrains the allowable mass of the companion.

Zoom In Zoom Out Reset image size

The planet orbiting JS 183 can therefore not yet be confirmed through RV detection of the stellar reflex motion. For a planet with a radius of about 0.32 R_J, we would expect a mass of roughly 15 R_⊕ (see, e.g., Chen & Kipping 2016). For a host star with M_* = 0.44 ± 0.04 M_⊙ and a circular orbit with a period of 10.134588 days, we would expect an RV semiamplitude of about 0.008 km s⁻¹, far below the ≈0.3 km s⁻¹ errors of our HIRES RV observations (Table 2). Although the object would be detectable using current state-of-the-art precision RV techniques having 1–2 m s⁻¹ errors (see, e.g., Plavchan et al. 2015), the star is too faint in practice for these methods to be applied. In such cases, validation of the planet interpretation is typically undertaken either empirically or statistically. Empirical methods use a suite of observations that may be sensitive to detection of certain kinds of blend scenarios (e.g., O'Donovan et al. 2006) that could mimic the observed transit signal. Statistical methods can estimate the probability of various contaminating eclipsing binary scenarios given the source location on the sky relative to the modeled galactic field star population and assumptions regarding the mass function, binary fraction, and binary mass ratios (e.g., Morton 2014).

As we were writing this paper, we became aware of work by Obermeier et al. (2016) and Mann et al. (2017), who independently discovered and followed up the same K2 time series data on EPIC 211916756. The thorough false-positive analysis presented in the Obermeier et al. (2016) paper need not be repeated here. From empirical considerations, no companions are detected in high spatial resolution imaging or in spectroscopy. From statistical considerations, a false-positive probability calculation rules out line-of-sight blended and bound hierarchical eclipsing binary systems. Remaining for consideration is the planet hypothesis.