---
author:
- 'T. Kaczmarek'
- 'C. Olczak'
- 'S. Pfalzner'
bibliography:
- 'reference.bib'
title: Evolution of the binary population in young dense star clusters
---

[ Field stars are not always single stars, but can often be found in bound double systems. Since binary frequencies in the birth places of stars, young embedded clusters, are sometimes even higher than on average the question arises of how binary stars form in young dense star clusters and how their properties evolve to those observed in the field population. ]{} [We assess, the influence of stellar dynamical interactions on the primordial binary population in young dense cluster environments.]{} [ We perform numerical N-body simulations of the Orion Nebula Cluster like star cluster models including primordial binary populations using the simulation code .]{} [ We find two remarkable results that have yet not been reported: The first is that the evolution of the binary frequency in young dense star clusters is independent predictably of its initial value. The time evolution of the normalized number of binary systems has a fundamental shape. The second main result is that the mass of the primary star is of vital importance to the evolution of the binary. The more massive a primary star, the lower the probability that the binary is destroyed by gravitational interactions. This results in a higher binary frequency for stars more massive than $2\Msun$ compared to the binary frequency of lower mass stars. The observed increase in the binary frequency with primary mass is therefore most likely not due to differences in the formation process but can be entirely explained as a dynamical effect. ]{} [ Our results allow us to draw conclusions about the past and the future number of binary systems in young dense star clusters and demonstrate that the present field stellar population has been influenced significantly by its natal environments.]{}

Introduction
============

Today the formation process of isolated single stars seems to be understood in principle. One remaining uncertainty, at least to some extent, concerns the formation of binary stars. There are several ways to form binary stars - the fragmentation of a collapsing cloud or core or the dynamical capture of a companion. However, determining the dominating process of binary formation requires knowledge of the parameters of the primordial binary population, i.e. the population as it appears just after star formation has finished.

Several surveys have been performed to determine the present-day parameters of binary populations in different kinds of environments. Observations of the field binary population have shown that the binary frequency is quite high, $60\%$ for G-type dwarf stars and $35 - 42\%$ for M-type dwarfs [@1992ApJ...396..178F; @1997AJ....113.2246R]. Both and [@1992ApJ...396..178F] determined the mass ratio and separation/period distributions of their samples. They found that the mass ratio distribution for G-type dwarfs increases as the mass ratios decreases while the mass ratio distribution for M-dwarfs is flat. The separation/period distributions for G-type and M-type dwarfs can be fitted by a Gaussian distribution in logarithmic space where both peak at nearly the same separation ($\sim 30 \AU$ for G-type and $9-30\AU$ for M-type dwarfs) with the FWHM of both distributions spanning over three orders of magnitude .

Another feature of binary populations is their dependence on the primary mass: observations indicate that the probability of finding a massive star in a binary system is higher than for low-mass stars. For OB stars, a binary frequency of $\sim 70\%$ is found [@2002AJ....124..404P], whereas for G-dwarfs a binary frequency of $\sim 60\%$ and for M-dwarfs a value as low as $\sim 31\%$ [@2003IAUS..211..311M] has been observed.

An understanding of the presently observed properties of evolved binary populations requires knowledge of the primordial binary population. To achieve this one has to look at sites where star formation has just happened, such as young star clusters. However, observations of binary populations in star clusters do not provide such a clear picture as the observations of the field binary population. While in dense star clusters such as the Orion Nebula Cluster (ONC) or NGC 2024 the observed binary frequency is similar to the field binary frequency , sparse clusters such as the Taurus-Auriga association have much higher binary frequencies . This leads to the question of what defines the binary frequency in the different environments. Possible explanations are differences in the temperatures of the parent clouds or the dynamical evolution within the cluster [e.g. @1999NewA....4..495K]. The latter possibility will be investigated in this study.

According to current knowledge, stars do not form in isolation, but in a clustered environment consisting of several tens to thousands of stars. Most of them possibly formed in clusters such as the ONC or even denser environments . The stellar density in ONC-like clusters is sufficiently high to potentially affect the evolution of the binary population via gravitational interactions.

The gravitational interactions that alter binaries in a star cluster are mostly three- and four-body interactions between binary and single stars. @1975MNRAS.173..729H [@1975IAUS...69...73H] and @1975AJ.....80..809H showed that three-body interactions are very effective in altering the parameters of binary systems. Both authors found that loosely bound binaries (also called soft) tend to become softer (meaning that their internal energy increases resulting in the widening of the binaries) owing to three-body interactions, while firmly bound (hard) binaries tend to become harder (Heggie-Hills law).

Four-body interactions (binary-binary encounters) become important in clusters that consist of a high number of binary stars [see for example Fig. 1 in @2010MNRAS.tmp.1645L]. In addition to three-body interactions, four-body interactions can form two “new” binary systems by exchanging the binary components, can destroy one or even both binary systems, or can lead to the formation of bound hierarchies of three or more stars. [@2003gnbs.book.....A; @1999JCoAM.109..407S], the simulation code used for this study, uses several techniques such as the regularization of three-body interactions [@1974CeMec..10..185A] or the chain method [@1993CeMDA..57..439M] to compute the evolution of these interactions correctly. However, in this study we only diagnose the evolution of the binary population.

Numerical investigations of the dynamical evolution of binaries in star clusters started in the 70’s with the development of the first N-body codes . The early works concentrated on the formation of binaries in star clusters with initially pure single stellar systems. Only in recent years simulations have been performed to investigate the evolution of a primordial binary population. For example, @1999NewA....4..495K simulated the evolution of ONC-like star cluster models where all stars initially were part of a binary system, i.e. a primordial binary frequency of $100\%$. They showed that the cluster evolution destroys wide binaries quite efficiently leading to a reduction in the initial binary fraction as the cluster evolves.

Here we also want to study the evolution of the primordial binary population of a young dense star cluster. In contrast to the previous studies, we investigate the effect of the initial binary frequency ($0-100\%$) and a log-uniform semi-major axis distribution similar to that observed in very young star forming regions [e.g. @2007AJ....134.2272R; @2008AJ....135.2526C]. Because observations indicate that the binary frequency increases with stellar mass, we also focus on the dependence of the binary evolution on stellar mass.

This paper is structured in the following way. In Sec. \[sec:definitions\] we define parameters that describe a binary population. Sec. \[sec:method\] covers the setup of the simulations, including the generation of the cluster models and initialization of the primordial binary population. Finally in Sec. \[sec:results\] we present the results obtained from the simulations.

Definitions {#sec:definitions}
===========

The most simple measure of the binary population of a star cluster is the number of binaries () within the cluster. Here, a pair of stars is called a binary system if the system is bound and the two stars are at the same time nearest neighbors. In this definition, a star that is part of a multiple system can be part of several “binary” systems simultaneously since two stars do not have to be mutually nearest neighbors. Hence triple or other higher order multiples are not handled separately by our diagnostic but are counted as multiple binary systems. In this study, we mostly use the normalized number of binary systems  defined as $$\label{eq:nbin_norm}
  \Nbinnorm (t) = \frac{\Nbin (t)}{\Nbin (0)},$$ where $\Nbin (0)$ denotes the initial and $\Nbin (t)$ the number of binary systems at time $t$ within the cluster. The usual measure of binary populations in star clusters is the binary frequency   $$\label{eq:binary_frequency}
  \binaryfrequency (t) = \frac{\Nbin (t)}{\Nsing (t)+\Nbin (t)},$$ where $\Nsing (t)$ denotes the actual number of observed single stars and $\Nbin (t)$ the actual number of binary systems. One major result of our work is that the normalized number of binary systems, Eq. \[eq:nbin\_norm\], is the appropriate measure of a binary population.

Method and cluster setup {#sec:method}
========================

Here the evolution of binary populations in young dense star clusters is studied using the parallized version  of the high-accuracy N-body code  [@2003gnbs.book.....A; @1999JCoAM.109..407S]. The study presented can be understood as a follow up of . The main difference to is that here the clusters contain a primordial binary population, with varying initial binary frequencies. In both studies, the ONC has been chosen as a model cluster because it is one of the most well-studied star forming regions for which many parameters are well constrained [see for example @1994ApJ...421..517P; @1997AJ....113.1733H; @1998ApJ...492..540H], in particular several surveys of its binary population exist. This allows a detailed comparison of the simulated binary populations with the observed ones.

Our simulations started with 4000 stellar systems, where stellar systems correspond to either single stars or binary systems. This means that for a simulation with a primordial binary frequency of 100% the dynamical evolution of 8000 stars is calculated. The masses $M$ of the 4000 stellar systems were sampled from the initial mass function given by @2002Sci...295...82K $$\label{eq:kroupa-imf}
  \xi(M) =
  \begin{cases}
    \; M^{-1.3} \quad & \textrm{if} \quad 0.08 \leq M \left[ \Msun \right]< 0.5, \\
    \; M^{-2.3} \quad & \textrm{if} \quad 0.5 \leq M  \left[ \Msun \right]< 1.0, \\
    \; M^{-2.3} \quad & \textrm{if} \quad 1.0 \leq M  \left[ \Msun \right]< \infty.
  \end{cases}$$ The minimum and maximum masses of the stellar systems have been chosen as $M_{\textrm{min}}=0.08 \Msun$, corresponding to the hydrogen-burning limit, and $M_{\textrm{max}}=50\Msun$ the mass of the most massive stellar system in the ONC - .

Observations of the stellar content in the ONC have shown that the stellar radial density profile of the cluster can be approximated by an isothermal sphere (i.e. $\rho_{\mathrm{present}} \propto 1/r^{2}$, see e.g. [@1988AJ.....95.1755J; @1997AJ....113.1733H]) with the exception of the cluster core ($R_{\mathrm{core}}\approx 0.2\pc$), which has a much flatter density profile ($\rho_{\mathrm{present,core}} \propto r^{-0.5}$, @2005MNRAS.358..742S). To account for this, we adopt an initial density profile of the form $$\label{eq:initial_radial_density_profile}
  \rho_{0} (r) = 
  \begin{cases}
    \rho_{0}\; (r/R_{\mathrm{core}} )^{-2.3}&,\quad r\in (0,R_{\mathrm{core}}] \\
    \rho_{0}\; (r/R_{\mathrm{core}} )^{-2.0}&,\quad r\in (R_{\mathrm{core}}, R]\\
    \qquad 0&, \quad r\in (R,\infty]
  \end{cases}$$ where $\rho_{0} = 3.1\times 10^{3}\pcdens$, $R_{\mathrm{core}} = 0.2\pc$, and $R=2.9\pc$. demonstrated that this initial radial density profile evolves towards the currently observed one after about $1\Myr$ of dynamical evolution, the estimated age of the ONC [@1997AJ....113.1733H].

Observations of the ONC indicate that the cluster are mass-segregated [e.g. @1997AJ....113.1733H]. Numerical simulations by @1998MNRAS.295..691B comparable to the cluster models studied here demonstrated that this mass segregation cannot be due to dynamical evolution of the cluster and therefore must be primordial. Following their suggestion, we initially place the most massive stellar system at the cluster centre and the three next most massive stars at random positions within a sphere of radius $R=0.6\Rhm$ centred on the cluster centre, where  is the half-mass radius of the cluster.

The velocities of the stellar systems are sampled from a Maxwellian velocity distribution in virial equilibrium. The resulting velocity dispersion of our model is lower than the observational value [$\sigma = 4.3\pm
0.5\kms$: @1988AJ.....95.1755J]. However, because of significant mass of the molecular cloud associated with the ONC that was neglected in our numerical model, the observed velocity dispersion is consistent without any assumption about a stellar system in dynamical equilibrium. In future investigations of the binary populations, we plan to include the effect of the gas associated with the cluster.

For the primordial binary population, we chose the model described in . They compared observations of the Scorpius OB II association with simulated observations of their binary population model to eliminate observational selection effects and biases that affect real observations of star clusters. Their analysis shows that the binary frequency in Scorpius OB II is above $70\%$ with a confidence of $3\sigma$. The closest agreement with observations is obtained for models with a binary frequency of $100\%$. The model of was chosen here since it is a well-studied binary population, close to its initial conditions. Scorpius OB II is a relatively young (about $5-20\Myr$) and probably too sparse region to have undergone significant dynamical evolution.

In contrast to the field population that exhibits a log-normal period distribution of G and M-dwarf binaries , the younger population in Sco OB II has a flat distribution of the semi-major axis in logarithmic space $$\label{eq:kouwenhoven07-oepiks-law}
  f_a(a) \propto a ^{-1} \qquad \Leftrightarrow \qquad f_{\log_{10} a} (\log_{10} a) = \mathrm{const.},$$ where $a$ is the semi-major axis of the binary system.

A flat distribution of the semi-major axes in logarithmic space corresponds to a similar distribution of periods $P$ with the relation $$\begin{aligned}
  \label{eq:kouwenhoven07-period-separation-distribution}
  \overline{ \log_{10} a } & =  \frac{2}{3}\overline{ \log_{10} P } - \frac{1}{3} \log\left( \frac{4\pi^2}{G\;M_T}\right), \\
  \sigma_{\log_{10} a} & = \frac{2}{3}\sigma_{\log_{10} P},\end{aligned}$$ where $G$ is the gravitational constant and $M_T$ the total mass of a binary system.

A flat binary period distribution has not only been observed in Sco OB II but also in other regions [e.g. @2008AJ....135.2526C; @2007AJ....134.2272R] and has been used in several theoretical studies involving the evolution of primordial binary frequencies [@1992MNRAS.257..513H; @2005MNRAS.358..572I; @2007ApJ...665..707H].

In addition to a log-uniform shape of the initial period distribution, other distributions have been used in numerical studies of the dynamical evolution of binary populations. For example, @1995MNRAS.277.1491K demonstrated that a dynamically evolved period distribution that is initially rising for short periods and nearly constant for long periods matches the observed period distribution of the field G-type dwarfs . In contrast, this cannot be achieved with an initially log-uniform period distribution because of the much larger number of dynamically stable short-period binaries. However, physical processes that were not taken into account by @1995MNRAS.277.1491K could significantly lower the number of short period binaries. For example “the orbital decay of embedded binary stars” [@2010MNRAS.402.1758S] is capable of reshaping a log-uniform binary period distribution within less than 1  towards a log-normal distribution (Hövel et al., in preparation). Consequently, as the effect of different initial period distributions on the evolution of the binary population in star clusters is basically unknown we will explicitly focus on this topic in a future investigation. In this context the present study can be understood as preparatory work.

Figure \[fig:initial\_semi-major\_axis\_distribution\] shows the initial semi-major axis distribution of our simulations. Although the semi-major axis distributions have been set up log-uniformly in the range $4.6\times 10^{-2} \AU\; - \; 4.6\times 10^4 \AU$, the distribution is only flat up to about $2\times 10^3\AU$ and declines afterwards to zero. The apparent lack of wide binary systems in the $2\times 10^3 - 4.6\times 10^4\AU$ range is due to the generating algorithm. In the code, binary systems are formed in isolation and then placed in the cluster environment. The wider a binary, the higher the probability that (at least) one of its components has a closer neighbour among the other cluster members. According to our definitions (see Sec. \[sec:definitions\]), this chance neighbour is now the potential companion and most probably not bound.

![Initial semi-major axis distribution for all simulated clusters. Crosses, stars, open and filled squares denote simulations with initially 30%, 50%, 75%, and 100% binaries, respectively.[]{data-label="fig:initial_semi-major_axis_distribution"}](number_of_bin_sys_vs_semi__all_binaries__0_0_Myr__light){width="\linewidth"}

also found the mass ratio distribution to have the form $$\label{eq:kouwenhoven07-mass-ratio-distribution}
  f_q(q)\propto q^{\gamma_q}, \qquad \gamma_q \sim -0.4,$$ where $q = \Msec / \Mprim$ is the mass-ratio of a binary with primary mass $\Mprim$ and secondary mass $\Msec$. Adopting this functional form in our models, the binary system masses were sampled from the @2002Sci...295...82K IMF and then divided using Eq. \[eq:kouwenhoven07-mass-ratio-distribution\]. This approach, however, does not reproduce the observed stellar IMF as can be seen in Fig. \[fig:comparison\_imf\_kroupa\] (thin lines), which displays the stellar IMFs produced in our simulations (filled squares, representing simulation with initially $100\%$ binaries). For comparison, the standard IMF was included by setting up simulations using random pairing of the binary components instead of using Eq. \[eq:kouwenhoven07-mass-ratio-distribution\] to pair the binary components. Random pairing has been shown to reproduce the stellar IMF and is here therefore used to represent the standard IMF. Using Eq. \[eq:kouwenhoven07-mass-ratio-distribution\] to pair the binary components generates too few low-mass stars in the range $0.08 \Msun - 0.16\Msun$ and too many intermediate mass stars in the mass range $0.15 - 0.5\Msun$, where the deviations become larger the higher the initial binary frequency is. The number of stars generated with masses exceeding $0.5\Msun$ reproduces the @2002Sci...295...82K IMF quite well, although their numbers are slightly under-predicted.

The thick lines in Fig. \[fig:comparison\_imf\_kroupa\] show the initial binary-system mass ($m_{\mathrm{Hy's} = \Mprim + \Msec}$) functions for the simulations used in this study and additionally the initial binary-system mass function if the binary components are drawn randomly from the stellar IMF. The binary system mass function in both cases is very similar. We do not expect the differences in the single-star IMFs to affect the evolution of the binary populations because it depends on the system masses as we show in Sec. \[sec:results\].

![Comparison of the stellar initial mass function (thin lines) and the initial binary-system mass function (thick lines). The open circles denote simulations that used random pairing to pair the binary components while the solid squares denotes simulations that used Eq. \[eq:kouwenhoven07-mass-ratio-distribution\].[]{data-label="fig:comparison_imf_kroupa"}](multiplot_relative_stellar_IMF_number_bin_sys_vs_mass){width="\linewidth"}

The eccentricity of the binary systems in Sco OBII was found by to be thermally distributed, in agreement with the eccentricities of long period binaries in the field. Following @2001ApJ...555..945K who illustrated that the dynamical evolution is not capable of changing the eccentricity distribution significantly, a thermal period distribution was applied to all of our simulations.

Observations show that the binary frequency varies among low-mass stars in the field ($35-42\%$) [@1992ApJ...396..178F; @1997AJ....113.2246R], Sun-like stars in the solar neighbourhood ($\sim 60\%$) , and the high binary frequency in low-mass star-forming regions such as Taurus ($\sim 100\%$) . We account for the large spread among the different environments by investigating cluster models initially with $30\%$, $50\%$, $75\%$, and $100\%$ binary frequency.

To improve the statistical significance of our results, for each initial binary frequency 100 simulations were performed. Only those that fulfilled the following conditions were considered in the subsequent diagnostic process:

-   after 1  of cluster evolution, the most massive stellar system in the simulations has to be within $1\pc$ of the cluster centre.

-   the cluster density profile after 1  of cluster evolution has to agree with the cluster density profiles observed by [@1997AJ....113.1733H] and @2002Msngr.109...28M.

Results {#sec:results}
=======

Observations of binary populations in different environments usually use the binary frequency as the main observable of a binary population since it can easily be compared between different environments and different observing campaigns. Theoretical studies of the evolution of binary populations are unaffected by from this shortcoming but usually also use this simple observable for reasons of comparison with observational data [see for example @1995MNRAS.277.1491K; @2001MNRAS.321..699K; @2003MNRAS.346..343K; @2005MNRAS.358..572I; @2007ApJ...665..707H; @2008MNRAS.388..307S]. If we also do this here, we obtain the result described by Fig. \[fig:evolution\_normalized\_bf\_and\_number\_binary\_systems\]a. It shows the evolution of the normalized binary frequency during the first three  of cluster evolution in the simulations. The normalization is performed by dividing the given binary frequencies by the initial binary frequency and allows direct comparison of the evolution of the binary frequencies for simulations with different initial binary frequencies.

As expected from previous studies [@1995MNRAS.277.1491K; @2003MNRAS.346..343K; @2009MNRAS.397.1577P], dynamical evolution reduces the binary frequency with time in all simulations. This can be explained by the destruction of wide binary systems by means of three body encounters [see @1995MNRAS.277.1491K; @1995MNRAS.277.1507K for more details]. Apart from this general trend, Fig. \[fig:evolution\_normalized\_bf\_and\_number\_binary\_systems\]a shows that the evolution of the binary population depends on the initial binary frequency in the sense that the higher the initial binary frequency, the higher the fraction of binaries that are destroyed during cluster evolution. This finding could lead to the conclusion that the evolution of binary populations is stronger in clusters with a higher initial binary frequency. An analogous dependence was found by @2008MNRAS.388..307S, who investigated the evolution of the binary population in globular star clusters using an analytical approach.

![Evolution of the normalized binary frequency, normalized number of binary systems and normalized number of single stars with time. Crosses, stars, open, and filled squares denote simulations with initially 30%, 50%, 75%, and 100% binaries, respectively.[]{data-label="fig:evolution_normalized_bf_and_number_binary_systems"}](multiplot_normalized_binary_frequency__number_of_binary_systems__number_of_single_stars_vs_time__all_binaries){width="\linewidth"}

In short, the evolution of the binary population in star clusters seems to depend on the initial binary frequency. However, because the calculation of the *binary* frequency depends on the number of *single* stars (see Eq. \[eq:binary\_frequency\]) it is probably not an unbiased measure of the evolution of the binary population. We thus introduce the normalized number of binary systems  (see Eq. \[eq:nbin\_norm\]) as the proper quantity and show its evolution in Fig. \[fig:evolution\_normalized\_bf\_and\_number\_binary\_systems\]b.

As expected from Fig. \[fig:evolution\_normalized\_bf\_and\_number\_binary\_systems\]a, the normalized number of binaries also decreases with time. However, in contrast to the normalized binary frequency the evolution of the normalized number of binary systems with time does *not* depend on the initial binary frequency. The difference among the data of only $1\%$ is of the order of the statistical error in the simulation. How can this difference in the evolution of the binary frequency and the binary population be explained? It is the evolution of the single star population in the star clusters that significantly affects the evolution of the binary frequency. Its evolution is show in Fig. \[fig:evolution\_normalized\_bf\_and\_number\_binary\_systems\]c.

In contrast to the binary frequency and the number of binaries, the number of single stars in the star cluster obviously increases with time because the dissolution of a binary system creates two new single stars. The evolution of the normalized number of single stars, however, depends even more strongly on the initial binary frequency than the evolution of the binary frequency.

We show by means of simple arithmetic why the evolution of the single star population depends on the initial binary frequency while the evolution of the binary population does not. First, we assume that the evolution of the number of binary systems can be described by a function $\alpha (t,\binaryfrequency,\dots)$ that may depend on the time of evolution, initial binary frequency, and other parameters, so that $\Nbin (t) = \alpha (t,\binaryfrequency,\dots) \Nbin(0)$, where $\Nbin (t)$ is the number of binary systems at time $t$. Dividing this by the initial number of binary systems we obtain the evolution of the normalized number of binary systems $$\label{eq:evolution_normalized_number_of_binary_systems}
\Nbinnorm (t) = \frac{\alpha (t,\binaryfrequency,\dots) \Nbin(0)}{\Nbin(0)} = \alpha (t,\binaryfrequency,\dots).$$ If we also assume, based on the findings of Fig. \[fig:evolution\_normalized\_bf\_and\_number\_binary\_systems\]b, that the evolution of the normalized number of binary systems does not depend on the initial binary frequency, the function $\alpha$ can be expressed as $\alpha (t,\binaryfrequency,\dots ) \approx \alpha (t)$. Because the dissolution of a binary system leads to the generation of two new single stars, the evolution of the number of single stars can be written as $$\label{eq:evolution_number_of_single_stars}
\Nsing (t) = \Nsing (0) +  2\left[ 1-\alpha (t,\binaryfrequency,\dots) \right] \Nbin (0),$$ where the first summand represents the initial number of single stars and the second summand gives the number of newly formed single stars. Applying the relation $\Nbin (0) = \frac{\binaryfrequency (0)}{1-\binaryfrequency (0)} \Nsing(0)$ and dividing by the initial number of single stars $\Nsing (0)$, we find that the evolution of the normalized number of binary systems is given by $$\label{eq:evolution_normalized_number_of_single_stars}
\Nsingnorm (t) = 1 + 2( 1 - \alpha (t) ) \frac{\binaryfrequency (0)}{1 - \binaryfrequency (0)}.$$ Hence, in contrast to the normalized number of *binary* systems (Eq. \[eq:evolution\_normalized\_number\_of\_binary\_systems\]), which depends only on $\alpha$, the normalized number of *single* stars depends explicitly on the initial binary frequency in the cluster. Consequently, using Eqs. \[eq:evolution\_normalized\_number\_of\_binary\_systems\] and \[eq:evolution\_normalized\_number\_of\_single\_stars\], we derive an analogous dependence for the normalized binary frequency $$\begin{aligned}
\label{eq:evolution_normalized_number_of_binary_frequency}
^{\binaryfrequency (t)} / _{\binaryfrequency (0)} &= \frac{\Nbin (t)}{\Nsing (t) + \Nbin (t)} \frac{1}{\binaryfrequency (0)} \notag\\
                    &= \frac{\alpha (t)}{\frac{1-\binaryfrequency (0)}{\binaryfrequency (0)} - \alpha (t) + 2} \frac{1}{\binaryfrequency (0)}\end{aligned}$$ In summary, our simple model (Eqs. \[eq:evolution\_normalized\_number\_of\_binary\_systems\]-\[eq:evolution\_normalized\_number\_of\_binary\_frequency\]) describes the evolution presented in Fig. \[fig:evolution\_normalized\_bf\_and\_number\_binary\_systems\] very well: it is the evolution of the single star population that is very sensitive to the initial binary frequency, while the binary population evolves independently of the initial binary frequency. The normalized number of binary systems (Eq. \[eq:nbin\_norm\]) introduced here is therefore the most simple measure of the evolution of a binary population *independent* of the initial binary frequency.

As a consequence, we can use Eq. \[eq:evolution\_normalized\_number\_of\_binary\_frequency\] to calculate the initial binary frequency in young ONC-like star clusters from only a few observables. Solving Eq. \[eq:evolution\_normalized\_number\_of\_binary\_frequency\] for $\binaryfrequency (0)$ reveals that $$\label{eq:initial_binary_frequency}
  \binaryfrequency (0) = \frac{\binaryfrequency (t) - \alpha(t) + 1}{\alpha (t)}.$$ This finding has a fundamental consequence: if we are able to determine the explicit form of the function $\alpha(t,
\binaryfrequency, \dots)$, we are able to predict the past and future number of binary systems from the present number of binaries in any cluster environment. If we again assume that $\alpha (t,\binaryfrequency, \dots) = \alpha (t)$, we find that $\alpha_{\mathrm{fit}}(t) = at^{b} + ct + 1$ with $a=-0.35 \pm 0.04, b=0.89 \pm 0.01 $, and $c =0.27 \pm 0.04 $, provides an excellent fit to the data in Fig. \[fig:evolution\_normalized\_bf\_and\_number\_binary\_systems\]b. However, we emphasize that our assumptions neglects the potential effects of other parameters such as the stellar density or distribution of binding energy, which will be investigated in a forthcoming paper.

Nonetheless, for illustrative purposes we assume that the ONC has cluster age of $1\Myr$ @1997AJ....113.1733H and obtain $\alpha (1\Myr) = 0.92$ using the given fit. Assuming a present-day binary frequency of $60\%$, motivated by observations that find a consistency between the binary populations of field G dwarfs and the ONC , we derive an initial binary frequency of 74%. This is clearly below the binary frequency observed in young sparse star-forming regions as Taurus-Auriga and suggests that the initial binary frequency depends on the star-forming environment as stated by .

We note that this result is based on our specific cluster model, which does not include all the physical processes present in real star clusters. Thus, further investigations may alter the determined initial binary frequency to some extent. For example, including the gas the young clusters are embedded in would probably require a higher initial density since star clusters expand significantly after gas expulsion [e.g. @2007MNRAS.380.1589B]. Higher densities should increase the fraction of dynamically destroyed binary systems during cluster evolution, resulting in larger values of $\alpha$. This means that the ONC could have been born with a higher initial binary frequency than estimated from our current simulations [as demonstrated by @2001MNRAS.321..699K].

We next use our homogeneous data set to investigate why the binary population evolves independently of the initial binary frequency. Figure \[fig:fraction\_number\_of\_bins\_vs\_semi\] shows the fraction of binary systems per semi-major axis bin that is destroyed after 1 of cluster evolution. It is evident that only binaries with semi-major axis wider than about $\approx
30\AU$ are destroyed during cluster evolution. This dependence on the semi-major axis was previously reported in other studies of the evolution of binary populations [e.g. @1995MNRAS.277.1507K] and can be explained by the outcome of three body encounters [for details see @1975MNRAS.173..729H]: to dissolve a binary system in a three body interaction, the kinetic energy of the intruding star must exceed the binding energy of the binary system, or, in other words, the velocity of the intruder must be larger than the orbital velocity of the binary system. With this, we estimate a lower limit to the semi-major axis of destroyed binaries in the cluster. Following @2003MNRAS.346..343K, we estimate the “thermal period” $\log_{10}
P_{\mathrm{th}} \approx 4.75$ at which destruction still is possible from $$\label{eq:thermal_period}
\log_{10}P_{\mathrm{th}} = 6.986 + \log_{10} m_{\mathrm{sys}} - 3\log_{10} \sigma_{3\mathrm{D}},$$ using the least massive equal-mass system $m_{\mathrm{sys}} = 0.16\Msun$ and assuming a typical relative velocity equal to the cluster velocity dispersion in our simulations $\sigma_{3D} \approx 3\kms$. Our result, $\log_{10} P_{\mathrm{th}} \approx
4.75\days$, corresponds to a minimum semi-major axis of about $53\AU$.

![Fraction of binary systems per semi-major axis bin that is destroyed after 1 of cluster evolution. Crosses, stars, open, and filled squares denote simulations with initially 30%, 50%, 75%, and 100% binaries, respectively.[]{data-label="fig:fraction_number_of_bins_vs_semi"}](individual_relative_deviations_of_number_of_bin_sys_vs_semi__all_binaries__1_0_Myr__light){width="\linewidth"}

This estimate is in good agreement with the results of our simulations in Fig. \[fig:fraction\_number\_of\_bins\_vs\_semi\]. However, in our simulations, there is no sharp transition between stable and destroyed systems because the thermal period $P_{th}$ depends strongly on the system mass (Eq. \[eq:thermal\_period\]). Since the periods of the binaries in the simulations have been sampled independently of the binary system masses the fraction of binaries destroyed by the dynamical evolution of the cluster decreases with binary mass. This is demonstrated in Fig. \[fig:fraction\_number\_bin\_sys\_vs\_sys\_mass\]. The fraction of destroyed binaries clearly decreases with binary mass.

The solid line in Fig. \[fig:fraction\_number\_bin\_sys\_vs\_sys\_mass\] represents a fit of the function $f( m_{\mathrm{sys}}) = a
\log_{10} m_{\mathrm{sys}} + b$ to the data of the simulations with initially $75\%$ binaries in the range from $m_{\mathrm{sys}}
= 0.16\Msun$ up to $10\Msun$. For higher masses, the results suffer from low number statistics and have therefore been excluded from the fit. We note that $f(m_{\mathrm{sys}})$ also fits the data of the other simulations demonstrating that the evolution of the binary population is independent of the initial binary frequency.

![Fraction of binaries per binary system mass bin that are destroyed after 1 of cluster evolution. The solid line is the fit of the function $f( m_{\mathrm{sys}}) = a \log_{10} m_{\mathrm{sys}} + b$ to the data with initially $75\%$ binaries. Crosses, stars, open, and filled squares denote simulations with initially 30%, 50%, 75%, and 100% binaries, respectively.[]{data-label="fig:fraction_number_bin_sys_vs_sys_mass"}](individual_relative_deviations_of_number_of_bin_sys_vs_sys_mass__all_binaries__1_0_Myr__light){width="\linewidth"}

The strong mass dependence of the binary destruction is further illustrated in Fig. \[fig:normalized\_bin\_freq\_vs\_time\_\_masses\]a. It shows the evolution of the normalized number of binaries with primary masses either exceeding $2\Msun$ (thin lines) and lower than $2\Msun$ (thick lines) in clusters with various initial binary frequencies. Clearly, a larger fraction of low-mass binaries is destroyed during the dynamical evolution of the cluster.

The boundary of $2\Msun$ between high and low-mass primaries was chosen in accordance with the study of , who performed a high spatial resolution infrared study of the binary population in the ONC. They found the frequency of binaries with high-mass primaries to be significantly higher. Figure \[fig:normalized\_bin\_freq\_vs\_time\_\_masses\]a implies that this finding can be explained by the dynamical evolution of the cluster and does not depend on the initial binary frequency.

However, a consistent comparison of our simulation and the observations by must focus on the frequency and not the number of binaries. Figure \[fig:normalized\_bin\_freq\_vs\_time\_\_masses\]b shows the evolution of the frequency of binaries with high and low-mass primaries using the same symbols as in Fig. \[fig:normalized\_bin\_freq\_vs\_time\_\_masses\]a. Clearly the binary frequencies of the high and low-mass primaries evolve quite differently. The former decreases much more slowly with time than the latter and does not depend on the initial binary frequency. In contrast, the evolution of binaries with low-mass primaries *does* depend on the initial binary frequency.

The differences in the evolution are explained by two mechanisms. First, binaries with low-mass primaries are preferentially destroyed by dynamical interactions in the cluster because of their lower binding energy (Fig. \[fig:normalized\_bin\_freq\_vs\_time\_\_masses\]a). Secondly, the destruction of a binary by a high-mass primary can result in the generation of two new single high-mass stars or a high and a low-mass star. The latter possibility means that the number of single low-mass stars increases, hence the frequency of binaries with low-mass primaries decreases relative to the former case. Simultaneously, the frequency of binaries with a high-mass primary does not decrease as fast as implied by Eq. \[eq:evolution\_normalized\_number\_of\_binary\_frequency\] because of the smaller increase in the number of single high-mass stars. Hence, the possible generation of a high- and a low-mass star in the destruction of a binary with a high-mass primary accelerates the reduction of the binary frequency of low-mass stars and decelerates the reduction of the binary frequency of high-mass stars.

We emphasize the different evolution of binaries with high and low-mass primaries in Fig. \[fig:normalized\_bin\_freq\_vs\_time\_\_masses\]c. Here we show the difference between the binary frequencies of high and low mass stars over time, $\Delta_{\mathrm{hl}} \binaryfrequency (t)$, relative to its initial value $\Delta_{\mathrm{hl}} (t)$. This normalization accounts for the possible initial difference of the binary frequencies of high and low-mass stars, so that all values are zero initially.

In all simulations, the difference of the binary frequency of high and low-mass stars increases with time because the binary frequency of low-mass stars decreases faster than the binary frequency of high-mass stars (Fig. \[fig:normalized\_bin\_freq\_vs\_time\_\_masses\]b). However, the amount of increase depends strongly on the initial binary frequency in the cluster with the lowest initial binary frequency ($30\%$) resulting in the smallest increase of less than $2\%$ after 3 while in simulations with the maximum initial binary frequency ($100\%$) the difference increases by more than $10\%$. Under the assumption that initially the frequency of binaries with low- and high-mass stars was the same in the ONC, the observational data from imply that the binary frequency in the ONC initially was very high with values well above $75\%$. If the initial binary frequencies were lower initially the observed difference must have been part of the initial conditions of the binary population.

![Panel a: Evolution of the normalized number of binary systems with primary masses exceeding $2\Msun$ (thin lines) and binaries with masses below $2\Msun$ (thick lines) with time. Panel b: Evolution of the normalized binary frequencies of binaries with high- (thin lines) and low-mass (thick lines) primaries with time. Panel c: Evolution of the difference of the binary frequency of high- and low-mass stars with time. The difference at each time-step has been reduced by the initial difference of the binary frequencies to improve the readability. Crosses, stars, open, and filled squares denote simulations with initially 30%, 50%, 75%, and 100% binaries, respectively.[]{data-label="fig:normalized_bin_freq_vs_time__masses"}](multiplot_relative_filters_number_of_binaries_binary_frequency_deviation_of_binary_frequency_vs_time){width="\linewidth"}

Conclusions {#sec:conclusions}
===========

We have studied the dynamical evolution of the binary population in young dense ONC-like star-clusters using two different measures: the (normalized) binary frequency (Eq. \[eq:binary\_frequency\]) and the normalized number of binary systems (Eq. \[eq:nbin\_norm\]). As expected, the evolution of the binary frequency depends heavily on the initial binary frequency: the higher the initial binary frequency in a cluster, the stronger the decrease in the binary frequency during cluster evolution. However, our main result is that the evolution of the normalized number of binaries is independent of the initial binary frequency. During the first  of cluster evolution, the number of binary systems decreases by about $6-9\%$ independent of the initial binary frequency in the clusters.

The reason for this discrepancy between the two measures is the dependence of the binary frequency on the population of single stars. Since the destruction of a binary system leads to the generation of two new single stars, the single star population can easily dominate the evolution of the binary frequency in a cluster, especially in the case of high initial binary frequencies (Eq. \[eq:evolution\_normalized\_number\_of\_binary\_frequency\]). Hence, the number of binaries – unbiased by the evolution of the single star population – represents the proper measure of the evolution of the binary population. We can therefore conclude that the evolution of the *binary population* in ONC-like stars clusters is *independent* of the initial binary frequency, which is in striking contrast to what would be inferred from the evolution of the binary frequency.

This independence of the number of binaries has allowed us to determine the initial binary frequency and predict the evolution of the binary population in young ONC-like stars clusters from present-day observational data. For example, taking the observed binary frequency of the ONC of $\sim 60\%$ and its age of $1\Myr$ [@1998ApJ...500..825P; @1997AJ....113.1733H], and estimating $\alpha (1\Myr) = 0.92$ from our fit to the simulations, Eq. \[eq:initial\_binary\_frequency\] provides an initial binary frequency of $74\%$. This binary frequency is significantly lower than the binary frequency observed in sparse star-forming regions such as Taurus-Auriga and implies that the initial binary frequency depends on the star-forming region. This finding agrees well with the results of , who claimed that the primordial binary frequency in star-forming regions depends on the temperature of the natal molecular cloud in which the stars are forming of. However, different initial conditions of the clusters (e.g. higher initial stellar densities) could possibly increase the amount of destroyed binary systems, which would further increase the estimated initial binary frequency [@2001MNRAS.321..699K]. It is not trivial to quantify these effects, so further numerical studies are required.

In addition, our simulations show that the evolution of the binary populations depends crucially on the masses of the primary stars: the lower the mass of the primary, the more likely that the binaries are destroyed. Over time, this leads to a higher binary frequency for stars with masses exceeding $2\Msun$ than for lower mass stars. This is even so in the case where the binary frequencies of high and low-mass stars were initially the same. This finding is in good agreement with the observations of , who found a significantly higher binary frequency for high-mass stars than for low-mass stars in the ONC. Our results suggest that the difference in the binary frequencies of high- and low-mass stars are probably not the outcome of the binary formation process but rather the result of the dynamical evolution with time. A possible test of this hypothesis would be to compare the maxima of the period/semi-major axis distributions of high- and low-mass stars in young dense star clusters. Since the probability of destroying a binary in a three body encounter depends on the binding energy and therewith on the separation and mass of the binary, the maximum of the period and semi-major-axis distributions for massive stars should be at longer periods and wider semi-major axis than for lower mass stars. However, this test is not easy to be carried out since observations of the binary population in young dense star clusters usually only cover a small range of periods/separations of the binaries, and it is often not even possible to determine the physical orbital parameters based on the observations.

The simulations presented here were performed for ONC-like clusters only and used the specific binary population model of . We aim to generalise the present results by performing simulations with different cluster and binary population models for an forthcoming paper.

Simulations were performed at the Jülich Supercomputer Centre, Research Centre Jülich, Project HKU14.
