Barium and related stars, and their white-dwarf companions. III. The masses of the white dwarfs

Masses are one of the most difficult stellar properties to measure. In the case of the white-dwarf companions of Barium stars, the situation is worse. These stars are dim, cool, and difficult to observe via direct methods. However, Ba stars were polluted by the Asymptotic Giant Branch progenitors of these WDs with matter rich in heavy elements, and the properties of their WD companions contain key information about binary interaction processes involving AGB stars and about the slow-neutron-capture(s)-process of nucleosynthesis. We aim to determine accurate and assumption-free masses for the WD companions of as many Ba stars as possible. We want to provide new observational constraints that can help us learn about the formation and evolution of these post-interaction binary systems and about the nucleosynthesis processes that took place in the interiors of their AGB progenitors. We combined archival radial-velocity data with Hipparcos and Gaia astrometry using the software package orvara, a code designed to simultaneously fit a single Keplerian model to any combination of these types of data using a parallel-tempering Markov chain Monte Carlo method. We adopted Gaussian priors for the Ba star masses and for the parallaxes, and assumed uninformative priors for the orbital elements and the WD masses. We determined new orbital inclinations and companion masses for 60 Ba star systems, including a couple of new orbits and several improved orbits for the longest-period systems. We also unravelled a triple system that was not known before and constrained the orbits and the masses of the two companions. (Continued in the manuscript)


Introduction
About half of the elements heavier than iron are synthesized by the slow neutron capture (s-) process of nucleosynthesis (e.g. Burbidge et al. 1957;Clayton et al. 1961;Käppeler et al. 2011). The main astrophysical site that meets the appropriate conditions for the s-process to operate is the helium-rich intershell in the interiors of thermally pulsing Asymptotic Giant Branch (tp-AGB) stars (e.g. Lugaro et al. 2003b;Cristallo et al. 2009;Karakas 2010;Käppeler et al. 2011). However, the overabundance of s-process elements on the surface of a star is not a unique feature of AGB stars. Barium (Ba) stars are an example of s-process enriched objects that have not reached the tp-AGB phase yet. They are known to form when an AGB companion pollutes them in a binary system (e.g. McClure et al. 1980;Mc-Clure 1984;Udry et al. 1998a;Jorissen et al. 1998). The mass donors in these systems evolved off the AGB long ago and are now dim white dwarfs (WD), while the accretors -the Ba stars -are observed on the main sequence (e.g. North & Duquennoy 1991;Jorissen & Boffin 1992;North et al. 1994North et al. , 2000Pereira 2005; Kong et al. 2018;Escorza et al. 2019b), the red-giant (e.g. Bidelman & Keenan 1951;McClure 1983;Udry et al. 1998b;Jorissen 2004;Escorza et al. 2017;Jorissen et al. 2019), and the AGB (as extrinsic S stars, e.g. Jorissen et al. 1998Jorissen et al. , 2019Shetye et al. 2020) phases.
Although their exact formation channel and the masstransfer mechanisms involved are not well understood (e.g. Tout & Eggleton 1988;Han et al. 1995;Soker 2000;Pols et al. 2003;Bonačić Marinović et al. 2008;Izzard et al. 2010;Dermine et al. 2013;Abate et al. 2018;Saladino & Pols 2019;Gao et al. 2023), our knowledge about the spectroscopic orbital parameters of Ba star systems and about the stellar properties of the Ba stars themselves is generally well established (e.g. Escorza et al. 2019b;Jorissen et al. 2019, and references therein). Additionally, the evolutionary link between dwarf and giant Ba stars is well accepted (e.g. ). However, not much is known about the WD companions. The mass-function distribution of Ba star systems is consistent with a narrow distribution of companion masses peaking at 0.6 M (e.g. Webbink 1986;McClure & Woodsworth 1990;Jorissen et al. 1998;Merle et al. 2016;Jorissen et al. 2019;Escorza et al. 2019a), but very few absolute masses have been determined, since there is normally no information about the orbital inclinations of these systems (a few exceptional cases were published by Pourbaix & Jorissen 2000;Escorza et al. 2019b;Jorissen et al. 2019, among others, by combining the orbital parameters of Ba stars with Hipparcos astro-Article number, page 1 of 46 arXiv:2301.04232v2 [astro-ph.SR] 18 Jan 2023 A&A proofs: manuscript no. main metric data). These WDs are cool, dim, and directly undetectable in most cases; although, Böhm-Vitense et al. (1984; Gray et al. (2011), among others detected UV excess flux attributable to the WD in a few Ba star systems.
The masses of the WD companions of Ba stars contain important information about the AGB progenitors and the nucleosynthesis processes that took place in their interiors, and they are important input for binary interaction models. Even though mixing and dilution processes such as thermohaline mixing (e.g. Proffitt & Michaud 1989;Charbonnel & Zahn 2007;Stancliffe et al. 2007;Stancliffe & Glebbeek 2008;Aoki et al. 2008), rotationally induced mixing (e.g. Denissenkov & Tout 2000), or atomic diffusion (e.g. Matrozis & Stancliffe 2016 might impact the final level of s-process abundance on Ba stars, correlations between these abundances and the WD mass can give us observational information about the efficiency of the s-process at different masses and metallicities and help us constrain AGB models (e.g. Cseh et al. 2022) and mass-transfer and dilution models (e.g. Stancliffe 2021). The ratio between the amount of heavy s-process elements (hs), such as Ba, La, or Ce, and light s-process elements (ls), such as Sr, Y, or Zr, on the surface of Ba stars suggests that the material accreted by these stars was synthesized by low-mass AGB stars (< 3 M ; Lugaro et al. 2003aLugaro et al. , 2012Lugaro et al. , 2016Cseh et al. 2018;Karinkuzhi et al. 2018), which still needs to be confirmed by measuring these WD masses. Additionally, Jorissen et al. (2019) suggested that WD companions of strong Ba giants (based on the Ba index introduced by Warner 1965) are more massive on average than the WD companions of mild Ba stars. However, most of their masses were determined under the assumption of a constant (or very narrow distribution of) Q = M 3 WD /(M Ba + M WD ) 2 as proposed by Webbink (1988) and McClure & Woodsworth (1990), so this trend still needs to be confirmed with assumption-free measurements of WD masses.
In the first two papers of this series, Jorissen et al. (2019) and Escorza et al. (2019b) collected old and new radial-velocity (RV) data to study the orbits of giant and dwarf Ba stars, respectively. Additionally, we used spectroscopically-determined stellar parameters and Gaia DR2 distances (Lindegren et al. 2018;Bailer-Jones et al. 2018) to locate these stars on the Hertzsprung-Russell diagram (HRD). By comparing their location on the HRD with STAREVOL evolutionary tracks (Siess et al. 2000;Siess 2006Siess , 2008 and following the methodology described in Escorza et al. (2017), we also determined accurate masses for the primary stars of these systems, the Ba stars. In this third article, we focus on the faint WD companions. We used the orvara software (Brandt et al. 2021c) to combine all the radialvelocity data available, the astrometric measurements from the Hipparcos mission (Perryman et al. 1997), the Gaia positions and proper motions (Lindegren et al. 2021), and the information in the Hipparcos-Gaia Catalogue of Accelerations (HGCA; Brandt 2018Brandt , 2021 to determine the astrometric orbital parameters of as many Ba star systems as possible (see Sect. 2 for the description of the sample), and then derive the mass of the secondary stars. All these data sets are described in Sect. 3. An important improvement with respect to what has been attempted before for these objects is that we use a joint astrometric-spectroscopic model (see Sect. 4) to find new best-fitting orbital parameters instead of relying only on RV data or imposing the spectroscopic solution on the astrometric data. Our results are presented in Sect. 5 and their implications are discussed in Sect. 6. We also discuss the feasibility of the direct detection of the WD companion for a subset of the longest-period systems in Sect. 7.

Target selection
For our methodology (see Sect. 4) to be applicable, a target must fulfil three requirements: (i) it must be part of the HGCA, (ii) we must have a good initial estimate of the mass of the primary star in the system, and (iii) the Hipparcos solution cannot not be more complex than the 5-parameter solution. As a starting point, we selected all the Ba stars from the samples studied by Jorissen et al. (2019), Escorza et al. (2019b) and North et al. 2020 that have Hipparcos identifiers. We excluded confirmed triple systems, stars whose Ba star nature was not certain or is under current investigation (e.g. Escorza et al. 2023), and a few systems that had an acceleration solution or an orbital solution in the Hipparcos data reduction (solution types,Sn,equal to 7,9 and 75). We ended up with 60 systems. Table 1 presents our target list. In addition to the most commonly used identifier, we include the Hipparcos identifier of each system and the Ba star type. We distinguish between pre-RGB, which are all the stars classified as dwarfs or subgiants by Escorza et al. (2019b) and North et al. (2020), and mBag or sBag which are stars classified as mild or strong Ba giants by Jorissen et al. (2019) based on their [La/Fe] and [Ce/Fe] values (as measured by Smith 1984;Allen & Barbuy 2006a,b;Pereira et al. 2011;Karinkuzhi & Goswami 2014Luck 2014;de Castro et al. 2016;Merle et al. 2016; Van der Swaelmen et al. 2017;Karinkuzhi et al. 2018;Jorissen et al. 2019) and on the Ba index introduced by Warner (1965). The table also lists the Ba star masses (M Ba ) that we used as a prior in our MCMC model (see Sect. 4) and the metallicity of the system, both values collected from Jorissen et al. (2019), Escorza et al. (2019b) or North et al. (2020) unless explicitly specified. For this work, we recomputed the primary masses for the ten systems that were part of the non-single-star (NSS) Gaia DR3 catalogues (Gaia Collaboration et al. 2022). We followed the exact same procedure followed and described in the mentioned papers and used the same STAREVOL grids of models (Siess et al. 2000;Siess & Arnould 2008;Escorza et al. 2017), but we used the NSS Gaia DR3 parallaxes to recalculate their luminosities and masses. Finally, the last column of Table 1 includes the sources where we found the archival RV data used in our analysis.

CORAVEL, HERMES and other radial-velocity data
The most important radial-velocity monitoring programs of Ba stars were carried out with the two CORAVEL spectrometers and with the HERMES high-resolution spectrograph. The CORAVEL spectrometers (Baranne et al. 1979) were installed on the 1-m Swiss telescope at the Haute-Provence Observatory and on the 1.54-m Danish telescope at ESO -La Silla, while HERMES (Raskin et al. 2011;Raskin & Van Winckel 2014) is mounted on the 1.2-m Flemish Mercator telescope at the Observatory El Roque de Los Muchachos.
The main results of these radial-velocity programs were published by Jorissen & Mayor (1988); Jorissen et al. (1998);Udry et al. (1998a,b); North et al. (2000); Gorlova et al. (2013); Jorissen et al. (2019); Escorza et al. (2019b) among others, and the strength of combining the two data sets, particularly for the longest-period systems, was discussed in the last two mentioned papers. Jorissen et al. (2019) and Escorza et al. (2019b) also described the data reduction process for the two instruments and the existence of a non-zero radial-velocity offset between the data sets due to the use of a different system of standard stars. Table 1: List of Ba star systems to which our methodology was applied. Column 1 lists the most commonly used identifiers, while column 2 lists the Hipparcos identifiers. Column 3 lists the Ba star type, which can be preRGB for stars classified as dwarfs or subgiants, or mBag or sBag for stars classified as mild or strong Ba giants, respectively. Column 4 lists the primary star masses and column 5, the metallicity of the system. These values were derived or collected by Escorza et al. (2019b) or Jorissen et al. (2019) for preRGB and giant systems, respectively, unless otherwise indicated. Finally, the last column gives the sources of the archival RV data we used.  Escorza et al. (2019b) determined a relation between the offset and B-V by comparing old and reprocessed CORAVEL data and calculated a fixed offset for each studied Ba star. For this work, we combined the two approaches. Where the RV data of a specific instrument spanned over a full orbit or more, we treated the offset as an additional free parameter that was optimized during the orbital fitting process. However, for systems with very few HERMES points or for some very long orbits, the offsets from Jorissen et al. (2019) or Escorza et al. (2019b) were adopted and fixed. This will be clearly indicated in the captions of each RV fit shown in Appendix A. Future monitoring with HERMES would remove the need for this assumption, allowing us to fit the offset term directly.
To complement the main CORAVEL and HERMES data, we collected additional radial-velocity measurements from other works and instruments, and the sources are listed in Table 1. An optimizable RV offset, such as the one described above between CORAVEL and HERMES, was considered for each data set.

Hipparcos astrometric data
The Hipparcos satellite ESA (1997), launched in 1989, was the first space mission with precision astrometry as its main goal. Between 1989 and 1993, Hipparcos measured the location and motion on the sky of more than 100,000 stars many times, to figure out their astrometric path. For each target in Table 1, we used the positions and the proper motions from the Hipparcos Cata-logue (Perryman et al. 1997). Additionally, we also queried the individual astrometric measurements from the re-reduction of the Hipparcos intermediate astrometric data (IAD;van Leeuwen 2007). The coordinates are expressed in the International Celestial Reference Frame (ICRF) at the 1991.25 epoch.
Since the code we are using is not yet prepared to deal with Hipparcos solutions more complex than the 5-parameter solutions, we excluded a few targets with acceleration or orbital solutions from the initial sample. Some of our remaining targets have a stochastic Hipparcos solution (Sn = 1). These represent cases where the residuals are significantly larger than expected, but since the proper motions and the IAD were obtained using a 5-parameter solution, we included them and gave them no special treatment.

Gaia astrometric data
The Gaia mission (Gaia Collaboration et al. 2016, 2018 was launched in 2013 as a successor of Hipparcos. For now, none of the Gaia Data Releases (DR) published individual astrometric measurements, so we queried the positions and proper motions published for our targets in the Early DR3 catalogue (Lindegren et al. 2021). In contrast with the Hipparcos data, these are expressed in the ICRF at the 2016 epoch. The Gaia EDR3 parallaxes were also queried and used as prior in the fit (see Sect. 4). Finally, in order to use an equivalent to epoch astrometry, we also used the Gaia Observation Forecast Tool (GOST 1 ). The GOST provides the predicted observations and scan angles for any Gaia source. We note that not all the planned observations will be used in the final astrometric solution, since some predicted scans might correspond to satellite dead times or might be unusable or rejected as outliers. For example, up to 20% of the observations predicted by GOST were excluded from the analysis published in Gaia DR2 (Brandt et al. 2021b).
Ten of the 60 targets presented in this study had a nonsingle-star (NSS) solution in Gaia DR3 (Gaia Collaboration et al. 2022). These targets are: HD 50264, HD 207585, HD 221531, HD 34654, HD 49841, HD 199939, HD 224621 and HD 87080, which had a non-single-star solution compatible with a combined astrometric and single lined spectroscopic model, and HD 44896 and HD 123585, which had a solution compatible with an astrometric binary. For these targets, we used the Gaia DR3 NSS parallax as priors, instead of the EDR3 value. Even though the Gaia DR3 NSS catalogue provided orbital inclinations for these 10 systems, we decided not to include an inclination prior in our calculations to first, treat all systems equally, and second, compare our independently determined inclinations with the new Gaia ones and validate our method.

The Hipparcos-Gaia Catalogue of Accelerations.
As an additional astrometric constraint, we used the difference in Hipparcos and Gaia proper motions via the Hipparcos-Gaia Catalogue of Accelerations (HGCA; Brandt 2018Brandt , 2021. This catalogue puts the Hipparcos, Gaia, and Hipparcos-Gaia (H-G) proper motions into the same reference frame to make them suitable for orbital fitting. The Hipparcos-Gaia proper motion is derived from the right ascension and declination measurements from the two missions and is by far the most precise due to the long time elapsed between them (proper motion uncertainties scale inversely with the time baseline). This acceleration in the inertial frame can be used to improve the dynamical parame-ters of the companion and to measure its mass because it breaks the mass-inclination degeneracy that RV data suffers from. We used the EDR3 version of the HGCA (Brandt 2021) for all our targets.
The EDR3 version of the HGCA also provides a χ 2 value between the two most precise proper motion measurements (normally EDR3 and H-G). This value is meant to find accelerating candidates for follow-up and if it is higher than ∼11.8 (Brandt 2021), the measured acceleration is considered significant and statistically different, by 3σ, from constant proper motion. In our case, since all our targets are known binaries, we do not need this χ 2 value to detect accelerators, but it can give us a hint about which systems are truly benefiting from the HGCA measurement. The queried HGCA χ 2 values are included in the last column of our result table (Table 2).

Orbital analysis with orvara
Orvara, developed by Brandt et al. (2021c), is designed to simultaneously fit a single Keplerian model to any combination of radial velocity data and relative and absolute astrometry. The combination of these different data sets, using Orvara or not, has recently proven to be very powerful to improve the accuracy of orbits and to measure precise companion masses, even for very long period systems where the observations only cover part of the orbit (e.g. De Kervella et al. 2020;Brandt et al. 2021c;Venner et al. 2021;Franson et al. 2022;Brandt et al. 2021a;Leclerc et al. 2022).
Orvara integrates the Hipparcos and Gaia intermediate astrometry package (htof; Brandt et al. 2021b) to fit the Hipparcos epoch astrometry and the times and scan angles of individual Gaia epochs. The code uses a parallel-tempering Markov chain Monte Carlo method (ptmcmc, Foreman-Mackey et al. 2013) and first fits the RV data. Orvara allows RV points from each instrument to have a different RV zero point, which we need at least for the CORAVEL-HERMES combination as discussed in Sect. 3.1. Then the absolute astrometry is included and fit for the five astrometric parameters (positions, α and δ, proper motions, µ α and µ δ , and parallax, ) using htof at each MCMC step. On top of the five astrometric parameters, we fitted the six Keplerian orbital elements (semimajor axis, a, eccentricity, e, time of periastron passage, T 0 , argument of periastron, ω, orbital inclination, i, and longitude of the ascending node, Ω), the masses of the two components (M Ba and M WD ), and a radial-velocity jitter per instrument to be added to the uncertainties. Note that while the difference between the Hipparcos and Gaia reference frames is taken into account in the HGCA, this is not the case for the IAD. However, the rotation difference in the proper motions is w = (−0.120, 0.173, 0.090) mas/yr (Fabricius et al. 2021). These values are very small compared to the amplitudes of the proper motion curves that we are measuring (of the order from a few to a couple of tens mas/yr, see Appendix A), and smaller than the residuals to these fits in most cases, so we did not take them into account.
For this work, we assumed uninformative priors for the orbital elements and for the WD mass, but we adopted Gaussian priors for the primary mass and for the parallax. For each target, we used the M Ba value given in Table 1 but using three times the error bar as sigma to be conservative and take into account systematic errors not accounted for in the statistical uncertainty. Concerning the parallax, the Gaia EDR3 value was used as prior for most targets, and the Gaia DR3 NSS parallax was used for the 10 targets with a NSS solution. We used 15 temperatures and for each temperature we use 100 walkers with 100,000 steps per walker. In a few cases, we needed to run twice as long or repeat the calculations using an educated starting position based on our knowledge about the systems from the RV-only fits published by Jorissen et al. (2019) or Escorza et al. (2019b), however, in most cases, the MCMC chains converged quite quickly. We discarded the first 300 recorded steps (the first 15000 overall, as we saved every 50) as the burn-in phase to produce the results presented in Sect. 5.
For more details about the computational implementation in orvara and htof and for case studies showing the performance of the code, we refer to the mentioned publications. Table 2 lists the obtained astrometric-spectroscopic orbital parameters, the best-fitting WD masses, the χ 2 of the best fit, and the HGCA χ 2 values discussed in Sect. 3.4. To make the table easier to read, we assume that the error bars we obtained from the MCMC fit are symmetric and listed only the largest value. This means that in some cases, the table lists an overestimated uncertainty in one of the two directions. The χ 2 values are an overall absolute astrometric χ 2 , computed adding the χ 2 for the Hipparcos proper motions (χ 2 H ), the χ 2 for the long-term Hipparcos-Gaia proper motions (χ 2 HG ), and the χ 2 for the Gaia proper motions (χ 2 G ). orvara uses RV jitter terms such that the reduced χ 2 of the RV fit is 1, so we did not take it into account to evaluate the goodness of the fit.

Results
The table is ordered based on the orbital period, with the systems with the longest periods first. This way, we can notice that all the systems with periods longer than ∼ 3 years have significant astrometric accelerations according to their HGCA χ 2 values, while most of the systems below that threshold do not. Finding such a clear threshold in a sample of confirmed binaries is an indication of the type of systems that the HGCA can help identify.
In addition to the table and in order to illustrate and discuss how the results that we get from orvara look like, we include the results for the main-sequence Ba star HD 2454 in Fig. 1. HD 2454 was first identified as a Ba dwarf by Tomkin et al. (1989), and North et al. (2000) confirmed its binarity even though they did not have enough data to estimate its orbital period. More recently, Gray et al. (2011) found direct evidence of the presence of a WD companion in the system thanks to the Galaxy Evolution Explorer (GALEX; Martin et al. 2005) UV observations and, since 2011, HD 2454 has been part of the longperiod binary monitoring program carried out with the HER-MES spectrograph (see Sect. 3.1). In spite of having almost three decades of RV data between the CORAVEL and HERMES measurements, Escorza et al. (2019b) were not able to constrain the orbit either. However, combining all these RV data points with the Hipparcos and Gaia information, we can finally estimate the orbital elements of HD 2454 as well as the mass of its WD companion. Fig. 1 shows, on the top left panel, the astrometric orbit of HD 2454, including the predicted position of the companion on the scheduled date of Gaia DR3. The best-fitting orbit is plotted as a black thick line, while 40 other well-fitting orbits are colour-coded as a function of the companion mass. On the top right panel, we show the RV curve of HD 2454. For this target we had CORAVEL (orange circles), SOPHIE (pink diamonds), and HERMES (green triangles) RV data. The plot shows that leaving the RV offsets between instruments completely free produces families of solutions with similar orbits and masses but different RV offsets (displaced vertically in the RV plot). This is especially noticeable in cases like this one, where no data sets covers even half an orbit. We want to note that even though we left the RV offsets free in most cases, we always made sure that the best-fitting solution required reasonable values and, especially in the CORAVEL-HERMES case, that these values were close to the values obtained by Jorissen et al. (2019) and Escorza et al. (2019b).
The two bottom panels of Fig. 1 show the fit to the proper motions in the right ascension (left) and declination (right) directions, as measured by Hipparcos (squared data point) and Gaia (circular data point). All the data sets included in the figures were fitted at the same time, and the plotted models are the same in all plots. Finally, Fig. 2 shows the one and two-dimensional projections of the posterior probability distributions of the masses of the two components in the system and a few orbital parameters (semi-major axis, eccentricity, and inclination) from the joint RV and astrometric MCMC computations. This corner plot shows that the two masses are correlated, and that the semimajor axis is also correlated with the total mass of the system. These correlations are even stronger for other targets.
We have included in Appendix A figures similar to Figs. 1 and 2 for all the targets in our sample. Additionally, an individual case of study of a Ba dwarf using the same method was presented in Escorza & De Rosa (2022).

Spectroscopic orbital parameters
Even though the main goal of this work was deriving the masses of the WD companions of all these Ba stars, an important additional result of this new method are the new orbits of HD 2454 and BD-11 o 3853, which could not be constrained before, as well as the improved orbits of a few other long-period systems. When comparing the orbital periods obtained using orvara to those presented in Jorissen et al. (2019), Escorza et al. (2019b) and North et al. (2020), which were obtained by fitting only the RV data, we get a very tight relation. The purely spectroscopic parameters and the new parameters are consistent with each other within error bars in almost all cases, and we discuss the exceptions below.

HD 218356
Our first orbital fit for this system converged to a period of more than 40 years, while the period published by Griffin (2006) and Jorissen et al. (2019) for HD 218356 was 111 days. No third object has been detected in this system in the past, but the mild s-process enhancement in the visible star has been flagged as surprising given the close orbit. We performed a three-body fit, setting strong constraints on the inner orbit using the published spectroscopic parameters, and we succeeded to recover the orbital parameters of two companions, confirming that HD 218356 is actually a triple system with a third companion in a much longer orbit than the published period. The orbital parameters of the system are included in Table 3 and the combined RV fit can be seen in Fig. 3. In order to test the significance of this detection, we compared the Akaike Information Criterion (AIC) of the two-and three-component models using the radvel package (Fulton et al. 2018). We found a ∆AIC of 439 favouring the three-component model. Given the masses of the two companions, we expect the WD that polluted the Ba star to be in the outer orbit. This would also explain the mild s-process enhancement reported for HD 218356. We included the corner plots with Table 2: Orbital elements and WD masses derived following the method described in Sect. 4 and listed in order of decreasing periods. The columns list, in order, the most commonly used identifier, the orbital period P, the eccentricity e, the time of periastron passage, the absolute semimajor axis of the orbit a, the argument of periastron of the visible star ω * , the longitudes of the ascending node Ω, the orbital inclination i, and the WD companion mass M WD . To keep the table cleaner, we assumed symmetric error bars, and included only the largest one of the two. The last two columns include the χ 2 of the best fitting model and the HGCA χ 2 value discussed in Sect. 3.4.  Table 3 for the remaining parameters.

Star ID
Article number, page 6 of 46 A. Escorza and R. J. De Rosa: Barium and related stars, and their white-dwarf companions III. The masses of the white dwarfs Table 2 continues  the parameters of both orbits in Appendix B. Only the outer orbit information is listed together with the other WD orbits in Table  2.

HD 201657
Our orbit fit for HD 201657 converged to twice the published orbital period and to a much more eccentric orbit. The astrometric data favours the longer orbit, and the RV data is not very constraining since we have only 15 CORAVEL points and one HERMES point. However, given the eccentricity-period diagram of Ba stars, the orbit published by Jorissen et al. (2019), the least eccentric of the two, is the most likely. We attempted to recover this orbit in order to check the quality of such a fit and calculate the WD companion mass by including an orbital eccentricity prior of 0.15 ± 0.15. We recovered Jorissen et al. (2019)'s orbital solution, although with a slightly higher χ 2 for the astrometric data. Since we considered this solution more likely for a Ba star, we listed this orbit in Table 2, but we show both fits and corner plots in Appendix C. More HERMES data would be essential to solve this case.

Astrometric orbital parameters
Finally, in addition to the new and improved orbital parameters, this method provided us with orbital inclinations for all these Ba star systems. Fig. 4 shows the distribution of the obtained cos(i) values. This distribution should be flat if we could assume our sample of binaries is randomly distributed on the sky, and even though we only have 60 systems, the distribution is compatible with a uniform one. We performed a Kolmogorov-Smirnov (KS) test (e.g. Press et al. 1986), and we obtained p-values higher than 0.8 when comparing our cos(i) distribution with uniform distributions of the same sample size.
The new orbital parameters are also compatible with the astrometric parameters published in Gaia DR3 for the ten targets available in their catalogue. Concerning the periods, all Gaia DR3 values are consistent with our values within 2σ. The largest difference is found for HD 221531, for which Gaia DR3 published a period of 1668 ± 135 days, about 260 days longer than our period. The Gaia DR3 time span is about 1000 days (Gaia Collaboration et al. 2022), while our data covers a few decades in most cases. Hence, we think that our method is more reliable to obtain the periods of long-period binaries. The eccentricities are compatible as well, without significant exceptions, and finally, we used the Thiele-Innes elements published in the Gaia DR3 catalogue and followed Halbwachs et al. (2022) to compute the orbital inclinations of these systems from the Gaia DR3 data. The Gaia DR3 inclinations are also compatible with the inclinations we obtained with our full RV+astrometric model within 1.2 times our σ. As discussed above, the HGCA is not very constraining for systems with periods below about 3 years, so while we think our method is better to determine the orbital periods of Ba stars, the Gaia DR3 inclinations are probably of better quality than ours for the shorter-period systems. When the epoch astrometry of the Gaia mission is published, we will be able to combine these data with all our other data sets and improve our results for the shortest period systems. Table 2 lists the masses we obtained for the companions to all the Ba stars in our sample, and Fig. 5 shows the distribution of these masses as a purple dashed histogram. Also in Fig. 5 we compare this new distribution to the distribution obtained by Jorissen et al. (2019) and Escorza et al. (2019a) for the same stars, which is drawn in black. The insert in the figure shows the cumulative frequency of the two distributions, including an envelope with the 1 − σ uncertainty for our distribution, which also envelopes the old distribution. We obtained a p-value of 0.010 on a KS test, which is not low enough to reject the null hypothesis. The two distributions are not statistically different. In Fig. 6, we plot the mass distributions of the companions to strong and mild Ba giants, separately. As mentioned in the introduction, this distinction is made based on the abundance ratios [La/Fe] and [Ce/Fe] and on the Ba index introduced by Warner (1965). We do not include the pre-RGB stars in this comparison, because the distinction between strong and mild enhancement has not been as clearly established as it has for the giants. We note that the WDs occupying the high-mass tail belong to systems with strong Ba giants. However, we performed a KS test, and we obtained a p-value of 0.45, meaning that we cannot reject that the two samples are drawn from the same distribution. The cumulative distributions plotted in the insert also show that  There are a few individual systems that appeared as clear outliers or that even have WDs with unphysical masses. These are briefly discussed below. There are two systems for which our simulations converged to very low WD masses. These are HD 18182 and HD 95241. The fit we achieved for the former is less than ideal (see Figs. A.14 and A.16), and even though the mass is small, taking the error bars into account, the value is compatible with an average WD in our sample. The CORAVEL RV data is not very constraining and the HERMES points, being of much higher quality, still fall on the same range of orbital phases, covering in total less than half of the orbit. Additionally, the Hipparcos and Gaia proper motions in the right ascension direction are very similar, not adding strong constraints to the fit either. This WD mass should be taken with caution.

White Dwarf masses
The fit for HD 95241, on the other hand, is significantly better. We used 97 RV points that cover very well the whole orbit (see Fig. A.21) and obtained clean and symmetric posterior distributions (see Fig. A.23). Of course, M Ba and M WD are very strongly correlated, so if the M Ba prior was incorrect, too small in this case, it would directly affect M WD . The mass of HD 95241 was determined by Escorza et al. (2019b) by comparing the location of the star on the HR diagram with STAREVOL (Siess et al. 2000;Siess & Arnould 2008) evolutionary tracks. The stellar parameters were determined from HERMES high-resolution highsignal-to-noise spectra and are in agreement with other studies (e.g. Takeda 2007;Soubiran et al. 2016). However, HD 95241 was flagged as a mild Ba dwarf by Edvardsson et al. (1993) having only a marginal overabundance of s-process elements with respect to iron. Other Ba dwarf candidates of their sample have been proven to be wrongly flagged. Most of them are likely single stars (Escorza et al. 2019b). It is possible that HD 95241 has a low-mass companion that is not a WD, and if it is a WD, its AGB progenitor was not massive enough to reach the thermally pulsing AGB phase and produce s-process elements. HD 95241 is likely not a Ba star and will be removed from further analysis.

The most massive WDs: HD 49641 and HD 31487
On the high-mass end of the distributions, there are two systems with WD masses clearly outlying from the initial mass distribution (M WD ≥ 1.2 M ). These are HD 49641, with M WD = 1.2 ± 0.4 M , and HD 31487, with M WD = 1.59 ± 0.22 M . The fit for HD 49641 is not very good, because the available RV data was scarce and old, so one should take this WD mass with caution, but the fits for HD 31487 seems reliable, including a clean result for the orbital projection on the sky (see Fig. 7). In order to try to explain this last mass, one could again try to invoke a wrong M Ba prior. We used the primary mass determined by Karinkuzhi et al. (2018). The primary mass listed by Jorissen et al. (2019) is not in agreement with Karinkuzhi et al. (2018)'s within error bars, but we decided to use the latter after studying their HR diagram (their Fig. 16). In any case, Jorissen et al.
(2019)'s mass is higher, and would result in a higher companion mass. Karinkuzhi et al. (2018)'s value seems reasonable given the location of the star on the HR diagram, and it is a very average value for giant Ba stars. Additionally, there is no big discrepancy between the parallaxes published in the different Gaia Data Releases. While a wrong parallax could have led to a wrong luminosity, hence mass, determination, we have no good reason to doubt this mass. From the posterior distributions and 1D projections shown in Fig Only with the dynamical information that we currently have, it is difficult to confirm that this 'massive companion' is a single object, and not a close pair formed, for example, by a faint mainsequence star and a WD (see van den Heuvel & Tauris 2020 for an example of such a situation). The strong s-process enhancement strongly suggests that there is a WD in the system, but since we cannot be certain of its mass, HD 31487 will be removed from further discussion.

Mass distributions
The mass distribution that we obtained for the WD companions of Ba stars is compatible with current estimates for field WD masses. The average mass of DA WDs (WDs with only Balmer lines in their spectra) is about 0.60 M , while that of DB WDs (WDs with no H or metals lines in their spectra, only helium lines) is 0.68 M (Kleinman et al. 2013). The weighted average of our mass distribution is 0.65 M , after removing the two targets mentioned in Sect. 5.3. There is a high-mass tail present in the mass distribution of WDs orbiting Ba giant that Jorissen et al. (2019) and Escorza et al. (2019a) already discussed (see also Fig. 5).
In order to evaluate if Q = M 3 WD /(M Ba + M WD ) 2 is constant, we computed this value for all our targets and present the average and the standard deviation for each one of the three subsamples separately in Table 4. The new distributions are marginally different to literature Q distributions (see Table 1 in Escorza et al. 2019a). We obtained p = 0.048 for the strong Ba giants, p = 0.035 for the mild Ba giants and p = 0.012 for the Ba dwarfs, when we performed KS tests. The main difference is that the new distributions are not as narrow as obtained in the past when modelling f (m) = Q sin 3 i, with f (m) being the spectroscopic mass function. In order to check if this is caused by the fact that the individual inclination uncertainties play a role now, while an inclination distribution was assumed in the past, we calculated new Q distributions removing the 10 and 20% systems with larger uncertainties. All the observed distributions are broader than the literature ones, but not significantly different.
In Table 4, we have also included the average and standard deviations of the current mass ratios of the three Ba star subsamples. The two subsamples of giants show closer values, with the

A comment on nucleosynthesis predictions
It is difficult to make a direct correlation between the WD companion mass and the s-process enhancement of the Ba star because many parameters, apart from the WD progenitor mass, strongly affect the final Ba star abundances and the unknowns are still stronger than the observational constraints (see Cseh et al. 2022, for a study of abundances in individual Ba giant systems). For example, the efficiency of the mass transfer and the dilution factor, the ratio between the accreted mass and the mass in the Ba star envelope over which this is mixed in, are major uncertainties in our understanding of the formation of Ba stars and will directly affect the final s-process enhancement (e.g. Stancliffe 2021). Of course, the efficiency of the s-process of nucleosynthesis in the interiors of AGB stars, which strongly depends on the mass and the metallicity of the star itself, is also a key parameter in order to explain a possible correlation between WD mass and Ba enhancement(e.g. Busso et al. 2001;Karakas & Lugaro 2016; Van der Swaelmen et al. 2017). Additionally, even the number of thermal pulses and third dredge-ups experienced by the AGB star before the mass transfer episode took place will have an effect on the final s-process abundance pattern (e.g. Shetye et al. 2018), as well as mixing and diffusion below the AGB star's convective envelope will (e.g. Goriely & Siess 2018). Standard stellar-evolution models do not predict solarmetallicity low-mass AGB stars undergo third dredge-ups (e.g. Cristallo et al. 2015;Karakas & Lugaro 2016). This limit can go down to 1 M at lower metallicities (e.g. Stancliffe et al. 2005;Lugaro et al. 2012;Fishlock et al. 2014). However, including different additional effects in the models can help, for example, Weiss & Ferguson (2009) showed that including some overshooting below the convective pulse, their models could make a 1 M AGB star undergo third dredge-ups. Additionally, Shetye et al. (2019Shetye et al. ( , 2021 found several low-mass AGB stars currently undergoing third dredge-ups and their models succeeded to reproduce the s-process overabundance including diffusive mixing at the bottom of the stellar envelope. Additionally, according to several studies, the AGB stars that polluted Ba stars need to have masses below 3 M to be able to reproduce their abundance ratios with models (Lugaro et al. 2003a(Lugaro et al. , 2012Cseh et al. 2018;Karinkuzhi et al. 2018). Figure 8 shows the relation between the metallicity (listed in Table 1) and the obtained WD masses (Table 2) for the preRGB Ba stars (orange circles), the strong Ba giants (blue squares) and the mild Ba giants (green triangles) in our sample. The figure shows an expected correlation between the Ba-type and the metallicity, caused by the fact that the efficiency of the s-process in AGB stars decreases as the metallicity increases (e.g. Cseh et al. 2018;Jorissen et al. 2019). However, there is no obvious correlation between the WD mass and the metallicity, even though the AGB mass directly affects the s-process efficiency as well. The least massive WDs are in systems with [Fe/H] < −0.1, in agreement with the models, and the most massive WDs accompany Ba giants of [Fe/H] between −0.4 and −0.2, with the three most massive WDs being in a strong Ba star systems.
Among our sample of 58 systems (after having removed (HD 95241 and HD 31487 from the WD sample), we do not find Ba stars with unexpectedly low mass companions. As discussed in Sect. 5.3, the companion mass for HD 18182 should be taken with caution, but all other Ba star systems have WDs of around or more massive than 0.5 M , meaning that their progenitors were AGB stars of around or more massive than 1 M . Note that to make such a statement, one needs to rely on initial-final mass relationships (IFMR). We used as a reference the relation published by Marigo et al. (2020Marigo et al. ( , 2022. Using the same relation, we can claim that a fraction of the AGB stars that polluted our sample of Ba stars were more massive than the expected 3 M limit, since we found that several WDs have masses around or higher than 0.8 M . This is the case even taking into account the kink that Marigo et al. (2020Marigo et al. ( , 2022 find for WDs of about 0.70 -0.75 M with carbon AGB progenitors. Most IFMRs (e.g. Weidemann 2000; Kalirai et al. 2008;Williams et al. 2009;Andrews et al. 2015;Cummings et al. 2016;El-Badry et al. 2018) flatten at around M WD ∼ 0.8 M , making stars with a wide range of initial masses accumulate at that WD mass. However, their progenitors are expected to have initial masses in the range between 3.5 and 5.5 M , hence more massive than what the Ba stars abundance ratios seem to indicate.
The presence of these massive WDs orbiting around both strongly and mildly polluted Ba stars presents important constraints, as well as an interesting challenge, for evolutionary and nucleosynthesis models. Future studies of these systems following the line presented by Stancliffe (2021) or Cseh et al. (2022), but using these new WD masses, might be able to tell us new things about AGB stars. We note that our error bars are significant and that these statements blur if we consider two or three sigma uncertainties. This will improve when we have NSS parallaxes to obtain more accurate Ba star masses and Gaia astrometric epochs to improve the RV+astrometry fit. Direct imaging observations could also help constrain the longest-period systems better (see Sect. 7).

Future observational prospects: direct imaging of these white dwarfs with SPHERE
The nearby (d 100 pc) Ba stars that are host to long-period (P 10 yr) companions are suitable candidates for high-contrast imaging observations to spatially resolve the companion. These observations would provide relative astrometric and photometric measurements between the WD and the Ba star host. A single measurement of the instantaneous angular separation between the components would constrain both the total semi-major axis, and thus the total system mass, and the inclination (unless the observation occurred when the companion was crossing the line of nodes). Photometric measurements of the companion could be used to estimate the bolometric luminosity of the companion which, in conjunction with the mass, can be compared to  (Table 1) and the obtained WD masses (Table 2) for the preRGB Ba stars (orange circles), the strong Ba giants (blue squares) and the mild Ba giants (green triangles) in our sample.
WD cooling models to estimate the age of the companion (see e.g. Bonavita et al. 2020 for the discovery and analysis of a WD companion around a K-type star with SPHERE and e.g. Gratton et al. 2021 for the study of a sample of Sirius-like systems, long-period main-sequence + white dwarf binaries).
We assessed the feasibility of spatially resolving the companion by comparing the predicted angular separation and flux ratio between the WD companion and the Ba host star to the performance of the VLT/SPHERE instrument (Beuzit et al. 2019). We filtered the sample to exclude systems with a median apoastron distance within 100 mas, calculated from the MCMC samples described in Section 4. For these systems, the companion would always be within the inner working angle of the instrument, and impossible to resolve with SPHERE. This filter resulted in a subsample of eight systems for which the companion will be at a projected separation of ρ > 100 mas at some point in its orbit.
The feasibility of direct detection also depends on the flux ratio between the WD companion and Ba star host. We estimate the H-band flux ratio for each MCMC sample using purehydrogen (DA) atmosphere mass-luminosity relations from Holberg & Bergeron (2006). We assigned an age to each MCMC sample at random from a uniform distribution between 10 6 and 10 10 yr to account for the unknown age of the WD. The model grid was linearly interpolated in (log t, M) to extract an absolute H-band magnitude. This was converted into a flux ratio using the parallax from the MCMC sample and the apparent H-band magnitude of the Ba star. We assumed the companion has negligible flux in the H-band relative to the Ba star, such that the catalogue H-band magnitude of the system can be entirely ascribed to the Ba star.
We converted the orbital elements to the angular separation between the WD and the Ba star host at the epochs 2023, 2024, and 2025 for each MCMC sample. The predicted angular separation and flux ratio for each sample was then compared to the SPHERE contrast curve given in Wahhaj et al. (2021). We accounted for the degradation in contrast performance for fainter stars (e.g. Jones et al. 2022) by scaling the contrast curve by the square root of the H-band flux ratio between the Ba star host and HR 8799, the star observed by Wahhaj et al. (2021) to measure the contrast curve.
The predicted separations in 2023 and H-band contrast for each of the eight systems are shown in Figure 9. There are six systems with a non-negligible probability of detection at this epoch; the others are too faint to be detected given the expected contrast curve. The majority of these systems exhibit a strong correlation between the separation at the 2023 epoch and the mass of the WD companion. This can partly be explained by the constraint provided by a direct measurement of the semi-major axis of the system, leading to a much more precise measurement of the total mass of the system.

Summary and conclusions
The WD companions of Ba stars contain important information about the formation of these chemically peculiar stars, about the binary interaction processes that these systems underwent in the past, and about the nucleosynthesis processes that took place inside their AGB progenitors. However, they are cool, dim, and generally not detected by direct methods, so they have not been studied in detail in the past. A few absolute masses had been determined before this work by combining the spectroscopic orbital parameters of these systems with Hipparcos astrometric data. However, most published masses for WD companions of Ba stars were computed by making assumptions on the relation between the masses of the two stellar components in these systems or on their orbital inclinations.
In this work, we used the software package orvara to combine radial-velocity data, Hipparcos and Gaia positions and proper motions through the Hipparcos-Gaia Catalogue of Accelerations, and astrometric epoch measurements from the Hipparcos mission, and determine the astrometric orbital parameters of 60 stars flagged as Ba dwarfs or giants. Using this method, we could constrain the orbits of two long-period systems that could not have been constrained before with RV data only, and we improved the orbital solution of a few other systems. Orbital inclinations were also determined for the first time for many of these systems, and finally, including a prior on the Ba star masses, we also derived the mass of the secondary stars in these systems. Finally, we discovered that HD 218356, one of the shortest period Ba star systems known, is actually a triple system. We determined the parameters of both the inner and outer orbits and the masses of the two components, and it is very likely that the WD companion that polluted HD 218356 is in the outer orbit, explaining the mild s-process enhancement of the Ba giant.
The WD mass distribution presented in this work includes all systems published by Jorissen et al. (2019), Escorza et al. (2019b) and North et al. (2020) that had a single-star Hipparcos solution and that were not confirmed triples. This mass distribution is compatible with field WD mass distributions and with those published before for Ba stars. The distribution extends to high WD masses, higher than expected by theoretical models of the s-process of nucleosynthesis that have focused on reproducing the abundance ratios measured on Ba star atmospheres. This work brings new observational constraints for these models and an interesting challenge to our understanding of the formation of Ba stars.
In order to look at Ba stars with new eyes, we plan future direct imaging observations of six of the longest-period systems with SPHERE. On the one hand, this data will provide us with a measurement of the instantaneous angular separation between the components of the system, partially breaking the total mass -semimajor axis correlation and helping us get more accurate masses. On the other hand, we will be able to estimate the bolometric luminosity of the WD, which combined with its mass, can be compared to WD cooling models to estimate the age of these systems.