Fundamental properties of the pre-main sequence eclipsing stars of MML 53 and the mass of the tertiary

TL;DR

This study provides a detailed analysis of the MML 53 eclipsing binary system, confirming the third star's bound nature, measuring stellar parameters with high precision, and discussing implications for early binary formation and stellar evolution models.

Contribution

First comprehensive analysis of MML 53, including the third star, with precise measurements of stellar masses, radii, and temperature, and insights into early binary formation mechanisms.

Findings

Masses measured within 1% accuracy

Radii inflated by 10-20% compared to models

Third star's mass confirmed around 0.7 Msun

Abstract

We present the most comprehensive analysis to date of the Upper Centaurus Lupus eclipsing binary MML 53 (2.097892 d), and for the first time, confirm the bound-nature of the third star (~9yr orbit). Our analysis uses new and archival spectra and time-series photometry. We determined the temperature of the primary star to be 4880+-100 K. The study of the close binary incorporated treatment of spots and dilution by the tertiary in the light curves, allowing for the robust measurement of the masses of the eclipsing components within 1% (M1=1.0400+-0.0067 and M2=0.8907+-0.0058 Msun), their radii within 4.5% (R1=1.283+-0.043 and R2=1.107+-0.049 Rsun), and the secondary temperature (4379+-100K). From the analysis of the eclipse timings, and the change in systemic velocity of the eclipsing binary and the radial velocities of the third star, we measured the mass of the outer companion to be 0.7 Msun within 20%. The age we derived from the evolution of the temperature ratio between the eclipsing components is fully consistent with previous estimates of the age of UCL (16+-2 Myr). At this age, the tightening of the MML 53 eclipsing binary has already occurred, thus supporting close-binary formation mechanisms that act early in their evolution. The eclipsing stars roughly follow the same theoretical isochrone, but appear to be inflated in radius (by 20% for the primary and 10% for the secondary). However, our primary radius measurement of is in full agreement with the independent measurement of the secondary of NP Per which has the same mass and a similar age. The eclipsing stars of MML 53 are found to be larger but not cooler than predicted by non-magnetic models, it is not clear what is the mechanism that is causing the radius inflation given that activity, spots and/or magnetic fields slowing their contraction, require the inflated stars to be cooler to remain in thermal equilibrium.

Tables

Table 1. Table 1: MML 53 WASP light curve data

HJD–2 450 000†	$\Delta$mag	${\sigma }_{\mathrm{mag}}$
3860.38987	0.0028	0.0723
3860.39021	-0.0028	0.0625
3860.39687	0.0076	0.0743
3860.39731	0.0028	0.0743
3860.40458	0.0011	0.0559
…

Table 2. Table 2: MML 53 FTS V-band light curve data

HJDUTC–2 450 000†	$\Delta$mag	${\sigma }_{\mathrm{mag}}$
5731.99983	-0.1372	0.0015
5732.00075	-0.1346	0.0015
5732.00169	-0.1331	0.0015
5732.00265	-0.1295	0.0015
5732.00360	-0.1287	0.0015
…

Table 3. Table 3: MML 53 CTIO light curve data

BJDTDB–2 450 000	$\Delta$mag	${\sigma }_{\mathrm{mag}}$	Filter
4970.47351	0.024	0.009	U
4970.48376	0.054	0.009	U
4970.49238	0.087	0.009	U
4970.49889	0.124	0.009	U
4970.50535	0.167	0.009	U
…

Table 4. Table 4: Radial velocity measurements

Instrument	BJDTDB	Primary	Secondary	Tertiary
	$-\mathrm{2\hspace{0.17em}450\hspace{0.17em}000}$	RV	RV	RV
		(kms-1)	(kms-1)	(kms-1)
FEROSa	3909.62035	-85.8	111.1	-3.5
UVES	5026.66177	96.3b	-108.6c	11.1d
	5040.51000	-76.7	91.0	7.8
	5042.47920	-90.7	108.3	10.4
	5044.61337	-88.5	104.6	9.1
	5045.62608	92.3	-107.6	10.1
	5049.49566	69.6	-80.3	10.3
	5060.55105	56.8	-61.6	9.0
	5061.50408	-74.5	87.7	8.2
	5062.47025	87.1	-102.7	11.1
	5081.48640	67.5	-75.0	11.2
	5081.49009	66.0	-74.0	11.6
	5083.48843	82.2	-96.3	11.5
	5084.48765	-87.5	102.1	10.3
	5092.49966	-68.4	81.6	8.9
CTIO 2010e	5449.48870	-92	103	18

Table 5. Table 5: Systemic EB and tertiary radial velocities

Instrument	BJDTDB	${\gamma }_{EB}$	RV3
+Epoch	$-\mathrm{2\hspace{0.17em}450\hspace{0.17em}000}$	($\mathrm{km}{s}^{-1}$)	($\mathrm{km}{s}^{-1}$)
FEROS 2006	3909.62035	$5.2\pm 0.8$	$-3.5\pm 1.1$
UVES 2009	5061.24102	$0.76\pm 0.15$	$10.1\pm 1.3$
CTIO 2009	4982.22138	$1.4\pm 0.9$	$11.0\pm 3.0$
CTIO 2010	5449.48870	$-2.0\pm 1.8$	$18.0\pm 2.5$

Table 6. Table 6: Stellar Surface spot parameters on the primary star

	Colatitude	Longitude	Radius	Temperature
	(deg)	(deg)	(rad)	factor
Spot 1	45	0	30	0.94
Spot 2	90	135	10	0.85

Table 7. Table 7: Constraint on level of third light from EB light curve model

	Minimum	Maximum
	Dilution	Dilution
$L_3/(L_1+L_2+L_3$)
in U	0.0	0.55
in B	0.0	0.56
in V	0.0	0.59
in RC	0.0	0.60
in IC	0.0	0.63
$i{(}^{\circ })$	81.4	90.0
${\chi }_{\mathrm{reduced}}^2$	2.2	2.5
$\Delta$ $M_1$†	3.3%
$\Delta$ $M_2$†	3.3%
$\Delta$ $R_1$†	1.1%
$\Delta$ $R_2$†	1.1%
$\Delta$ $T_{\mathrm{eff},2}$/$T_{\mathrm{eff},1}$†	2.4%
† Difference in value between models with minimum and maximum dilution

Table 8. Table 8: Times of light curve minima for primary eclipse

Dataset	Time of minimum
	(BJD${}_{\mathrm{TDB}}{}^{}-\mathrm{2\hspace{0.17em}450\hspace{0.17em}000}$)
WASP 2006	3911.1182 $\pm$ 0.0007
WASP 2007	4227.9021 $\pm$ 0.0004
WASP 2008	4563.5651 $\pm$ 0.0005
CTIO May 2009	4972.6509 $\pm$ 0.0002
WASP 2011	5692.2197 $\pm$ 0.0004
FTS-20110619	5732.0764 $\pm$ 0.0001
FTS-20110709	5750.9574 $\pm$ 0.0001
WASP Feb 2012	6038.3701 $\pm$ 0.0002
WASP 2013	6418.0912 $\pm$ 0.0002

Table 9. Table 9: Orbital solution of EB and tertiary

Table 10. Table 10: Flux of each stellar component relative to the total flux

	V-band	R-band	I-band
$L_1/(L_1+L_2+L_3$)	0.65	0.61	0.55
$L_2/(L_1+L_2+L_3$)	0.24	0.26	0.28
$L_3/(L_1+L_2+L_3$)	0.11	0.13	0.17

Table 11. Table 11: Physical properties of MML53 eclipsing stars and their orbit

Parameter		Value	Units
Orbital period	$P_{\mathrm{orb}}$‡	2.097892 $\pm$ 0.000005	days
Eccentricity	$e$	0.	(fixed)
Mass ratio	$q_{\mathrm{EB}}$ = $M_2$/$M_1$	0.8565 $\pm$ 0.0034
Systemic velocity	${\gamma }_{\mathrm{EB}}$ †	0.77 $\pm$ 0.15	$\mathrm{km}{s}^{-1}$
Semi-major axis	$a\sin i$	8.492 $\pm$ 0.017	R⊙
	$a$	8.584 $\pm$ 0.018	R⊙
	$a$	0.03992 $\pm$ 0.00008	au
Binary total mass	$M_B$	1.9307 $\pm$ 0.0119	M⊙
Inclination	$i$	81.61 $\pm$ 0.12	∘
Sum of fractional radii	$r_1+r_2$	0.2784 $\pm$ 0.0027
Temperature ratio	$T_{\mathrm{eff},2}$/$T_{\mathrm{eff},1}$	0.8972 $\pm$ 0.0018
Primary mass	$M_1$	1.0400 $\pm$ 0.0067	M⊙
Primary radius	$R_1$	1.283 $\pm$ 0.043	R⊙
Primary surface gravity	$\log g_1$	4.24 $\pm$ 0.03	dex (cgs)
Primary temperature	$T_{\mathrm{eff},1}$	4880 $\pm$ 100	K
Secondary mass	$M_2$	0.8907 $\pm$ 0.0058	M⊙
Secondary radius	$R_2$	1.107 $\pm$ 0.049	R⊙
Secondary surface gravity	$\log g_2$	4.30 $\pm$ 0.04	dex (cgs)
Secondary temperature	$T_{\mathrm{eff},2}$	4379 $\pm$ 100	K
‡ The orbital period was adopted from the analysis of Hebb et al. (2010).
† Systemic velocity of the EB components at the time of the UVES RV observations.

Full text

1 Introduction

With the growing number of transiting planet surveys and follow-up radial velocity data, the detection and characterization of eclipsing binary (EB) stars has seen a resurgence over the last decade. Eclipsing binaries that are also double-lined spectroscopic systems have long provided crucial observational constraints for stellar evolution models by allowing the direct measurements of the masses and radii of the components, and also importantly, a measure of their temperatures (Andersen 1991; Torres et al. 2010; Stassun et al. 2014). Well-constrained stellar masses and radii are especially important for understanding pre-main-sequence (PMS) stars, as these are rapidly evolving systems which have not yet fully contracted.

Only ten years ago, the known pre-main sequence EBs with precisely measured properties were all members of the Orion nebula cluster, with ages between 1 and 2 Myr probing the youngest and earliest stages of stellar evolution, and of the Orion OB1 group, with an older stellar population of $\sim$10 Myr (e.g., Mathieu et al. 2007). It has been over the last few years that new EB systems in other young clusters and associations have been discovered, observed and carefully analyzed. Currently in the literature, there are only 14 known EB systems where the eclipsing components have directly measured masses ($M_{\star}$ $\lesssim$ 1.4 M⊙) and are on the pre-main sequence. Of these, there are seven EB systems belonging to the Orion star formation complex: ASAS J052821+0338.5 (Stempels et al. 2008); RX J0529.4+0041 (Covino et al. 2001, 2004); V1174 Ori (Stassun et al. 2004); Parenago 1802 (Stassun et al. 2008; Cargile et al. 2008; Gómez Maqueo Chew et al. 2012); Parenago 2017 (Morales-Calderón et al. 2012); JW 380 (Irwin et al. 2007), and 2MASS J05352184–0546085 (Stassun et al. 2006, 2007; Gómez Maqueo Chew et al. 2009). There are five EBs that are members of the Scorpius-Centaurus OB complex: HD 144548 (Kiraga 2012; Alonso et al. 2015); MML 53 (Hebb et al. 2010, 2011); UScoCTIO 5 (Kraus et al. 2015; David et al. 2016); EPIC 203710387 (Lodieu et al. 2015; David et al. 2016), and EPIC 203868608 (David et al. 2016). And there is one known EB in NGC 2264, CoRoT 223992193 (Gillen et al. 2014, 2017), and one EB in the Perseus star-forming complex, NP Per (Perova et al. 1966; Lacy et al. 2016). Other pre-main sequence EB candidates have been identified but their fundamental properties have not yet been measured (e.g., van Eyken et al. 2011; Morales-Calderón et al. 2012). Given the large spread of measured masses (from brown dwarfs of $\sim$0.02 M⊙ up to $\sim$1.4 M⊙ stars) and radii ($\sim$0.25 to 2.4 R⊙) of the known eclipsing objects (see also Fig. 11) and the spread in ages ($\sim$1–17 Myr) of their star-formation regions, it is clear that analyses of these systems provides strong empirical constraints on models of pre-main sequence evolution of low-mass stars and brown dwarfs.

MML 53, the first pre-main sequence EB discovered outside of the Orion star forming region and the subject of this paper, is an interesting pre-main sequence EB. Its young, pre-main sequence nature has been comfirmed by numerous observations measuring the X-ray emission, H$\alpha$ emission, and Li I $\lambda$6708 absorption from the system (Wichmann et al. 1997; Mamajek et al. 2002; Torres et al. 2006; White et al. 2007; Hebb et al. 2010), and it has long been known as a spatial and kinematic member of the $16\pm 2$ Myr old Upper Centaurus Lupus (UCL) subgroup of the Scorpius-Centaurus OB association (Mamajek et al. 2002; Pecaut & Mamajek 2016; Pecaut et al. 2012).

The eclipsing nature of the system was first discovered by Hebb et al. (2010) in data obtained as part of the WASP transiting planet survey (Pollacco et al. 2006). Analysis of the 2006–2008 WASP light curve combined with additional radial velocity measurements taken in 2009 with the 1.5m telescope at Cerro Tololo Inter-American Observatory (CTIO) found the EB to be composed of a $1.0$ M*⊙* and a 0.86 M*⊙* pair of stars in an 2.09 day eclipsing orbit (Hebb et al. 2011). Features from a third unresolved star were also detected in the spectra. As we show in this paper, variations in the radial velocity of the tertiary component compared to the systemic radial velocity of the binary confirm MML 53 is a gravitationally bound triple system. In addition, Hebb et al. (2010) detected small changes of $\sim$3 minutes in the epoch of the eclipses from 2006–2008. Subsequent WASP observations described in this work have confirmed these variations, which are attributed to light travel time effects. As the EB in MML 53 orbits the third star over the timescale of about a decade, the distance from the Earth to the EB changes causing the epoch of the eclipse minima to vary with its orbital position.

In summary, MML 53 is a 16 Myr old, hierarchical triple system consisting of a close eclipsing binary and a lower mass tertiary component that has recently been spatially resolved (Schaefer et al. 2018). Due to its unique age among the known pre-main sequence EBs, precise measurements of the fundamental properties of its component stars have the potential to test a previously unconstrained part of parameter space in the theoretical stellar evolution models. In this paper, we present precise fundamental properties of the eclipsing components derived by incorporating new, high-quality spectroscopic and photometric observations of the MML 53 system into a comprehensive eclipsing binary model (§3.5), which accounts for stellar surface spots and the effect of light from the tertiary star in the light curve modeling. This paper also presents the first constraints on the mass and orbital parameters of the third stellar component through a combined analysis of long term variations in the systemtic radial velocity of the EB and corresponding changes in its measured eclipses times.

2 Observations

In this section, we describe the new and archival observations utilized in our analysis of the MML 53 system for deriving the fundamental properties of its eclipsing stars and the orbital parameters of the bound tertiary component.

2.1 Photometric data

2.1.1 WASP photometry

Hebb et al. (2010) described the photometric time series data obtained on MML 53 between 2006–2008 as part of the WASP transiting planet survey (Pollacco et al. 2006). Subsequently, MML 53 was observed again in the field-of-view of the WASP-South telescope between 2011–2013. All nights showing full or partial eclipses were extracted from the full WASP light curve and used to measure the epoch of minimum light of the eclipsing pair for each year of data between 2006–2013 as described in §3.3.1. A total of 10328 photometric data points were obtained between February–August 2011; 10648 were obtained between February–June 2012 with an additional 23952 points in July and August 2012; and 76754 data points were observed using three cameras in an intensive campaign of this field between February–August 2013. All WASP photometric data of MML 53 are provided in Table 1 (in full in online version111Tables 1–3 are only available in full in electronic form at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsweb.u-strasbg.fr/cgi-bin/qcat?J/A+A/).

These data were processed and removed of systematics with the standard WASP pipeline (Collier Cameron et al. 2006) resulting in 135078 brightness measurements obtained over this time period. The typical photometric precision of the early 2006–2012 data is $\sim 7$ mmag, which is measured by the standard deviation of the data points in the out-of-eclipse phases. After changing the WASP-South lenses from 200mm to 85mm in July 2012, the updated observing strategy lead to more observed data points with a lower precision of $\sim 20$ mmag. Starspot modulations are present in these data from which we can measure a rotational period (§3.2). However, the starspot modulation can affect the derived time of the eclipse epochs if they are not modeled correctly or removed. Therefore, each night of data was rectified by fitting a first or second order polynomial baseline to the out-of-eclipse data and subtracting the model values from the observed magnitudes at all times. Individual out-of-eclipse data points were rejected at this stage if they deviate by more than 5$\sigma$ from the polynomial baseline. The resulting phase-folded primary and secondary eclipses derived from the rectified light curves are shown in Fig. 1 for each year of WASP data. We omitted the July–August 2012 data since all eclipses occurred while the sun was up due to the near integer day period of the EB. All data were phase-folded with the, $P_{\mathrm{orb}}$ = 2.097892 $\pm$ 0.000005 days, and time of minimum light in BJDTDB units, T${}_{0}=2454972.650850$, derived from the detailed modeling of the 2009 CTIO data. Noticeable shifts in the time of eclipse minima as compared to T0 and due to the tertiary component are visible in these data. The EB model light curves (described in §3.3.1) are over plotted in red on all rectified light curves. The times of minimum light that change from year to year are defined with these EB models and are used to constrain the parameters of the tertiary’s orbit in §3.3.

2.1.2 Faulkes Telescope South photometry

Two primary eclipses of MML 53 were observed on 19 June 2011 and 9 July 2011 in order to continue tracking the eclipse timing variations. The data were obtained with the Spectral Camera on the 2-m Faulkes Telescope South (FTS) through the Las Cumbres Observatory Global Telescope Network (LCOGT). We observed the target with 60 second exposure times repeatedly for approximately 5.5 hours in the Johnson V-band filter. We employed the $2\times 2$ binning mode for faster readout time and defocused the camera by 0.3 mm to avoid saturation. The data were processed in the standard way with the LCOGT imaging data pipeline (BANZI) 222https://lco.global/observatory/data/BANZAIpipeline/, which includes bad pixel masking, bias and dark frame subtraction, and flat-field division of each individual science frame with the best available calibration images. The pipeline also performs source extraction and astrometry. The $5^{\prime}\times 5^{\prime}$ field-of-view of the instrument contained eight bright comparison stars that were used in deriving the differential magnitudes with a photometric precision of $\sim 2$ mmag. The phase-folded FTS light curves are shown in Fig. 2, and the data are given in Table 2.

Attempts were made to get additional eclipse photometry in the 2012 and 2013 seasons, but poor weather prevented such observations.

2.1.3 CTIO photometry

MML 53 was observed between May 18 and June 08, 2009 with the CTIO-1m telescope and Y4K-Cam camera. The detector consists of a 4K$\times$4K array of $15\mu$ pixels placed at Cassegrain focus giving a $0.3^{\prime\prime}$/pixel platescale. Thus the entire array projects to a $20^{\prime}\times 20^{\prime}$ field of view. The observed signal is fed into four amplifiers causing the raw images to have a quandrant effect with the readnoise between 11-12 $e^{-}$ and gain of 1.45-1.52 $e^{-}$/ADU, depending on the amplifier. The detector has a readout time of 51 seconds and a 71k-electron well depth before non-linearity sets in. This converts to a saturation of 40,000 counts/pixel in $1\times 1$ binning mode.

Throughout each observing night, MML 53 and the surrounding field were monitored in the standard Kron-Cousins optical filter set (UBVRcIc) alternating continuously between all five filters. Exposure times were chosen to maximize the flux in the target star and nearby reference stars while keeping the peak pixel value in MML 53 below $40,000$ counts. The telescope was defocused to allow for longer exposure times to build up signal in the fainter reference stars without saturating MML 53. We adopted an exposure time of 7 seconds for the V, RC, and IC–band observations and longer exposures of 45 seconds and 90 seconds in the B and U band filters, respectively, where the detector is less sensitive. We achieved an overall light curve cadence of approximately 8 minutes in each filter accounting for the exposure times, the read out time, and other overheads, like filter changes. Since the orbital period of the eclipsing binary is very close to $\sim$2 days, three consecutive primary eclipses and three secondary eclipses were observed during the first six (6) clear nights of the observing campaign. Four of these six nights were photometric. On the other two nights, thin clouds were visible, but it did not affect the overall observing cadence or photometric precision of the data. Due to poor weather, MML 53 was observed sparsely for the next seven nights (2009-05-24 to 2009-05-30). The weather improved for the final week of the observing campaign, but all the eclipses occurred during the day, so these data only sample the out-of-eclipse variation and allow the characterization of the stellar spots (§3.2.1). The eclipsing binary analysis described below uses only the first six (6) nights of data in which the eclipses occur. These data are presented in Table 3.

Flat field and bias calibration frames necessary for processing the images were obtained during each observing night. Sets of 11 bias frames were taken at the beginning and end of each night, and single frames were observed periodically throughout each night. Eleven dome flats were observed per night in all five filters, and twilight flats (3-4 per filter) were obtained on the few photometric nights in the beginning of the run. Due to the relatively small number of twilight flats obtained in each filter, the dome flats were used for the flat-field calibration correction.

The images were processed in a standard way using routines written by L. Hebb in the IDL programming language. Each of the four amplifiers was processed independently. All object and calibration frames were first overscan corrected (by subtracting a line-by-line median overscan value), bias subtracted and then trimmed. Stacked bias images were created by averaging all bias frames observed each night and subtracted from all science and flat-field frames. All dome flats observed during the first six nights of the observing campaign were averaged into a single dome flat in each filter and then applied to the trimmed and bias-corrected science images.

Souce detection and aperture photometry were performed on all processed science images using the Cambridge Astronomical Survey Unit catalog extraction software (Irwin & Lewis 2001). The software has been compared with SExtractor (Bertin & Arnouts 1996) and found to be very similar in the completeness, astrometry and photometry tests 333https://www.ast.cam.ac.uk/ioa/research/vdfs/docs/reports/simul/index.html. This photometry software was applied to all processed images of MML 53. Adopting conservative parameters to define the detection threshold, the target star and dozens of fainter stars in the field were detected in each image. Aperture photometry was performed on all detected stars using a 4 pixel radius circular aperture, which was selected to match the typical seeing over the first six nights of the observing run. The same aperture was used on all nights of data. Eight bright, non-variable reference stars were selected from the many detected stars and used to perform differential photometry on the target star. In each image, the flux from all reference stars was summed into a single super comparison star that was divided by the aperature flux from MML 53 and converted to a differential magnitude. The resulting phase folded differential photometry light curves of MML 53 obtained from the first six nights of the observing run are shown in Fig. 3 (VRCIC) and Fig. 4 (UB).

2.2 Spectroscopic data

The spectroscopic data presented in this paper are used to model the short term radial velocity variations of the primary and secondary EB components and to track the long term secular variations in velocity as the EB and the tertiary orbit their common center of mass.

2.2.1 UVES spectra

MML 53 was observed fourteen times between 14 July 2009 and 18 September 2009 with the UVES spectrograph on the ESO Very Large Telescope (program ID 383.C-080). The observations were obtained with the dichroic mode on the instrument, with the blue arm centered at 3900Å, and the red arm at 5800Å. Data in the blue arm are of poor quality and were not considered in this paper. We adopted a slit width of 0.6*′′* which allows for achieving a resolution, R$\sim 60,000$ at the red end of the spectrum. With exposure times of 3 minutes per observation, we achieved a signal-to-noise of 100 on the V$\sim 10.8$ star. The data were processed with the REDUCE package (Piskunov & Valenti 2002), which uses advanced order-tracing and slit-modeling techniques to reconstruct and extract the stellar spectrum.

2.2.2 FEROS spectrum

One high resolution ($R\sim 50,000$) spectrum of MML 53 covering a wavelength range between 3765 and 8862 Å was found in the European Southern Observatory (ESO) archive. The spectrum was obtained on 23 June 2006 using the FEROS échelle spectrograph on the 2.2m MPG/ESO telescope. This spectrum was presented in Hebb et al. (2010) where it was used to confirm the presence of the tertiary star and measure the radial velocity values of all three components. The details of the radial velocity analysis can be found in that paper, but the measured velocities are $-85.8,111.1$, and $-3.5$ km s*-1*, for the primary secondary and tertiary, respectively. Based on our experience, we adopted an uncertainty of 1.1 km s*-1* for the an individual radial velocity (RV) measurement from this instrument. Using the mass ratio derived from the final EB analysis (§11) and these primary and secondary star RV measurements, we determineds the systemic radial velocity for the EB to be $5.2\pm 0.8$ km s*-1* at the time of this observation. We report this value in the Table 5 and use it in the binary-tertiary analysis (§3.3).

2.2.3 CTIO spectra

A series of thirteen spectra of MML 53 were obtained in queue mode between 18 May 2009 and 12 June 2009 with the SMARTS 1.5m échelle spectrograph444See http://www.ctio.noao.edu/~atokovin/echelle/index.htm. at the Cerro Tololo Inter-American Observatory (CTIO). We also observed a single spectrum with the same instrument the following season on 09 September 2010 to continue monitoring the radial velocity variations in the tertiary star. A detailed description of the processing and analysis of the 2009 data are presented in Hebb et al. (2011), which we summarize briefly here since it is the same for the newly presented 2010 spectrum.

The bench-mounted spectrograph has a fixed cross-disperser and échelle grating, but accommodates a variety of slit widths that allow for resolutions of 25,000–40,000. In order to maximize the signal-to-noise in these observations, we obtained $3\times 600$s exposures each night with a large slit width of 140 $\mu m$  which translates into a signal-to-noise of $S/N\sim$25 per resolution element and a resolution of R$\sim$25,000, which is sufficient to identify and resolve the three individual components of MML 53. The spectral images taken on each night were processed in the standard way with overscan subtraction, 2-D bias subtraction, trimming, and flat-fielding before the three individual images were median combined. The spectra were extracted from each processed science frame and then wavelength calibrated with nightly ThAr lamp exposures using standard échelle data processing routines in IRAF555IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation (Tody 1993).. A single radial velocity standard was obtained on each night of the science observations and processed in an identical fashion. In 2009, a single spectrum of HD 81797 was used as the radial velocity template in the cross-correlation analysis, and in the 2010 season, a single 60 second exposure of HD 223807 was observed for the same reason. This star has a radial velocity of $-15.83$ $\mathrm{km\,s}^{-1}$(Nidever et al. 2002).

A cross-correlation analysis using the IRAF routine fxcor was performed on the calibrated spectra obtained in 2009 to measure twelve and ten independent radial velocities for the primary and secondary components, respectively. This analysis, presented in Hebb et al. (2011), resulted in measurements of the mass ratio and the systemic radial velocity of the EB of $+1.4\pm 0.9$ km s*-1*. During this time, the radial velocity of the tertiary star was also measured using fxcor in five spectra obtained near quadrature. The average radial velocity of the tertiary derived from these spectra is $+11.0\pm 3.0$ km s*-1*. In this paper, we analyzed the reduced spectrum from 2010 as described in §3.1 and derived the radial velocity of all three components.

3 Analysis and results

The various analysis steps to characterize this system are not independent. First, the determination of the mass of the tertiary body described in Sect. 3.3 depends on knowing the sum of the masses in the eclipsing binary ($M_{\rm B}$ = $M_{\mathrm{1}}$ + $M_{\mathrm{2}}$) which is derived from the EB model described in Sect. 3.5. However, the final EB model solution depends on knowing the value of the third light, and the third light depends on the mass of the tertiary. Furthermore, the spectral disentangling (described in Sect. 3.1) requires the relative luminosity ratios of the three unresolved components of the system which depends on the EB model and the binary-tertiary model. Finally, the spectral synthesis necessary to determine the temperature of the primary star requires knowledge of the gravity of the primary and secondary stars determined from the EB model. Therefore, the analysis steps described below were performed in an iterative manner until all solutions were consistent with each other, and the derived properties are measured robustly. Below, we describe the details of each analysis step and the final results derived from it during the final iteration.

3.1 Spectroscopic analysis

The high-resolution UVES spectra obtained in 2009 spectrally resolve all three components of MML 53, and cover a range of orbital phases of the system. This allowed us to determine the individual radial velocities of the components, as well as to perform a spectroscopic analysis of the three stars in this system.

Using the method of least-squares deconvolution (LSD, see Donati et al. 1997; Kochukhov et al. 2010) we determined from each spectrum combined-average line profiles of MML 53, concentrating on the region between 5500Å and 6500Å which is populated by a large number of narrow absorption lines. We recovered three profiles one for each star, which were separated in velocity-space. To each of the recovered profiles we then fitted a three-component model consisting of three rotational profiles, calculated by disk-integrated radiative transfer, which includes limb-darkening and non-rotational broadening processes such as micro- and macroturbulence. In this analysis, we used a macroturbulence value of 1.2 $\mathrm{km\,s}^{-1}$, and a microturbulence value of 1.0 $\mathrm{km\,s}^{-1}$, as appropriate for pre-main sequence stars (e.g., Padgett 1996). This allowed us to determine radial and projected rotational velocities for each component. The recovered radial velocities were measured relative to synthetic spectra, which are based on laboratory wavelengths from VALD (Piskunov et al. 1995; Kupka et al. 1999), and are presented in Table 4, and used in Sects. 3.2 and 3.5. The measured rotational velocities ($\upsilon\sin i$) are $30.6\pm 1.0$ $\mathrm{km\,s}^{-1}$, $26.6\pm 1.3$ $\mathrm{km\,s}^{-1}$, and $25.8\pm 3.2$ $\mathrm{km\,s}^{-1}$ for the primary, secondary and tertiary components, respectively. The errors in the measured $\upsilon\sin i$ were derived from the standard deviation of the $\upsilon\sin i$ derived from each of the UVES spectra (0.3, 0.8, and 3 $\mathrm{km\,s}^{-1}$, for the primary, secondary and tertiary respectively) to which we added in quadrature 1.0 $\mathrm{km\,s}^{-1}$uncertainty to account for systematic errors due to the continuum normalization and/or macroturbulence and microturbulence values used in the analysis.

We also analyzed the reduced CTIO spectrum from 2010 using LSD and find the radial velocity of all three components to be $-92$, $103$, and $18$ km s*-1*, for the primary secondary and tertiary, respectively. An uncertainty of $2.5$ km s*-1* was adopted for these measurements. Using the mass ratio derived from the final EB analysis (§11) and these primary and secondary star RV measurements, we determine the systemic radial velocity for the EB to be $-2.0\pm 1.8$ km s*-1* at the time of this observation. We report this value in the Table 5 and use it in the binary-tertiary analysis (§3.3).

We also applied the method of spectral disentangling (see Bagnuolo & Gies 1991) to our set of UVES spectra. This technique inverts the relation that each observed spectrum is a linear combination of the spectrum of each of the three stellar components. Using the radial velocities recovered above as input we numerically reconstructed the spectra of the three components as described in Stempels & Hebb (2011). Since this technique requires an assumption of the relative luminosity ratio of the three components, we adopt the values presented in Table 10 for this parameter. The recovered spectra are typical of young K-type dwarfs, with Li i 6708Å absorption clearly present. Both the secondary and the tertiary components are affected by narrow emission in H$\alpha$ and He i 5876Å; no emission was present in the spectrum of the primary star, although the H$\alpha$ line appeared to be filled-in. Because of the composite spectra, a more detailed analysis of the activity of the MML 53 stars is beyond the scope of this paper.

Once extracted, we calculated synthetic spectra for each component with the SME software package (Valenti & Piskunov 1996; Piskunov & Valenti 2017), using MARCS model atmospheres (Gustafsson et al. 2008) and atomic and molecular line lists from the VALD database (Piskunov et al. 1995; Kupka et al. 1999). While it is possible to estimate the surface gravity ($\log g$) from the spectrum, this parameter is much better determined from the masses and radii recovered from the EB modeling (see §3.5), and these values were therefore used as an input value when determining the effective surface temperatures $T_{\mathrm{eff}}$ for the eclipsing stars. From this analysis we recover for the primary, secondary, and tertiary $T_{\mathrm{eff}}$ = 4880 $\pm$ 100 K, $T_{\mathrm{eff}}$ = 4482 $\pm$ 100 K, and $T_{\mathrm{eff}}$ = 4500 $\pm$ 250 K, respectively. In the calculation of synthetic spectra, the metallicity was assumed to be solar, and the micro- and macroturbulence were estimated to be $1.0$ and $1.2$ $\mathrm{km\,s}^{-1}$.

The agreement of observed and disentangled spectra with synthetic spectra based on these parameters is illustrated in Fig. 6. In our analysis we find that the overall agreement between the synthetic and observed spectra is excellent for the primary star. Thus, we are highly confident in the primary star parameters derived from this synthesis. However, for the secondary star, we required to include an additional continuum source caused by magnetic activity, also referred to as “veiling”. Without the veiling corresponding to 20% of the light of the secondary star, the depth of calculated Na D lines are consistently too deep, while the shape (which is highly temperature-dependent) corresponds well to the observed spectrum. Also, H$\alpha$ emission suggests significant levels of magnetic activity is present in both the secondary and tertiary spectra. This unknown veiling quantity can have a moderate affect on the derived parameters. The tertiary disentangled spectrum is also in good agreement with the observation, however the properties are more uncertain because the luminosity of the tertiary is much lower than the other two components (as reflected in the uncertainty in the derived effective temperature).

Fortunately, only the primary star temperature and $\upsilon\sin i$ are necessary (along with the EB model solution) to derive all individual properties of the binary components, which is the aim of this paper.

3.2 Preliminary eclipsing binary model

For the analysis of the eclipsing components, we utilized the information from previous studies of the MML 53 system and adopt the assumptions described in this section.

The orbital period of the eclipsing binary was adopted from the careful determination presented in the discovery paper from the detailed analysis of the times of the eclipses (Hebb et al. 2010).

Utilizing all the available radial velocity data and light curves, we explored the possibility of a non-circular orbit. None of the solutions that allowed for a non-zero eccentricity converged. Moreover, the sinusoidal shape of the RV curves and the fact that the primary eclipse occurs at phase 0.0 and the secondary at phase 0.5 are robust indicators that the orbit is circular. Thus, for the rest of this analysis we adopt a circular orbit for the eclipsing binary (with eccentricity $e$ = 0.0).

We applied a Lomb-Scargle periodogram to each independent WASP light curve after removing the primary and secondary eclipses. Searching for periods between 0.5–30 d resulted in seven of the nine light curves having the strongest peak close to the same rotation period. We averaged the seven independent periods to determine a rotation period of 2.091 $\pm$ 0.013 d from the spot modulation present in the WASP light curves. Fully consistent with the orbital period, the components are found to be rotationally synchronized to their orbital motion, as expected for a circular orbit given that the tidal circularization timescale is expected to be longer than the synchronization timescale (e.g, Zahn 1977; Mazeh 2008).

Given that the spin-orbit alignment timescale is of the order of the synchronization timescale, we also assume that the stellar spin axes are aligned with the plane of the eclipsing binary orbit. Calculating the condition producing spin-orbit misalignment in the inner binary due to a tertiary component from Anderson et al. (2017), we find that in the case of MML 53 the eclipsing components of the binary are not likely to be misaligned.

Given that the UVES RV curves were taken over a period of $\sim$66 days (corresponding to $\sim$2% of the tertiary orbit) and that the peak-to-peak RV variation of the center of mass velocity of the EB ($\gamma_{\rm EB}$) is $\sim$10 $\mathrm{km\,s}^{-1}$, we considered the change in $\gamma_{\rm EB}$ to be negligible over the timespan that the UVES observations encompass.

Because the Baraffe et al. (2015) stellar evolutionary models do not provide constraints in the U and B broadband filters, and thus, there are no constraints on the third light in those bands, we did not use the UB light curves to derive physical properties of the eclipsing components.

Only in this preliminary model, we adopted the effective temperature of the primary component from the previous spectroscopic determination ($T_{\mathrm{eff,1}}$= 4890 K; Hebb et al. 2010). The primary temperature was updated from the analysis in §3.1 for the final EB model (§3.5).

Only in this preliminary model, the level of the third light was based on the relative heights of the CCF peaks from Hebb et al. (2010), and a coeval, lower mass star predicted from stellar models (e.g., Baraffe et al. 2015). Thus, we utilized a third light that corresponds to 9%, 11% and 18% of the total light for the VRCIC pass bands, respectively, as the dilution in the light curves due to the tertiary. For the final EB model (§3.5), we used the levels of third light derived in Sect. 3.4.

Given the above, we first fitted with PHOEBE (Prša & Zwitter 2005), the available CTIO light curves to derive the time of mid-transit. We then fitted the two radial velocity curves to derive the EB parameters that are fully defined by the RV curves, namely: $a\sin{i}$, mass ratio $q_{\rm EB}$, and systemic velocity $\gamma_{\rm EB}$ at the time of the UVES observations. These values are reported at the top of Table 11, and remained fixed for the rest of the analysis. The radial velocity curves and the best-fit RV model are shown in Fig. 5. We then reached a solution manually and utilizing the PHOEBE Levenberg–Marquardt (LM) minimization algorithm to determine a model that fits the VRCIC light curves well. The solution was attained by varying the inclination, the potentials, and the secondary temperature. At each step the limb-darkening coefficients were interpolated for each passband. This provides estimates for the radii and thus surface gravity of the eclipsing components, and the secondary temperature that were used in the determination of the spot properties below (see §3.2.1). Iteratively, such that we derived a consistent solution for the RV curves and the light curves, we also refine the time of mid-transit, which is derived to be $2454972.65085\pm 0.00016$ days (BJDTDB).

3.2.1 Stellar surface spots in the CTIO 2009 photometry

The presence of stellar surface spots is most evident from the out-of-eclipse phases of the light curves. We measured a peak-to-peak amplitude of 0.14 mag in the U-band, 0.09 mag in the B-band, 0.06 mag in the V-band, 0.04 mag in the RC-band, and 0.03 mag in the IC-band. These measured amplitudes are at least $5\times$ larger than the corresponding median photometric uncertainty of the CTIO light curves. The amplitude of the variation in the out-of-eclipse light curves of MML 53 increases with decreasing observed wavelength, as expected for stellar surface spots (e.g., Bouvier & Bertout 1989; Gómez Maqueo Chew et al. 2009). Additionally, the eclipsing binary is a detached system, meaning that the components are not interacting (i.e., there is no mass transfer). The eclipsing components are far enough away from each other that reflection effects have amplitudes that are smaller than the photometric precision of each light curve ($<$ 5 mmag; Wilson 1990), and little deformation of the stars occurs ($|r_{equator,i}-r_{pole,i}|<1\%$). Finally, the asymmetry of the light curves before and after the secondary eclipse (see Figs. 3 and 4 from phase 0.4 to 0.6) indicates that there is a non-homogeneous distribution of surface brightness in the combined light from the stellar disks (unlike ellipsoidal variation). The result of these qualitative observations leads to the conclusion that the observed deviation of the out-of-eclipse phases from the relatively flat light curve is most likely due to stellar surface spots.

The depth and shape of the eclipses are affected by the presence of spots, and consequently so are the derived radii (e.g., Covino et al. 2004; Morales et al. 2010; Windmiller et al. 2010) and the temperature ratio (Gómez Maqueo Chew et al. 2009). In order to characterize this in-homogeneity of surface brightness and derive robust physical properties, we attempted to model the light curve features with cooler surface spots. The spot parameters are degenerate, and we have very little constraint given our data sets on their properties. We have some information about the longitude of the spots based on the position of the deepest modulation in the out-of-eclipse light curve, and about the temperature of the spots relative to the stellar temperature from the relative depth of the spot modulation as a function of wavelength. However, the spot size, temperature, and latitude are all highly degenerate and multiple combinations of parameters can easily produce the same light curve variation. Despite the degeneracy, in order to study the stellar properties, we need to only adopt a single set of parameters that represent the light curve, thus minimizing the effect of the spots on the derived bulk properties of the eclipsing stars (i.e., radius and temperature).

The light curves show two clear regions where independent spots are affecting the light curve. This prompted us to adopt a two-spot model. The apparent dip in brightness due the spots is greatest at the primary eclipse causing us to place one spot on the side of the primary star that faces the secondary star (i.e., defined as longitude 0∘ in PHOEBE). In addition, to model the region around the secondary eclipse, we placed a second spot on the primary star 135∘ in longitude from the first spot. These longitudes remained fixed for the rest of the light curve modeling. We iterated on the spot positions with the LM solver in PHOEBE and found equally good fits for these positions within $\pm$2∘ in the longitude, so we adopted these values exactly for our spot longitudes.

The duration of the spot modulation around the primary eclipse covers a large fraction of the total orbital phase of the light curve. To fit the large feature that encompasses from about phases $-$0.3 to 0.3 (Figs. 3 and 4), we could adopt a larger spot at the equator or a smaller spot at a latitude that is closer to the pole. In order to fit this feature, we needed to choose a relatively large spot at a non-equatorial latitude to create the large duration feature. As mentioned above, the size of the spot and its latitude are somewhat degenerate, so many combinations of parameters result in equally adequate fits to any given part of the light curve. Therefore, to model the large feature, we adopted an angular radius of 30∘ and a latitude of 45∘. In addition, the light curve is best fitted when the spot is not occulted by the secondary star—again causing us to choose a non-equatorial spot.

The light curves are well fitted by two spots, both located on the primary stellar surface. The placement and sizes of the two spots on the stellar surface were optimized to match the shape of the asymmetries in the observed light curves about both eclipses. Because the spot temperature is highly degenerate with the spot size (e.g., Gómez Maqueo Chew et al. 2009), the temperature factor (i.e., spot to stellar surface temperature ratio) of each spot was fitted for any given level of third light, because for a fixed size it determines the amplitude of the effect in the light curves due to spots. The best-fit spot parameters for the adopted values of third light (see §3.4) are given in Table 6.

Other observed evidence that supports magnetic activity and thus the presence of surface spots are: the measured activity indicators in the stellar spectra (e.g., H-$\alpha$ is measured in emission; Hebb et al. 2010, and references therein), the blue-excess in the level of third light (§3.5; e.g., Gómez Maqueo Chew et al. 2012; Gillen et al. 2017), and the observed rotational modulation in the WASP light curves.

3.2.2 Effect of third light on eclipsing component properties

From the previous analyses, we are certain that there is dilution in the light curves because of the presence of an unresolved third component of MML 53, which at the distance of UCL (140 $\pm$ 2 pc; de Zeeuw et al. 1999) is a few tens of milli-arcseconds in angular separation from the EB (see §3.3 and Schaefer et al. 2018). In this section, we explore the effects of the third light level on the physical properties of the eclipsing components by modeling the eclipsing binary with varying levels of third light. The largest differences in the EB physical parameters come from the comparison of (a) the case in which there is no dilution in the light curves (i.e., the third light represents 0% of the total light of the system), and (b) the case in which the inclination of the EB is 90∘(i.e., the model eclipses are deepest and thus the dilution by the third light has to be the highest to match the observed eclipse depth). We modeled these two cases with PHOEBE and show our results in Table 7. The temperature factor of each of the two spots (see §3.2.1) was modified to fit the observed amplitude in the light curves, depending on the level of third light. These two extreme cases show that the largest uncertainty in the masses of the eclipsing components due to the level of third light is $\sim$3%. In the case of the radii, this uncertainty is $\sim$1%. We also find that in the case of maximum dilution, the level of third light required to match the observed eclipse depths does not significantly decreases toward the bluer bands as would be expected for a lower-mass tertiary, as is suggested by the height of the peaks of the CCF of the combined spectra and as is determined in §3.3. In fact, the third light level required is relatively flat in all passbands, which indicates an excess in the level of third light in the bluer bands in the case of a lower-mass tertiary.

We consider the uncertainty on the physical properties of the eclipsing stars due to the amount of third light to be much smaller than these values because: (a) the fit to the observed light curves at the two extremes is worse than our best fit model (best $\chi^{2}_{\rm reduced}\approx$1.5 and Table 7); (b) both cases are not physical, because we know that the tertiary exists and it is diluting the light curves, and it is a lower-mass star gravitationally bound to the system; and (c) we do have constraints on the kind of star that is diluting the light curves (see §3.3), even if the amount of light we adopt is model dependent.

3.3 Binary-tertiary model

The tertiary is visible in the spectrum as an independent component, but variations in its radial velocity compared to the systemic radial velocity of the binary confirm that the tertiary is part of a triple system. In addition, the timing of the eclipses of the binary components vary periodically due to the binary’s orbital motion around the tertiary star (Hebb et al. 2010). We used the all of the available light curves to investigate eclipse timing variations caused by the wide tertiary star.

3.3.1 Eclipse timings

The rectified light curves derived from the 2006-2013 WASP data and the two individual eclipse events obtained with the FTS telescope in 2011 are used to measure variations in the time of minimum light for the MML 53 EB due to its motion around a common center of mass with the tertiary companion.

We fitted the rectified light curves from individual seasons of data (including all observed primary and secondary eclipses) using the fast, analytic EB modeling code EBOP (Popper & Etzel 1981; Southworth et al. 2007). This program treats the stars as detached, nearly spherical geometric shapes in order to derive the orbital parameters (i.e., period, epoch, eccentricity) of the binary and some eclipse parameters that are directly related to the shape of the light curve (i.e., sum of the stellar radii, surface brightness ratio), but it does not provide direct physical properties of the stars (i.e., individual temperatures and stellar radii). Since we require only the time of minimum light and not a full EB model solution from these data, this program is suitable for the analysis.

While analyzing the light curves with EBOP, we used the Levenberg-Marquart fitting option and allowed the time of minimum light to be a free parameter, but the orbital period of the binary, the stellar masses, the eccentricity, and the secondary parameters (limb darkening, gravity brightening, reflection coefficients, and third light contribution) remained fixed to the values derived in the final EB model from PHOEBE. The relative sum of the stellar radii, the surface brightness ratio, and the inclination angle were allowed to vary in each case in order to provide sufficient freedom to find the best fitting model light curve while accounting for small variations in the relative eclipse depths due to the filter, the contribution of third light to that filter, and starspots.

The best fitting model light curves are overplotted in red on the phase-folded, rectified light curves shown in Figs. 1 and 2, and the final minimum eclipse times derived from these fits are reported in Table 8. The 2013 WASP data consists of three independent light curves from different WASP cameras that cover the same time period. Each light curve was fitted independently with EBOP, and the weighted average of the three eclipse time measurements is reported in the table. The uncertainty on the 2013 measurement is the standard error of the mean of these three values, but all other reported uncertainties come directly from of the final EBOP results. All final epoch times are converted to BJDTDB time units using the routines provided online by Jason Eastman (Eastman et al. 2010). We note that the WASP data obtained on this target in the summer of 2012 is not shown in the figure or used in the analysis because no primary or secondary eclipse minima were observed during that time period.

3.3.2 Orbital solution of binary-tertiary system

We combine the eclipsing timing epochs in Table 8 with the radial velocity measurements described in §3.1 to constrain the orbital parameters of the binary-tertiary system. We also incorporate into the fits the recent measurements described in Schaefer et al. (2018) of the angular separation of the binary and tertiary along with the angular position change of the components in the plane of the sky. To do this, we developed a program that solves for the orbital parameters of a two-body Keplerian system by treating the MML 53 eclipsing binary as a single mass, $M_{B}=M_{1}+M_{2}$, in orbit with the tertiary star, $M_{3}$, around a common center of mass. This is justified since the separation between the binary components is less than 1% of the separation between the binary and the tertiary. Incorporating the effect of both binary components on the motion of the tertiary is needlessly complicated given the final uncertainties on the orbital parameters of greatest interest.

There are seven independent parameters that are used to define the orbital motion the binary-tertiary system: the orbital period, $P_{3}$; the mass ratio, $q_{3}=M_{3}/(M_{B}$); the eccentricity, $e_{3}$; the argument of periastron, $\omega_{3}$; the tilt of the orbital plane from the observer’s line-of-sight, $\theta_{3}$; the systemic velocity of the triple system, $\gamma$; and the time of periastron, $t_{\rm peri}$. In addition, we assume the mass of the binary, $M_{B}=1.9307$ M⊙ is known, and we used Kepler’s Law to derive the orbital separation, $a_{3}$, between the binary and the tertiary components. These parameters were used in the model generating engine of our program to produce synthetic three-dimensional velocities and positions as a function of time for the binary and tertiary components. In order to identify the optimum values of these parameters that best reproduce the observations, the model generating engine was wrapped by an affine invariant Markov chain Monte Carlo (MCMC) sampler that was integrated into the program itself (Foreman-Mackey et al. 2013).

The program allows the MCMC algorithm to explore the parameter space from a random starting point for all parameters. For each MCMC trial, it uses the known $M_{B}$ along with the adopted $P_{3}$, $e_{3}$, $q_{3}$, and $t_{\rm peri}$ values for that trial to find the true anomaly of the binary and tertiary components in their eccentric orbit. The true anomaly combined with the $\gamma$ velocity, and the orbital separation, $a_{3}$, were used to determine the positions and velocities of the binary and tertiary masses in the ellipse reference frame (with the orbital angular momentum axis pointed perpendicular to the observers line of sight). The ellipse frame was then rotated by $\omega_{3}$ and tilted down from the z-axis by $\theta_{3}$ to achieve the final three-dimensional positions and velocities of the components at each time of an observation. The $\omega_{3}$ was measured from the center of mass to the tertiary star orbit, as for visual binaries. To calculate the $\chi^{2}$ value of that trial set of parameters, the program incorporates a comparison between the synthetic line-of-sight velocity of both components to the observed radial velocities at each time. It also uses the line-of-sight position of the binary component relative to the center of mass divided by the speed of light to compare to the observed eclipse timing offsets. Furthermore, the distance between the EB and tertiary in the plane of the sky was compared to the measured angular separation multiplied by the Gaia (data release 2) distance of 130.2 pc and added to the $\chi^{2}$ value. Lastly, to incorporate the angular motion of the system around the common center of mass in the plane of the sky, which is reported in Schaefer et al. (2018), we first solved for the angular offset which minimizes the difference between the measured angles and those in the model. This is necessary because the orientation of the model angles must match the reference point defined in the observations. This angular offset corresponds to the longitude of the ascending node, $\Omega_{\rm ascending}$, which is the angle from the reference direction north to the line connecting the center of mass and the orbital plane of the tertiary component when it crosses the plane of the sky in the direction away from the observer. This angular offset was then applied to the model angles in the plane of the sky before comparing them to the measurements. For completeness, we also report the angular semi-major axis between the visual components, $a_{\rm angular}$. According to our model, the maximum separation occurs close to the first Schaefer et al. (2018) measurement obtained in 2014.

We ran the model generating engine in the affine invariant MCMC sampler with 1000 walkers. In applying this algorithm, we took special care with certain parameters. Due to the degeneracies in the system, we only allowed $\theta_{3}$ to vary from $0-\pi/2$. We adopt $\sqrt{e_{3}}\cos{\omega_{3}}$ and $\sqrt{e_{3}}\sin{\omega_{3}}$ as sampling parameters and converted these values to $e_{3}$ and $\omega_{3}$ in order to avoid the Lucy-Sweeney bias (Lucy & Sweeney 1971; Eastman et al. 2013). The mass ratio, $q_{3}$ can vary outside of the range from $0.0-1.0$ when updating its value at each step, so we rejected all values greater than 1.0 after the MCMC is complete since these values are unphysical in our model and in our understanding of the MML 53 system.

We allowed each walker to run for 30,000 trial steps of which 30-33% were accepted resulting in a final distribution of 9000-10,000 accepted steps per walker. After examining the output MCMC file of accepted parameters, we chose to remove the first 300 accepted steps from each walker as it constitutes the burn-in phase. Furthermore, only 836 of the initial 1000 walkers converge to a single solution at the global minimum of the $\chi^{2}$ space. The remaining walkers appear to get stuck in local minima at much higher $\chi^{2}$ values with much lower acceptance rates. We removed these walkers from the final distribution, but considered the remaining walkers to be converged (shown in Fig. 8). The final best fitting model has a $\chi^{2}_{min}=26.88$. The best fitting parameters and their $1\sigma$ uncertainties are shown in Table 9. The mass of the tertiary, $M_{3}=0.72^{+0.16}_{-0.09}$ M⊙, is derived by multiplying $M_{B}=1.9307$ M⊙and the newly derived mass ratio $q_{3}=0.37^{+0.08}_{-0.05}$. In Fig. 7, we show the best fitting model radial velocity curves of the EB and third star system, and eclipse timing curves compared to the observations.

3.4 Third light determination

Light from the unresolved tertiary star causes both the primary and secondary eclipses of the EB to appear shallower than they should. This directly affects the derived orbital inclination angle of the EB, which indirectly influences the individual masses and radii. Thus, in order to derive accurate fundamental properties of the primary and secondary components from the eclipsing light curve, quantitative values of the light contributed by the tertitary star in each filter must be incorporated into the EB model. The third light values applied here were derived from the latest theoretical stellar evolution models from Baraffe et al. (2015) based on the mass and age of the tertiary star. These isochrones provide VRCIC band absolute magnitudes as a function of mass and age that are interpolated to find the flux contribution of the third star relative to primary and secondary components in each filter.

An initial guess for the third light was derived from the relative height of the cross-correlation peaks in the 2006 FEROS spectrum. This third light value is used to perform a preliminary EB model which results in masses and luminosities for the primary and secondary components as described in §3.2. This preliminary value of third light does not affect the final values of third light used nor the EB physical properties. To derive the third light used in the final EB model, the primary and secondary star masses were combined with the binary-tertiary mass ratio, $q_{3}$, to find the mass of the tertiary star. With preliminary masses for all three components known, we implemented a quadratic interpolation of the published isochrones that are less than $50$ Myrs old at each of the component masses to derive a series of theoretical VRCICband absolute magnitude values for each star as a function of age.

At this point, we could interpolate the mass tracks at the independent age measured for the Upper Centaurus Lupus cluster ($\sim 16$ Myr) to find the relative luminosities of the stellar components, however we chose instead to find self-consistent values that match both the luminosities derived from the EB model and a single theoretical isochrone for all three stars. The EB model provides a measured value for the primary-to-secondary flux ratio in each filter based on the temperatures and radii of each star. We used this value as a constraint and convert the theoretical absolute magnitude values into an array of primary to secondary star flux ratios ($L_{2}/L_{1}$) in the VRCICbands as a function of age. We then interpolated the models at the measured $L_{2}/L_{1}$ value from the EB model in each filter.

This provides an independent age estimate for the system, which is the same for all three filters within 1 Myr. Finally, we interpolated the VRCICband absolute magnitude values for each component at that age and calculate the tertiary star’s relative contribution to the overall light of the system in each filter.

This is an iterative process in which the EB model, the binary-tertiary model, and the third light calculation are performed in consecutive order until the masses and relative light contributions of the three components have converged within the uncertainties. The relative VRCIC band fluxes determined for each star after several iterations are shown below including the third light contribution ($L_{3}/(L_{1}+L_{2}+L_{3}$)) that is used in the final EB model. As a consistency check, we measured the relative surface brightness from the CCF of the FEROS spectrum fitting a three-Gaussian model obtaining 0.59:0.23:0.18, which are fully consistent with the values presented in Table 10.

3.5 From EB model: masses, inclination, sum of the radii and temperature ratio

Utilizing the eclipsing binary tool PHOEBE (Prša & Zwitter 2005), we modeled both the UVES radial velocity measurements for the primary and secondary components, and the CTIO light curves. We only fitted the VRCIC light curves to derive the physical properties of the eclipsing components to limit the uncertainty introduced by the level of third light (3.4). The VRCIC light curves are shown in Fig. 3 and compared to the best-fit model described in this section. Not only are the UB bands not included in the evolutionary models that determine the amount of expected dilution due to the tertiary, but we find that the U-band has an additional third light contribution than would be expected for a less massive star. Thus, the UB bands were not used to determine the best-fit model, but are shown for reference in Fig. 4, and are included in this paper to distribute the full CTIO photometric dataset to the community.

Adopting the spot sizes and placements from §3.2.1, the level of third light from §3.4, and the $T_{\mathrm{eff,1}}$ (= 4880 K) from the spectral disentangling (§3.1), we randomly sampled 80,000 times the following parameters (and ranges): orbital inclination (79.5 $\leq$ $i$ $\leq$ 84.1∘); the primary potential (5.9 $\leq$ $\Omega_{1}$ $\leq$ 10.0); the secondary potential (5.7 $\leq$ $\Omega_{2}$ $\leq$ 9.7), and the secondary temperature (4245 $\leq$ T_{\mathrm{eff,2}}$$\leq 4545 K). For each combination of these parameters, we fitted the temperature factor of each stellar spot and the overall luminosity for each light curve, interpolated the limb-darkening coefficients for each band, and calculated the $\chi^{2}$ of each model.

From the resulting multidimensional $\chi^{2}$-space and considering the detached and circular orbit of the eclipsing binary, we obtained confidence levels (shown in Fig. 9) for the properties that are derived directly from the light curves, namely: the inclination angle, $i$; the temperature ratio, $T_{\mathrm{eff,2}}$/$T_{\mathrm{eff,1}}$from the relative depth of the eclipses, and the sum of the fractional radii, $r_{1}+r_{2}$ = ($R_{\mathrm{1}}$+$R_{\mathrm{2}}$)/a from the duration of the eclipses. In the case of MML 53 because the eclipsing binary orbit is circular and the eclipses are V-shaped, we are able to constrain the sum of the fractional radii from the light curves and not the radius ratio. Given this degeneracy, although we sampled the primary and secondary potentials for the parameter-space exploration, we do not report the individual values as they are not significant, as is well known for EBs in circular orbits (Kallrath & Milone 2009). The best-fit model to the VRCIC light curves (see Fig. 3) was identified for having the lowest $\chi^{2}$, with a corresponding reduced-$\chi^{2}\approx$ 1.5.

In the case of the U and B light curve models, we adopted the VRCIC best-fit solution and fitted the level of third light to match the amplitude of the variation due to spots and the depth of the eclipses. We find that the B-band is well fitted with a dilution by the tertiary of 8% of the total luminosity of the three-body system, while the U-band requires a 15% dilution. Given that the tertiary is less massive (and thus, redder) than the eclipsing components, we would expect the level of dilution due to the tertiary to decrease toward bluer wavelengths. An additional blue component in the third light levels could be due to accretion on to the tertiary, as has been identified in at least another PMS eclipsing binary Par 1802 (Gómez Maqueo Chew et al. 2012). Similarly, a $u$-band excess has been shown by the eclipsing components of CoRoT 223992193 (Gillen et al. 2017). Other young, single stars have been observed to have UV excess and optical spectra accretion features (e.g., Findeisen & Hillenbrand 2010). Figure 4 shows the model and observed light curves in the U and B passbands. We have not included the residuals to the models because these light curves are not used to derive the best-fit model.

3.6 Derivation of the semi-major axis, the individual radii and the secondary temperature

Based on the resulting parameters and their associated uncertainties from the EB model to both the RVs and VRCIC light curves (§3.5), we calculated the physical properties of interest, namely the semi-major axis, the primary and secondary radius, and the effective temperature of the secondary component. All values are summarised in Table 11.

To determine the physical scale of the orbit, we utilized the parameters derived directly from the RV curves ($a\sin{i}$ and $q$) and their formal uncertainties to the fit together with the values and uncertainties from the confidence levels of the quantities that depend solely from the light curves to derive the physical properties of the eclipsing components, their orbit, and their corresponding uncertainties. Specifically, we derived the semi-major axis of the eclipsing orbit from the definition of $a\sin{i}$ and the measured $i$. Once $a$ was determined, we derived the individual masses and the total mass from $q$, orbital period and $a$ through the equations of Keplerian motion.

Given the circular orbit of the MML 53 EB, the grazing nature of its eclipses and the contamination of the photometry by the third star, we require an external constraint in order to derive the individual radii of the eclipsing components (e.g., Kopal 1959; Stassun et al. 2004, 2014). However, in the case of MML 53, the flux ratio derived from the CCF is uncertain, and thus, it was not utilized in this analysis to derive individual radii as the external constraint. The primary radius was instead determined from the measurement of the $\upsilon\sin i$ of the primary component from the LSD analysis of the high-resolution UVES spectra (§3.1), the inclination $i$ from the EB model (§3.5), and that the primary star is synchronized and its spin-axis is aligned with its orbital motion (§3.2), given that by definition $\upsilon\sin i$1 = $2\pi\leavevmode\nobreak\ R_{\mathrm{1}}\leavevmode\nobreak\ \sin{i}/P_{\rm rot,1}$.

The secondary radius was derived from the sum of the fractional radii (Fig. 9, bottom panel), $a$, and the primary radius. With the individual radii and masses of the eclipsing components, the surface gravities for the eclipsing components ($\log g_{\mathrm{1}}$ and $\log g_{\mathrm{2}}$) were derived readily utilizing the fundamental constants from Prša et al. (2016). As a consistency check, we calculated the secondary radius from the LSD measurement of its $\upsilon\sin i$ and find it to be in agreement within the uncertainties with the secondary radius reported in Table 11.

The secondary temperature was derived from the spectroscopically-determined primary temperature $T_{\mathrm{eff,1}}$ (§3.1), and the temperature ratio resulting from the EB model (Fig. 9, top panel). We present our measurements for the physical properties of the eclipsing stars derived from the best-fit eclipsing binary model, adjusting both the RVs and VRCIC light curves in Table 11.

4 Summary and discussion

MML 53 is a gravitationally bound hierarchical triple system where all three components are in the pre-main sequence. It consists of a close eclipsing binary composed of a 1.0400 M⊙ primary star and a $0.8907$ M⊙ secondary star in a $P_{\rm EB}\sim 2.09$ day orbit, and a distant, lower mass ($\sim$0.7 M⊙) star in a longer period ($\sim$8.5 yr) orbit. The masses of the eclipsing components have been determined with $<$1% precision, which allows the mass of the tertiary to be measured to $\sim$20%. Additionally, our analysis of the EB allows us to measure the radii of its components to be 1.237 and 1.153 R⊙ for the primary and secondary stars with a precision of 2.7% and 3.5%, respectively. We also measure the individual temperature of the primary from the spectral analysis to be 4880 $\pm$ 100 K (2% precision), and derived from the temperature ratio that of the secondary star to be 4380 $\pm$ 100 K (2% precision). Although MML 53 shows all the axes of complexity for pre-main sequence EBs (higher multiplicity, spots, possible accretion), we are able to measure precisely the individual properties of the eclipsing stars and constrain for the first time the mass of the tertiary.

Furthermore, our analysis permits the even more precise, direct measurement of two EB quantities: (1) the sum of the fractional radii (0.9% precision) from the duration of the eclipses, and (2) the temperature ratio (0.2% precision) from the relative depth of the eclipses in each passband. These two quantities derived from our eclipsing binary model are robust measurements, and are independent from the light contamination by the third star. Thus, the comparison of these two direct measurements to those predicted by theoretical evolutionary models is important. Their evolution predicted by the Baraffe et al. (2015) models is shown in Fig. 10 and describes the evolution of the eclipsing stars. The primary has reached the Henyey track (Henyey et al. 1965), heating up and slowing down its contraction; whereas the secondary star continues contracting at roughly the same temperature along the Hayashi track (Hayashi 1961). The errors in the temperature ratio and sum of the radii predicted by the models come from the uncertainty in the measured masses (Table 11); the errors in the measured temperature ratio and sum of the radii were derived from the contour maps from the EB modeling (Fig. 9). Assuming, as is standard for non-interacting, close binaries, that the eclipsing stars are coeval, Fig. 10 shows that that the two ages derived independently from the theoretical evolutionary models compared to our measurements of the temperature ratio ($\sim$17 Myr; top panel) and the sum of the radii ($\sim$9 Myr; bottom panel) are not in mutual agreement.

In-detail analyses (Pecaut & Mamajek 2016; Pecaut et al. 2012) present the mean age of the UCL subgroup (and of MML 53 itself) to be $16\pm 2$ Myr based on a comparison between 14 high-mass, turn-off stars to the theoretical models for rotating stars from Ekström et al. (2012), and on F-type and G-type members compared to pre-main sequence theoretical stellar evolution models (Dotter et al. 2008; Baraffe et al. 2015; Tognelli et al. 2011; Chen et al. 2014). An observed large spread in the ages of individual members of the subgroup is thought to be partly due to observational uncertainties and unresolved multiplicity, but also to an intrinsic spread of ages within the subgroup. MML 53 is in the region of the association that is close to the average age, and is not in a part of the cluster that is thought to be much younger (high galactic longitudes) or older (edges of region). Additionally, the age derived from the temperature ratio (top panel; Fig. 10) is consistent with these age estimates for UCL; thus we adopt an age of $16\pm 2$ Myr for this specific object. Because the radii that we measured for the two eclipsing stars of MML 53 match to a much younger age (bottom panel in Fig. 10), the eclipsing components appear to be inflated.

The presence of a third independent component in the spectrum of MML 53 had been identified from its discovery, as had the changes in the timing of the eclipses of the binary (Hebb et al. 2010). However, it is only recently that the tertiary star has been resolved (Schaefer et al. 2018) at its widest separation from the EB. It is in this paper that we have confirmed that the tertiary is a bound component of the MML 53 system and have constrained its mass dynamically. The $\sim$20% uncertainty in the tertiary mass is a conservative estimate from the binary-tertiary model (§3.3), since the constraint from the optical spectra showing a lower luminosity for the tertiary has not been incorporated in the tertiary mass determination.

After carefully accounting for the third star, its effects on the measurement of the individual properties of the eclipsing components of MML 53 were minimized, allowing for a meaningful comparison between these measurements, the theoretical evolutionary models of Baraffe et al. (2015) and other direct measurements of pre-main sequence stars in double-lined, eclipsing binaries (see Figs. 11 and 12). The radii of both eclipsing components of MML 53 are consistent with a single isochrone. However, they match the younger 10 Myr track (Fig. 11). The radius of the primary star is inflated by 10% with respect to the radius predicted by stellar models interpolated at 16 Myrs for a 1.04 M⊙ star, and the secondary radius is inflated by 15% as predicted for a 0.89 M⊙ star.

At 16 $\pm$ 2 Myr, MML 53 is very similar in age to the EB NP Per ($\sim$17 Myr; Lacy et al. 2016), and importantly, the primary star of MML 53 and the secondary of NP Per have the same mass of 1.04 M⊙, and as shown in the mass–radius diagram (Fig. 11), both stars have the same radius, measured independently, that appears to be too large, as predicted by the Baraffe models. More generally, the stars of a given eclipsing system, except for Par 1802 (blue, downward triangles at 0.4 M⊙) and NP Per (green squares), also fall on the same isochrone on the mass–radius diagram, showcasing that the theoretical models can describe the overall behavior of young, low-mass stars. However, the direct measurements of the EBs masses and radii appear to be younger than the ages derived by independent methods (e.g., from the study of the young associations: turn-off stars, mass of arrival at the main-sequence, lithium abundance). The apparent radius inflation of these young stars (including the eclipsing components of MML 53) could be due to physical processes not included in the evolutionary tracks, for example, the inhibition of convection due to magnetic fields that slow down contraction (Feiden 2016; MacDonald & Mullan 2017), causing the lowest-mass young stars to appear younger.

Stars that are younger than 10 Myr and have masses lower than 1.0 M⊙ have been at approximately the same temperature throughout their evolution, as they contract along the Hayashi track. This behavior is evident in the overlap of the theoretical isochrones and EB measurements shown in the mass–temperature diagram (Fig. 12). In the case of MML 53 at 16 Myr and with masses straddling 1.0 M⊙, its primary component has begun increasing in temperature, while for the lower-mass secondary the tracks predict a small change in temperature ($\sim$100 K, within our errors) from 1 to 16 Myr. Importantly, our measurements of the individual temperatures follow the slope of the 16 Myr isochrone, even if slightly above it. Both stars are slightly hotter than predicted by the 16 Myr track. It is puzzling that the MML 53 eclipsing stars do not appear to be cooler, as would be expected if the inflated radii were due to magnetic activity, by inhibition of convection and/or spots (e.g., Feiden 2016; Somers & Pinsonneault 2015), in order for the stars to remain in thermal equilibrium. In the mass–temperature plane, the primary of MML 53 and the secondary component of NP Per are not consistent within one sigma, falling above and below the 16 Myr Baraffe model, respectively. Their individual temperatures would be consistent with each other, and with the 16 Myr model within two sigma. It is not surprising that in general the scatter in the mass–temperature diagram is larger than in the mass–radius diagram, as individual temperatures are harder to measure from the combined spectra of (at least) two stars.

Some studies of triple systems composed of a close binary bound to an outer third star suggest that the formation of close binaries may be the result of tidal tightening of the inner binary due to the third star (e.g., Tokovinin et al. 2006; Naoz & Fabrycky 2014). However at 16 Myr, the primary and secondary stars of MML 53 are already in a close binary ($\sim$0.03 au). Our results of MML 53 are more in agreement with recent population synthesis models suggesting that the mechanisms causing the tightening of the close binary orbit occur in most cases early in the stars’ evolution (Moe & Kratter 2018).

Acknowledgements.

YGMC was supported in part by UNAM-PAPIIT IN-107518. LHH was supported by grants NSF-AAG-1009810 and NSF-AAG-1312453. We would like to thank G. Feiden, R. Barnes and D. Fleming for useful discussions. Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO program ID 383.C-080. This research has made use of the services of the ESO Science Archive Facility. This work makes use of observations from the LCOGT network. Based on observations at Cerro Tololo Inter-American Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. The CTIO observations made use of 1m and 1.5m telescopes operated by the SMARTS Consortium. WASP-South is hosted by the South African Astronomical Observatory and we are grateful for their ongoing support and assistance. Funding for WASP comes from consortium universities and from the UK’s Science and Technology Facilities Council.