Accurate Low-Mass Stellar Models of KOI-126
aa r X i v : . [ a s t r o - ph . S R ] S e p Accepted for Publication to ApJL
Preprint typeset using L A TEX style emulateapj v. 11/10/09
ACCURATE LOW-MASS STELLAR MODELS OF KOI-126
Gregory A. Feiden , Brian Chaboyer , and Aaron Dotter Accepted for Publication to ApJL
ABSTRACTThe recent discovery of an eclipsing hierarchical triple system with two low-mass stars in a close orbit(KOI-126) by Carter et al. (2011) appeared to reinforce the evidence that theoretical stellar evolutionmodels are not able to reproduce the observational mass-radius relation for low-mass stars. We presenta set of stellar models for the three stars in the KOI-126 system that show excellent agreement with theobserved radii. This agreement appears to be due to the equation of state implemented by our code.A significant dispersion in the observed mass-radius relation for fully convective stars is demonstrated;indicative of the influence of physics currently not incorporated in standard stellar evolution models.We also predict apsidal motion constants for the two M-dwarf companions. These values should beobservationally determined to within 1% by the end of the
Kepler mission.
Subject headings:
Binaries: eclipsing — Stars: evolution — Stars: low-mass — Starspots INTRODUCTION
Double-lined eclipsing binaries (hereafter DEBs) arepowerful tools for testing stellar evolution models. Awealth of information can be gleaned from observationsof DEBs, including precise masses and radii of the com-ponent stars along with the apsidal motion of the system,if the orbit is sufficiently eccentric (see Torres et al. 2010,for a review). To date, observations of DEB systemswith at least one low-mass component (below 0.8 M ⊙ )have painted a grim picture for stellar evolution models.The radii predicted by models are systematically 5 - 10%smaller than those determined from observations, whilethe effective temperatures derived from models are 5%cooler (e.g., Ribas et al. 2008; Morales et al. 2008, 2009;Torres et al. 2010, and references therein).The aforementioned discrepancies have been at-tributed to the effects of large scale magnetic fields(Ribas 2006; Chabrier et al. 2007), as most of the low-mass DEBs discovered are close binary systems withshort orbital periods (generally less than 3 days). Tidalinteractions act to spin-up the stars as they progress to-wards complete tidal synchronization. Short rotationperiods result from this process, amplifying the stellarmagnetic field strength. Strong magnetic fields inhibitthe efficiency of convective energy transport and increasethe total surface coverage of starspots, effectively low-ering the temperature at the stellar surface. In orderto conserve flux, the stellar radius is forced to inflate(e.g., Gough & Tayler 1966; Mullan & MacDonald 2001;Chabrier et al. 2007).Recently, Carter et al. (2011, henceforth C11) discov-ered a triply eclipsing hierarchical triple system, KOI-126, in the Kepler data set which contains two low-mass stars (KOI-126 B, C). The two low-mass starsare in a tight 1.77 day orbit itself orbiting a moremassive primary star (KOI-126 A) on a fairly eccen-tric path about every 34 days. The authors were ableto derive fundamental parameters for all three stars: Department of Physics and Astronomy, Dartmouth College,6127 Wilder Laboratory, Hanover, NH 03755 Space Telescope Science Institute, 3700 San Martin Dr., Bal-timore, MD 21218 M A = 1.347 ± M ⊙ , R A = 2.0254 ± R ⊙ , M B = 0.2413 ± M ⊙ , R B = 0.2543 ± R ⊙ , M C = 0.2127 ± M ⊙ , and R C = 0.2318 ± R ⊙ . Spectroscopy of the more massive primary indicatedit has an effective temperature of 5875 ±
100 K with asuper-solar metallicity ([Fe/H] = +0.15 ± ± isochrones(Demarque et al. 2004). When C11 compared the ob-served mass-radius relation to 1, 2, 4, and 5 GyrBaraffe et al. (1998, BCAH98) isochrones for [Fe/H] ≤ § § §
4, we discussthe implications of this study with regards to the onlyother well known low-mass eclipsing system, CM Dra. MODELS
We constructed individual stellar evolution modelsfor each KOI-126 component along with a series oftheoretical isochrones using DSEP . The physics in-corporated in the models have been described previ- http://stellar.dartmouth.edu/ ∼ models/ Feiden et al.ously (Chaboyer et al. 2001; Bjork & Chaboyer 2006;Dotter et al. 2007, 2008), but we shall provide a briefsummary.Our models include the effects of helium and heavy el-ement diffusion following the prescription of Thoul et al.(1994), though for fully convective stars, diffusion physicsare unimportant. The opacities utilized by DSEP arethe high-temperature OPAL opacities (Iglesias & Rogers1996) and low-temperature opacities of Ferguson et al.(2005). Surface boundary conditions were defined us-ing the PHOENIX model atmospheres (Hauschildt et al.1999a,b) and were attached to the interior model at T = T eff by interpolating in model atmosphere tables.Attaching the atmospheres to T = T eff makes the valueof the convective mixing length used in the atmospherecode inconsequential (Baraffe et al. 1997).Above 0.8 M ⊙ , DSEP uses a general ideal gasequation of state with the Debye-H¨uckel correction(Chaboyer & Kim 1995). In the low-mass regime, DSEPemploys the FreeEOS in the EOS4 configuration, se-lected for its treatment of arbitrary heavy element abun-dances and its inclusion of the H +2 molecule. FreeEOShas also been shown to be valid for modeling stars moremassive than 0.1 M ⊙ (Irwin 2007). Convective core over-shoot is included using the method of Demarque et al.(2004). Rotation was not considered.The only modification made to the underlying physicsin DSEP is related to the partial inhibition of elementdiffusion. We have introduced turbulent diffusion asdescribed by Richard et al. (2005). Turbulent diffusionmodifies the atomic diffusion coefficient and acts to ex-tend the mixing region below the convection zone. Themagnitude of the turbulent diffusion coefficient is tied toan adjustable reference temperature, T , and varies withdensity via: D T = ωD He ( T ) (cid:20) ρρ ( T ) (cid:21) − (1)where ω characterizes the relative strength of turbulentdiffusion and D He ( T ) and ρ ( T ) are the helium diffusioncoefficient and density at the prescribed reference tem-perature, respectively. Proffitt & Michaud (1991) moti-vate the ρ − dependence in order to reproduce the solarberyllium abundance, which appears to be unchangedover time. Thus, any non-standard mixing in the Sunmust be localized to a narrow region below the solar con-vection zone. We select ω = 400 and leave it fixed, inconcordance with Richard et al. (2005). The referencetemperature used primarily in our models is T = 10 ,which was found to best reproduce the observed abun-dance trends of NGC 6397 (Korn et al. 2007).A solar calibration model was generated to determinethe appropriate initial mass fractions of helium ( Y init )and metals ( Z init ), given the solar heavy element compo-sition of Grevesse & Sauval (1998), as well as to calibratethe convective mixing length, ( α MLT = ℓ/H P ). At thesolar age (4.57 Gyr; Bahcall et al. 2005) we were able toreproduce the solar radius, solar luminosity, radius of theconvective boundary, and ( Z/X ) ⊙ with α MLT = 1.938, Y init = 0.27491, and Z init = 0.01884. All of the modelsutilized in this study were calculated using the solar cal- By Alan Irwin: http://freeeos.sourceforge.net R ad i u s ( R ⊙ ) Age (Gyr)1.347 M ⊙ , [Fe/H] = +0.151.315 M ⊙ , [Fe/H] = +0.231.379 M ⊙ , [Fe/H] = +0.07 Fig. 1.—
Stellar evolution tracks used to constrain the age ofKOI-126 A. The dark band signifies the radius constraints imposedby the observations. ibration as a reference and were assumed to be coeval.Super-solar metallicity models with [Fe/H] = +0.15 weregenerated with Y init = 0.28419 and Z init = 0.02469. RESULTS
Stellar Age
The age of the system was constrained by evolving a1.347 M ⊙ model, with [Fe/H] = +0.15 and solar cali-brated α MLT , and matching the model radius with theobserved radius of KOI-126 A (Figure 1). The age wederive for the system is 4.1 ± Mass-Radius Relation
The primary results of this Letter are demonstratedin Figure 2. DSEP accurately reproduces the ob-served radius for each of the low-mass stars at 4.1Gyr with [Fe/H] = +0.15. We find that for massesM = 0.2410 M ⊙ and M = 0.2130 M ⊙ the pre-dicted radii from DSEP are R = 0.2544 R ⊙ andR = 0.2312 R ⊙ , indicating a relative error between themodel and observed radii of less than 0.3%. At solarmetallicity, the models predict radii approximately 1%smaller than those predicted by the super-solar metallic-ity models. The predicted radii are robust. Artificiallyreducing the mixing length ( α MLT = 1 .
00) and fittingthe atmosphere to a deeper point in the stellar envelope( τ = 100) both produced radius changes under 0.5%.Solar metallicity isochrones from DSEP display radiiapproximately 1% larger than radii predicted byBCAH98. The difference is likely a consequence of theequation of state (EOS) utilized by each group. Whereasthe EOS used by BCAH98 is calculated for a pure hydro-gen/helium plasma, FreeEOS calculates the EOS for anccurate Low-Mass Stellar Models of KOI-126 3 R ad i u s ( R ⊙ ) Mass (M ⊙ ) CM DraKOI-126DSEP [Fe/H] = 0.0DSEP [Fe/H] = +0.15BCAH98 [Fe/H] = 0.0 Fig. 2.—
Mass-radius relationship as defined by CM Dra (red– points) and KOI-126 (black – points). Overlaid are 4.1 Gyrtheoretical isochrones from DSEP with [Fe/H] = 0.0 (black – solid)and [Fe/H] = 0.15 (blue – dash) and from BCAH98 with [Fe/H] =0.0 (red – dash-dot). arbitrary metal abundance. Our models can be more reli-ably calculated above solar metallicity. To test the effectof the EOS, we ran DSEP with the Saumon et al. (1995,SCVH) EOS as an attempt to mimic the BCAH98 modelradii predictions. Our models using the SCVH EOS pro-duced radii within 0.5% of the BCAH98 models at thesame mass and composition, illustrating the importanceof the EOS.
Relative Fluxes
The effective temperature of a star is another mea-sure with which to compare our models. While C11found the effective temperature of KOI-126 A to be5,875 ±
100 K from spectroscopy (which our modelsmatch), they were not able to determine the effectivetemperatures of the two M-dwarfs. However, the dy-namical/photometric model utilized by C11 yielded theflux for each M-dwarf relative to the primary. Themodel results were: f B /f A = (3 . ± . × − and f C /f A = (2 . ± . × − .To compare our stellar models, we derived a syntheticcolor- T eff transformation for the Kepler bandpass usingthe spectral response function provided in the
Kepler In-strument Handbook (Van Cleve & Caldwell 2009). Withmodel fluxes determined using the PHOENIX model at-mospheres, it was possible to generate synthetic colorsfor each of the model components of KOI-126. We found f B /f A = 4 . × − and f C /f A = 3 . × − , a 6 σ and 3 σ difference, respectively. We also derived relativefluxes utilizing the empirical relation given on the Kepler
Guest Observer website . The empirical relation relies onthe SDSS g and r magnitudes to derive the approximate Kepler magnitude. SDSS magnitudes for our modelswere calculated using a synthetic color- T eff transforma-tion (Dotter et al. 2008), allowing for an estimate of the Kepler magnitudes. This semi-empirical transformationyielded: f B /f A = 3 . × − and f C /f A = 2 . × − ,consistent with the relative fluxes derived by C11. See http://keplergo.arc.nasa.gov/CalibrationZeropoint.shtml
In an effort to decrease the relative flux for each M-dwarf in the purely theoretical transformation, we var-ied the parameters of the primary star within the givenobservational constraints and changed the amount ofconvective core overshoot (CCO). When considering thelargest mass (1.379 M ⊙ ), lowest metallicity (+0.07), anda relatively high amount of CCO, the relative fluxes werereduced to within 2 σ of the accepted values. As withthe age determination, the tight observational constraintimposed on the radius of the primary made the range ofrelevant radii insignificant in deriving our results. Apsidal Motion Constant
By the end of the nominal
Kepler mission, C11 pre-dict they will know the apsidal motion constant of thetwo low-mass stars with a relative precision of nearly 1%.From our low-mass stellar models, we were able to pre-dict an apsidal motion constant for each star, quantifyingthe degree to which the mass within the star is centrallyconcentrated. To do this we followed the prescription ofKopal (1978) and solved Radau’s equation with j = 2: r dη dr + 6 ρ ( r ) h ρ i ( η + 1) + η ( η −
1) = 6 (2)where η ( r ) is the value of a particular solution toRadau’s equation related to the stellar deviation fromsphericity, ρ ( r ) is the density of the plasma, and h ρ i isthe average density. We used a 4 th order Runge-Kuttaintegration method to solve Radau’s equation to obtaina particular solution at the surface of each star. Havingdetermined η ( R ), we were able to determine the apsidalmotion constant via: k = 3 − η ( R )4 + 2 η ( R ) (3)As a check on our apsidal motion constant integrator, wegenerated polytropic models characterized by the poly-tropic constant n = 1.0 and 1.5 and compared the resultsof our code with the results of Brooker & Olle (1955).The apsidal motion constants derived from the inte-rior structure of our models for KOI-126 B and C are k = 0.1499 and k = 0.1512, respectively. This suggestsour models are slightly less centrally condensed than ann = 1.5 polytrope which is characterized by k = 0.1433.In fact, the run of density and pressure for our modelsindicate our models are best described by a polytropewith n ∼ k values are negligible. Employing the formulaeof Stothers (1974), we find that rotation affects our k values at the 0.02% level for stars rotating with a periodof 1.7 days. DISCUSSION
Low-mass stars below the fully convective boundaryare a wonderful tool to test basic physics. Their low massaffords theorists a stable, long lived ( > yr) labora-tory with which to test the physics incorporated in themodels. These fully convective stars have relatively sim-ple structures, and uncertainties in the opacities, surfaceboundary conditions, and the treatment of convection Feiden et al. TABLE 1Apsidal motion constants (k ) for fully convective stars with [Fe/H] = -0.50, 0.0, and +0.15. Values are quoted at 7different ages (in Gyr). Mass ( M ⊙ ) [Fe/H] 0.5 1.0 2.0 3.0 4.0 5.0 10.00.20 -0.50 0.1523 0.1523 0.1522 0.1522 0.1522 0.1522 0.1521 · · · · · · +0.15 0.1518 0.1517 0.1517 0.1517 0.1517 0.1517 0.15160.25 -0.50 0.1499 0.1498 0.1498 0.1498 0.1498 0.1497 0.1496 · · · · · · +0.15 0.1496 0.1496 0.1496 0.1495 0.1495 0.1495 0.14930.30 -0.50 0.1482 0.1482 0.1481 0.1481 0.1481 0.1480 0.1479 · · · · · · +0.15 0.1479 0.1479 0.1479 0.1478 0.1478 0.1477 0.1476 have relatively small effects on the predicted propertiesof the models (Dotter 2007, § T eff transformationsare of interest and might not provide a fully accuratetransformation to the observational plane. In the low-mass regime, opacities are complicated by the formationof molecules and the peak of the stellar spectrum is nearthe cut-off of the Kepler response function. Interest-ingly, the semi-empirical transformation yields relativefluxes which are entirely consistent with the photometricmodels of C11. However, we must be cautious with thisresult as the systematic errors are not well constrainedfor the transformation from the SDSS magnitudes to the
Kepler magnitude. It must also be noted that the flux ofthe primary star is sensitive to the details of CCO. Weobserved that increasing the amount of CCO brought thepurely theoretical fluxes closer to the observational val-ues. More investigation will be required to accuratelydiagnose the discrepancies.Finally, we note that the determination of the apsidalconstant will provide a crucial test of our stellar evolutionmodels. In particular, it will test the EOS, which directlydetermines the run of density necessary for the computa-tion of the apsidal motion constant. Morales et al. foundthat using the BCAH98 models, k was approximately0.11 whereas our models predict a larger value of approx-imately 0.15. The difference between these two values isdirectly attributable to the EOS, which determines therun of density within a stellar model. If C11 are able toaccurately derive the apsidal motion constant to within1%, it will provide a stringent benchmark against whichto test the interior physics of low-mass stellar evolutionmodels. SUMMARY
In their discovery paper, C11 reported that the triplyeclipsing hierarchical triple KOI-126 appeared to sup-port the mounting evidence that current standard low-mass stellar models are unable to reproduce the observedmass-radius relation. However, we have generated stellarmodels and theoretical isochrone tracks using the Dart-ccurate Low-Mass Stellar Models of KOI-126 5mouth Stellar Evolution Program and find our modelradii agree with the observations. Combining the KOI-126 measurements with previous observations of the low-mass binary system CM Dra, we find that the dispersionin the observed fully convective mass-radius relation issignificant and stands in contrast to theoretical predic-tions. The fact that CM Dra, a system with sub-solarmetallicity, lies on the super-solar side of the theoreti-cal mass-radius relation is indicative of physics currentlynot incorporated in standard stellar models. We predictthe apsidal motion constant for each of the KOI-126 low-mass stars and find k ≃ Kepler mission. This will provide a crucial test for ourmodels and will provide a stringent constraint againstwhich to test all current and future low-mass standardstellar evolution models.We are grateful for the National Science Foundation(NSF) grant AST-0908345, which supported this work.Gratitude is also given to the anonymous referees whosecomments and suggestions have improved the quality ofthis manuscript. We wish to thank F. Allard for illumi-nating some of the particulars of the PHOENIX atmo-sphere code.mission. This will provide a crucial test for ourmodels and will provide a stringent constraint againstwhich to test all current and future low-mass standardstellar evolution models.We are grateful for the National Science Foundation(NSF) grant AST-0908345, which supported this work.Gratitude is also given to the anonymous referees whosecomments and suggestions have improved the quality ofthis manuscript. We wish to thank F. Allard for illumi-nating some of the particulars of the PHOENIX atmo-sphere code.