Inferring binary and trinary stellar populations in photometric and astrometric surveys
Axel Widmark Boris Leistedt David W. Hogg
astro-ph.SR astro-ph.GA



















































Axel Widmark
axel.widmark@fysik.su.se

0000-0001-5686-3743]Axel Widmark
Department of Physics, AlbaNova, Stockholm University, SE-106 91, Stockholm, Sweden
The Oskar Klein Centre for Cosmoparticle Physics, AlbaNova, SE-106 91 Stockholm, Sweden

0000-0002-3962-9274]Boris Leistedt
Center for Cosmology and Particle Physics, Department of Physics,  
   New York University, 726 Broadway, New York, NY 10003, USA
NASA Einstein Fellow

0000-0003-2866-9403]David W. Hogg
Center for Cosmology and Particle Physics, Department of Physics,  
   New York University, 726 Broadway, New York, NY 10003, USA
Center for Data Science, New York University, 60 Fifth Avenue, New York, NY 10011, USA
Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA
  
















Multiple stellar systems are ubiquitous in the Milky Way, but are often unresolved and seen as single objects in spectroscopic, photometric, and astrometric surveys.
Yet, modeling them is essential for developing a full understanding of large surveys such as Gaia, and connecting them to stellar and Galactic models. In this paper we address this problem by jointly fitting the Gaia and 2MASS photometric and astrometric data using a data-driven Bayesian hierarchical model that includes populations of binary and trinary systems. This allows us to classify observations into singles, binaries, and trinaries, in a robust and efficient manner, without resorting to external models. We are able to identify multiple systems and, in some cases, make strong predictions for the properties of its unresolved stars. We will be able to compare such predictions with Gaia Data Release 4, which will contain astrometric identification and analysis of binary systems.























§ INTRODUCTION

A large fraction of stars, possibly a majority, exist in binaries and higher multiple stellar systems <cit.>. In stellar surveys, such systems are resolved only if they are sufficiently close. When unresolved, they are reported as single objects where the respective stellar fluxes are added into a total magnitude. Even though a multiple stellar system cannot be resolved by pure angular separation, there are other techniques to do so. With multi-epoch photometry it can be possible to observe eclipses, depending on the spatial orientation of the system. For very tight binaries with high rotational periods, spectroscopy with high time resolution can be used to infer the rotation from a periodic widening/narrowing of stellar emission lines, as applied to white dwarf binaries in the Supernova Progenitor Survey <cit.>. Also, without such time variations, clues are given by spectral or photometric information. For example, the sum of two black body spectra of different temperatures is not itself a black body spectrum. With good stellar models it is possible to disentangle such systems, as discussed in <cit.>. In this article, we use a data-driven model to infer multiple stellar systems, meaning that we do not rely on stellar models. The strength of this approach is its generality, in that we fit our model to the data with minimal a priori assumptions about the stellar population. In the era of large stellar surveys such as Gaia <cit.>, it should be possible to simultaneously model and fit for single and multiple systems without the need of external models, as we demonstrate in the present work. While other identification methods are observationally expensive and often require special purpose surveys, this method relies on readily available photometry.

In this paper we focus on the identifying multiple stellar systems in the first data release of the Gaia survey <cit.>. The aim is to infer what stellar objects in the Gaia catalogue are binary or trinary from two-band photometry and an astrometric distance measurement. Combining two stars from the main sequence into an unresolved binary does not result in a stellar system that is itself consistent with the main sequence distribution, at least not where the main sequence is sufficiently narrow in color-magnitude space. For example, combining two identical stars will double the flux, giving identical color but a magnitude difference of Δ M ≃ 0.75.

We use a simple hierarchical Bayesian model in a small window of color-magnitude space, with cuts on observed parallax (giving a volume limited sample). We assume a main sequence of single stellar systems that is Gaussian in width and infer a population of binary and trinary stellar systems. This is a powerful way to identify multiple stellar systems and in some cases even infer the color and magnitude of the individual stars, which in turn could be turned into constraints on the individual masses and temperatures using stellar models.

In the longer term, we hope to include this formalism in a larger framework of data-driven models, as used in <cit.>. In both these sources, the authors use Gaussian Mixture Models. With this long term goal in mind, we describe the population of stellar systems in color-magnitude space in terms of a mixture of multivariate Gaussians.

By the end of its mission, Gaia is predicted to have photometry errors in the order of milli-magnitudes. However, for now the photometry is nowhere near that level of precision, and does not provide color information, which is why we in this article use photometric information from the 2MASS survey.

This article is outlined as follows. In section <ref>, we describe the catalogue and the sample of data we use. In section <ref>, we present the statistical model for the stellar population. In section <ref>, we show how to infer properties of component stars in a multiple stellar system. In section <ref>, we present our results. In section <ref>, we discuss and conclude.







§ TGAS/2MASS CROSS MATCH


In this work we use a catalogue of astrometric data from the Tycho-2 Gaia Astrometric Solution <cit.> cross matched with photometric data from the Two Micron All Sky Survey (2MASS). The cross match is taken from <cit.>.




 §.§ Observables and errors


The observables we take into account are parallax, ϖ, and apparent magnitudes in J and K bands, m_J and m_K. A hatted quantity (ϖ̂) refers to the observed value, while a non-hatted quantity (ϖ) refers to its true value. We operate in the space of color and absolute magnitude in the J-band, where the true and observed quantities are


    c     = m_J-m_K, 
    ĉ    = m̂_J-m̂_K, 
    
        M     = m_J - 5[log_10(arcs/ϖ)-1], 
    M̂    = m̂_J - 5[log_10(arcs/ϖ̂)-1].

In this article, the only absolute magnitude considered is that of the J-band, so for the sake of brevity we write the absolute magnitude without the J-index.

We propagate the error on the parallax to first order, giving a error covariance matrix over c and M


    Σ_c,M^(i) = 
        [               σ_J^2+σ_K^2                                               σ_J^2;                     σ_J^2                           σ_J^2 + (5σ_ϖ/ϖ̂log 10)^2 ],

with an implicit index i that denotes the stellar object. This first order propagation is motivated by a high parallax significance, ϖ̂/σ_ϖ. Objects in our sample have a median parallax significance of 25 and a minimum value of 5. See section <ref> below for details on how the sample is constructed.

Because the population is modeled as a Gaussian mixture, it is very convenient to have errors also in the form of a covariance matrix. Computing the posterior becomes more efficient, since the convolution of two bivariate Gaussians has a simple analytical form.

We do not consider dust extinction corrections to the photometric magnitudes, as these effects are expected to be negligible in the color–magnitude box we consider given our volume limited sample. To check this, we calculated the dust corrections using the <cit.> three-dimensional dust map, evaluated at the point distance estimate 1/ϖ̂ for each object. We find that the median dust corrections to the J and K band magnitudes are 0.0068 and 0.0029, respectively, while the 90th percentiles are 0.0412 and 0.0174. These corrections are small compared to the width of the components of our Gaussian mixtures, and we neglect them in order to simplify our analysis.



 §.§ Data cuts


We make cuts to the data according to the following criteria:


	
  * an observed parallax ϖ̂>5 mas (corresponding to a distance <200 pc),
    
  * an observed color ĉ in interval 0.5–0.8,
    
  * an observed absolute magnitude M̂ in interval 3–6,
    
  * no NaN-valued observables,
    
  * no excessive photometric noise, defined as σ_J^2+σ_K^2<0.1^2 (this criterion removes 11 objects from our sample).


This gives a total of 23,112 individual objects. The full color-magnitude diagram of the TGAS/2MASS cross-match is shown in figure <ref>, also marking the smaller window of parameter space that we operate in. These cuts on color and magnitude limits our single star population to K-type dwarfs <cit.>.







§ MODEL


In this section we outline the hierarchical model that includes binaries and trinaries in our window of color-magnitude space. The distribution of single stars is modeled by Gaussians generated from population parameters α, β_1, β_2, and γ. These Gaussians are centered on a straight line, according to


    M = β_1 + β_2 (c-0.5),

with a spacing of 0.05 in color between Gaussians. The covariance matrix describing the Gaussians have one eigenvector pointing along the line described by equation (<ref>), which is set to width 0.05 in c, and one eigenvector perpendicular to said line of width γ.

The amplitude, i.e. the number of stars in each Gaussian, is proportional to


    n_singles(c) ∝ e^-α c.


In other words, the parameter α describes how the number density changes with color; β_1 and β_2 describe the slope of the distribution; γ describes the thickness of the main sequence distribution of single stars. The Gaussian mixture that describe the population of single stars has parameters


    c̅_s    = 0.05s 
    M̅_s    = β_1 + β_2(0.05s-0.5) 
    Σ_(s)   = [   0.05^2+β_2^2(β_2^2+1)^-1γ^2                                0.05^2β_2+β_2(β_2^2+1)^-1γ^2;  0.05^2β_2+β_2(β_2^2+1)^-1γ^2                               0.05^2β_2^2 + (β_2^2+1)^-1γ^2 ],

where s is the index of the Gaussian and the bar in c̅ and M̅ indicate that they are mean values of the respective Gaussians. The integer on the covariance matrix Σ_(s) is in parenthesis, as it is identical for all Gaussians of the singles population. The two eigenvectors are manifest as the two separate terms of the respective matrix element. Because we limit our sample to an interval in color, c∈[0.5,0.8], the Gaussian mixture of the singles population includes Gaussian with a mean color in region c̅∈[0.4,0.9], giving a total of 11 Gaussians.
Despite the simplicity of this model, we expect from stellar models and existing data that it will provide an accurate description of the population of single stars.

In addition to aforementioned four population parameters, there are three more population parameters that describe the binary and trinary populations: η, which parameterizes the pairing mechanism of multiple stellar systems (how favored are systems where the component stars have similar properties); and the fraction of binary and trinary systems, f_b and f_t, corresponding to the (relative) prior probabilities of an object being an unresolved binary or trinary system as opposed to a single star.



 §.§ Gaussian mixture and binary population


We modeled the population of single star objects as a sum of Gaussians in color-magnitude space. The population of binary and trinary stellar systems are sums of the single population Gaussians. Adding up the magnitude of the objects making an unresolved multiple system follows the formula


    M_sum = -5/2log_10( ∑_i 10^ -2/5M_i ).

Because this is not a straight forward addition, the sum of two Gaussian distributions in color-magnitude space is not exactly a Gaussian itself. However, the actual result is close to a Gaussian distribution and can be approximated as such to an adequate degree of accuracy. Given two Gaussians with means {c̅_i,M̅_i} and covariance matrices Σ_i, where i={1,2}, a first order approximation to the sum is a Gaussian distribution, where the means are added together according to equation (<ref>), and the summed covariance matrix is given by


    Σ_sum = 𝒥_1 Σ_1 𝒥_1^⊤ +  𝒥_2 Σ_2 𝒥_2^⊤.

Here, 𝒥_i={1,2} are Jacobians of the form


    𝒥_i     = [ ∂ c_sum/∂ c_i|_c_i=c̅_i                         ∂ c_sum/∂ M_i|_M_i=M̅_i; ∂ M_sum/∂ c_i|_c_i=c̅_i                         ∂ M_sum/∂ M_i|_M_i=M̅_i ] = 
        = [                             g(M̅_i-c̅_i)/ g(M̅_1-c̅_1)+g(M̅_2-c̅_2)                                                                     g(M̅_i)/ g(M̅_1)+g(M̅_2) -  g(M̅_i-c̅_i)/ g(M̅_1-c̅_1)+g(M̅_2-c̅_2);                                                                   0                                                                                                                g(M̅_i)/ g(M̅_1)+g(M̅_2) ],

written in terms of the function g(x)=10^-2x/5.

While this is only propagating the Gaussians to first order, it does work very well, as is evident from figure <ref>. The approximation works well for adding Gaussians that are sufficiently small, such that higher order terms does not come into significant effect.




For the Gaussian mixtures of the binary and trinary populations, the amplitude of the Gaussian is proportional to exp(-αc̅)exp(ηΔ_c̅), where Δ_c̅=|c̅_A-c̅_B| for binaries and Δ_c̅=|c̅_A-c̅_C| for trinaries, which is the difference in color between the highest and lowest Gaussian mean color values. The first factor of the amplitude is directly proportional to the number density of single stars of the same color. The seconds factor comes from the fact that multiple stellar systems with similar properties might be favored or disfavored, which we parameterize with population parameter η. We would expect that binaries with stars of approximately equal mass, and similar color, are more likely than binaries with a significant mass difference <cit.>.

As in the case of the Gaussian mixture that describes the singles population, we constrain ourselves to include Gaussians that have a mean color in the interval c̅∈[0.4,0.9]. A Gaussian in the mixture that describes the binary or trinary population is built from two or three parent Gaussians belonging to the singles population. The parent Gaussians do not necessarily have a mean color within the interval c̅∈[0.4,0.9], even though the composite stellar object does so. For example, it is possible to have a very dim component star, of large color value c, such that its contribution to the total color is small. For this reason we assume the distribution of singles to be extended analytically beyond this window according to the model. As the hydrogen burning limit is reached around c=1.1 <cit.>, the outermost single star parent Gaussian has mean color c̅=1.05. This makes the distribution of single stars fall off smoothly at c=1.1.

The total number of Gaussians in the binary Gaussian mixtures is 89. For the case of the trinary population, the number of Gaussians is very large, and unnecessarily so. Therefore we do not use all possible combinations of parent stars. Rather, the three parents are given by three single stars with indices (s,s+2k,s+2l), where s,k,l are integers fulfilling that l≥ k≥ 0. This gives a total of 193 Gaussians for the trinary Gaussian mixture. Although we are thinning the number of trinary Gaussians, the total number of trinary stellar systems are unaffected; this is accounted for by a normalization factor, as we shall see in section <ref>.




 §.§ Hierarchical model


There is a total of 7 population parameters, capsuled in a vector , as is listed in table <ref>. Also visible in the table are the parameters of the Gaussian mixture (color and magnitude mean values and a 2 × 2 covariance matrix), the stellar object parameters (intrinsic color and magnitude), and the data of an observed object (observed color and magnitude and a 2 × 2 covariance error matrix). When writing the object parameters or data without an index i, we refer to the stellar parameters or data of all objects, like = {_1,_2,...,_N }. The index i always denotes a stellar object, while the index s, b, and t always denotes a Gaussian in the Gaussian mixture of the single, binary, or trinary population.



The posterior on the population parameters  and stellar parameters  is written


    Pr(,|𝐝) = 
        Pr() ∏_i Pr(𝐝_i | _i)Pr(_i | ).

The prior demands that γ, f_b, and f_t are positive and that f_b+f_t<1, and is otherwise flat in all parameters, like


    Pr() = θ(f_b) θ(f_t) θ(1-f_b-f_t) Pr(α, β_1, β_2, γ, η),

where θ is the Heaviside step function and Pr(α, β_1, β_2, γ, η) is uniform in some arbitrarily large volume in parameter space.

The posterior density of the population parameters, with stellar parameters marginalized, is


    Pr(|𝐝) = 
        Pr() ∏_i ∫d_i Pr(𝐝_i | _i)Pr(_i | ).

Because both the errors are Gaussian and the model is a sum of Gaussians, the integral over an individual object can be written as a sum over bivariate Gaussian distributions, like


    ∫d_i Pr(𝐝_i | _i)Pr(_i | ) =  
    1-f_b-f_t/N̅_single   ∑_single
    Gaussians
     sexp(-αc̅_s)𝒩(_i-_s,Σ_i+Σ_s)  
    
        + f_b/N̅_binary   ∑_binary
    Gaussians
     bexp(-αc̅_b)exp(ηΔ_c̅,b)𝒩(_i-_b,Σ_i+Σ_b)  
    
        + f_t/N̅_trinary   ∑_trinary
    Gaussians
     texp(-αc̅_t)exp(ηΔ_c̅,t)𝒩(_i-_t,Σ_i+Σ_t).

The mean color c̅_s and α sets the amplitude of the sth Gaussian of the single population, and for binaries (trinaries) there is the additional factor given by η and Δ_c̅,b (Δ_c̅,t). The quantities N̅_single, N̅_binary, N̅_trinary are normalizations equal to the integrated distribution for the population in our color-magnitude sample window. In the case of the population of single stars (and equivalently in the two other cases), we have that


    N̅_single = ∑_single
    Gaussians
     sexp(-αc̅_s)1/4[erf( 0.8-c̅_s/√(2(σ_c,s^2+σ̃_c^2)))
        +erf( c̅_s-0.5/√(2(σ_c,s^2+σ̃_c^2))) ] 
    ×[erf( 6.0-M̅_s/√(2(σ_M,s^2+σ̃_M^2)))
        +erf( M̅_s-3.0/√(2(σ_M,s^2+σ̃_M^2))) ],

which is a sum over all Gaussians in the mixture, of their respective amplitudes times the probability that a star of that Gaussian will be included in the window of parameter space. The quantities σ_c,s and σ_M,s are the width of the sth Gaussian in color and magnitude, i.e. the diagonal elements of Σ_s. The quantities σ̃_c and σ̃_M are the median errors to the color and magnitude data, accounting for the fact that observational errors can cause an object that is actually inside the allowed region of parameter space to be mistaken for being outside, and vice versa. For our sample, they have values σ̃_c=0.032 and σ̃_M=0.090. We have made some approximations here: we take the median values instead of integrating over the full distribution of errors (where this distribution could even be a function of c and M); furthermore, we do not account for the non-diagonal elements of the covariance matrices. Both of these factors are small and have minuscule bearing on our end result. The simplification of taking the median error is only applied when calculating the normalization factor N̅; the individual errors of each object is accounted for in all other factors of equation (<ref>).

The probability of a stellar object x being binary, given population parameters, is given by


    Pr( x  is binary | 𝐝_x, ) =
        f_b/N̅_binary∑_bexp(-αc̅_b)exp(ηΔ_c̅,b)𝒩(_x-_b,Σ_x+Σ_b)∫d_x Pr(𝐝_x | _x)Pr(_x | ).

The denominator is a normalization as given by equation (<ref>). The probability of an object being single or trinary is expressed in an analogous manner. The probability of being single, binary and trinary adds up to unity and can be expressed using a color code according to figure <ref>.








§ DISENTANGLING BINARY AND TRINARY SYSTEMS

Given our model, where the main sequence is narrow and Gaussian in width, it is not only possible to infer that a stellar object is a binary or a trinary system; one can also constrain the properties of the component stars. In this section we only make the binary case explicit, but analogous equations can be straightforwardly derived for the trinary case.

Given population parameters  and assuming that the stellar object with data 𝐝_x = {_x,Σ_x} is indeed a binary, the posterior on the component stars' properties, _A={c_A,M_A} and _B={c_B,M_B}, is proportional to


    Pr(_A,_B|,𝐝_x) ∝𝒩(_x-_A+B,Σ_x) Pr(_A,_B | )

where _A+B is the color and absolute magnitude of single stars A and B combined, following equation (<ref>).

The population of binaries is modeled with a Gaussian mixture, where each Gaussian has two parent single star Gaussians. However, to get the probability density of the two component stars, Pr(_A,_B | ), we have to look at individual stars. The number density of binaries decays with the color of the binary (not the colors of its component stars), according to the same exponential law as the single star population. The probability density of the two component stars becomes


    Pr(_A,_B | ) ∝
    exp(-α c_A+B)
        exp(η |c_A-c_B|)
        ∑_s𝒩(_A-_s,Σ_s)
        ∑_s'𝒩(_B-_s',Σ_s')


When inferring the properties of individual objects, we do so in the regime of Empirical Bayes, in the sense that we reuse the information on the inferred population parameters as if it was independent of the object in question[<https://en.wikipedia.org/wiki/Empirical_Bayes_method>]. This is formally incorrect but since we are working with a sample of so many objects, removing one object from the sample has a negligible effect on the population parameter posterior.






§ RESULTS





The posterior density for the 7 population parameters is visible in figure <ref>. Those were obtained by histograming samples obtained with a Monte-Carlo Markov Chain, using a Metropolis random walk of 40,000 steps, thinned by a factor 10. The step length was calibrated in a thorough burn-in phase.

The posterior density is well constrained. The value for α is largely due to selection effects, as the true abundance of stars does not decrease with color. The other parameters of the main sequence slope (β_1, β_2, and γ) and have reasonable values and are strongly constrained. The parameters of binary and trinary populations (η, f_b, and f_t) also well constrained, although these values are heavily influenced by selection effects.



In figure <ref> we see a scatter plot of a thinned set of data points, the highest posterior model without noise, and the model with noise. In the model panels, the black/red/blue correspond to single/binary/trinary systems. In the data panel, the color corresponds to the posterior probability of the object type, marginalized over the population parameters. The color scale is continuous according to figure <ref>, so some points are for example purple, signifying that binary and trinary classifications are roughly equally likely. In the inferred model, the three populations form three separate bands with very little overlap, although observational errors blurs this structure. Objects are classified as single/binary/trinary mainly as a function of brightness with respect to the main sequence, although neighboring data points can differ somewhat depending on their respective measurement errors. An object with an observed magnitude in the binary band can be consistent with the single population if the errors are large enough. Because the number density of the singles population is higher, such an object will be inferred to have a lower intrinsic absolute magnitude. Equivalently, an object in the trinary band with large errors will tend towards the binary population. In the model with noise, shown in the right panel, the points are given an additional scatter corresponding to the noise of the data. In order to account for a magnitude dependence to the noise, we have binned the sample region into six bins in absolute magnitude, such that the noise in the right panel is drawn at random from the corresponding bin in the data panel.



Figure <ref> shows a summary statistic of the inferred object types. A majority of objects are very likely single stars, while strongly inferred binaries are less numerous and trinaries even more so. While very strongly inferred singles and binaries should indeed be interpreted as such, it is possible that the objects strongly inferred to be trinaries are in fact higher multiples. Including quaternary systems in this analysis would probably make an even better fit for the objects furthest from the main sequence.





In figures <ref> and <ref>, we see the posterior on the color and magnitude of the component stars of a multiple stellar system, where the posterior is marginalized over the population parameters. These Monte-Carlo Markov Chains were found using the Emcee implementation <cit.> of the affine-invariant ensemble sampling method <cit.>. Both figures are for the same object, but in figure <ref> it is assumed that this object is a binary system, while in figure <ref> it is assumed that it is trinary. This object[Gaia DR1 source ID: 6220743838828320384 – 2MASS identification number: 14230246-3029168] is chosen as it has an equal probability between these two possibilities: it has a 50.0 % probability of being binary and a 50.0 % probability of being trinary. Its observed color and magnitude are ĉ=0.512 and M̂=3.54. In both these figures, the posteriors are constrained by the criterion that colors are in ascending order, such that c_A < c_B (< c_C). Clearly, the color and magnitude of the brighter component star is much better constrained than those of the dimmer component star(s).



In figure <ref>, we turn to the objects of our catalog  that are inferred to be binaries with a probability greater 0.5, and show the maximum posterior values of the component star colors for objects. It is clear that objects for which the difference in color between the component stars is large, seen towards the top of the figure, are also more consistent with being single stars, as seen by the dark red scatter points. Conversely, objects where the difference in color between the component stars is small, towards the bottom of the panel, are also more consistent with being trinaries, as seen by the purple scatter points.




§ DISCUSSION AND CONCLUSION


We developed a flexible data-driven model of the color-magnitude diagram simultaneously fitting for populations of single stars as well as unresolved binaries and trinaries for the first time. We applied this model to photometric and astrometric data from the Gaia TGAS cross-matched with the 2MASS survey, and found that this model can explain the scatter and distribution of objects in the data. The same method can be applied to any data set where the distances to the stellar objects are known to adequate accuracy, such as an open cluster or a nearby local group galaxy.

The model itself relies on describing the population of single stars with Gaussian distributions positioned on a line, with amplitude exponentially decreasing as a function of color. 
With this model at hand, the populations of binaries and trinaries can be predicted, and we showed that despite the non-linear transformation they are also well approximated by Gaussians. Therefore, the full model is parametrized by the slope and intercept of the line supporting the Gaussians describing single stars, their covariances, the decay rate of their amplitudes as a function of color, as well as the fractional probabilities of objects being binaries or trinaries. This approach provides an accurate description of the data and is suitable for efficient inference of its parameters.

With the parameters of this model constrained from the TGAS-2MASS data, we derived probabilities for individual objects to be single, binary, or triple systems, as well as making predictions for their individual colors and magnitudes (typically strongest for the system's brightest star).
Those results are robust, with all sources of uncertainties fully included, or equivalently, all parameters marginalized over.
They minimally rely on assumptions about the underlying physics, not depending on rigid stellar models but rather on a flexible data-driven model providing a good description of the data. For simplicity, we did not include higher multiples in our analysis, as a population of singles, binaries, and trinaries describes the overall population very well. However, the systems strongly inferred to be trinary could very well be quaternary system instead. In the inferred model, visible as the center panel of figure <ref>, the three populations form separate bands. With more precise parallax measurements such bands would be clearly discernible, making the inference on the model more powerful and robust.

One of the main limitations of this analysis is the absence of a model for the selection function (i.e., the detection and noise properties of the objects in our sample). We know for a fact that selection effects are very severe, coming from the TGAS catalogue and the complex scanning law of the Gaia survey, and also the TGAS/2MASS cross correlation. Selection effects are dependent on position on the sky, apparent magnitude, and proper motion, among other effects, and are difficult to account for. Furthermore, the chance of confusing two stars as one object is a function of distance (a binary systems might be resolved if nearby but not if far away). For those reasons, we cannot fully interpret f_b and f_t as actual fractions of binary and trinary systems in the solar neighborhood. Rather, they are the fraction of unresolved binary and trinary systems in our sample, where selection effects are included.

Future Gaia Data Releases will allow us to improve this model and resolve some of those limitations. 
The selection function will be simpler, and Gaia colors will be available, therefore reducing the need for external data like 2MASS.
We also intend to extend this work in the following ways: include higher multiples, model other regions of color-magnitude space, use a less restrictive cut in parallax, develop more sophisticated models of the data involving Gaussian Mixtures and stellar models. A major improvement would be to account for how the resolving power decreases with distance, which could give a more robust result and constrain the distribution of orbital separation in multiple stellar systems.



We would like to thank our referee, Maxwell Moe, for a constructive peer review that was of great benefit to this article. This work was supported by collaborative visits funded by the Cosmology and Astroparticle Student and Postdoc Exchange Network (CASPEN). BL was supported by NASA through the Einstein Postdoctoral Fellowship (award number PF6-170154). DWH was partially supported by the NSF (AST-1517237) and the Moore-Sloan Data Science Environment at NYU.





This work has made use of data from the European Space Agency (ESA) mission Gaia[<http://www.cosmos.esa.int/gaia>], processed by the Gaia Data Processing and Analysis Consortium[<http://www.cosmos.esa.int/web/gaia/dpac/consortium>] (DPAC). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.
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 research utilized the following open-source Python packages: Matplotlib <cit.>, Emcee <cit.>. This work additionally used the Gaia science archive (https://gea.esac.esa.int/archive/) and the NASA Astrophysics Data System.






















