© ESO 2020

1. Introduction

Hot subdwarf-B (sdB) type stars are core-helium-burning stars that have lost the majority of their hydrogen-rich envelope (M_env ≈ 0.01 M_⊙). Most of these systems have masses very close to the core-helium-flash mass of 0.47 M_⊙ (Heber 2009, 2016). Although many different formation channels have been proposed, it is currently accepted that all sdB stars form solely due to binary interactions (Heber 2016).

Hot subdwarf binaries are excellent probes of binary evolution. Since the masses of the sdB stars in these binaries are typically close to the core-helium-flash mass, many of them must have formed near the end of the red giant phase, which in turn constrains their progenitor properties. Furthermore, the long-period sdBs with main-sequence (MS) companions, which are the subject of this study, are believed to have undergone only one stable mass transfer (MT) phase, and thus may be used to constrain the properties of mass transfer. Finally, these binaries are typically double-lined spectroscopic sources, which allows one to study both the sdBs and their cool companions.

There are three main formation channels for sdB stars. In the common envelope (CE) ejection channel, the sdB star forms as the core of an evolved red giant which has lost most of its hydrogen-rich envelope in a common envelope episode (Paczynski et al. 1976; Han et al. 2002). The CE ejection channel produces sdB binaries with short orbital periods (i.e. hours to a few tens of days) and white dwarf (WD) or MS companions. Conversely, in the stable Roche-lobe overflow (RLOF) channel, a red giant loses its envelope through stable mass transfer (Han et al. 2000, 2002). The stable RLOF channel produces wide sdB binaries with orbital periods of the order of years, up to about 1600 days. Finally, in the WD merger channel, the sdB star forms as a result of a merger of two white dwarfs (Webbink 1984). The WD merger channel creates single sdB stars with masses that can potentially be higher than the canonical value. The outcomes of the different sdB formation channels were studied through binary population synthesis by Han et al. (2003) and, more recently, by Chen et al. (2013). In this study, we focus on long-period sdB binaries with MS companions which formed through the stable RLOF channel and ignite degenerately through a helium flash.

Composite long-period sdB binaries have been predicted by Lamontagne et al. (2000), Green et al. (2001), Aznar Cuadrado & Jeffery (2001) and Reed & Stiening (2004), but the first solved long-period systems based on spectroscopic monitoring programmes were published only a few years ago (Østensen & Van Winckel 2011; Deca et al. 2012). Recently, the orbital properties of 23 wide sdB binaries have been published (see, for example, Vos et al. 2019, for an overview).

These observations have already led to improvements in the binary interactions models for stable RLOF from red giants. First binary population studies of long-period sdB binaries predicted orbital periods of the order of a few hundred days (Han et al. 2002). The periods of the first observed systems with solved orbits, however, were significantly longer. These longer periods were explained with a more detailed treatment of the angular momentum loss and the inclusion of atmospheric RLOF (Chen et al. 2013). A second discovery was that almost all long-period sdBs had eccentric orbits, which contradicted the tidal circularisation models that predict that all sdB progenitor binaries should circularise long before they start transferring mass. Using the stellar evolution code Modules for Experiments in Stellar Astrophysics (MESA, Paxton et al. 2011, 2013, 2015, 2018, 2019), Vos et al. (2015) proposed two possible mechanisms that could explain the observed eccentricities: phase-dependent mass loss and the interaction of a circumbinary disc with the binary.

The most recent discovery was the strong relation between the mass ratio and the orbital period in these wide sdB binaries. This relation was attributed to the stability of RLOF (Vos et al. 2019), but up to now, this has not been explained by binary population synthesis models. In this study, we show that the observed relation between the mass ratio and the orbital period in long-period sdB binaries can be fully explained by the chemical evolution of our Galaxy.

Galactic evolution has a strong effect on the period distribution of long-period sdBs. The Galactic metallicity, [Fe/H], evolves with time, and at the same time metallicity has a strong effect on the sizes of the red giant progenitors of sdB stars. Stars more massive than 1.5 M_⊙ formed relatively recently and, therefore, produce sdB binaries at solar metallicity. Stars of about 0.9 M_⊙ formed early in the Galactic history and, therefore, produce sdB binaries at sub-solar metallicities of about −0.4 (Robin et al. 2003). Compared to solar metallicities, the red giant radii of such stars are about 30% smaller (Choi et al. 2016), and the resulting sdB periods are, correspondingly, shorter. Therefore, to model the P  −  q relation of the whole population of sdB binaries, one needs to model how the Galactic metallicity and star formation rate evolved with time.

In this paper, we perform a population synthesis of composite sdB binaries by using the binary stellar evolution code MESA, by adopting a conventional model of red-giant mass transfer and by accounting for the Galactic chemical evolution by using the Besançon model (Robin et al. 2003, 2014; Czekaj et al. 2014). Compared to typical population studies which ignore the Galactic evolution, we have been able to reproduce the observed P  −  q distribution with minimal assumptions and without explicitly fine-tuning any parameters.

2. Observations

The current sample of long-period hot subdwarf binaries contains 25 systems with solved orbital periods, of which 23 systems have a mass ratio. This sample contains six systems that were part of the long-term monitoring programme with the Mercator telescope (Vos et al. 2012, 2013, 2017). Two systems were solved by Barlow et al. (2013) based on observations with the Hobby-Eberly telescope. For the original long-period sdB binary solved by Deca et al. (2012), we used updated orbital parameters from Deca et al. (2018), and 11 binaries are part of the UVES survey (Vos et al. 2019). Furthermore, the sample contains one system that has an orbital period determined from lightcurves (Otani et al. 2018), but no mass ratio can be determined this way. Recently, also two long-period binaries with sdO components were analysed (Molina et al., in prep.). The two sdO stars have masses consistent with the canonical sdB mass, and are likely in a later evolutionary stage: He-shell burning. As the orbital parameters are not affected by the evolution from He-core burning to He-shell burning, these systems are also included in this sample. Most of these systems were selected based on GALEX colours, and the selection criteria of the UVES and Mercator sample are discussed in Vos et al. (2018a).

Before comparing the observations with models, it is important to check if the selection criteria of the sample and detection limits of the telescopes and instruments have any impact on the observed correlations. These effects have been studied in detail by Vos et al. (2019), where they found that neither the target selection nor the sensitivity of the telescope could be responsible for the observed correlations in orbital period with mass ratio. The sensitivities of the instruments, UVES and HERMES, are also sufficiently high to detect systems with orbital period significantly longer than those observed in this sample. The sample that we analyse here can, therefore, be considered to have no bias for the orbital features that we aim to explain. At the same time, as the observational sample is magnitude-limited to the bright systems in the solar neighbourhood, we expect it to contain significantly more systems from the Galactic thin disc, rather than the sparse and less massive thick disc or the halo. The sample contains no systems from the distantly-located bulge as that region was explicitly avoided.

The wide sdB binary sample shows two main correlations. The first one is the correlation between the orbital period and eccentricity, and the second one is the correlation between the orbital period and the mass ratio. The observed eccentricity of these systems is unexpected as they should have tidally circularised before the onset of mass transfer. This unexplained eccentricity is also observed in other evolved binaries, such as, for example, barium stars, post-AGB stars and He-WD binaries (see Vos et al. 2017, for an overview). Vos et al. (2015) found that phase-dependent mass-loss and circumbinary discs can explain the observed eccentricity in wide sdB binaries. However, while their models allowed for the observed parameters, they could not predict the observed correlation with the orbital period.

In this work, we focus on the relation between the orbital period and the mass ratio. In the correlation, longer orbital periods (P) correspond to smaller mass ratios q = M_sdB/M_comp. Assuming that the sdB binaries formed through degenerate ignition from low-mass progenitors (M ≲ 2.0 M_⊙) and that the sdB mass is close to canonical (0.47 M_⊙), this P  −  q relation is equivalent to saying that longer orbital periods correspond to larger companion masses. The plots showing the observed orbital periods and mass ratios are presented in Sect. 6.

An interesting notion in the observations is that there is a secondary group of systems that follow the same P  −  q and P  −  e correlations as the main group but at shorter orbital periods. In observable features, there is no difference known between hot subdwarfs in the main group and the secondary group. There are only a few systems found in the secondary group, which makes a statistical comparison difficult. Up to now, the reason for the difference in both groups is unknown, although we discuss the possible reasons in Sect. 7.7.

3. Binary evolution model

Long-period sdB binaries result from stable mass transfer from evolved red giants onto main-sequence stars. To result in an sdB object, the naked helium core of the red giant has to ignite helium and have less than 0.01 M_⊙ of hydrogen remaining in its atmosphere (Heber 2016).

We adopt a standard observationally-motivated formation scenario for long-period sdB binaries and schematically show it in Fig. 1. Additionally, in Fig. 2, we show a typical evolutionary MESA track leading to a long-period sdB binary, starting from a 1.2 M_⊙ donor star. We consider binaries with primary masses between 0.7 M_⊙ and 2.0 M_⊙ and initial periods between 100 and 800 days. These are primaries which may have left the main sequence by the present time, ignited their helium core degenerately (Hurley et al. 2000) and initiated mass transfer near the end of the red-giant branch (Han et al. 2002; Chen et al. 2013). Initially eccentric, these binaries circularise before the mass transfer starts (Vos et al. 2015; Beck & Vos, in prep.). As the red giant expands, it fills its Roche lobe and initiates the mass transfer. For a broad range of initial binary mass ratios, q_i = M_primary/M_comp ≲ 2.0, the mass transfer proceeds stably at mass transfer rates below 10⁻³ M_⊙ yr⁻¹, resulting in long-period binaries (Tauris et al. 2000; Pavlovskii & Ivanova 2015). The exact stability criterion for mass transfer also depends on the evolutionary stage of the red giant and the mechanism of mass-loss in the binary.

Fig. 1. Formation channel for the main P − q branch of wide-period composite sdB binaries. Initially, a 0.7 M⊙–2.0 M⊙ main-sequence primary (MS) evolves to become a developed red giant (RG), which initiates mass transfer through Roche lobe overflow. The mass transfer rate gradually grows and exceeds 10−5 M⊙ yr−1, at which point the accretor cannot accept the incoming mass and ejects it in an outflow (shown schematically). Even though the red giant loses all its envelope, its degenerate core ignites helium and turns into a luminous sdB star. The companion main-sequence star accretes only a small amount of mass of ≲0.05 M⊙, but, nevertheless, becomes polluted (MS*). The radius of the RG is strongly sensitive to the metallicity of the primary, which in turn correlates with its age and mass due to Galactic evolution. This way, the Galactic age-metallicity correlation is an important factor determining the final periods of long-period sdB binaries.

Fig. 2. Example MESA track for a binary with a 1.2 M⊙ donor star, a 0.85 M⊙ companion, an initial period Pi = 350 days and a metallicity [Fe/H] = −0.15, which produces a hot composite sdB binary. Top panel: track for the donor star in the HR diagram. The mass loss phase is marked in red; the thin red line shows the phase when Ṁ > 10−9 M⊙ yr−1, thick red line corresponds to Ṁ > 10−6 M⊙ yr−1 and the moment of maximum mass loss is marked with a horizontal red marker. The first He flash is marked with a blue diamond and the core He burning phase is shown in a thick blue line. The phase between the He flash and the core He burning is very short compared to the core He burning phase. Bottom panel: orbital period and mass loss rate during the mass loss phase. The same markings are used as in the HR diagram: the thin black line shows the phase with Ṁ > 10−9 M⊙ yr−1, thick black line corresponds to Ṁ > 10−6 M⊙ yr−1 and maximum mass loss is marked with a vertical black marker.

Observations of companions of long-period sdB binaries indicate that they accrete only a small, up to a few times 0.01 M_⊙, amount of material, that is to say they form through non-conservative mass transfer (Vos et al. 2018a). Therefore, we adopt the simplest non-conservative mass transfer scenario consistent with these observations. As the red giant donor expands, progressively larger parts of its envelope end up outside the Roche lobe. Gravitationally unbound from the donor, these layers flow through the L₁ point and form an accretion disc. Initially, mass transfer rates are low, and we assume that most of the transferred mass is conservatively accreted on the companion. As the mass transfer rates reach Ṁ ≳ 10⁻⁵ − 10⁻⁶ M_⊙ yr⁻¹, the companion accretes a few 0.01 M_⊙ and gets critically spun up by accretion (Popham & Narayan 1991). Subsequently, the companion cannot accrete any further material as the material cannot efficiently remove the excess angular momentum in order to land on the companion (Popham & Narayan 1991; Paczynski 1991; Deschamps et al. 2013). Additionally, at similar mass transfer rates, the gravitational energy released due to accretion cannot be efficiently radiated by the accretor, making it swell, which is also expected to lead to mass loss (Kippenhahn & Meyer-Hofmeister 1977; Pols & Marinus 1994; Toonen et al. 2012). As we discuss in Sect. 5, both effects occur at similar mass transfer rates, but our MESA simulations suggest that typically over-spinning takes place first.

We assume that mass gets lost through a non-collimated jet-like outflow with the angular momentum of the accretor, see for example Tauris & van den Heuvel (2006) and Shiber & Soker (2018):

$$\begin{array}{c}\hfill {\displaystyle \frac{{\dot{J}}_z}{J_z}}={\displaystyle \frac{{(M_{\mathrm{RG}}/M_{\mathrm{comp}})}^2}{1+(M_{\mathrm{RG}}/M_{\mathrm{comp}})}}{\displaystyle \frac{{\dot{M}}_{\mathrm{RG}}}{M_{\mathrm{RG}}}}.\end{array}$$(1)

As follows from Chen et al. (2013) and as we discuss in detail in Sect. 7.7, the exact choice of the angular momentum loss prescription should only have a small effect on the final period-mass ratio correlation of sdB binaries. Depending on the binary parameters, mass transfer rates may reach values of about 10⁻² M_⊙ yr⁻¹. In this case, we assume that binaries enter a common envelope stage, for example Pavlovskii & Ivanova (2015), and produce short-period systems.

The final periods of wide sdB binaries are sensitive to their metallicity, for example Chen et al. (2013). At the same time, Galactic metallicity has evolved with time. The oldest stars generally have metallicity of about 0.4 dex lower than the presently-formed stars (Edvardsson et al. 1993). We illustrate the importance of metallicity through a simple analytic argument here, and present a detailed demonstration based on MESA runs further in Sect. 6.4. Let us consider two binaries: Binary 1 has primary mass of 1.0 M_⊙ and metallicity [Fe/H] = −0.4, Binary 2 has the same primary mass of 1.0 M_⊙, but its metallicity is solar ([Fe/H] = 0). Binary 1 represents an average system (for its primary mass) which formed about 8.4 Gyr ago in the Galaxy and reached the tip of the red-giant branch (RGB) today. The primary then has a radius of about 135 R_⊙ at the tip of the RGB (Choi et al. 2016). Binary 2 represents the same system as Binary 1, assuming it was born at solar metallicity, a much less typical metallicity for the Galaxy at the time the binary was born. Due to its overestimated metallicity, Binary 2 would have a radius of 165 R_⊙ at the tip of the RGB.

We can now estimate the effect of metallicity on the final periods of these binaries. For both Binary 1 and Binary 2 we assume a 0.7 M_⊙ companion which, according to our MESA simulations described later, could lead to the production of an sdB from the primaries in both cases. Due to the difference in red giant radii, the initial periods of these binaries differ by 35 per cent. We can assume that the ensuing mass transfer is fully non-conservative (Vos et al. 2018a) and that the system loses mass with the angular momentum of the accretor. The final periods after the mass transfer phase are then given by Soberman et al. (1997):

$$\begin{array}{c}\hfill {\displaystyle \frac{P_{\mathrm{sdB}}}{P_{\mathrm{init}}}}={\left({\displaystyle \frac{M_{\mathrm{RG}}}{M_{\mathrm{sdB}}}}\right)}^3·{\left({\displaystyle \frac{M_{\mathrm{RG}}+M_{\mathrm{comp}}}{M_{\mathrm{sdB}}+M_{\mathrm{comp}}}}\right)}^2·e^{3({\displaystyle \frac{M_{\mathrm{sdB}}}{M_{\mathrm{comp}}}}-{\displaystyle \frac{M_{\mathrm{RG}}}{M_{\mathrm{comp}}}})}\end{array}$$(2)

Assuming a canonical sdB mass of 0.47 M_⊙ (Heber 2009, 2016), both binaries reach similar final mass ratios of q of about 0.65. Their periods, however, are 1100 d and 1500 d, correspondingly. While the binary with the correctly initialised metallicity ends up on the main branch of the P  −  q relation, the period of the binary with solar metallicity is inconsistent with the main branch of the P  −  q relation.

The argument above, consistent with the detailed MESA analysis in Sect. 6, shows that the choice of initial metallicities can affect the final periods by several hundred days. In contrast, since mass transfer is non-conservative in our model, the choice of metallicity should not have any significant direct effect on the final mass ratios q. It should also be noted that Eq. (2) applies to real systems only if we know whether the initial binary does indeed produce an sdB-MS binary. In particular, as shown by Chen et al. (2013), the final periods of sdBs depend on core mass and the radius of the red giant donor at the end, rather than at the beginning of mass transfer. However, since the radius of red giants even after mass transfer is sensitive to metallicity, the same argument still holds. Generally, in order to produce an sdB, the primary must ignite the core and have lost most of its hydrogen envelope. These conditions may be fully captured only through detailed binary simulations we describe further.

4. Galaxy evolution model

An accurate evolutionary model for the Galactic metallicity is needed to generate the present-day population of long-period sdB binaries. The lowest-mass stars reaching the red giant branch at present formed 10 Gyr ago at low metallicities of about −0.4 (Edvardsson et al. 1993). The more massive stars of 2.0 M_⊙ formed in the last few Gyr at solar-like metallicities. Since the metallicity has a strong effect on the radius of red giants, as we show in Sect. 3, the history of metallicity also has a strong effect on the shape of the P  −  q relation for long-period sdB binaries. In contrast, typical binary population synthesis studies, of sdB binaries and other types of binaries in the Galaxy, have generated populations either by assuming a constant or a uniformly distributed metallicity (see e.g. Izzard et al. 2009; Hamers et al. 2013; Wijnen et al. 2015; Stanway & Eldridge 2018). Although see also for example Lamberts et al. (2019), Boco et al. (2019) and Olejak et al. (2020), where the metallicity evolution is accounted for.

We model the evolution of the Galactic metallicity by using the Besanćon model of the Galaxy (Robin et al. 2003). The Besanćon model describes the evolution of star formation and metallicity in the Galaxy and is calibrated by large photometric surveys of the field.

The model represents the Galaxy through 10 population bins: 7 bins for the thin disc and one bin each for the bulge, thick disc and the halo, all of which had a relatively short star-formation history. Each population bin is associated with its stellar mass fraction and mean and 1-σ range of metallicities.

We summarise the parameters of the Galactic bins in our study in Table 1. In the table, the mass fractions are expressed in terms of the Galactic stellar mass, for which we use the value of M_MW = 6.43 ⋅ 10¹⁰ M_⊙ from McMillan (2011), consistent with the recent study by Licquia & Newman (2015). Thin disc parameters, and the stellar mass fraction of the thick disc are based on the original study by Robin et al. (2003). We obtain mass fractions for the thin disc bins by integrating the Galactic profiles from Robin et al. (2003). For the mass fraction of the Galactic bulge and the age of the thick disc, we use more recent data from Robin et al. (2012, 2014). Stars are assumed to form at a constant rate within each stellar bin, and the initial mass function (IMF) is based on Kroupa & Haywood v6 IMF from Czekaj et al. (2014).

5. Simulations

5.1. Modifications to the MESA code

In this work, we have used the stellar and binary evolution code MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019), version r10390. In order to follow through the evolution described in Sect. 3, we have placed a limit on the mass accretion of the secondary based on its rotational velocity. Our MESA models have shown that the companion reaches critical rotational velocity long before it starts swelling due to the thermal energy released from accretion. By this point, the donor loses at most ≲0.03 M_⊙, which is close to the estimates of the mass accreted by the companions of the observed sdBs, for example Vos et al. (2018a). Therefore, the companion is allowed to accrete all mass lost from the donor star until it reaches its critical rotational velocity. At that point, the remainder of the material is considered to be lost from the system with the angular momentum of the accretor. In terms of stability, a system is considered stable as long as the donor mass-loss rate does not exceed 10⁻² M_⊙ yr⁻¹, which is a common criterion in binary evolution studies, for example Pavlovskii & Ivanova (2015). This is the case of all the long-period hot subdwarf systems produced in our model sets, and we additionally kept track of the peak mass transfer rates reached in our runs. To be clear, we have not used the fixed mass-loss fractions provided in MESA, but dynamically calculate the amount of mass accreted by the companion and lost from the system at each timestep. To extract parameters of interest from the MESA runs we have used the Neural Network assisted Population Synthesis code¹ (NNaPS, Vos et al. 2020).

5.2. Extracting observational properties from the MESA runs

The outcomes of MESA simulations must be identified with sdB binaries. The term hot subdwarf star was originally a purely spectroscopic classification. More recently, it is mostly used to describe only the core He-burning stars that match the spectroscopic parameters. Therefore, in this study, we identify the MESA outcomes as sdB binaries based on their position on the HR diagram. Specifically, we require the star to be core He-burning, and base the selection criterion on the average effective temperature during the core He-burning phase. If the star reaches an effective temperature between 15 000 and 40 000 K during its core He-burning phase, it is considered a hot subdwarf. We also only considered hot subdwarfs that ignite He under degenerate conditions.

Furthermore, we make a division between sdA-type stars with T_eff <  20 000 K and sdB-type stars with higher effective temperatures. For the hot subdwarfs that result from degenerate He-core ignition, the average effective temperature and surface gravity during the He-core burning phase are strongly related. This relation may be seen in Fig. 3, which shows the He-core burning systems resulting from our runs. This figure also shows that the sdBs seem to group together at roughly 28 000 K and at 34 000 K. This is a known effect that can be attributed to the treatment of convection during the He-flashes in MESA. The group with the hottest sdBs with temperatures above 34 000 K ignite late on the WD cooling track and convectively mix H and He in the outer layers, igniting most of its hydrogen in a H-flash following shortly after the He-flash, and this way leaving a hot core with almost no H remaining. The systems with lower temperatures experience He-flash early on without igniting hydrogen, this way leaving cooler H atmospheres. See for example Xiong et al. (2017) and Wu et al. (2018) for a discussion. The remaining envelope mass during the core He burning phase for sdB stars ranges from 0.0005 to 0.004 M_⊙ while sdA stars have a slightly larger envelope masses ranging from 0.003 to 0.018 M_⊙.

Fig. 3. Average effective temperature and surface gravity during the core He burning phase for the hot subdwarfs in the MESA sample. There is a strong relation between log g and Teff for hot subdwarfs formed from progenitors with degenerate He cores. Hot subdwarfs with Teff <  20 000 K are classified as sdA stars while the hotter subdwarfs are classified as sdB stars.

Having identified hot subdwarf stars in the MESA models, we select those binaries which would have been identified spectroscopically as hot composite sdB binaries, thus modelling the observed sample. In other words, we apply the same selection criteria to the model sample as the ones used to select the observational composite hot subdwarf sample. Firstly, the observed sample consists only of systems from the thin and thick disc, and potentially from the halo. The systems in the bulge would be too faint for the current observing programmes. The second criterion is that the system should be recognizable as a composite binary. Namely, it should be possible to identify both the sdB and the cool companion in the spectrum. In the observational sample, systems are selected based on visual and IR photometry. However, if it is not possible to recognise both components in the optical spectrum, the systems are rejected from the observational sample.

The Gaia catalogue (Gaia Collaboration 2016, 2017) provides the most comprehensive database or stellar positions, magnitudes and proper motions potentially suitable for selecting large samples of candidate long-period composite sdB binaries. It has proven very useful to build catalogues (e.g. Geier et al. 2019; Pelisoli & Vos 2019) and is therefore ideal as a base for an observability criterion. We have calculated synthetic colours and spectra for different sdB+MS combinations and found that a criterion based on the sdB flux contribution in the Gaia-G band, as described below, serves as a good indicator for composite binaries. For systems where the sdB flux contribution is between 20 and 90%, both components will most likely be recognizable in the optical spectrum. These tests were done with solar metallicity spectra for the companion, and low rotational velocities (v sin i = 10 km ⋅ s⁻¹). Changes in metallicity and rotational velocity can change these results, but we are confident that for this purpose the Gaia-G band selection criterion is sufficient. An example of two sdB+MS systems at both ends of the visibility selection criterion is given in Fig. 4. The systems that fulfil this flux-based selection criterion will be referred to as composite binaries, while the systems that do not fulfil this criterion will be called single-lined binaries (regardless of which component is visible).

Fig. 4. Spectral energy distribution for an sdB+MS binary with 90% light contribution of the sdB in the Gaia-G band in the top panel, and with 20% contribution in the bottom panel. The binary SED is plotted in black while the contributions of the sdB and MS component are shown respectively in blue and green. In black-dotted line, the transmission curves for the Gaia-G, 2MASS-J and 2MASS-Ks filter are shown from left to right.

5.3. Main set of runs

We summarise the main and the additional sets of runs performed in this study in Table 2. The main model set of this article is the model set of all possible and likely sdB progenitors accounting for the Galactic evolution of metallicity and star formation rates; Model set 1. This set consists of 2060 MESA simulations based on a synthesised population described as follows.

First, we initialised a larger population of systems representing all the binaries in the Galaxy (in practice, about 10⁴ times larger then the number of MESA simulations). We used the Kroupa-Haywood v6 initial mass function (Kroupa 2008; Haywood et al. 1997) as implemented in the Besançon Galactic model by Czekaj et al. (2014). Therein, the IMF, dN/dm ∼ m^−α, is set by slope α equal to 1.3 for m between 0.09 M_⊙ and 0.5 M_⊙, to 1.8 for m between 0.5 M_⊙ and 1.53 M_⊙ and to 2.3 for m between 1.53 M_⊙ and 120 M_⊙. After drawing the mass of the primary, we drew its Galactic bin by using bin mass fractions in Table 1 as weights. Then we initialised the age by drawing uniformly in the age range of the Galactic bin and chose the metallicity by drawing a Gaussian variable with the parameters from Table 1. We drew the companion mass assuming a uniform distribution in M_comp/M_primary (Raghavan et al. 2010) and the orbital period P assuming a flat distribution in log P in the range between 1 and 10⁴ days (Abt 1983). We assumed circular binaries since red giants efficiently circularise before overflowing the Roche lobe (Vos et al. 2015).

For the actual MESA simulations, we selected from the large population described above only those binaries which may, in principle, result in long-period composite sdB binaries forming at present epoch. In particular, we considered only binaries with primary masses above 0.7 M_⊙, since lower-mass stars remain on the main sequence, and below 2.0 M_⊙, since higher-mass stars ignite the helium core non-degenerately (Hurley et al. 2000) and are thus expected to lead to a different population of binaries. Furthermore, we selected the binaries with M_primary/M_comp between 1 and 3 since higher mass ratios lead to unstable mass transfer on the giant branch, for example Tauris et al. (2000) and Pavlovskii & Ivanova (2015). We selected the period to be in the range between 0.1 and 1.4 of the period corresponding to Roche-lobe overflow at the tip of the red giant branch as determined from single evolutionary tracks from MIST database (Choi et al. 2016). The lower limit of 0.1 corresponds to Roche lobe overflow on the sub-giant branch, at which point the helium cores do not have enough mass to ignite. The higher limit of 1.4 was set to ensure that the binary MESA runs lead to Roche lobe overflow on the red giant branch in agreement with the single track-based estimates from the MIST database. Finally, we selected the binaries which become sdB stars at the present time. Specifically, we selected the binaries in which the primary reaches the tip of the red giant branch (based on single evolutionary tracks from MIST) not earlier than 300 Myr before the present day and not later than 700 Myr after the present day. If the tip was reached earlier than 300 Myr before the current time then, even allowing for possible 100 Myr difference between MIST and actual MESA tracks, the sdB star would have already reached the end of its lifetime and not be observed as an sdB star at present. If the tip was reached more than 700 Myr after present, then even allowing the possible difference of 100 Myr between MIST and MESA, the primary would still be on the main sequence at the present time and not able to produce an sdB. This way, through the selection process, the initial sample of 3.0  ⋅  10⁷ general Galactic binaries was reduced to a sample of 2060 systems which may potentially produce sdB binaries today. By examining the results of our MESA simulations, we ensured that all the produced sdB binaries originated from the systems fully enclosed within the parameter space we selected here.

A large fraction of the systems in the main set of runs (Set 1) are members of the thin disc as it is the most massive component of the Galaxy, whereas the thick disc and the halo, correspondingly, contribute little or do not contribute at all to the systems in the figure. To investigate the P  −  q relation caused by the different components of the Galaxy in more detail (specifically, the thick disc, halo and bulge), we created separate model sets focused on those populations. These additional sets are numbered 1.1, 1.2 and 1.3 for, respectively, the thick disc, bulge and the halo and are summarised in Table 2.

To estimate the effect of using the Galactic model compared to more simple metallicity prescriptions, we simulated 1200 binaries at fixed metallicities of 0, −1, and uniformly distributed metallicities between −1 and 0.25, also focusing on the sdB systems forming in the Galaxy today and keeping all the other parameters the same. These runs correspond to runs 2, 3 and 4 in Table 2.

To demonstrate the importance of using the MESA code for producing sdB binaries as opposed to using synthetic codes which do not evolve stellar structure, we have re-simulated the main set of runs with the BSE code (Hurley et al. 2002), which corresponds to Set 5 in Table 2. To reproduce the evolution described in Sect. 3, we also had to modify the BSE code. The default implementation of red giant mass transfer in BSE significantly underestimates the critical initial mass ratio leading to stable mass transfer, for example Tauris et al. (2000) and Pavlovskii & Ivanova (2015). Therefore, we replaced the stability criterion based on simplified polytropic models by the stability condition q_i = M_primary/M_comp ≤ 2.0, to represent the typical critical mass ratios we observe in MESA runs, as we present further in Sect. 6. Furthermore, we have modified the code to produce mass loss with the angular momentum of the accretor. Finally, we have removed the limiter on the red giant mass transfer rate which is set to be the thermal-timescale mass transfer rate of about 10⁻⁶ M_⊙ yr⁻¹, since in MESA runs red giants reach mass transfer rates well above that value.

Finally, to explore in detail the effects of different initial binary parameters on the final periods and mass ratios, we have simulated additional 121 binaries using the MESA code, which corresponds to Set 6 and which we present in detail in Sect. 6.4.

To obtain the Galactic formation rates from our study, we used the total stellar Galactic mass of 6.43 ⋅ 10¹⁰ M_⊙ (McMillan 2011), the overall binary fraction of 0.45 and the close binary fraction (sub-fraction of the binaries which are in period range of 1–10⁴ d) of 0.25, suitable for the mass range and the typical metallicities considered here (Moe et al. 2019). With these parameters set, we could identify how many stars in the Galaxy produce one sdB-progenitor binary and this way, by considering the sdB formation rate in our simulations, to re-normalise it to the sdB formation rate in our Galaxy. We calculated the present-day number counts of different sub-groups of sdB binaries by using their mean harmonic lifetimes represented by their core-helium burning phase. We ignored the effects of anti-correlation of the close binary fraction and metallicity in this study (Moe et al. 2019) to identify the effect of purely-Galactic metallicity evolution.

6. Results

Our main set of 2060 runs (model set 1) has produced 149 sdB and 40 sdA binaries, of which 89 sdB and 28 sdA binaries would have been identified as long-period composite subdwarf binaries in actual observations following the criteria described in Sect. 5.2.

6.1. P  −  q relation

Our models can explain the main branch of the P  −  q relation in long-period sdB binaries by using a simple Galactic chemical evolution model and without any free parameters, which is the most important result of this study. In Fig. 5, we show the P  −  q relation of the models obtained for several different metallicity distributions from Table 2. In these plots, the blue circles show the produced systems which would have been observed as hot composite sdB binaries, that is fulfil the visibility criteria described in Sect. 5.2. The real observations are shown in grey rectangles.

Fig. 5. Period-mass ratio relation of the sdB binaries produced by the different MESA samples. The grey squares show the real observed systems, while the blue circles show the MESA models which would have been observationally identified as composite long-period sdB binaries, i.e. fulfil the visibility criteria in Sect. 5.2. The observed systems, located within the nearest 1 kpc, belong to the solar neighbourhood. Most of the simulated sdB binaries in the Galactic population belong to the thin disc, which is the dominant Galactic component (see Table 2).

The model set with fixed metallicity and constant star-forming rate (model sets 2 and 3) can only cover a small part of the observed P  −  q relation. This result is expected as the RGB radius is strongly dependent on the metallicity, and as the system needs to interact when the RGB core is close to the He ignition mass, the period range at which the interaction can take place and form an sdB is small. The sample with a sub-solar metallicity of [Fe/H] = −1.0 has orbital periods around 500 days while the solar metallicity sample reaches orbital periods around 1200 days. The final mass ratio range is similar for both samples and is between q = M_sdB/M_comp = 0.3 and 0.8. This result is consistent with the idea that the main effect of metallicity is on the final periods rather than mass ratios. We demonstrate this effect in more detail in Sect. 6.4. The sample with a sub-solar metallicity of −1.0 covers most of the systems in the secondary branch of the P  −  q relation.

The model set with a uniform metallicity distribution (model set 4) can cover all observed orbital periods, including those of the secondary branch. However, this sample can not, in any meaningful way, explain the observed P  −  q correlation of the main group. As the metallicity can freely vary between [Fe/H] = −1.0 and 0.25, sdB binaries with orbital periods between ∼500 and ∼1400 days can be created, as expected from model sets 2 and 3.

The last model set in Fig. 5, which follows the Galactic model (model set 1), covers the P  −  q correlation in the main group very well, predicting systems with high mass ratios at short orbital periods and low mass ratios at long orbital periods. The relation causing this pattern is shown in Fig. 6. This figure shows the sdB progenitor mass as a function of the initial metallicity of the system. The systems with a lower-mass sdB progenitor formed when the Galaxy was younger and consequently have a lower metallicity. These systems need to have a shorter orbital period to interact when the sdB progenitor has a sufficiently massive core and will produce shorter orbital period systems. Since lower-mass progenitors also have lower-mass companions, the resulting sdB binaries have high final mass ratios q = M_sdB/M_comp and end up on the top of the P  −  q relation. Higher-mass massive progenitors, on the opposite, have higher metallicities, hence higher final periods, and higher-mass companions, and hence lower final mass ratios. As a result, higher-mass progenitors lead to sdB binaries at the bottom of the P  −  q relation. For a given sdB progenitor mass and orbital period, there is only a small range in mass ratios and initial periods that will be stable but still have sufficient mass loss and helium core size to result in an sdB star. We further explore the effects of all the initial parameters in detail in Sect. 6.4.

Fig. 6. Initial metallicity as a function of the sdB progenitor mass for all systems in the model sample following the Galactic model (sample 4). Systems that fulfil the criteria to be recognizable as a composite sdB binary are plotted in filled blue circles, while those where only one component is visible are shown in open grey circles.

The different Galactic populations outlined in Table 1 can populate different parts of the P  −  q diagram. This is shown in Fig. 7 for the four main Galaxy components: the thin disc, the thick disc, the halo and the bulge. Since the thin disc is the most massive component of the Galaxy, most modelled systems are located in the thin disc. The Galactic bulge, which is an old high-metallicity population, will produce sdB binaries at long orbital periods, to the right of the currently observed P  −  q relation, as well as sdB binaries which overlap with the currently observed P  −  q relation. Systems from the thick disc have a lower metallicity and can be found at slightly shorter orbital periods. The halo contributes so little to the Galactic mass (Table 1) that there are no sdB binaries, observable as composite systems, formed from the halo population in the main sample. As follows from the halo set of runs (Set 1.3), members of the halo population have short periods below 500 d and are located to the left of the thick disc population. While no Galactic population in our model covers the systems in the secondary branch in the P  −  q relation, Figure 5, we discuss the possible origins for these objects in Sect. 7.7.

Fig. 7. Schematic plot showing where wide sdBs formed from the different Galactic sub-populations are located in the period – mass ratio diagram. Thin-disc objects are shown in blue, thick-disc objects are shown in orange, bulge objects are shown in green and halo objects in purple. The observed systems are plotted in grey squares.

We have also verified if synthetic stellar population codes, such as BSE, can reproduce the P  −  q relation of sdB binaries. By running the same Galactic population as in model set 1 through the BSE code (model set 5), we found that even when modified to have the same angular momentum-loss prescription and similar stability criteria as MESA, all the sdB-binaries acquired orbital periods below 400 days. This result agrees with the original models by Han et al. (2003). As also discussed in Chen et al. (2013), these binaries fall into a region of the P  −  q plane which is quite far away from the observed systems. Our study shows that the reason why the BSE code does not reproduce the P  −  q relation is not the angular momentum loss prescription. Indeed, in contrast to our MESA results, even with the same angular momentum loss prescription as in MESA, no binaries with initial periods longer than 300 days have produced sdB objects in the BSE sample. In other words, in BSE, sdB binaries were produced by progenitors with incorrectly-short periods. The too short periods of these progenitors explain the too short periods of the resulting sdBs. In contrast, the objects which would have been progenitors in MESA did not produce any sdB binaries in BSE. Therefore, we conclude that the reason why BSE produces too short periods of sdBs is the incorrect criteria by which the BSE code defines the formation of stripped helium stars. The BSE code is based on interpolating over the models of single stars, whereas the formation of sdB-binaries is strongly sensitive to the details of binary mass transfer, for example Chen et al. (2013). Therefore, the conditions, by which the BSE code defines whether a stripped He star has formed, are not consistent with the more detailed MESA simulations.

6.2. M_sdB – Metallicity – P_final correlations

Since the maximum red giant radius strongly depends on metallicity, the final orbital periods of sdB binaries must also correlate with their initial metallicity. We show this relation in Fig. 8. While the periods and metallicities have been known to correlate for hot subdwarfs, (e.g. Chen et al. 2013), and other post-mass-transfer systems, such as low-mass and extremely-low-mass WDs, (e.g. Istrate et al. 2016), this relation has never been depicted for a MESA-based synthesised population which follows the Galactic metallicity and star formation history. Neither has it been compared to actual observations.

Fig. 8. Relation between the orbital period and the metallicity for the sdB binaries (left) and He-WD binaries (right) in our sample. The relation is a result of the donor radius on the RGB being strongly sensitive to its metallicity. The period-metallicity relation also depends on the sdB/He-WD mass, which is shown by the colour scale. The systems that are recognizable as composite sdBs are shown in filled circles while the single-lined ones are shown in grey open circles. For the He-WDs no distinction is made between composite and single-lined systems.

The correlation of the final orbital periods with the initial metallicity may potentially be used to test the model. Moreover, the correlation is stronger for sdB binaries compared to other post-mass transfer systems, which is because, for a given metallicity, the red giant radius is the dominant parameter determining the final period (as we discuss in Sect. 6.4). Indeed, the range of possible sdB masses is very narrow (in all our models, between 0.44 and 0.49 M_⊙) and the range of initial periods and mass ratios leading to sdB binaries is also fairly narrow (at most, 40%). Therefore, the relation between the period and initial metallicity of sdB binaries is very tight and has a very small spread, as shown in the left panel of Fig. 8.

Only 12 wide composite sdB binaries have measurements of their metallicities available in the literature (Vos et al. 2017, 2018a, 2019). These measurements have, on average, large errors as the spectra used to derive the metallicities have low signal-to-noise ratios. In Fig. 9, the observed period-metallicity relation is shown together with the predictions from our models. The observations match very well with the predictions, and this way provide independent verification of our model.

Fig. 9. Observed versus modelled orbital period and metallicity correlation for long-period sdB binaries. The black squares show all the observed wide sdB binaries for which such measurements are available. The results from our MESA model following the Galactic model (model set 4) are shown in blue circles. Only model systems that fulfil the observability criteria are shown.

For comparison, we show the same relation for the He WD outcomes in our Galactic sample in the right panel of Fig. 8. These are systems in which the donor star lost too much mass to ignite He, and which end as He WDs. They range in mass between 0.30 and 0.46 M_⊙. The He WDs show the same relation of longer orbital periods at higher metallicity for systems with the same final donor mass. However, as the mass range for He WDs and their companions is much larger than for sdBs, and as the mass of a He WD is not easy to determine accurately, there is no simple relation that can be tested from the observations of He WDs.

The mass of an sdB star formed from a low-mass progenitor igniting He under degenerate conditions is close to the canonical mass of 0.47 M_⊙ but varies with orbital period and initial metallicity. The MESA models here can be used to derive a fitting equation of the sdB mass as a function of the orbital period and metallicity. An equation that is linear in these two parameters can reproduce the sdB mass with an error less than 0.1%:

$$\begin{array}{c}\hfill M_{\mathrm{sdB}}=0.395+0.005·\frac{P_{\mathrm{orb}}}{100\mathrm{d}}-0.058·[\mathrm{Fe}/\mathrm{H}]\end{array}$$(3)

A similar equation with an error less than 0.2% can be derived for the He-WDs in our sample:

$$\begin{array}{c}\hfill M_{\mathrm{WD}}=0.335+0.012·\frac{P_{\mathrm{orb}}}{100\mathrm{d}}-0.048·[\mathrm{Fe}/\mathrm{H}]\end{array}$$(4)

It is important to notice that Eq. (4) applies to the Galactic binaries in which the WD has formed recently (within the last few 100 Myr) and is only valid for the parameters range covered in our models. This range corresponds to the same range that is shown in Fig. 8.

6.3. Linking to the progenitor properties

Our models offer an interesting possibility of connecting the properties of the observed sdB systems to their progenitors. Such an inference is possible because several present-day observables are strongly correlated to the initial parameters of the binaries. The two strongest correlations are between the final mass ratio and the sdB progenitor mass, and between the final mass ratio and the initial orbital period. These relations are shown in Figs. 10 and 11. These correlations can be applied to the data from ongoing large spectroscopic surveys, such as APOGEE, Gaia-ESO, WEAVE and 4MOST. These surveys, together with Gaia DR3, will provide a large sample of wide main-sequence binaries with periods and metallicities, for example Price-Whelan et al. (2020) and Hayden et al. (2015). These data sets can be used to verify our models further and to connect the sdB populations to their progenitors.

Fig. 10. Correlation between the mass ratios in the observable composite sdB binaries (filled circles) and the initial sdB progenitor masses in our models. The colour shows the initial mass ratios and the empty circles show the sdB binaries, which would not be observationally identified as composite systems. One may see that the observed mass ratios constrain, to a certain extent, the initial progenitor masses of sdB binaries.

Fig. 11. Correlation between the mass ratios in the observable composite sdB binaries (filled circles) and the initial orbital periods of the systems in our models. The colour shows the initial metallicity and the empty circles show the sdB binaries, which would not be observationally identified as composite systems. Given the rather tight correlation, one may constrain the initial period and, in certain cases, the progenitor metallicity from the present-day mass ratios.

In the first figure, it may be seen that systems with a larger final mass ratio have a lower sdB progenitor mass. As may be seen from the figure, the spread in this relation is caused by the initial mass ratio. By examining the spread, we see that, for a given sdB progenitor mass, there is a limited range in companion masses that will lead to an sdB star. If the companion mass is outside of that range, a He-WD or a horizontal-branch (HB) star is formed instead.

The second strong correlation is between the initial orbital period and the final mass ratio. This correlation shows that systems with larger observed mass ratios are produced from systems with larger initial orbital periods. For a given observed mass ratio, the initial orbital period also varies with the metallicity of the system. When both the mass ratio and the metallicity can be determined, the initial orbital period can be derived quite precisely.

Finally, as we show in Fig. 12, our models imply that the companion masses in the observed sdB binaries should correlate with their metallicities. This correlation may, in principle, be used to constrain our model, or, alternatively, the metallicity history of the Galaxy. The advantage of this correlation is that these properties do not require the determination of the orbital parameters, which makes it easier to obtain a large sample in less time.

Fig. 12. Metallicity of the companion stars as a function of the companion mass after at the binary sdB formation time. The box plot shows the average metallicity with an orange line, the 1-sigma range with a box and the 3-sigma range with whiskers. It may be seen that the present-day metallicity and companion masses of sdB binaries are correlated, according to our model.

6.4. Sensitivity to initial parameters

To investigate the sensitivity of the final parameters to the initial parameters, we ran a dedicated set of models by varying each initial parameter separately (Set 6 in Table 2). The results of these runs are shown in Fig. 13. We chose models with three representative initial donor masses of 1.0, 1.2 and 1.8 M_⊙, each of which leads to a long-period hot composite sdB binary. For each of these three masses, we varied the initial mass ratio, metallicity and initial orbital, in equal steps, in the range that would still produce a hot subdwarf star after the mass-loss phase. If the initial mass ratio is too high, too much mass is lost and the donor ends up as a He-WD. If the mass ratio is too low, too little mass is lost, and the donor ends up as a horizontal branch star. Similarly, if the initial period is too short, the red giant core at the end of mass transfer is too small in mass to ignite helium, whereas if the period is too large, too much envelope remains at the end of mass transfer to produce an sdB star. The range in metallicity that can produce a hot subdwarf star is large. However, due to Galactic evolution, only a small range in metallicity is possible for each initial donor mass. In Fig. 13 we, therefore, show only the one-sigma range of possible metallicities for each initial donor mass.

Fig. 13. Effect of the initial parameters: orbital period (Pi), mass ratio (qi), metallicity ([Fe/H]) and sdB progenitor mass on the final orbital periods and mass ratios of sdB binaries. We consider primary masses of 1.0, 1.2 and 1.8 M⊙, indicated in boldface below each set of tracks. Systems that result in hot composite sdB stars are marked with blue circles, composite sdA stars are marked with purple diamonds, while systems resulting in He-WDs or horizontal branch stars are marked with green triangles and orange squares respectively. The locations of the main and second branches of the observed P  −  q relation are indicated with the red dashed ellipses. The initial mass ratios and orbital periods are varied in equal steps in the range that still produces hot composite subdwarf stars for the given initial donor mass, while, for the metallicity, we show the one-sigma range of the Galactic metallicity corresponding to the time when the progenitors of each given mass formed (Table 1). This figure clearly shows that for a given initial donor mass, only a small range in the final mass ratios can be reached by varying the other initial parameters. Varying the metallicity has almost no effect on the final mass ratio, but strongly changes the final orbital period. The correlation between the final mass ratio and the initial donor mass combined with the link between metallicity and initial donor mass due to Galactic evolution then explains the observed P  −  q relation. See Sect. 6.4 for details.

From these models, it is clear that the dominant factor determining the final mass ratio is the initial donor mass. This dependency may be understood as follows. Lighter donors have on average lower-mass companions, whereas heavier donors have on average higher-mass companions. In contrast, the final sdB masses only weakly depend on the primary mass. Similarly, since companions accrete only small amounts of mass, as we discuss in Sect. 3, their final masses remains close to their initial masses. Therefore, lighter donors with lower-mass companions lead to larger final mass ratios q = M_sdB/M_comp, and heavier donors with higher-mass companions lead to lower final mass ratios. This relation between the initial donor mass and the final mass ratios is also shown in Fig. 10.

For a given initial donor mass, the final mass ratio will only change little by varying the initial mass ratio or the initial orbital period. For example, for the system with an initial donor mass of 1.2 M_⊙, the final mass ratio can vary from 0.48 to 0.55 by changing the initial mass ratio or the orbital period. Similarly, changing the metallicity of the system only has a minimal effect on the final mass ratio. It is also important to notice that systems with a lower initial donor mass have a higher final mass ratio and vice versa.

The initial metallicity affects the final orbital period significantly. The one-sigma range in metallicity can change the final orbital period by 200 – 300 days. Since heavier donors have higher metallicities and thus create sdB binaries with longer final orbital periods than the lighter donors, there is a correlation between the initial donor mass and final orbital period. Combining this and the earlier correlation between the primary mass and the final mass ratio then explains the observed period-mass ratio relation.

Long-period sdA binaries behave similarly to sdB binaries. The main difference between both types of systems is that sdA binaries have more envelope mass left and are therefore cooler (see Fig. 3). In Fig. 13, the sdA binaries are located on the right ends of the hot subdwarf series, just before they become horizontal branch stars. Therefore, the relation between orbital period and mass ratio that exists for sdB binaries also exists for sdA binaries. As these systems have more envelope mass left, their P  −  q relation is shifted to longer final orbital periods.

6.5. Mass distribution

The mass distribution of the sdB components in the Galactic sample is shown in Fig. 14. Mass distributions of the sdB stars in composite binaries and single-lined systems are shown, respectively, in blue and grey. Both distributions peak at similar masses, and have a similar spread. The mass distribution of the sdBs in composite systems peaks at 0.465 M_⊙ and has a spread of about 0.01 M_⊙. These results are consistent with the earlier results of Han et al. (2002), and the observed distributions of for example Fontaine et al. (2012) and Schneider et al (2020).

Fig. 14. Mass distribution of the sdB stars in the Galactic sample. The distribution of the composite systems is shown with the blue line, while that of sdBs in single-lined systems is shown in grey. The bins were selected using Knuths rule (Knuth 2006). Both distributions are similar and consistent with the empirically determined masses, for example Fontaine et al. (2012) and Schneider et al (2020).

7. Discussion

We have shown that the main branch of the observed P  −  q relation in long-period composite sdB binaries may be fully explained, without free parameters, through a standard binary interaction prescription implemented in the MESA code and a standard model for the chemical history of the Galaxy. In this section, we show that both the MESA code with the binary interaction model and the model for the Galactic evolution have been essential ingredients in our model. Our method can likely improve the modelling of other post-mass transfer binaries containing hot subdwarf and HB stars, white dwarfs and, possibly, compact objects. Furthermore, by applying the criteria for the observational selection of long-period composite sdB binaries, we have shown that our model predicts several observational correlations. These correlations may be used to probe the progenitor properties of sdB binaries, as well as to test the model. Finally, our study offers an exciting possibility of mutually constraining the models of the Galactic chemical evolution and the observed intrinsic properties of sdB binaries, as well as of other exotic binaries in the field.

7.1. The importance of Galactic evolution

It is well known that metallicity has an important effect on the evolution of stars. For binary systems, the initial metallicity determines, to a large extent, the red giant radius of the stars and thus dictates when the interaction phase starts on the RGB. Therefore, metallicity has an important effect on the orbital period distribution of post-interaction systems. However, metallicity is often ignored in binary population synthesis studies. Typical binary population studies use a fixed (usually solar) or, in the best case, a uniform metallicity distribution (see e.g. Izzard et al. 2009; Hamers et al. 2013; Wijnen et al. 2015; Stanway & Eldridge 2018).

Using a fixed or uniform metallicity distribution has a benefit of simplifying the models and has its merits in some cases. In a realistic setting, however, stellar metallicity is coupled to the stellar age due to Galactic evolution. This coupling is particularly important when dealing with low-mass stars that have long MS lifetimes. If such stars reach the RGB, they must have formed in sub-solar metallicity environments, and are thus sensitive to metallicity changes in the history of the Galaxy. In our models, stars with masses smaller than about 1.4 − 1.6 M_⊙ are most sensitive to the metallicity increase with time, see for example Fig. 6. Similar sensitivity will likely manifest itself in properties of other Galactic populations which are experiencing or have recently experienced mass transfer on the RGB from a low-mass donor. These include short-period sdB binaries (Heber 2016), symbiotic binaries (Munari 2019), extremely low-mass WD binaries (Pelisoli & Vos 2019) and binary pulsars with a He-WD companion (Istrate et al. 2014), among others.

The importance of the cosmological metallicity history has been recently recognised in population studies of binary compact objects, such as black holes or neutron stars, due to its effect on massive stellar evolution, for example Lamberts et al. (2018), Boco et al. (2019), Chruslinska et al. (2019), Neijssel et al. (2019) and Toffano et al. (2019). At the same time, only very few Milky Way studies have considered a realistic metallicity history so far, for example Lamberts et al. (2019) and Olejak et al. (2020), focusing on gravitational wave sources. We expect that using simple, realistic Galactic metallicity models, as adopted in this study, will be common in modelling other Galactic populations in the near future.

7.2. The importance of using MESA for modelling sdB binaries

Modelling sdB binaries requires the use of a binary evolution code such as MESA, rather than traditional binary population synthesis codes. Such codes are based on fits to evolutionary tracks of single stars in isolation and are not designed to model a realistic response of red giants to mass loss. In comparison, binary evolution codes continuously solve the stellar structure equations for both stars, this way accurately evolving red giants through the envelope loss and the subsequent He-flash phases which lead to the formation of sdB stars. Furthermore, population synthesis codes make use of simplified analytic prescriptions when determining outcomes of binary mass transfer, this way limiting, potentially, the correctness of the results. For example, the BSE code determines the stability of red giant mass transfer by using polytropic models of red giants. In its standard implementation, BSE predicts that most of the sdB models calculated in this article would have been unstable and would have entered a common-envelope phase.

In Sect. 6, we showed that even when modified to reproduce similar stability criteria as observed in the MESA code and to follow the same binary interaction model, BSE code still fails to produce sdB binaries at the observed orbital periods. This issue is related to the incorrect initialisation of the He star stage in the BSE code compared to MESA. The reason for the difference may be because the BSE fits for initialising He stars were designed under the old assumption of polytropic models, which lead to incorrect stability criteria for red giant mass transfer. Nevertheless, with the above caveats in mind, traditional BPS codes may certainly benefit significantly from using Galactic metallicity distributions.

We expect that binary evolution codes such as MESA will be increasingly more important in the future. While such codes are more computationally expensive than population synthesis codes (tens of hours versus milliseconds to model a binary on a single modern CPU), population samples of a few thousand objects can be efficiently studied with binary evolution codes. If the initial parameters are chosen carefully, a sample of a few thousand objects can be statistically significant. As the computational power increases over time, binary evolution codes will allow for modelling progressively more rare populations and also provide useful calibrations to synthetic codes.

7.3. Observational tests and predictions

The models in this article match the main P  −  q branch of the observed wide sdB binaries very well. The model sets that use the same binary interaction model but with a simplified metallicity distribution (fixed or uniform [Fe/H]) can match the observed range of mass ratios, and in the case of uniform [Fe/H], also the range of orbital periods. However, none of these models can explain the correlation between the orbital period and mass ratio. It is now clear that this correlation is mainly an effect of the metallicity evolution in the Galaxy.

Our models also predict a very strong correlation between the orbital period, sdB core mass and metallicity. This correlation is due to the RGB radius dependence on metallicity, which is one of its main determining parameters. Of these quantities, the orbital period and metallicity can be observed directly. Even though only a few systems have a published value for the metallicity, and usually with a large error, the observations match the predicted period-metallicity relation very well (see Fig. 9). This is the first time the observed period-metallicity correlation of long-period sdB binaries has been quantitatively reproduced by a population synthesis study.

We note that although the correlation of the sdB core mass and metallicity has been discussed in Chen et al. (2013), the period range obtained in their study (200–1050 d) does not match the range of the observed periods (800–1400 d). Based on our modelling, we infer that there are three primary sources of the difference. Firstly, the authors used [Fe/H] from −1 to 0 for the whole range of masses they considered, thus introducing systems with low metallicity of −1 and high masses of 1.6 M_⊙, which are very rare in the Galaxy and result in short final periods, cf. Figure 6. Secondly, the maximal initial mass for the progenitors was set to 1.6 M_⊙. As we show, higher-mass (1.9 M_⊙) solar-metallicity progenitors would explain the observations of systems with periods above 1050 d. Finally, the authors used several prescriptions for the angular momentum and mass loss to construct their fits, thus introducing a scatter when compared to the observed populations. With these differences accounted for, the Chen et al. (2013) model would probably result in a correlation which matches the observations more accurately.

Element diffusion in the atmospheres of sdB companions, although not included in our models, may potentially decrease their surface metallicities compared to the initial values by as much as 0.20 dex (Dotter et al. 2017). However, we expect that the material recently accreted from the red giants, chemically homogenised by convection, should have restored the initial surface metallicities of the companions. The fact that there is no visible offset between our model and the observations in Fig. 9 supports this idea.

A second possible observational confirmation of our model are the wide sdB binaries in the Galactic bulge. The bulge population is old and metal-rich, and the wide sdB binaries formed in this population are predicted to have on average longer orbital periods and a different slope in the P  −  q distribution than their lower-metallicity counterparts in the Galactic thin disc (see the green versus the blue populations in Fig. 7). Following the reasoning in Sect. 3, the narrower range of final mass ratios q in the P  −  q plane may be explained by a narrower range of masses of red giants in the bulge and hence a narrower range of masses of their companions compared to the thin disc population. A different slope of the correlation compared to that of the thin disc is related to the bulge being a mono-age population with a wider metallicity spread. Observations of wide sdB binaries in the bulge are difficult as the sources are faint and located in crowded fields, and there are currently no wide sdB binaries known in the bulge.

Both of the previous tests require the orbital period of the wide sdB binaries to be known. As these systems can have periods above several years, obtaining spectroscopic observations covering the entire orbit is a very time-consuming process. One may try targeting low-metallicity candidates with shorter periods, by selecting them kinematically from the thick disc or the halo, for example through the combined use of Gaia and LAMOST data. However, our models predict that the metallicities and the masses of the companion stars should also be correlated. This relation is shown in Fig. 12, and it does not require knowing the orbital period. The metallicity for this correlation can be obtained from a single high-resolution, high S/N spectrum. The companion mass can be determined from the surface gravity and radius or by fitting single stellar evolution models to the observed effective temperature, surface gravity, metallicity and radius (see e.g. Vos et al. 2018b; Maxted et al. 2015). The radius can be obtained using the Gaia parallax, while the other parameters can be obtained from one high-resolution, high-S/N spectrum. This method only requires one observation per target system and can thus easily be used to create a large test sample for the model. We do not expect any short-period sdB binaries to pollute the sample since there are no short period sdB binaries with FGK-type companions known (see e.g. Kawka et al. 2015). Similarly, we estimate chance alignments between sdBs and background or foreground MS stars to be very rare. An alternative way of obtaining the periods of sdB binaries is through high-cadence photometric variability studies, for example Otani et al. (2018), although this method is also time-consuming.

7.4. Gaia selection criteria

Based on the effective temperature, surface gravity and radius of the MESA models, we can calculate a theoretical colour-magnitude diagram for the wide sdB binaries. The spectroscopic parameters are taken to be time-averages over the core He-burning phase. TMAP (Werner et al. 2003) and Kurucz atmosphere models (Kurucz 1979) were used respectively for the sdB and the cool companion to calculate the synthetic magnitudes and colours. We show the resulting diagram in Fig. 15. The majority of the systems are located under the main sequence and can be easily selected from a colour cut in the Gaia bands. Some wide sdB binaries, however, overlap with the main sequence. These are systems with strong MS companions that outshine the sdB and are not detectable as composite binaries (shown in open circles) as well as systems that are spectroscopically recognizable as composite sdBs. The latter systems which are not apparent from a cut in the Gaia magnitudes could still be detectable if UV photometry from, for example, GALEX (Martin & GALEX Science Team 2003) is included. There are a few single-lined sdB binaries located at the RGB. These systems had an initial mass ratio very close to unity, and the companion has ascended the RGB during the sdB phase of the primary. The single-lined sdB+MS binaries dominated by the sdB star are located on the blue end of the diagram. These systems would photometrically and spectroscopically not look like a composite system but could be detected from RV variations.

Fig. 15. Gaia colour-magnitude diagram. Full blue circles show composite-spectrum sdB+MS binaries and empty blue circles show the single-lined sdB+MS binaries. sdA+MS binaries are similarly coded in green circles. The systems that will form an sdB+MS or sdA+MS binary are shown with orange squares at the moment when the interaction phase starts. As a comparison, the Gaia colour-magnitude diagram of the HIPPARCOS sample is shown in grey.

We can define a cut in the Gaia colour-magnitude diagram that includes all composite sdB binaries. The upper and lower limits on the absolute G-band magnitude are given in the following equations:

$$\begin{array}{cc}\hfill M_G& \le -2.06{(G_{\mathit{BP}}-G_{\mathit{RP}})}^2-1.35(G_{\mathit{BP}}-G_{\mathit{RP}})+5.10,\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill M_G& \ge -1.8(G_{\mathit{BP}}-G_{\mathit{RP}})+2.05.\hfill \end{array}$$(6)

The cut in the G_BP − G_RP colour is given by:

$$\begin{array}{c}\hfill (G_{\mathit{BP}}-G_{\mathit{RP}})\ge -0.08M_G-0.07\end{array}$$(7)

To exclude the main sequence stars, we also apply an absolute magnitude-dependent extra colour cut, following Geier et al. (2019), which is shown in red dotted line in Fig. 15:

$$\begin{array}{cc}\hfill (G_{\mathit{BP}}-G_{\mathit{RP}})& \le (M_G-1.84)/5.6,1.0<M_G,\le 3.8\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill (G_{\mathit{BP}}-G_{\mathit{RP}})& \le (M_G+1.83)/16.0,M_G>3.8.\hfill \end{array}$$(9)

After applying this main sequence Gaia cut on our simulated systems, 59 sdB binaries and 14 sdA binaries remain.

We show the progenitors of the sdB+MS binaries at the moment they start the mass transfer phase (Ṁ = 10⁻¹⁰ M_⊙ yr⁻¹) in Fig. 15 in orange squares. They are all grouped just underneath the red clump but are spread broader than the typical RGB track. These predictions make it possible to search for progenitor systems just before they start the RLOF phase. Comparing the population of progenitors with the population of the long-period sdB binaries may lead to further insights into the details of mass transfer in these systems, as well as the nature of the period-eccentricity relation observed in long-period systems. The location underneath the red clump strongly reduces the number of potential candidates and makes this a more feasible observational project. The progenitors of sdB+MS binaries may potentially already be present in the current datasets, for example Price-Whelan et al. (2020), which provide orbits of close binaries in the Galaxy.

7.5. Galactic formation rates and local observations

Our main sample of runs encloses the parameter space of all the systems formed through the Galactic history which may produce long-period composite sdB binaries degenerately at present day. Since our initial population is based on the Galactic star formation history, our set of MESA runs can be identified with the present-day formation rates of sdBs binaries, as we describe in detail in Sect. 5.3. By applying the method to the different Galactic components, we obtain the rates shown in Table 3. Specifically, the present-day formation rate of all the observable wide composite sdB binaries on the main P  −  q branch, contained in the Gaia cut (Fig. 15), is 1.2 ⋅ 10⁻³ yr⁻¹, which corresponds to 2.1 ⋅ 10⁵ systems currently present in the Galaxy. The primary sources of uncertainty for the numbers in Table 3 are the Galactic stellar mass, the binary fraction and, for the 500 pc sample, the Poisson noise, which altogether may lead to an uncertainty of a factor of a few. The uncertain stellar mass of the Galaxy may be factored out from the values in Table 3. In order to re-scale the values to a new Galactic stellar mass M_⋆, new or a new average star formation rate R_⋆, new over the Galactic age T_MW, the values should be multiplied by a factor M_⋆, new/(6.43 ⋅ 10¹⁰ M_⊙) or (R_⋆, new/(4.59 M_⊙ yr⁻¹)) ⋅ (T_MW/(14 Gyr)), correspondingly.

Despite the above uncertainty, our present-day formation rates are by 10–50 times lower than the estimate given in Han et al. (2003). One major factor for the difference comes from the fact that Han et al. (2003) used too high star formation rate for the Galaxy. Specifically, they assumed that there is one binary with the primary more massive than 0.8 M_⊙ forming in the thin disc and the thick disc per year (2 systems per year in total) at a constant rate over 15 Gyr. In comparison, the Besanćon model we use in this study predicts that such binaries form in the whole Galaxy at a rate of 0.17 yr⁻¹, that is about 10 times slower than assumed in Han et al. (2003). Part of this difference may come from the uncertain stellar mass of the Galaxy and part of it may be since Han et al. (2003) stated that they used an upper estimate on star formation rates. Additionally, we selected the systems based on them being observationally identifiable as sdB binaries, which may produce another factor of 2, as may be seen from Table 2. With these differences in mind, the rates broadly agree. We also have found the present-day formation rate for composite observable sdA binaries to be 0.27 kyr⁻¹, their total Galactic number being 0.43 ⋅ 10⁵, subject to a similar uncertainty range of a factor of a few, as the sdB production rates. Furthermore, assuming that the mass transfer phase of the red giant progenitor systems (shown with orange squares in Fig. 15) lasts on average about 10³–10⁴ yr, typical for the MESA runs, and that the sdB phase lasts about 100 Myr, we expect between several and a few tens of sdB progenitors currently in the mass transferring phase in the Galaxy.

We calculate the local number density of wide composite sdB binaries in Table 3 by using the local stellar number densities for the solar neighbourhood for each Galactic component from Robin et al. (2003) and by assuming that the metallicity history of the stars in the solar neighbourhood is representative of Milky Way as a whole. We obtain the total local density of long-period sdB binaries of 150 kpc⁻³. The thin disc population makes up for almost 97% of the locally observed systems, and 3% and 0.005% correspond to the thick disc and the halo, respectively. The negligible fraction of the halo objects is partially due to the stellar halo containing a small fraction of the Galactic mass and partially due to the halo being much more sparse compared to the disc population. By integrating the thin, thick disc and the halo density profiles from Jurić et al. (2008) and McMillan (2011), we obtain the expected number counts of the systems in the nearest 500 pc and 1000 pc. We observe approximately the same fraction of thick disc objects and also a negligible fraction of halo objects in these samples. We note that while the absolute rates and numbers in Table 3 are uncertain within a factor of a few, the relative rates and number counts of sdBs in different Galactic components are much more accurate and directly reflect our knowledge about the stellar mass fractions in each of these components.

It is important to note that, while the Besanćon model was constructed to represent the Milky Way Galaxy as a whole (Robin et al. 2003), it is still, nevertheless, calibrated to provide the best representation of the observed stars in the solar neighbourhood. Our observed sample of sdB binaries and the locations of wide-period sdB binary candidates also belong to the solar neighbourhood and, therefore, should be well modelled by the age-metallicity correlation provided by the Besanćon model, also in agreement with Casagrande et al. (2011) and Bensby et al. (2014). The Milky Way, as a whole, has a broader spread in the age-metallicity correlation compared to the solar neighbourhood, for example Feuillet et al. (2019). While this broader spread should not have much effect on the overall Galactic rates estimated in Table 3, we expect the P  −  q relation of sdB binaries to also become somewhat broader as larger parts of the Galaxy get included in the observed sample of sdB binaries.

7.6. Constraints from/on the Galactic evolution

We expect the metallicities of long-period sdB binaries to correlate with their vertical velocity dispersions in the Galaxy. Lower-metallicity binaries, for example, have formed earlier in the Galactic history and have experienced more dynamical interactions and thus have acquired higher vertical velocity dispersions compared to higher-metallicity longer-period sdBs which formed more recently. Based on the kinematics of white dwarfs, the vertical velocity dispersion for the oldest, lowest-metallicity sdB binaries in the thin disc is expected to be about 30 km s⁻¹ (Seabroke & Gilmore 2007) and their vertical scale height in the Galaxy is expected to be about 200 pc (Tremblay et al. 2016). The highest-metallicity sdB binaries, instead, originate from recently-formed stars and are expected to have vertical velocity dispersions of 10 km s⁻¹ and scale heights of about 75 pc.

The radial distributions of sdB binaries in the Galaxy depend on the radial metallicity structure of the Galactic disc and, implicitly, on the radial migration history of the Galaxy and, in particular, the churning and blurring processes in the Galactic disc, for example Frankel et al. (2019) and Feuillet et al. (2019). The fact that the Milky Way has a radial metallicity gradient of −0.07 dex kpc⁻¹, for example Robin et al. (2003), suggests that lower-metallicity and shorter-period binaries should at present, on average, reside at larger galactocentric radii than higher-metallicity longer-period binaries. While the long-period sdB binary samples are too small to constrain the spatial distributions of age and metallicity in the Galaxy, it is conceivable that larger samples of exotic binaries, including short-period sdB binaries, binary white dwarfs, symbiotic binaries, cataclysmic variables and others, may in a similar manner eventually prove complementary to the current Galactic chemical evolution studies.

For example, long-period sdB binaries may be used to put constraints on the thick disc properties. The location on the P  −  q diagram of the sdB binaries coming from the thick disc depends rather sensitively on the age assumed for the thick disc. For example, the age of 11 Gyr, assumed in Robin et al. (2003) places the thick disc binaries visibly outside of the observed P  −  q relation, whereas the age of 10 Gyr from the more recent study by Robin et al. (2014) lies closer to the relation, as seen in Fig. 7. As may be inferred from Table 3, constructing a complete sample of sdB stars within 1 kpc will also allow putting constraints on the Galactic mass fraction in the thick disc. For example, if the thick disc mass fraction was higher than 10–20%, the 1 kpc sample would show an overdensity of sdB binaries in the corresponding region in Fig. 7. We also do not expect any noticeable contribution for halo sdB stars in our observed sample, as shown in the previous section, which may be relevant to other types of sdB binaries (Luo et al. 2019).

Our study suggests that the P  −  q relation and, generally, the binary sdB population, should be different in environments with other metallicity histories compared to the Milky Way. This concerns, in particular, dwarf galaxies, other Milky Way-like galaxies, wherein the sdB populations may potentially be probed photometrically, and especially mono-age populations such as old open and globular clusters, wherein any P  −  q relation would have a different nature from that of the Galactic population. Similarly, our study implies that the sdB population has been rather different at different times in Galactic history. In the extreme example of the early universe, a formation channel similar to the one studied here would involve binaries formed from massive low-metallicity stars leading to bright, massive helium stars and possibly contributing to the reionisation of the universe (Götberg et al. 2017, 2020).

7.7. Constraining stellar evolution with sdB binaries

Successful match to the observations shown in Figs. 5 and 9, lends strong support to our model of red giant mass transfer onto main sequence stars. Since all the sdB binaries in the main branch of the P  −  q relation are the outcomes of stable mass transfer, our study serves as an additional observational confirmation that evolved red giants lead to stable mass transfer for initial mass ratios M_RG/M_comp of up to about 1.8, for primaries in mass range between 0.7 and 2.0 M_⊙. This result is in contrast to polytropic models used, for example, in the BSE code, which predict that systems that lead to M_sdB/M_comp <  0.5 should be unstable for mass ratios M_RG/M_comp above 1, that is always unstable in our scenario. This way, long-period sdB binaries can be used to constrain the stability of mass transfer from giant donors observationally.

Our models also validate that the red giant radii predicted by the MESA code with the default choice of the mixing length parameter are correct up to at least ∼10%. If the actual red giant radii were smaller or larger, our model, which does not contain any fitted parameters, would lead to proportionally shorter or longer periods for all the modelled sdB binaries, which would not match the observations shown in Fig. 5.

As the observed sample becomes larger, our models can potentially provide support to the observations of anti-correlation of metallicity and the close binary fraction in the Galaxy (Moe et al. 2019). The close binary fraction at metallicities of −0.6 is expected to be about 50% larger than at metallicities of −0.2 (Moe et al. 2019). The dependence of binary fraction on metallicity would be reflected in the relative number fractions of low-metallicity systems towards the short period end of the main branch and of the higher-metallicity systems towards the long-period end of the main branch in Fig. 5.

We have used the most standard binary and Galactic evolution prescriptions with a minimal number of model parameters. Since the determining quantities for the positions of sdB binaries on the P  −  q branch are the primary masses and the metallicities of the progenitors, we expect that the Galactic studies well calibrate these parameters. The main uncertainties of our model come from the less certain stellar mass fraction in the Galaxy, the variation of binary fractions and period distributions with metallicity and the inherent uncertainties in the standard binary mass transfer models. All these uncertainties should be reflected in the total rates of systems rather than the locations on the P  −  q plane.

Our modelling also reinforces the idea that sdB binaries are excellent tools for constraining binary mass transfer models. While we have aimed at using the most canonical binary mass transfer model, our simulations also suggest what the effects of the assumptions about binary mass transfer may be. In particular, if we allowed the accreting star to gain significant amounts of transferred material, the resulting mass ratios would be also be significantly reduced, producing a mismatch with the observed P  −  q relation. Since observations of the companion stars also indicate that they accrete only a few times 0.01 M_⊙ of material at most (Vos et al. 2018b), our study suggests that low-mass stars below about 1.5 M_⊙ with metallicities above −0.5 are generally inefficient at accreting material at mass transfer rates above 10⁻⁵ M_⊙ yr⁻¹, which are the typical rates for sdB progenitors.

This study does not constrain the correct angular momentum loss prescription. As discussed, for example by Rappaport et al. (1995) and Chen et al. (2013), the formation of sdB is strongly coupled to the core mass at the end of red giant mass transfer, which in turn couples to the red giant radius at the end of mass transfer and therefore to the final orbital period. In other words, the fact the red giant core ignites to produce an sdB is coupled to the period of the final binary, rather than to the angular momentum loss history which brought the binary to that period. However, since the red giant radius at the end of mass transfer is also sensitive to [Fe/H], the final periods are directly affected by metalicity.

Further studies of sdB binaries may help to constrain the angular momentum loss prescription. As may be seen from Eq. (2), angular momentum prescription directly affects which window of initial periods may lead to the formation of sdB binaries. Furthermore, angular momentum loss prescription dictates which fraction of binaries undergo stable mass transfer. These two effects may have a direct impact on the number of systems which can lead to sdB binaries. One can also potentially probe the angular momentum loss prescription as a function of progenitor mass, by studying how the number of observed systems varies along the P  −  q branch, although the number distribution of sdB-binaries along the P  −  q branch is degenerate with the Galactic star formation history. Finally, we point out that using the standard stability criteria of mass transfer was sufficient to explain the P  −  q relation of long-period sdB binaries. However, the choice of stability criteria may be further verified and constrained by studying through detailed MESA modelling of the populations of both short- and long-period composite sdB binaries simultaneously.

Since we have been able to explain the main branch of the P  −  q relation through degenerate ignition, in the sample of all the binaries with the primary masses between 0.7 and 2.0 M_⊙ formed through the Galactic history, the origins of the second P  −  q branch remain unclear. It seems unlikely to us that more massive stars could produce the second branch through non-degenerate ignition. The most massive progenitors on the main branch lead to q_final ≈ 0.3. Along the main branch, increasing the primary mass leads to lower final mass ratios. Naively extrapolating, the more massive non-degenerately igniting stars should end up having q_final ≲ 0.3, which is further enhanced by the fact that the final cores masses in the non-degenerate case are typically lower than in the degenerate case (unless the primary mass is larger than a few solar masses). However, a detailed binary evolution modelling of the non-degenerate channel is needed to make any conclusive predictions. It may also be possible that the second-branch objects are members of a peculiar Galactic sub-population with a metallicity of about [Fe/H] ≈ −1 and a broad range of ages between 5 and 10 Gyr not included in the Besanćon model, for example Chiappini et al. (2015), Martig et al. (2015) and Hekker & Johnson (2019). However, it seems unlikely that a separate unidentified population comparable in size to that of the thin disc would be present in the solar neighbourhood. Further insights into the Galactic membership of the second branch can be gained from a kinematic analysis of the long-period sdB population. Triple stellar evolution (Toonen et al. 2016) is unlikely to lead to a tight correlation seen in the second branch. Another interesting possibility could be that binary mass transfer prescriptions are also sensitive to metallicity. For example, if mass transfer were to become conservative at low metallicities, the binaries of the second branch could potentially be explained through the thick disc population, receiving their short periods due to low metallicity and the low final mass ratios due to the companion accreting additional amounts of mass. In any case, we expect that understanding the origins of the second branch will better inform us either about the Galactic populations or the binary evolution. Detecting possible systems in between the main and the second branches of the P  −  q relation may put important constraints on some of the above scenarios.

8. Conclusions

In this article, we have studied the relation between the orbital periods and mass ratios observed in long-period sdB+MS binaries. This relation was discovered by Vos et al. (2018a), where it was attributed to the stability of RLOF on the red giant branch. In this work, we have shown that while the stability of RLOF will influence the range in initial mass ratios that can form an sdB star, it alone is not responsible for the strong correlation between P and q. It is in fact the metallicity history of our Galaxy that causes the P  −  q relation. The slow increase over time of the average metallicity in the disc, where most of the sdBs originate, causes a correlation between sdB progenitor mass and metallicity. Due to the strong dependence of the red giant radius on the metallicity, this correlation causes low-mass progenitors with lower metallicities to form sdBs at shorter orbital periods than their heavier counterparts. The fact that the sdB mass range is very narrow, combined with the narrow range of mass ratios in the original binaries which lead to sdBs, creates the strong correlation between orbital period and mass ratio. A standard binary interaction model combined with the Besanćon Galaxy evolution model can match the main branch of the observed P  −  q relation very well. At the same time, while there are several formation scenarios possible for the second branch of the relation, its origins remain unexplained.

Our models show a strong correlation between the final orbital period and the initial metallicity of sdB+MS binaries, which is in very good agreement with the observations. The good match with the observations serves as further support for our model. This correlation also depends on the mass of the final product and holds both for sdB+MS and He-WD+MS systems. We provide a fitting formula for these correlations in Eqs. (3) and (4). Furthermore, there are strong correlations between initial and observed parameters that can aid in linking observed systems to their progenitors.

Our study also shows that it is possible to conduct a small but statistically significant population study with the MESA code. As we showed in Sect. 6, in the case of sdB binaries it is, in fact, necessary to use MESA compared to synthetic population synthesis codes such as BSE. Since presently, such codes do not capture the correct conditions for the formation of sdB stars. Using our models, we have determined the current formation rate for the long-period composite sdB binaries on the main branch of the P  −  q relation in the Galaxy to be 1.2 ⋅ 10⁻³yr⁻¹, and we can expect to observe about 240 such systems within 1 kpc. Moreover, we expect the population of wide composite sdB binaries in the bulge to be different from the currently observed sample, as we show in Fig. 7. Our study also provides strong support for the standard mass-loss model we used for mass transfer from low-mass red giants onto main sequence stars.

Our study has been the first to demonstrate how the chemical history of the Galaxy manifests itself in an observed binary population. We expect that using chemical evolution models of the Galaxy will have a significant impact in studying the observed binary populations, the production of which involves at least one low-mass (M_init <  1.6 M_⊙) component that has evolved off the main sequence. Such binary populations include hot subdwarfs, (extremely) low-mass white dwarfs, cataclysmic variables, ultra-compact X-ray binaries, symbiotic binaries, among others. In the longer term, we expect that the accurately measured properties of exotic binaries will provide new constraints on chemical evolution studies of the Galaxy.

Acknowledgments

This work was supported by a fellowship for postdoctoral researchers from the Alexander von Humboldt Foundation awarded to JV. We would like to thank Christian Sahlholdt, Diane Feuillet, Silvia Toonen, Jennifer Johnson, Chris Belczynski, Ashley Ruiter, Neige Frankel, Alina Istrate, Alejandra Romero, Ingrid Pelisoli, Paul McMillan, Astrid Lamberts, Melvyn B. Davies and Ross Church for helpful comments and discussions at different stages of this work.

References

All Tables

All Figures

Fig. 1. Formation channel for the main P − q branch of wide-period composite sdB binaries. Initially, a 0.7 M⊙–2.0 M⊙ main-sequence primary (MS) evolves to become a developed red giant (RG), which initiates mass transfer through Roche lobe overflow. The mass transfer rate gradually grows and exceeds 10−5 M⊙ yr−1, at which point the accretor cannot accept the incoming mass and ejects it in an outflow (shown schematically). Even though the red giant loses all its envelope, its degenerate core ignites helium and turns into a luminous sdB star. The companion main-sequence star accretes only a small amount of mass of ≲0.05 M⊙, but, nevertheless, becomes polluted (MS*). The radius of the RG is strongly sensitive to the metallicity of the primary, which in turn correlates with its age and mass due to Galactic evolution. This way, the Galactic age-metallicity correlation is an important factor determining the final periods of long-period sdB binaries.

In the text

Fig. 2. Example MESA track for a binary with a 1.2 M⊙ donor star, a 0.85 M⊙ companion, an initial period Pi = 350 days and a metallicity [Fe/H] = −0.15, which produces a hot composite sdB binary. Top panel: track for the donor star in the HR diagram. The mass loss phase is marked in red; the thin red line shows the phase when Ṁ > 10−9 M⊙ yr−1, thick red line corresponds to Ṁ > 10−6 M⊙ yr−1 and the moment of maximum mass loss is marked with a horizontal red marker. The first He flash is marked with a blue diamond and the core He burning phase is shown in a thick blue line. The phase between the He flash and the core He burning is very short compared to the core He burning phase. Bottom panel: orbital period and mass loss rate during the mass loss phase. The same markings are used as in the HR diagram: the thin black line shows the phase with Ṁ > 10−9 M⊙ yr−1, thick black line corresponds to Ṁ > 10−6 M⊙ yr−1 and maximum mass loss is marked with a vertical black marker.

In the text

	Fig. 3. Average effective temperature and surface gravity during the core He burning phase for the hot subdwarfs in the MESA sample. There is a strong relation between log g and Teff for hot subdwarfs formed from progenitors with degenerate He cores. Hot subdwarfs with Teff <  20 000 K are classified as sdA stars while the hotter subdwarfs are classified as sdB stars.
In the text

Fig. 4. Spectral energy distribution for an sdB+MS binary with 90% light contribution of the sdB in the Gaia-G band in the top panel, and with 20% contribution in the bottom panel. The binary SED is plotted in black while the contributions of the sdB and MS component are shown respectively in blue and green. In black-dotted line, the transmission curves for the Gaia-G, 2MASS-J and 2MASS-Ks filter are shown from left to right.

In the text

Fig. 5. Period-mass ratio relation of the sdB binaries produced by the different MESA samples. The grey squares show the real observed systems, while the blue circles show the MESA models which would have been observationally identified as composite long-period sdB binaries, i.e. fulfil the visibility criteria in Sect. 5.2. The observed systems, located within the nearest 1 kpc, belong to the solar neighbourhood. Most of the simulated sdB binaries in the Galactic population belong to the thin disc, which is the dominant Galactic component (see Table 2).

In the text

	Fig. 6. Initial metallicity as a function of the sdB progenitor mass for all systems in the model sample following the Galactic model (sample 4). Systems that fulfil the criteria to be recognizable as a composite sdB binary are plotted in filled blue circles, while those where only one component is visible are shown in open grey circles.
In the text

	Fig. 7. Schematic plot showing where wide sdBs formed from the different Galactic sub-populations are located in the period – mass ratio diagram. Thin-disc objects are shown in blue, thick-disc objects are shown in orange, bulge objects are shown in green and halo objects in purple. The observed systems are plotted in grey squares.
In the text

Fig. 8. Relation between the orbital period and the metallicity for the sdB binaries (left) and He-WD binaries (right) in our sample. The relation is a result of the donor radius on the RGB being strongly sensitive to its metallicity. The period-metallicity relation also depends on the sdB/He-WD mass, which is shown by the colour scale. The systems that are recognizable as composite sdBs are shown in filled circles while the single-lined ones are shown in grey open circles. For the He-WDs no distinction is made between composite and single-lined systems.

In the text

	Fig. 9. Observed versus modelled orbital period and metallicity correlation for long-period sdB binaries. The black squares show all the observed wide sdB binaries for which such measurements are available. The results from our MESA model following the Galactic model (model set 4) are shown in blue circles. Only model systems that fulfil the observability criteria are shown.
In the text

Fig. 10. Correlation between the mass ratios in the observable composite sdB binaries (filled circles) and the initial sdB progenitor masses in our models. The colour shows the initial mass ratios and the empty circles show the sdB binaries, which would not be observationally identified as composite systems. One may see that the observed mass ratios constrain, to a certain extent, the initial progenitor masses of sdB binaries.

In the text

Fig. 11. Correlation between the mass ratios in the observable composite sdB binaries (filled circles) and the initial orbital periods of the systems in our models. The colour shows the initial metallicity and the empty circles show the sdB binaries, which would not be observationally identified as composite systems. Given the rather tight correlation, one may constrain the initial period and, in certain cases, the progenitor metallicity from the present-day mass ratios.

In the text

	Fig. 12. Metallicity of the companion stars as a function of the companion mass after at the binary sdB formation time. The box plot shows the average metallicity with an orange line, the 1-sigma range with a box and the 3-sigma range with whiskers. It may be seen that the present-day metallicity and companion masses of sdB binaries are correlated, according to our model.
In the text

Fig. 13. Effect of the initial parameters: orbital period (Pi), mass ratio (qi), metallicity ([Fe/H]) and sdB progenitor mass on the final orbital periods and mass ratios of sdB binaries. We consider primary masses of 1.0, 1.2 and 1.8 M⊙, indicated in boldface below each set of tracks. Systems that result in hot composite sdB stars are marked with blue circles, composite sdA stars are marked with purple diamonds, while systems resulting in He-WDs or horizontal branch stars are marked with green triangles and orange squares respectively. The locations of the main and second branches of the observed P  −  q relation are indicated with the red dashed ellipses. The initial mass ratios and orbital periods are varied in equal steps in the range that still produces hot composite subdwarf stars for the given initial donor mass, while, for the metallicity, we show the one-sigma range of the Galactic metallicity corresponding to the time when the progenitors of each given mass formed (Table 1). This figure clearly shows that for a given initial donor mass, only a small range in the final mass ratios can be reached by varying the other initial parameters. Varying the metallicity has almost no effect on the final mass ratio, but strongly changes the final orbital period. The correlation between the final mass ratio and the initial donor mass combined with the link between metallicity and initial donor mass due to Galactic evolution then explains the observed P  −  q relation. See Sect. 6.4 for details.

In the text

Fig. 14. Mass distribution of the sdB stars in the Galactic sample. The distribution of the composite systems is shown with the blue line, while that of sdBs in single-lined systems is shown in grey. The bins were selected using Knuths rule (Knuth 2006). Both distributions are similar and consistent with the empirically determined masses, for example Fontaine et al. (2012) and Schneider et al (2020).

In the text

Fig. 15. Gaia colour-magnitude diagram. Full blue circles show composite-spectrum sdB+MS binaries and empty blue circles show the single-lined sdB+MS binaries. sdA+MS binaries are similarly coded in green circles. The systems that will form an sdB+MS or sdA+MS binary are shown with orange squares at the moment when the interaction phase starts. As a comparison, the Gaia colour-magnitude diagram of the HIPPARCOS sample is shown in grey.

In the text