The mass discrepancy problem in massive stars: an inferential perspective

Introduction

The determination of stellar masses is fundamental to stellar astrophysics. Stellar masses can be inferred dynamically through spectroscopic and eclipsing binary observations or derived from stellar evolution models — see review by [11]. The resultant spectroscopic, and in the case of binary systems also the dynamical masses, are typically discrepant with the inferred masses from the evolutionary models, resulting in the mass discrepancy problem [18]. In massive stars, this discrepancy has been partially attributed to the underestimation of the convective core mass because of unconstrained convective-boundary mixing (CBM; [32, 22]) processes, in addition to limitations in atmospheric modelling [e.g. 19]. Since massive stars have convective cores and radiative envelopes during the main sequence, CBM can increase the convective core mass and prolong their lifetimes (see [12]). Enhanced CBM has been reported to decrease the mass discrepancy in massive binary systems [51, 19].

Moreover, accounting for other mechanisms, such as rotation [51, 52], has been proposed to reduce the mass discrepancy. [50] argued that the mass discrepancy may yet be solved from a proper calibration of all observed parameters, stimulating a debate on whether the discrepancy reflects the need for improved physics in structure and evolution models or better treatment of measurements and modelling tools. Therefore, accurate measurement of the convective core masses and the efficiency of interior mixing processes of massive stars may hold the key to resolving the mass discrepancy problem.

Name	Component	$M \ (\mathrm{M}_\odot)$	$T_\mathrm{eff} \ (K)$	$\log L \ (\mathrm{L}_\odot)$	$\log g \ ($cm$\,$s$^{-2}$)	Reference
CW Cephei	A	12.951(52)	28300(460)	4.243(29)	4.066(6)
CW Cephei	B	11.877(49)	27431(434)	4.119(28)	4.099(6)
HD 254346	A	15.52(4)	28000*	4.57(4)	3.79(3)
HD 254346	B	10.77(5)	16300(420)	3.48(5)	3.79(3)
HD 74455	A	16.30(3)	29700*	4.36(4)	4.08(3)
HD 74455	B	10.32(4)	23300(450)	3.79(4)	4.14(4)
KX Vel	A	13.68(3)	26000*	4.451(35)	3.738(32)
KX Vel	B	8.28(4)	22100(220)	3.416(59)	4.272(41)
$\varphi$Ori	A	17.69(5)	25600*	3.59(3)	4.06(3)
$\varphi$Ori	B	11.61(5)	19900(420)	3.99(4)	4.28(4)
V453 Cygni	A	13.96(23)	28800(500)	4.666(31)	3.708(18)
V453 Cygni	B	11.10(18)	27700(600)	4.163(39)	4.044(9)
V1295 Sco	A	17.26(2)	30500*	4.724(25)	3.843(7)
V1295 Sco	B	12.66(3)	28100(390)	4.227(33)	4.064(9)

Eclipsing binaries with pulsating components are among the best candidates for the study of the mass discrepancy, as they offer the opportunity to apply either binary dynamical constraints or asteroseismic constraints or both — see review by [28]. Masses from eclipsing binary modelling can be measured to an excellent accuracy of better than 1% [20, 15, 11]. Hence, such systems are benchmarks for measuring the mixing efficiency and offer great potential to resolve the mass discrepancy. For example, eclipsing binary systems have been used to help calibrate CBM in intermediate-mass stars [42, 48, 49, 43], and demonstrate that CBM increases with mass up from about $1 \, \mathrm{M}_{\odot}$ to about $2 \, \mathrm{M}_{\odot}$ and then flattens, as also found from pressure-mode asteroseismology of single stars [9]. [8] studied the core masses for the range $[3,9]\,$M$_\odot$ from gravito-inertial modes and found an increasing trend from their sample of 26 Slowly Pulsating B (SPB) stars. [56] revisited eclipsing binaries in the mass range of $1.25 - 4.2\, \mathrm{M}_{\odot}$ and could not recover the previously reported stellar mass dependence. They discussed that there is substantial uncertainty in the derived mass dependence owing to model degeneracies, and suggested that the errors in previously reported CBM trends in mass could be underestimated. Moreover, the mass of the convective core is strongly dependent on the evolutionary stage, as illustrated from gravity-mode asteroseismology [6, 7]. To tackle the effect of mass on the value of CBM properly, one must thus use methods that are capable of high-precision age-dating via estimation of the central hydrogen mass fraction still left in the burning convective core.

The mass dependence of CBM has also been reported in the literature for stars at the high mass end [e.g. 46], with larger CBM having been reported to partially resolve the mass discrepancy. Using a sample of 21 O-type stars in eclipsing binaries with masses $14-56$ M$_{\odot}$, including those in the Large and Small Magellanic Clouds, [50] obtained good agreement between evolutionary model and dynamical masses, suggesting that the discrepancy might have been resolved. However, [51] studied two O-type eclipsing binaries in the Large Magellanic Cloud and found them to be under-massive by 11% or over-luminous by 0.2 dex. The apparent differences observed in the reported possible agreement (or disagreement) of the dynamical masses with the evolutionary masses could be reflecting possible systematics or methods in constraining the mixing parameters rather than only the missing physics not accounted for by the mixing parameters. In particular, [19] conducted a detailed analysis of CBM to reproduce the dynamical masses of 11 eclipsing binary stars with masses of $4.58-17.07$ M$_{\odot}$ derived from binary modelling based on TESS data. They concluded that enhanced CBM partly resolves the mass discrepancy at a given metallicity, while also showing that a change in metallicity significantly affects the mass of a star and its position in the HR diagram. Further observational mapping of the mass discrepancy is thus warranted (e.g. [47]).

Previous studies have largely focused on improving stellar evolutionary physics through enhanced CBM, rotation effects, or additional mixing processes for many reasons [5], including to reduce the evolutionary mass discrepancy. Comparatively less attention has been paid to the influence of the statistical inference approach itself. In particular, when dynamical masses are incorporated directly into Bayesian likelihood functions, their exceptionally small uncertainties may dominate the inference, making it difficult to distinguish genuine improvements arising from stellar physics from those introduced by the construction of the likelihood. The primary objective of this work is to investigate the mass discrepancy using multiple Bayesian inferential formulations designed to separate methodological effects from physical effects of enhanced CBM to the inferred mass discrepancy, while simultaneously deriving the core properties such as the convective core mass (M$_{\rm cc}$), central hydrogen mass fraction (X$_{\rm c}$) and mixing parameters. To this end, we formulate several complementary likelihood constructions, introduce Bayesian multi-objective model selection [e.g. 62, 63, 61] to balance statistical fit quality against mass-discrepancy reduction, and compare inference under both single-star and coeval binary evolutionary models. The resulting inferential strategies enable us to derive the convective core mass and mixing parameters, as well as to quantify how much of the apparent reduction in the mass discrepancy originates from improved mixing parameters and how much reflects the inferential methodology itself.

We describe the sample in § 2, and the evolutionary models in § 3. The analysis framework is described in § 4. In § 5 we describe the results, and we discuss the results and conclude in § 6.

Sample Description

To evaluate the mass discrepancy, we select a small sample of double-lined spectroscopic (i.e. SB2) eclipsing binaries of spectral type B, which contain at least one $\beta$ Cep component. We specifically focus on binaries with a $\beta$ Cep component as these are common among massive stars (i.e. $M \gtrsim 8$ M$_{\odot}$; [31, 14]). The $\beta$ Cep stars pulsate in low-order pressure or gravity modes that allow for robust asteroseismic inferences [33, 30, 35, 29]. In the current study, we do not exploit any asteroseismic constraints but limit ourselves to the binary observables. This has the advantage that our current work can be applied to non-pulsating binaries as well. We leave the inclusion of stricter additional inferences from identified pulsation modes to future work.

To ensure homogeneity and consistency, only systems with complete published parameters are considered. We assembled a sample of seven such pulsating eclipsing SB2 systems, which have component masses that range between $8-20$ M$_\odot$. Consequently, we have a total of 14 stars in the sample. Table 1 shows the list of systems in our sample, which includes their reported dynamical masses, effective temperature, luminosity, and surface gravity. Any system missing any of these four observational parameters is considered to have incomplete parameters and is not included in the sample.

We note that CW Cephei and V453 Cygni have been previously modelled in specific studies of the mass discrepancy using a $\chi^{2}$ merit function [19]. Our larger sample and more sophisticated statistical inference treatment in the covered mass regime are the major novelties of this work.

Stellar Models

In stellar evolution models, the total time evolution of the mass fraction $X_{i}$ of a chemical element $i$ is given as [23, 22]:

where $r$ is the radius coordinate, $m_{i}$ is the mass of the $i^{\rm th}$ chemical element, R$_{\rm ji}$ and R$_{\rm ik}$ denote the nuclear reaction rates producing and destroying species $i$, respectively, $\rho$ is the local density, $D_\mathrm{mix}$ is the turbulent diffusivity in units of $\text{cm}^2\,\text{s}^{-1}$ [24, 22], and $\mathcal{M}_i^{\text{micro}}$ accounts for extra microscopic atomic diffusion effects such as gravitational settling and radiative levitation (see [3, 4, 13]).

The term $D_\mathrm{mix}$ in Eqn. ((1)) assumes that CBM can be described diffusively and is used to capture the combined contributions from different physical processes such as convection, magnetic fields, rotation, and waves in different regions of the star that are not included otherwise. Hence for a slowly rotating, non-magnetic massive star, the main contribution to $D_\mathrm{mix}$ is from the convective core, the CBM region, and mixing in the radiative envelope. Therefore, $D_\mathrm{mix}$ is defined following [22], such that:

where $D_\mathrm{CBM}(r)$ accounts for the mixing in the convective boundary region. In fast rotating and/or magnetic stars, the mixing contributions from rotation and/or magnetic fields are usually added to Eqn. ((2)) as extra terms. In the case of close binary systems, where binary interaction is important, the effect of tides on mixing could also be captured in $D_\mathrm{mix}(r)$ as an additional term.

There are generally two commonly used prescriptions for diffusive CBM: an exponentially decaying coefficient or a constant one described by a step function. In the case of an exponential description, the mixing efficiency decays exponentially with distance from the convective boundary, with the diffusion profile [24, 26, 22] defined as:

where $f_\mathrm{CBM}$ is the exponential CBM scaling parameter, $D_{0}$ is the convective mixing efficiency at $r_{0}$ and $H_{p,0}$ is the pressure scale height at $r_{0}$. At $r_{\rm cc}$, where the convective boundary occurs, the mixing length theory (MLT; [27]) predicts the convective velocities of rising convective parcels of gas, and hence the associated convective mixing efficiency, to vanish. The mixing length, $l_\mathrm{mix}$, which is the distance a convective fluid element travels before mixing with the surrounding, is typically implemented as a product of the local pressure scale height, $H_{p}$, resulting in $l_\mathrm{mix}=\alpha_{\rm MLT}H_{p}$ [25], such that $\alpha_{\rm MLT}$ is the mixing length parameter. Hence, the exponential CBM is calibrated at $r_0 = r_{\text{cc}} - f_0 H_{p,\text{cc}}$ where $H_{p,\text{cc}}$ is the pressure scale height at the convective boundary and $f_{0}$, in the case of MESA models, is fixed to values between $0.002-0.005$ [22]. The temperature gradient in the overshoot zone is assumed to be the radiative one, given as:

We use this description of convective overshooting in our work.

To determine the properties of the seven pulsating SB2 binary systems in this work, we used precomputed grids of stellar evolutionary models. We used the specific MESA grid calculated by [64], which was computed based on a revised version of [30] using the Modules for Experiment in Stellar Astrophysics [MESA: 37, 38, 39, 40, 41, 36]. The varied parameters include mass, central hydrogen mass fraction (X$_{\rm c}$), CBM using a diffusive exponential overshooting ($f_{\rm CBM}$), additional diffusive mixing in the envelope defined as $D_{\rm mix} = D_{\rm mix,0} \frac{\rho_0}{\rho}$ according to the hydrodynamical simulations by [65], where $\rho$ the local density and $\rho_0$ the density at the base of the envelope. A solar metallicity of $Z=0.014$ was used. The mass ranges from 7 to 30 M$_{\odot}$ in logarithmic steps of $\Delta\log \rm M = 0.015$, $f_{\rm CBM}$ from 0.005 to 0.035 in steps of 0.005, $\log D_{\rm mix,0}$ from 1 to 5 (or 2 to 6 when $M > 13\,$M$_\odot$) in steps of 1, and $X_{\rm c}$ from 0.7 to 0.0001 with the resolution increasing as the $X_{\rm c}$ decreases to ensure the end of the main sequence is well sampled.

In the rest of this paper, we only use this stellar model grid, which was designed specifically for the application of $\beta\,$Cep sample asteroseismology by [64]. This implies that conclusions are only valid for this particular choice of the input physics that went into this grid. Comparisons among grids of models with a variety of input physics in terms of fit quality, as in [2] or [32] are not treated in our work.

Analysis

Each of the 14 stars in our sample of seven pulsating SB2 eclipsing binaries is assumed to have negligible tidal deformation. First, all stars are treated as effectively single stars, such that they are assumed to have evolved independently of their companion. As a second case, we treat both components in each binary as coeval. Comparison of these two cases provides insight of how the mass discrepancy depends on the assumption of independent and coeval evolution for both components of the binaries.

Interpolation framework

The grid of MESA evolutionary models is discrete. As a result, to evaluate the models within the grid as if they were continuous functions, an interpolation scheme was applied. We use a two-step piecewise cubic hermite interpolating polynomial (PCHIP) [57]; first along the evolutionary state ($X_{\rm c}$) in single-star case or age in coeval case, and second across models for fixed values of $f_{\rm CBM}$ and $\log D_{\rm mix,0}$. PCHIP preserves monotonicity and avoids spurious oscillations between grid points [58]. This is particularly relevant for the chosen evolutionary coordinates (i.e. $X_{\rm c}$ and age), which evolve monotonically during a star’s main-sequence within the main sequence lifetime. Hence, its use preserves the local shape of the evolutionary tracks, ensures numerical stability of likelihood evaluations, and prevents interpolation-induced artefacts from influencing the statistical results.

In the single-star case, the interpolation was performed using $X_{\rm c}$ as the evolutionary coordinate. The tracks were sorted in decreasing $X_{\rm c}$ and PCHIP interpolators were built along each track for the parameters $\log T_{\rm eff}$, $\log L$, $\log g$, age, and $M_{\rm cc}$, after removing possible duplicate values of $X_{\rm c}$.

In the coeval binary case, the components of a binary system are assumed to have the same age (and metallicity) since they formed from the same birth cloud. To enforce this coeval behaviour, the interpolation procedure was modified and the stellar age instead of $X_{\rm c}$, was adopted as the evolutionary coordinate. The models were sorted by age and PCHIP interpolators were constructed as functions of age for the parameters $\log T_{\rm eff}$, $\log L$, $\log g$, $X_{\rm c}$, and $M_{\rm cc}$ after removing duplicate values of age. For a proposed common binary age, each track was evaluated at that age. The resulting quantities were then interpolated across stellar mass using PCHIP interpolation, for fixed $\log D_{\rm mix,0}$ and $f_{\rm CBM}$.

Statistical Framework

Using observational constraints derived from the eclipsing binary models, we conducted Bayesian inferential analysis to determine mixing parameter values that best fit the observables and possibly resolve the mass discrepancy. Because covariance matrices for the published binary observables are unavailable, we adopt the standard approximation of independent Gaussian uncertainties for the observables and a simplified Gaussian likelihood function. Thus, the likelihood function, $\mathcal{L}$ based on $\chi^2$ follows

where $D$ are the observed data, $\theta$ represent the model parameters, $\mathcal{M}$ is the adopted stellar evolutionary hypothesis (e.g. single-star evolution, coeval binary evolution, or a specific mixing parameter), and

with $O_i, \sigma_i$ and $M_i$ the observed quantity, its error, and the theoretical prediction for the quantity, respectively. To estimate the free parameters, we use a log-likelihood merit function,

Using Bayes’ theorem, the posterior probability distribution is thus defined as:

where $P(\theta \mid \mathcal{M})$ is the prior and $\mathcal{Z}(D \mid \mathcal{M})$ is the posterior evidence. Uniform priors were adopted for the sampled parameters (mass, X$_{\rm c}$, age) over the ranges explored in the model grid. The mass and $X_{\rm c}$ are sampled in the case of single-star models and the mass and age are sampled in the case of coeval binary models, while $f_{\rm CBM}$ and $\log D_{\rm mix,0}$ are treated as discrete hyper-models¹ explored by outer loops in both cases. Binary observables such as $T_{\rm eff}$, $\log g$, $L$, and mass are considered in the likelihood in a controlled experimental manner described in the likelihood framework given in § 4.3. The analysis was conducted through nested sampling [e.g. 59] with dynesty² [60, 66], which is a Python package for estimating Bayesian posteriors and evidences using nested sampling.

For each combination of mixing parameters ($f_{\rm CBM}$, $\log D_{\rm mix,0}$), nested sampling yields the posterior distribution of the sampled stellar parameters together with the Bayesian evidence, $\ln{\mathcal{Z}}$. The posterior distributions are used to infer stellar properties such as the stellar age, M$_{\rm cc}$, evolutionary mass, and the mass discrepancy ($\Delta M$) between the mass of the evolutionary model and the dynamical mass, while $\ln{\mathcal{Z}}$ provides a quantitative measure of the statistical support for each set of mixing parameters after marginalising over the sampled space of the other free parameters. The models are then compared with each other to identify the most likely one based on their Bayesian evidence.

To ensure numerical convergence of the nested-sampling procedure, preliminary calculations were performed with live points ($N_{\rm live}$)³ of 150, 300, 350, and 500, together with multiple independent random seeds. Convergence was assessed from the stability of both the Bayesian evidence and the posterior-derived quantities relevant to the mass discrepancy analysis. Evidence estimates obtained at successive $N_{\rm live}$ values were consistent within their nested-sampling uncertainties, while the posterior medians and credible intervals showed no systematic changes with increasing $N_{\rm live}$. Therefore, increasing the number of live points beyond $N_{\rm live}=300$ (or 500 for the more complex model formulations) produced no material change in the inferred parameters or model-comparison conclusions. We consequently adopted the lowest $N_{\rm live}$ of 300 (or 500 for the coeval case) for which these quantities were numerically stable, thereby avoiding additional computational expense without a corresponding material gain in inference precision.

Likelihood framework

Using the description in § 4.2, we built seven different models for each star. These models are defined based on the quantities that are included in the likelihood function, the stellar case, and the parameters used for statistical inferences. The different statistical formulations described below represent different statistical experiments because they are based on different sets of observables and conditions. Consequently, Bayesian evidences are compared only among models evaluated under the same likelihood function definition. For example, the reference model (RM) and Pareto model (PM) use only $T_{\rm eff}$, $\log \rm L$ and $\log \rm g$ in the likelihood, whereas the mass model (MM), Pareto mass model (PMM) and coeval cases additionally include the dynamical mass. Since these models condition on different data vectors, the values of their evidence ($\ln \mathcal Z$) are not interpreted as directly comparable across the different likelihood frameworks. Instead, the evidence is used to rank competing mixing parameter values within each likelihood definition.

Reference model (RM):

The reference model is the model with minimal mixing ($f_{\rm CBM}= 0.005$ and $\log D_{\rm mix,0} =1.0$). $T_{\rm eff}$, $\log g$ and $L$ are included in the likelihood function to derive the reference model. This gives the baseline evolutionary mass at minimal mixing and a supposedly underestimated core mass with maximum possible discrepancy.

Mass model (MM):

In this model, the dynamical mass is included in the likelihood function in addition to those in the reference model. Due to the accurate measurement of the dynamical masses, the relative uncertainties are usually very small compared to those of other measured quantities. As a result, the likelihood is usually dominated by the mass and the model with the least $\Delta M$ most often gives the maximum evidence in this case.

Pareto model (PM):

Here, we remove the dynamical mass from the likelihood similar to the RM, but compute models for all the possible mixing parameter values allowed in the grid. Instead of picking the model with the maximum posterior evidence, we treat the selection of the preferred mixing parameter values as a Bayesian multi-objective model-selection problem. For each combination of the mixing parameters, the nested-sampling analysis yields the Bayesian evidence, $\ln \mathcal Z$, which measures the overall statistical support of the stellar model given the observables included in the likelihood. In addition, the posterior samples are used to evaluate the probability that the inferred evolutionary masses of the stellar components are simultaneously consistent with their independently measured dynamical masses within their combined observational and model uncertainties. We denote this quantity by $P_{\rm consistent}$.

The Bayesian evidence and $P_{\rm consistent}$ quantify complementary aspects of model quality. In simpler terms, the evidence measures the statistical agreement between the stellar evolutionary model and the observational constraints, whereas $P_{\rm consistent}$ measures the probability that the model successfully reproduces the dynamical masses of the components. Since the model that fits the observables may not always be the one that reduces the mass discrepancy, we formulate the model-selection problem as a Bayesian multi-objective Pareto-based optimization [e.g. 62, 63, 61].

Since the dynamical masses are excluded from the likelihood function, they play no role in determining either the posterior distribution or the Bayesian evidence. Consequently, the inferred evolutionary masses constitute genuine posterior-predictive quantities, and their agreement with the independently measured dynamical masses provides an external test of the adopted mixing parameters. The quantity $P_{\rm consistent}$ therefore measures the predictive success of each mixing parameter in this case, rather than the quality of a fitted observable.

A model is considered Pareto optimal if no alternative model provides concurrently both stronger statistical support, as demonstrated by the evidence, and better agreement with the dynamical mass [see e.g. 62, 63, 61]. The collection of such models forms the Pareto frontier and represents the set of best compromise solutions between the overall statistical fit quality and simultaneous reproduction of the binary masses. Figure 1 illustrates this Pareto frontier in the two-dimensional objective space ($P_{\rm consistent}$, $\ln \rm \mathcal{Z}$) for the binary CW Cep. Each point represents one set of mixing parameter values explored by the Bayesian analysis.

When a single representative solution is required for presentation or subsequent diagnostic analyses in systems where the Pareto-optimal models are more than one, we adopt a Bayesian decision rule. We first retain all Pareto-optimal models whose evidences satisfy $\Delta \ln \mathcal{Z} = \ln \mathcal{Z}_{\rm max} - \ln \mathcal{Z} \leq 1$ which correspond to models with comparable statistical support. Among these statistically competitive Pareto-optimal models, we select the solution with the largest value of $P_{\rm consistent}$. This procedure avoids introducing an arbitrary scalar merit function or weighting between competing objectives while yielding a representative model for subsequent analysis.

The component-wise mass discrepancies, $\Delta M_{1}$ and $\Delta M_{2}$, are reported separately throughout this work as diagnostic quantities describing the residual evolutionary–dynamical mass offsets of the individual stellar components, but they are not themselves used as optimization objectives. In summary, the PM formulation separates model fitting from model prediction: the atmospheric observables constrain the stellar models through the likelihood function, whereas the dynamical masses provide an independent posterior-predictive assessment of the inferred mixing parameter values.

Pareto mass model (PMM):

The PMM formulation follows the same multi-objective optimisation strategy as the PM model but differs only in the inclusion of the dynamical masses in the likelihood function. Consequently, the Bayesian evidence reflects the statistical support for the complete set of fitted observables, while the residual mass discrepancy is interpreted as the remaining mismatch after fitting the dynamical masses rather than as an independent model prediction. The Pareto-based procedure is therefore used to balance the overall statistical support for the fitted model with the remaining residual mass discrepancy, thereby distinguishing models that achieve similar statistical support but differ in their ability to minimise the fitted mass residual.

Models RM, MM, PM, and PMM are computed for the single star case for each of the targets in the sample. However, the components of eclipsing binaries are reported to be coeval. As a result, enforcing coeval age for both components adds additional constraints and better captures the evolution of the components of eclipsing binaries. In subsequent models built for the targets, the coeval age was enforced and the PMM strategy was adopted.

Coeval fixed mixing model $\text{(Coeval\_FM)}$:

Here, we aim to test whether a given set of mixing parameter values can simultaneously reproduce both components of a given binary system in our sample. We sample the masses of both components and the coeval age. In this model, a set of $f_{\rm CBM}$ and $\log D_{\rm mix,0}$ are obtained that are expected to simultaneously reproduce both components.

Coeval varying mixing model $\text{(Coeval\_VM)}$:

This model is similar to the coeval fixed mixing model, except that the same set of $f_{\rm CBM}$ and $\log D_{\rm mix,0}$ for both components is not enforced. In this case, the mixing parameters for both components are independent and are expected to be independently determined by the respective masses of the components under the coeval assumption. For example, CBM is reported to be dependent on the mass [e.g. 42, 43, 64, 67]. Thus, allowing different mixing parameters for the components clearly captures the effect of mass on the mixing parameters and shows whether enhanced independent mixing for both components is enough to solve the mass discrepancy in binary systems.

Coeval mass-dependent mixing model $\text{(Coeval\_MD)}$:

Here we assume varying mixing parameters for both components like the coeval varying mixing model except that the values of CBM for both components are not independently free, but coupled to each other by a mass-dependent relation. A common envelope mixing is shared by both components in this formulation.

Enhanced mixing has been reported to be required for high-mass stars [e.g. 44, 45, 19] and to partially resolve the mass discrepancy problem [19]. Motivated by the previous observational evidence, and to investigate whether CBM varies systematically with mass in our current sample, we introduced an empirical first-order Taylor-like mass-dependent parameterisation of:

where $f_{\rm CBM}(M)$ is the CBM parameter assigned to a star of mass $M$, and $f_{\rm CBM,ref}$ is the reference value of CBM parameter evaluated at the reference mass $M_{\rm ref}$. Here, $M_{\rm ref}$ is chosen as the midpoint of the binary component masses and given as

We also define

as the local gradient at the reference mass with units of CBM per solar mass ($\mathrm{M}_\odot^{-1}$). A gradient of $s_f=0$ reflects the $\text{coeval\_FM}$ model, where the same mixing parameter values successfully reproduce both components. The choice of reference mass as the mean of the binary masses is to symmetrically parametrise the mixing relation within the binary system.

Because the evolutionary grid is computed only for discrete values of $f_{\rm CBM}$, the inferred values of $f_{\rm CBM} (M)$ are restricted to the available grid points. Consequently, the gradient $s_{\rm f}$ is not interpreted as a precise continuous derivative of CBM with mass. Instead, it is used here as an empirical diagnostic indicating whether the preferred grid-level value of CBM increases, decreases, or remains constant between the two binary components. A finer grid or continuously sampled $f_{\rm CBM}$ would be required to infer a more continuous-like mass dependence.

Results

Influence of Bayesian inferential formulations

Results for each of the different statistical formulations are obtained. As expected, their comparison demonstrates that the derived results are systematically influenced by the adopted likelihood constructions. When the dynamical mass is included in the likelihood (like in the PMM formulation), its strict demand and comparatively small uncertainty exerts a strong influence on the posterior distribution, resulting in a significantly reduced mass discrepancy. In contrast, when the mass is excluded from the likelihood function (like in PM formulation), the inferred solutions are determined by the combined constraints imposed by the effective temperature, luminosity and surface gravity, allowing the residual mass discrepancy to be evaluated independently.

The Pareto-based model selection procedure further shows that the statistically preferred models are not always those that maximize the evidence or minimise the residual mass discrepancy. Instead, the preferred solutions balance overall statistical support demonstrated by the statistical evidence with agreement between evolutionary and dynamical masses demonstrated by the posterior probability consistent with zero discrepancy. Figure 1 illustrates the Pareto trade-off between the Bayesian evidence and discrepancy reduction.

Physical implications for mixing

Across the explored grid of evolution models, increasing the values of the mixing parameters generally increases the convective core mass and reduces the inferred evolutionary mass discrepancy, thus corroborating previous findings in the literature [45, 19]. However, the dependence is not strictly monotonic. Excessively increasing the values of the mixing parameters may also shift the evolutionary tracks away from the observed positions in the Hertzsprung–Russell (HR) diagram or fail to satisfy the remaining observational constraints. Consequently, the Bayesian evidence does not invariably favour the strongest available mixing.

As clearly visible in Figure 2 and Figure 3, which summarise the behaviour of the preferred mixing parameters on the residual mass discrepancy and the corresponding convective-core masses, both CBM and envelope mixing have a significant impact on the measured mass discrepancy and the convective core mass, as already found by [19]. However, the effect of envelope mixing appears to become comparably stronger as its value increases.

In addition, while the inferred convective core masses correlate with larger values of mixing and anti-correlates with the mass discrepancy, stars with similar core masses do not always exhibit identical discrepancies. This indicates that the inferred stellar properties depend jointly on the adopted mixing parameters, the evolutionary state, and the observational constraints entering the inference. An average of $3-9\%$ of the relative convective core mass growth arises from the increase in the values of the mixing parameters for the chosen input physics in the adopted model grid.

Coeval binary inference

The coeval analyses provide an additional physically motivated constraint because both components of each binary system are required to share a common age. This reduces the available parameter space and therefore provides a more stringent test of the adopted mixing parameters. By comparing PMM and Coeval_VM, which share similar likelihood function conditions and independent mixing parameters for both components, we find that enforcing coeval evolution reduces the observed mass discrepancy by a median of only 0.7% and an average of $4\%$ for both components. The additional constraint imposed by the common age barely improves the optimisation procedure. The slight reduction of the mass discrepancy for the adopted input physics of the grid is shown in Fig. Figure 4, where we overplot the Kernel Density Estimate (KDE) of the sparse histogram to place the result into context given our small sample size.

For other coeval models such as Coeval_FM and Coeval_MD, which share slightly different mixing parameterisation with PMM (although with similar likelihood conditions), the comparisons yielded inconclusive results. Comparing the left panel of Figure 3 to Figure 2 demonstrates the impact of assuming coeval evolution on the mass discrepancy for the case of CW Cep.

The coeval_FM formulation reproduces several binary systems with a common mixing parameter value for both components fairly well, whereas other systems require different preferred mixing efficiencies for the two components. It is not clear why some components of comparably similar mass difference between the components may be reproduced by similar mixing parameter values while others are not, as our sample is too small to derive conclusions at population level.

Allowing component-dependent mixing (like in the case of Coeval_VM models) improves the agreement for most binaries in the sample because there are more degrees of freedom. Figure  Figure 5 shows an example of the convective core mass relation with $f_{\rm CBM}$ at the best fit $\log D_{\rm mix,0}$, coeval age, and mass of the Coeval_MD model for CW Cep. The figure illustrates the local gradient for this particular binary when we rely on the chosen input physics of the adopted grid.

Figure 6 illustrates the mass-dependent CBM explored from our Bayesian analyses over-plotted with the preferred CBM for both components in the case of HD 254346. While the figure suggests that more massive components require more efficient mixing, the apparent step-like behaviour is just a consequence of the discrete sampling of $f_{\rm CBM}$ in the underlying MESA grid rather than evidence for a discontinuous physical relation. The values of $s_{\rm f}$ ($0-2\times10^{-3}~M^{-1}_{\rm \odot}$) obtained for the Coeval_MD models are only weakly constrained and should only be interpreted qualitatively.

The derived parameters of the targets are listed in Table Table 2 for the single case model and in Table Table 3 for the coeval binary case model. We emphasise that systems whose components were not successfully reproduced by a coeval model were also not well reproduced in the single star case, pointing towards other physical processes at play that are not captured by the input physics in our MESA model grid.

RM
Star	CW Cep	HD 254346	HD 74455	KX Vel	V453 Cygni	V1295 Sco
$M_{1}\, (\mathrm{M}_{\odot})$	13.42(29)	15.89(43)	14.93(42)	14.13(33)	17.13(37)	18.39(36)
$M_{2}\, (\mathrm{M}_{\odot})$	12.37(21)	-	9.37(21)	7.74(15)	12.67(32)	13.34(23)
$M_\mathrm{cc, 1} \, (\mathrm{M}_{\odot})$	4.07(15)	4.20(21)	4.95(27)	3.25(15)	4.43(17)	5.63(19)
$M_\mathrm{cc, 2} \, (\mathrm{M}_{\odot})$	3.67(11)	-	2.33(13)	2.11(7)	3.62(16)	4.02(12)
$\frac{M_\mathrm{cc, 1}}{M_{1}}$	0.30	0.26	0.33	0.23	0.26	0.31
$\frac{M_\mathrm{cc, 2}}{M_{2}}$	0.29	-	0.25	0.27	0.29	0.30
$X_{\rm c,1}$		0.26(2)	0.53(3)	0.206(22)	0.197(3)	0.316(6)
$X_{\rm c,2}$	0.53(1)	-	0.50(3)	0.655(29)	0.475(9)	0.501(9)
$\mathrm{age}_1\,$(Myr)	6.2(3)	8.5(4)	4.70(69)	10.67(44)	8.37(24)	6.54(18)
$\mathrm{age}_2\,$(Myr)	6.3(3)	-	11.5(1.6)	5.3(2.4)	7.52(41)	6.37(29)
MM
Star	CW Cep	HD 254346	HD 74455	KX Vel	V453 Cygni	V1295 Sco
$M_{1}\, (\mathrm{M}_{\odot})$	12.959(49)	15.25(4)	16.29(3)	13.68(3)	14.40(19)	17.26(2)
$M_{2}\, (\mathrm{M}_{\odot})$	11.89(5)	-	10.29(4)	8.25(4)	11.42(14)	12.67(3)
$M_\mathrm{cc, 1} \, (\mathrm{M}_{\odot})$	4.27(3)	4.59(8)	5.78(15)	3.59(6)	4.89(9)	6.56(1)
$M_\mathrm{cc, 2} \, (\mathrm{M}_{\odot})$	3.82(2)	-	2.75(9)	2.22(8)	3.51(7)	4.16(2)
$\frac{M_\mathrm{cc, 1}}{M_{1}}$	0.33	0.30	0.35	0.26	0.34	0.38
$\frac{M_\mathrm{cc, 2}}{M_{2}}$	0.32	-	0.27	0.27	0.31	0.33
$X_{\rm c,1}$	0.529(5)	0.301(18)	0.556(29)	0.256(17)	0.294(2)	0.383(5)
$X_{\rm c, 2}$	0.553(5)	-	0.523(35)	0.619(31)	0.503(7)	0.528(7)
$f_{\rm CBM, 1}$	0.035	0.005	0.005	0.02	0.035	0.035
$f_{\rm CBM, 2}$	0.035	-	0.005	0.005	0.035	0.035
$\log D_\mathrm{mix, 1}\, (\mathrm{cm}^{2}\,\mathrm{s}^{-1})$	2.0	6.0	3.0	2.0	6.0	6.0
$\log D_\mathrm{mix, 2}\, (\mathrm{cm}^{2}\,\mathrm{s}^{-1})$	5.0	-	3.0	2.0	5.0	5.0
$\mathrm{age}_1\,$(Myr)	7.46(17)	10.11(29)	3.74(61)	11.83(26)	12.96(34)	8.21(8)
$\mathrm{age}_2\,$(Myr)	7.91(23)	-	8.91(1.3)	7.58(3)	10.67(41)	7.98(26)
PM
Star	CW Cep	HD 254346	HD 74455	KX Vel	V453 Cygni	V1295 Sco
$M_{1}\, (\mathrm{M}_{\odot})$	13.13(29)	15.27(46)	14.95(42)	13.65(33)	15.60629	17.32(36)
$M_{2}\, (\mathrm{M}_{\odot})$	12.11(22)	-	9.37(21)	7.74(15)	12.29792	12.98(25)
$M_\mathrm{cc, 1}\, (\mathrm{M}_{\odot})$	4.39(15)	4.80(24)	4.95(28)	3.74(16)	5.48(19)	6.59(20)
$M_\mathrm{cc, 2}\, (\mathrm{M}_{\odot})$	3.93(11)	-	2.33(13)	2.11(7)	3.94(16)	4.32(13)
$\frac{M_\mathrm{cc, 1}}{M_{1}}$	0.33	0.31	0.33	0.27	0.35	0.38
$\frac{M_\mathrm{cc, 2}}{M_{2}}$	0.32	-	0.25	0.27	0.32	0.33
$X_\mathrm{c, 1}$	0.533(5)	0.321(17)	0.534(29)	0.269(18)	0.297(2)	0.383(4)
$X_\mathrm{c, 2}$	0.556(5)	-	0.501(35)	0.655(29)	0.510(8)	0.531(7)
$f_{\rm CBM, 1}$	0.035	0.03	0.005	0.025	0.035	0.035
$f_{\rm CBM, 2}$	0.035	-	0.005	0.005	0.035	0.035
$\log D_\mathrm{mix, 1}\, (\mathrm{cm}^{2}\,\mathrm{s}^{-1})$	5	4	3	3	6	6
$\log D_\mathrm{mix2}\, (\mathrm{cm}^{2}\,\mathrm{s}^{-1})$	5	-	1	1	5	5
age$_{1}$( Myr)	7.36(39)	9.76(50)	4.70(69)	12.08(55)	11.21(48)	8.16(29)
age$_{2}$( Myr)	7.55(38)	-	11.52(1.63)	5.34(2.40)	9.07(55)	7.57(38)
PMM
Star	CW Cep	HD 254346	HD 74455	KX Vel	V453 Cygni	V1295 Sco
$M_{1}\, (\mathrm{M}_{\odot})$	12.96(5)	15.25(4)	16.29(3)	13.68(3)	14.40(19)	17.26(2)
$M_{2}\, (\mathrm{M}_{\odot})$	11.89(5)	-	10.29(4)	8.25(4)	11.42(14)	12.673)
$M_\mathrm{cc, 1}\, (\mathrm{M}_{\odot})$	4.28(3)	4.59(7)	5.78(15)	3.77(6)	4.89(9)	6.56(1)
$M_\mathrm{cc, 2}\, (\mathrm{M}_{\odot})$	3.78(2)	-	2.74(11)	2.22(8)	3.51(7)	4.12(2)
$\frac{M_\mathrm{cc, 1}}{M_{1}}$	0.33	0.30	0.35	0.28	0.34	0.38
$\frac{M_\mathrm{cc, 2}}{M_{2}}$	0.32	-	0.27	0.27	0.31	0.33
$X_\mathrm{c, 1}$	0.529(5)	0.301(19)	0.556(29)	0.269(18)	0.294(2)	0.383(5)
$X_\mathrm{c, 2}$	0.551(5)	-	0.523(36)	0.619(31)	0.503(7)	0.527(8)
$f_{\rm CBM, 1}$	0.035	0.005	0.005	0.025	0.035	0.035
$f_{\rm CBM, 2}$	0.03	-	0.005	0.005	0.035	0.035
$\log D_\mathrm{mix, 1}\, (\mathrm{cm}^{2}\,\mathrm{s}^{-1})$	4	6	3	2	6	6
$\log D_\mathrm{mix, 2}\, (\mathrm{cm}^{2}\,\mathrm{s}^{-1})$	5	-	2	2	5	2
age$_{1}$( Myr)	7.47(18)	10.11(29)	3.74(61)	12.02(28)	12.96(34)	8.21(7)
age$_{2}$( Myr)	7.76(22)	-	8.92(1.4)	7.39(2.15)	10.67(41)	7.84(26)

Coeval_FM
Star	CW Cep	HD 254346	HD 74455	KX Vel	V453 Cygni	V1295 Sco
$M_{1}\, (\mathrm{M}_{\odot})$	12.949(53)	15.24(4)	16.29(3)	13.68(3)	14.79(17)	17.26(2)
$M_{2}\, (\mathrm{M}_{\odot})$	11.901(46)	10.67(5)*	10.30(4)	8.25(4)	11.26(11)	12.67(3)
$M_\mathrm{cc, 1}\, (\mathrm{M}_{\odot})$	4.29(7)	3.69(8)	5.66(13)	3.10(7)	4.60(90	6.18(2)
$M_\mathrm{cc, 2}\, (\mathrm{M}_{\odot})$	3.83(13)	2.79(2)	3.08(4)	2.09(1)	3.42(5)	4.17(1)
$\frac{M_\mathrm{cc, 1}}{M_{1}}$	0.33	0.24	0.35	0.23	0.31	0.36
$\frac{M_\mathrm{cc, 2}}{M_{2}}$	0.32	0.26	0.30	0.25	0.30	0.33
$X_\mathrm{c, 1}$	0.528(9)	0.205(19)	0.532(26)	0.211(21)	0.280(2)	0.371(4)
$X_\mathrm{c, 2}$	0.558(9)	0.484(7)	0.632(11)	0.565(4)	0.495(5)	0.534(2)
$f_\mathrm{\rm CBM}$	0.035	0.005	0.005	0.005	0.035	0.035
$\log D_\mathrm{mix,0}\, (\mathrm{cm}^{2}\,\mathrm{s}^{-1})$	5	2	1	4	5	5
age (Myr)	7.69(15)	9.64(20)	4.24(51)	11.14(26)	11.29(21)	7.76(7)
Coeval_VM
Star	CW Cep	HD 254346	HD 74455	KX Vel	V453 Cygni	V1295 Sco
$M_{1}\, (\mathrm{M}_{\odot})$	12.957(46)	15.23(4)	16.29(3)	13.68(29)	14.79(16)	17.26(2)
$M_{2}\, (\mathrm{M}_{\odot})$	11.891(46)	10.68(5)*	10.30(4)	8.25(4)	11.28(11)	12.66(3)
$M_\mathrm{cc, 1}\, (\mathrm{M}_{\odot})$	4.27(2)	4.78(5)	5.96(9)	4.06(5)	4.55(8)	6.16(2)
$M_\mathrm{cc, 2}\, (\mathrm{M}_{\odot})$	3.69(2)	2.69(2)	3.09(4)	2.03(1)	3.44(5)	4.05(1)
$\frac{M_\mathrm{cc, 1}}{M_{1}}$	0.33	0.31	0.37	0.29	0.31	0.36
$\frac{M_\mathrm{cc, 2}}{M_{2}}$	0.31	0.25	0.30	0.25	0.30	0.32
$X_\mathrm{c, 1}$	0.529(4)	0.284(1)	0.551(25)	0.297(16)	0.278(2)	0.371(5)
$X_\mathrm{c, 2}$	0.546(4)	0.452(7)	0.630(11)	0.543(5)	0.497(5)	0.522(2)
$f_{\rm CBM, 1}$	0.035	0.035	0.02	0.035	0.035	0.035
$f_{\rm CBM, 2}$	0.025	0.005	0.005	0.005	0.035	0.025
$\log D_\mathrm{mix, 1}\, (\mathrm{cm}^{2}\,\mathrm{s}^{-1})$	2	4	2	2	3	2
$\log D_\mathrm{mix, 2}\, (\mathrm{cm}^{2}\,\mathrm{s}^{-1})$	2	3	2	1	5	5
age (Myr)	7.45(14)	10.62(21)	4.30(56)	12.27(28)	11.18(19)	7.69(7)
Coeval_MD
Star	CW Cep	HD 254346	HD 74455	KX Vel	V453 Cygni	V1295 Sco
$M_{1}\, (\mathrm{M}_{\odot})$	12.945(49)	15.23(4)	16.29(3)	13.68(3)	14.78(16)	17.26(2)
$M_{2}\, (\mathrm{M}_{\odot})$	11.898(46)	10.67(5)*	10.30(4)	8.25(4)	11.25(11)	12.67(3)
$M_\mathrm{cc, 1}\, (\mathrm{M}_{\odot})$	4.29(3)	4.10(6)	5.77(11)	3.29(7)	4.59(8)	6.18(2)
$M_\mathrm{cc, 2}\, (\mathrm{M}_{\odot})$	3.83(2)	2.75(2)	3.08(4)	2.07(1)	3.42(4)	4.17(1)
$\frac{M_\mathrm{cc, 1}}{M_{1}}$	0.33	0.27	0.35	0.24	0.31	0.36
$\frac{M_\mathrm{cc, 2}}{M_{2}}$	0.32	0.26	0.30	0.25	0.30	0.33
$X_\mathrm{c, 1}$	0.529(4)	0.240(17)	0.539(25)	0.232(19)	0.280(2)	0.371(5)
$X_\mathrm{c, 2}$	0.558(3)	0.474(7)	0.631(11)	0.559(5)	0.495(5)	0.534(2)
$f_{\rm CBM, 1}$	0.035	0.015	0.01	0.01	0.035	0.035
$f_{\rm CBM, 2}$	0.035	0.005	0.005	0.005	0.035	0.035
$\log D_\mathrm{mix,0}\, (\mathrm{cm}^{2}\,\mathrm{s}^{-1})$	5	3	1	2	5	5
age (Myr)	7.69(14)	9.97(21)	4.26(52)	11.31(26)	11.30(19)	7.77(7)

Separating methodological and physical contributions

We now compare the results from the various Bayesian fitting procedures (14 different models) in an attempt to distinguish results arising from the statistical formulation (7 different likelihood functions, for each of the single and binary treatment of the components) from those associated with different values of mixing in the adopted MESA model grid. Because the analyses of the likelihood functions differ primarily in the treatment of the dynamical mass and the demand of equal age (or not), their comparison may provide information on the methodological contribution to the inferred mass discrepancy for the input physics and the free parameter ranges adopted in the MESA grid.

In order to compare the 14 models we define

where $\Delta M_{\rm ehmix}$ is the discrepancy between the dynamical mass and the mass of the selected best model for the case of enhanced mixing. For $\Delta M_{\rm min}$ we consider two options: (i) we take the discrepancy between the dynamical mass and the evolutionary mass of the reference models with the absolute minimal amount of core overshooting and envelope mixing in the grid. Hence for this case, $f_{\rm CBM}=0.005$ and $\log D_{\rm mix,0}=1.0$ are fixed values and $\Delta M_{\rm min}$ only changes due to the variations in the likelihood functions; (ii) rather than using the RM as a common reference, $\Delta M_{\min}$ is the mass discrepancy between the dynamical mass and the evolutionary mass obtained from the other model configurations, evaluated at the minimum mixing parameters, $f_{\rm CBM}=0.005$ and $\log D_{\rm mix,0}=1.0$, while retaining the same likelihood formulation as adopted for their corresponding enhanced-mixing solutions. In the case of (i), $\eta$ indicates the need of extra mixing compared to models without it to resolve the mass discrepancy, where zero means no need for extra mixing. For option (ii), $\eta$ quantifies the extent to which enhanced mixing reduces the mass discrepancy across the different model configurations. Both options use all the variants of the likelihood functions considered in this work.

The outcomes of $\eta$ in these two cases are shown in Figure 7. For case (i) in the left panel we find overwhelming evidence that extra mixing is required, irrespective of the chosen likelihood formulation. The median efficiency of the mass discrepancy reduction $\eta$ is 97% for both the cases with only core overshooting or envelope mixing enhanced (’Evol_CBMoraD’), or both types of mixing together (’Evol_CBMD’). This is in agreement with conclusions in the literature presented in the introductory section of our paper. For case (ii) the KDEs shift to median $\eta$ values of 42% for ’model_CBMoraD’ and 37% for ’model_CBMD’ models. For case (i), $\eta$ is affected by both the enhanced mixing parameters and the variety in the likelihood functions, while for case (ii) the difference mainly comes from the mixing parameters alone, since the variety in the likelihood treatments should affect the numerator and denominator of $\eta$ in a similar way. Comparing both scenarios hints at approximately half of the apparent reduction in the inferred mass discrepancy is attributable to the need of extra mixing and half to the methods adopted in the construction of the likelihood functions, when considering the chosen input physics and parameter ranges of our grid of MESA modes and while evaluating our small sample of eclipsing binaries.

Discussion and conclusions

One of the principal outcomes of this study is that the inferred mass discrepancy depends not only on the adopted stellar physics and the range of free parameters associated with it but also on the construction of the statistical inference problem. By comparing otherwise identical analyses with and without dynamical masses included in the likelihood function, we explicitly isolate the methodological contribution to the inferred discrepancy. This demonstrates that likelihood construction constitutes a source of systematic uncertainty that should be considered when fitting stellar evolution models. This conclusion is complementary to the one by [19], who stressed the importance of the spectroscopic parameter estimation from atmosphere models, where variations in the microturbulent velocity and/or turbulent pressure are often ignored in the likelihood functions.

The mass discrepancy has traditionally been associated with inadequacies in the adopted physics of stellar interiors, with enhanced mixing shown to reduce its magnitude in several studies [e.g. 45, 19, 30, 67]. Consequently, discussion of the mass discrepancy has largely focused on identifying the missing stellar physics. Our results suggest that this interpretation may be incomplete. By comparing alternative Bayesian inferential formulations applied to the same stellar evolutionary grid, we find that part of the inferred reduction in the mass discrepancy depends on the statistical formulation of the inference problem itself. A precisely measured mass in the likelihood comparable to other observables with larger relative uncertainties may reduce the inferred discrepancy independently of changes in the underlying stellar physics. Our results reinforce the need to separate methodological from physical contributions to resolve the mass discrepancy.

Despite our conclusions, mixing remains one of the principal uncertainties in the evolution of intermediate- and high-mass stars. We stress that we have only used one type of mixing profile to describe core overshooting and to mimic envelope mixing, respectively, among many available choices for these two mixing phenomena [cf. 34 for a broader discussion when performing asteroseismic modelling]. CBM inferred from eclipsing binaries and gravity-mode asteroseismology exhibits substantial star-to-star variation, with reported values depending on the adopted stellar models, evolutionary diagnostics, and fitting methodology [e.g. 42, 48, 49, 6, 33, 32, 12, 64]. Whether this variation reflects genuine physical diversity or limitations of current stellar models remains an open question.

Our analysis supports the growing evidence that enhanced mixing reduces the evolutionary mass discrepancy. This is expected as stronger CBM increases the convective core mass and thus prolongs the main-sequence lifetime and increases the luminosity of a star of a given mass at a given metallicity [54, 53, 12]. Consequently, models with stronger mixing can reproduce the observed positions of the star in the HR diagram using lower evolutionary masses, thereby bringing them into closer agreement with the independently measured dynamical masses. In addition, increasing mixing parameters also alters the effective temperature, surface gravity, and main-sequence lifetime of stellar models [55, 33]. This is also corroborated by the findings in this work as shown in the HR diagram in Fig. Figure 8, where the evolutionary tracks of both the statistically preferred mixing and the minimal mixing under the coeval_VM formulation are plotted for the sample. However, mixing parameters cannot be infinitely increased as excessively large mixing parameters shift the evolutionary tracks away from the observed positions in the HR diagram or require implausible common ages to reproduce both binary components simultaneously.

In addition, the observed mass discrepancy anti-correlates with the convective core mass corroborating the long established fact in the literature that underestimation of the core mass results in more mass discrepancy. The coeval binary analyses reinforce these conclusions. Several systems are also satisfactorily reproduced using common mixing parameter values for both components, whereas others favour different mixing efficiencies for the two components. Similar star-to-star variations have been reported in binary and asteroseismic studies [e.g. 42, 19], indicating that mixing cannot be well modelled by assuming that the parameters describing them do not change along the evolution. Results by [6] also show that the shrinking convective cores of intermediate-mass stars can nowadays be measured adequately along the entire main sequence from asteroseismic observables. Keeping this in mind, it is unlikely that the internal mixing can be described by a single universal time-independent parameter [cf. 1 for the same conclusion from numerical simulations of envelope wave mixing as a function of age]. Our empirical mass-dependent CBM follows the same qualitative trend of higher CMB for higher masses, but should be regarded as exploratory because of the small sample size and limited evolutionary stages of our considered binaries.

Our results further show that models with the highest Bayesian evidence do not always correspond to those with the smallest $\Delta M$, especially if the mass is not enforced. This demonstrates that resolving the mass discrepancy alone is insufficient as successful models are expected to simultaneously reproduce the observed temperatures, luminosities, surface gravities, and coevality constraints. The Pareto-optimal solutions introduced here represent a balance between statistical fit quality and reduction in mass-discrepancy. This corroborates the growing evidence that the best stellar models are not necessarily those with the strongest amount of mixing, but rather those that achieve a physically self-consistent compromise among multiple observables [42, 34].

Bayesian inference has become an increasingly important tool for interpreting modern stellar observations. This is because it naturally combines heterogeneous observables with physically motivated prior information while providing rigorous uncertainty estimates. Most applications, however, employ a single likelihood formulation and implicitly assume that the inferred stellar parameters are insensitive to the choice of observational constraints entering the likelihood. The present analysis shows that this assumption may not always be appropriate. The alternative inferential formulations explored here to study the mass discrepancy produce systematically different posterior distributions of the discrepancy despite employing just one grid of stellar evolutionary models. Future Bayesian analyses involving a variety of stellar evolution model grids with different input physics will surely benefit from explicitly examining the sensitivity of assumed and inferred stellar parameters to the adopted likelihood formulations.

This study demonstrates that the inferred evolutionary mass discrepancy depends not only on stellar evolution physics but also on the statistical formulation of the inference problem. Comparison of the alternative likelihood formulations indicates that about half of the apparent reduction in mass discrepancy may originate from the statistical treatment of the observational constraints rather than solely from changes in the mixing parameters adopted in the used grid of stellar models. The resulting central hydrogen mass fractions for the stars in the coeval cases, ranging from 0.632(11) to 0.205(19), show that the sample spans the early to late main-sequence evolutionary stages. Their coeval ages range from 4.25(51) to 12.27(28) Myr and are determined with precisions of $2$ – $13\%$ for the adopted input physics. The methodology presented here could be extended to larger samples and to future analyses incorporating rotation, binary interaction, asteroseismic constraints, or more sophisticated stellar evolution calculations with time-dependent mixing parameters.

Acknowledgements

This work was supported by the Polish National Science Foundation (NCN) under grant nr. 2021/43/B/ST9/02972 led by GH. DMB gratefully acknowledges UK Research and Innovation (UKRI) in the form of a Frontier Research grant under the UK government’s ERC Horizon Europe funding guarantee (SYMPHONY; PI Bowman; grant number: EP/Y031059/1), and a Royal Society University Research Fellowship (PI Bowman; grant number: URF\R1\231631). MV and CA acknowledge financial support from the European Research Council (ERC) under the Horizon Europe programme (Synergy Grant agreement N$^{\circ}$101071505: 4D-STAR). While partially funded by the European Union, views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. CA further acknowledges financial support from the Flemish Government under the long-term structural Methusalem funding program by means of the project SOUL: Stellar evOlution in fUll gLory, grant METH/24/012 at KU Leuven, Belgium.