New low-mass eclipsing binary systems in Praesepe discovered by K2

TL;DR

This paper reports the discovery and detailed analysis of four low-mass eclipsing binary systems in the Praesepe cluster, introducing a new modeling method and testing stellar evolution models with high-precision measurements.

Contribution

It presents a new Gaussian process eclipsing binary model, GP-EBOP, and a method for simultaneous temperature and distance determination, along with precise measurements of low-mass binaries in Praesepe.

Findings

Three systems are detached and double-lined, with precise mass and radius measurements.

The observed stellar parameters are more consistent with PARSEC models than BHAC15.

Discovery of the first brown dwarf transiting a mid M-dwarf in a young open cluster.

Abstract

We present the discovery of four low-mass ($M<0.6$ $M_\odot$) eclipsing binary (EB) systems in the sub-Gyr old Praesepe open cluster using Kepler/K2 time-series photometry and Keck/HIRES spectroscopy. We present a new Gaussian process eclipsing binary model, GP-EBOP, as well as a method of simultaneously determining effective temperatures and distances for EBs. Three of the reported systems (AD 3814, AD 2615 and AD 1508) are detached and double-lined, and precise solutions are presented for the first two. We determine masses and radii to 1-3% precision for AD 3814 and to 5-6% for AD 2615. Together with effective temperatures determined to $\sim$50 K precision, we test the PARSEC v1.2 and BHAC15 stellar evolution models. Our EB parameters are more consistent with the PARSEC models, primarily because the BHAC15 temperature scale is hotter than our data over the mid M-dwarf mass range probed. Both ADs 3814 and 2615, which have orbital periods of 6.0 and 11.6 days, are circularized but not synchronized. This suggests that either synchronization proceeds more slowly in fully convective stars than the theory of equilibrium tides predicts or magnetic braking is currently playing a more important role than tidal forces in the spin evolution of these binaries. The fourth system (AD 3116) comprises a brown dwarf transiting a mid M-dwarf, which is the first such system discovered in a sub-Gyr open cluster. Finally, these new discoveries increase the number of characterized EBs in sub-Gyr open clusters by 20% (40%) below $M<1.5$ $M_{\odot}$ ($M<0.6$ $M_{\odot}$).

Tables

Table 1. Table 1: Names, coordinates, properties and membership information for the four newly identified EBs.

Property	Units	AD 3814	AD 2615	AD 3116	AD 1508	Refs.
EPIC		211972086	212002525	211946007	212009427
2MASS		J08504984+1948364	J08394203+2017450	J08423943+1924520	J08312987+2024374
Other names		…	…	HSHJ 430	…	1
RA	J2000.0	08:50:49.84	08:39:42.03	08:42:39.43	08:31:29.87
Dec	J2000.0	+19:48:36.4	+20:17:45.0	+19:24:51.9	+20:24:37.5
$u$	AB	$21.009\pm 0.093$	$21.747\pm 0.185$	$22.290\pm 0.190$	$18.102\pm 0.014$	2
$g$	AB	$18.769\pm 0.008$	$19.416\pm 0.012$	$19.646\pm 0.014$	$15.540\pm 0.004$	2
$r$	AB	$17.299\pm 0.006$	$17.905\pm 0.007$	$18.206\pm 0.007$	$14.151\pm 0.004$	2
$i$	AB	$15.803\pm 0.005$	$16.324\pm 0.004$	$16.675\pm 0.005$	$13.700\pm 0.001$	2
$z$	AB	$14.999\pm 0.005$	$15.456\pm 0.006$	$15.845\pm 0.006$	$12.905\pm 0.004$	2
V	Vega	17.80	18.46	18.73	14.79	3
$J$	Vega	$13.529\pm 0.026$	$14.027\pm 0.021$	$14.348\pm 0.032$	$11.674\pm 0.022$	4
$H$	Vega	$12.911\pm 0.024$	$13.456\pm 0.026$	$13.769\pm 0.037$	$10.949\pm 0.023$	4
$K_s$	Vega	$12.651\pm 0.022$	$13.136\pm 0.034$	$13.499\pm 0.043$	$10.767\pm 0.020$	4
WISE 1	Vega	$12.478\pm 0.024$	$12.938\pm 0.024$	$13.299\pm 0.029$	$10.677\pm 0.023$	4
WISE 2	Vega	$12.291\pm 0.026$	$12.773\pm 0.031$	$13.096\pm 0.039$	$10.638\pm 0.021$	4
Spectral type	M sub-type	$3.4\pm 0.1$	$4.0\pm 0.1$	$3.9\pm 0.3$	$0.1\pm 0.1$	5
H$\alpha$ emission	Å	2.4–3.5, …	3.0–4.3, 10.7	3.1–5.2, 4.6	2.0–2.1, …	6,7
RA proper motion, ${\mu }_{\alpha }$	mas yr-1	-37.5	-39.3	-37.5	-37.3	5
Dec proper motion, ${\mu }_{\delta }$	mas yr-1	-14.1	-11.6	-8.2	-16.7	5
Membership probability	%	97.9	99.7	99.1	98.3	5

Table 2. Table 2: Radial velocities derived from Keck/HIRES spectra for ADs 3814, 2615, 3116 and 1508 ( top  to bottom ).

Epoch			S/N	RV   (km s-1)
UT date	BJD	Phase *	7500 Å	Primary	Secondary
……………………………………………….   AD 3814   ……………………………………………….
2015 12 24	2457381.15090	0.607	16	$54.08\pm 0.77$	$-6.83\pm 0.93$
2015 12 29	2457386.14539	0.437	16	$21.42\pm 0.76$	$58.70\pm 1.10$
2016 02 02	2457420.89940	0.214	16	$0.10\pm 0.75$	$95.82\pm 1.13$
2016 02 03	2457421.92652	0.385	15	$12.91\pm 0.76$	$77.18\pm 0.93$
2016 05 17	2457525.80479	0.652	14	$60.96\pm 0.37$	$-19.29\pm 1.15$
2016 12 22	2457744.96970	0.084	15	$15.72\pm 0.35$	$65.38\pm 0.56$
2016 12 26	2457748.97551	0.750	13	$68.56\pm 1.04$	$-28.91\pm 1.18$
2017 01 13	2457766.85771	0.723	12	$67.17\pm 0.40$	$-29.35\pm 0.56$
……………………………………………….   AD 2615   ……………………………………………….
2015 12 29	2457386.16741	0.039	13	$26.03\pm 0.86$	$43.50\pm 0.77$
2016 05 17	2457525.78402	0.059	13	$20.45\pm 0.76$	$45.69\pm 0.80$
2016 05 20	2457528.78074	0.317	14	$-1.28\pm 0.60$	$65.74\pm 0.60$
2016 10 14	2457676.07100	0.997	10	$35.22\pm 0.29$
2016 12 22	2457745.03382	0.935	13	$49.63\pm 0.48$	$20.66\pm 0.41$
2017 01 13	2457766.90914	0.818	5	$72.09\pm 0.53$	$5.63\pm 0.60$
……………………………………………….   AD 3116   ……………………………………………….
2016 02 02	2457420.92116	0.102	12	$26.28\pm 0.82$	—
2016 02 03	2457421.90606	0.599	12	$40.47\pm 0.83$	—
2016 05 17	2457525.76250	0.978	12	$39.75\pm 0.59$	—
2016 05 20	2457528.75747	0.488	13	$27.46\pm 0.54$	—
2016 10 14	2457676.09435	0.796	13	$55.96\pm 0.30$	—
2016 12 22	2457744.98886	0.542	12	$31.91\pm 0.44$	—
2017 01 13	2457766.87986	0.583	6	$37.23\pm 0.77$	—
……………………………………………….   AD 1508   ……………………………………………….
2016 12 22	2457745.047527495	0.971	40	$50.62\pm 1.40$	$16.36\pm 1.57$
2016 12 26	2457748.953315984	0.479	30	$21.66\pm 2.79$	$42.37\pm 3.23$
2017 01 13	2457766.842635326	0.970	40	$52.56\pm 2.19$	$18.69\pm 1.82$

Table 3. Table 3: Spectroscopic light ratios and quadratic limb darkening priors applied in the GP–EBOP modeling for ADs 3814, 2615, 3116 and 1508.

System	Spectroscopic light ratio		Limb darkening coefficients and assumed model atmosphere parameters *
	BJD	$l_{\sec }/l_{\mathrm{pri}}$	Component	$\mu$	${\mu }^{\prime }$	$T_{\mathrm{eff}}$ (K)	$\log g$ (cgs)
AD 3814	245 7766.9	${0.41}_{-0.19}^{+0.25}$	Pri	$0.46\pm 0.13$	$0.21\pm 0.46$	$3200\pm 200$	$4.9\pm 0.1$
			Sec	$0.49\pm 0.24$	$0.23\pm 0.76$	$3100\pm 200$	$5.0\pm 0.1$
AD 2615	245 7766.9	${1.13}_{-0.20}^{+0.24}$	Pri & Sec	$0.46\pm 0.13$	$0.21\pm 0.46$	$3200\pm 200$	$4.9\pm 0.1$
AD 3116	—	—	Pri	$0.46\pm 0.13$	$0.21\pm 0.46$	$3200\pm 200$	$4.9\pm 0.1$
			Sec	$0.68\pm 0.17$	$0.17\pm 0.46$	$2500\pm 200$	$5.0\pm 0.1$
AD 1508	245 7745.0	${0.63}_{-0.26}^{+0.41}$	Pri & Sec	$0.47\pm 0.14$	$0.20\pm 0.31$	$3700\pm 200$	$4.8\pm 0.1$

Table 4. Table 4: Effective temperatures and distance values for each EB estimated from SED modeling and the empirical relations of Mann et al. ( 2015 ) and David et al. ( 2016a ) .

Method ∗	Model †	$T_{\mathrm{eff}}$  ‡		Distance
		Primary	Secondary
		(K)	(K)	(pc)
…………………………………….. AD 3814 ………………………………..
SED	PHOENIX	$3193\pm 17$	$3085\pm 21$	${168.8}_{-7.3}^{+6.1}$
SED	BT-SETTL	$3230\pm 36$	$3121\pm 35$	$172.1\pm 9.3$
ER	M15	$3241\pm 76$	$3013\pm 79$	$172\pm 12$
ER	D16	3251	3023
SED	Combined	${3211}_{-36}^{+54}$	${3103}_{-39}^{+53}$	${170.4}_{-8.9}^{+11.0}$
…………………………………….. AD 2615 ………………………………..
SED	PHOENIX	$3132\pm 21$	$3112\pm 20$	$177.2\pm 7.9$
SED	BT-SETTL	$3172\pm 37$	$3150\pm 37$	$181\pm 11$
ER	M15	$3187\pm 75$	$3145\pm 90$	$177\pm 15$
ER	D16	3197	3156
SED	Combined	${3152}_{-40}^{+57}$	${3131}_{-38}^{+56}$	${179}_{-10}^{+13}$
…………………………………….. AD 3116 ………………………………..
SED	BT-SETTL	$3184\pm 29$	$1639\pm 248$	—
ER	M15	$3237\pm 74$	$880\pm 217$	$183\pm 14$
ER	D16	3236	880
…………………………………….. AD 1508 ………………………………..
SED	PHOENIX	$3754\pm 78$	$3679\pm 121$	$164\pm 22$
SED	BT-SETTL	$3779\pm 87$	$3706\pm 117$	$167\pm 23$
ER	M15	$3738\pm 76$	$3639\pm 284$	$156\pm 28$
ER	D16	3746	3649
SED	Combined	${3767}_{-85}^{+99}$	${3693}_{-135}^{+122}$	${166}_{-23}^{+25}$

Table 5. Table 5: Published double-lined eclipsing binary systems in sub-Gyr open clusters where both components are below 1.5 M ⊙ , ordered by ascending primary mass.

Name	$M_{\mathrm{pri}}$	$M_{\sec }$	$R_{\mathrm{pri}}$	$R_{\sec }$	Cluster a	Age	Year	Refs.
	($M_{\odot }$)	($M_{\odot }$)	($R_{\odot }$)	($R_{\odot }$)		(Myr)
EPIC 203868608	$0.02216\pm 0.00045$	$0.02462\pm 0.00055$	$0.2823\pm 0.0051$	$0.2551\pm 0.0036$	Upper Sco	5–10	2016	1
2MJ0535-05	$0.0572\pm 0.0033$	$0.0366\pm 0.0022$	$0.690\pm 0.011$	$0.540\pm 0.009$	ONC	1–2	2006	2,3
EPIC 203710387	$0.1183\pm 0.0028$	$0.1076\pm 0.0031$	$0.417\pm 0.010$	$0.450\pm 0.012$	Upper Sco	5–10	2015	4,1
JW 380	$0.262\pm 0.025$	$0.151\pm 0.013$	$1.189\pm 0.175$	$0.897\pm 0.170$	ONC	1–2	2007	5
HCG 76	$0.2768\pm 0.0072$	$0.3020\pm 0.0073$	$0.319\pm 0.036$	$0.34\pm 0.11$	Pleiades	125	2016	6
UScoCTIO 5	$0.3336\pm 0.0022$	$0.3200\pm 0.0022$	$0.862\pm 0.012$	$0.852\pm 0.013$	Upper Sco	5–10	2015	7,1
Par 1802	$0.391\pm 0.032$	$0.385\pm 0.032$	$1.73\pm 0.015$	$1.62\pm 0.015$	ONC	1–2	2008	8,9
MHO 9	$0.41\pm 0.18$	$0.172\pm 0.069$	$0.46\pm 0.11$	$0.321\pm 0.060$	Pleiades	125	2016	6
2MJ0446+19	$0.47\pm 0.05$	$0.19\pm 0.02$	$0.56\pm 0.02$	$0.21\pm 0.01$	NGC 1647	150	2006	10
CoRoT 223992193	$0.668\pm 0.012$	$0.4953\pm 0.0073$	$1.295\pm 0.040$	$1.107\pm 0.044$	NGC 2264	3–6	2014	11
MML 53	$0.994\pm 0.030$	$0.857\pm 0.026$	$2.201\pm 0.071$ b		UCL	15	2010	12,13
HD144548	$0.984\pm 0.007$	$0.944\pm 0.017$	$1.319\pm 0.010$	$1.330\pm 0.010$	Upper Sco	5–10	2015	14
	$1.44\pm 0.04$ c		$2.41\pm 0.03$ c
V1174 Ori	$1.006\pm 0.013$	$0.7271\pm 0.0096$	$1.338\pm 0.011$	$1.063\pm 0.011$	Ori OB 1c	5–10	2004	15
V818 Tau	$1.06\pm 0.01$	$0.90\pm 0.02$	$0.76\pm 0.01$	$0.77\pm 0.01$	Hyades	600–800	2002	16
RXJ 0529.4$+$0041A	$1.27\pm 0.01$	$0.93\pm 0.01$	$1.44\pm 0.10$	$1.35\pm 0.10$	Ori OB 1a	7–13	2000	17,18,19
NP Per	$1.3207\pm 0.0087$	$1.0456\pm 0.0046$	$1.372\pm 0.013$	$1.229\pm 0.013$	Per OB 2	6–15	2016	20
ASAS J0528$+$03	$1.375\pm 0.028$	$1.329\pm 0.020$	$1.83\pm 0.07$	$1.73\pm 0.07$	Ori OB 1a	7–13	2008	21
Praesepe systems published in this paper
AD 3814	$0.3813\pm 0.0074$	$0.2022\pm 0.0045$	$0.3610\pm 0.0033$	$0.2256\pm 0.0063$	Praesepe	600–800	2017	this work
AD 2615	$0.212\pm 0.012$	$0.255\pm 0.013$	$0.233\pm 0.013$	$0.267\pm 0.014$	Praesepe	600–800	2017	this work
AD 1508	$0.45\pm 0.19$	$0.53\pm 0.21$	$0.548\pm 0.099$	$0.45\pm 0.10$	Praesepe	600–800	2017	this work

Table 6. Table 6: Fitted and derived parameters of the models applied to AD 3814, AD 2615, AD 3116 and AD 1508.

Parameter	Symbol	Unit	Value
			AD 3814	AD 2615	AD 3116	AD 1508
Eclipse parameters
Sum of radii	$\left(R_{\mathrm{pri}}+R_{\sec }\right)/a$		${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.05044}}_{-0.00055}^{+0.00069}$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.02979}\pm 0.00034$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.0845}}_{-0.0053}^{+0.0066}$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.1774}}_{-0.0076}^{+0.0066}$
Radius ratio	$R_{\sec }/R_{\mathrm{pri}}$		${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.624}}_{-0.010}^{+0.017}$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}1.15}\pm 0.11$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.3599}}_{-0.0128}^{+0.0094}$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.83}\pm 0.24$
Orbital inclination	$i$	∘	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}89.177}}_{-0.064}^{+0.051}$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}88.996}\pm 0.013$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}88.41}}_{-0.42}^{+0.49}$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}80.54}}_{-0.39}^{+0.46}$
Orbital period	$P$	days	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}6.015717}\pm 0.000013$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}11.615254}\pm 0.000073$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}1.9827960}\pm 0.0000060$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}1.5568370}}_{-0.0000090}^{+0.0000100}$
Time of eclipse center	$T_{\mathrm{prim}}$	BJD	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}2457178.982842}\pm 0.000059$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}2457176.63998}\pm 0.00019$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}2457178.817792}\pm 0.000080$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}2457147.26784}\pm 0.00026$
	$\sqrt{e}\cos \omega$		$-{0.0301}_{-0.0057}^{+0.0103}$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.0337}}_{-0.0128}^{+0.0067}$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.364}}_{-0.026}^{+0.016}$	$-{0.0081}_{-0.0094}^{+0.0069}$
	$\sqrt{e}\sin \omega$		$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.031}\pm 0.034$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.020}\pm 0.052$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.04}\pm 0.14$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.010}}_{-0.041}^{+0.049}$
Central surface brightness ratio	$J_{\mathrm{K2}}$		$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.748}\pm 0.034$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.950}\pm 0.060$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.0051}}_{-0.0036}^{+0.0049}$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.90}}_{-0.21}^{+0.30}$
Primary linear LDC *	$u_{\mathrm{pri}\mathrm{K2}}$		$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.54}\pm 0.20$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.51}\pm 0.13$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.66}\pm 0.16$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.48}\pm 0.14$
Primary non-linear LDC *	$u_{\mathrm{pri}\mathrm{K2}}^{\prime }$		$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.24}\pm 0.29$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.31}\pm 0.22$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.03}\pm 0.23$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.22}\pm 0.24$
Secondary linear LDC *	$u_{\sec \mathrm{K2}}$		$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.39}\pm 0.11$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.41}\pm 0.13$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.43}\pm 0.13$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.44}\pm 0.14$
Secondary non-linear LDC *	$u_{\sec \mathrm{K2}}^{\prime }$		$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.12}\pm 0.21$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.04}}_{-0.18}^{+0.25}$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.12}\pm 0.26$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.12}\pm 0.23$
Out-of-eclipse variability parameters
Amplitude	$A_{\mathrm{K2}}$	%	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.00785}}_{-0.00057}^{+0.00077}$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.0179}\pm 0.0023$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.0103}}_{-0.0020}^{+0.0033}$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.136}}_{-0.136}^{+0.085}$
Timescale of SqExp term	$l_{\mathrm{SE}\mathrm{K2}}$	days	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}8.55}}_{-0.36}^{+0.43}$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}16.97}\pm 0.97$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}7.32}}_{-0.67}^{+0.90}$	—
Scale factor of ExpSine2 term	${\Gamma }_{\mathrm{ESS}\mathrm{K2}}$	days	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}11.5}\pm 4.0$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}9.56}}_{-0.82}^{+0.96}$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.55}}_{-0.20}^{+0.27}$	—
Period of ExpSine2 term	$P_{\mathrm{ESS}\mathrm{K2}}$	days	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}7.375}}_{-0.069}^{+0.059}$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}12.150}}_{-0.062}^{+0.074}$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}2.252}\pm 0.020$	—
Timescale of Matern32 term	$l_{\mathrm{M32}\mathrm{K2}}$	days	—	—	—	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}223.6}\pm 2.3$
White noise scale factor	${\sigma }_{\mathrm{K2}}$		$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}1.901}\pm 0.039$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}1.448}\pm 0.02$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}1.363}\pm 0.019$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}1.329}}_{-0.023}^{+0.031}$
Radial velocity parameters
Systemic velocity	$V_{\mathrm{sys}}$	km s-1	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}33.60}\pm 0.24$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}34.91}\pm 0.39$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}34.93}}_{-0.53}^{+0.61}$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}33.1}\pm 1.7$
Primary RV semi-amplitude	$K_{\mathrm{pri}}$	km s-1	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}33.90}\pm 0.39$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}39.86}}_{-0.88}^{+0.80}$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}18.66}}_{-1.00}^{+0.95}$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}98}\pm 15$
Secondary RV semi-amplitude	$K_{\sec }$	km s-1	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}63.93}\pm 0.49$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}33.12}}_{-0.89}^{+0.83}$	—	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}84}\pm 14$
HIRES jitter term	${\sigma }_{\mathrm{HIRES}}$	km s-1	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.50}}_{-0.30}^{+0.37}$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.95}}_{-0.35}^{+0.48}$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.93}}_{-0.62}^{+1.18}$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}1.6}}_{-1.2}^{+2.5}$
Fundamental parameters
Primary mass	$M_{\mathrm{pri}}$	M⊙	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.3813}\pm 0.0074$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.212}\pm 0.012$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.276}\pm 0.020$ a	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.45}}_{-0.14}^{+0.19}$
Secondary mass	$M_{\sec }$	M⊙	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.2022}\pm 0.0045$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.255}\pm 0.013$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.0517}\pm 0.0041{\left(54.2\pm 4.3\right)}^{b,c}$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.53}}_{-0.16}^{+0.22}$
Primary radius	$R_{\mathrm{pri}}$	R⊙	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.3610}\pm 0.0033$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.233}\pm 0.013$	$0.29\pm 0.08$ d	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.549}}_{-0.082}^{+0.099}$
Secondary radius	$R_{\sec }$	R⊙	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.2256}}_{-0.0049}^{+0.0063}$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.267}\pm 0.014$	$0.10\pm 0.03{\left(1.02\pm 0.28\right)}^{c,e}$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.454}}_{-0.101}^{+0.094}$
Primary effective temperature	$T_{\mathrm{pri}}$	K	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}3211}}_{-36}^{+54}$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}3152}}_{-40}^{+57}$	$3191\pm 27$	${3767}_{-85}^{+99}$
Secondary effective temperature	$T_{\sec }$	K	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}3103}}_{-39}^{+53}$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}3131}}_{-38}^{+56}$	${1669}_{-258}^{+244}$	${3693}_{-135}^{+122}$
Mass sum	$M_{\mathrm{pri}}+M_{\sec }$	M⊙	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.583}\pm 0.011$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.468}\pm 0.023$	—	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.98}}_{-0.29}^{+0.38}$
Radius sum	$R_{\mathrm{pri}}+R_{\sec }$	R⊙	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.5868}}_{-0.0073}^{+0.0084}$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.4991}}_{-0.0102}^{+0.0096}$	—	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}1.00}\pm 0.13$
Semi-major axis	$a$	R⊙	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}11.630}\pm 0.073$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}16.75}\pm 0.28$	—	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}5.67}\pm 0.65$
Eccentricity	$e$		${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.00194}}_{-0.00057}^{+0.00253}$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.00254}}_{-0.00078}^{+0.00406}$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.146}}_{-0.016}^{+0.024}$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.00108}}_{-0.00078}^{+0.00347}$
Longitude of periastron	$\omega$	∘	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}116.0}\pm 39.0$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}27.0}\pm 69.0$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}5}\pm 20$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}91.0}\pm 29.0$
Primary surface gravity	$\log g_{\mathrm{pri}}$	(cm s-2)	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}4.9040}}_{-0.0064}^{+0.0073}$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}5.031}\pm 0.048$	—	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}4.61}\pm 0.13$
Secondary surface gravity	$\log g_{\sec }$	(cm s-2)	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}5.037}}_{-0.026}^{+0.019}$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}4.993}}_{-0.035}^{+0.042}$	—	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}4.84}\pm 0.20$
Primary synchronized velocity	$V_{\mathrm{pri}\mathrm{sync}}$	km s-1	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}3.036}\pm 0.028$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}1.014}\pm 0.058$	—	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}17.8}\pm 3.2$
Secondary synchronized velocity	$V_{\sec \mathrm{sync}}$	km s-1	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}1.898}}_{-0.041}^{+0.053}$	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}1.162}}_{-0.059}^{+0.050}$	—	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}14.7}\pm 3.3$
Synchronization timescale	$t_{\mathrm{sync}}$	Myr	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}27.29}\pm 0.49$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}152.6}\pm 4.7$	—	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.0510}}_{-0.0090}^{+0.0120}$
Circularization timescale	$t_{\mathrm{circ}}$	Gyr	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}17.30}\pm 0.10$	$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}467.6}\pm 1.5$	—	${\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0.01040}}_{-0.00013}^{+0.00033}$

Equations

$$\mathbf{F}\sim \mathcal{N}(\mu ,\mathbf{K}).$$

$${\mathbf{K}}_{ij}=k(t_i,t_j)+k_w(i,j)$$

$$k\mathrm{QPESS}(t_i,t_j)=A^2\exp \left[-\Gamma {\sin }^2\left(\frac{\pi \left|t_i-t_j\right|}{P}\right)-\frac{{\left(t_i-t_j\right)}^2}{2l^2}\right].$$

$$k_w(i,j)={\sigma }^2{\delta }_{ij}$$

$$k_{M32}(t_i,t_j)=A^2\left(1+\frac{\sqrt{3}\left|t_i-t_j\right|}{l}\right)\exp \left(-\frac{\sqrt{3}\left|t_i-t_j\right|}{l}\right)$$

Full text

binaries: eclipsing – binaries: spectroscopic – stars: fundamental parameters – stars: low-mass – brown dwarfs

††facilities: Kepler/K2, Keck/HIRES, SDSS, 2MASS, WISE.††software: astropy (Astropy Collaboration et al., 2013), emcee (Foreman-Mackey et al., 2013), george (Ambikasaran et al., 2014).

1 Introduction

Stellar evolution theory underpins much of observational astrophysics, yet significant uncertainties remain at low masses ($M\lesssim 0.8$ $M_{\odot}$) and young ages ($t\lesssim 1$ Gyr). Unfortunately, this mass and age range is also where observational constraints are scarce. The fundamental goal of stellar evolution theory is to accurately predict the observables (radius, temperature and luminosity) for a star of given mass, age and metallicity. The evolutionary pathway of a star is governed primarily by its mass, which is accessible only through study of gravitational interactions such as in binary or higher order multiple star systems. For eclipsing binaries (EBs) the ratio of the radii of the stars is attainable. Eclipsing binaries are particularly important objects if they are also detected as double-lined systems in spectra, as the individual masses and radii of both stars can be extracted from the combined light curves and radial velocity curves of the system. Radii can also be measured directly using interferometric techniques, but only for the brightest of nearby stars. When the inferred mass and radius values reach a precision of a few percent or less, they provide one of the strongest observational tests of stellar evolution theory available (e.g. Torres et al., 2010; Stassun et al., 2014).

Open clusters are fruitful astrophysical laboratories given that their members share broad coevality, composition and distance. The detection of multiple EBs in a given cluster, with each member of each pair sharing the same age and metallicity but spanning a range of masses, offers a particularly strong test of stellar evolution theory. The pursuit of EB parameters, among other science goals, has motivated numerous programmes to target open clusters via time-series photometry, e.g. the ground-based Monitor, PTF Orion and YETI projects (Aigrain et al., 2007; van Eyken et al., 2011; Neuhäuser et al., 2011), and space-based observations with CoRoT and Spitzer (Gillen et al., 2014; Morales-Calderón et al., 2012). Furthermore, since March 2014, the re-purposed Kepler mission, K2 (Howell et al., 2014), has targeted a number of star forming regions and young (sub-Gyr) open clusters across the ecliptic for $\sim$80 days each. To date, the nearby $\rho$ Ophiuchi star forming region and Upper Scorpius young OB association ($\sim$1 and 5-10 Myr, respectively) were observed in Campaign 2, as were the Pleiades and Hyades open clusters ($\sim$125, 600–800 Myr, respectively) in Campaign 4, and Praesepe (600–800 Myr) in Campaign 5.

The Praesepe open cluster, also known as the Beehive cluster or M44, was targeted by K2 in Campaign 5 (April–July 2015). Praesepe is a relatively nearby, metal-rich, several hundred Myr cluster hosting $>$1000 high probability members ($>$80%) and more than 100 candidate members ($>$50% probability) (Kraus & Hillenbrand, 2007; Rebull et al., 2017). Given its richness and proximity, Praesepe is a well-studied benchmark cluster. The parallaxes of bright Praesepe members in the GAIA DR1 suggest a distance of $182.8\pm 1.7\pm 14$ pc, where the first error represents the uncertainty on the cluster center determination and the second reflects the observed radial spread of high probability members on the sky (Gaia Collaboration et al., 2017). This is in agreement with the commonly quoted Hipparcos distance to the cluster, $181.5\pm 6.0$ pc (van Leeuwen, 2009). The cluster has a low reddening along the line of sight of $E(B-V)=0.027\pm 0.004$ (Taylor, 2006). Metallicity estimates typically fall within the range [Fe/H] $\sim$ 0.12–0.16 (e.g. Boesgaard et al., 2013; Yang et al., 2015; Netopil et al., 2016) but can be as high as [Fe/H] = $0.27\pm 0.10$ (e.g. Pace et al., 2008). The age of Praesepe is estimated in the range $\sim$600–900 Myr (e.g. Adams et al., 2002) with traditional estimates typically falling at the lower end, often through association with the Hyades (e.g. Salaris et al., 2004). More recently, however, Brandt & Huang (2015b) included stellar rotation to conclude that the upper main sequences of both Praesepe and the Hyades were consistently well-fit at an age of $\sim$750–800 Myr. The age of Praesepe is further discussed in section 6.2.3.

The binary fraction within the cluster has been extensively studied. Pinfield et al. (2003) noted that binaries in Praesepe appear to favor similar-mass systems. Boudreault et al. (2012) focused on the low-mass population, finding binary frequencies of: $25.6\pm 3.0$% between 0.2 $<$ M $<$ 0.45 $M_{\odot}$, $19.6\pm 3.0$% between 0.1 $<$ M $<$ 0.2 $M_{\odot}$ and $23.2\pm 5.6$% between 0.07 $<$ M $<$ 0.1 $M_{\odot}$. Wang et al. (2014) analyzed the full Praesepe membership to find a binary occurrence rate of 20–40%. Furthermore, a significant population of binaries and higher order systems were identified by Khalaj & Baumgardt (2013), who propose a binary fraction of $35\pm 5$% in the mass range 0.6–2.2 $M_{\odot}$, assuming mass-dependent pairing of primary stars following the results of recent star formation simulations (e.g. Bate, 2009).

This paper presents the characterization of four high-probability, low-mass eclipsing binary members of Praesepe. §2 describes the sources and previous literature characterization. In §3 we detail the photometric and spectroscopic observations. In §4 we present a modified eclipsing binary model for detached systems, GP–EBOP, and describe the light curve and radial velocity analyzes. We then present the results for each system in §5. In §6 we present an updated method to simultaneously determine the effective temperatures of both stars as well as the distance to an EB system, before discussing these new EBs in the context of calibrating stellar evolution models, and informing tidal evolution theory in close binaries. Finally, we conclude in §7.

2 New Eclipsing Binaries Among Praesepe Members

Half a dozen deep proper motion surveys of Praesepe have been published since 2000 (Adams et al., 2002; Kraus & Hillenbrand, 2007; Baker et al., 2010; Boudreault et al., 2012; Khalaj & Baumgardt, 2013; Wang et al., 2014). Three of our four EBs are considered Praesepe members in at least four of those six studies (AD 3814, 2615 and 3116). Our fourth EB (AD 1508) is identified as a Praesepe member in only two of those studies.

In the top panel of Figure 1, we show where these four objects fall in a V vs. V-Ks color-magnitude diagram, where we have derived V-Ks estimates based on a conversion (Rebull et al., 2017) from G-Ks, where G is the star’s magnitude in the Gaia DR1 catalog. All four stars have photometry consistent with Praesepe membership. AD 1508 is the earliest type (brightest) of the four; it is located well above the single star main-sequence locus, suggesting that it is a nearly equal mass binary. AD 3116 and 3814 are located nearly on the single star main-sequence locus, and so their binary companions are presumably very low mass. AD 2615 is displaced about 0.4 mag above the single star locus, and so is likely to have an intermediate mass binary companion.

Three of the four stars have published spectral types: AD 3814 - M5 (West et al., 2011); AD 2615 - M4.0 (Adams et al., 2002), M5 (West et al., 2011); and AD 3116 - M4.5 (Adams et al., 2002), M3.9 (Kafka & Honeycutt 2006. These spectral types are broadly consistent with their V-Ks colors. All four systems have spectral types estimated from photometry (Kraus & Hillenbrand, 2007): AD 3814 - M$3.4\pm 0.1$; AD 2615 - M$4.0\pm 0.1$; AD 3116 - M$3.9\pm 0.1$; and AD 1508 - M$0.1\pm 0.1$. As these form a homogeneous set for our EBs, we adopt these spectral types here. For each system, properties extracted from the literature are reported in Table 1.

In the bottom panel of Figure 1 we show the Praesepe V vs. V-Ks color vs. rotation period diagram and indicate our four systems in red. Given the wide spread in rotation periods for mid–M dwarfs, ADs 3814, 2615 and 3116 all lie along the single star trend, but the early–M dwarf AD 1508 lies far below the single star trend with a short rotation period.

3 Observations

3.1 Photometry

We proposed targets for the K2 Campaign 5 observations, which included Praesepe, as part of the K2 Young Suns Survey (PI Stauffer). Targets were collated through merging various proper motion surveys (Klein Wassink, 1927; Jones & Cudworth, 1983; Jones & Stauffer, 1991; Kraus & Hillenbrand, 2007; Wang et al., 2014) with published BVRI photometry (Mermilliod et al., 1990; Stauffer, 1982, and references therein). The K2FOV tool was used to select targets falling ‘on silicon’ and we further limited our proposal to stars with spectral type later than F0 (i.e. possessing outer convective envelopes) and brighter than $R<17$. This gave 477 high-probability Praesepe targets in total. In addition to our proposed systems we also investigated light curves of Praesepe candidates from other K2 programs.

The K2 observations of Praesepe spanned 27 April – 10 July 2015 and the FoV centerd on 08:40:38 +16:49:47. Given the typical 30-minute cadence of Kepler observations, this resulted in $\sim$3300 data points for each target. Short cadence (1 min) observations are also possible for a small number of targets but all systems presented here were observed in standard long cadence mode. We discuss our method to reduce the K2 photometry in §3.1.1. For objects showing the signatures of eclipses in the K2 time series photometry, we cross-referenced the EPIC identifiers with literature information in order to determine basic system properties (see §3.1.2) and to identify which systems to pursue with high dispersion spectroscopy (see §3.2).

3.1.1 K2 data detrending and eclipse detection

We started from the Simple Aperture Photometry (SAP) light curves, which were made available at the Mikulski Archive for Space Telescopes (MAST) as part of K2 Data Release 7111See https://keplerscience.arc.nasa.gov/

k2-data-release-notes.html#k2-campaign-5 for details.. We used the k2sc pipeline (Aigrain et al., 2016) to correct the light curves for systematics caused by the quasi-periodic rolling motion of the spacecraft, while preserving the intrinsic variability of the target stars. k2sc works by modeling the SAP flux as the sum of two smooth, random functions: one depending on the star’s position on the detector, and one depending on time, plus white noise. The position component represents instrumental systematics associated with the satellite’s pointing variations (mainly intra- and inter-pixel sensitivity variations), while the time component represents the star’s intrinsic variability, plus any long-term instrumental effects not accounted for by the position component. Both components are modeled using Gaussian Process (GP) regression (see §4.1 for further details and references on GPs). While both components are initially treated as aperiodic, a quasi-periodic GP is automatically used for the time component if the light curve shows any evidence of periodic behavior after a first pass treatment with default parameters.

A careful treatment of outliers ensures that k2sc mostly preserves short-duration events such as planetary transits or stellar eclipses. However, once the eclipses were identified (by visual examination) in the four systems discussed in the present paper, their light curves were re-processed using k2sc’s periodic mask option. This option enables the user to supply the period, epoch and duration of the eclipses, and any in-eclipse points are then ignored when training the GP model. In effect, we are using the k2sc GP model to interpolate in both flux and position space to the times affected by the eclipses, thereby providing a model prediction for the total system flux across each eclipse. In our analysis we use the k2sc light curve that has been detrended for instrument systematics but which retains the stellar variability component. This allows us to simultaneously model both the stellar variability and eclipses (see §4).

3.1.2 Estimation of primary star properties from broadband colors

We estimated primary effective temperatures and masses using broadband color relations and absolute magnitudes presented in Table 1, respectively.

Effective temperatures ($T_{\rm eff}$) were estimated using the empirical color-$T_{\rm eff}$ relations presented in Mann et al. (2015) (their eq. 6) and David et al. (2016b) (their eq. 1, which is derived from fitting polynomials to the color and temperature data presented in Pecaut & Mamajek (2013) for dwarf stars, and is valid for $0.3<V-K_{s}<7.0$). These predict primary effective temperatures of $\sim$3250, 3190, 3240 and 3750 K for ADs 3814, 2615, 3116 and 1508, respectively. In §6.1 we directly determine the effective temperatures of both stars in each EB through modeling their spectral energy distributions (SEDs) and compare our $T_{\rm eff}$ values to these empirical predictions in Table 4.

We estimated primary masses from absolute K band magnitudes using the semi-empirical relation of Mann et al. (2015) (their eq. 10) and the empirical relation of Benedict et al. (2016) (their eq. 11). For this, we converted apparent to absolute magnitudes assuming a cluster distance of $182.8\pm 14$ pc (Gaia Collaboration et al., 2017) and assumed a reddening along the line of sight of $E(B-V)=0.027\pm 0.004$ (Taylor, 2006). These two relations predict primary masses of: $\sim$0.43, 0.34, 0.28 and 0.72 $M_{\odot}$ for ADs 3814, 2615, 3116 and 1508, respectively. For AD 1508 we used only the Mann et al. (2015) mass prediction as this system lies outside the validity range ($0.1\lesssim M\lesssim 0.6$ $M_{\odot}$) of the Benedict relation.

We note that these predictions are for single stars and hence are not appropriate for binary systems unless the system magnitudes are dominated by the primary component. Furthermore, these empirical relations are approximations only and are estimated from systems that typically do not contain as high a metallicity as Praesepe ([Fe/H] $\sim$ 0.1–0.27). Nonetheless, they serve to highlight the expected temperature and mass regimes of the systems to be analyzed.

3.2 Spectroscopy

We obtained high resolution spectra for each of the identified eclipsing binary systems using the Keck HIRES spectrograph (Vogt et al., 1994). The observations were taken between 2015 December and 2017 January, with the exact epochs along with estimated signal-to-noise ratios and measured radial velocities given in Table 2. The spectra cover the wavelength range $\approx$4800–9200 Å at a spectral resolution of $R>36,000$, and were reduced using the $makee$ software written by Tom Barlow. We measured radial velocities using the cross correlation techniques within the $fxcor$ task in $IRAF$, with absolute reference to between 3 and 5 (depending on the night) late type radial velocity standards. The standards and their approximate spectral types include: GJ 514 (M0.5), HD 95650 (M1), LHS 3433B (M2), Gl821 (M2), GJ 408 (M2.5), GJ 176 (M2.5), GJ 109 (M3.5), GJ 402 (M4), Gl 876 (M4), GJ 105B (M4.5), GJ 388 (M4.5), GJ411 (M4.5), GJ 406 (M6.5), with the reference velocities generally taken from Nidever et al. (2002). Telluric-free spectral regions were selected over between 6 and 19 orders (depending on the signal-to-noise of the target spectrum) for cross correlation function fitting. Depending on the velocity separation of the peaks, they were fit either singly or simultaneously, and depending on the signal-to-noise of the spectrum, the fitting function was either Gaussian or parabolic. Errors in the quoted radial velocities were determined from the empirical scatter among the measured orders and reference stars for each observation, with some hand editing to remove extreme outliers deriving from particularly poor measurements. In general, the scatter among the measurements that is quoted as the radial velocity error, is smaller than or comparable to the mean among the errors in the individual measurements over the orders and reference stars included in the quoted radial velocity value. This gives us some confidence that we are accurately representing the random errors in our methods.

AD 3814, AD 2615, and AD 1508 are detected as double-lined systems, with measurable radial velocities for each component at nearly all epochs. AD 3116, however, presented only a single line set, which we attribute to the primary. In the double-lined systems, the CCF peak height ratios were used to approximate the light ratio between the two components, which was then applied as a prior in the light curve modeling (see §4).

In addition to the radial velocities, H$\alpha$ equivalent width measurements were made for each EB using the $splot$ task in $IRAF$. The values quoted in Table 1 represent the combined system and the range records the variability over the various epochs of observation in Table 2.

4 Analysis with the GP–EBOP model

Both young and low-mass stars typically display photometric and spectroscopic modulation arising from the longitudinal inhomogeneity of active regions on the stellar surface, with activity timescales a strong function of stellar mass. In close binaries ($P\lesssim 15$ days), activity levels are generally observed to be higher than in their single star counterparts. This variability is important to properly account for when analyzing the observed stellar eclipses since it can subtly modify the detailed shape of individual eclipses. Ideally, therefore, we would model the stellar variability at the same time as fitting for the eclipses and, in doing so, propagate any uncertainties in the variability modeling through into the posterior distributions for the EB parameters. This approach motivated the development of a new eclipsing binary model, GP–EBOP, which we use here to characterize the new Praesepe EBs by simultaneously modeling the K2 light curves and Keck/HIRES radial velocity measurements, accounting for activity-induced effects. The method is distinct from those that account for stellar variability by detrending first and then modeling eclipses second.

4.1 GP–EBOP

GP–EBOP comprises a central eclipsing binary (EBOP) model coupled with a Gaussian process (GP) model, which has an MCMC (Markov chain Monte Carlo) wrapper. It can be used to model both eclipsing binary systems and transiting planets: we use it here in its first capacity but note its tested ability to model planet transits (e.g. Pepper et al., 2017). Below we briefly describe the main components of the model:

• •

EB component. The EB model is a modified version of the (JKT)EBOP family of models, which was first presented in Irwin et al. (2011). Each star is modeled as a sphere when computing light curves from the eclipses and as a biaxial spheroid for the calculation of reflection and ellipsoidal effects. This model is able to compute light ratios and radial velocities, and can correct for the “classical” light travel time across the system.

Differing from previous EBOP-based models, this implementation uses the analytic method of Mandel & Agol (2002) for the quadratic limb darkening law. GP–EBOP utilizes the LDtk toolkit (Parviainen & Aigrain, 2015), which allows uncertainties in the stellar parameters (effective temperature, surface gravity and metallicity) to be propagated through the PHOENIX stellar atmosphere models (Husser et al., 2013) and into priors on the limb darkening coefficients. Limb darkening parameterization within the fitting process follows the triangular sampling method of Kipping (2013).

• •

GP component. The GP model utilizes the george package222http://dan.iel.fm/george (Ambikasaran et al., 2014) and is used to model the out-of-eclipse (OOE) photometric data. A detailed description of Gaussian process regression is beyond the scope of this paper but the interested reader is referred to Roberts et al. (2012) for a gentle introduction, Rasmussen & Williams (2006) for a more detailed entry, Aigrain et al. (2012) for application to stellar light curves and Gillen et al. (2014) for application to eclipsing binary light curves and cross-correlation functions.

A simple way to view GPs is to think of them as a means of modeling a light curve by parameterizing the covariance between pairs of flux measurements, rather than explicitly specifying a functional form of model to fit the data. In this Gaussian process model, the joint distribution of the observed flux measurements is taken to be a multivariate Gaussian, whose covariance matrix is populated through a covariance function that depends on the observation times. As such, a GP is distribution over functions. When the parameters of a GP (called hyperparameters) are varied, we step through function space rather than the more familiar parameter space of conventional methods.

Crucially for our application, the power of GP regression is that we obtain an uncertainty on the prediction for the OOE variability across each eclipse, which we can then propagate through into our posterior distributions for the EB parameters.

• •

MCMC wrapper. GP–EBOP explores the posterior parameter space using the the Affine Invariant MCMC method, as implemented in emcee (Foreman-Mackey et al., 2013).

4.2 Light curves

The K2 light curves are a timeseries of flux measurements. GP–EBOP models the light curves by assuming the joint distribution over the flux measurements ${\bf F}$ is given by a multivariate Gaussian whose mean function $\mu$ is an eclipse model and whose covariance matrix ${\bf K}$ is described by a Gaussian process:

[TABLE]

The elements of the covariance matrix ${\bf K}$ are given by:

[TABLE]

where the first term represents the specific kernel chosen to describe to the out-of-eclipse variations and the second term describes the white noise component.

Figure 2 shows the raw light curves of the four new EBs and Figure 3 shows these phase-folded on the photometric variability period. The OOE light curves of all four systems presented here display evolving starspot modulation with characteristic amplitudes, periods and evolutionary timescales. To model these smoothly evolving data, therefore, we chose a GP with a quasi-periodic Exponential Sine Squared kernel (hereafter QPESS). This is a periodic kernel that is allowed to evolve over time, i.e. mimicking evolving starspot modulation. The QPESS kernel has the required flexibility to explain the large-scale flux variations in the OOE light curves. It is given by:

[TABLE]

The first exponential describes the periodic component and the second the evolution of the periodic signal. $A$ is the characteristic amplitude of the variations, $\Gamma$ is the scale of the correlations, $P$ the period of the oscillations and $l$ the evolutionary timescale. $t_{i}$ and $t_{j}$ represent example times of two flux measurements within the time series. The resulting periods (Table 6) differ from those reported by Rebull et al. (2017) (based on Lomb-Scargle techniques) at about the $\sim$1% level. The white noise term is given by:

[TABLE]

where $\sigma$ is the standard deviation and $\delta_{ij}$ is the Kronecker delta function. Within GP–EBOP the white noise term is incorporated via a multiplicative scale factor on the observational uncertainties, as george adds these scaled uncertainties in quadrature to the diagonal of the covariance matrix.

We model the k2sc light curves that have been detrended for instrument systematics but which still contain stellar activity variations. After visual inspection of the SAP and PDC k2sc light curves we opted to work with the PDC versions as these display lower point-to-point scatter.

As can be seen in Figure 2, numerous stellar flares are present throughout the light curves. Flares were treated in two ways depending on whether or not they affected the stellar eclipses. Those that did not were automatically removed using the following method: the light curve was smoothed using a running median filter, which was followed by running sigma cuts to identify flares. The data before and after the flare peak was removed until the light curve returned to the smoothed light curve value. Flares affecting the stellar eclipses were treated more carefully: as even a detailed modeling would not correct the photometry to a precision required to include in our eclipse modeling, we opted to conservatively mask out the affected data via visual inspection. The resulting light curves, which were modeled in our analysis, can be seen in Figures 7, 11 and 15 for ADs 3814, 2615 and 3116, respectively. The light curve of AD 1508 was treated slightly differently as only a preliminary solution is presented here (see §5.4 for details).

The full light curves (eclipses and out-of-eclipse variability) and radial velocity variations were simultaneously modeled by GP–EBOP  stepping through the parameter space 50,000 times with each of 144 ‘walkers’. The first 25,000 steps were discarded as burn-in and the remainder of each chain was thinned following inspection of the autocorrelation lengths for each parameter. To account for the $\sim$30 minute cadence of K2 observations, GP–EBOP was supersampled at 1 minute cadence and numerically integrated to the K2 sampling for model evaluation. The uncertainties on the limb darkening coefficients were inflated by a factor of 30, above the uncertainties derived from the PHOENIX models. This inflation factor was determined by comparing quadratic limb darkening coefficients of LDtk, Claret et al. (2012) and Sing (2010) for common $T_{\rm eff}$  $\log g$ and metallicity values in a representative range for our EBs across the Kepler bandpass. We used the spread in their predictions, and applied a further increase to account for systematic uncertainties in M-dwarf model atmospheres, to determine our inflation factor. Reflection effects and gravity brightening were not included in the modeling. The former is accounted for by the GP model and the latter makes no significant difference to the model posterior distributions, which we tested by performing model runs with different gravity brightening exponent ($\beta$) values. We note that Alencar & Vaz (1997) found that $\beta$ ranges between 0.2 and 0.4 for stars with temperatures between $3700\leqslant T\leqslant 7000$ K and that the typical Lucy (1967) value of $\beta=0.32$ best describes stars with $T=6500$ K.

4.3 Radial velocities

The Keck/HIRES RVs were modeled using Keplerian orbits simultaneously with the K2 light curves. Spectroscopic light ratios (available for three of the four systems presented here) were estimated from cross-correlation peak heights and applied as priors on the light curve model component. This can help break the well-known degeneracy between the radius and surface brightness ratios, which can often be a limiting factor in the individual radius estimates for near equal-mass EBs.

An RV jitter term, incorporated in GP–EBOP, was used to allow the uncertainties on the Keck/HIRES RV measurements to be scaled, if necessary. This helps account for additional variations arising from e.g. stellar activity and instrument systematics. This jitter term is added in quadrature to the observational uncertainties. When RVs from multiple instruments are obtained, GP–EBOP can scale the uncertainties for each instrument individually and account for offsets between different instrument RV zero points.

5 Results

The K2 light curves and Keck/HIRES radial velocity measurements of the four new EBs (ADs 3814, 2615, 3116 and 1508) were modeled with GP–EBOP; the results for each system are discussed in turn below. Throughout our analysis we define the primary as the star that, when occulted, gives the deepest eclipse, and the secondary as the occulting star. We note that these adjectives do not necessarily imply that the primary star is the more massive or brighter star, as we find to be the case with AD 2615.

5.1 AD 3814

AD 3814 has been extensively studied in the literature. The M3.4 spectral type, broadband photometric magnitudes and colors, and proper motion give AD 3814 a high probability of cluster membership. Figure 7 shows the K2 light curve used in the modeling after flares were removed. Three eclipses were masked in the flare removal process (see section 4.2): two secondary eclipses at rBJD333rBJD = BJD $-$ 2454833.$\sim$2315 and 2361, and one primary eclipse at rBJD$\sim$2364. The red line and pink shaded region indicate the mean and 2$\sigma$ uncertainty of the posterior GP–EBOP eclipse model, which is able to reproduce both the eclipses and the slowly evolving starspot modulation.

Detrending with respect to the GP component and phase-folding on the binary orbital period allows us to inspect the shape of the eclipses in detail. These are shown in Figure 8, where the top panel displays the full phase-folded light curve and the bottom panels show zooms around primary and secondary eclipses (left and right, respectively). There is clear evidence of increased scatter in the residuals across each eclipse, which is presumably due to uncorrected differential starspot effects. Starspots on the background star will have a differential effect on the eclipse shape, with the eclipse being shallower if starspots on the background star are preferentially occulted by the foreground star and deeper if the unspotted photosphere is preferentially occulted. As the timescale for such differential effects are much faster than the typical starspot modulation observed out of eclipse, the QPESS kernel will struggle to account for this effect given its covariance properties, which constrain it to smooth variations. Instead, the GP will opt to inflate its uncertainty due to the increased scatter. One could theoretically include an additional kernel within the GP model to try and account for such differential effects across eclipses, but this is beyond the scope of the current work.

The 8 Keck/HIRES RVs were modeled simultaneously with the K2 light curve. The resulting phase-folded RV orbit is shown in Figure 9 (primary in red and secondary in blue). The colored lines and shaded regions indicate the median and 2$\sigma$ uncertainties on the posterior orbits of the two stars, which are well-fit to the observed RVs. The systemic velocity of the system is $V_{\rm sys}=33.60\pm 0.24$ km s*-1* (dashed gray line), which is consistent with the recessional velocity of the cluster $V_{\rm rec}$ $\sim$33–35 km s*-1* (e.g. van Leeuwen, 2009; Quinn et al., 2012; Yang et al., 2015) and hence provides further evidence of cluster membership. We note that the residuals in the phase-folded RV plot display an interesting structure. Inspection of the RV residuals in time, however, does not suggest any long term trend indicative of a tertiary companion, which is consistent with the lack of a detectable tertiary peak in the cross-correlation function. Possible explanations for the residuals are issues with the absolute radial velocity calibration, the RV stability of the reference standards, or the precise placement of the target star in the center of the slit. GP–EBOP attempts to account for this unknown noise component by including an additional jitter term that acts to scale the observational uncertainties. We note that if the origin of this noise component were known, it may be possible to model directly within the fit, but this is beyond the scope of the present analysis.

Figure 10 depicts the system, to scale, at both primary and secondary eclipse, indicating the geometry responsible for the observed eclipses and RV variations. The model parameters and 1 sigma uncertainties for AD 3814 are presented in the first results column of Table 6. The light curve and RV modeling with GP–EBOP yields masses and radii for each star: the primary and secondary masses are $0.3813\pm 0.0074$ and $0.2022\pm 0.0045$ M*⊙* with corresponding radii of $0.3610\pm 0.0033$ and $0.2256\,^{+0.0063}_{-0.0049}$ R*⊙*. The masses of both components are constrained to 2% and the primary and secondary radii to 1% and 3%, respectively. The fundamental parameters are compatible with the estimated M3.4$\pm$0.1 spectral type and the primary mass estimate from section 3.1.2. The masses, radii and effective temperatures (derived in §6.1) of AD 3814 are compared to the current suite of stellar evolution models in section 6.2.

We applied a prior on the system light ratio and priors on the quadratic limb darkening coefficients (see Table 3). The light ratio was determined from the cross-correlation peak height ratio in a HIRES spectrum taken close to quadrature, which is acceptable as the HIRES spectral range is a reasonable match to the K2 bandpass. We note that the degeneracy between the surface brightness and radius ratios is not apparent in our posteriors, although it is not expected to be significant in this system given the mass and brightness ratios.

We conclude by noting that this system would benefit from a more detailed modeling of the individual eclipses, incorporating a full starspot model, to assess whether the large-scale underlying starspot distribution can be reconstructed from the eclipses, which track different longitudes on the stellar surfaces over the K2 run.

5.2 AD 2615

AD 2615 is an M4.0 high probability member of Praesepe. The analysis presented here is consistent with the photometric, spectroscopic and membership information from previous studies. The light curve of AD2615 that was used in the modeling is shown in Figure 11. One secondary eclipse, at rBJD$\sim$2367, was masked following the flare removal process (see section 4.2). The red line and pink shaded region represent the mean and 2$\sigma$ uncertainty of the posterior GP–EBOP eclipse model. As with AD 3814, the model is able to capture both the stellar eclipses and the evolving starspot modulation. The model’s predictive power can be seen before and after the light curve, where it is able to predict the form of the evolving modulation pattern, given the covariance properties of the data; this also drives the motivated prediction and uncertainty across each eclipse.

Figure 12 shows the phase-folded light curve, which has been detrended with respect to the GP component. The eclipse model is an acceptable fit to the data. There is no clear evidence for increased scatter in the residuals, which suggests that the geometry of the eclipses does not preferentially track bright or dark regions on the stellar surfaces, perhaps because the underlying starspot distribution in AD 2615 is more homogeneous than in AD 3814.

Figure 13 shows the phase-folded RV orbit (red for primary and blue for secondary). The 5 HIRES RVs of both stars are well-fit by the Keplerian model. The 2$\sigma$ uncertainties on the orbits (red and blue shaded regions) increase around quadrature, as expected. The systemic velocity of $V_{\rm sys}=34.91\pm 0.39$ km s*-1* (dashed gray line) is compatible with the cluster’s recessional velocity, providing further kinematic evidence of cluster membership. We note that a sixth RV observation was conducted but lay too close to primary eclipse to disentangle the two stellar components and hence was not used in the fit. In principle, we could determine an upper limit on the separation of the two stars and use this as an additional constraint in the modeling. However, at phase = 0.997, the solution is already tightly constrained and hence this upper limit would not place useful constraints on our existing solution. We further note that spectral disentangling may offer an interesting alternative route of RV determination for this system, which could utilize this sixth observation. While traditional spectral disentangling techniques require many high SNR spectra, powerful new techniques are emerging designed for fewer and lower SNR spectra (e.g. Czekala et al., 2017). It would be interesting to compare the standard CCF-based RV determination with these new spectral disentangling techniques, but this is beyond the scope of the present paper.

Figure 14 depicts the system, to scale, at primary and secondary eclipse, showing the configuration responsible for the observed eclipses. The medians and 1$\sigma$ uncertainties of the GP–EBOP model posteriors are reported in Table 6 (second results column). The primary and secondary masses are $0.212\pm 0.012$ and $0.255\pm 0.013$ M*⊙, with corresponding radii of $0.233\pm 0.013$ and $0.267\pm 0.014$ R⊙*. We remind the reader that we define the primary star as the star which, when occulted, gives the deeper eclipse, but that this does not necessarily mean it is the more massive or brighter of the two components, as indeed is the case in this system. The masses and radii are constrained to 6% for the primary and 5% for the secondary. This system would benefit from additional RVs around quadrature to increase the precision of the mass determination. The fundamental parameters are compatible with the estimated M4.0 spectral type but the mass of either component is lower than the estimate from the system’s absolute K-band magnitude (Section 3.1.2), presumably because this system is a near equal mass binary and so both stars contribute significantly to the K-band flux, resulting in an overestimated single-star mass. The masses, radii, and effective temperatures (c.f. §6.1) are compared to stellar evolution models in section 6.2.

We applied priors on the system light ratio and stellar limb darkening coefficients (see Table 3). Even though the system is near-equal mass and brightness, our spectroscopic light ratio was able to break the degeneracy between the surface brightness and radius ratios, which can be a limiting factor in determining radii in such systems.

5.3 AD 3116

AD 3116 is an M3.9 high probability member of Praesepe. The system sits at the bottom of the cluster sequence (see Figure 1) suggesting the secondary component contributes little optical light to the total system flux and hence is comparatively low-mass.

Analysis of the K2 light curve and 7 HIRES spectra reveals the system to be single-lined with eclipses visible only on the primary component, consistent with its position in color-magnitude space. Secondary spectroscopic lines could not be detected, even after dividing the two spectra and looking for similar but weaker patterns in the CCF, which suggests the secondary contributes very little ($<20-35$%) to the system’s optical light. Given the lack of a detectable secondary eclipse and secondary RV orbit, the data alone are not able to constrain the solution precisely. There exist two families of solutions: one consisting of a small secondary that fully transits and the other a larger secondary on a grazing trajectory. The primary RV orbit requires the secondary to be eclipsed, and hence both models find a negligible surface brightness ratio in the Kepler band to remain consistent with the lack of a detectable secondary eclipse. For the solution comprising a large ($R_{\rm sec}/R_{\rm pri}\gtrsim 1$), grazing secondary, this would require an unusual object possessing a very low temperature given its radius. Inspection of the system mass function revealed that the secondary lay in the brown dwarf regime ($M_{\rm sec}\sim$55 $M_{\rm Jup}$), which further supported the solution comprising a small, fully-transiting secondary. We tested the reliability of the primary RV solution to see if individual RVs close to the systemic velocity (i.e. which could be biased by low-level secondary light) may be affecting the eccentricity of the RV orbit and hence the system parameters. We removed all bar the three RVs closest to quadrature and, as expected, the model converged again on a solution requiring the secondary to be eclipsed. This, combined with the small primary RV semi-amplitude and lack of secondary spectroscopic lines, rules out a scenario where the secondary is of comparable size and brightness to the primary but there is no secondary eclipse due to the eccentricity of the orbit. All available information and tests pointed towards a very low-mass, small and cool secondary component.

We therefore chose to place loose uniform priors on the radius ratio and surface brightness ratio to encourage the solution towards a physically sensible secondary component. These priors were $0.0\leqslant R_{\rm sec}/R_{\rm pri}\leqslant 0.46$ and $0.0\leqslant J_{\rm K2}\leqslant 0.25$ which, given the expected primary star properties (c.f. §3.1.2) and secondary star mass estimate, act simply to exclude physically implausible solutions and do not act to constrain the remaining physically plausible solutions. We performed further tests allowing $R_{\rm sec}/R_{\rm pri}$ and $J_{\rm K2}$ to extend up to 0.75 and 0.35, respectively, but find consistent posterior values.

The model fit is shown in Figures 15–18, whose descriptions are the same as for ADs 3814 and 2615 in sections 5.1 and 5.2 above. The model is a good fit to the primary eclipse and large-scale evolving starspot structure in the K2 light curve (Figure 15). We note that two primary eclipses, at rBJD$\sim$2319 and 2350, were masked following the flare removal process (see section 4.2). Figure 16 shows the phase-folded and GP-detrended light curve: the primary eclipse is well-fit, although there is a modest increase in the residual scatter, which is larger than in AD 2615 but smaller than in AD 3814. The RV data suggests a moderately eccentric orbit ($e\sim 0.15$; see Figure 17) with a systemic velocity of $V_{\rm sys}=34.93\,^{+0.61}_{-0.53}$ km s*-1* (dashed gray line). This is consistent with the cluster’s recessional velocity, providing additional kinematic evidence of cluster membership.

Using the empirical relations of Benedict et al. (2016), and assuming the van Leeuwen (2009) cluster distance of $181.5\pm 6.0$ pc, the $K_{\rm s}$ magnitude of AD 3116 implies a primary mass of $M_{\rm pri}=0.28\pm 0.02$ $M_{\odot}$, where the uncertainty arises equally from the empirical relation scatter and our assumed 0.1 mag uncertainty on the quoted $K_{\rm s}$ value. We checked this value using the empirical relations of Mann et al. (2015) and find $M_{\rm pri}\sim 0.29$ $M_{\odot}$, consistent with the Benedict et al. value. Taking the Benedict value, the mass function from our final solution then yields $M_{\rm sec}=54.2\pm 4.3$ $M_{\rm Jup}$. This is one of only $\sim$20 known transiting brown dwarfs (e.g. Csizmadia, 2016; Nowak et al., 2016; Bayliss et al., 2017) and the primary component is one of only three M-dwarfs known to host a transiting brown dwarf. Furthermore, this is only the second known transiting brown dwarf in an open cluster (i.e. where the age is well-constrained), and the first younger than a Gyr.

Figure 18 shows the system geometry at primary and secondary eclipses. That the brown dwarf is fully occulted yet shows no detected signature in the K2 band theoretically allows us to place an upper limit on the optical reflected light and hence albedo of the object. Using the Mann et al. (2015) empirical relations to estimate the primary radius, and hence secondary radius and semi-major axis from our light curve modeling, we can estimate the system scale. This then allows us to compute the angle on the sky that the brown dwarf subtends as seen from the primary. With $R_{\rm sec}\sim 0.11$ $R_{\odot}$ and $a\sim 4.7$ $R_{\odot}$, the brown dwarf intercepts $\sim$0.007% of the visible light from the primary star. Therefore, even if the brown dwarf reflected all incident flux (i.e. an albedo of 1), we would not detect a drop in flux in the K2 light curve when the brown dwarf is occulted.

We applied priors on the limb darkening coefficients (see Table 3). The secondary temperature was set to be as low as the PHOENIX models allow but is likely still too high (see Table 4). However, as the secondary gives no detectable eclipse it makes no significant difference to the presented solution. Given the system is single-lined we did not place a prior on the system light ratio in the K2 band.

5.4 AD 1508

AD 1508 is a high probability M0.1 member of Praesepe, which sits high above the cluster sequence (see Figure 1), suggesting a near-equal mass system. The preliminary analysis presented here is consistent with this picture. The K2 light curve of AD 1508 (see Figures 2 and 3; bottom panels) is dominated by evolving starspot modulation at the few percent level. Very shallow grazing eclipses are also present with a depth of less than 1%. We obtained only three RVs for this system, which unfortunately fall close to primary and secondary eclipses (see Table 2). Given this, and the shallow eclipses, a precise solution is not possible. Instead, we provide our initial analysis and offer the system to the community for further pursuit.

The K2 light curve and three Keck/HIRES RVs were simultaneously modeled with GP–EBOP. However, given the preliminary nature of the modeling, and unlike the other three systems, we opted to simplify the light curve analysis by performing an initial detrending of the starspot modulation and then modeled the residuals with GP–EBOP to analyze the stellar eclipses. To do this, the out-of-eclipse light curve was flattened through two iterations of a cubic basis spline with knots every 2 hours and rejection of 0.5-$\sigma$ outliers. Figure 19 shows the resulting detrended light curve that was modeled with GP–EBOP. Low-level (likely systematic) residual variations are present, which show a relatively rough behavior. Accordingly, we chose a Matern-3/2 kernel for the GP component, which is given by:

[TABLE]

where $A$ is the amplitude and $l$ the characteristic timescale of the variations.

Detrending with respect to the GP component and phase-folding on the orbital period, as shown in Figure 20, we see that the eclipses are well-fit by the model. There is no significant evidence of increased scatter across the eclipses. We note that the light curve of AD 1508 appears noisy in comparison to the other systems discussed here, even though it is significantly brighter. This is simply because the plot scales in Figures 19 and 20 are small as the eclipses are shallow and the starspot modulation has already been detrended for. It is not a reflection of the true noise level in this system: the point-to-point scatter of all systems discussed here decreases with system brightness, as expected.

The phase-folded RV orbit is shown in Figure 21 which, given only three RVs at non-optimal phases, is not well-constrained. This is reflected in the large 2$\sigma$ uncertainties on the posterior orbits (red and blue for the primary and secondary stars, respectively). Nonetheless, the systemic velocity is relatively well-constrained at $V_{\rm sys}=33.1\pm 1.7$ km s*-1*, which is consistent with the cluster recessional velocity and hence provides further kinematic evidence of Praesepe membership. Figure 22 shows the system, to scale, at primary and secondary eclipse. The shallow eclipses simply result from the very grazing trajectory of the stellar orbits, as viewed from Earth.

The median and 1$\sigma$ uncertainties resulting from our preliminary analysis are reported in Table 6 (fourth results column). Given the available data, significant uncertainties exist in the derived masses and radii. The primary and secondary masses are $0.45\,^{+0.19}_{-0.14}$ and $0.53\,^{+0.22}_{-0.16}$ M*⊙* with corresponding radii of $0.549\,^{+0.099}_{-0.082}$ and $0.454\,^{+0.094}_{-0.101}$ R*⊙*. The solution is currently limited by the lack of RV constraints and future analysis would benefit from additional RV measurements, especially around quadrature. Nonetheless, the fundamental parameters are compatible with the estimated M0.1$\pm$0.1 spectral type and the primary mass estimate from section 3.1.2. Given the existing uncertainties we do not compare this system to stellar evolution models in section 6.2.

We applied priors on the system light ratio and limb darkening coefficients (see Table 3). Although large uncertainties remain, the spectroscopic light ratio was able to break the degeneracy between the surface brightness and radius ratios, which can be a limiting factor in determining individual radii in near-equal mass and brightness systems.

6 Discussion

The direct determination of fundamental stellar parameters offers an opportunity to test stellar evolution models. The fundamental predictions of these models are the radius and $T_{\rm eff}$ for a star of given mass and metallicity as a function of age. Ideally, therefore, we would be able to determine the mass, radius and $T_{\rm eff}$ of both stars as, together, these offer a particularly strong test of stellar evolution theory. However, while the masses and radii of stars in EBs naturally fall out of the joint light curve and radial velocity modeling, estimating effective temperatures is more challenging. In §6.1 we present a method of simultaneously estimating the effective temperature of both stars, and the distance to the system in a manner that makes full and correct use of the light and radial velocity constraints. We then compare our $T_{\rm eff}$’s and distances to empirical $T_{\rm eff}$ relations and to previous distance estimates to Praesepe. In §6.2 we compare our masses, radii and $T_{\rm eff}$’s to the predictions of stellar evolution models for individual systems and also place the newly characterized EBs in the context of other known low mass EBs and briefly discuss the constraints that can be placed on the age of Praesepe. Through this model comparison, and in §6.3 where we comment on the synchronization of the new EBs, we discuss several astrophysical implications of our findings.

6.1 Simultaneous determination of effective temperatures and distance from the spectral energy distribution

The standard method of estimating $T_{\rm eff}$ is the following: 1. estimate the primary star $T_{\rm eff}$ from either system colors adopting empirical single-star relations, or use (typically) low resolution spectra to infer a combined spectral type (SpT) and convert this into a primary star $T_{\rm eff}$. 2. estimate the secondary star $T_{\rm eff}$ from the primary $T_{\rm eff}$ and the light curve surface brightness (and hence temperature) ratio. There are a number of issues with this approach: empirical color-$T_{\rm eff}$ and SpT-$T_{\rm eff}$ relations for single stars are not necessarily applicable for all binary systems and the temperature ratio estimated from the light curve is specific to that band, i.e. $(T_{\rm sec}/T_{\rm pri})_{\rm\,band}$; it is not a $T_{\rm eff}$ ratio.

A more direct approach would be to model the system’s spectra, but to do so would require high SNR (signal-to-noise) data, which would normally require the co-adding of spectra. While feasible for single star systems, this is not possible for binaries as there are two varying components. One approach would be to disentangle the spectra into their individual components and model these directly to estimate $T_{\rm eff}$ of each star (e.g. Czekala et al., 2015, 2017). However, while powerful, this approach is both time and computationally intensive, and the distance to the system remains unknown (unless the spectra are also flux calibrated).

A method of simultaneously determining $T_{\rm eff}$ of both stars, and the distance to the system, is to model the system’s spectral energy distribution (SED). This approach is not computationally intensive, does not rely on empirical single-star relations and readily incorporates priors from the joint light curve and RV modeling. Importantly, with respect to the last point, it correctly interprets the band-specific surface brightness ratio from the light curve modeling. Therefore, we simultaneously estimate $T_{\rm eff}$’s and the distance to ADs 3814, 2615, 3116 and 1508 using the following method:

SEDs were constructed using broadband magnitudes readily available in the literature. We obtained SDSS $ugriz$ magnitudes from the Sloan Digital Sky Survey Data Release 13, and 2MASS JHKs and WISE data from the NASA/IPAC Infrared Science Archive. These are reported in Table 1 along with their formal measurement uncertainties. 2. 2.

Model grids of both BT-SETTL (Allard et al., 2012) and PHOENIX v2 model spectra (Husser et al., 2013) were convolved with commonly available bandpasses ($ugriz$, UBVRI, 2MASS JHKs, Spitzer/IRAC, WISE and Kepler) to create a model grid of bandpass fluxes. 3. 3.

Each SED was modeled by interpolating the model grids in $T_{\rm eff}$–$\log g$ space. We opted to fix the metallicity at Z=0.0 given the cluster [Fe/H] value but note it is possible to include in the interpolation. 4. 4.

The parameters of the fit were the $T_{\rm eff}$, radius and $\log g$ of each star, the distance to the system, the interstellar extinction and the uncertainty scale factor ($T_{\rm pri}$, $T_{\rm sec}$, $R_{\rm pri}$, $R_{\rm sec}$, $\log g_{\rm pri}$, $\log g_{\rm sec}$, $d_{\rm sys}$, $A_{\rm v}$ and $\sigma_{s}$). The radii and $\log g$’s have priors from the joint light curve and RV solution, $A_{\rm v}$ had a prior determined for the cluster (Taylor, 2006), and the temperatures, distance and uncertainty scale factor had uninformative priors. The uncertainties on the magnitudes were initially set by adding the observed variability level to the formal measurement errors in quadrature and a further inflation term ($\sigma_{s}$) was then fit for. 5. 5.

The posterior parameter space was explored using emcee with 50,000 steps and 196 ‘walkers’. Convergence was assessed using the Gelman–Rubin diagnostic plus examination of individual sections of the chains. A conservative burn-in was estimated comprising the first 25,000 steps for all systems and parameter distributions were derived from the remainder after thinning each chain based on the autocorrelation lengths of each parameter. 6. 6.

This method also gives the option of placing additional priors in the modeling. For example, one can place a prior, from the light curve modeling, on the surface brightness ratio between the two stars in the band observed, rather than incorrectly placing a $T_{\rm eff}$ ratio constraint. In the case of single-lined systems, radius ratio constraints and surface brightness upper limits can also be placed.

Both BT-SETTL and PHOENIX v2 model spectra are able to reproduce the broadband magnitudes of ADs 3814, 2615, 3116 and 1508. We note, however, that the BT-SETTL models consistently underpredict the optical $r$ band fluxes, whereas the PHOENIX v2 models predict higher red-optical fluxes in agreement with the data for all sources. Accordingly, in Figure 4 we show the PHOENIX v2 model fits to the observed broadband magnitudes of ADs 3814, 2615 and 1508 reported in Table 1 (for AD 3116 we show the BT-SETTL fit as the PHOENIX models do not extend to low enough temperatures to explain the secondary brown dwarf component). The $T_{\rm eff}$ and distance values derived from our SED fitting procedure with both the BT-SETTL and PHOENIX v2 models are reported in Table 4 along with the empirical relation predictions of Mann et al. (2015) and David et al. (2016a). We discuss the effective temperature and distance estimates in the following two sections.

6.1.1 Effective temperatures

We find that the BT-SETTL model temperatures are typically $\sim$40 K hotter than the PHOENIX v2 values, although both sets of temperatures agree to within 1$\sigma$. They are also both in agreement with the temperatures predicted by empirical relations. We note that both sets of empirical relations used the BT-SETTL models to calibrate their temperature scale and hence caution should be applied when interpreting the slightly closer agreement between the empirical relations and BT-SETTL SED temperatures than with the PHOENIX v2 values.

Given the slight offset between the BT-SETTL and PHOENIX temperatures we opted to combine the two predictions for each star as our final $T_{\rm eff}$ values. These are reported in Table 4 as the “combined” model and are: $T_{\rm pri}=3211\,^{+54}_{-36}$ K and $T_{\rm sec}=3103\,^{+53}_{-39}$ K for AD 3814; $T_{\rm pri}=3152\,^{+57}_{-40}$ K and $T_{\rm sec}=3131\,^{+56}_{-38}$ K for AD 2615; and $T_{\rm pri}=3767\,^{+99}_{-85}$ K and $T_{\rm sec}=3693\,^{+122}_{-135}$ K for AD 1508. For AD 3116 we used only the BT-SETTL models given the expected temperature of the brown dwarf secondary.

While both SED modeling and empirical relations yield consistent results, the SED modeling constraints are significantly tighter (even combining both sets of results), which is perhaps unsurprising given they are system-specific and capitalize on the joint light curve and RV modeling constraints. Furthermore, interpreting the temperature ratio from the light curve modeling as a genuine $T_{\rm eff}$ ratio is incorrect in all cases where the bandpass observed does not cover the majority of the integrated spectra of both EB components, and the system is not equal mass. For both ADs 3814 and 2615, using the Kepler bandpass temperature ratio as a $T_{\rm eff}$ ratio (as required when using empirical relations) results in a steeper temperature scale than the light curve modeling results actually imply, i.e. the secondary is predicted to be cooler than expected relative to the primary temperature. This effect is most noticeable in AD 3814 given the larger mass ratio in this system.

6.1.2 Distance to Praesepe

Literature distance estimates to Praesepe range from $\sim$160–190 pc with the more recent determinations clustering around $\sim$175–185 pc (Mermilliod et al., 1990; Reglero & Fabregat, 1991; Gatewood & de Jonge, 1994; Percival et al., 2003; An et al., 2007; van Leeuwen, 2009; Gaia Collaboration et al., 2017). Gaia DR1 parallaxes imply a distance of $182.8\pm 1.7\pm 14$ (the two uncertainties are the error on the cluster center determination and the observed spread of cluster members on the sky; Gaia Collaboration et al., 2017). Our distance estimates for ADs 3814, 2615 and 1508 are $170.4\,^{+11.0}_{-8.9}$, $179\,^{+13}_{-10}$ and $166\,^{+25}_{-23}$ pc, respectively, which are all in agreement with the Gaia parallax distance. As AD 3116 is single-lined, we do not have precise radii and surface gravities, so we placed a prior on the distance to the system of $d_{\rm sys}=182.8\pm 14$ pc, and hence do not quote a distance for this system as we essentially recover our prior.

Empirical bolometric corrections (BCs) are available for M-dwarfs (e.g. Mann et al., 2015). Combining these with our calculated radii gives the system bolometric flux, which can be converted to absolute bandpass magnitudes using the derived BCs and compared to apparent magnitudes to estimate the distance using the distance modulus (see M15 distances in Table 4). We note that these are also in agreement with both our distances and the Gaia cluster value.

6.2 Comparison with stellar evolution models

6.2.1 The Newly characterized EBs

With precise masses, radii and effective temperatures for both stars in ADs 3814 and 2615 we can test the predictions of stellar evolution theory for low mass stars at the beginning of the main sequence phase of evolution. Figure 5 compares the fundamental parameters of ADs 3814 and 2615 to the PARSEC v1.2 (Bressan et al., 2012; Chen et al., 2014) and BHAC15 (Baraffe et al., 2015) models. Praesepe is slightly metal-rich ([Fe/H]$\approx$0.14) but the closest BHAC15 models in metallicity are solar composition. Therefore, we compare our results with both the solar metallicity PARSEC and BHAC15 models (Figure 5, top row) and also compare to the PARSEC models at Praesepe metallicity (Figure 5, bottom row). In the mass-radius plane (left panels) the PARSEC models (solid lines) predict slightly larger radii than BHAC15 (dashed lines) for a given mass, but both models are able to explain the two components of each system with a single isochrone at the 1$\sigma$ level (for PARSEC this is true for both solar and Praesepe metallicities). This agreement is encouraging as the masses of AD 3814 are constrained to 2% for both components and the primary and secondary radii to 1% and 3%, respectively. The uncertainties on the masses and radii of AD 2615 are slightly larger, given the system is fainter, and there are fewer eclipses and RVs, but the masses and radii are still both constrained to 6% for the primary and 5% for the secondary.

We note that both systems are young (sub-Gyr) and display modest H$\alpha$ emission. Therefore, compared to old M dwarfs these Praesepe stars are expected to have relatively strong magnetic fields and high spot coverage. Higher activity levels are thought to result in stars with lower effective temperatures and inflated radii (e.g. Chabrier et al., 2007; MacDonald & Mullan, 2014), and this is often seen in observations (e.g. Feiden & Chaboyer, 2012). Stars in EBs with longer orbital periods appear to show better agreement with the models, but those that do show disagreement tend to be fully convective. This might suggest that for stars with radiative cores and convective outer envelopes, disagreements with models are driven by rotation and magnetic activity but comparisons for fully convective stars are subject to other errors (Feiden, 2015). That these two fully (or almost fully) convective EB systems are active and have relatively short (6–12 day) periods yet agree well with the radius predictions of non-magnetic models presents a further challenge to stellar evolution theory.

While the masses and radii appear to be in agreement, including $T_{\rm eff}$ complicates the picture. We next compare our results in the $T_{\rm eff}$–$\log g$ plane. The surface gravity, $\log g$, combines the mass and radius information, which agree well for both models, and hence this parameter should also be well explained. In the middle column of Figure 5 we see significant discrepancies between the data and models, which points towards problems in the model $T_{\rm eff}$ scales. The models substantially diverge in their $T_{\rm eff}$ predictions, with the BHAC15 models being hotter by $\sim$200–250 K across the mid M-dwarf range, and the PARSEC models being perhaps 10–25K cooler than the data. We note that this is also seen in the mass–$T_{\rm eff}$ and radius-$T_{\rm eff}$ planes (not shown here). Both sets of models essentially predict the same $T_{\rm eff}$ independent of age for $t\gtrsim$400 Myr out to 10 Gyr. Our SED analysis yields $T_{\rm eff}$ values that are in closer agreement to the PARSEC models than BHAC15, but both models predict a steeper $T_{\rm eff}$ scale than the data suggest (note that a steeper model $T_{\rm eff}$ scale manifests as a shallower gradient in $T_{\rm eff}$-$\log g$ space, as observed). One option is that the model $T_{\rm eff}$ scales are too steep for mid M-dwarfs but it could also be that additional phenomena, not included in the models, are responsible for the observed slope difference. Both ADs 3814 and 2615 display starpot modulation in the K2 light curves. As neither PARSEC nor BHAC15 include the effect of magnetic fields and starspots it could be that some of the discrepancy arises from these phenomena rather than the model $T_{\rm eff}$ scale being too steep per se.

Although the primary component of AD 3814 agrees with the PARSEC Praesepe metallicity models, the secondary lies above the relation. We can take the primary star as an example to explore the required spot coverage and contrast ratio needed to bring its computed $T_{\rm eff}$ onto the same expected isochrone as the secondary component. We note that this scenario would require the PARSEC $T_{\rm eff}$ scale to be underpredicting the true unspotted $T_{\rm eff}$ but this is plausible so we continue with the exercise nonetheless. Assuming a spot-to-unspotted photospheric temperature ratio of 0.8 (e.g. Grankin, 1998) would require $\sim$25% spot coverage. To bring the primary and secondary components within 1$\sigma$ would only require a 10% spot coverage on the primary. We note, however, that the radius posterior medians sit just below the zero age main sequence predicted by the PARSEC models and invoking starspots to redress the $T_{\rm eff}$ slope differences would imply a corresponding decrease in the radii for these stars without spots.

To bring the primary and secondary components of AD 3814 into agreement with the BHAC15 models would require spot coverages of 30–40% on each star. While high, this is consistent with observations of active late-type stars, especially those in close binaries (e.g. O’Neal et al., 2004). We note that the BHAC15 models track a steeper path in $\log g$-$T_{\rm eff}$ space beyond 3400 K (corresponding to a shallower $T_{\rm eff}$ scale). Simply shifting the BHAC15 models cooler by 250 K would bring them into agreement with all four stars. This is not possible with the PARSEC models, so it remains a valid option that the PARSEC model temperature scale is too steep over the mass range probed ($\sim$0.2–0.4 $M_{\odot}$). However, more precisely characterized M-dwarf binaries are required to confirm this tentative statement.

The radii and effective temperatures combine to determine the luminosity of a star. Stellar evolution models are typically found to underpredict the radii and overpredict the effective temperatures of active low-mass stars; however, these combine to essentially recover the correct luminosity. The right column of Figure 5 shows the radiative $T_{\rm eff}$-luminosity relation. As expected, the BHAC15 models appear to underpredict the luminosity because the model $T_{\rm eff}$ is too high. The PARSEC models are in better agreement: they are able to follow the general trend of the data and explain the primary component of AD 3814 and the secondary of AD 2615, but the other components are slightly discrepant at the $\sim$1.5$\sigma$ level.

6.2.2 Updated mass–radius relation for low-mass EBs

Figure 6 shows the mass-radius relation for detached double-lined eclipsing binaries below $M<1.5$ $M_{\odot}$. Field EBs are shown in gray while members of young open clusters – including our newly discovered systems reported here – are colored by cluster (see figure caption for color scheme). The fundamental parameters of the known cluster EBs with both components below 1.5 $M_{\odot}$ are reported in Table 5. The three double-lined Praesepe systems reported here make a significant contribution to known cluster EBs, increasing the total number below 1.5 $M_{\odot}$ by almost 20% (and increasing the known double-lined M-dwarf EB population by 30%). Furthermore, ADs 3814 and 2615 add precise constraints for stellar evolution models at the zero-to-early age main sequence for low-mass stars.

6.2.3 Age of Praesepe

As briefly discussed in the introduction, the age of Praesepe has been debated in recent years. It has typically been estimated at $\sim$600–650 Myr by isochrone fitting, often through association with the Hyades (e.g. Perryman et al., 1998; Salaris et al., 2004; Fossati et al., 2008). However, Brandt & Huang (2015b) found that including rotation in stellar models implied an age of $790\pm 60$ Myr (2$\sigma$ uncertainty) for Praesepe, which is in agreement with their Hyades age of $\sim$750–800 Myr (Brandt & Huang, 2015a). This older age estimate arises from the fact that rotation results in longer main sequence lifetimes and hence older ages for post-turnoff populations. This result was corroborated by David & Hillenbrand (2015), who also include the effect of stellar rotation in their comparison between stellar atmospheric parameters (derived from Strömgren photometry) and theoretical isochrones.

Somewhat orthogonal to the ages inferred from radiative properties such as $L$ and $T$, the ages of EB systems can be determined through comparison of their masses and radii with stellar evolution models (see section 6.2.1). Unfortunately, over the mass range probed by our EBs, the several hundred Myr Praesepe sits roughly at the zero age main sequence. As M-dwarf evolution is slow, their increase in radius as they evolve through their first several Gyr on the main sequence is correspondingly small. Therefore, using our masses and radii to independently estimate the age of Praesepe would carry significant uncertainty and would not provide useful input to the current 600 vs. 800 Myr age discussion.

6.3 Circularization and synchronization

6.3.1 Tidal circularization

In this section we compare our findings for the new EB systems to the expectations for tidal circularization and spin-orbit synchronization at the age of Praesepe. The binaries presented here are particularly valuable benchmarks for studies of tidal dissipation timescales in close binaries, as they are at or near the beginning of their main sequence evolution. Zahn & Bouchet (1989) posited that essentially all tidal circularization should occur during the PMS phase, when stars are larger and have deeper convective envelopes. If this theory were correct, all late-type main sequence binaries with periods less than $\sim$8 days should be circularized. Binaries with longer orbital periods would retain their primordial eccentricities and experience negligible tidal circularization after the PMS phase.

However, Meibom & Mathieu (2005) used observations of binaries in coeval stellar populations to clearly show that tidal dissipation proceeds to circularize orbits well after the PMS stage (see their Fig. 9). While standard equilibrium tide theory (Zahn, 1989; Claret & Cunha, 1997) and dynamical tide theory (Witte & Savonije, 2002) do predict exactly this trend, binaries are generally observed to circularize more quickly than theory predicts (i.e. tidal dissipation is a more efficient process than expected). The binary population of Praesepe and the Hyades is a conspicuous outlier to this trend, indicating agreement with theory but significant tension with observations of all other well-characterized clusters. However, Zahn & Bouchet (1989) cautioned that two short-period eccentric binaries in Praesepe and Hyades (KW 181 and VB 121) are single-lined systems, in which the secondaries could possibly be white dwarfs, meaning that the standard theory of tidal dissipation would not apply. Ignoring these two systems, those authors estimated binaries with periods below 8.5–11.9 days should be circularized by the age of the Hyades, and by extension Praesepe. Our findings for AD 3814 and AD 2615 corroborate the notion that the circularization period for Praesepe is larger than previously measured, and to our knowledge AD 2615 is the longest period circular binary in either Praesepe or the Hyades. Revisiting the analysis of Meibom & Mathieu (2005) including these two systems would bring the observations for Praesepe into better agreement with those of other clusters, in the sense that binaries of a given age are observed to be circular out to longer periods than theory predicts.

As for AD 3116, tidal dissipation proceeds differently for extreme mass ratio systems (Ogilvie, 2014), and so we caution against drawing conclusions based on its relatively high eccentricity ($e=0.15$) given its short orbital period of $<2$ days. In fact, the recently discovered transiting brown dwarf in the significantly older Ruprecht 147 cluster similarly exhibits a relatively high eccentricity and short orbital period (Nowak et al., 2016).

Finally, we note that the transition between circular and eccentric binaries in a coeval stellar population (as demarcated by either the “cutoff period”, i.e. the longest period circular binary, or preferably by the “tidal circularization period”) can in principle be used to estimate the age of the stellar population (Mathieu & Mazeh, 1988). Given sufficient data and a well-calibrated relation amongst clusters, the method could also be extended to close binaries in the field to provide an upper limit in age if the binary is eccentric, or a lower limit if it is circular.

6.3.2 Spin-orbit synchronization

The theoretical outcome of tidal evolution within a binary system is a circular orbit and a state of double synchronous rotation with spin axes aligned to the orbital angular momentum vector. However, as noted by Ogilvie (2014), this theoretical prediction has never been observationally verified for a binary star system. This is in part due to the difficulty of measuring stellar rotation, particularly for both components of a binary, and the need for an eclipsing system to precisely measure obliquities.

Binaries for which the rotation period of one or more component can be measured, particularly within coeval stellar populations, are thus critical benchmarks for tidal synchronization studies. For the four binaries discussed here, one appears to be nearly synchronized (AD 1508) while the other three appear to be rotating subsynchronously (i.e. at a frequency lower than the orbital frequency). This observation is based on the measured $P_{\mathrm{spot}}/P_{\mathrm{orb}}$ ratios of 1.25, 1.08, and 1.14 for ADs 3814, 2615, and 3116, respectively. On the surface, this is surprising given that 1) the expected synchronization timescales are much smaller than the cluster age, and 2) tidal synchronization is expected to occur more quickly than circularization in close binaries (Zahn, 1977; Hut, 1981) and two of the subsynchronous binaries are on nearly circular orbits (ADs 3814 and 2615).

It is important to note that photometric variations indicate a star’s surface rotation rate, but the spin of interior layers is not measured and known only to the extent to which there is reason to believe the interior is coupled to the surface. For the binaries with mass ratios near unity, both stars are contributing to the observed brightness modulations, but in the absence of multiple distinct peaks in a periodogram we infer the modulation period to indicate the rotation of the primary. Notably, surface differential rotation can lead to configurations in which the spin of equatorial regions is synchronized with the orbit while higher latitudes may be rotating more slowly. Such a scenario was suggested to explain observations of the late-type EB HII 2407 in the Pleiades (David et al., 2015). Indeed, there is observational evidence (Barnes et al., 2015) and theoretical motivation (Schuessler & Solanki, 1992; Granzer et al., 2000) for polar spots on rapidly rotating, fully convective stars.

However, unlike the Pleiades EB, the binaries presented here exhibit much larger discrepancies between the rotation and orbital periods. The measured rates of differential rotation have been observed to decrease strongly with stellar temperature (Barnes et al., 2005; Collier Cameron, 2007). Using the empirical formula of Collier Cameron (2007), the expected rates of differential rotation for the stars considered here are all below $10^{-3}$ rad d*-1*. If we assume the orbits are synchronized at the equator and that polar spots are responsible for the measured rotation periods, then the implied rates of differential rotation for ADs 3814, 2615, and 3116 would be 0.21, 0.04, and 0.38 rad d*-1*, respectively. These values are significantly higher than the differential rotation rates measured for fully convective stars (Morin et al., 2008; Reinhold et al., 2013; Davenport et al., 2015). Our observations therefore indicate either: 1) tidal synchronization proceeds more slowly in fully convective stars than the theory of equilibrium tides predicts, 2) magnetic braking is currently playing a more important role in the spin evolution of these binaries than tidal forces, or 3) differential rotation in fully convective stars can be much more important than previously appreciated. We consider the last explanation to be the least plausible.

Subsynchronous rotation has previously been observed for short period binaries in the younger M35 and M34 clusters, aged $\sim$150 Myr and $\sim$250 Myr, respectively (Meibom et al., 2006). As those authors noted, this result is in direct contradiction with expectations of tidal evolution on the main sequence which predicts binaries with periods near or below the circularization period (which AD 3814 and AD 2615 apparently are) to be rotating pseudosynchronously (synchronized with the instantaneous orbital angular velocity at periastron) or slightly supersynchronous.

We conclude by noting that current theories of tidal evolution carry significant and under-explored uncertainties. In particular, theory for solar-type and early-type stars is more developed than that for fully convective stars. Tidal dissipation is expected to be more efficient, and thus circularization more rapid, in stars with convective outer layers (Zahn, 1975), which is supported observationally (Van Eylen et al., 2016).

7 Conclusions

We have presented photometric timeseries data from Kepler/K2 and follow-up high dispersion spectroscopy from Keck/HIRES in order to characterize four new EB systems in the sub-Gyr old Praesepe cluster. These new discoveries increase the number of characterized EBs below 1.5 $M_{\odot}$ in sub-Gyr open clusters by 20%, and add 40% of the cluster EB population with masses $M\lesssim 0.6\,M_{\odot}$.

We analyze these low-mass EBs with GP–EBOP, a new multi-purpose Gaussian process eclipsing binary and transiting exoplanet model, to determine model-independent stellar masses and radii. We present an updated method of simultaneously determining the effective temperatures of both stars as well as the distance to an EB by modeling the system’s spectral energy distribution. This approach capitalizes on the posterior constraints from the joint light curve and RV modeling to break existing degeneracies and also correctly interprets the light curve model’s band-specific surface brightness ratio, rather than using it to approximate an effective temperature ratio.

We determine the masses of AD 3814 to 2% precision and the primary and secondary radii to 1% and 3%, respectively. The masses and radii of AD 2615 are both determined to 6% precision for the primary and to 5% for the secondary. Together with effective temperatures determined to a typical precision of $\pm$50 K, we test the PARSEC v1.2 and BHAC15 stellar evolution models. Overall, the EB parameters are most consistent with the PARSEC models, primarily because the BHAC15 temperature scale is too hot over the mass–age range probed. Both the PARSEC and BHAC15 models are able to explain the masses and radii of ADs 3814 and 2615 with a single isochrone in the range $\sim$400–1000 Myr, but predicting $T_{\rm eff}$ proves more challenging. Our SED-derived $T_{\rm eff}$ values, which are consistent with those derived from empirical M-dwarf relations, are better matched to the PARSEC models. We find that the BHAC15 models predict temperatures $T_{\rm eff}$ $\sim$100–300 K hotter than our data, whereas the PARSEC models lie in the correct $T_{\rm eff}$ range. However, both models predict a steeper $T_{\rm eff}$ track over the mass range $M\sim 0.2-0.4$ $M_{\odot}$ than our data suggest. More M-dwarf EBs with precise $T_{\rm eff}$ values on the main sequence are required to confirm this tentative statement. Our luminosities are in agreement with the PARSEC model predictions but we find that the BHAC15 models overpredict this parameter primarily due to their high $T_{\rm eff}$ values. While both ADs 3814 and 2615 possess precise solutions, we note that AD 3814 would benefit from a more detailed modeling of the individual eclipses (especially incorporating a full starspot model), and AD 2615 would benefit from additional RVs to tighten the existing solution.

We present a preliminary solution for a third detached double-lined system, AD 1508. The K2 light curve displays clear, but shallow, eclipses on both stars and the three Keck/HIRES RVs we obtained show the two stars not to be rapid rotators. This system is therefore amenable to precise characterization but would require further RV measurements throughout the orbital phase and may also benefit from targeted eclipse monitoring with moderate-aperture ground-based telescopes.

The final system, AD 3116, comprises a mid M-dwarf primary star with a transiting brown dwarf companion (M$$\sim54 $M_{\rm Jup}$). There are only $\sim$20 transiting brown dwarf systems known: AD 3116 is one of only three systems where the primary is an M-dwarf, and is only the second transiting brown dwarf system discovered in an open cluster (and the first younger than a Gyr). It will therefore be a favorable target for future transiting brown dwarf studies.

Finally, we find that ADs 3814 and 2615, which have orbital periods of 6.0 and 11.6 days, are circularized but not synchronized, with at least one component rotating sub-synchronously. This contradicts the expectations of tidal evolution, which would predict synchronization to proceed faster than circularization in these systems and for it to have been achieved by the age of Praesepe. Our observations therefore suggest that either tidal synchronization proceeds more slowly in fully convective stars than the theory of equilibrium tides predicts, or magnetic braking is currently playing a more important role in the spin evolution of these binaries than tidal forces.

We thank Pierre Maxted for interesting discussions and John Southworth for help compiling Table 5. This paper includes data collected by the Kepler/K2 mission. Funding for the K2 mission of Kepler is provided by the NASA Science Mission directorate. Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) under support by the NASA Office of Space Science via grant NNX09AF08G and by other grants and contracts. Some of the data presented herein were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W.M. Keck Foundation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. Finally, we would like to thank the anonymous referee for their careful reading of the manuscript and helpful suggestions for improvement.