{'fraction_non_alphanumeric': 0.03266525532714718, 'fraction_numerical': 0.033575371465880535, 'mean_word_length': 4.3282732447817835, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The dependence of convective core overshooting on stellar mass: reality check, and additional evidence

Overshooting from the convective cores of stars more massive than about 1.2 M(Sun) has a profound impact on their subsequent evolution. And yet, the formulation of the overshooting mechanism in current stellar evolution models has a free parameter (f[ov] in the diffusive approximation) that remains poorly constrained by observations, affecting the determination of astrophysically important quantities such as stellar ages. In an earlier series of papers we assembled a sample of 37 well-measured detached eclipsing binaries to calibrate the dependence of f[ov] on stellar mass, showing that it increases sharply up to a mass of roughly 2 M(Sun), and remains constant thereafter out to at least 4.4 M(Sun). Recent claims have challenged the utility of eclipsing binaries for this purpose, on the basis that the uncertainties in f[ov] from the model fits are typically too large to be useful, casting doubt on a dependence of overshooting on mass. Here we reexamine those claims and show them to be too pessimistic, mainly because they did not account for all available constraints --- both observational and theoretical --- in assessing the true uncertainties. We also take the opportunity to add semi-empirical f[ov] determinations for 14 additional binaries to our previous sample, and to update the values for 8 others. All are consistent with, and strengthen our previous conclusions, supporting a dependence of f[ov] on mass that is now based on estimates for a total of 51 binary systems (102 stars).


INTRODUCTION
In recent years there has been a resurgence of interest in the phenomenon of convective core overshooting in stars. Empirical constraints of different kinds have been brought to bear on the problem of calibrating the free parameters of the two most often used prescriptions for overshooting, which are α ov for the classical stepfunction implementation, and f ov for the diffusive approximation (Freytag et al. 1996;Herwig et al. 1997). In the classical formulation the extension of the core beyond the boundary set by the Schwarzschild criterion is d ov = α ov H p , whereas in the alternate prescription the mixing is modeled as a diffusive process with a diffusion coefficient at a radial distance r from the boundary given by D(r) = D 0 exp (−2r/f ov H p ), in which D 0 is the coefficient inside the boundary and H p the pressure scale height. In a series of papers over the last three years (Claret & Torres 2016, hereafter Paper I, Paper II, and Paper III) we employed a total of 37 double-lined eclipsing binaries (DLEBs) with wellmeasured masses, radii, effective temperatures, and in some cases metallicities, to examine the dependence of the strength of overshooting on stellar mass. The measurements were compared against current stellar evolution models to infer semi-empirical values for the overshooting parameter for each star.
It was found that α ov and f ov increase sharply starting at ∼1.2 M ⊙ to a maximum at ∼2 M ⊙ , after which they show little or no change up to the 4.4 M ⊙ limit of our sample. The peak value for f ov is about 0.016. The mass behavior was found to be independent of the element mixtures we tested (Grevesse & Sauval 1998;Asplund et al. 2009), and was qualitatively the same for α ov and f ov , except for a scaling factor between the two of α ov /f ov ≈ 11.4 (Paper II).
In recent work Constantino & Baraffe (2018) (hereafter CB2018) have revisited the issue with the goal of assessing the degree to which observations for DLEBs constrain the overshooting parameter. They focused specifically on f ov , and selected a subsample of eight of our binaries for their experiments, which they considered to be representative of the range of stellar masses and evolutionary states of the parent population. Their modeling led them to claim that in most cases the range of f ov values permitted by the observations is very large (often the full interval they explored, typically f ov = 0.000-0.040), such that it is very difficult to discern any dependence at all on stellar mass. They concluded therefore that the method of using DLEBs for this type of calibration is of limited utility. Overall they found that a constant value of the overshooting parameter around 0.013 or 0.014 seems adequate to fit all eight of the systems they studied, given that they have such large uncertainties. Values of f ov of this order would be consistent with the degree of overshooting we found in our previous work for masses larger than about 2 M ⊙ (Paper II, Paper III), but not for lower mass stars.
Intrigued by the CB2018 claims that seem to conflict with our own experience, as a first goal we set out here to examine their results more closely, and to repeat their experiments following their grid-based approach to the extent possible. We depart from their basic procedure only in that we use the same stellar evolution models and physical ingredients that we employed in our previ-ous series of papers, to ensure consistency with our earlier results. To our surprise we find significant disagreements with CB2018 for most of the systems, which in some cases may be explained by the differences in the models. Furthermore, we believe their assessments are overly pessimistic in part because they did not consider other available constraints, some empirical and some rooted in theory, that significantly reduce the formal uncertainties in f ov for some systems. In particular, we review physical arguments that set a maximum size for the convective cores of stars in the low-and intermediate-mass regime ( 2 M ⊙ ), effectively making values of f ov as large as some of those proposed by CB2018 unrealistic. We believe these arguments, along with our own grid experiments for the eight binaries studied by CB2018, show that in many cases DLEB observations do indeed have enough discriminating power to discern a variation of f ov with mass.
As a second goal of this paper we take the opportunity to report semi-empirical determinations of f ov for an additional set of eclipsing binaries we have identified that are suitable for this type of analysis. As we will show, those results are consistent with and strengthen the general trend found earlier.
We have organized our paper as follows. In Section 2 we summarize our procedures for exploring the range of acceptable values of f ov for the systems studied by CB2018. Descriptions of these grid experiments are presented and discussed individually in Section 3 along with a comparison against the results CB2018, highlighting the disagreements. Our conclusions regarding the claims by CB2018 are given in Section 4. In Section 5 we expand our previous sample of 37 binaries and determine semi-empirical values of f ov for another 14 DLEBs. We also revise the values for 8 of them, based on improvements to their stellar parameters that have appeared recently in the literature. This section also reports our updated f ov vs. mass relation in a new rendering that more fairly displays the uncertainties in our determinations. We conclude with a discussion in Section 6. Appendix A describes in some detail important physical arguments suggesting that the degree of overshooting cannot be as large as claimed by CB2018 for stars of low and intermediate mass.

GRID SEARCH METHODOLOGY
CB2018 based their study on eight DLEBs listed in their Table 1. Three have component masses smaller than 2 M ⊙ and are therefore of particular interest (AY Cam, HD 187669, and BK Peg), and the other five are more massive. In most cases the primary and secondary stars are fairly similar in terms of their mass. For each binary they explored a range of values of the overshooting parameter f ov , the standard mixing-length parameter α MLT , and the system metallicity [Fe/H], all of which they report in the same table. They assumed both α MLT and f ov to be identical for the two components. The stellar evolution models they used are based on the MONSTAR code (Campbell & Lattanzio 2008), with physical ingredients described therein. In particular, these models result in a mixing length parameter for the Sun of α MLT⊙ = 1.60, somewhat lower than ours (see below).
To maintain consistency with the results of Paper II and Paper III, our stellar evolution tracks used to fit the binaries were calculated here with the Modules for Experiments in Stellar Astrophysics package (MESA; Paxton et al. 2011Paxton et al. , 2013Paxton et al. , 2015 version 7385, with the physical ingredients given in our previous work, and ignoring the effects of rotation. The solar-calibrated value of the mixing length parameter for these models and our setup is α MLT⊙ = 1.84. All our calculations have included microscopic diffusion. CB2018 did not explicitly say whether they took this effect into account, and we are unable to determine this from an examination of the literature sources for their stellar evolution code. We return to this issue later (Subsection 3.9). The element mixture adopted here is that of Asplund et al. (2009), which is the same as used by CB2018. The corresponding mass fraction of metals for the Sun is Z ⊙ = 0.0134. The helium abundance in our model grids follows the enrichment law Y = 0.249 + 1.67Z, as in our earlier work. CB2018 used fixed values of Y = 0.25 or 0.26.
In general we have chosen to explore the same ranges in f ov , α MLT , and [Fe/H] as did CB2018, though in a few instances we expanded them somewhat. Our grids for each system are uniform and complete in f ov and α MLT . CB2018 did not specify the step sizes they used for f ov or α MLT , and from their description it does not appear that their grids are uniformly spaced or complete. Instead, in most cases we presume they sampled the range manually at a few representative values of f ov and α MLT , as well as [Fe/H]. For the metallicity we typically explored values within the range allowed by the measured abundances, when available, or else we used the intervals listed by CB2018.
We note also that while three of the binaries in their Table 1 show a range for the mixing length parameter, the other five appear to have been treated with a single value, usually α MLT⊙ = 1.60, which is the mixing length for the Sun according to the MONSTAR models of CB2018. Holding α MLT fixed at a value appropriate for the Sun carries the risk of biasing the fits for stars that are significantly different from the Sun. Particularly in those cases, but also in others, we have expanded the range of α MLT values to account for the theoretical expectation, based on 3D simulations, that α MLT may depend on the evolutionary state of the star or on its metallicity (see, e.g., Magic et al. 2015). 3 Our ranges typically extend toward larger values than CB2018, based on the difference seen in the mixing length parameter for the Sun between MESA and MONSTAR. The larger range of α MLT values we explore will tend to allow a larger spread of f ov values. Doing this is therefore more conservative when it comes to assigning uncertainties to the overshooting parameter.
CB2018 considered a fit to be acceptable when the model predictions for the two components computed at their measured masses are consistent with the measured radii and temperatures within their reported uncertainties. We used the same criterion here. As in our previous work, to be conservative we accounted for possible errors in the models themselves by allowing an age difference of up to 5% between the components, which again has the effect of allowing a larger range in both f ov and α MLT . CB2018 did not explicitly say whether they did something similar, or whether they required the two stars to have precisely the same age in all cases. A further "common sense" rule we have adopted concerns systems for which the observations admit solutions with one of the components in a rapid phase of evolution, and also solutions in a slower phase with a different f ov and/or α MLT . As the star can only be in one evolutionary state at a time, in these situations we have chosen the slower phase because that scenario is more likely a priori. CB2018 allowed both types of solutions at the same time, which results in larger ranges of permitted f ov values.
As we will see below, for several systems it is also the case that f ov values for the primary are better constrained than the secondary because it is a more evolved star. When the mass ratio is close to unity, it is not unreasonable to expect the two stars will also have similar degrees of overshooting, which can effectively constrain the otherwise large range that the secondary on its own would permit. In the next section we show how these reasonable assumptions can be applied to some of the CB2018 systems.

RESULTS OF OUR GRID EXPERIMENTS
We examine each of the eight DLEBs studied by CB2018 in order of decreasing primary mass. The figures accompanying each system below are all drawn to the same scale to permit a direct comparison. A red diagonal line is shown to indicate the f ov intervals that CB2018 claim are allowable, which are always identical for the two components. Note that while a first glance at their Table 1 may suggest that all primary/secondary f ov combinations within their reported ranges give acceptable fits, strictly speaking they have only explored the diagonal of those 2D regions, as shown in the figures that follow. For our own determinations we display the 2D regions to illustrate the sometimes weaker constraint on f ov for the less evolved secondaries, and to identify which combinations of unequal f ov values provide acceptable fits.

OGLE-LMC-ECL-CEP-0227
Our grid of evolutionary tracks for this system covers the f ov and α MLT intervals 0.010-0.020 and 1.80-2.20, respectively, with step sizes of 0.002 and 0.10. Both ranges extend beyond those considered by CB2018, who used a single α MLT value of 2.00. To our knowledge there is no available spectroscopic metallicity for this binary. CB2018 adopted [Fe/H] = −1.00, for which we find no acceptable solutions with MESA that accomodate both component radii and both temperatures within the observational errors, for ages differing by no more than 5%. We do find good fits for the metal fraction of Z = 0.0022 used by Claret & Torres (2017), corresponding to [Fe/H] = −0.78 with the Asplund et al. (2009) element mixture, so we adopted that value here. The acceptable fits are shown in Figure 1a (dots and crosshatched area), and restrict f ov more than indicated by CB2018 despite the freedom we allowed for α MLT and the fact that we did not constrain the primary and secondary values of f ov to be identical. As we cannot make a direct comparison using the same metallicity they adopted, we have at least made it for the same α MLT Figure 1. (a) Values of the overshooting parameter for the primary and secondary of OGLE-LMC-ECL-CEP-0227 that give acceptable fits to the observations within our 5% limit on the age difference (dots). The full ranges we find are shown by the crosshatched area. The red diagonal line in this and subsequent figures marks the full ranges in fov reported by CB2018. Labels indicate the adopted system metallicity and the range of values explored in our grids for fov and α MLT . (b) Same as panel (a), restricted to α MLT = 2.00. value they chose (2.00). Under these conditions we find a single satisfactory solution, indicated in the right panel.

LMC-562.05−9009
For this binary our grids were computed over the f ov range of 0.008-0.016, with a step size of 0.002, and a range in α MLT of 2.20-3.00, in intervals of 0.10. Again these are both slightly larger than those of CB2018. This is another system lacking a spectroscopic metallicity measurement, so we adopted a value of [Fe/H] = −0.73 (Z = 0.0025), close to the CB2018 choice. The solutions reported by CB2018 for this binary are quite similar to ours, and are shown in Figure 2a. One significant difference, however, is that they allowed values of α MLT as high as 3.00, which exceeds the maximum values predicted by the models of Magic et al. (2015), based on full 3D radiative hydrodynamic calculations. In the right panel of the figure we have restricted α MLT to more moderate values between 2.20 and 2.40, which results in a smaller range of allowed values of f ov than reported by CB2018.

χ 2 Hya
This particularly interesting binary system features rather unequal masses (3.605 ± 0.078 M ⊙ and 2.632 ± 0.049 M ⊙ ; Torres et al. 2010). CB2018 used the same f ov values for the primary and secondary. Our models for each star were computed over the same f ov range as theirs (0.000-0.050) in steps of 0.005, and the α MLT values we explored are between 1.60 and 2.20 with our usual step size of 0.10. CB2018 considered a single mixing length parameter of 1.60 for both components (the MONSTAR α MLT⊙ ). No metallicity determination is available for χ 2 Hya, so we investigated the [Fe/H] range explored by CB2018 between −0.15 and the solar composition. The MESA models give no acceptable solutions at either extreme, but we do find good fits for an intermediate value [Fe/H] = −0.05, shown in Figure 3a. At primary f ov values of 0.035 the model predictions are barely within observational errors for that star, and only for mixing length parameters α MLT ≥ 2.00, which are considerably higher than the value adopted by CB2018. For f ov ≤ 0.010 the primary is in a rapid evolutionary phase; we disregard those solutions in favor of those at slower, more likely phases that are permitted at higher values of f ov . The full range of values that the observations allow for the primary and secondary is indicated by the crosshatched region.
If we consider only the single mixing length value adopted by CB2018 (1.60), the acceptable solutions span smaller ranges in f ov for both components (Figure 3b), indicating stronger constraints on overshooting than reported by CB2018 under similar conditions. However, given that our solar value of α MLT is considerably larger than theirs, for a more fair comparison with CB2018 we show in this latter figure the good fits for α MLT = 1.80 (square symbols), which is the closest in our grid to the MESA solar-calibrated value of α MLT⊙ = 1.84. We find slightly more extended ranges in f ov for both stars than with α MLT = 1.60, but still appreciably smaller than indicated by CB2018.
χ 2 Hya is special in that it has a well measured temper- ature ratio of T eff,2 /T eff,1 = 0.945 ± 0.009 from the work of Clausen & Nordström (1978), which sets a strong additional constraint. This quantity is closely related to the difference in the eclipse depths, and is usually much more accurate than the individual absolute temperatures. We find that adding the requirement that the predicted temperature ratio be within the observational error restricts f ov even further. This is shown in Figure 3c for grids over the same α MLT range as before, and in Figure 3d for fixed α MLT values of 1.60 and 1.80.
Regarding the secondary component of χ 2 Hya, it is worth pointing out that even though models formally allow them, overshooting parameters as low as zero or 0.005 are probably unrealistic for a star of this mass (M = 2.632 M ⊙ ). To our knowledge there is no empirical evidence of such small values of f ov for similar stars (see, e.g., Schröder et al. 1997;Stancliffe et al. 2015;Torres et al. 2015;Moravveji et al. 2015Moravveji et al. , 2016Valle et al. 2017). With this consideration, the uncertainty in f ov would be reduced even further.

OGLE-LMC-ECL-26122
For this system our grids cover a broader interval in f ov than explored by CB2018 (0.005-0.024, step size = 0.002), and the same range in α MLT (1.90-2.20, step size = 0.10). Pietrzyński et al. (2013) have reported a spectroscopic metallicity estimate for the binary of [Fe/H] = −0.15 ± 0.10. It is worth noting, however, that CB2018 chose to use considerably lower values reaching [Fe/H] = −0.50. 4 Restricting ourselves to the measured range only, we find no satisfactory solutions at the upper limit of [Fe/H] = −0.05, but do find them at the lower limit, [Fe/H] = −0.25. These fits are shown in Figure 4 with filled circles, and indicate the observations provide very strong constraints on the primary value of f ov , which can only run between 0.020 to 0.024. For the secondary component, values smaller than about 0.014 place it in a rapid phase of evolution that we consider less likely than allowable slower phases corresponding to a higher degree of overshooting. These more likely fits are indicated by the crosshatched area in the fig-4 CB2018 stated that this value "is only marginally more metalpoor than the best fit from Claret & Torres (2017), Z = 0.0050". The value we reported in Paper II is actually Z = 0.0070, corresponding to [Fe/H] = −0.28, which is much closer to the lower limit of the measured value. It is possible that CB2018 mistakenly quoted the Z value from our earlier Paper I (0.0050), which is based on a different element mixture. That lower value was found in Paper II to be biased too low due to the use of an outdated value for the primordial helium abundance.
ure. Once again the relatively similar masses of the components (3.593 ± 0.055 M ⊙ and 3.411 ± 0.047 M ⊙ ; Pietrzyński et al. 2013) would enable us to restrict the secondary interval even further by requiring its f ov to be similar to that of the primary, as CB2018 assumed.
In conclusion, we find that the range of f ov values permitted by the observations of OGLE-LMC-ECL-26122 is considerably smaller than proposed by CB2018.

SZ Cen
The grids were computed for f ov values between 0.0125 and 0.0300 in steps of 0.0025, and α MLT values between 1.60 and 2.20 every 0.10. CB2018 used only their solarcalibrated mixing length parameter of α MLT⊙ = 1.60, and considered a very narrow metallicity range between [Fe/H] = −0.25 and −0.20. There is no measured metallicity available for this binary. The acceptable solutions we find for [Fe/H] = −0.25 are shown in Figure 5a, where we include those with f ov = 0.030 for the secondary even though they barely satisfy the observations within observational errors (and only for α MLT well above 1.60). Selecting only the fits with α MLT = 1.60, to match what CB2018 did, results in a reduced range of permitted f ov values (Figure 5b, circles), which is smaller than reported by CB2018. Solutions with α MLT = 1.80, close to our own solar-calibrated value, are shown with squares. Hatched regions cover the full range for the primary and secondary in each case.
SZ Cen is another system with a highly precise measurement of the temperature ratio: T eff,2 /T eff,1 = 1.035± 0.003 (Grønbech et al. 1977). If we apply this additional observational constraint under the same conditions as before, we obtain the viable solutions shown in Figures 5c (full range of α MLT ) and 5d (restricted α MLT ). In the latter panel there are no acceptable fits with α MLT = 1.60, and only a single one with α MLT = 1.80, showing how much better f ov can be determined.
The results for the other end of the metallicity range, [Fe/H] = −0.20, are shown in Figure 5e-5h using the same parameters as in Figure 5a-5d. The allowed f ov ranges are quite similar for the two compositions.

AY Cam
Our grids of stellar evolution tracks span the full interval of f ov values considered by CB2018 (0.000-0.040), and were calculated with a step of 0.005. Our range in α MLT is 1.60-2.00 in steps of 0.10, whereas CB2018 used only the fixed value α MLT⊙ = 1.60. AY Cam lacks a spectroscopic metallicity estimate, and the [Fe/H] range considered by CB2018 is fairly small. We find no acceptable fits at their lower limit (solar composition) despite scanning a larger range in α MLT , but do find solutions for the upper end, [Fe/H] = +0.10. These fits may be visualized in the left panel of Figure 6, in which each acceptable solution is marked with a dot and may correspond to a different mixing length parameter within the limits we explored. Solutions with f ov < 0.005 place the primary in a rapid stage of evolution, so we have chosen to accept only values of 0.005 or larger, corresponding to slower phases. The full ranges for the primary and secondary are indicated by the crosshatched region.
The right panel of Figure 6 displays the allowed f ov ranges when we hold the mixing length parameter fixed at the value of α MLT⊙ = 1.60 used by CB2018 for both components. This shrinks the allowed region somewhat, as expected, which now becomes considerably smaller than reported by CB2018, particularly for the primary. As before, the permitted interval for the secondary is larger than the primary because it is relatively unevolved. As done previously for χ 2 Hya, we considered also solutions with α MLT = 1.80, near our value for the Sun. In this case the primary f ov values are confined between 0.010 and 0.020 (squares). The sections of parameter  space allowed by these two choices of α MLT are shown by the hatched areas.
CB2018 constrained the f ov values for the AY Cam components to be the same, although the mass difference is actually non-negligible (1.905 ± 0.040 M ⊙ and 1.709 ± 0.036 M ⊙ , q = 0.8972 ± 0.0032; Williamon et al. 2004). If we were to apply the same constraint and force the two stars to have the same f ov , the ranges permitted by the observations at a fixed value of α MLT (as per CB2018) would be even smaller. Thus, under similar conditions as CB2018, we find the degree of overshooting to be considerably better defined for this system than indicated by those authors.

HD 187669
In this case we computed models over the full range in f ov specified by CB2018 (0.000-0.040, with a step of 0.005), and allowed α MLT to vary between 1.60 and 2.20 every 0.10. CB2018 held the mixing length parameter fixed at their solar-calibrated value of α MLT⊙ = 1.60 for both stars. We adopted the same metallicity as CB2018, [Fe/H] = −0.25, which is the measured value for this object (He lminiak et al. 2015). The models that match the observations within errors are marked in Figure 7a with circles. We find no viable solutions for the α MLT value adopted by CB2018, but do find them for α MLT = 1.80. These may be seen in Figure 7b.
For this binary star the full range of viable f ov values (in both panels of the figure) is as large as that found by CB2018 for both components, i.e., 0.000-0.040. Note also that our grids show a tendency for the good fits to cluster along the diagonal of the area, consistent with the fact that the two components of HD 187669 have essentially identical masses (1.505±0.004 M ⊙ and 1.504± 0.004 M ⊙ ; He lminiak et al. 2015).

BK Peg
Our grids in f ov , with a fine step size of 0.002, span the same range as CB2018 (0.000-0.040). For α MLT we explored values of 1.20-2.20 in steps of 0.10, which is much wider than the interval of 1.23-1.35 used by CB2018. The measured metallicity of BK Peg is [Fe/H] = −0.12 ± 0.07 (Clausen et al. 2010). CB2018 explored models in the [Fe/H] range between −0.06 and +0.05 that goes beyond the observational errors. We explored an even wider range between −0.13 and +0.05, and find acceptable solutions only for the upper limit, [Fe/H] = +0.05. We display them in Figure 8a, where it can be seen that the f ov values for the primary are better constrained than reported by CB2018. Restricting α MLT for both stars to the small interval 1.23-1.35 results in no viable solutions, except combinations with one star in that range and the other far outside it (α MLT > 1.70). Mixing length parameters as low as those considered by CB2018 are rather unusual. In Figure 8b we show the good fits with α MLT = 1.80 for both stars. The parameter space for f ov becomes marginally smaller than before in the case of the primary star. We point out, finally, that CB2018 constrained the f ov values for the BK Peg components to be identical despite noting that the mass difference is actually non-negligible (1.414 ± 0.007 M ⊙ and 1.257 ± 0.005 M ⊙ , q = 0.889 ± 0.002; Clausen et al. 2010), and that this could imply the "correct" values of f ov should be different if f ov depends on mass. Here we have not made that assumption.

Summary of results
The figures in the preceding subsections (right panels) show that for most systems we obtain allowable ranges in f ov for the primary and/or secondary that are consistently smaller than reported by CB2018, even though we did not constrain the overshooting to be the same for the two stars, as they did. Doing so would shrink our permitted f ov ranges even more (see below). We believe that in some cases (LMC-562.05-9009, BK Peg) the CB2018 intervals for f ov are larger partly because they allowed the mixing length parameters to reach values beyond what is typically considered realistic for similar stars, whereas we find perfectly acceptable fits without the need for such extreme values. In the case of OGLE-LMC-ECL-26122 the metallicities allowed by CB2018 also deviated considerably from the measured composition, whereas we have always kept the models within the observational errors in [Fe/H], when available, or explored the same range as CB2018 otherwise. The systems for which we tend to agree the most are ones in which CB2018 did not fix the mixing length parameter to their solar-calibrated value (OGLE-LMC-ECL-CEP-0227, LMC-562.05-9009, OGLE-LMC-ECL-26122).
While all our calculations have included the effects of microscopic diffusion, it is not clear whether or not CB2018 did the same in their work. The influence of diffusion can be very important, particularly for low and intermediate mass stars. It enhances the heavy element content of the core, increasing its size as the star evolves. Thus, to some degree it plays a similar role as overshooting, and can affect the determination of f ov (see, e.g., Michaud et al. 2004;Deheuvels et al. 2016). If CB2018 did not include diffusion, then there would be a tendency for them to allow higher values of f ov in their fits, such as they report. Additionally, diffusion alters the surface abundances (i.e., the abundances accessible to observation) to a degree that depends on metallicity and evolutionary state (see, e.g., Dotter et al. 2017;Deal et al. 2018). Enforcing the measured (surface) abundances when fitting models while at the same time ignoring diffusion can therefore bias the results, including the inferred f ov and α MLT values.
We emphasize here again that our calculations do not take into account the influence of stellar rotation, which can also lead to extra mixing and is therefore degenerate with the effect of overshooting, to some extent (see, e.g., Ekström et al. 2012). Including rotation in the modeling would tend to result in lower values of f ov . However, we do not expect rotation to change the overall shape of any dependence there may be between the overshooting parameter and mass.

DOUBLE-LINED ECLIPSING BINARIES
If we now restrict f ov (along with α MLT ) to be the same for the two components of each binary to place our results on the same footing as those of CB2018, the f ov intervals become smaller in several cases. Figure 9 displays a comparison of our f ov uncertainties (gray shaded rectangles) with theirs (open rectangles) for all eight DLEBs in their sample. In one case (OGLE-LMC-ECL-CEP-0227) our grid search returns a single acceptable value of f ov . On the other hand, for LMC-562.05-9009, SZ Cen, and HD 187669 we find good fits over essentially the same ranges as CB2018. There are, however, additional considerations that should not be ignored and that rule out some of the fits, making the formal uncertainties in f ov smaller.
One of those considerations is the tight constraint on the temperature ratio that is available for χ 2 Hya and SZ Cen, which CB2018 did not use. The impact of this was discussed in the previous section, and is indicated again in Figure 9 with an arrow for these two systems. It shows that for χ 2 Hya the viable range in f ov is reduced Figure 9. Comparison between the acceptable ranges for the overshooting parameter from CB2018 (white rectangles) and from our own grids (gray rectangles) assuming the same values of fov for the primary and secondary, and under closely similar conditions regarding α MLT . Systems are in the same order as discussed previously, and primary and secondary masses in solar units are indicated on the right. The hatched rectangles for χ 2 Hya and SZ Cen represent the much reduced ranges we obtain by adding the constraint from the temperature ratio (see previous section). For the three lower-mass systems AY Cam, HD 187669, and BK Peg the hatched areas result from considering the maximum plausible values of fov from theory (Roxburgh 1992) (see Appendix A).
to one tenth the size reported by CB2018, and for SZ Cen it is confined to a single point from our grid.
In five of the eight DLEBs studied by CB2018 the authors claimed the observations are formally consistent with values of f ov exceeding 0.025, and even as high as 0.050 in the case of χ 2 Hya. Aside from the fact that there is no credible empirical evidence for such a large degree of overshooting as this for stars in the mass range of this sample (∼1.2-4 M ⊙ ), we believe a more serious issue is that these extreme f ov values run counter to common sense limits grounded in theory. For stars under about 2.0 M ⊙ that have small convective cores it is well known that typical values of the overshooting parameter coupled with a direct dependence between the expansion of the core and the pressure scale height would lead to core sizes that become unrealistically large for such stars because of the sharp increase in the pressure scale height (see, e.g., Wuchterl & Feuchtinger 1998;Woo & Demarque 2001;Pietrinferni et al. 2004;Demarque et al. 2004). A more detailed discussion of this issue and its consequences may be found in the Appendix. Roxburgh (1992) established that for stars with small cores the extension due to overshoot cannot exceed about 18% of the size of the classical core set by the Schwarzschild criterion (see also Woo & Demarque 2001). This maximum core size effectively limits f ov to values much smaller than those allowed by CB2018, and matters the most for the three lowest-mass stars in the sample, AY Cam, HD 187669, and BK Peg. These are precisely the ones that led CB2018 to claim a lack of evidence for a dependence of f ov on mass. Referring back to Figures 6-8 for these three systems, we have highlighted with large open circles the grid points that satisfy the Roxburgh upper limit. The corresponding ranges allowed for f ov have also been transcribed to Figure 9, and are indicated with an arrow. They are between four and eight times smaller than the CB2018 uncertainties.
Thus, based on our grid experiments and the discussion above, we feel that the general conclusion reached by CB2018 regarding the lack of utility of DLEBs for estimating f ov is overly pessimistic, and perhaps somewhat misleading. From their work it is apparent that their approach to the issue of estimating error bars for f ov was strictly statistical in nature, with no other consideration. The situation is more nuanced, as we have shown, and there are other empirical as well as theoretical constraints that should be taken into account if the uncertainties are to be physically realistic. We have been mindful of these constraints in our previous work, even if not stated explicitly, particularly of the theoretical upper limit on the size of the convective cores for lower-mass stars that seems fairly obvious to us.
However, it is clear -and in this we agree with CB2018 -that the constraint on f ov is much weaker for some of the stars, particularly the secondaries of χ 2 Hya, AY Cam, and BK Peg because they are less evolved, as well as both components in HD 187669 for the reasons mentioned by CB2018. For the unevolved stars this is hardly unexpected, of course, and we have indicated as much in our previous papers. We concur also that while the strength of the constraint on f ov can vary significantly from system to system, there is value in an approach such as ours that draws on a large sample of binaries, adding together whatever information each system is able to contribute for the purpose of investigating how the overshooting parameter may depend on mass.
We point out also that in our previous studies the typical error bars adopted for the f ov estimates from the eclipsing binaries were 0.003 for giants and 0.004 for dwarfs, which we still believe to be quite reasonable when all constraints are taken into account, as described above. Nevertheless, based on the much more detailed examination of the eight systems in the present study as well as additional experiments with new binaries presented below, we consider it prudent to adjust our earlier error estimates to 0.004 for giants and 0.006 for unevolved stars so as to be more conservative. We adopt these uncertainties also in the following.

NEW fov DETERMINATIONS
Since our most recent study the properties of several of the DLEBs studied in Paper II have been updated by Graczyk et al. (2018), in some cases significantly. Additionally, we have identified 13 more well-measured detached binary systems (mass and radius uncertainties less than ∼2%) in which one or both components are sufficiently evolved to be useful for estimating f ov . Most have metallicity estimates. Almost all the new binaries are giants, and in five cases one or both stars are less massive than 2 M ⊙ , and are therefore especially valuable for investigating the dependence of f ov on mass. We have also added χ 2 Hya to the list, which we investigated in detail above as part of the CB2018 sample. This system had been dropped from our previous studies because we had been unable to obtain a good fit to the observations at sufficiently similar ages for the two components using a different set of models (see Paper I), but we have now been able to do so with MESA.
The revised physical parameters for the eight DLEBs from Paper II and for the 14 new binaries including χ 2 Hya are collected in Table 1. Each block of the table has the systems sorted in decreasing order of the primary mass. In two of the new systems, KIC 10031808 and KIC 9246715 (He lminiak et al. 2019), the spectroscopic abundance analysis also yielded a measure of [α/Fe], which we have accounted for below in fitting the evolutionary tracks.
Our new determinations of f ov and α MLT were carried out with the same MESA models described in Section 2, following the same procedure as in our previous studies. Briefly, we made use of coarse grids of evolutionary tracks calculated for the measured component masses and spanning a range of f ov and α MLT values. These grids served to guide subsequent manual adjustments toward the final estimates of the overshooting and mixing length parameters for each star, making use of all observational and theoretical constraints. The values reported here are those giving the smallest chi-squared value, with the main observables typically being the absolute radii and temperatures. The ages for the two components in each system were allowed to differ by no more than 5%. Adopted uncertainties in f ov are 0.004 for giants and 0.006 for dwarfs (see above). The α MLT values have typical errors of 0.20. For full details we refer the reader to our earlier studies (e.g., Section 3 of Paper II, or Section 3 of Paper III). The results are presented in Table 2, along with the best-fit heavy element abundance Z and the mean age of each system.
With the present determinations along with those reported in Paper II (updated here for 8 cases; see Table 1) and Paper III we now have f ov estimates for a total of 51 DLEBs (102 stars), all based on the same MESA models with the same element mixture (Asplund et al. 2009).
The new and revised f ov values are consistent with the results from our earlier studies, and support the existence of a dependence with mass. A representation of this relation that provides a sense of the overall uncertainties is shown in Figure 10 in the form of a "heatmap", in which each f ov measurement was replaced with a distribution of  10,000 random draws from a Gaussian distribution with σ = 0.004 or 0.006 depending on whether the star is evolved or not. To achieve smoothness in the horizontal direction we also perturbed the masses in the same way with σ = 10% of their individual values. Each pixel in the diagram is colored according to the number of random draws it contains. The actual f ov measurements are shown as well.

DISCUSSION
The challenging problem of calibrating the convective core overshooting parameter has been approached in different ways over the last two or three decades using a variety of empirical data. These include color-magnitude diagrams of open clusters and associations, binary stars, rotational velocities of massive stars, and more recently also asteroseismic constraints. Due in large part to ob-servational limitations (measurement errors, sample size, etc.), many of these determinations have reported estimates of a mean or typical value of α ov or f ov for a given population of stars, but have generally lacked the sensitivity to explore changes as a function of mass or other stellar properties. As a result, there is conflicting evidence in the literature: some studies support a variation of the overshooting parameter with mass, while others find that a constant value is perfectly adequate.
For example, Schaller et al. (1992) adopted a fixed overshooting parameter (α ov = 0.2 in the step-function approximation) for stars of all masses, on the basis of a comparison of their models and earlier ones with the morphology of the turnoff of several dozen open clusters (see also Maeder & Mermilliod 1981;Mermilliod & Maeder 1986). Demarque et al. (2004) also used clusters, but  Ekström et al. (2012) and Georgy et al. (2013) from an investigation of the width of the terminalage main-sequence for stars in the range 1.35-9 M ⊙ . On the other hand, in their recent paper presenting a large grid of MESA models Choi et al. (2016) took a different view and adopted a fixed value of f ov = 0.016 for stars of all masses, determined from isochrone fits to the colormagnitude diagram of M67. On the asteroseismic side Aerts (2015) reported seeing no significant trend of the strength of overshooting with mass for a set of mostly OB stars, whereas the study of Deheuvels et al. (2016) did find a hint of a possible increase with mass over the range ∼1.1-1.5 M ⊙ , although their sample was small. While it seems undisputed that overshooting is required in most (if not all) stars above about 2 M ⊙ (Moravveji et al. 2015(Moravveji et al. , 2016Valle et al. 2017, and others), the evidence is just as compelling that other less massive stars such as the components of the eclipsing binary AI Phe (M = 1.23 and 1.19 M ⊙ ) can only be fit properly without overshooting (Higl & Weiss 2017). This imples some sort of transition in the values of f ov or α ov between high-and low-mass stars. In fact, in a sign of a general tendency we see in the literature toward acceptance of such a transition, many current series of stellar evolution models (but not all) already incorporate a dependence of overshooting on mass (e.g., Demarque et al. 2004;Pietrinferni et al. 2004;VandenBerg et al. 2006;Mowlavi et al. 2012;Bressan et al. 2012;Georgy et al. 2013;Spada et al. 2017;Hidalgo et al. 2018), although the adopted functional form of the trend has so far been largely arbitrary.
Early studies based on binaries (Schröder et al. 1997;Ribas et al. 2000;Claret 2007) also suggested a variation of overshooting with mass, but were generally afflicted by small samples and/or large uncertainties. Other more recent binary studies such as those of Meng & Zhang (2014) and Stancliffe et al. (2015) have not found clear signs of a change, but also dealt with few systems. The present investigation, which continues and expands our previous work in Paper II and Paper III, uses by far the largest sample of binaries (51) with accurately measured parameters and presents strong evidence ( Figure 10) for a gradual change in f ov as a function of mass, with a break point near 2 M ⊙ .
Given the history of claims and counterclaims regarding this issue, the study of CB2018 comes as no surprise and is in fact welcome, as it has motivated us to reexamine our own work more critically. They contend, based on their model fits for a subset of 8 representative DLEBs they selected from our previous samples, that the formal uncertainties in the fitted f ov values are so large as to make the more unevolved binaries essentially useless for calibrating overshooting. Our replication of those experiments with the MESA models under the same conditions as theirs, to the extent possible, finds considerably smaller f ov uncertainties in many cases, especially when taking into account the constraint from the temperature ratio available for some systems, which CB2018 ignored. Moreover, the implausibly large f ov values allowed by them in several cases are the result of not taking into account physical arguments that limit the maximum over-shooting distance in stars with small convective cores, specifically, the Roxburgh criterion. This is particularly relevant below ∼2 M ⊙ , which is precisely the regime in which we (and others) have found that f ov varies with mass (see also Figure 13). The purely statistical approach of CB2018 for estimating f ov is of course expected to lead to larger formal uncertainties, but is too simplistic if the goal is to extract overshooting parameters and uncertainties that make physical sense. Their assumption of a single f ov value in binaries with components of appreciably different mass (χ 2 Hya, AY Cam, BK Peg) also appears questionable to us, even more so if the goal is to detect a trend of f ov with mass.
Beyond the notion that there is a maximum permissible extension of the convective core due to overshooting, our current understanding of the phenomenon does not yet allow us to predict the value of the overshooting parameter from first principles with any certainty. Nevertheless, we note that recent theoretical work by Jermyn et al. (2018) has led them to postulate that the parameter depends on mass in much the same way as indicated by the semi-empirical evidence we have presented here. Their Figure 8 for α ov versus mass bears a striking resemblance to our Figure 10 for f ov (once the scale factor α ov /f ov ≈ 11.4 is accounted for; Claret & Torres 2017), except that their peak value of α ov seems about 25% higher than indicated by the observations. Through similar theoretical arguments they find an analogous mass dependence for the diffusive parameter f ov (see their Figure 9), although in this case their scale is a factor of ∼5.5 too large.
We mention in closing that a similar trend of overshooting with mass can be inferred theoretically (and independently of the work by Jermyn et al. 2018) from the homology arguments presented by Claret & Torres (2017). This follows from Equation (7) of their Appendix, which gives an expression for the fractional increase in the mass of the convective core. Details of this derivation will be published elsewhere. This basic physical argument makes it immediately clear that the formal error bars for f ov reported by CB2018 should be reduced significantly at the upper end, and can hardly reach values as large as 0.040 for relatively low-mass stars with small cores. This conclusion is quite aside from the fact that there is no empirical evidence supporting such large values of f ov for stars under about 2 M ⊙ , as pointed out earlier.
Both the step-function and the diffusive formulations of overshooting are largely ad hoc in nature and lack any predictive power to set α ov or f ov from first principles. A valuable attempt in this direction was carried out more than 40 years ago by Roxburgh (1978), who established a criterion to estimate the overshooting distance based on an integral involving the difference between the luminosities due to nuclear and radiative processes. The integrand evaluated at the lower and upper limits of the convective plus overshoot region (r 1 and r 2 , respectively) changes sign at some intermediate radius r ov , permitting a definition of the overshooting distance as r 2 − r ov .
As a device for estimating the overshooting distance itself the Roxburgh criterion was initially called into question by Baker & Kuhfuss (1987), but has since been shown to be valid as an upper limit to that distance (see, e.g., Kuhfuss 1987;Zahn 1991;Roxburgh 1992;Canuto 1997). In particular, Roxburgh (1992) estimated that for stars with small cores the maximum overshooting distance is approximately 0.18 times the size of the classical core, independently of the details of energy generation and opacity. Woo & Demarque (2001) arrived at a similar estimate. For useful discussions about the validity of the Roxburgh criterion, and more generally about the modeling of convection, we refer the reader to the reviews by Maeder (2009) and Kupka & Muthsam (2017).
The Roxburgh (1992) upper limit to the extension of the core is shown in bottom panel of Figure 12 as a dashed line. This naturally sets a maximum value of f ov for stars with small convective cores. We illustrate this more directly in Figure 13, which reproduces our previous estimates of f ov from Claret & Torres (2017. The Roxburgh limit is seen to be consistent with the semi-empirical determinations of f ov . This upper limit from theory is most relevant for the three binaries in the CB2018 sample with the lowest mass: AY Cam, HD 187699, and BK Peg. For these three systems we have marked in Figures 6-8 with an open circle the solutions with f ov values for the primary and secondary that satisfy the Roxburgh criterion (or slightly exceed it, to be conservative). This constraint reduces the range of allowable values quite significantly, and should not be ignored. Figure 13. Previously published estimates of fov from Paper II and Paper III (37 binary systems) as a function of stellar mass, along with the upper limit (dashed line) expected from the criterion of Roxburgh (1978) for stars with small convective cores (see also the bottom panel of Figure 12). The solid line, taken from Paper III, is an approximate representation of the observed trend drawn to guide the eye.