© ESO 2020

1 Introduction

Stars form in smaller or larger clusters. Some of them remain today as open clusters or loosely bound stellar associations, but the vast majority of the star-forming regions have dissolved and spread their stars into the Milky Way field. One major aim of the GALactic Archaeology with HERMES survey (GALAH; De Silva et al. 2015) is to quantify this process and to trace individual stars back to the regions where they originally formed (Freeman & Bland-Hawthorn 2002; Bland-Hawthorn et al. 2010). A main tool for this is the so-called chemical tagging, where elemental abundances are used to trace stars with a common origin (examples where chemical tagging has been applied, developed, and tested include De Silva et al. 2007; Mitschang et al. 2014; Blanco-Cuaresma et al. 2015; Quillen et al. 2015; Liu et al. 2015; Ting et al. 2016; Kos et al. 2018; Andrews et al. 2019; Simpson et al. 2019). To this end, efficient and robust pipelines for the analysis of large numbers of single stars have been developed (Ness et al. 2015; García Pérez et al. 2016; Buder et al. 2018; Ting et al. 2019).

However, we know from studies of the stellar populations near to the Sun that a large fraction of stars is in binaries. To understand star formation, we need to know the properties of the binary population. Until now, the binary population in the field (i.e. not in stellar clusters, where photometric displacement in the cluster Hertzsprung-Russel, HR, diagram points to binarity) has been very difficult to study, and we are uncertain about the percentage of stars that are in binary systems. The fraction could be as low as 20% or as high as 80%, and might also vary with stellar type (see the review by Duchêne & Kraus 2013). The prime targets of the large spectroscopic surveys, FGK stars, can have a binary fraction as high as 50% (Raghavan et al. 2010). Additionally, binary population statistics are important for surveys targeting exoplanets, and they in turn contribute to understanding of binarity through serendipitous detections (e.g. Bluhm et al. 2016; Wittenmyer et al. 2020).

Input catalogues for spectroscopic surveys are built from all-sky photometric surveys (see e.g. Majewski et al. 2017; Stonkutė et al. 2016; Chiappini et al. 2019; Bensby et al. 2019), and they therefore include many as yet undiscovered spectroscopic non-eclipsing binary systems. This opens new unbiased possibilities for studying binary stars and their population statistics across the whole Milky Way and not in the solar neighbourhood alone. Because the Gaia-ESO Survey (Gilmore et al. 2012) reaches faint magnitudes (down to V ~ 19), it provides a good example of how we can extend the knowledge about the binary population far beyond the solar neighbourhood. Merle et al. (2017) found 185 confirmed double-lined binaries and 5 triple systems in the field populations and across the Gaia-ESO magnitude range (amounting to ~0.7 % of the analysed internal data-release 4 sample).

All of this means that to fully understand the stellar populations in the Milky Way, we also need to be able to analyse the spectra from binary systems that are observed with these large surveys. Most analyses of spectra belonging to binary systems have been made in dedicated studies with multi-epoch spectroscopy (see e.g. the collection of spectroscopic binary orbits by Pourbaix et al. 2004). Large spectroscopic surveys re-observe only a small and random fraction of their targets, primarily for internal consistency checks, and they offer limited opportunities for multi-epoch radial velocity studies (see e.g. Birko et al. 2019). Apache Point Observatory Galactic Evolution Experiment (APOGEE; Majewski et al. 2017), and to some extent, the RAdial Velocity Experiment (RAVE; Steinmetz et al. 2006) and Gaia-ESO Survey, are the only surveys that have multi-epoch spectroscopy that is rich enough for a time-series analysis (El-Badry et al. 2018; Skinner et al. 2018b; Matijevič et al. 2011; Merle et al. 2017). Because many surveys do not have multi-epoch spectroscopy, it becomes important to develop robust tools for analysing large sets of double-lined stellar spectra. In this study we develop such a method. Our approach is to first detect double-lined spectra and then use all available information for the associated stellar systems, such as spectra, photometry, and parallaxes, in a combined analysis. We develop a pipeline based on Bayesian inference that explicitly takes uncertainties in the observed data into account. The framework allows for the derivation of a large number of stellar properties. In this work we focus on deriving the stellar parameters, extinction, and overall metallicities. The inclusion of the parallax enables the derivation of radii of both stars in the binary. We apply our general framework to the GALAH survey.

The GALAH survey is an on-going spectroscopic survey, observing a million stars in the magnitude range 12–14 in V across the Southern sky. Data products from the analysis of the single-star spectra include elemental abundances for as many as 32 different elements (Buder et al. 2018). The selection of targets is simple: a magnitude cut is used, and the Galactic plane is avoided. The survey thus provides excellent opportunities for the observation of binary stars. Because the stars are relatively bright and the dwarfs therefore are nearby, they have excellent Gaia parallaxes (Gaia Collaboration 2018), which we exploit in the analysis.

Another non-trivial task is identifying the potential binary systems in a large survey data-set before their subsequent analysis (see e.g. Merle et al. 2017). Several approaches are possible, and we show that a combination of methods yields a larger catch of binaries. Westart with using the t-distributed stochastic neighbour embedding (t-SNE) classification to detect binaries (Traven et al. 2017), but also employ a method that uses cross-correlation with a spectral template (e.g. Merle et al. 2017). The two approaches are highly complementary.

This paper is organised as follows: Sect. 2 introduces the data, Sect. 3 describes the comprehensive framework that we develop, Sect. 4 describes our validation of the method; Sect. 5 describes and discusses our detection of binary stars in the GALAH data-set, Sect. 6 describes the cuts that we applied to obtain our final sample, Sect. 7 describes our results in some detail, Sect. 8 provides the discussion, and Sect. 9 contains our conclusions. A catalogue of the binaries and their derived stellar parameters is provided in Appendix A (the full catalogue is made publicly available at the CDS), and additional ancillary figures are shown in Appendix D.

2 Data

We here analyse a specific data-set: the extended GALAH data-set. This consists of stellar spectra from the GALAH survey (reduced as explained in Kos et al. 2017), apparent magnitudes from a variety of photometric catalogues (AAVSO Photometric All Sky Survey – APASS; Henden et al. 2016, Gaia DR2; Gaia Collaboration 2018, Two Micron All Sky Survey – 2MASS; Skrutskie et al. 2006, Wide-field Infrared Survey Explorer – WISE; Wright et al. 2010), and the parallax measurements from Gaia DR2. We applied a zero-point correction to Gaia DR2 parallaxes of +0.029 mas. This correction is based on the estimate that the parallaxes are on the whole too small, following from the Gaia investigation of quasars and validation solutions (see Lindegren et al. 2018).

Our analysis was performed on a per-spectrum basis (some targets have spectra from different epochs), therefore we chose the identifier sobject_id (see Table 5 in Buder et al. 2018) as a unique reference to each observed spectrum in the GALAH survey, accompanied by corresponding information that is included in the extended GALAH data-set. Hereafter, we refer to all investigated entities identified by sobject_id as GALAH objects, which means that repeated observations of the same target are treated independently and are compared only at the level of final products.

GALAH spectra are recorded in three visual bands and in one infra-red band (see Table 1), with a resolution of R ≈ 28 000 and a typical signal-to-noise ratio (S/N) per resolution element ~100 (Buder et al. 2018). The stellar spectra used in this study come from the GALAH internal data release 2 (iDR2). We only included spectra that successfully passed the reduction pipeline up to the point where the normalised spectrum, shifted to the rest frame, is produced (for more details, see Kos et al. 2017 and Buder et al. 2018). GALAH iDR2 also includes sources targeted by the GALAH pilot (Duong et al. 2018), K2-HERMES (Wittenmyer et al. 2018), and TESS-HERMES (Sharma et al. 2018) surveys, which are reduced using the same reduction pipeline as the main GALAH survey.

3 Method

In order to exploit all available information about the binary sources, we have designed our binary pipeline to incorporate constraints from spectroscopic, photometric, and astrometric measurements. The binary pipeline presented here is specific to the data-set we used, but it generalises to include all types of observational data and prior information. All observational data are fit simultaneously, which allows us to avoid converging to local minima. To provide a probability distribution function (PDF) for all parameters of interest, we employ Bayesian inference in tandem with a Monte Carlo Markov chain (MCMC) sampler (Gregory 2005; Hogg et al. 2010; Foreman-Mackey et al. 2013; Sharma 2017). This approach allows us to (a) define and incorporate a generative model of our data, (b) apply prior information about the model parameters, (c) explore a relatively large parameter space, and (d) estimate the model parameters and their statistical uncertainties.

3.1 Bayesian scheme

In order to estimate the model parameters and their uncertainties within the binary pipeline, we used Bayes’s theorem (see e.g. Hogg et al. 2010; Sharma 2017). We write the probability distribution of the model parameters θ (see Table 2) given the data X and their uncertainties σ_X, p(θ|X, σ_X), as $$p(\theta |X,{\sigma }_X)=\frac{p(X|\theta ,{\sigma }_X)\cdot p(\theta )}{p(X)},$$(1)

where $\mathcal{L}\equiv p(X|\theta ,{\sigma }_X)$ is the likelihood of the parameters or the probability of the observed data given the model, $\mathcal{P}\equiv p(\theta )$ is the prior knowledge about the model parameters, and p(X) is a constant (known as the evidence) that ensures that the posterior probability distribution p(θ|X, σ_X) is normalised to unity.

We assume Gaussian uncertainties for all measured quantities, and express the likelihood as the product of conditional probabilitiesfor each data point in the following form: $$\mathcal{L}={\displaystyle \prod _{i=1}^{N}p}(X_i|\theta ,{\sigma }_{X,i})={\displaystyle \prod _{i=1}^N\frac{1}{\sqrt{2\pi {\sigma }_{X,i}^2}}}\exp \left(-\frac{{\left[X_i-\mathcal{M}{(\theta )}_i\right]}^2}{2{\sigma }_{X,i}^2}\right),$$(2)

where N is the number of data points and $\mathcal{M}{(\theta )}_i$ represents the model value for a corresponding data point i given some independent variable (e.g. the wavelength in the case of spectroscopic data).

Similarly, when we have prior informative knowledge about the model parameters θ in the form of their value and standard deviation, σ, the prior can be written as $$\mathcal{P}={\displaystyle \prod _{j=1}^{M}p}({\theta }_j)={\displaystyle \prod _{j=1}^M\frac{1}{\sqrt{2\pi {\sigma }_{\theta ,j}^2}}}\exp \left(-\frac{{\left[{\theta }_j-{\theta }_{j,0}\right]}^2}{2{\sigma }_{\theta ,j}^2}\right),$$(3)

where M is the number of model parameters, while θ_j,0 and σ_θ,j are the informative value and standard deviation for a given parameter θ_j.

In this scenario, we seek to maximise $\mathcal{L}\cdot \mathcal{P}$, which is for computational convenience commonly implemented as maximising the natural logarithm of the same expression, yielding $$\ln \left(\mathcal{L}\cdot \mathcal{P}\right)=K-{\sum }_{i=1}^N\frac{{[X_i-\mathcal{M}{(\theta )}_i]}^2}{2{\sigma }_{X,i}^2}-{\sum }_{j=1}^M\frac{{[{\theta }_j-{\theta }_{j,0}]}^2}{2{\sigma }_{\theta ,j}^2},$$(4)

where K contains all the terms that remain constant in the process of inference or the exploration of the parameter (θ) space (the evidence term is independent of θ). Equation (4) shows that maximising $\ln \left(\mathcal{L}\cdot \mathcal{P}\right)$ is in our case identical to minimising a sum of different χ² terms.

3.2 Binary star model and prior information

In accordance with our philosophy of simultaneously applying information from the spectroscopic, photometric, and astrometric domain to characterise binary stars, our model of the data and prior information represent three distinct entities: (1) a model for the spectrum of the binary system, (2) a model for apparent magnitudes of the binary system, and (3) a prior on the parallax of the binary system. When we take this into account and remove the constant terms from Eq. (4), the objective function that we aim to minimise in our Bayesian inference becomes $$\begin{array}{l}{\chi }_{\text{obj}}^2={\chi }_{\text{spec}}^2+{\chi }_{\text{phot}}^2+{\chi }_{\text{astro}}^2,\\ {\chi }_{\text{spec}}^2=\frac{1}{2}{\sum }_{i=1}^{N_{\text{spec}}}\frac{{[X_{\text{spec},i}-{\mathcal{M}}_{\text{spec},i}(\theta )]}^2}{{\sigma }_{X_{\text{spec},i}}^2},\\ {\chi }_{\text{phot}}^2=\frac{1}{2}{\sum }_{i=1}^{N_{\text{phot}}}\frac{{[X_{\text{phot},i}-{\mathcal{M}}_{\text{phot},i}(\theta )]}^2}{{\sigma }_{X_{\text{phot},i}}^2},\\ {\chi }_{\text{astro}}^2=\frac{1}{2}\frac{{[\varpi -{\varpi }_0]}^2}{{\sigma }_{\varpi }^2}.\end{array}$$

Subscripts_spec,_phot, and_astro refer to the spectroscopic, photometric, and astrometric parts of the objective function and also indicate the source and type of data (for details, see Sects. 3.2.2 and 3.2.3). N_spec and N_phot represent the number of data points, X_spec,i and X_phot,i are individual values of data points, ${\mathcal{M}}_{\text{spec},i}(\theta )$ and ${\mathcal{M}}_{\text{phot},i}(\theta )$ are the corresponding values from the spectroscopic and photometric model, and ${\sigma }_{X_{\text{spec},i}}$ and ${\sigma }_{X_{\text{phot},i}}$ represent all pertaining uncertainties. The parallax ϖ is one of the model parameters θ, and ϖ₀ and σ_ϖ represent the prior information on the parallax and its standard deviation obtained from Gaia DR2. In an ideal situation, our models would be accurate representations of the physical systems we study, and all uncertainties pertaining to observed data would be properly accounted for. This is almost never the case, therefore the uncertainties ${\sigma }_{X_{\text{spec}}}$, ${\sigma }_{X_{\text{phot}}}$, and σ_ϖ must be scaled appropriately if they are under- or overestimated. The under- or overestimation of uncertainties in observational data, along with any systematic biases, can be estimated in the Bayesian scheme with additional θ parameters. However, doing so simultaneously for different data-sets while also including uncertainties of the spectroscopic and photometric model can be challenging. We plan to address this problem in a future study. For now, ${\sigma }_{X_{\text{spec}}}$, ${\sigma }_{X_{\text{phot}}}$, and σ_ϖ are included in the binary pipeline as reported in the respective sources of observational data.

3.2.1 Priors

All θ parameters except for the parallax are assigned uniform or flat priors, that is, we only fixed their allowed range. The values that we chose are given in Table 3. In the absence of strong prior information for these parameters, the flat priors used here do not hinder our inference because of the strongly identifiable likelihood. The prior on parallax is Gaussian, with the mean and standard deviation given by the Gaia DR2 parallax value and its uncertainty.

We set a generous prior for the absolute radius of a star because this quantity directly depends on the assumed distance to the binary system given by the inverse of the sampled parallax. Because matching the synthetic and observed apparent magnitudes effectively determines the angular size of the star (2 ⋅ R∕d) through luminosity and distance consideration, a potentially erroneous value of the parallax might produce unrealistically small or large radii, even though the other parameters are well determined. Therefore, we allow for such (rare) exceptions with a generous prior on the absolute radius and advise the user to confirm the quality estimates of the parallax before using the values of absolute radii (see e.g. the RUWE parameter in Table A.1).

3.2.2 Photometric model

Synthetic apparent magnitudes are obtained by convolving a model of the spectral energy distribution (SED) of the stars in the binary system with the desired passband transmission functions, taking the radius of each star and the parallax of the system into account. The SED of the binary is constructed by summing the stellar SEDs, after which reddening is applied. These computed apparent magnitudes form ${\mathcal{M}}_{\text{phot},i}(\theta )$ in Eq. (7). In addition, the SEDs are used to compute η_n as the luminosity ratio of the two stars (contribution of light from the secondary compared to the primary star) in a binary system for each of the four GALAH bands, $${\eta }_n={10}^{0.4[M_{p,n}-M_{s,n}]},n\in (1,2,3,4),$$(9)

where M_p,n and M_s,n are absolutemagnitudes of the primary and the secondary star in a given GALAH band (see Table 1).

We impose a common [Fe/H] for both stars in a binary (see the discussion in Sect. 8.3). Because both stars are on the main-sequence (MS), turn-off, or sub-giant branch, the more massive star always has the larger radius. Therefore we define the primary star as the one with the larger radius. We prefer this definition to the traditional one that is based on the luminosity in a given passband because when a star evolves away from the MS, its temperature drops and its radius increases, which may lead to a configuration in which one of the stars contributes more light at the red end of the spectrum, while the other star dominates in the blue. Figure 1 depicts such a binary system, where the secondary star contributes more light in the GALAH blue band, and less in the others. In the rare case of detecting binary systems harbouring two giants, the primary star by our definition might no longer be the more massive one because of the complex evolution on the giant branch. Results for such systems should therefore be handled with care.

The extended GALAH data-set includes observed apparent magnitudes from Gaia DR2 and three other photometric surveys (APASS, 2MASS, and WISE) that span the wavelength range from the visible to the infrared regime. These magnitudes and their uncertainties are converted into fluxes for all computations in the binary pipeline ($X_{\text{phot},i},{\sigma }_{X_{\text{phot},i}}$ in Eq. (7)). The synthetic fluxes (${\mathcal{M}}_{\text{phot},i}(\theta )$ in Eq. (7)) are based on the SEDs generated from the ATLAS9 stellar atmosphere models (Castelli & Kurucz 2004), using the Python packages PyPhot (Fouesneau 2019) and PySynPhot (STScI Development Team 2013) for their construction and convolution with given passband transmission functions.

The ATLAS9 stellar atmosphere models were chosen because they are incorporated in PySynPhot and because they agree excellently with the MARCS models for FGK stars, as discussed in Gustafsson et al. (2008). The latter is especially important because our aim is that the underlying physical models that we used both in the SME spectroscopic fit of the training set (see The Cannon model below) and in the photometric model described here are the same or very similar because some of the parameters (e.g. the effective temperature) are constrained by the spectroscopic and the photometric fit. Different extinction curves are available for PySynPhot, and we adopted mwavg, which is modelled for the Milky Way diffuse medium (standard R(V) = 3.1 extinction law by Cardelli et al. 1989).

The modelled SEDs (with T_eff, logg, and [Fe/H] as input) were convolved with the following ten photometric passbands¹ : B_J, V_J, G_BP, G, G_RP, J, H, K_S, W₁, and W₂ from the APASS, Gaia DR2, 2MASS, and WISE photometric catalogues, respectively. These passbands were chosen because corresponding observational data for GALAH stars is available, whose intrinsic quality is also discussed in Čotar et al. (2019), and becausetheir transmission curves are readily available in the PyPhot package.

Fig. 1 Spectral energy distributions for an example of a binary system (Teff,1 = 5018 K, Teff,2 = 5715 K, R1 = 1.9 R⊙, R2 = 1.3 R⊙, E(B− V) = 0.13). The legend explains which SEDs are plotted (P is the primary star, and S is the secondary star). The four GALAH spectral bands are indicated by vertical shaded regions, while the width of the transmission functions of photometric bands is indicated on the top. The black points are observed integrated fluxes in the following order of photometric bands: BJ, GBP, VJ, G, GRP, J, H, KS, W1, and W2 (see Sect. 3.2.2), and the short green lines indicate the synthetic integrated fluxes at the same position. The error bars of observed fluxes are exaggerated by a factor of 10 for clarity. Despite this, the uncertainty for the Gaia G value is smaller than the point size.

3.2.3 Spectroscopic model

In the spectroscopic model (Eq. (6)) we used spectral templates derived from observed spectra. This brings a number of advantages compared to using synthetic spectra: the majority of spectral lines are accounted for, the effects of the instrument are embedded automatically, the resolution is identical, and the spectral parameters directly relate to the main survey results. The downside is that some of the properties are less well determined, such as the exact continuum placement, but the continuum normalisation is effectively the same for the observed and modelled spectra (templates). Additionally, some parameters may be under-sampled, for instance, the higher rotational velocities, which are common to stars in close binary systems but not to their single counterparts.

The main GALAH pipeline incorporates The Cannon (Ness et al. 2015), a data-driven interpolation method for propagating astrophysical information within the same observed data-set of stellar spectra. For consistency with the main pipeline, we also adopt The Cannon as a generative model for the GALAH spectra. The Cannon model is based on the assumption that for a given spectrum n, we can approximate its flux value at some wavelength (f_n,λ) by a linear combination of the label terms pertaining to that spectrum, l_n, the corresponding coefficients for those terms, (Θ_λ), and the noise at that wavelength, summarised as $$\begin{array}{l}{f}_{n,\lambda }={\Theta }_{\lambda }^T\cdot l_n+{\text{noise}}_{\lambda }\\ l_n=(T_{\text{eff},n},\log g_n,{[\text{Fe/H}]}_n,T_{\text{eff},n}^2,\mathrm{...}),\end{array}$$

where the label terms can be any parameters and their combinations to an arbitrary order, and the noise term accounts for the observational noise (σ_n,λ) as well as for the inability of the model to reproduce the data (s_λ). The Θ_λ and s_λ finally represent The Cannon model, and they are obtained by fitting Eq. (10) to a training set (more details are provided in Ness et al. 2015).

For consistency with the main GALAH pipeline, we employed the same training set as was used to train The Cannon model there (see Sect. 3.2 in Buder et al. 2018). It consists of 10 605 spectra of dwarf and giant stars, but here we excluded 150 spectra that were recognised as binary star candidates in Sect. 5, and are thus left with 10 455 spectra, to which we refer hereafter as the training set. The parameters for the training set were estimated using the detailed spectrum synthesis code Spectroscopy Made Easy (SME; Piskunov & Valenti 2017), and the spectra were selected to be of the highest quality and to sufficiently cover the Hertzsprung-Russell diagram (for details, see Buder et al. 2018).

We produce spectral templates by feeding The Cannon model with a desired label vector (a set of parameters as describedin Table 2). These templates thus become the model part of Eq. (6) – ${\mathcal{M}}_{\text{spec},i}(\theta )$. The scatter of The Cannon model (s_λ) is added in quadrature to the uncertainty of observed flux values of GALAH spectra ${\sigma }_{X_{\text{spec},i}}$, and X_spec,i contains flux values from all fourGALAH spectral bands (see Table 1), which contribute approximately equally to the final number of spectroscopic data points N_spec of 12 230.

In the construction of spectral templates, we set the A_Ks parameter of The Cannon model to 0, unlike in the main GALAH pipeline, where it is fitted. The A_Ks parameter encapsulates the spectral signature of diffuse interstellar bands (see Sect. 3.3.1 in Buder et al. 2018), and if not nullified, would produce duplicated interstellar absorption lines in the modelled spectrum of a binary star because we combine two templates at different RV shifts. This parameter is therefore not listed in Table 2, and although it is related to interstellar absorption, it carries different information than E(B−V), which describes the influence of extinction on the overall SED.

It is difficult to determine the S/N of templates constructed by The Cannon model, but we estimate that it is much higher than the median S/N of the training set spectra (~ 200 per resolution element in the GALAH spectral band 2) because The Cannon effectively averages over them. The S/N of templates is also parameter dependent because it will change with the number of training set spectra in a given bin of the parameter space. The quality of the fit will depend on the S/N of templates, but also on the ability of The Cannon to accurately reproduce the training set spectra, which in turn also depends on the complexity of The Cannon model. Most of the literature employs The Cannon model with labels in the second order, arguing that the level of complexity thus obtained is sufficient. We follow this here, largely because we wish to be consistent with the main GALAH pipeline, but we nevertheless performed a brief investigation of proceeding to the third and fourth order in labels that can be found in Appendix B.

The final synthetic spectrum of the stars in a binary system is constructed by summing together the model spectra for the primary and secondary star, using the luminosity ratio η_n. Before summation, both the primary and secondary spectra are shifted according to their radial velocities (RVs) given by θ.

We identified a specific problem when The Cannon model was used to produce both spectral components of the synthetic binary template. In the regime of higher v_broad (e.g. > ~20 km s⁻¹) and with a certain combination of other parameters, all within the prior ranges given in Table 3, the model can produce double lines in a single-component spectrum. This leads to a rejection of solutions for genuine double-lined spectra at the low ΔV_r end of the distribution of binary systems, with the primary component mimicking double lines with small velocity separations and the secondary component ending up aligned with the primary and thus not satisfying the minimum ΔV_r for a reliable solution.

We have checked the relatively small number of spectra of the training set with high values of v_broad and cannot safely reject any of them as binaries, which might resolve this problem. We therefore imposed a more strict prior in the analysis of binary spectra by changing the upper limit of the flat prior on v_broad to 15 km s⁻¹ (which still includes the majority of the training set v_broad values). This limit therefore represents the highest v_broad of the modelled spectrum, and can limit the accuracy of the spectroscopic fit in the hotter regime of the HR diagram or in the case of close binaries with synchronous rotation (that are normally spun up), but it does not significantly influence the majority of our results. Because our solution is not ideal, further attention tothis problem is required in future work, while in the current catalogue of results, values of v_broad close to 15 km s⁻¹ should be considered with care.

3.2.4 Posterior probability distributions of model parameters

The model parameters θ and their uncertainties σ_θ are obtained by sampling their posterior probability distributions based on Eq. (5) and using the emcee code, an implementation of the affine-invariant MCMC ensemble sampler (Foreman-Mackey et al. 2013, version 2.2.1). This scheme offers a way to explore a given parameter space and very elegantly integrates (marginalises) posterior probability distributions.

We report the results for the model parameters θ by the statistics of their posterior probability distributions, each of which contains 25 600 samples (64 walkers times last 400 steps of the MCMC chain). We calculated the mean value and extracted the 16th, 50th (the median), and 84th percentile for each posterior distribution. Additionally, we provide the mode as the centre of the largest bin by binning the posterior distribution in 30 equal bins.

Because the MCMC scheme automatically provides a sampling of posterior probability distributions, which can be carried forward to subsequent inferences (useful for propagating uncertainty), or standing in as a prior, we used this feature in subsequent steps of our method, as discussed in Sect. 3.2.5. We list the model parameters in Table 2. The additional hidden parameter η_n is listed separately at the bottom of Table 2 because it is not directly sampled from the posterior probability distribution, but it is determined by the photometric model and used by the spectroscopic model.

3.2.5 MCMC sampling process

By trial and error, we developed an MCMC procedure for posterior sampling that produces converged solutions close to global optima while minimising computation time. The MCMC process for each system was initialised with 1024 walkers in order to cover the relatively high dimensional parameter space. The initial parameter values for walkers were assigned randomly from the prior distributions (see Table 3). First of all, 250 steps were used to obtain priors on RVs of both binary components, using only the spectroscopic fit. Hereafter, the full model was used, starting with 100 steps, after which the upper 10% of best solutions were kept, and from then on, the MCMC walk continued using 64 walkers. In this second part, 750 steps were made, after which the upper 25% of best solutions were kept, followed by 500 more steps. At this point, the convergence of solutions was checked using a custom convergence criterion, and if convergence was not achieved, the MCMC walk continued in batches, each of which contained 500 steps and a convergence check. Taking only the upper best solutions of subsequent batches potentially excludes multimodal solutions, but close to convergence, the results are not highly multimodal. The process stopped either with achieved convergence or a maximum of ten repeated batches.

The convergence criterion consists of two custom rules for comparison of the scatter and slope of the MCMC chain in subsequent batches, which proved to be efficient and faster than standard convergence diagnostics such as auto-correlation. The first rule checks that the median values of ${\chi }_{\text{obj}}^2$ at the start and end of the batch do not differ by more than one standard deviation of ${\chi }_{\text{obj}}^2$ at the end. The second rule effectively checks that there are no multiple chains by comparing 10% of the best and 10% of the worst ${\chi }_{\text{obj}}^2$ at the end of the batch, imposing that their average standard deviation multiplied by 20 is always larger than the difference of their medians.

The MCMC setup used in this work takes on average 8 CPU hours to arrive at a converged solution for each observed spectrum. Evaluating ${\chi }_{\text{phot}}^2$ is approximately ten times slower than evaluating ${\chi }_{\text{spec}}^2$; the computation of SEDs using PySynPhot is by far the most expensive calculation.

4 Validation of the method

We first validated our method for deriving the properties of stars in a binary system by applying it to single stars. This makes the process of validation transparent and comprehensible because the single-star model is less complex, and also lets us benchmark selected derived values against external sources. The photometric and spectroscopic models change such that the synthetic apparent magnitudes are derived using only one SED, and the synthetic spectra are composed of only one component. We used different sets of GALAH objects for two different validation tests: The Cannon validation sample (Sect. 4.1), and the benchmark validation sample (Sects. 4.2 and 4.3).

4.1 Spectroscopic validation

We verified that the binary pipeline, when using only the spectroscopic fit, produced similar results to the main GALAH pipeline. We expect this to be the case because (a) the generative model for stellar spectra is produced by the same technique (The Cannon), (b) the training set for The Cannon is the same as was used in the main GALAH pipeline, and (c) the objective function for comparing the model to data (see Eq. (6)) is essentially the same as in the main GALAH pipeline. The only significant difference is the process of minimisation of the objective function. The test was made on The Cannon validation sample, which is a subset of ~2000 random stars from the training set. This sample covers the HR diagram in the same way as the training set, but more sparsely. It contains dwarf as well as giant stars, and their positions on the sky are indicated in Fig. 13.

Figures 2 and B.1 show the comparison of SME values (see Buder et al. 2018) to the results of both the binary pipeline (Fig. 2) and the main GALAH pipeline (Fig. B.1), confirming their close agreement. We observe some trends and biases in these figures, such as the inconsistency for the lowest metallicities in the bottom panel of Fig.2, which are a consequence of both the flexibility of the The Cannon model and the adequacy of the training set,where it is difficult to ensure equal representation of spectra in all bins of the modelled parameter space (see Sect. 3.3.1 in Buder et al. 2018).

4.2 Photometric validation

In this validation check and in Sect. 4.3, we used the benchmark validation sample, which is a subset of GALAH objects withexternal information about the temperature from the infrared flux method (IRFM; Casagrande et al. 2019) and radii (asteroseismology; Stello et al. 2017; Sharma et al. 2019). This sample of 620 giant stars is therefore suitable for benchmark testing, and is dominated by stars from K2 campaigns C1 and C6 (Sharma et al. 2019, see also Fig. 13). Compared to the previous spectroscopic validation, this sample has almost no overlap with the training set and hence no SME parameters to compare to.

We took the results from the binary pipeline, using only the photometric fit (Eq. (7)), and compared them to the IRFM temperature values and reddening values from SFD (Schlegel, Finkbeiner, David) reddening maps (Schlegel et al. 1998) and 3D reddening maps based on Pan-STARRS 1 photometry (Green et al. 2018; hereafter the Bayestar maps). This comparison is shown in Fig. 3, where we display results with and without the Gaia magnitudes in the photometric fit. We distinguish the Gaia magnitudes from other photometric passbands used in this work because their uncertainties are approximately an order of magnitude lower. They therefore affect the fit more strongly.

Figure 3 shows some important problems. The most striking is a bias towards higher temperatures as compared to both the GALAH DR2 values and those from the IRFM. An additional test that we performed using a larger validation sample (without the need for external information about the radii) also shows strong grid effects, which diminish considerably, however, when we exclude Gaia magnitudes. The grid effect has a spacing of about 250 K, and it might originate from the SED model in combinationwith small uncertainties of Gaia data. Excluding Gaia magnitudes from the fit also slightly reduces the bias towards higher temperature and reddening values. In any case, this effect is caused by the degeneracy between the two parameters, where an increase in the modelled temperature is compensated for by the increased extinction (and therefore reddening). More details of this comparison, which is summarised in Fig. 3, are shown in Fig. E.1, where thex-axis of the panels spans the literature values to which we compare our results. The grid effect shown in Fig. E.1 is reflected in the bimodality of distributions in Fig. 3.

The photometric validation test implies that our photometric data alone cannot constrain the shape of the SED to an extent where we could produce reliable estimates of the interstellar reddening and temperature of the stars. The photometric fit therefore has to be complemented by another constraint on the temperature (or reddening), which in this work is provided by the spectroscopic fit.

Fig. 2 Comparison of SME and The Cannon labels for The Cannon validation sample. The Cannon values produced by the binary pipeline are subtracted from the SME values.

Fig. 3 Comparison of Teff (orange) and E(B−V) (green) parameters as derived by the binary pipeline (using only the photometric fit) to the literature values (indicated on the horizontal axis). The literature values are subtracted from those given by the binary pipeline. The lighter shaded distributions correspond to results from the binary pipeline when the Gaia magnitudes are excluded from the fit.

4.3 Overall validation for single stars

We now investigate the results by fusing the spectroscopic and photometric part of the binary pipeline and the prior on parallax. By including the photometric fit in the binary pipeline, we can determine the radii and reddening for our target stars. The results for stellar radii do not depend significantly on the degeneracy between the temperature and extinction. However, this degeneracy is detrimental to determining E(B−V), as illustrated in Fig. 3. It arises from the interplay between the effects of temperature and extinction on the modelled SED. By increasing the temperature, the blackbody SED experiences a relatively larger increase of flux in the blue part of the spectrum. This can be compensated for by an increase in reddening (extinction), which decreases the flux more strongly at shorter wavelengths. Thus the approximate shape of the SED and the apparent synthetic magnitudes are conserved, while the modelled radius is largely unaffected.

Figure 4 demonstrates that the bias towards higher temperatures and reddening, which was observed in Fig. 3, is now largely overcome with the aid of the spectroscopic part of the fit in the binary pipeline. Figure 4a, which compares the temperature values, is consistent with Fig. 16 from Buder et al. (2018), indicating a good agreement for the temperature parameter between the binary pipeline and the main GALAH pipeline. The reddening parameter shown in Figs. 4b and c is still somewhat overestimated, however, and in the comparison with SFD reddening maps, there seems to be a better agreement at higher extinction regions than in the case of Bayestar maps. Nevertheless, we might expect the Bayestar maps to be more realistic because they are in 3D, and thus it would be plausible that the binary pipeline overestimates reddening in low and high extinction regions, as is indicated in Fig. 4b.

Figure 5 shows a comparison of radii derived by the binary pipeline, Gaia DR2, and those computed from asteroseismic parameters ν_max and Δν using scaling relations adopted from Huber et al. (2017, see their Eq. (3)). We applied a Δν correction following Sharma et al. (2016). The constants ν_max⊙ and Δν_⊙ are taken from Huber et al. (2017), whereas T_eff,⊙ is set to 5772 K (Mamajek et al. 2015) and the T_eff is the one computed by the binary pipeline. The Gaia DR2 radii are taken as is, computed based on luminosity (obtained using temperature-only dependent bolometric correction) and effective temperature as determined by Gaia (for more details, see Andrae et al. 2018). We did not correct Gaia radii for extinction, the parallax zero-point, or the difference in temperaturescales employed by Gaia and this work. The comparison shown in panels c and d in Fig. 5 should therefore be taken as approximate because thorough corrections have to be applied for a more faithful benchmarking against Gaia radii, such as the one performed by Zinn et al. (2019) for seismic radii of Kepler field stars.

Because we did not correct for Gaia DR2 radii, both the binary pipeline and seismic values show an offset towards larger radii in our crudecomparison. In the case of seismic – Gaia comparison, the discrepancy is thus larger than the reported agreement for Kepler field stars (Zinn et al. 2019). However, the agreement between the binary pipeline and seismic values is encouraging, showing a negligible bias and a median relative difference in radii of 5.7%. The stars with the largest absolute differences (ΔR) in all panels of Fig. 5 also have larger relative parallax uncertainties (σ_ϖ ∕ϖ; up to 40%).

In summary, the validation of the binary pipeline on single stars indicates that we can achieve results that are consistent with selected literature sources to ~45 K in effective temperature, ~6% in stellar radii, and ~40% in reddening (measured in magnitudes), keeping in mind that this holds true for a subset of giant stars that constitute the benchmark validation sample. We do not speculate on how this comparison translates into accuracy because different methods can be more reliable for different sub-populations of investigated objects. However, in the case of binary stars, where two components constrain the overall shape of the SED or the combined spectrum, we expect less degeneracy between model parameters and thus lower uncertainties on temperatures, radii, or reddening, for instance.

Fig. 4 Temperature and reddening for the benchmark validation sample derived by the binary pipeline compared to the literature values from the IRFM method and Bayestar and SFD reddening maps (see also Fig. 3). The colour-coding is by absolute Galactic latitude.

Fig. 5 Comparison of radii for the benchmark validation sample derived by the binary pipeline, Gaia DR2, and computed from asteroseismic parameters (more details in Sect. 4.3). The colour-coding is by relative parallax uncertainty except in panel b, which is a copy of panel a, but colour-coded to indicate stars with less reliable seismic information (νmax < 40; S. Sharma, priv. comm.). The more reliable part of the sample is used to calculate the bias and scatter as reported in panels a, c, and d. The colour scale in panel a is clipped for clarity, while panels b and c have a smaller range of σϖ ∕ϖ because some Gaia DR2 values for radii were lacking.

4.4 Validation of the binary pipeline on artificial binary stars

In order to estimate the performance of the binary pipeline in producing reliable parameters for SB2 systems, we constructed an artificial sample of 2000 binary systems based on single stars from the training set. This sample was produced in a simplified approach such that the properties of both stars roughly represented realistic combinations of stars in binary systems. A more thorough consideration of stellar properties based on the mass, mass ratio, and evolutionary stage of the binary system, for example, is beyond scope of this validation, and we are additionally limited by the choice of stars in the training set.

The simplified approach of combining the primary and secondary star differs in the way we treat dwarf and giant systems, and we only produced dwarf-dwarf and giant-giant combinations (we followed the dwarf-giant distinction given by Zwitter et al. 2018). For every star inthe training set (the primary), we imposed three conditions for possible secondary candidates: (1) a matching metallicity (|Δ[Fe∕H]| < 0.02), (2) the surface gravity (logg) of the secondary must be higher than or equal to that of the primary, and (3) T_eff2 ≤ T_eff1 for dwarf-dwarf pairs and T_eff2 ≥ T_eff1 for giant-giant pairs. We then randomly selected five candidates for the secondary among the training set stars that satisfied the three conditions.

To obtain the necessary properties for combining primary and secondary spectra, we assigned a mass ratio to each pair. For giant-giant pairs we assumed a mass ratio of one because both stars need to have approximately equal mass in order to simultaneously evolve to the giant stage, whereas for dwarfs, we adopted an empirical relation between mass and effective temperature ($m_2/m_1={[T_{\text{eff2}}/T_{\text{eff1}}]}^{1/0.7}$) based on the study by Eker et al. (2015). The primary was then assigned a radius of 1 R_⊙, while the radius of the secondary was computed from the previously determined mass ratio and the measured surface gravity of both stars. We set a parallax of 1 mas for all systems, after which artificial observed magnitudes, together with luminosity ratios for combining spectra of both components, were computed using the same tools as presented in Sect. 3. When the luminosity ratio in any GALAH spectral band was below 0.1, the system was rejected. The E(B−V) value for each system was then drawn randomly from an exponential distribution (scale = 0.15) that closely reflected the overall distribution of reddening for GALAH targets.

The primary star in each pair was kept at rest frame, whereas we shifted the secondary by randomly sampling its RV from an exponential distribution (scale = 50 km s⁻¹) with a random sign. Then we combined the spectra of the two components. Finally, we randomly chose 1500 out of ~ 140k produced pairs with 15% giant-giant and the remaining dwarf-dwarf systems.

The results of the comparison between the input (SME parameters, additional assigned properties) and the parameters recovered by the binary pipeline are characterised in Table 4. The overall distributions largely resemble the uncertainty of the spectroscopic model, as illustrated in Figs. 2 and B.1, with slightly amplified biases that indicate the degeneracy between the primary and secondary star. The scatter of [Fe/H] in this validation is smaller than that in Fig. 2 and might result from the constraint on this parameter by both components.

5 Identification of multiple star candidates

A number of techniques to identify stellar binaries from their integrated spectra have been developed, mostly based on the analysis of the cross-correlation function (CCF). We complemented this technique with a sample of binaries detected by applying a dimensionality reduction algorithm to classify the spectrum. We identified binaries using the procedure described by Traven et al. (2017) and updated in Buder et al. (2018), and confirmed and extended the sample with a CCF analysis following Merle et al. (2017).

5.1 Detection through a dimensionality reduction

The detection and selection of multiple star candidates was first performed with the help of a dimensionality reduction technique, the t-SNE, which is very efficient at grouping similar objects in a visually manageable 2D map (Fig. C.1). The selection process includes a clustering algorithm, DBSCAN (Ester et al. 1996), that alleviates the detection of boundaries for differentmorphological groups of spectra in the t-SNE map, but it can also be complemented by manual selection of smaller regions or individual spectra. We refer to Buder et al. (2018) for more details on the iterative semi-automatic classification procedure.

We selected the spectra that were flagged as binary stars in Buder et al. (2018), but we also added a few hundred spectra that were classified as binaries in a previous study (Traven et al. 2017), which are not flagged as such in the new classification. This already implies that our classification procedure is not perfect (see discussion in Traven et al. 2017), but the majority of the double-linedbinary stars are detected. A group of binary stars that is generally harder to isolate in the map are the most massive and therefore hotter stars, where the spectral features caused by the high temperature (strong, broad lines) dominate the binary signature (double lines). The hotter binary stars are better accounted for in the second detection procedure discussed in the following section.

5.2 Detection through the cross-correlation function

The second selection of multiple star candidates was performed with the CCF technique and automated identification of multiple RV components developed by Merle et al. (2017). Because the observed stars are mainly of FGK spectral type (and most of the GALAH binaries are nearby dwarf systems), we used a single spectral template for the CCF computation following Matijevič et al. (2010). We used the MARCS model atmosphere (Gustafsson et al. 2008) of the Sun (T_eff = 5777 K, log g = 4.44 and [Fe/H] = 0) produced by the 1D local thermal equilibirum (LTE) radiative transfer code Turbospectrum (Plez 2012).

The most prominent absorption lines in GALAH spectra are the hydrogen lines Hα and Hβ, which are usually either saturated or very strong in late-type stars. They produce a very broad peak in the CCF that can hide close secondary components. We show an example in Fig. 6, where we plot CCFs that are computed with (dashed lines) and without (plain lines) Hα (band 3) and Hβ (band 1). By excluding Hα and Hβ from the template, the typical width of a single RV component CCF is similar in all four bands of the GALAH spectra because bands 2 and 4 do not include saturated or strong absorption lines. Hence we chose to exclude the hydrogen lines when we computed the CCF.

The solar templates (one for each band) used in the computation of the CCFs were calculated by Turbospectrum with a wavelength step equal to half of the original sampling of GALAH spectra (0.020, 0.025, 0.030, and 0.035 Å), corresponding to resolving powers larger than 220 000. The templates were not convolved to the resolution of the instrument. Each CCF was computed over a RV range of 700 km s⁻¹ centred on the zero velocity with a velocity step of 0.5 km s⁻¹. The number and positions of RV components were computed using the code detection of extrema, DOE² (see Merle et al. 2017). For each CCF, the number of components were obtained using the first three derivatives of the CCF. The use of the first derivative does not allow us to automatically identify two close components when a local minimum between the two components is missing. However, the minima of the second derivative provide the positions of the CCF peaks. The ascending zeros of the third derivative were used to identify these positions. When the number of components were identified, we performed a Gaussian fit of the core of each component. Merle et al. (2017) provide more details on the process. An example of an SB2 candidate for which the first derivative of the CCF is not sufficient to detect the second component is given in Fig. 7. The RVs obtained by the Gaussian fits are slightly different than the zeros of the third derivative because the successive derivatives are smoothed (by a Gaussian kernel with σ =8 km s⁻¹) to decrease the numerical noise.

Fig. 6 Effect of excluding Hα (in band 3) and Hβ (in band 1) in the solar template on the width of the CCF. The decrease in CCF peak width reveals a second close component that is difficult to detect when Hα and Hβ are includedin the template.

5.3 SB2 candidate sample from combined detections

We combined the SB2 spectra detected with t-SNE and CCF and obtained a final list of 19 773 unique candidates, where the intersection between t-SNE and CCF detections amounts to 7556 objects. The CCF approach reproduces the majority of t-SNE identifications in the region of isolated groups of binary stars (blue points in Fig. 8), but some candidates well inside these groups and most of those in the transition regions are missed. On the other hand, the CCF approach additionally contributes 6841 candidates from different parts of the parameter space (red points in Fig. 8). These identifications with CCF are, as expected and clearly shown in Fig. 8, mostly outliers or edge members in groups of spectra from distinct regions of the parameter space.

The colour-magnitude diagram in Fig. 9 displays all t-SNE and CCF detected SB2 candidates against the background of all GALAH spectra considered here. We can distinguish the binary main-sequence, which is offset to brighter absolute magnitudes. The over-density of detections just below this binary main-sequence shown in Fig. 9b points to false positives due to a smooth transition from well-distinguished double lines to a mixture of unresolved binary and single-star spectra. Figure 9a also shows that CCF extends the detection towards cooler and hotter temperatures than t-SNE.

Fig. 7 Identification of the RV components in a CCF of a close SB2 candidate. The top panel shows the CCF. The three panels below correspond to the three successive derivatives. The vertical black lines show the RV components detected by taking the ascending zeros of the third derivative (bottom panel). The vertical red lines give the positions of RV components obtained by fitting the multi-peaked CCF.

6 Selecting the final sample of SB2s

In this section we define the final sample of binary stars for our investigation. We consider different aspects of the results produced bythe binary pipeline, and summarise the four criteria that we used to define the final sample (Table 5). It is important to recognise that these selection criteria by themselves produce a selection effect, adding to others, for example, from the identification process described above. After applying the selection criteria to the sample of 19 773 analysed objects, we are left with12 760. From this final sample, we show two examples of the fit performed by the binary pipeline in Fig. 10.

6.1 Selection criteria

The main aim of the selection presented here is to reject physically or technically erroneous solutions. We note that the cuts based on these selection criteria significantly reduce the number of binary stars in our sample. This is mostly due to false identifications and the conservative nature of our selection criteria. Below we discuss the four criteria listed in Table 5 in more detail.

Fig. 8 t-SNEprojection map of 587 153 spectra as presented in Buder et al. (2018). The points (spectra) are colour-coded by all t-SNE detected SB2s (blue) and the unique additional detections by CCF (red). The count of spectra for each group is given in the legend.

Criterion 1

The convergence of the solution to the global minimum in the binary pipeline, or at least to a minimum as defined by our custom convergence criterion (Sect. 3.2.5), can fail for several reasons. By letting the MCMC sampling process run much longer, more objects would arrive at convergence, but it is not guaranteed that these would really be good solutions because convergence issues are often connected to data of poor quality, multiple minima, and thus degenerate solutions, and other issues that might make a lengthy convergence doubtful in any case. When we apply each criterion from Table 5 on the full sample instead of in sequence, we find that the overlap of objects rejected by criterion 1 and the others is large, indicating that criterion 1 is connected to specific issues that the rest of them address.

Fig. 9 Colour-magnitude diagram of SB2 identifications. (a) All identifications with CCF, (b) all identifications with t-SNE. The grey points are all GALAH objects considered in this work. The absolute magnitudes and colours are taken from the Gaia DR2 catalogue, and no reddening corrections have been applied.

Fig. 10 Example fits of binary stars from the final sample. Top and bottom panels: binary system at low- and high-velocity separation of double lines, respectively. The spectral region around Hβ and Hα is shown in the left and middle, respectively, while a region of a few well-isolated lines is shown on the right. Some of the derived θ parameters are indicated for each binary system, and the bottom example shows that the primary star has already evolved from the main-sequence. The primary and secondary component of the modelled binary spectrum are shown separately at a continuum level in accordance with their contribution. The black error bars in the inset plot of the photometricfit represent the relative uncertainty of the observed fluxes in the indicated photometric passbands, while the red points represent the relative offset of synthetic fluxes. The error bars of Gaia photometry are multiplied by a factor of 10 for clarity.

Criterion 2

The restriction on the difference of RVs of the two stars, or the velocity separation of double lines, is based on an estimate of the limit to which GALAH spectra can be used to resolve almost blended double lines. It is difficult to quantify the smallest reliable velocity difference where we can be certain that double lines are indeed present, and it depends on the flux ratio, stellar type, stellar rotation, varying resolution across the spectral range, and other effects. By visual inspection of the ${\chi }_{\text{spec}}^2$ fits and the degeneracies of stellar parameters in the case of severely blended lines, we set a conservative limiting value of 15 km s⁻¹. An example of partially blended lines at this velocity separation is shown in Fig. 10. Nevertheless, some line profiles that seem double-lined but heavily blended can originate from spotted surfaces or pulsation (either radial or non-radial), which is present in stars from many places of the HR diagram, causing diverse line asymmetries.

Criterion 3

In the case of spectra where the lines of the secondary star are very weak compared to those of the primary star, the spectroscopic fit of the secondary component becomes unreliable. For a main-sequence binary with the Sun as the primary, η ~ 0.1 translates into q ~ 0.5. It is also possible that our method fits noisy features of observed spectra with the model of the secondary star when the actual lines of the secondary star are comparable to the noise level or when the spectrum actually belongs to a single star. We avoided such unreliable solutions by the criterion on the luminosity ratio (η_n > 0.1), which has to be met in all of the four GALAH spectroscopic bands. In principle, the condition on the most extreme luminosity ratio is S/N dependent, meaning that in spectra with higher S/N, we can fit the lines of the secondary star more reliably, but we defer these considerations to future studies.

Criterion 4

We expect an approximate agreement between the velocity separation of double lines (ΔV_r) as determined by the binary pipeline and the one produced by the CCF detection method (see Sect. 5.2). Upon inspecting the difference between the two values, we decided to use the large discrepancy between the two (>20 km s⁻¹) as a criterion for rejection of poor solutions, which are mostly due to binary pipeline fit of the secondary component at arbitrarily large velocity separations in the case of noisy or otherwise peculiar spectra, where the binary signature is less evident.

6.2 Ensemble properties

As an additional confirmation of our results and to verify the adequacy of the selection criteria, we investigated ensemble distributions of model parameters θ derived by the binary pipeline and their interdependence for our final sample. First of all, Fig. 11 shows some expected correlations such as the connection between the temperature and surface gravity for individual stars (illustrating the main sequence) and the strong connection between temperature and microturbulence (see Eqs. (1) and (2) in Buder et al. 2018). The strong positive correlation between the temperatures of the primary and the secondary star stems from the constraint on the luminosity ratio, meaning that we cannot observe a hot primary in combination with a very cool secondary if we are to detect double lines in the spectrum. For the same reason, the radii of both stars in a binary system are also positively correlated.

We also note the negative correlation between the temperature of the two stars and [Fe/H]. Concerning this anti-correlation, we checked the GALAH DR2 sample of reliable results for single dwarf stars (following the dwarf-giant distinction given by Zwitter et al. 2018), and we found a similar trend, with a measured Pearson coefficient of − 0.33. This anti-correlation might be explained by the shortcomings of 1D hydrostatic models of stellar structure, which can overestimate the temperature of stellar atmospheres, especially in the lower metallicity regime (Bergemann et al. 2012 and references therein). Alternatively, this might be a feature of the underlying stellar population because when we consider equal-mass stars of different metallicities on the main-sequence, the more metal-poor stars have a smaller radius and higher surface temperature than their metal-rich counterparts.

7 Results

In this section we present the results of our analysis of binary stars. The catalogue contains the final sample as defined in Sect. 6, and its contents are described in Appendix A (Table A.1).

7.1 Kiel diagram of binary stars

Figure 12 shows the Kiel diagram for our final sample of binary stars. The positions of the primary stars agree well with the background distribution of the parameter values for single stars from the GALAH DR2 published results. As expected, the secondary stars are overall cooler and of higher surface gravity than the primary stars. However, the measured T_eff and log g values of secondary stars have larger uncertainties than those of primary stars (compare Fig. D.1), which causes the secondaries to be scattered more around their true position in the different parts of the Kiel diagram.

The poorly determined logg values of many secondary stars, whose radii indicate a main-sequence membership but which are instead positioned on the giant branch in Fig. 12, cannot be simply filtered out by applying cuts on the RV separation of double lines or the luminosity ratio. This issue might be alleviated if the spectra used here included more diagnostic lines for log g determination such as Mg b triplet or the Ca I 6162 Å lines (Ruchti et al. 2016, and references therein). It may also be that this is a chance alignment of a foreground dwarf (primary) and a background giant (secondary), in which case the binary pipeline would declare a small radius for the giant.

In general, we do not expect to find many double-lined binaries among giants for several reasons, for example, the constraint on similar luminosity and high enough orbital velocity. However, based on investigating all systems with a primary star situated above the dashed line in Fig. 12, which marks the distinction between dwarfs and giants as defined in Zwitter et al. (2018), we confirm their true binary nature in approximately half of the cases, whereas for others, the double lines in spectra are either not conclusive or the fit is clearly problematic. No general rule for distinguishing between the true giant binary systems and anomalous ones in this part of the Kiel diagram exists, therefore we advise caution in using such solutions.

7.2 Spatial distribution of binary stars

The volume of the Galaxy probed by the GALAH survey is described in Buder et al. (2018), while Figs. 13 and 14 place our sample of binary stars in the perspective of the survey footprint. The sky map shows a fairly homogeneous location of binary stars compared to all GALAH DR2 stars. The cumulative distributions of the distances in Fig. 14 demonstrate that the binary stars are observed to larger distances than their single-dwarf counterparts because of the magnitude limit of the survey and because binary stars are intrinsically brighter.

Fig. 11 Ensemble of derived parameters for the final sample of analysed binary systems. Distributions of θ parametervalues for 12 760 analysed objects are shown with diagonal panels, whereas the other panels indicate correlations between these model parameters. The associated uncertainty distributions are displayed in Fig. D.1. For clarity, we show the logarithm (base 10) of values for R1, R2, and E(B−V). The parameter units are given in Table 2.

7.3 Extinction

Following the same procedure as in Sect. 4, we compared E(B − V) values of the final sample of binary stars with the SFD and Bayestar reddening maps. Figure 15 shows a very good agreement with the Bayestar values for our sample, and the bias for the SFD values is systematic, which is expected for mostly nearby binary systems because SFD integrates extinction along the line of sight. We emphasise that in comparison to the results of the binary pipeline for single stars shown in Fig. 4, the reddening parameter E(B−V) appears to be better constrained in the case of binary stars because the photometric fit has less freedom for compensating the change in temperature by a corresponding change in E(B−V).

7.4 Population statistics

We first examine the metallicity distribution function (MDF) of our final sample of binary stars. Figure 16 shows the difference between the GALAH dwarf, giant, and binary stars. As expected, the giant stars have overall lower metallicities with a larger spread because they are located in a larger volume that reaches farther away from the Galactic plane than dwarf stars. The binary stars exhibit a narrow MDF reminiscent of the MDF for dwarf stars, but with lower metallicity values similar to those of giants. Several aspects need to be considered in trying to explain this, and we discuss them in Sect. 8.

The ΔV_r distribution in Fig. 17a shows a strong decline of binary stars with increasing RV separation. This is expected in general because for larger ΔV_r, the orbital period of the binary system must be short and the inclination high (edge-on view), and we have to catch the system at just the right time (close to quarter phase). Otherwise, these factors will greatly reduce the line-of-sight projected velocities of the two stars in the binary system and we will not see it.

The ratio of temperatures can indicate the evolutionary stage of stars in a binary system. Following our definition of the primary and secondary star, which states that the primary is the larger (Fig. 17c) and hence more massive star, and assuming that both are on the main-sequence, the condition T₁ ≥ T₂ is true. It can also be true for systems in which one or both stars are already evolved, but T₁ < T₂ indicates that at least the primary star has already evolved, and therefore the T₁ < T₂ condition sets the lower limit of binary stars with at least one component already evolved from the main-sequence. Figure 17b shows that there are ~25% of such binary stars in our final sample. When we assume that the majority of the remaining ~75% of systems (and both stars in them) lie on the main-sequence, and using standard relations (i.e. luminosity ∝ R²T⁴ ∝ M^3.8), we estimate that we probe systems of binary stars with a mass ratio down to 0.5.

The relatively large fraction (~25%) of binary stars with at least one evolved component can be partly attributed to the Malmquist bias because evolved stars are intrinsically brighter than their counterparts on the main-sequence and are therefore over-represented in a magnitude-limited sample of stars. On the other hand, it is also possible that the temperature estimate for the primary and secondary star has a different bias or systematics, which might result in a shift of the division that is indicated by Fig. 17b.

Combining the ratio of temperatures and radii in Fig. 17d, we note the high density of points that represents the main-sequence, while the largest relative difference in radii is observed for systems with already evolved stars. There is a clear indication that a diagonal line connects the highest ratio of radii and the highest ratio of temperatures, which suggests a forbidden region in the upper right part of the plot. This is due to the luminosity ratio restriction that we imposed in Sect. 6 (criterion 4).

Fig. 12 Kiel diagram of Teff and log g parameters as derived by the binary pipeline for the final sample of binary stars. The dashed line denotes a dwarf-giant distinction given by Zwitter et al. (2018), and the black edge of circles indicates systems where the primary lies above the line. The size of the circles relates to the radius of each star as derived by the binary pipeline, and the contours correspond to the density of the GALAH DR2 results for single stars (see Buder et al. 2018; we only used their results with flag_cannon = 0).

Fig. 13 Footprint of GALAH DR2 stars (black contours; single dwarf stars according to Zwitter et al. 2018), the final sample of binary stars (see Sect. 6), The Cannon validation sample, and the benchmark validation sample (see Sect. 4). GALAH stars include sources targeted by the GALAH pilot (Duong et al. 2018), K2-HERMES (Wittenmyer et al. 2018), and TESS-HERMES (Sharma et al. 2018) surveys, which undergo the same reduction pipeline as the main GALAH survey.

Fig. 14 Cumulative distribution of the distances towards sources based on Gaia DR2 parallaxes. The binary star sample is the sample defined in Sect. 6, while the GALAH DR2 single-star samples are those with flag_cannon = 0 (Buder et al. 2018) and the dwarf-giant distinction following Zwitter et al. (2018).

Fig. 15 Comparison of E(B−V) for the finalsample of binary stars derived by the binary pipeline to (a) Bayestar and (b) SFD reddening maps. The colour-coding is by absolute Galactic latitude.

Fig. 16 Cumulative distribution function of ([Fe/H]) for different samples of stars as indicated in the legend. The binary star sample is the one defined in Sect. 6, while the GALAH DR2 single-star samples are those with flag_cannon = 0 (Buder et al. 2018) and the dwarf-giant distinction following Zwitter et al. (2018).

7.5 Note on the accuracy of the results

All model parameters θ reported in the catalogue are accompanied by the corresponding uncertainties as given by our Bayesian scheme (see Sect. 3.2.4). The distributions of these formal uncertainties are displayed in Fig. D.1. However, as discussed in Sect. 3.2, our method ensures that these uncertainties are realistic only when all observational and model uncertainties are properly accounted for.

While it would be possible to account for the reliability or uncertainty of the spectroscopic model given its internal precision (see Figs. 2, B.1, and D.1) and the uncertainty of the parameters of the training set spectra, the evaluation of the reliability of the photometric model represents a greater challenge. One reason for this is the lack of flux-calibrated observed spectra that would span the FGK parameter space considered in this work and could be compared to model SEDs. Another obstacle is the assessment of validity of passband transmission curves because they are defined for different observing locations, different reference systems, with varying availability and quality of information in the literature (for more information on this, see Moro & Munari 2000; Fiorucci & Munari 2003). Additionally, the photometric observations, which cover a wide spectral range from visible to near-IR, suffer more from acomplex interplay of diverse extinction effects than the spectroscopic data in the form of normalised spectra.

The method we presented here incorporates only observational uncertainties as reported in the source catalogues, but it is general enough to allow inclusion of model uncertainties or systematic biases in observational data. We provide columns of additional information in our catalogue that indicate the performance of our method and quality of the results, and the user is encouraged to take these into account (see also Sect. 3.2.1). They are explained in Appendix A. Additionally, some indication of the performance of our method is given in the following sections.

7.5.1 Internal consistency check with repeat observations

In order to better evaluate the internal consistency of our method, we compared the derived binary parameters for systems that have repeat observations in GALAH. Distinct GALAH objects (as defined in Sect. 2) can represent repeat observations of the same star (or binary system), of which we have 341 binary systems with two repeats, followed by 36, 3, 1, and 2 systems with three, four, six, and seven repeats, respectively.

Although there are physical reasons for a possible change of parameters between repeat observations due to a different phase of the binary system (e.g. eclipses or reflection effects), and different results can be obtained due to a different quality of repeat observations (spectra), we expect that for the majority of cases, these differences are small.

We compute the maximum difference between repeat observations for each θ parameter and list the quantiles of the corresponding distributions in Table 6. This shows, for example, that for the majority of binary systems with repeat observations, the agreement of parameter values that are constrained by the spectroscopic fit is comparable to the indicated uncertainty of the spectroscopic model (see Fig. 2).

It is expected that the RV of the secondary component will change between repeats as a result of a change in binary phase, whereas the RV of the primary is expected to remain constant because spectra are predominantly shifted to the rest frame of the strongest spectral component in the main GALAH pipeline. Both these assumptions are confirmed by inspecting Table 6. However, exceptions are possible in the case of equally massive components.

7.5.2 Testing the co-evolution (same-age) assumption

We performed a cursory examination of the ages for the two components of our binary systems. The age estimates discussed here are derived following a Bayesian fitting of PARSEC isochrones introduced in Howes et al. (2019), which has also recently been used by Sahlholdt et al. (2019) to derive benchmark ages of Gaia benchmark stars. The input to the method are θ parameters T_eff[1,2], [Fe∕H], and bolometric luminosities of the two components. The luminosities are computed using T_eff[1,2] and R_[1,2], while the uncertainties for temperature and radii are compiled by summing the formal uncertainties and the median value from Table 6 in quadrature. The uncertainties of the luminosities are propagated in a Monte Carlo fashion given the uncertainties for temperature and radii. For each of the two stars in our final sample of binary systems, the method produces a PDF of the age estimate, and we extracted only the mode and the confidence limits of the PDF (see section 3.6 of Jørgensen & Lindegren 2005, confidence limits correspond to 1 σ if the PDF is Gaussian).

We used the average of the confidence limits to compute the relative uncertainty on the age estimate, and retained only those with a value below 10%. As expected, these are mainly stars from the turn-off and sub-giant branch. The comparison of primary and secondary ages for this selection of binary systems is shown in Fig. 18.

The ages of the primary and secondary star mostly agree given the individual age uncertainties, which is an encouraging confirmation of the validity of our results. However, we observe some outliers and a trend towards secondary stars being older than the primaries. This asymmetry can be explained by the potential systematic uncertainty in luminosity of the two components (i.e. one not reflected in the age uncertainties) and by considering how ages are computed. When we take two stars in a binary system whose true position is on the main-sequence, they will both move into the turn-off region when we increase their luminosity. Their ages can be reliably determined there, and those results will appear in Fig. 18. However, in this case, the secondary will always fall on an isochrone of greater age than the primary. On the other hand, when we decrease the luminosity of both stars, they will move to a region of densely packed isochrones with very poor age determination and are then excluded from the selection shown in Fig. 18.

In our modelled binary systems, the scatter and a potential bias in luminosity towards greater values can come from several sources: (a) uncertainty of parameters used for calculating luminosity, (b) uncertainty of the parallax for binary sources (see e.g. the recent study by Belokurov et al. 2020), (c) underestimated Gaia DR2 global parallax zero-point correction and possible unknown local variations, and (d) third-light contributions (e.g. foreground or background contamination, triple or higher order multiple systems). The latter is also supported by our close inspection of outliers in Fig. 18, where we identify a few resolved, nearly equal-mass triple systems, in which case the luminosity of an assumed binary system can be erroneously boosted by a factor of ~ 1.5. In this study, we apply a correction of +0.029 mas to the Gaia DR2 parallax values, which slightly decreases derived luminosities, however, some authors (see e.g. Fig. 14 in Chan & Bovy 2020) suggest larger corrections. Additionally, the complex morphology of isochrones can lead to several evolutionary stages being assigned to the same star, which can limit the accuracy of the age estimates presented here to ±20% (see e.g. discussions in Lebreton et al. 2014; Sahlholdt et al. 2019).

Fig. 17 Cumulative distribution functions of the RV separation (ΔVr) between the two stars in a binary system, ratio of temperatures (T1∕T2), ratio of radii (R1∕R2), and the relation between T1∕T2 and R1 ∕R2 in the fourth panel. The distribution of ΔVr is cut at the lower level of 15 km s−1 as explained in Sect. 6, while R1∕R2 starts at unity by construction because we define the primary as the star with a larger radius.

Fig. 18 Age comparison between the two components for binary systems with most reliable age estimates. The black line represents a one-to-one relation.

8 Discussion

8.1 Identification of binary stars

The FGK binary stars investigated here represent systems with relatively short orbital periods and other properties that enable the detection of multiple lines in spectra, for instance, similar luminosities. We applied conservative selection criteria to obtain our final sample of binary stars, but the true number of binaries in the data is much greater. Our detection methods might be further improved in order to extend the sample towards longer orbital periods and lower mass ratios (see also the methods in El-Badry et al. 2018 and Merle et al. 2020). Epoch photometry can be used to help identify binary star candidates in addition to spectroscopic identification. Eclipsing binaries, although much fewer in numbers and biased towards short-period systems, can be efficiently discovered by large-scale photometric surveys with high enough cadence. Facilities such as Gaia, LSST (Ivezić et al. 2019), and TESS (Ricker et al. 2015) promise to deliver orders-of-magnitude larger binary star samples than we currently have, but it will be a challenge to secure high-quality spectra for even a fraction of these detections in order to perform an analysis such as presented in this work.

Another detection approach is possible by considering high-quality astrometry, which enables the comparison of an astrometric parallax with the luminosity parallax. The latter can be inferred by knowing the absolute magnitude of the stellar system, which can be deduced based on the standard stellar parameters (T_eff, log g, [Fe/H]) measured from e.g. assumed single-star spectra. The binarity would be revealed in case of a significant mismatch between the two parallaxes, and this approach is discussed and attempted in Xiang et al. (2019). Similarly, Čotar et al. (2019) identified objects that are more luminous than predicted by their spectral type and Gaia distance, and specifically focused on double and triple Sun-like stars. Alternatively, comparing the results of spectral analysis that incorporate either a single or a binary star model can be used to distinguish between the two realities based on which model fits the data better, as demonstrated by El-Badry et al. (2018). The approaches described in this paragraph demand that detailed spectral analysis is performed on the whole data-set before attempting to identify binarity, however.

When we consider relatively wider binary systems, a detection method based on comparing 6D phase-space information combined with chemistry can be undertaken to infer binarity (see e.g. Simpson et al. 2019). However, when the angular separation between the two stars is sufficient, such a study would then transition to separate observations rather than observing light of the two stars in a single spectrum.

One potential contamination not considered here is chance alignments. Two stars can be positioned at angular distances smaller than the fibre aperture and thus be treated as a single target, although they can be spatially completely unrelated and thus can falsely pretend to represent a binary star. This might occur in GALAH, although efforts were made in the construction of the input catalogue (based on 2MASS) to avoid such contamination. A pilot study using a model of the Galaxy and observed fields of the GALAH survey suggests that these chance alignments are possible, especially towards lower Galactic latitudes, but more work is needed to properly account for them.

8.2 Undetected binarity and the effect on stellar parameters

In the context of Galactic archaeology, it is important to detect binary spectra, whether they result from chance alignments or actual gravitationally bound stars. Large-scale spectroscopic surveys should be aware of the influence of binarity on their analysis because some recent studies have shown that neglecting the proper treatment of binary spectra introduces systematic biases in the derived stellar parameters and elemental abundances that are not negligible.

Schlesinger et al. (2010) reported the uncertainty from binarity in SEGUE (Yanny et al. 2009) targets for spectra of S∕N ≥ 25 to be ~80 K in temperature and ~0.1 dex in [Fe/H], with a smaller impact on logg, while a study by El-Badry et al. (2017) estimated typical systematic errors of 300 K in temperature, 0.1 dex in [Fe/H], and 0.1 dex in logg for APOGEE-like spectra of solar-type stars. A proper investigation of this effect is therefore essential in order to make the findings of current and future projects dedicated to Galactic archaeology more reliable (see also Xiang et al. 2019).

Fig. 19 Mean metallicity of GALAH DR2 single stars binned by distance (top), distance from the Galactic plane (middle), and absolute Galactic latitude (bottom). The GALAH DR2 single-star samples are those with flag_cannon = 0 (Buder et al. 2018), and the dwarf-giant distinction following Zwitter et al. (2018). The sawtooth-shaped part of the lines indicates very few data points in the corresponding bins.

8.3 Systematically lower metallicity of the final sample

Figure 16 clearly shows an overall lower metallicity of the final sample of binary systems investigated in this work (median [Fe/H] ~0.2 dex lower compared to GALAH single dwarfs). To some extent, as indicated in Fig. 14, binary stars probe a more distant volume of space than single dwarf stars because their luminosity is intrinsically higher. Some change towards lower metallicities can be expected because of this: the trend of decreasing metallicity with distance (as well as with height above Galactic plane) is observed for GALAH stars (see Fig. 19). Furthermore, the distinction between GALAH single giant and dwarf stars adopted here (Zwitter et al. 2018), and the fact that there is a significant number of evolved systems in our binary star sample (see Figs. 12 and 17), make a comparison between single dwarf stars and binary stars more difficult. These binary systems with evolved and giant stars push the distance limit even further.

We plan to further investigate the comparison of metallicity values for single versus binary GALAH stars. The metallicity comparison can be performed per each small volume of physical space where both single and binary stars are located, while also examining the orbits of stars to determine their common origin. This analysis is defered to a future study because it will greatly benefit from the next Gaia data release (eDR3), which will bring more accurate 6D phase-space information.

Another aspect of lower metallicity values of our binary stars has to do with the fitting of the model to data for single as compared to binary stars. In general, there is a potential bias in the spectroscopically derived metallicity because the fit suffers to some extent from the degeneracy between temperature and metallicity, for example. However, in the case of binary stars, this degeneracy should be alleviated by the assumed common metallicity, which is now constrained by two sets of spectral lines. We expect this to be the case even more so when the two stars in a binary system have significantly different temperatures. On the other hand, when we consider binary systems where the mass ratio of the two stars or their age is such that they are found in significantly different evolutionary stages, we would have to consider effects of atomic diffusion and radiative levitation that can challenge the common metallicity assumption.

It is also possible that our binary identification techniques are unable to identify a significant portion of metal-rich binary spectra. In order to investigate this, we studied the shape of the CCF for a few dozen representative binary spectra that were created by summing GALAH DR2 single-star spectra. These binary spectra were split into three metallicity bins ([Fe/H] < −0.5, −0.5 < [Fe/H] < −0.1, [Fe/H] > −0.1) defined with the help of Fig. 20, and the two stellar components put at ΔV_r of 15 km s⁻¹ and higher, but no obvious diagnostic issue of the resulting CCFs was found for any of the metallicity/ΔV_r combinations.

Finally, the change to lower metallicities might be attributed to the widely observed negative correlation between metallicity and binary fraction for short-period binary systems. From the observational side, this has recently been addressed by Moe et al. (2019), who examined five different surveys and thoroughly accounted for their underlying selection biases. They measured the intrinsic occurrence rate of close solar-type binaries and showed that all five observational samples exhibit a quantitatively consistent anti-correlation between the close binary star fraction (a < 10 AU; this type of binaries also make up our final sample) and metallicity. Furthermore, hydrodynamical simulations performed by Bate (2019)agreed with this observed trend, with the explanation that the greater cooling rates at high gas densities that are due to the lower opacities at low metallicities increase the fragmentation on small spatial scales, leading to a higher ratio of close binary stars.

Fig. 20 Relationship between [Fe/H] and distance from the Galactic plane for our final binary star sample as defined in Sect. 6 and GALAH DR2 single-star sample with flag_cannon = 0 (Buder et al. 2018), where GALAH single stars are colour-coded by temperature and the binary stars are plotted in red.

8.4 Improvements of the model and the fitting procedure

The spectroscopic model for the data we presented here is based on a data-driven technique, The Cannon, which is trained on a training set analysed by a model-driven method, SME. We argue that this combination might be ideal because model-driven techniques are considered to be theoretically accurate although generally more expensive (and thus more appropriate for smaller data-sets), while data-driven ones are fast and can account for systematics in the data, but can be limited by their complexity and availability of good-quality and sufficient training data. Selecting good training data is a complex process and subject of ongoing investigation, where unsupervised classification methods and data visualisation (see Sect. 5.1) can help greatly in determining the types and number of distinct objects (spectra) in a given data-set.

The spectral flux modelling can also be improved with readily available techniques that would allow for a higher complexity or flexibilityof the spectroscopic model, such as the neural-networks-based The Payne (developed as a model emulator by Ting et al. 2019, see their Fig. 1) and its data-driven counterparts (see e.g. Xiang et al. 2019). Additionally, there are other complex physical aspects of flux modelling that might be crucial for our work, such as proper treatment of the limb darkening and rotational velocity that it depends on, the effects of which can for example change the fitted T_eff. Improvements in this direction can be achieved through higher complexity of the spectroscopic model, a better and larger training set, or even by post-processing of modelled spectra.

The photometric fit is strongly influenced by the Gaia photometry, but because of the large width of the three Gaia passbands, these are mostly important for constraining the total flux rather than the overall shape of the SED, and thus inclusion of more photometric data points (passbands) would benefit the accuracy of results. Furthermore, we here relied on the PySynPhot package to provide model-driven SEDs for our photometric model. However, in the future, a sufficient collection of observed SEDs might become available in order to train a data-driven model. Additionally, construction of SEDs provided by PySynPhot represents the majority of computational time in our method, thus it is possible that a method such as The Cannon or The Payne might perform better, having been trained either on model-driven or data-driven training data.

Passband transmission functions are another aspect of the photometric model where uncertainties are sometimes not known or difficult to estimate. The revised passbands for Gaia DR2 are an example of a step forward in this context (Weiler 2018). Moreover, synthetic photometry involving photometric systems that have been defined from the ground and at different locations, or photometry that is given in the Vega system, are often only accurate to about 5% (STScI Development Team 2013), which mostly pertains to the overall zero-point comparison between synthetic and observed magnitudes, rather than the shape of the SED. This might explain the discrepancy between the radii derived here for single stars and the seismic values (see Sect. 4).

From a purely technical point of view, the convergence of our method to the best possible solution, or in other words, to the global minimum, is not strictly guaranteed. We chose the scheme of MCMC sampling of the posterior distribution function in order to explore the whole model parameter space. However, we implement some custom functionalities to make it more time efficient, and thus it can happen that in some cases, the convergence is not optimal. Solutions potentially include the use of Hamiltonian Monte Carlo (HMC; Neal 2011), which effectively reduces the number of Markov chain samples required for the same result, and other methods from the statistical modelling and data analysis platform (STAN; Carpenter et al. 2017), such as variational inference for a much quicker first estimate of parameter distributions.

8.5 Comparison to previous works

This paper presents a complementary effort to previous studies of binary stars performed in the context of large all-sky surveys (Holmberg et al. 2009; Matijevič et al. 2010; Gao et al. 2014; Merle et al. 2017; Skinner et al. 2018a; El-Badry et al. 2018; Kounkel et al. 2019). The binary star samples investigated in those studies were connected to diverse observational strategies, instrumentation used by the corresponding surveys (e.g. the Geneva-Copenhagen Survey – GCS, Holmberg et al. 2009; RAVE, Steinmetz et al. 2006; Sloan Digital Sky Survey – SDSS, Blanton et al. 2017), the methods, and derived parameters, which makes them difficult to directly compare to our sample. Purely in terms of the number of identified and analysed binary systems, our sample of 12 760 SB2s is a significant increase, for example, from Matijevič et al. (2010, 2011), Merle et al. (2017), or El-Badry et al. (2018, ~ 3000 binary stars), mainly by virtue of the size and properties of the GALAH data-set. Additionally, our sample includes several dozen SB2 giant systems, which otherwise are a rare occurrence, with only a handful known in the literature (see e.g. Merle et al. 2017).

The philosophy and approach of combining all available observational data in order to extract the maximum amount of information from studied objects have previously been proposed and developed for example by Schönrich & Bergemann (2014), and recentlyby Cargile et al. (2019). In contrast to these studies of single stars, we avoid using models of stellar evolution (e.g. isochrone fitting) in order to obtain the properties of both stars in a binary system with the least amount of model-dependent assumptionsas possible. However, the framework that we develop can be easily extended to include more theoretical background if needed. We also recognise the work by Stassun & Torres (2016) on benchmarking astrometric parallaxes as a parallel effort to combine different observational data in the context of binary stars, albeit with a different prior or model configuration and desired outcomes.

8.6 Gaia revolution

In the context of statistics of binary stars, it is important to mention Gaia DR3, which is scheduled for the second half of 2021, in which solutions for many astrometric and spectroscopic binaries will be published, along with about one million eclipsing binaries that are expected to have been identified during processing of Gaia DR2. Gaia data in synergy with other data-sets of binary stars, such as the one presented in this work, have the potential to truly revolutionise the field of multiplicity statistics across the full range ofprimary masses, mass ratios, orbital periods, eccentricities, chemistry, ages and evolutionary states, and environments.

9 Conclusions

We presented a general comprehensive method for the analysis of double-lined FGK binary stars, validated and applied on a large sample of identified SB2 systems. We identified binary stars in the GALAH spectroscopic sample using two different approaches: a sophisticated method based on the cross-correlation function, and a classification approach based on the dimensionality reduction technique t-SNE. By combining candidates from these two identifications, we obtained a sample of 19773 candidate SB2s out of 587 153 spectra from the second GALAH internal data release.

To characterise the binary candidates, we combined photometric, spectroscopic, and astrometric data in a Bayesian framework and used an MCMC method to obtain the best-fit model parameters. We validated our method on a benchmark sample of 620 single stars to confirm that our derived parameters are consistent with those from external sources, such as the effective temperature from the IRFM method or stellar radius from asteroseismology. Following the analysis of candidate binary systems and an application of conservative selection criteria, we present a catalogue of 12 760 reliable detections. We derived standard stellar parameters for the two stars in a binary system, their radii, and interstellar reddening.

Based on the parameters derived here, we conclude that our catalogue represents a population of close binary systems (a < 10 AU) with mass ratios 0.5 ≤ q ≤ 1 that can be detected through double lines in their spectra. The results additionally indicate an overall lower metallicity of binary stars compared to single dwarf stars that are observed in the same magnitude-limited sample of the GALAH survey. Some of the parameters are used to construct cumulative distributions shown in Fig. 16, which can provide further insight into the statistical properties of the analysed population of binary stars.

We aim to improve and employ our method for the future GALAH releases, and to apply it to forthcoming stellar spectroscopic surveys such as WEAVE (Dalton et al. 2012) and 4MOST (de Jong et al. 2012).

Acknowledgements

We acknowledge helpful conversations with Paul McMillan and Thomas Bensby. We are grateful to the referee for the very useful comments which also encouraged us to investigate certain aspects of our analysis that improved the understanding and validity of our results. We acknowledge the use of data reduction infrastructure developed within the GALAH Survey collaboration. The GALAH survey web site is www.galah-survey.org. This work is based on data acquired through the Australian Astronomical Observatory, under programmes: A/2014A/25, A/2015A/19, A2017A/18 (The GALAH survey); A/2015A/03, A/2015B/19, A/2016A/22, A/2016B/12, A/2017A/14 (The K2-HERMES K2-follow-up program); A/2016B/10 (The HERMES-TESS program); A/2015B/01 (Accurate physical parameters of Kepler K2 planet search targets); S/2015A/012 (Planets in clusters with K2). We acknowledge the traditional owners of the land on which the AAT stands, the Gamilaraay people, and pay our respects to elders past and present. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www. cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/ dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research was made possible through the use of the AAVSO Photometric All-Sky Survey (APASS), funded by the Robert Martin Ayers Sciences Fund and NSF AST-1412587. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. This research has made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This work was supported by the Swedish strategic research programme eSSENCE. G.T. and S.F. were supported by the project grant “The New Milky Way” from the Knut and Alice Wallenberg foundation and by the grant 2016-03412 from the Swedish Research Council. T.M. is supported by a grant from the Fondation ULB. Y.S.T. is grateful to be supported by the NASA Hubble Fellowship grant HST-HF2-51425.001 awarded by the Space Telescope Science Institute. G.T., K.Č., and T.Z. acknowledge financial support of the Slovenian Research Agency (research core funding No. P1-0188 and project N1-0040). The computations were partially performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at LUNARC partially funded by the Swedish Research Council through grant agreement no. 2016-07213 and partially funded by the Royal Fysiographic Society of Lund. JDS acknowledges the support of the Australian Research Council through Discovery Project grant DP180101791. SB acknowledges funds from the Alexander von Humboldt Foundation in the framework of the Sofja Kovalevskaja Award endowed by the Federal Ministry of Education and Research. This research has been supported by the Australian Research Council (grants DP150100250 and DP160103747). This research was partly supported by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. This research has made use of NASA’s Astrophysics Data System and the CDS services (Strasbourg, France).

Appendix A The catalogue

The binary star parameters derived in this work are presented in the catalogue, which is described in Table A.1. In order to get the reduced (perdata point) χ² values of the spectroscopic and photometric fit separately, the columns chi2_binary_pipeline_spec and chi2_binary_pipeline_phot should be divided by the fixed number of spectroscopic data points (N_spec = 12 230) and the number of available photometric data points (obtained by inspecting NULL values of photometric data in Table A.1), respectively. These can then be used as quality indicators or the formal goodness of fit. If all the uncertainties of our analysis scheme were properly accounted for, these reduced χ² values would be close to 1.

These values, however, can be large for several reasons such as strong oversubtraction of the background light in the spectral reduction, unrecognised high values of noise in spectra, improper continuum normalisation and discordant wavelength solution in spectra, erroneous photometric data, underestimated observational uncertainties, issues with the spectroscopic or photometric model, and failed convergence of the solution. Conversely, the reduced χ² values can also be too low, for example, because the observational uncertainties are overestimated. For example, the median of ${\chi }_{\text{spec}}^2$ values for our results is ~0.8, which implies that our spectroscopic observational uncertainties, as derived in Kos et al. (2017), might be slightly overestimated. One reason for this might also be correlations between different contributions to the observational uncertainty budget.

It is beyond scope of this work to determine which of these issues influence our reduced χ² values in a specific case, but we aim to conduct this proper analysis in the future. The user of the data presented in Table A.1 should be aware of the possible underlying factors.

The reported radius of the stars in binary systems is directly related to the value of parallax, therefore we encourage users to investigate the quality of the parallax measurement when the radii values derived here are used. This can be aided by the RUWE parameter from the Gaia DR2 catalogue (included in Table A.1).

Appendix B Higher order of The Cannon

Figure B.1 shows a comparison of SME derived labels and the labels reproduced by The Cannon interpolation for the training set, and demonstrates that augmenting the labels to the third order improves The Cannon performance. This trend continues with the fourth order as expected because of the increased flexibility (complexity) of the model, but a more thorough investigation of The Cannon of higher orders is required to determine the limits of over-fitting, for instance, or issues with extrapolation in a model that is too complex.

We onlyshow plots of the three basic stellar parameters in Fig. B.1, but six parameters were used in generating The Cannon model (T_eff, log g, [Fe/H], v_mic, v_broad, and A_Ks), in accordance with the main GALAH pipeline. The values of the bias and standard deviation in Fig. B.1 (left column) arean indication of the internal precision of The Cannon model we used, and they agree with those reported in Buder et al. (2018; see their Fig. 8) for all six parameters. We note that these values should always be appropriately incorporated in the presentation of final results obtained by inference that includes The Cannon method.

Fig. B.1 Comparison of SME and The Cannon label values for the training set (single-star spectra). Panels a–c: similar to Fig. 8 in Buder et al. (2018), but we also include a comparison to the model generated with labels in the third order shown in panels d–f.

Appendix C t-SNE map

Figure C.1 shows the t-SNE map (presented in Buder et al. 2018) that was used to identify binary star candidates in addition to the CCF method. Only a few dozen triple stars are flagged, and they are hardly visible next to the larger group of binary starsto the left, whereas additional smaller groups of binaries are located in the top left corner in the figure.

Fig. C.1 t-SNE projection map of 587 153 GALAH spectra. The points (spectra) are colour-coded by classification category. The count of spectra for each category is given in the legend. Black points denote normal, well-behaved spectra.

Appendix D Formal uncertainties of the binary pipeline

When the results from our catalogue of binary star solutions are used, formal uncertainties (distributions shown in in Fig. D.1) should be taken in consideration together with validation tests in Sects. 4 and 7.5.1.

Fig. D.1 Ensemble of uncertainties of θ parameters for the final sample of 12 760 analysed binary systems (see also Fig. 11). The uncertainty distributions are shown with diagonal panels, and the other panels indicate their correlations. The uncertainty units are those of the corresponding parameters (see Table 2). The tails of the distributions are truncated for clarity.

Appendix E Photometric validation tests

Figure E.1 shows results of the photometric validation for the benchmark validation sample of stars. The grid effects are especially clearly visible in the temperature comparison, where Gaia magnitudes are included in the fit.

Fig. E.1 Comparison of Teff and E(B−V) derived here to the literature values (see text) for the benchmark validation sample. Panels e–h: results for the special case where we exclude Gaia magnitudes from the fit. The colour-coding is by Teff (panels c, d, g, and h) and the formal uncertainty on Teff (panels a, b,e, and f). The label “GALAH” stands for GALAH DR2 published values.

References

All Tables

All Figures

Fig. 1 Spectral energy distributions for an example of a binary system (Teff,1 = 5018 K, Teff,2 = 5715 K, R1 = 1.9 R⊙, R2 = 1.3 R⊙, E(B− V) = 0.13). The legend explains which SEDs are plotted (P is the primary star, and S is the secondary star). The four GALAH spectral bands are indicated by vertical shaded regions, while the width of the transmission functions of photometric bands is indicated on the top. The black points are observed integrated fluxes in the following order of photometric bands: BJ, GBP, VJ, G, GRP, J, H, KS, W1, and W2 (see Sect. 3.2.2), and the short green lines indicate the synthetic integrated fluxes at the same position. The error bars of observed fluxes are exaggerated by a factor of 10 for clarity. Despite this, the uncertainty for the Gaia G value is smaller than the point size.

In the text

	Fig. 2 Comparison of SME and The Cannon labels for The Cannon validation sample. The Cannon values produced by the binary pipeline are subtracted from the SME values.
In the text

Fig. 3 Comparison of Teff (orange) and E(B−V) (green) parameters as derived by the binary pipeline (using only the photometric fit) to the literature values (indicated on the horizontal axis). The literature values are subtracted from those given by the binary pipeline. The lighter shaded distributions correspond to results from the binary pipeline when the Gaia magnitudes are excluded from the fit.

In the text

	Fig. 4 Temperature and reddening for the benchmark validation sample derived by the binary pipeline compared to the literature values from the IRFM method and Bayestar and SFD reddening maps (see also Fig. 3). The colour-coding is by absolute Galactic latitude.
In the text

Fig. 5 Comparison of radii for the benchmark validation sample derived by the binary pipeline, Gaia DR2, and computed from asteroseismic parameters (more details in Sect. 4.3). The colour-coding is by relative parallax uncertainty except in panel b, which is a copy of panel a, but colour-coded to indicate stars with less reliable seismic information (νmax < 40; S. Sharma, priv. comm.). The more reliable part of the sample is used to calculate the bias and scatter as reported in panels a, c, and d. The colour scale in panel a is clipped for clarity, while panels b and c have a smaller range of σϖ ∕ϖ because some Gaia DR2 values for radii were lacking.

In the text

	Fig. 6 Effect of excluding Hα (in band 3) and Hβ (in band 1) in the solar template on the width of the CCF. The decrease in CCF peak width reveals a second close component that is difficult to detect when Hα and Hβ are includedin the template.
In the text

Fig. 7 Identification of the RV components in a CCF of a close SB2 candidate. The top panel shows the CCF. The three panels below correspond to the three successive derivatives. The vertical black lines show the RV components detected by taking the ascending zeros of the third derivative (bottom panel). The vertical red lines give the positions of RV components obtained by fitting the multi-peaked CCF.

In the text

	Fig. 8 t-SNEprojection map of 587 153 spectra as presented in Buder et al. (2018). The points (spectra) are colour-coded by all t-SNE detected SB2s (blue) and the unique additional detections by CCF (red). The count of spectra for each group is given in the legend.
In the text

	Fig. 9 Colour-magnitude diagram of SB2 identifications. (a) All identifications with CCF, (b) all identifications with t-SNE. The grey points are all GALAH objects considered in this work. The absolute magnitudes and colours are taken from the Gaia DR2 catalogue, and no reddening corrections have been applied.
In the text

Fig. 10 Example fits of binary stars from the final sample. Top and bottom panels: binary system at low- and high-velocity separation of double lines, respectively. The spectral region around Hβ and Hα is shown in the left and middle, respectively, while a region of a few well-isolated lines is shown on the right. Some of the derived θ parameters are indicated for each binary system, and the bottom example shows that the primary star has already evolved from the main-sequence. The primary and secondary component of the modelled binary spectrum are shown separately at a continuum level in accordance with their contribution. The black error bars in the inset plot of the photometricfit represent the relative uncertainty of the observed fluxes in the indicated photometric passbands, while the red points represent the relative offset of synthetic fluxes. The error bars of Gaia photometry are multiplied by a factor of 10 for clarity.

In the text

Fig. 11 Ensemble of derived parameters for the final sample of analysed binary systems. Distributions of θ parametervalues for 12 760 analysed objects are shown with diagonal panels, whereas the other panels indicate correlations between these model parameters. The associated uncertainty distributions are displayed in Fig. D.1. For clarity, we show the logarithm (base 10) of values for R1, R2, and E(B−V). The parameter units are given in Table 2.

In the text

Fig. 12 Kiel diagram of Teff and log g parameters as derived by the binary pipeline for the final sample of binary stars. The dashed line denotes a dwarf-giant distinction given by Zwitter et al. (2018), and the black edge of circles indicates systems where the primary lies above the line. The size of the circles relates to the radius of each star as derived by the binary pipeline, and the contours correspond to the density of the GALAH DR2 results for single stars (see Buder et al. 2018; we only used their results with flag_cannon = 0).

In the text

Fig. 13 Footprint of GALAH DR2 stars (black contours; single dwarf stars according to Zwitter et al. 2018), the final sample of binary stars (see Sect. 6), The Cannon validation sample, and the benchmark validation sample (see Sect. 4). GALAH stars include sources targeted by the GALAH pilot (Duong et al. 2018), K2-HERMES (Wittenmyer et al. 2018), and TESS-HERMES (Sharma et al. 2018) surveys, which undergo the same reduction pipeline as the main GALAH survey.

In the text

	Fig. 14 Cumulative distribution of the distances towards sources based on Gaia DR2 parallaxes. The binary star sample is the sample defined in Sect. 6, while the GALAH DR2 single-star samples are those with flag_cannon = 0 (Buder et al. 2018) and the dwarf-giant distinction following Zwitter et al. (2018).
In the text

	Fig. 15 Comparison of E(B−V) for the finalsample of binary stars derived by the binary pipeline to (a) Bayestar and (b) SFD reddening maps. The colour-coding is by absolute Galactic latitude.
In the text

	Fig. 16 Cumulative distribution function of ([Fe/H]) for different samples of stars as indicated in the legend. The binary star sample is the one defined in Sect. 6, while the GALAH DR2 single-star samples are those with flag_cannon = 0 (Buder et al. 2018) and the dwarf-giant distinction following Zwitter et al. (2018).
In the text

Fig. 17 Cumulative distribution functions of the RV separation (ΔVr) between the two stars in a binary system, ratio of temperatures (T1∕T2), ratio of radii (R1∕R2), and the relation between T1∕T2 and R1 ∕R2 in the fourth panel. The distribution of ΔVr is cut at the lower level of 15 km s−1 as explained in Sect. 6, while R1∕R2 starts at unity by construction because we define the primary as the star with a larger radius.

In the text

	Fig. 18 Age comparison between the two components for binary systems with most reliable age estimates. The black line represents a one-to-one relation.
In the text

Fig. 19 Mean metallicity of GALAH DR2 single stars binned by distance (top), distance from the Galactic plane (middle), and absolute Galactic latitude (bottom). The GALAH DR2 single-star samples are those with flag_cannon = 0 (Buder et al. 2018), and the dwarf-giant distinction following Zwitter et al. (2018). The sawtooth-shaped part of the lines indicates very few data points in the corresponding bins.

In the text

	Fig. 20 Relationship between [Fe/H] and distance from the Galactic plane for our final binary star sample as defined in Sect. 6 and GALAH DR2 single-star sample with flag_cannon = 0 (Buder et al. 2018), where GALAH single stars are colour-coded by temperature and the binary stars are plotted in red.
In the text

	Fig. B.1 Comparison of SME and The Cannon label values for the training set (single-star spectra). Panels a–c: similar to Fig. 8 in Buder et al. (2018), but we also include a comparison to the model generated with labels in the third order shown in panels d–f.
In the text

	Fig. C.1 t-SNE projection map of 587 153 GALAH spectra. The points (spectra) are colour-coded by classification category. The count of spectra for each category is given in the legend. Black points denote normal, well-behaved spectra.
In the text

	Fig. D.1 Ensemble of uncertainties of θ parameters for the final sample of 12 760 analysed binary systems (see also Fig. 11). The uncertainty distributions are shown with diagonal panels, and the other panels indicate their correlations. The uncertainty units are those of the corresponding parameters (see Table 2). The tails of the distributions are truncated for clarity.
In the text

	Fig. E.1 Comparison of Teff and E(B−V) derived here to the literature values (see text) for the benchmark validation sample. Panels e–h: results for the special case where we exclude Gaia magnitudes from the fit. The colour-coding is by Teff (panels c, d, g, and h) and the formal uncertainty on Teff (panels a, b,e, and f). The label “GALAH” stands for GALAH DR2 published values.
In the text