© The Authors 2025

1. Introduction

A significant fraction of stellar systems in our Galaxy are formed by binary stars (e.g. Moe & Di Stefano 2017; Niu et al. 2021; Torres et al. 2022). The physical processes involved in the evolution of such systems are more varied and, at the same time, often more complex than those describing the evolution of their single counterparts. Angular momentum transfer, magnetic braking, gravitational wave emission, rejuvenation and ageing, orbital synchronization, and, especially, mass transfer episodes, are some examples of such physical processes characteristic of binary evolution.

Unfortunately, many aspects of these processes remain poorly understood. Precise observational data coupled with comprehensive statistical analyses are essential to constrain certain theoretical models. Of particular interest are those binary systems formed by at least a white dwarf. These objects represent the most common remnant among stellar evolution and their physical properties are reasonably well understood from a theoretical point of view (see Althaus et al. 2010; Isern et al. 2022, for recent reviews). This understanding enables the derivation of precise evolutionary cooling models, thereby establishing white dwarfs as reliable cosmochronometers (Fontaine et al. 2001). The photometric magnitudes of white dwarfs are thus correlated with their cooling age. Hence, with an appropriate initial-to-final mass relationship we can derive the mass of the progenitor stars, their lifetime in the main sequence, and ultimately, the total age of the system. This information, when coupled with other properties of the non-evolved companion, such as their metallicity, kinematics, or activity, can lead to important age-relation estimates (e.g Rebassa-Mansergas et al. 2021, 2023; Raddi et al. 2022).

However, during a mass transfer episode, the inference of the initial progenitor’s parameters can be especially challenging. This is particularly true when dynamically unstable mass transfer occurs, thus leading to mass overflow through the Roche lobe. This generally results in a common envelope phase through which the binary separation decreases dramatically. When formed of a white dwarf and a main-sequence star, the resulting system typically exhibits a short orbital period ranging from a few hours to a few days (Rebassa-Mansergas et al. 2008; Nebot Gómez-Morán et al. 2011). These systems are of special interest as precursors of cataclysmic variables (e.g. Pala et al. 2022), super-soft X-ray sources (e.g. Parsons et al. 2015) and double degenerates (e.g. Napiwotzki et al. 2020). All these exotic objects may later result in a type Ia supernovae (e.g. Wang & Han 2012).

Several attempts to reconstruct the progenitors of these kinds of systems, that is, the initial conditions of the parent stars, can be found in the literature (e.g. Nelemans & Tout 2005; Zorotovic et al. 2010; Davis et al. 2011; Hernandez et al. 2021, and references therein). However, given the irreversible nature of many stellar processes (particularly the common-envelope phase), reversing this process to precisely reconstruct the complete past history solely from final observable parameters is a tough task. An approach to avoid these issues is to adopt fixed values for certain parameters or relying on a grid of predefined models. However, this limits the exploration of the initial parameter space, which, in our opinion, is one of the major drawbacks of the methodologies used so far.

Detailed binary evolution codes, such as MESA (Paxton et al. 2015) or BINSTAR (Siess et al. 2013), allow for precise tracking of the physical processes governing the evolution of binary systems. However, these tools are generally time-consuming and can only simulate stellar evolution forward in time based on initial conditions. On the other hand, binary population synthesis allows for faster, although less precise, simulations of a sample of binary stars from a given set of conditions. Hence, techniques to effectively identify the progenitors of a specific set of binary stars using population synthesis are required. Recent advances, such as the dart_board code (Andrews et al. 2018), integrate Markov Chain Monte Carlo methods with rapid binary evolution models, improving efficiency and allowing for flexible modeling of both populations and individual binaries, while naturally incorporating observational uncertainties and additional constraints.

In this paper, we aim to develop an algorithm based on a population synthesis code capable of identifying the initial parameters that, through binary evolution, result in a set of final parameters consistent with the observable ones. No restrictions are imposed on the initial parameters, allowing us to identify possible correlations among these parameters. Time efficiency and the ability to be applied to a wide range of binary systems are also two desirable characteristics of our algorithm. Finally, the algorithm will be applied to a set 30 white dwarf plus main-sequence (WDMS) eclipsing binaries identified by the Zwicky Transient Facility survey (Brown et al. 2023). The accurate estimation of stellar parameters in these systems represents a clear advantage over previously used samples subjected to substantially larger observational uncertainties and will result in a restricted sample of compatible progenitors. Consequently, this restriction will refine the conditions of theoretical models and shed light on the evolution of these systems.

2. Population synthesis code

In this work, the MRBIN population synthesis code is used for the simulation of binary systems. It was initially based on the binary stellar evolution code (BSE) developed by Hurley et al. (2002) but incorporated several updates provided by Camacho et al. (2014), Cojocaru et al. (2017) and Canals et al. (2018). The code is a multipurpose binary population synthesis simulator, but in particular it has been widely used for the study of the white dwarf population in binary systems (see, for instance, Torres et al. 2022, and references therein).

Our code belongs to a family of synthetic stellar evolution codes (other examples include, binary_c-python Hendriks & Izzard (2023); StarTrack Belczynski et al. (2008); SeBa Toonen et al. (2012), and references therein). These algorithms stand out for being relatively fast, as they rely on a set of adjustments of analytical formulas that replace precise and highly time-consuming calculations such as detailed stellar-structure and binary-evolution codes, for example MESA (Paxton et al. 2015) or BINSTAR (Siess et al. 2013). Despite this, the error introduced in the approximation does not exceed 5% and they also enable a broader parameter analysis space compared to the more detailed evolution code counterparts (see Fragos et al. 2023, for a recent and detailed analysis).

Consequently, the use of MRBIN is justified as an initial proxy for the progenitor’s solution. A detailed analysis of each system using a binary evolution code is beyond the scope of the present study and deferred to future work.

Several input parameters are needed to be considered for the evolution of a binary system. Among them, the most relevant are the mass of the primary star (the initially more massive main-sequence star), the mass of the secondary star, the eccentricity, the separation, and orbital period, the efficiency of the common-envelope phase, the metallicity, and the age when the system was formed (in our analysis we assumed that the formation of both stars are coeval). Additionally, the code incorporates a set of prescriptions that account for other aspects of binary evolution such as tidal evolution, stellar winds, gravitational radiation, magnetic braking, and angular momentum losses. All these processes introduce an additional set of parameters as described in Hurley et al. (2002), which are kept constant throughout our simulations.

For the treatment of the common-envelope phase we used the α prescription (Webbink 1984; Tout et al. 1997). This formalism, essentially based on the energy conservation principle, introduces the α_CE parameter as a measure of the efficiency of the orbital energy to unbind the common envelope of the system:

$$\begin{array}{c}\hfill E_{\mathrm{bin}}={\alpha }_{\mathrm{CE}}\mathrm{\Delta }{E}_{\mathrm{orb}},\end{array}$$(1)

where E_bin is the binding energy of the envelope, and ΔE_orb is the change in orbital energy. An approximation in a compact form for the binding energy is commonly employed:

$$\begin{array}{c}\hfill E_{\mathrm{bin}}=-\frac{G{M}_{\mathrm{donor}}{M}_{\mathrm{env}}}{\lambda R_1},\end{array}$$(2)

where M_donor, M_env, and R₁ are the total mass, the envelope mass, and the radius, respectively, of the primary star at the beginning of the common-envelope phase, and λ is the binding energy parameter (de Kool et al. 1987; Hurley et al. 2002). This parameter accounts for the approximation of the previous equation with respect to the exact value of the gravitational binding energy. The λ parameter is approximated by a series of values that mainly depend on the type of star and its envelope mass (see Claeys et al. 2014, and references therein). In addition to the binding energy, we can take into account the internal energy, mainly due to the hydrogen recombination energy. However, it has been estimated that this internal energy is negligible in the case of post-common envelope WDMS binaries (Zorotovic et al. 2010; Camacho et al. 2014), as the systems in our study sample. Consequently, we have adopted a null internal energy contribution in our model.

Finally, our code incorporates some updates with respect to BSE version. For instance, we used the PARSEC evolutionary tracks (Bressan et al. 2012), providing an updated version of the mass-radius and, thus, effective temperature-luminosity relation for the main-sequence stars. Analogously, the white dwarf evolution is followed according to La Plata cooling tracks (Althaus et al. 2015, 2021; Camisassa et al. 2016, 2017, 2019) and the corresponding photometry by using Koester’s (Koester 2010) hydrogen-rich atmospheric models. The set of La Plata models incorporate the most updated physical processes in white dwarf evolution, such as crystallization, phase separation and sedimentation of major species such as Ne22, as well as covering the full range of chemical core compositions (helium, carbon-oxygen, and oxygen-neon).

3. Methodology: inverse population synthesis

Reconstructing the past history of post-common envelope WDMS binaries involves determining the parameters of their progenitors: their masses, separations, eccentricities, ages, and ultimately, a set of initial parameters. In the following sections, we outline the strategy presented in our work, which is based on the population synthesis code, MRBIN, described in the previous section.

3.1. The inverse population synthesis algorithm

In Figure 1 we show a schematic representation of the problem. Σ is the space containing all possible initial values, while Ω represents the space containing all possible output values after binary evolution. Our goal is to find the subset S ⊂ Σ, such that a system, x_in = (x₁,x₂,…,x_n)_in, that belongs to S, after binary evolution leads to the final set of parameters, y_out = (y₁,y₂,…,y_m)_out, that is at a distance, d, smaller than a certain σ from the set of observed parameters, ${\mathit{y}}_{\mathrm{obs}}={(y_1,y_2,\dots ,y_m)}_{\mathrm{obs}}$. As an initial assumption we can adopt a Euclidean normalized distance between the output point and the observed point:

Fig. 1. Schematic representation of the problem. The space that contains all possible initial values (Σ) is projected, after binary evolution, onto the output space (Ω). Only those points within S lead to a final set of parameters, called O, that are at a distance, d < σ, from a certain set of observed parameters. See text for details.

$$\begin{array}{c}\hfill d=\sqrt{{\left(\frac{y_{1,\mathrm{out}}-y_{1,\mathrm{obs}}}{y_{1,\mathrm{obs}}}\right)}^2+{\left(\frac{y_{2,\mathrm{out}}-y_{2,\mathrm{obs}}}{y_{2,\mathrm{obs}}}\right)}^2+\dots +{\left(\frac{y_{m,\mathrm{out}}-y_{m,\mathrm{obs}}}{y_{m,\mathrm{obs}}}\right)}^2}.\end{array}$$(3)

The goal of the algorithm is to find S in an efficient and computationally low-cost way. In order to do that the algorithm performs a random walk in the initial space (Σ) to identify all possible combinations of parameters that lead to a desired result in the observable space (Ω). That is, by applying random displacements, the algorithm explores a set of initial parameters which are evolved using the MRBIN population synthesis code into a set of output parameters. Next, we calculate the distance (Eq. 3) of this output point to a given point in the observable space. This distance is compared with an arbitrary σ, which represents the allowed error threshold (e.g., 1%, 5%, 10%, etc.). If the result is within the allowed error range (i.e., d ≤ σ), the set of initial parameters is considered valid and saved. Otherwise, the point is discarded and the algorithm returns to the previous point. Then, a new set of random displacements is applied and the process is repeated.

To maximize the efficiency of the algorithm (that is, to find the border that delimits the solution space, S, in the minimum number of steps), adaptive steps between points can be used. This means that if the initial generated point is considered valid, the next step taken from that point will be larger than the previous one. However, if the new point exceeds the maximum error allowed, the algorithm will return to the previous point and generate another one using a smaller step size. In our case, we adopt a scaling factor, ϵ, of 1.25; that is, we increase or decrease each new random displacement by 25%.

In Figure 2 we plot an example of the random walk process in a 2-dimensional space. It starts from point 1 (green circle) inside the solution space S. It makes a first step of length Δ_1 → 2 to point 2, that is, a random displacement in each of the n-dimension of the input space Σ. The new point is evolved and checked to be valid (i.e., d ≤ σ), so that the algorithm continues from it. The next step is now increased by the scaling factor, ϵ = 1.25, so that Δ_2 → 3 is on average a 25% larger than before. The algorithm falls to point 3, which is also valid. The next step is even larger, but in this case point 4 is invalid (i.e., d > σ), so the algorithm has to return to point 3 and do a new step from there, now a 25% smaller. After two valid points (5 and 6), the same happens with point 7, implying that the algorithm goes back to point 6 and continues in another direction with a smaller step. Points 8, 9, and 10 are consecutively valid, so that larger and larger steps are taken.

Fig. 2. Schematic representation of our adaptive random walk algorithm. See text for details.

To minimize the time required for the algorithm to explore the solution space, we initiate a new random walk from any arbitrary solution point that has been found thus far. For instance, point 9 is adopted as the new starting point for the random walk, with points 2′, 3′, and 4′ marked in blue. Since this new random walk operates independently from the previous one, both can be executed in parallel. Moreover, since this process can be applied iteratively, it leads to a significant reduction in the overall computation time.

Finally, if the initial guess does not belong to the solution space, S, a similar adaptive random walk process is applied until we locate a solution within S. In our example of Fig. 2 that would be the case when starting with point 1″ (red circle). If the new point, 2″, is at a distance, d, larger than the previous one, we return to the original point since we find a new one, point 3″, whose distance is smaller. This process is repeated until we finally find a point, 4″, that belongs to the solution space, S. It then continues as a random walk as previously described.

3.2. Validation of the algorithm

We have tested the feasibility of our algorithm with the assistance of our population synthesis code, MRBIN. Since, for synthetic systems, we know both the initial and final values, we can thus verify the capability of our algorithm to find solutions given a set of final values that would play the role of observed values. Different combinations of input and output parameters have been tested, as well as different error thresholds (σ).

In all cases, we have been able to verify that the algorithm based on adaptive random walks can find a stable set of solutions. That is, the distribution of the initial parameters is sufficiently smooth, and they are independent of the algorithm’s starting point. This is accomplished as long as the number of generated points is sufficiently high (on the order of 10⁶), and a minimum of 100 random walks are generated.

For illustrative purposes we have chosen here three representative cases that correspond to WDMS systems with final periods of the order of a few hours, ten hours, and days (P = 1.16,  15.50, and 41.73 h, respectively). The output space consists of five variables (these five variables are the ones that we can get observationally, see Section 4): the masses and effective temperatures of the WDMS, and the orbital period of the system. The free parameters correspond to six variables of the input space: the initial masses of the stars, the initial period and eccentricity, the common envelope efficiency parameter, and the age of the system. The metallicity has been fixed to the solar value, and for the three systems shown here we adopted α_CE = 0.3.

Results are displayed in Figure A.1 using a corner plot, which represents the marginal distribution of each of our considered input parameters (top histograms) and density maps of possible projections between any two parameters (rest of sub-panels). Furthermore, correlation matrices are also generated in order to identify any potential relationships between parameters. Each term of the matrix, Pearson correlation coefficient, is calculated as:

$$\begin{array}{c}\hfill r={\displaystyle \frac{\sum ((X-\overline{X})·(Y-\overline{Y}))}{\sqrt{\sum {(X-\overline{X})}^2}·\sqrt{\sum {(Y-\overline{Y})}^2}}},\end{array}$$(4)

where X and Y represent the respective values of two parameters, while $\overline{X}$ and $\overline{Y}$ denote the means of the X and Y variables, respectively. We recall here that correlation coefficients range between −1 and 1 in the off-diagonal elements of a correlation matrix. A value of −1 indicates a perfect negative linear dependence, +1 indicates a perfect positive linear dependence, and 0 indicates no linear dependence at all.

The analysis of the results presented in Fig. A.1 reveals a similar behavior of the space parameter regardless of the final period of the systems. In all cases, a smooth distribution of the input parameters is found, indicative that the algorithm has reached a stable set of solutions. Moreover, a positive linear dependence, r ∼ 0.9, is presented between the common envelope efficiency, α_CE, and the mass of the white dwarf progenitor, M₁. Similarly, a negative dependence is found between the age of the system and M₁ and α_CE. These dependencies reflect the fact that a specific final period can be achieved by a more massive primary progenitor in a shorter time when using a larger α_CE parameter. Conversely, when the primary progenitor is less massive, a smaller α_CE is required, but a longer age is needed.

The previously mentioned high correlation reflects a strong degeneracy in the solution space, S. Different combinations of the initial parameters (see left panels of Fig. A.1) reproduce a set of similar final parameters within a certain σ threshold (1% in the case of the systems represented in Fig. A.1). Among the set of possible solutions, x_s ∈ S, we have adopted the most representative or optimal one, x_opt, as the most probable one. We estimated the optimal solution as that which maximizes the joint probability distribution,

$$\begin{array}{c}\hfill p({\mathit{x}}_{\mathrm{opt}})=max{\displaystyle \prod _i}{f}_{X_i}(x_i)|\mathit{x}\in S,\end{array}$$(5)

where f_Xi(x_i) is the marginal distribution of the input parameter X_i. The fact that the probability is constrained to those points, x, belonging to the solution space, S, implicitly incorporates the possible correlations between parameters. Moreover, we can estimate the error in the optimal solution from those points, x_err, whose likelihood-based confidence region is within 1σ from the probability of the optimal solution, x_opt, that is, $p({\mathit{x}}_{\mathrm{err}})\ge p({\mathit{x}}_{\mathrm{opt}})/\sqrt{e}$.

In Table 1, we present the synthetic input parameters for the three systems analyzed here, along with the corresponding values for the optimal solution (the one that maximizes the joint probability, Eq. (5)) for the three threshold distances: 0.1%, 1%, and 10%, respectively. The analysis of Table 1 reveals that the input parameters found by our algorithm are either perfectly close or quite close to the synthetic ones for threshold distances of 0.1% and 1%, respectively. However, this is not the case when the threshold distance is increased to 10%. In that situation, some of the input parameters, such as the period or eccentricity, clearly diverge from the synthetic ones. Even though, it is important to emphasize here that the solution found by our algorithm in that case is still completely valid. When enlarging the threshold distance, the space of possible solutions, S, is also enlarged. In that situation, the synthetic input values used as a reference for our validation systems do not necessarily represent the ones that maximize the joint probability distribution. Given the high correlation between some of the input parameters, solutions with a higher probability than the one derived from the synthetic values can be found.

Nevertheless, some caveats and warnings regarding the method should be addressed. First, provided that a sufficiently large number of random walkers are generated, we can reasonably assume that the method recovers the entire solution space consistent with the observational constraints. Second, observational uncertainties must be small enough so that the corresponding solution space, S, is also narrow. In this limited region, we assume that prior distributions (such as the initial mass function) can be approximated as approximately constant. Thus, the impact of including full priors would introduce only second-order corrections to the recovered parameter distributions.

As long as these two conditions are met, we can approximate the optimal solution by maximizing the product of the marginal distributions, thereby minimizing the effects of possible correlations between parameters. While this approach has both advantages and limitations, one of our broader goals is to ultimately derive an estimate of the priors from the collective properties of the final set of solutions (see Section 5.2).

In this sense, our method is not strictly a Markov Chain Monte Carlo method, as described in Andrews et al. (2018), since it does not apply a probability criterion based on the posterior distribution to guide the random walks. Instead, our approach is closer to a density estimation technique, aiming to approximate the probability density function rather than merely sampling from it. In conclusion, the algorithm presented here for inverse population synthesis reconstruction can be considered as a feasible tool for exploring the entire space of initial variables and efficiently finding the set of solutions and their possible correlations, within a reasonable low computational cost.

4. The observed sample

In this work we analyze the sample presented in Brown et al. (2023). In their work, they conducted a photometric follow-up of 43 eclipsing WDMS binaries detected by the Zwicky Transient Facility (ZTF) survey¹, performing an extensive photometric campaign, and acquiring precise light curves for each object. Of the 43 systems that they followed-up, 9 do not have measured parameters because they harbor magnetic white dwarfs. In addition, the light-curve fits of 4 objects suggest the white dwarf companions to be brown dwarfs, so that the best fit secondary masses can only be regarded as upper limits. Taking this into account we will apply the algorithm to the remaining 30 systems. After fitting the best model to each observed system, Brown et al. (2023) found a set of stellar parameters for each component of the system. In Table B.1 we show the most relevant stellar parameters obtained in that work and the corresponding observational errors. Among the parameters available in their study, we adopted the more representative ones, including the observed masses and effective temperatures of both components, as well as the orbital period of the system. Typical observational errors for these parameters are of the order of 5%.

To accurately reconstruct the past histories of these systems, it is also necessary to incorporate the metallicities of the stars, a parameter that is not provided by Brown et al. (2023). We utilized the PARSEC isochrones (Bressan et al. 2012) for this purpose. Using the observed values of effective temperature and mass, we calculate the metallicities by performing multivariable interpolation on the corresponding isochrones. This procedure has been applied to our observed sample, obtaining a metallicity proxy value for all systems except for ZTF J140537.34+103919.0. Due to its lowest mass companion (reported as just above the hydrogen-burning limit, Brown et al. 2023), proper determination of its metallicity is not feasible. Therefore, a solar value of Z = 0.014 has been adopted. The list of metallicities derived is presented in the right column of Table B.1. An average value of Z = 0.0145 with a dispersion of σ_Z = 0.003 has been found, thus indicating the thin disk origin of the vast majority of this sample and the perfect agreement with metallicity estimates of WDMS binaries (Rebassa-Mansergas et al. 2021).

5. Results

5.1. Reconstructing individual binary systems

Our algorithm has been applied to each of the 30 systems from Brown et al. (2023) presented in Table B.1. We remind the reader that the observed space parameter consists of the masses and effective temperatures of both stars, and the orbital period of the system. For the input space parameter (that is, our free parameters) we adopted a six-dimensional space, consisting in the masses of the progenitor stars, the initial eccentricity and orbital period, the common envelope efficiency, and the age when the system was born. For each system we adopted the metallicity presented in Table B.1 as a constant value along the simulation.

In Figure C.1 we present the distribution of input parameters found for three representative cases and the associated corner plots (left panels), together with the corresponding correlation matrices (right panels): ZTF J041016.82-083419.5, ZTF J063808.71+091027.4, and ZTF J122009.98+082155.0. We remark that the distribution found for each parameter (histogram panels of each corner plot) is reasonable smooth, thus indicating that a stable solution has been found.

A closer look reveals that, while the progenitor masses have a symmetric, near Gaussian distribution, that is not the case for the rest of parameters. In these cases, the asymmetric distribution implies that the average value (mean) does not agree with the most probable value (mode). Inspection of the correlation matrices also revealed a strong correlation among some input parameters. That is the case of the common-envelope efficiency parameter, α_CE, and the mass of the primary, M₁, where a positive correlation exits (r ∼ 0.9 ∼ 0.95) for most of the systems. Similarly, a negative or anticorrelation arises between the age and M₁ and between the age and α_CE (r ∼ −0.9 ∼ −0.95). These facts imply that the initial parameters are not independent variables. Therefore, if we were to adopt a solution set comprising the most probable or average value for each individual parameter, it would likely yield an incompatible solution with the observed values. Consequently, we have adopted the set of input values that correspond to the optimal solution, that is the solution that maximizes the joint probability distribution of the initial parameters (see Section 3.2), for the remainder of the discussion.

In Table C.1 we show the values corresponding to the optimal solution of our analyzed systems. Also, for sake of completeness, we include in the Appendix C Tables C.2 and C.3, corresponding to the mean and mode values, respectively. However, we should remark, as previously stated, that the initial conditions derived from the mean and mode values for each individual initial parameter, given the strong correlation among some of them, may not necessarily represent a solution of the binary system.

For the set of optimal solutions we computed the distance to the observed values following Equation (3). In Table C.1 these distances are shown as a percentage under column d_min. Likewise, we display in the last column of the table, the observational error, σ_obs; calculated as the sum in quadrature from the error estimates presented in Table B.1 for M₁, M₂, T_eff, 1, T_eff, 2, and, as a result of error propagation, 2.5% for the period. This σ_obs is adopted as our threshold error in the computation of solutions. An initial analysis revealed that for only 11 systems, the set of solutions provide a distance below the observed error threshold.

To investigate this result in more detail, we have depicted in Figure 3 the output set of values for our optimal solution, (M₁, M₂, T_eff, 1, T_eff, 2, P)_out, with respect to the corresponding observed values, (M₁, M₂, T_eff, 1, T_eff, 2, P)_obs. As a visual aid, we marked (solid black line) the 1:1 relationship between parameters. A practically perfect agreement is achieved for all parameters except for T_eff, 2. In that case, the obtained solution is systematically hotter than the observed value. This issue was reported in Brown et al. (2023, see Fig. 5.).

Fig. 3. Parameter comparison of output values for the solution that provides the minimum distance, (M1, M2, Teff, 1, Teff, 2, P)out, with respect to the observed values, (M1, M2, Teff, 1, Teff, 2, P)obs. As it can be seen, all parameters perfectly agree except Teff, 2, which shows a clear shift between the output versus observed values. See text for details.

From this analysis we conclude than the main source of error arises from the discrepancy in the effective temperature of the secondary. We can account for this problem by adopting an ad hoc shift of 310 K (dashed black line) in the observed effective temperature estimates or, equivalently, in theoretical model. The new distances for Equation (3) we would obtain in this case are written in Table C.1 under the column d′_min. With the exception of two cases, all the others now exhibit a distance below the observed threshold error. These two (ZTF J125620.57+211725.8 and ZTF J140537.34+103919.0) and have been reported in Brown et al. (2023) as systems of special interest.

ZTF J125620.57+211725.8 was originally identified in Rebassa-Mansergas et al. (2021) as an unresolved binary system formed by a extremely low mass (below 0.2 M_⊙) white dwarf. On the other hand, the mass estimate of the white dwarf by Brown et al. (2023) was notably larger, specifically 0.48 ± 0.01 M_⊙. Our inverse population synthesis study reveals that is not possible to find a set of progenitor parameters that reproduce the observed fitted values presented in Brown et al. (2023). Thus, we reinforce the idea of conducting spectroscopic follow-up alongside a detailed analysis of binary evolution for this particular system.

Regarding ZTF J140537.34+103919.0, this system was reported in Brown et al. (2023) as a binary with a possible sub-stellar companion, whose secondary star exhibits an effective temperature hotter than expected for its mass. The same authors suggested the possibility of it being a triple system to explain their observed peculiarities. Nevertheless, a thorough analysis of this system is deserved to determine the progenitors and understand their past history.

Finally, it is important to note that, although the data quality of the sample analyzed here is sufficient to yield observational uncertainties on the order of a few percent, this level of uncertainty slightly surpasses the ideal conditions (see Section 3.2) required to accurately recover the input parameters. Consequently, some caution should be exercised in the interpretation of the results presented here and in the forthcoming sections.

5.2. Global properties of the progenitors

Once each individual system was analyzed, we studied the properties of the entire sample (28 systems that have compatible solutions with the observed errors; see Table C.1). Recall that we adopted as progenitor solutions the set of parameters that resulted in the maximum joint probability within the observed parameter space. In Figure 4 we display a corner-plot for the progenitor parameters (left panel) and the respective correlation matrices (right panel). As input parameters we have used the masses of the progenitor stars, M₁ and M₂, the eccentricity and orbital period, e and P, the common envelope efficiency parameter, α_CE, the age of the system, T_age, and we have also added the metallicity, Z, as computed in Section 4.

Fig. 4. Distribution of parameters for the entire sample.

Some correlations between the parameters are detected. First, a clear anticorrelation, r = −0.801, arises between the total age and the mass M₁. As expected from stellar evolution theory, primary low-mass progenitors (roughly around 1 M_⊙), spend more time until they overfill their Roche-lobes. This natural effect is also increased if we consider that we have estimated low-metallicities for those low-mass primaries.

Second, a mild anticorrelation, r = −0.285, is shown between the α_CE parameter and the mass of the secondary. For a better analysis, we plot in Figure 5α_CE versus M₂ with its corresponding error bars and marked as a dashed red line the linear decreasing trend. The confidence interval for the slope is [−1.656,0.359]  (M_⊙)⁻¹ at a 95% confidence level, indicating that a null or even positive slope cannot be entirely excluded. However, the relation found suggests that the largest values of the common-envelope efficiency parameter correspond to low-mass secondaries, although not all of these systems present a large α_CE. On the other hand, more massive secondaries present low-values of α_CE. Moreover, the analysis of the distribution of the common-envelope efficiency parameter reveals that there is not a unique value (on the contrary there is a wide range from α_CE = 0.087 to 0.977) capable of reproducing the entire sample. The average value of the common-envelope efficiency for the whole sample result to be ⟨α_CE⟩ = 0.588, which is slightly larger than in previous analysis of post-common envelope WDMS binaries that assigned a low efficiency parameter, α_CE ∼ 0.2 − 0.4, (e.g. Zorotovic et al. 2010; Toonen & Nelemans 2013; Camacho et al. 2014). This is partially due to the fact that some of the systems with the lowest secondary-mass require a higher efficiency, in some cases α_CE is nearly as high as 1. For instance, this is the case for ZTF J164441.18+243428.2 and ZTF J102653.47-101330.3, both with M₂ < 0.11 M_⊙ and with α_CE = 0.932 and 0.942, respectively. This result also suggests that for low-mass secondaries, the inclusion of other sources of energy, in addition to the orbital energy, such as the thermal envelope energy or internal energy due to dissociation and ionization energy of the envelope, may be considered (Rebassa-Mansergas et al. 2012).

Fig. 5. Common envelope efficiency parameter, αCE, as a function of the mass of the secondary. The decreasing linear trend is marked by the dashed red line.

The aforementioned suggestion seems to be in line with previous analysis by De Marco et al. (2011), Davis et al. (2011). These authors hypothesized that lower mass companions require more time than a stellar dynamical time to spiral into the core of the giant. This prolonged process enables the giant to use its internal thermal energy to facilitate the unbinding of its envelope. This effect should be more evident in the case of brown dwarf companions instead of M dwarfs. However, in a recent analysis of a sample of close detached white dwarf plus brown dwarf binaries Zorotovic & Schreiber (2022) indicated that these systems can be reconstructed with a low-efficiency parameter without the need to include other sources of internal energy for expelling the envelope. However, it should be noted that this result is achieved by imposing the condition that the age of the analyzed systems should be above a minimum age derived from the literature. If this condition is omitted, larger values of the efficiency parameter are recovered, indicating in that case the contribution of a certain amount of internal energy.

It is worth mentioning here that our algorithm does not impose any restrictions on the input parameters, which include the total age of the system. Furthermore, since our observational sample includes a precise estimation of the temperatures of both stars, the algorithm is thus able to properly constrain the total age of the system.

Continuing our analysis of the output parameters, we also observed a certain anticorrelation between the common envelope efficiency and the age of the system. Large values of α_CE are linked to young ages. However, the opposite is not entirely true, as low efficiency values are associated with both young and old systems. This fact implies that, at least for the sample analyzed here, there is no clear anticorrelation between the common envelope efficiency and the mass ratio, q = M₂/M₁, as proposed by De Marco et al. (2011).

Another noteworthy result of the set of solutions found by our algorithm concerns to the eccentricity distribution. In Figure 6 we plot the initial eccentricity distribution for the systems analyzed here. For visual analysis we depict the canonical thermal distribution (dashed red line), f = 2e. We should be aware that a larger phase space at higher eccentricities is available for the random walkers (for instance, binary systems circularized at periastron upon Roche lobe overflow, allowing for a much wider range of orbital separations). However, given that the optimal solution corresponding to the product of the marginal distributions is maximized, individual biases can, to some extent, be compensated. For instance, binaries with shorter initial periods are naturally more prone to Roche lobe overflow and thus tend to interact earlier (as shown in the marginal distributions of the corner plots in Figure A.1). These effects therefore compensate for the larger phase space available at higher eccentricities when determining the optimal solution for the system. We conclude that, as long as the observational errors are small, the optimal solution can be considered a reasonable estimate, and therefore the distribution obtained from the analyzed sample reflects the properties of that sample. Hence, despite the fact that our sample cannot be considered complete by far, and the existence of some initially already circularized systems (∼10%), the eccentricity distribution found in this work sample seems to be compatible with the thermal one. A Kolmogorov-Smirnov test yields a p-value of 0.2054, failing to reject the null hypothesis that the data follows the thermal distribution. This result suggests that close binary systems exhibit an initially thermalized distribution of eccentricities, rather than a uniform or even superthermal distribution found in wide binaries (e.g., Geller et al. 2019; Hwang et al. 2022).

Fig. 6. Initial eccentricity distribution for the systems analyzed in this work. Despite the scarcity of the data, the distribution seems to be consistent with the canonical thermal distribution (red line), f = 2e, although some initially already circularized systems appear in the sample.

Finally, our study of the sample of eclipsing binary systems from Brown et al. (2023) has allowed us to estimate the metallicities (see Section 4) and, through our reconstruction algorithm, the ages of such binary systems. Hence, we can analyze the resulting age-metallicity distribution for objects in our sample. In Figure 7 we show the distribution obtained (black dots) along with a best-fit linear trend (dashed red line). For visual comparison we also plot the age-metallicity values (gray dots) for the sample of wide WDMS binaries from Rebassa-Mansergas et al. (2021). The dispersion in metallicities presented in our sample for younger ages (up to ∼3 Gyr) is in agreement with that presented in Rebassa-Mansergas et al. (2021). However, a noticeable trend towards lower metallicities is shown in our sample. Although the 95% confidence interval for the slope, [−3 870,799]  (yr⁻¹), reveals that it is not entirely incompatible with a null slope. In any case, the lack of older high-metallicity objects in our sample, in contrast to the sample of Rebassa-Mansergas et al. (2021), where this dispersion in metallicities is consistently maintained across all ages, can be attributed to the shorter main-sequence lifetime of these objects compared to their low-metallicity counterparts. Hence, they can enter the common-envelope phase earlier and spend more time cooling, becoming too faint to be detected, or even lead to a merger of the binary system. A complete volume-limited sample would provide valuable insight into this issue.

Fig. 7. Age-metallicity relationship for the objects analyzed in this work (black points). For visual comparison we also show the linear trend of our sample (dashed red line) and the sample of wide WDMS binaries from (Rebassa-Mansergas et al. 2021, gray dots, RM+21).

6. Conclusions

We have developed a binary stellar reconstruction algorithm based on inverse population synthesis techniques. The algorithm, employing adaptive random walks, efficiently explores the input parameter space and identifies a set of solutions. This set consists of initial parameters that, when the binary system evolves, result in output parameters compatible with the observed ones within a certain error threshold.

The algorithm was validated with the aid of MRBIN, a population synthesis code devoted to binary evolution and based on BSE code. We verified that our algorithm can identify a stable solution after computing approximately 10⁶ points per system and around 100 random walks in parallel, representing a reduced computational time. At the same time, this solution is obtained without the need to impose any restriction on the input variables.

Once validated, we applied our algorithm to a set of eclipsing WDMS binaries detected by the ZTF survey. The high-precision parameters of this set of observed systems allows us to obtain an accurate estimate of the progenitor parameters without the need to impose any external constraints. Out of the 30 systems analyzed, we found a valid solution, meaning that the distance between the output parameters and the observed values is below the error threshold, for 28 of them.

The input parameter space for those 28 systems is defined by the masses of the progenitor stars, M₁ and M₂, the eccentricity and orbital period, e and P, the common envelope efficiency parameter, α_CE, and the age of the system, T_age. We have also made a priori estimates of the metallicity, Z, computed through interpolation of PARSEC sequences, and kept fixed throughout the application of the algorithm. The output parameters correspond to the masses and effective temperatures of the observed objects, as well as their orbital periods: (M₁, M₂, T_eff, 1, T_eff, 2, P)_obs. The output parameter values provided by our algorithm perfectly match the observed ones, except for T_eff, 2. In this case, a shift of 310 K is adopted due to a calibration issue in the synthetic evolutionary tracks.

The analysis of the global properties has revealed significant trends and correlations among certain parameters. They can be summarized in the following points:

• A high degeneracy exists in the input parameter space, although the excellent precision of the observational data to which we apply our algorithm allows us to constrain the input values to a reasonable extent.

• Not a universal value of the common envelope efficiency, α_CE, allows for the reconstruction of all systems; instead, a wide range of values is needed. In fact, the average α_CE value we find is ≃0.6, which is slightly higher than the 0.3–0.4 values generally assumed in previous works.

• The α_CE parameter exhibits a mild anticorrelation with the secondary mass. Although no internal energy is required for the systems analyzed here, larger values found for the lowest secondary masses suggest that internal energy may be considered in this regime. Furthermore, the extrapolation of this result also suggests that the inclusion of this extra source of energy can be present and may even play a major role in systems with brown-dwarf companions.

• Initial eccentricities are compatible with a thermal distribution, although a small fraction (∼10%) of systems are initially circularized.

• The age-metallicity relation found in this work resembles that obtained for wide binary systems, although in our sample of post-common-envelope binaries, there is a lack of old and high-metallicity objects, probably due to an observational bias against these objects.

Although preliminary, the results presented in this work highlight important trends. A more detailed and robust analysis, incorporating a complete sample along with the effects of selection biases and prior assumptions, will be necessary to fully confirm and refine these findings.

Acknowledgments

We acknowledge the positive and valuable comments from our anonymous referee. We acknowledge support from MINECO under the PID2020-117252GB-I00 and PID2023-148661NB-I00 grant and by the AGAUR/Generalitat de Catalunya grant SGR-386/2021.

References

Appendix A: Validation input space parameters

Fig. A.1. Inverse population synthesis validation test: corner plots (left panels) and correlation matrices (right panels) corresponding to the three synthetic binary systems with final periods P = 1.16, P = 15.5, and P = 41.73 hours (top, middle, and bottom panels, respectively). Corner plots represent the set of solutions found by our algorithm within a 1% error threshold.

Appendix B: Observed stellar parameters

Appendix C: Stellar input parameters

Fig. C.1. Input space parameter obtained through our inverse population synthesis algorithm for three observed cases: ZTF J041016.82-083419.5, ZTF J063808.71+091027.4, and ZTF J122009.98+082155.0 (top, middle, and bottom panels, respectively.

Tables C.2 and C.3 correspond to the mean (average) and mode (most probable value) solution for the stellar input parameters of the sample analyzed in this work. It should be recalled that, due to the correlation between the parameters, the set of input values for a certain system does not guarantee that it should be a solution. That is, after evolving the system with these initial set of parameters, the output values can be in disagreement with the observed values (see Section 5).

All Tables

All Figures

	Fig. 1. Schematic representation of the problem. The space that contains all possible initial values (Σ) is projected, after binary evolution, onto the output space (Ω). Only those points within S lead to a final set of parameters, called O, that are at a distance, d < σ, from a certain set of observed parameters. See text for details.
In the text

	Fig. 2. Schematic representation of our adaptive random walk algorithm. See text for details.
In the text

	Fig. 3. Parameter comparison of output values for the solution that provides the minimum distance, (M1, M2, Teff, 1, Teff, 2, P)out, with respect to the observed values, (M1, M2, Teff, 1, Teff, 2, P)obs. As it can be seen, all parameters perfectly agree except Teff, 2, which shows a clear shift between the output versus observed values. See text for details.
In the text

	Fig. 4. Distribution of parameters for the entire sample.
In the text

	Fig. 5. Common envelope efficiency parameter, αCE, as a function of the mass of the secondary. The decreasing linear trend is marked by the dashed red line.
In the text

	Fig. 6. Initial eccentricity distribution for the systems analyzed in this work. Despite the scarcity of the data, the distribution seems to be consistent with the canonical thermal distribution (red line), f = 2e, although some initially already circularized systems appear in the sample.
In the text

	Fig. 7. Age-metallicity relationship for the objects analyzed in this work (black points). For visual comparison we also show the linear trend of our sample (dashed red line) and the sample of wide WDMS binaries from (Rebassa-Mansergas et al. 2021, gray dots, RM+21).
In the text

	Fig. A.1. Inverse population synthesis validation test: corner plots (left panels) and correlation matrices (right panels) corresponding to the three synthetic binary systems with final periods P = 1.16, P = 15.5, and P = 41.73 hours (top, middle, and bottom panels, respectively). Corner plots represent the set of solutions found by our algorithm within a 1% error threshold.
In the text

	Fig. C.1. Input space parameter obtained through our inverse population synthesis algorithm for three observed cases: ZTF J041016.82-083419.5, ZTF J063808.71+091027.4, and ZTF J122009.98+082155.0 (top, middle, and bottom panels, respectively.
In the text