text:
# Accretion and Dynamical Interactions in Small-N Star-Forming Clusters: I. N=5 case

## 1 Introduction

The advent of infrared imaging has in recent years brought about the realization that virtually all stars are formed in clusters (e.g. Clarke et al. 2000). On a smaller scale, most field stars are known to be binaries (e.g. Duquennoy & Mayor 1991) or higher order multiples (e.g. Tokovinin & Smekhov 2002), and there is convincing evidence that the multiplicity among young stars is even larger (Reipurth & Zinnecker 1993; Ghez et al. 1993; Köhler & Leinert 1998). Even at earlier stages, Reipurth (2000) reports on a binary frequency as high as $``$ 80 % for sources driving giant Herbig-Haro flows. Additionally, the evidences for coevality of young binary systems (e.g. White & Ghez 2001) further suggest that clustering at the scales of molecular cloud cores must be a common property of star-forming regions. As shown by Sterzik & Durisen (1998), dynamical processes occurring within these small stellar groups at birth are likely to affect the properties of the resulting stellar objects. Similarly, gas dynamics are equally important at this stage: circumstellar discs are likely to influence the dynamics of close encounters (McDonald & Clarke 1995) and differential accretion from a common envelope can imprint a large range of dynamical masses on the stars involved (Bonnell et al. 1997).

Additionally, from a theoretical point of view, the fragmentation of a cloud core into a small aggregate of stars can help to attenuate the angular momentum problem in star formation (Mestel & Spitzer 1956; Bodenheimer 1995). As suggested by Larson (2002), the formation of a binary or multiple system can redistribute the excess angular momentum of the infalling envelope into the orbital motion of the stellar companions. Other mechanisms such as transport processes in discs and outflows might not be as efficient.

An often used argument against the paradigm of small N clustering at birth is based on the sub-millimetre results of Motte et al. (1998), who showed that the mass function of dense cores in Ophiucus is very similar to the stellar IMF. This important result has been widely interpreted as favouring a one to one mapping from core to star, leaving rather limited role for multiple star formation. Since direct evidence for such substructure is hard to address observationally (small N clusters disintegrate over a few crossing times), it is interesting to explore the consequences that this hypothesis can have on the initial mass function and the orbital parameters of binary stars, among other observables. The present paper will focus on small aggregates of N=5 protostars set in an initially uniform background of molecular gas, where the potential of the cluster is initially dominated by the gaseous component. Accretion and dynamical interactions take place naturally in this scenario and provide a diversity of mechanisms by means of which close binary stars, multiple systems and brown dwarfs can be formed, and their properties statistically compared with observational data.

The structure of the paper is as follows. In Section 2 the progress in our understanding of the fragmentation of gas clouds and the disintegration of multiple systems is reviewed, Section 3 describes the calculations performed, with a caveat on the role of circumstellar discs. Section 4 gives a detailed account of the evolution of the clusters. In Section 5 the results per given core mass are presented. Section 6 introduces the convolution with a core mass function. The properties of the brown dwarfs formed are discussed in Section 7. Finally, our conclusions are given in Section 8.

## 2 Previous Work

### 2.1 The disintegration of small N clusters

The earliest studies of the disintegration of small clusters treated only point mass dynamical interactions (van Albada 1968a,b). Subsequently, a number of authors have investigated the pure N-body problem using a variety of numerical codes (e.g. Harrington 1974, 1975; Mikkola and Valtonen 1986). In general, the outcome of such simulations is the break up of the cluster after $``$ N crossing times, the energy for the break up being released by the dynamical formation of a bound subsystem, usually a binary star containing the two most massive cluster members. As pointed out by McDonald & Clarke (1993), such an outcome has immediate consequences for binary pairing statistics, implying a predilection for binary membership amongst massive stars and a dynamical bias against brown dwarfs as binary companions. Another obvious implication of this formation mode is that the kinematics of young stars are likely to reflect the hardness of the three body encounters that ejected them from the cluster: for a cluster containing stars of typical mass M which suffer such an encounter at radius $`r_{close}`$, the typical ejection velocity is $`(GM/r_{\mathrm{close}})^{1/2}1.5r_{c100}^{1/2}\mathrm{kms}^1`$, where $`r_{c100}`$ is normalised to an encounter at $`100`$ A.U and a mass of 1M. Sterzik & Durisen (1995) were the first authors to link ejections from compact clusters to the generation of Runaway T Tauri stars (RATTs), whilst Armitage & Clarke (1997) pointed out that the ejection process would be likely to prune the circumstellar discs of such stars and thus render them less likely to be Classical T Tauri stars. To date, the most accurate and statistically complete study of the pure N-body disintegration problem is contained in the series of papers by Sterzik & Durisen and collaborators (Sterzik & Durisen 1995, 1998; Durisen et al. 2001), which analyse the decay channels and resulting kinematics and binary properties of a large ensemble of Nbody simulations. This set of simulations provides an invaluable benchmark against which to compare hydrodynamical simulations.

It is evident that in the case of young stars with discs, encounters at separations less than about 100 AU will not in any case be dominated by point mass gravity, but will be further modified by dissipative interactions with circumstellar discs. Clarke & Pringle (1991) and McDonald & Clarke (1995) modelled the effect of discs by adding parameterised drag terms during close encounters and found that the possibility of energy loss via the disc weakened the dynamical bias in binary pairing characteristics - in other words, lower mass stars ended up in binaries with a higher probability, especially as binary secondaries, compared with the results of pure Nbody simulations.

In all the above calculations, the mass of each star is assigned at the outset, so that although they yield predictions for binarity and stellar kinematics, they can say nothing about the form of the IMF. Such pre-assignment of masses corresponds to the case that most of the mass of the parent core is contained in non-linear density perturbations at the onset of gravitational instability. The opposite extreme would be a model in which gravitationally unstable fluctuations contain initially only a small fraction of the core’s mass, with the remainder being distributed throughout the core. In this latter case, the gravitationally unstable seeds acquire most of their mass through subsequent accretion, which is competitive in the sense that all seeds attempt to feed from the same reservoir. In such calculations, the final masses of stars are determined by the same dynamical processes that create the binary pairs and lead to stellar ejections. Such ideas were first explored in the simulations of Bonnell et al. (1997), who showed that due to the inequitable nature of the competitive accretion process, a large dynamic range of final stellar masses will result even in the case that the seeds all start off with equal masses. More recently, Reipurth & Clarke (2001) have specifically linked the cluster formation scenario to the production of brown dwarfs, arguing that the low mass end of the mass spectrum was populated by seeds that had lost out in the competitive accretion process and were prematurely ejected.

### 2.2 Hydrodynamics: The fragmentation of molecular cloud cores

The favoured mechanism for the formation of most binary stellar systems is the fragmentation of a collapsing molecular cloud. Fragmentation has been studied numerically for $``$ 20 years (Boss & Bodenheimer 1979; Boss 1986; Bonnell et al 1991; Nelson & Papaloizou 1993; Bonnell 1994; Burkert & Bodenheimer 1993; Bate, Bonnell & Price 1995; Truelove et al. 1998). These calculations appear to show that it is possible to form binaries with similar properties to those that are observed in fragmentation calculations. However, these simulations lack predictive power since the results depend sensitively on the initial conditions, which are poorly constrained. Besides, following the calculation significantly beyond the point at which fragmentation occurs is extremely computationally expensive. Bate & Bonnell (1997) quantified how the properties of a binary system are affected by the accretion of a small amount of gas from an infalling gaseous envelope. They found that the effects depend primarily on the specific angular momentum of the gas and the binary’s mass ratio. Generally, accretion of gas with low specific angular momentum decreases the mass ratio and separation of the binary, while accretion of gas with high specific angular momentum increases the separation and drives the mass ratio towards unity (see also Artymowicz 1983, Bate 2000). From these results, they predicted that closer binaries should have mass ratios that are biased toward equal masses compared to wider systems, since the gas which falls on to a closer system is likely to have more specific angular momentum, relative to the binary’s, than for a wider one.

A more direct approach has recently been conducted by Bate et al. (2002a,b), in which a gas cloud which is subject to a supersonic turbulent velocity field is modelled. The divergence-free velocity field rapidly generates a richly non-linear density structure in the gas (Ostriker, Stone & Gammie 2001). Dense cores are continually created and destroyed and, eventually, some of them become Jeans unstable and collapse (e.g. see Padoan & Nordlund 2002 for details). However, these cores possess internal sub-structure and fragment further until the opacity limit for fragmentation (a few Jupiter masses) is reached. Thereafter, the stellar seeds grow in mass by accretion and interact strongly with other collapsed fragments and protostellar discs in their vicinity. This system readily demonstrates the formation of small N ensembles in which the sort of behavior described in the last section is observed to occur. Currently, however, the accuracy of the statistical properties of the stars that can be determined from such a one-off calculation are limited by the small number of objects ($``$ 50) formed. Additionally, the computational expense involved in this large-scale simulations makes it prohibitive to explore different initial conditions. Modeling the evolution of many small-N ensembles allows us to obtain more accurate statistics and to investigate the dependence of the results on the initial conditions imposed.

## 3 The simulations

The calculations reported in this paper have been performed using a three-dimensional smoothed particle hydrodynamics (SPH) code based on a version originally developed by Benz (Benz 1990; Benz et al. 1990). The smoothing lengths of particles vary in time and space, such that the number of neighbours for each particle remains approximately constant at Nneigh = 50. We use the standard form of artificial viscosity (Monaghan & Gingold 1983) with strength parameters $`\alpha _v`$ = 1 and $`\beta _v`$ = 2.

The stellar seeds are modelled by sink particles (Bate, Bonnell & Price 1995), which interact only via gravity. Sink particles are characterized by a constant sink radius Rsink such that any gas matter which falls into it and is bound to the sink particle is accreted. The gravitational acceleration between two sink particles is Newtonian for r $``$ Rsink, but is smoothed within this radius using spline softening (Benz 1990). The maximum acceleration occurs at r $``$ Rsink/4; therefore this is the minimum binary separation Rmin that can be resolved. The sink radius has been chosen to be $`10^3`$ of the initial radius of the cloud, so that Rmin is $``$ 1 AU. The gas component is modelled using $`10^4`$ SPH particles with equal masses.

### 3.1 Initial conditions

One hundred simulations of these small N ensembles has been performed. Each simulation consists of an isothermal spherical gas cloud, with initial constant density. Five stellar seeds (sink particles) were randomly placed inside a sphere of radius R equal to one fourth of the initial cloud radius Rcl. The choice of N=5 is motivated by two reasons: first, we want to study the effect of dynamical interactions in unstable multiple systems (therefore N must be greater than 2) and second, we would like to allow for the possibility of binary-binary encounters to take place (thus, at least N=4). The precise choice of N=5, rather than 4 or 6 is not particularly relevant, since the results are not sensitive to small variations in the number of seeds. Results from clusters with a larger number of seeds, say N=10,20, will be provided in a future paper. Regarding the initial density distribution, the choice of a homogeneous density profile has been made following the properties of observed cores: e.g. Alves et al. (2001) found that the Bok globule B68 could be fitted remarkably well by a Bonnor-Ebert distribution, which is characterised by a flat inner profile; many other observations support this same result (see André et al. 2000 for a review). In addition, theoretical work by Bate (2000) has shown that cloud cores with initially constant density profiles are more likely to produce binary systems with properties more similar to those that are observed. A steeper initial density profile would enhance dramatically the effects of competitive accretion, consecuently producing stellar objects with too unequal masses (see Section 4 for details on how competitive accretion works.)

The initial mass of the cloud Mcl is a scale-free quantity and therefore is set to 1 M by default, so that the final stellar mass gives directly the cloud mass fraction that was accreted by that object. Obviously, the final real masses will be given by the product of that mass fraction times the mass of the core, taken randomly from a certain core mass probability distribution.

The gas component comprises 90 % of the initial mass, and is initially in virial equilibrium with its own gravitational potential. This condition determines the radius of the cloud to be equal to:

$$R_{\mathrm{cl}}=\frac{2}{5}GM_g\frac{\mu }{R_gT}$$
(1)
where $`G`$ is the gravitational constant, $`\mu `$ is the mean molecular weight (the gas is assumed to consist of pure molecular hydrogen and therefore $`\mu `$ is set to 2), $`R_g`$ is the gas constant, $`T`$ is the temperature, and $`M_g`$ is the gas mass.

Correspondingly, the protostars are also set in virial equilibrium with the gas. This condition determines the initial velocity dispersion of the seeds. Note that virial equilibrium means that, in most simulations, if only the gravitational force between seeds is computed, the majority of the protostars are bound to each other, i.e., they form initially non-hierarchical high-order multiples. The absolute value of the binding energy of the pairs, however, is very low, of the order of several $`\times `$ 10-5 times the binding energy of the resulting tightest pair. The cluster is bounded by a small external pressure to prevent the escape of gas particles at the initial stages when the gas thermal energy is comparable with its gravitational energy. The gas is initially stationary and thus the total angular momentum is equal to zero. The initial $`\beta `$ parameter (ratio of the rotational energy of the protostars to the potential energy of the gravitational interaction between the protostars and the gas cloud) of the stellar system has a mean value of 0.03, i.e., the initial angular momentum of the stellar seeds is not negligible but certainly low. Thus, we favour models in which the stars lack rotational support against the pull of the gas once this collapses towards the centre.

The numerical values of all the cluster fundamental parameters are given in Table 1.

An efficiency of 100 % is assumed in all calculations: gas is allowed to be accreted until all the cloud mass is in the form of stars. Although this is not realistic, the lack of firm constraints on star-formation efficiency of molecular cloud cores has deterred us from switching the gas accretion off at any arbitrary time, thus allowing the investigation of this process at all stages of the cluster evolution.

An unstable multiple is likely to remain in the cluster once all the gas has been exhausted. Thereafter, the calculation is followed using a pure N-body algorithm. This has been done for approximately 80 % of the simulations. The N-body algorithm used works under the same principles as the SPH code: adaptive time-steps evolved using a Runge-Kutta-Fehlberg integrator which provides a spatial resolution limited by the softening radius $`ϵ`$ $``$ Rsink . No regularization procedure has been applied.

Each of these simulations required an average of $``$ 450 CPU hours on Sun UltraSPARC 1/5/10 workstations for the hydrodynamical calculations. The N-body procedure needed an average of $``$ 350 CPU hours, the standard deviation of this distribution being much larger than in the hydrodynamical case due to the inequitable timescales involved in the random process of dynamical decay. To test the dependency of the results on the numerical resolution, calculations were also performed with 5 $`\times `$ 104 particles. The results did not show any qualitative difference from the 104 particles run.

These simulations do not include magnetic fields nor radiative transfer. An isothermal equation of state (T = 10 K) is used throughout the calculations to resemble the approximate balance between compressional heating and radiative cooling characteristic of the low-density regions of molecular clouds. No feedback from the stars has been modelled. This may be an adequate choice for low-mass star formation where strong stellar winds are not expected.

### 3.2 Circumstellar discs

Protostellar discs are an intrinsic byproduct of the star formation process. Disc formation occurs whenever any non-zero angular momentum is present in the initial gas from which the stars form. These discs can be quite large and can thus affect the dynamics of the stellar interactions (McDonald & Clarke 1995; Hall et al. 1996). The calculations reported here do not include any discs in the initial conditions. The reason for this is twofold: first, to place discs at the onset of the simulations would require a large number of SPH particles, which would subsequently slow down significantly the calculations and therefore make the choice of a statistical approach virtually impossible. Secondly, the initial conditions in the gas cloud do not provide enough vorticity to generate large protostellar discs; initial velocity perturbations such as turbulent fields would do this but, again, would require a computationally prohibitive number of gas particles. Global rotational motion in the cloud is known to generate vorticity and will be explored elsewhere.

The role of protostellar discs is to provide close encounters with an additional dissipative term which can affect the pairing characteristics of binary stars. Clarke & Pringle (1991) and McDonald & Clarke (1995) included discs in their model and found that the possibility of energy loss via disc weakened the dynamical bias in binary pairs – in other words, more lower mass stars ended up in binaries, especially as binary secondaries, compared with the results of pure N-body simulations. On the other hand, the inclusion of gaseous discs introduces a mechanism for the formation of stellar objects via fragmentation induced by gravitational instability (Bonnell 1994; Whitworth et al. 1995; Burkert, Bate & Bodenheimer 1997; Bate et al 2002a,b), process that could set a natural scale for the initial mean separation of some of the interacting fragments smaller than that assumed here. Additionally, circumbinary discs can have an important effect on the evolution of the binary orbital elements: the binary may transfer orbital angular momentum to the material in the disc via gravitational torques (Artymowicz et al. 1991; Bate & Bonnell 1997); this loss of angular momentum reduces the separation of the binary.

The effects of circumstellar disc on the dynamics of these small N clusters will be explored in the future.

## 4 The evolution of Small-N clusters

### 4.1 Onset of instabilities, cloud collapse and the formation of a binary system

As each calculation begins, the protostellar seeds start accreting gas from the common envelope. The initial gas density is constant and there are no pressure gradients to oppose the cluster self-gravity, which therefore begins to collapse. The stellar motions generate small perturbations in the gas density, producing local departures from the expected r-2 profile. Such perturbations arise from the accretion of gas that is in the vicinity of the stars: under-dense regions are created and, in turn, replenished by the infall of gas from larger distances. This replenishment is dependent not only on the gravitational attraction of the individual stars but mostly on the overall gravitational potential of the cluster.

The presence of gas in the cluster has two main effects on the evolution of the system’s energy budget (see Bonnell et al. 1997). First, there is mass loading, where the gas is accreted on to the individual stars, increasing their masses without changing their momenta and hence decreasing their kinetic energies. Second, the presence of gas acts as a drag source on the stars. The motion of the seeds through the cluster accelerates and condenses the gas behind the stars, which subsequently experiences isothermal shocks that radiate the additional kinetic energy away. More importantly, the cloud collapse proceeds isothermally and therefore the gas gravitational energy increases enormously (in absolute value), dragging in the stars, until enough accretion on to the seeds turns the stellar component into the main focus of gravitational attraction. These two processes contribute to making the cluster more bound, increasing the chances of forming very tight binaries.

Once the collapse of the cloud drives the stars towards the cluster centre, where the density is larger, accretion can proceed more efficiently. The accretion rate of a star depends predominantly on its position in the cluster. Those stars closest to the centre have the benefit of the cluster potential in attracting gas. The stars that, at the time of collapse, stay at the outside of the cluster never accrete appreciable amounts of gas. This differential accretion process will be decisive in providing a large range of final stellar masses. As soon as one object accretes a significant fraction of the gas reservoir, it starts to dominate the potential of the cluster and, therefore, turns into the main component of an unstable multiple stellar system. Several possibilities are thereafter allowed by the chaotic nature of these simulations.

First, two seeds may start close enough so as to soon form a binary, therefore enhancing accretion onto the system and out-weighting the effect of competitive accretion between the two components. At this early stages, the seed masses are still close to their initial values and therefore, the binary mass ratio is close to unity. Since no initial rotation is present in the cloud, the specific angular momentum of the infalling gas about the binary rotation axis will be in general lower than the specific angular momentum of the binary. This makes the binary mass ratio depart slowly from unity and the binary semi-major axis progressively decrease (see Section 2.2). Subsequent dynamical interactions can change significantly the semi-major axis of the binary, produce ejection of lower-mass seeds and break up the stellar system into a hierarchical multiple, but not change the binary mass ratio appreciably. The evolution of the separation and mass ratio of the two inner components of a binary and a triple system are shown in Fig. 1. Note how the semi-major axis of the binary changes abruptly when an ejection takes place, whilst the mass ratio decreases slowly due to the differential gas accretion of the two components.

Second, it may happen that only one object accretes a significant amount of gas before any close binary system is formed. In this situation, the other seeds are attracted towards the most massive member and acquire very eccentric orbits. The high eccentricity enhances the differential accretion as only at pericentre is the density high enough to make gas accretion efficient. The result of this process is to drive the mass ratio of the binary to small values, although never too far from unity (binaries with mass ratios smaller than 0.2 are rare).

Of course, most of the simulations performed display an accretion history that lies somewhere in between the two extreme cases described above. Yet the common outcome is that most of the cluster mass ends up in the binary system, as it settles in the gas-rich centre very early on in the cluster evolution, and that the mass ratio of this binary depends substantially on how soon the pairing occurs relative to the accretion history of the seeds.

In summary, the consequences of gas accretion are, among others, to produce an unstable multiple system frequently dominated by a binary star whose mass ratio has, to a significant degree, been set. From this point of the calculation onwards (with about 70 % of the gas matter left) dynamical interactions play a major role.

### 4.2 Dynamical interactions: ejections, binary hardening and hierarchical multiples

Dynamical interactions become important when the cluster evolution is close to its first free-fall time. The central binary has accreted a significant fraction of the available gas and begins to dominate the cluster potential. The other three seeds are attracted to its vicinity and very close encounters soon take place. The effect of these encounters is twofold: first, repeated three-body encounters between the central binary and the incoming object contribute to harden the binary, therefore decreasing abruptly its semi-major axis, by a factor ranging from 2 to 10, and increasing its eccentricity. The sudden growth of eccentricity can have a small long-term effect on the mass ratio evolution of the central binary, as explained in Section 4.1.

Secondly, the incoming object, with no exception a lower-mass member of the cluster, is ejected. Ejection means that either this object stabilizes into a hierarchical orbit (at a distance much larger than the new semi major axis of the binary), or it attains an outward speed larger than the cluster present escape velocity, therefore leaving the cluster forever. In the latter situation, the ejection halts any further accretion, therefore setting the final mass of the escaper and its kinematical properties. An object which remains bound to the cluster in a large orbit may also be deprived of further gas accretion if the gas is already too centrally-condensed.

In general, these simulations show that at least one object and sometimes two are ejected and escape from the cluster at an early stage, when about 70 % of the gas matter is still present. These escapers leave with a very low fraction of the cluster mass and, on average, need only a very soft kick to escape the potential well. The typical ejection velocities of these objects are $``$ 1 km s-1.

On the other hand, it is common that a non-hierarchical triple system survives the cluster break-up once all the gas has been exhausted. The decay of this unstable multiple is a stochastic process that needs to be followed numerically using pure N-body algorithms. Sometimes the triple system disintegrates into a hardened binary and a single. The single star, frequently more massive than those ejected at previous times, receives also a stronger kick as the distance of close approach can be as small as $``$ 1 AU (i.e. as low as it is allowed by the numerical resolution imposed). In other runs, the system stabilizes into a hierarchical triple. To consider when this is the case, the two-step criterion proposed by Eggleton & Kiseleva (1995) has been used. First, the triple system is partitioned into two pairs; this can be done in three different ways (e.g. \[1-2\]-3, \[1-3\]-2, \[2-3\]-1). For each such configuration, two instantaneously Keplerian orbits can be computed: the inner orbit of two of the bodies in their centre of mass (CM) frame, and the outer orbit of the third body around the CM of the inner component in the CM frame of the whole system. For the system to be a instantaneous hierarchical triple, both orbits must be bound (negative energy) at least in one of the configurations, and the ratio of the outer to inner periods must be greater than unity . If there is more than one configuration in which both orbits are bound, that with the largest ratio of outer to inner periods is selected. Second, if at a given time the triple system is found to fulfill the previous criterion, the ratio of the outer pericentre distance to the inner apocentre distance is computed. If this ratio is greater than a critical value Y$`{}_{}{}^{\mathrm{min}}{}_{0}{}^{}`$:

$$Y_0^{\mathrm{min}}1+\frac{3.7}{q_{\mathrm{out}}^{1/3}}+\frac{2.2}{1+q_{\mathrm{out}}^{1/3}}+\frac{1.4}{q_{\mathrm{in}}^{1/3}}\frac{q_{\mathrm{out}}^{1/3}1}{q_{\mathrm{out}}^{1/3}+1}$$
(2)
which depends on the inner $`q_{\mathrm{in}}`$ and outer $`q_{\mathrm{out}}`$ mass ratios, then the system is said to be stable for a time 10n $`\times `$ the triple outer period, where $`n`$ is larger than 4 in $``$ 95 % of the cases. This latter criterion guarantees that a substantial fraction of the triples included in our results have survival times greater than 100 Myr.

Finally, approximately 15 % of the simulations produce hierarchical quadruple systems, the only escaper being invariably ejected at very early times. Quadruples have been identified using a generalized version of the criterion described in the previous paragraph. In very few cases (2 %), two mutually unbound binaries are produced, plus one single star. One of these simulations is particularly striking as, contrary to the rule, the total mass in the two binaries is much smaller than that of the single.

## 5 Statistical results per core mass

### 5.1 Mass function

Table 2 shows the frequency of binary, triple and quadruple systems found in these simulations. The frequency $`f`$ is defined as the fraction of binaries/triples/quadruples over the total number of systems, where a system can be a single, a binary, etc…(e.g. if out of five stars, two are in a binary, the binary fraction would be $`\frac{1}{4}`$). 47 % of the simulations produce one binary star and three single stars. The binary mass ratio is in the range 0.5-1 and the three escapers are of much lower mass. In a minority of cases only two or one stars (36 % and 15 %, respectively) are ejected as singles, with the third (and fourth) being retained at large separations. Once again, the binary contains most of the mass, and the multiple companion(s) are, like the single stars, of very low mass. This pattern of mass acquisition is reflected in the mass distribution of stars produced for a given parent core mass (see Figure 2). This Stellar Mass Fraction Probability Distribution (hereafter SMFPD) is bimodal: one peak corresponding to the binary members and the other to the low mass singles and wide companions.

### 5.2 Binary mass ratio

Figure 3 depicts the dependence of binary mass ratio $`q`$ on primary mass. All systems have been decomposed in pairs, following the hierarchical distribution. Thus, e.g. a triple will contribute to the diagram as two binaries: first, the inner two bodies (diamonds), and second, the pair formed by the centre of mass particle of the inner binary and the tertiary companion (asterisk). Quadruples are represented by squares. Figure 3 illustrates two striking trends that can readily be related to the processes of gas accretion and dynamical interactions described in Section 4. Firstly, there is a relative paucity of binaries with low mass primaries – essentially there are no binary primaries with masses less than a factor of 0.5 the minimum core mass assumed. This is just because the fraction of the core mass that goes into the single stars/outer multiple companions is very small and so the binary primary contains a large fraction of the mass of its parent core. Secondly, there are very few extreme mass ratio systems and, where they exist, the low mass component is almost always the wide outermost companion.

The reason for this second trend is rather different from the reason that extreme mass ratio systems are also rare in pure N-body simulations. Here, the main effect is not so much that low mass companions are exchanged out of binaries by more massive ones but that, instead, once two seeds form a binary, which very often occurs early on in their accretion histories ($`q`$ therefore being not far from unity), the combined gravitational field enhances accretion on to both components, which subsequently accrete similar amounts of gas.

Likewise, the reason that the binary contains most of the core mass is not the result of repeated exchange interactions, but instead that, once two seeds form a binary, they constitute the focus for further accretion, and thereby outstrip the other cluster members in their mass acquisition.

Although these results are roughly compatible with the mass ratio distribution derived by Mazeh et al. (1992) for close G dwarf binaries, they clearly do not agree with the mass ratio distribution for G dwarf binaries as a whole: the results of Duquennoy & Mayor (1991) show a rise with decreasing mass ratio, at least for mass ratios as low as $``$ 0.2.

This deficit can be somewhat alleviated if some of the outer pairs in multiples systems are added to the distribution. As Figure 4 shows, the rise with increasing mass ratio characteristic of close binaries (solid and dot-dashed line) contrast with the steep rise in the distribution at very low mass ratio when all the pairs are included (e.g. dashed line for triples). This finding is in line with the general feature of the simulations: low mass companions are much more likely to be the outlier in triple systems than the secondaries in the central binary. A prediction of this result is that where wide low mass companions have been detected, their primaries might, on closer examination, turn out to be binaries (Gizis et al. 2001).

As mentioned in Section 3.2, the presence of circumstellar discs, neglected in the present simulations, might contribute to increase the fraction of low mass ratio systems by acting as a source of energy dissipation in stellar encounters, or by fragmenting into one companion when orbiting a single star. Future experiments (Delgado-Donate et al. in prep.) will focus on larger N ensembles in order to discover whether ejections due to binary-binary interactions can give additional insight into the formation of extreme mass ratio systems. Likewise, the evolution of N=2 systems (studied by Bate 2000 for circular systems), in which competitive accretion is the dominant process and there is no place for dynamical interactions, will be explored.

### 5.3 Semi-major axis distribution

Figure 5 depicts the semi-major axis dependence on primary mass. As it can readily be seen, most binaries (diamonds) are rather tight, with a mean semi-major axis $`a`$ 5 AU. On the other hand, wide systems (asterisk for triples, and squares for quadruples) correspond, with few exceptions, to triple and quadruple stars. This result can be immediately related to the pattern of dynamical decay outlined in Section 4.2. Given that a multiple system will only survive if its orbital elements adopt particular ranges of parameter space (so that it is sufficiently hierarchical to be dynamically stable; see e.g. Eggleton & Kiseleva 1995), triples and quadruples must display a semi-major axis distribution whose mean is $``$ 10 times larger than that of the binaries comprising the inner components.

The minimum semi-major axis for binaries is set by the numerical resolution achieved in these calculations, as parameterized by Rmin. On the other hand, the low mean value of the semi-major axis distribution is determined by the number and strength of those dynamical interactions that lead to an ejection. Therefore, $`a`$ is strongly anti-correlated with the number of seeds N=5 used throughout these calculations. The initial centrally concentrated spatial distribution of the protostars (the initial separation between seeds is $``$ 103 AU, 10 $`\times `$ smaller than the core radius), is also responsible of the formation of very close binaries. A more distributed initial configuration would hace produced a larger fraction of wide binaries.

From these results, together with those described in Section 5.2, it is clear that there exists a relation between the mass ratio, semi-major axis and multiplicity of the system. In this kind of simulations, wide binaries (semi-major axis $`>`$ 103 AU) could be formed if there is no orbital decay due to dynamical interactions. For this to happen the binary components must accrete mass in relative isolation and therefore, companions at larger separations are not likely to be present. Since there is no other mechanism that can reduce the orbit substantially (the collapse of the core induces a semi-major axis shrinkage of $``$ an order of magnitude), the two stars remain too distant to compete equally for the gas reservoir, consequently driving the mass ratio quickly to very low values.

### 5.4 Eccentricity

Among binary stars (i.e. excluding all outer pairs in multiples), the eccentricity can take a whole range of values from 0.1 to 1, as shown by Figure 6. The thick solid line histogram includes all the pairs, whereas the thin solid line refers to pure binaries, the dashed line to the inner binaries in multiples, and the dot-dashed line correspond to triple and quadruple systems. It is apparent that there is a preference for high eccentricity values, following a probability distribution with a functional form f(e) $``$ e. Although the eccentricity is the orbital parameter most sensitive to star-disc interactions (high eccentricities would be damped by the accretion of material from the disc), this result is in good agreement with the expected distribution of eccentricities for close young binaries (in which tidal circularization has not had time to operate yet; Ambartsumian 1937), and long period binaries (Duquennoy & Mayor 1991).

Multiple systems are characterized by larger eccentricities, as it is expected for objects which have undergone several close triple encounters before settling into a stable hierarchical configuration. This trend is even more accentuated in quadruples than in triples. Thus, one would expect a multiple system to display an increasingly larger eccentricity as one moves along the hierarchy, as well as a progressively smaller mass ratio. This pattern is reflected in the diagram shown in Figure 6: the eccentricity distribution of inner binaries (dashed line) is approximately flat whilst that of outer pairs (dot-dashed line) peaks at very high eccentricities.

### 5.5 Kinematics

The break-up of these small aggregates imprints a clear kinematic signature on the objects involved. Binary stars attain centre of mass velocities that, in average, are a factor of 10 smaller than the typical ejection speed of single stars (see Figure 7; singles are denoted by crosses and multiples as in previous figures; the dashed and dotted lines indicate the multiples and singles mean velocity, respectively). This distinction between the velocity distribution of singles and binaries is in essence due to the fact that single stars are all ejected from the core, whereas the binary (or multiple) star remains close to the core centre of mass. A notable difference between the N-body case (e.g. Sterzik & Durisen 1998) and these hydrodynamical simulations is the relationship between the typical ejection speeds and the initial parameters of the core. In the dissipationless calculations, the final velocities of ejected stars are of the same order as the virial velocity of the initial cluster. In the simulations reported here there are also a few objects which are ejected early on with such velocities, where the initial virial velocity of the core is $``$ 0.2 km s-1. However, most of the stars are ejected with velocities much greater than this: about half attain speeds greater than 5 times the initial virial velocity (the corresponding number for Sterzik & Durisen is less than 1 %), and there is also a significant minority that attain velocities $``$ 4 km s-1 (dot-dashed line in Figure 7). This may be readily understood inasmuch as in gas dynamical simulations the protostellar system shrinks as the gas is accreted, so that later interactions are much closer and produce correspondingly larger velocity kicks.

Another noteworthy feature of the kinematics of ejected objects is shown in Figure 8. The velocity of singles has been averaged in ten logarithmic bins (squares) spanning the whole range of masses (per given core mass). The best fit to these statistical averages does not display any significant correlation between ejection speed and escaper mass, within the error bars (given by Poisson noise). This result arises naturally from the dynamics of the accreting cluster. Low mass objects can be ejected at all stages of the cluster evolution since they can remain bound to the cluster at large distances from the gas-rich centre; thus, their speed pattern will reflect the different degrees of impulse needed to overcome the cluster potential as it gets stronger. This modus operandi is illustrated by the large spread in velocities shown in Figure 8. Note that the apparently lower dispersion in ejection velocities for higher mass singles (mass $`>`$ 0.06 M) results from the smaller number of objects located in this mass interval, as demonstrated by a Kolmogorov-Smirnov test applied to the speed distribution.

Consequently, the escaper’s speed depends mainly on the binding energy of the cluster at the time of ejection (see Bate, Bonnell & Bromm 2003 for a similar conclusion), which in turn would correlate roughly with the total mass in the stellar component at that time. This correlation can be seen in Figure 9.

It is unclear however whether this predicted difference in the velocity of young single and binary stars would in practice be detectable in star forming regions, since the velocity offset between the two (a few km s -1) is comparable with the core-core velocity dispersion in such regions. Additionally, the presence of various binaries in the same core might spread the velocity dispersion of the binaries via ejection of some of the pairs, therefore making it less distinct from the population of escapers (Bate, Bonnell & Bromm 2003).

## 6 The convolution

In order to relate the findings described in Section 5 to the resultant IMF, it is necessary to convolve the final stellar masses per given core mass (or SMFPD) with a core mass function (henceforth CMF). This means that the final masses of the protostars and related quantities are determined by two processes: first, the seeds gain mass according to the acquisition pattern found to operate in these cores; and second, the star obtains its final mass once the initial mass of its parent core is determined. In this operation we are assuming that the division of mass is a scale-free process. This would certainly be the case either if there were no systematic dependence on core mass of dimensionless variables such as e.g. the Jeans number (as controlled by equation (1)), or else if the results depicted here are insensitive to those variables.

Mathematically, this convolution differs from the usual one in the sense that two numbers (the stellar mass fraction of a certain seed and the mass of the core in which this star is embedded) are combined by a multiplicative instead of additive operation. That is, the final mass of a stellar object m will be given by the product r $`\times `$ Mc, where r is the mass accreted by the seed in units of a given core mass, and Mc is the mass of that core.

Unfortunately, the choice of a core mass function is not unique since the process that leads to the formation of dense star-forming cores is only partially understood. Several functional forms for this CMF are possible:

1. A power-law with index of $``$2, as expected from the equipartition of bulk kinetic energy and gravitational energy in a pressureless fluid.
2. A log-normal mass function, as derived from the density probability distribution function of Jeans unstable clumps in computer simulations of compressible magnetohydrodynamic turbulence (Padoan et al 1997; Klessen 2001).
3. A power-law with Salpeter or nearly-Salpeter slope ($`\alpha `$ = $``$2.35), as found by Motte et al. (1998) in the mass distribution of dense cores in Ophiucus.
4. A combination of power-law and log-normal distributions, as proposed by Padoan & Nordlund (2002) based on turbulent fragmentation calculations.

In order to fully understand the effects that a multiplicative convolution has on the CMF, we must note that a SMFPD described by a delta function would conserve the shape and slope of the CMF, and simply displace it to the position of the delta. If instead of a delta function we use a narrow box function, the CMF is additionally stretched to the new limits given by \[ r∗,min $`\times `$ Mc,min , r∗,max $`\times `$ Mc,max \]. Also, the convolution of two power-laws – e.g. a SMFPD of slope $`\alpha `$ confined between 0 and 1, and a CMF of slope $`\beta `$ with an arbitrary larger upper cutoff (where $`\alpha `$ is assumed to be smaller than $`\beta `$) – would, in general, give rise to a new function characterized by a break at a position given by the product of the lowest CMF value and the highest SMFPD value. This break would split the new function into two regions: at lower values, a power-law of slope $`\alpha `$, and at higher values, a power-law of slope $`\beta `$. Note that $`\alpha `$ and $`\beta `$ would exchange their range of influence if we choose $`\alpha `$ to be larger than $`\beta `$; that is, the steepest slope dominates at larger values.

### 6.1 Initial Mass Function

Figure 10 shows the result of a convolution of the bimodal SMFPD with a CMF taken to be a power-law of index $``$2.35 (Salpeter IMF) in the range of core masses 0.25 to 10 M. The choice of the power index and upper cutoffs are motivated by the millimetre continuum observations of Motte et al. (1998) and the H13CO+ survey of Onishi et al. (2002), respectively. The lower cutoff comes from the minimum mass a core can have so that the initial seed mass lies above the opacity limit for fragmentation ($``$ 5 Jupiter masses; Low & Lynden-Bell 1976; Rees 1976; Boss 1989; Bate 1998, 2002), given the initial conditions imposed – 10% of the mass is in stellar form at the start of the calculations. Clearly, the resulting IMF essentially follows the core mass function down to a stellar mass close to the minimum core mass. The reason for this correspondence is that for a steep core mass function, the majority of stars of a given mass belong to the high mass peak of low mass cores, rather than the low mass peak of the less numerous higher mass cores. Thus the number of stars of given mass is governed by the frequency of cores whose mass exceeds the stellar mass by a rather modest factor (less than a factor two) and the IMF therefore follows the core mass function.

Similarly, convolutions performed using log-normal mass functions as CMF, in which the characteristic mass and standard deviation adopt typical values, and the cutoffs are set at the same values chosen for the power-law CMF discussed above, produce invariably the same outcome. That is, the resulting IMF resembles the original log-normal mass function and simply differs in its width, which stretches to the values set by the products of minimum and maximum core and stellar masses. Also, as it may be expected from the previous results, the functional form proposed by Padoan & Nordlund (2002) for the CMF (in the original, IMF) is equally preserved.

Thus we find that although star formation in this case is far from being a one to one mapping between core and star mass, the IMF of the stars is controlled by, and nearly identical to, the core mass function. We note that the observed similarity between the IMF and core mass function cannot therefore be used as an argument against multiple fragmentation and competitive accretion inside star-forming cores.

### 6.2 Binary Properties

The convolution of core masses with a CMF has an effect on the resulting properties of binary stars. Although the binary mass ratio and eccentricity is independent of the individual masses of the binary components, the dependence of these parameters on the mass of the primary is modified by the CMF assumed. Additionally, if the initial mass of a core and its radius are assumed to be correlated in some way (as it is the case when e.g. the initial cloud conditions are set according to equation (1)), the semi-major axis distribution is also affected by the choice of the CMF. Equation (1) also states that the final velocity dispersion is independent of the initial mass assumed for each individual core.

Assuming (1) to hold, and making use of the same CMF as in Section 6.1, we have derived the scaled primary masses and semi-major axis of all the pairs. Fig 11 shows the semi-major axis dependence on primary mass. The linear relation between primary mass and semi-major axis imposed by equation (1) is reflected in the diagram, as well as the effect of the convolution, which has smeared out the distribution of points initially clustered around the original mean value. From this diagram, it can also be seen that the process of convolution with a CMF results in a large fraction of binaries with low mass primaries – binary primaries have a mass $``$ 0.5 $`\times `$ the core mass assumed, but there are many more cores with masses below 1 M than above.

The histogram of periods is shown in Fig 12. It can be fitted by an approximate Gaussian distribution which peaks at $``$ 103 days. Patience et al. (2000) find a similar value for binaries in dense environments (open clusters), whereas Duquennoy & Mayor (1991) show an observational maximum at $``$ 104 days for field G type stars (comparable to that of T Tauri stars). This discrepancy might be due to the number of stellar objects imposed as initial conditions in the model. Systems where repeated strong interactions are not so likely might produce longer period binaries – for this to be a solution to the lack of wide binaries in these simulations, the fraction of cores with two initial seeds must be comparable or larger than that of cores possesing more substructure. The addition to the gas of some initial angular momentum might as well enhance the fraction of binaries with longer periods.

In summary, the choice of an observationally-based CMF results in a set of binary masses and orbital parameters that are in reasonable agreement with observations. The formation of wide binaries, as suggested above, might be accomplished in a scenario where dynamical interactions were not important, as in cores where only two or three collapsed fragments are present. The origin of very low mass ratio binaries, however, remains obscure: both pure N-body calculations (e.g. Sterzik & Durisen 1998) and the present hydrodynamic-N-body models fail to produce as large a fraction of extreme mass ratio binaries as observational data suggest (the Duquennoy & Mayor mass ratio distribution peaks at a value of 0.2). Although some of these double stars might turn out to have internal substructure, in the sense of the inner body being in fact a binary star (e.g. Tokovinin & Smekhov 2002), – this would be immediately accounted for in our simulations –, it is obvious that very low mass ratio binaries constitute a significant population and its origin needs to be further investigated.

## 7 Brown Dwarfs

Once a convolution with a suitable CMF is performed (e.g. a power-law with nearly Salpeter slope within the observed mass range), the calculations show that the population of escapers and outer companions is abundant in brown dwarfs.

As described in Section 4.2, low mass objects are naturally produced when dynamical interactions within an unstable accreting multiple are taken into account. These interactions provide a simple dynamical explanation for the formation of sub-stellar objects, which start with a mass of the order of the opacity limit of fragmentation and are prevented from reaching the hydrogen burning limit due to the limited amount of gas they can accrete before they are ejected from the common envelope (Reipurth & Clarke 2001; Bate et al 2002a).

Figure 13 shows the mass ratio of multiple systems in terms of the mass of the primary component. As explained in Section 5.2, all multiples have been partitioned in two-body systems (symbol code as in previous figures). Below the solid line are located those pairs which posses a brown dwarf companion; binary brown dwarfs occupy the region to the left of the dashed line. It can be readily seen that brown dwarfs produced by dynamical interactions are infrequent as binary companions. On the other hand, among hierarchical triples and quadruples the outermost object is very often a brown dwarf, with a probability higher than 90 % for low-mass primaries – this fraction is certainly an upper limit, since some of the triple and quadruple systems included in this study might not be stable enough to survive the pre-main-sequence stage. In this sense, no brown dwarf desert at wide separations should be apparent if the inner components were not resolved (see Gizis et al. 2001). Brown dwarfs can be wide binary companions, but in this case the primary is predicted either to be a binary system or to have a low mass (less than $``$ 0.3 M). Binary brown dwarfs are equally improbable in this scenario, with the exception of some scarce, tight, nearly equal mass pairs.

The above results on brown dwarf binary pairings represent a critical area for confronting dynamical decay models with observations. It is currently unclear how the binary fraction changes in the case of very low mass primaries, since the detection of brown dwarf binaries is in its infancy (see, for example, Reid et al. 2001). If the binary fraction were indeed lower for brown dwarfs, it could readily be explained by such models, whereas a high incidence of brown dwarf binaries would be explicable only if the core mass function extended to sufficiently low masses (i.e. close to or below the hydrogen burning mass limit). Thus the detection of a large population of binary brown dwarfs would imply that most brown dwarfs arose as the majority component of low mass cores, rather than being the low mass ejectae from much more massive cores. Clearly, a firmer estimate of the brown dwarf binary fraction would place important constraints on formation models.

It is also worth noting that, in these models with minimum core mass of 0.25 M, the fraction of brown dwarfs that remain in bound systems is approximately 25% (75% are ejected as singles). This is in clear contrast to objects above the hydrogen burning limit, which display a boundedness fraction close to 90%.

The kinematics of sub-stellar objects are also likely to reflect their origin. Although there is no appreciable dependence of the final velocity of the singles on mass, brown dwarfs acquire as a class higher velocities than stellar objects, including binary and multiple systems (see Figure 14; symbols and lines as in Figure 7). With velocities of a few km s-1, it is possible that ejected brown dwarfs would be detectable in the vicinity of Class O objects, as suggested by Reipurth & Clarke (2001). Relatively few brown dwarfs have velocities in excess of $`5`$ km s-1, however, so that the bulk of brown dwarfs in the Pleiades would have been retained, given that the escape velocity in the young Pleiades is likely to have exceeded this value. On the other hand, smaller stellar associations such as Taurus-Auriga could have already lost a significant fraction of its sub-stellar population and therefore, would show a relative excess of binary stars. The findings of Briceño et al. (2002) are in line with this argument: their survey indicates that Taurus has $``$ 2 $`\times `$ fewer brown dwarfs at 0.02-0.08 M than the Trapezium, whereas it is well known (Ghez et al. 1993; Leinert et al. 1993; Simon et al. 1995, Köhler & Leinert 1998) that the binary frequency among T Tauri stars in Taurus is enhanced by a factor of two compared to solar type stars on the main sequence.

## 8 Conclusions

We have modelled the decay of non-hierarchical gas rich systems comprising 5 seed protostars in order to study the mass function, binarity and kinematics of the resulting population. This study can be seen as complementary to that of Bate, Bonnell & Bromm (2002a,b; 2003) who have modelled similar processes in the context of the calculation of the global evolution of a molecular cloud. Although the approach used here has the disadvantage that its initial conditions are artificial, and that it evidently cannot treat any possible interactions between star forming cores, the computational economy effected in not treating the whole cloud ensures a) that one can obtain statistically significant populations of stars for a given initial condition and hence isolate the factors that determine the properties of the resulting stars and b) that one can pursue the evolution to completion, i.e. to the point where all the gas has been accreted and where the system has decayed into a hierarchical multiple system. This latter feature allows us, in contrast to the calculations of Bate et al., to explore the properties of wider binaries, which remain undecayed in the Bate et al simulations.

Our chief conclusions are as follows:

1. The mass function of the resulting stars is primarily determined by the mass function of parent cores, being comparable to this core mass function down to stellar masses that are similar to the minimum core mass, and declining at masses below this. The reason for this is that the division of mass for given core mass is essentially bimodal - there is a high mass peak corresponding to a binary pair and a lower mass peak corresponding to ejected singles or wide tertiary companions. When this is convolved with a core mass function that declines steeply towards high masses, the majority of stars of given mass belong to the high mass peak of low mass cores rather than the low mass peak of rare, high mass cores. Consequently, the stellar mass function follows the core mass function. Thus the observed similarity between the core mass function and the stellar IMF (Motte et al. 1998) cannot of itself be used to disprove the hypothesis that most stars arise in small non-hierarchical multiples.
2. Stars much less massive than the minimum core mass are unlikely to be binary primaries. This can be understood in terms of the evolutionary pattern described above. Thus if a large fraction of brown dwarfs turn out to be binary primaries, this must imply that the parent core distribution extends well into the substellar regime. In this case, although brown dwarfs would form in non-hierarchical multiples, the brown dwarf binaries would not be ejected from these systems but would represent the highest mass objects in their parent cores. Thus a firm evaluation of the binary fraction amongst brown dwarfs will place important constraints on the properties of their parent cores. \[There is the additional possibility that star-disc interactions can produce tight brown dwarf-brown dwarf binaries that would survive ejection. The lack of discs in the simulations reported here do not allow us to model this process. The simulations of Bate et al. (which do include discs) indicate that this is not common however\].
3. Binary pairs with extreme mass ratios ($`q<0.2`$) are not produced in the simulations. More extreme mass ratios only occur in the case of low mass tertiary companions in a wide orbit around a central tight binary. Such a configuration is the natural consequence of the decay of non-hierarchical systems, in which low mass stars may either be ejected or else end up in wide hierarchical orbits. We point out that it is an expectation of dynamical decay models that where brown dwarfs are found as distant companions to normal stars (Gizis et al 2001; Kirkpatrick et al. 2001), the central object should normally be a close binary. Over all, the distribution of binary mass ratios (Figure 4) is well matched to the observations of Mazeh et al (1992) for close G dwarf pairs but is probably deficient in low mass companions when compared with the statistics of G dwarf binaries as a whole (Duquennoy and Mayor 1991).
4. Single stars are ejected with typical velocities of a few km s-1 (Figure 7) irrespective of their mass. The final velocities of binary stars with respect to their parent cores is typically an order of magnitude lower. In principle, this signature of dynamical decay could be sought in star forming regions, although may be masked by the relative motion of the parent cores (typically at a few km s-1). Since the binary fraction is lower among low mass stars and brown dwarfs than in higher mass stars (see (ii) above), this depedence of ejection speed on binarity means that in practice brown dwarfs are born with a higher velocity dispersion (with respect to their parent cores) than higher mass stars. With typical velocities of a few km s-1, the majority of brown dwarfs would be expected to be retained in open clusters.

## Acknowledgments

We thank the referee for comments that helped to improve the paper, and Jim Pringle, Cristiano Porciani and Bo Reipurth for their useful input. EJDD is grateful for support to the European Union Research Training Network The Formation and Evolution of Young Stellar Clusters.


format: text

text:
# THE DISK-JET CONNECTION IN MICROQUASARS AND AGN

## 1 Introduction

Collimated jets are associated with many astrophysical objects, ranging from active galactic nuclei (AGNs) to young stellar objects (YSOs). While the precise mechanisms for acceleration and collimation are still unknown, there is a growing consensus that in most, if not all cases, jet formation involves an accretion disk threaded by a large-scale magnetic field. In the case of YSOs, there is direct observational evidence (in the form of high-resolution images from HST) linking jets to the centers of accretion disks (e.g. Burrows et al. 1996). In the case of black hole X-ray transients (“microquasars”) the evidence is less direct, but is based on multi-wavelength observations (e.g. Mirabel & Rodriguez 1994).

In the present Letter, we propose a new interpretation of microquasar and AGN observations, in the context of a general model that connects disks and jets through the structure of the poloidal magnetic field.

## 2 Proposed Model

Belloni et al. (2000) have undertaken a thorough analysis of the X-ray behavior of GRS 1915+105 as seen by RXTE. They categorize the behavior by separating it into 12 different classes, and note that it would have been possible to even split these classes into many more. However, the main conclusion of the analysis is that the variability can be viewed as a set of transitions between three basic states. These states (denoted by A, B and C) are separated into three regions of the X-ray color-color diagram. State C represents a low luminosity and hard spectral index. This state corresponds to the low luminosity parts of the strongly variable behavior (the dips and quiescent phases), and can also be associated with about half of the time when the source is observed to be in a constant low state. The spectrum in State C can be modeled as being predominantly a power law with a low energy cut off due to absorption. There is also a thermalized contribution that can be interpreted as coming from an accretion disk, but, if so, one has to assume that the disk emission is truncated at around ten or twenty inner disk radii (Belloni et al. 1997a). States A and B show softer spectra and can be modeled as corresponding to thermalized emission from an accretion disk with an inner disk radius of around 20 km. State B is more luminous and hotter than State A, and can be interpreted as having a higher accretion rate.

Belloni et al. (1997a,b; 2000) interpret the transition between State C and States A/B as being caused by the disappearance and reappearance of the inner accretion disk due to a disk instability mechanism of the kind thought to occur in dwarf novae and other cataclysmic variables. If such an explanation is to hold, then the instability needs to be of the “inside-out” kind, in which a cooling front starts at the inner disk edge, progresses out to a radius of some tens of inner radii and then moves inwards again (Smak 1984). Thus most of the disk remains in the high state (on the upper branch of the putative $`S`$-curve in the effective temperature-surface density plane) and just the inner region jumps between the upper and lower branches in the standard limit cycle fashion. Consequently, during the low/hard state (C) the inner part of the accretion disk is either absent or has a much reduced viscosity, rendering it hard to detect. For the semi-regular behavior reported in the earlier papers of Belloni et al. (1997a,b) this appeared to be a reasonable hypothesis, and the time-scales of the behavior were found to be approximately right.

However, the wide variety in the nature of the variability of the source reported more recently (Belloni et al. 2000) renders this explanation more problematic. In particular, for the source to be in State C for an extended period of time it would be necessary for the transition front between the two stable parts of the $`S`$-curve to be stationary, which can be shown not to be possible (Smak 1984; Papaloizou & Pringle 1985). Instead, we propose here that during State C, the inner accretion disk is still present, and is still accreting at essentially the same rate as before, but, rather, the energy liberated in the accretion is converted efficiently into magnetic energy and is emitted in the form of a magnetically dominated outflow or jet.

In all models for the transport of angular momentum in accretion disks (i.e. for the viscosity) which involve magnetic fields, the first step in the process is to convert shear energy in the disk into magnetic energy, by using a poloidal magnetic flux to generate a toroidal component (Shakura & Sunyaev 1973; Eardley & Lightman 1975; Galeev, Rosner, & Vaiana 1979; Pudritz, 1981; Tout & Pringle 1992). What then happens to that magnetic energy is still rather uncertain. Some of the toroidal flux must be used as a source to generate a poloidal component and so to close the generation loop for the dynamo, generally assumed to occur through the magneto-rotational instability (Balbus & Hawley 1991). But most of the energy, which comes ultimately from the gravitational potential energy of the accreting material, must be lost from the disk. There are three basic possibilities, which are not mutually exclusive. The first, and the one usually assumed to occur, is that the magnetic energy is dissipated within the disk. This appears to happen in the numerical dynamo simulations (Balbus, Hawley, & Stone 1996) and in this case we have a standard accretion disk that emits thermalized radiation from its surface (Shakura & Sunyaev 1973). In this case we expect the disk to be radiating at essentially all radii. The second possibility, discussed initially in a number of contexts (e.g. Lynden-Bell 1969; Mészáros, Meyer, & Pringle 1977; Ionson & Kuperus 1984), is that the energy might be dissipated locally, but instead of within the disk, in the low density regions above and below the disk in the form of a chromosphere or corona. Whereas the thermalized emission from the disk is radiated at or around the local black body temperature, such a corona can reach temperatures comparable to the local virial temperature, and can radiate either thermally or by Compton scattering of low energy photons. Such a model is now regarded as the standard model for emission from AGN disks (e.g. Liang 1979; Haardt & Maraschi 1991) and for the low/hard state of some of the X-ray binaries (for example Cyg X-1; Done et al. 1992). In this case, however, unless the radiation process is collimated in some way (for example by Compton scattering of soft photons in an outwardly moving bulk flow; Malzac, Beloborodov, & Poutanen 2001) then we still expect to see at least some emission from the inner disk regions, for example Fe K$`\alpha `$ emission caused by irradiation of the disk by coronal heating. The third possibility is that most of the energy is not dissipated locally, but is released in the form of a bulk flow (as a mixture of kinetic energy and Poynting flux) either as a wind or as a well-collimated jet (Blandford & Payne 1982). In this case, it is conceivable that the inner disk radiates at a severely reduced level, if at all, with essentially all of the accretion energy being lost in a wind or jet. We should emphasize that it is not all the *mass* that is ejected into the jet (although some is), but most of the accretion *energy*.

## 3 The Two States

The idea that a large fraction of the accretion energy might be emitted from the disk in the form of an outflow or jet is not new. Blandford & Payne (1982) give an idealized calculation in which, in principle, all of the accretion energy can be emitted as a collimated outflow (see also Ogilvie & Livio 2001). Considerations of the models for the dynamo process within the disk have led some authors to suggest that the conversion of a substantial part of the accretion energy to a magnetically dominated outflow might be a standard possibility (Kato & Horiuchi 1986; Tout & Pringle 1992). There is also already some circumstantial evidence that the inner regions of some accretion disks in close binary systems are under-luminous. In particular, a comparison of the temperature profiles (Rutten, van Paradijs, & Tinbergen 1992), and of the mean ultraviolet ($`\lambda \lambda 120`$–200 nm) spectra (la Dous 1991) of the nova-like variables (which are thought to be steady standard accretion disks) with those of dwarf novae in outburst (quasi-steady disks of comparable brightness), suggests that the inner parts of the nova-like variable disks are under-luminous. HST eclipse observations of the nova-like variable UX UMa (Knigge et al. 1998) confirm this possibility. These observations, when taken together, give rise to the idea that the inner parts of such accretion disks can exist in (at least) two states. Most of the magnetic field generated by the dynamo process is small-scale, and tangled within the disk. If the field remains small-scale then most of the dissipation takes place within the disk and the disk is radiative. But if a larger-scale field is generated, perhaps by an inverse cascade above and below the disk (Tout & Pringle 1996) then the possibility arises of the more globally structured, poloidal field being able to drive an energetically significant outflow. It is important to note that while current MHD simulations are generally incapable of following in detail the process of large-scale poloidal field generation, the most recent numerical results (Turner, Bodenheimer & Rozyczka 1999; Kudoh, Matsumoto & Shibata 2002) give some credence to the jet production process we propose. The concept that an accretion disk in a nova-like variable can switch as time proceeds from being initially a standard radiative disk to being one in which it releases less of its energy thermally finds some support in the observations of the anti-dwarf nova TT Ari before and after a short minimum (Tout, Pringle, & la Dous 1993). Recent computations on driving outflows in the nova-like variables IX Vel and V3885 Sgr also indicate the need for magnetic driving (Hartley et al. 2002).

## 4 Application to GRS 1915+105: The Jet/Disk Connection

By analogy with the above considerations, we suggest that in GRS 1915+105 States A/B correspond to a state of the disk in which a dynamo-generated field is mainly small-scale with no global poloidal component. The accretion energy is liberated locally, some within the body of the disk to generate the thermal component of the spectrum, and some in a corona to produce the power-law component. We propose, however, that in State C, where the inner parts of the accretion disk appear to vanish, most of the accretion energy is not dissipated locally, but is emitted in the form of a relativistic jet.

If this is the case, then it should be possible to identify consequences of the jet emission. Fender (2001) argues that in addition to the ‘micro-quasar’ GRS 1915+105, four persistent (Cygnus X$``$1, GX 339$``$4, GRS 1758$``$258 and 1E 1740.7$``$2942) and three transient (GS 2023+38, GRO J0422+32 and GS 1354$``$64) candidate black hole X-ray binary systems display persistent flat or inverted radio spectra occurring during the extended periods in the low/hard state which corresponds to State C of Belloni et al. (2000). He suggests that this spectral component can be identified as synchrotron emission from a conical, partially absorbed jet (in part, by analogy with similar behavior in AGN; see, for example, Blandford & Königl 1979), and concludes that, if so, the jet power must be a substantial fraction of the accretion luminosity.

It might be objected that, in contrast to the AGN case, although some of the X-ray binaries show direct evidence of jet-like structures (typically on a scale of parsecs), most do not. However, for a jet to be visible it is necessary that it radiates, and for it to radiate it must interact, either internally, or with its environment. There are three basic zones in which such interactions can occur. The first is at the base of the jet, where the jet production mechanism is unlikely to be 100 percent efficient, and is likely to produce a measure of dissipation. This region is typically of a size of the order of a few Schwarzschild radii for a black hole accretor, if the jet is to be relativistic, and so too small to be directly observable. It is likely however that the dissipation in the jet production process is the cause of the compact self-absorbed synchrotron emission. The second zone is external material which the jet might run into. This interaction zone is clearly the cause of the extended radio emission seen in radio galaxies, but it seems plausible that most X-ray binaries are not surrounded by a sufficiently dense cocoon of material. The third interaction zone is internal to the jet itself and occurs when variability of the jet leads to internal shocks. It seems likely that this process is the cause of the observed correlation between X-ray variability and emission.

Mirabel et al. (1998) report on simultaneous X-ray, infra-red and radio observations of GRS 1919+105 during which the source oscillated between States B and C on a time-scale of about 30 mins (Class $`\beta `$ of Belloni et al. 2001). They associate the emission of IR/radio emitting plasma with the times of recovery from the X-ray dips; that is, with the times the source changes from a period in State C to one in which the behavior is one of rapid oscillation between States A/B and C. Fender et al. (1999) report on high-resolution MERLIN radio images of multiple ejections from GRS 1919+105. The superluminal proper motions of the radio blobs can be traced back to a time when the source had just emerged from a three-week extended phase in the low/hard X-ray state and was undergoing oscillations between States A/B and State C. Feroci et al. (1999) report on the detection of a large isolated radio flare observed with the Ryle Radio Telescope, which coincided with a period of oscillatory X-ray behavior seen by BeppoSAX. Similarly, Fender et al. (2002) again report simultaneous radio (RATAN$``$600 and ACTA) and X-ray (RXTE) observations in which large radio flares are seen to emerge at the end of week-long plateau (State C) X-ray phases, when again the X-ray behavior is one of rapid oscillation between States A/B and State C.

Taken together, these observational data can be interpreted as indicating the presence of a steady radio jet during the extended phases in State C, with abrupt, enhanced jet emission (flares and the ejection of blobs) being caused by enhanced variability of the jet emission mechanism as the disk changes in a stuttering fashion from putting most of its energy into a jet to dissipating most of its energy locally. In this context, we draw attention to simultaneous monitoring of the AGN 3C120 at both X-ray (RXTE) and radio (VLBA) wavelengths reported by Marscher et al. (2002). Superluminal knots are seen to emerge along the radio jet in this object. By tracing the blobs back in time to when they emerge from the nucleus, Marscher et al. conclude that each time the X-ray flux dips, and the X-ray spectral index becomes harder, a superluminal knot emerges at the site of the radio core.

## 5 Modeling the Variability

Our suggested model for the behavior of GRS 1919+105 depends on the ability of the central regions of the accretion disk to switch between two basic states. We speculate briefly in this section as to how and why this might come about. We note that a suggestion of a switching behavior which results in a variation of the jet speeds has been discussed by Meier et al. (1997). We should also point out that while GRS 1919+105 may be regarded as somewhat peculiar, when comparing with other black hole systems, it is also the source on which the most detailed data are available.

We note first that the launching of a jet probably requires a poloidal magnetic field of magnitude $`B_p`$ to exist over a scale of radius $`R`$, where $`R`$ is of order the inner disk radius, and where $`B_p`$ is much smaller than the field strengths $`B_{\mathrm{disk}}`$ associated with the dynamo-driven magnetic disk viscosity (Pringle 1993). This comes about by estimating the jet luminosity as being due to the local poloidal energy density moving outwards at the local dynamical velocity $`v_\varphi `$, i.e.

$$L_{\mathrm{jet}}B_p^2R^2v_\varphi ,$$
(1)
and the accretion disk luminosity (or the rate at which accretion energy is released in the disk) as

$$L_{\mathrm{acc}}\dot{M}_{\mathrm{acc}}v_\varphi ^2,$$
(2)
where $`\dot{M}_{\mathrm{acc}}`$ is the accretion rate.

Writing the accretion rate as $`\dot{M}_{\mathrm{acc}}\nu \rho _dH`$, where $`\nu \alpha c_sH`$ is the viscosity, $`\rho _d`$ the mean disk density, $`c_s`$ the disk sound speed, and $`H`$ the disk thickness (assumed $`R`$), and assuming for the dynamo viscosity that $`\alpha B_{\mathrm{disk}}^2/\rho _dc_s^2`$, we find that (Pringle 1993; Tout & Pringle 1996)

$$\frac{B_p^2}{B_{\mathrm{disk}}^2}\frac{L_{\mathrm{jet}}}{L_{\mathrm{acc}}}\frac{H}{R}.$$
(3)
This calculation is confirmed by Merloni & Fabian (2002). Thus, if the jet takes the bulk of the accretion luminosity (i.e. $`L_{\mathrm{jet}}L_{\mathrm{disk}}`$), we require $`B_pB_{\mathrm{disc}}(H/R)^{1/2}`$.

Numerical simulations of the dynamo process in accretion disks, while demonstrating the feasibility of the existence of a self-regulating dynamo operating within an accretion disk, are unable to undertake the kind of global computations required to simulate the generation of a global poloidal magnetic field (although see initial steps in this direction by Kudoh et al. 2002). Thus, for the time being we must resort to physical estimation rather than to detailed modeling. Using a simple model for an inverse cascade process, involving reconnection of magnetic loops above and below the disk, Tout & Pringle (1996) demonstrated that such a process might in principle lead to large-scale fields being generated by the dynamo.111It should be noted that previous suggestions that such a poloidal field might be generated by the dragging in of some global ambient field by the accretion disk are now generally discounted (Lubow, Papaloizou, & Pringle 1994a). The particular scaling derived by Tout & Pringle (1996) indicated that $`L_{\mathrm{jet}}/L_{\mathrm{acc}}H/R`$. However, Livio (1997) generalized their analysis and found that it is possible to have $`L_{\mathrm{jet}}L_{\mathrm{acc}}`$ provided that the distribution $`n(l)`$ of the lengths of the magnetic loops is of the form $`n(l)l^{3/2}`$. Thus, while these estimates indicate that the suggestions we put forward in this paper are not physically unreasonable, as mentioned above, calculations of the details of this process, including the time-scales on which changes might occur, are beyond current computing power.

Instead, we suggest a highly simplified picture of how the relevant time-scales might be estimated. We consider that the dynamo process operates locally, so that each annulus of the disk, of radial size $`H`$, can be considered to operate independently. This seems a reasonable conclusion to draw from numerical simulations of disk dynamos. This radial size also corresponds to the distance a magnetic field can diffuse in the disk using a typical time-scale for the dynamo cycle of $`t_{\mathrm{dyn}}(\mathrm{a}\mathrm{few})\mathrm{\Omega }^1`$, where $`\mathrm{\Omega }`$ is the local dynamical frequency, and an effective magnetic diffusivity comparable to the local effective viscosity (i.e. a local effective Prandtl number of order unity). If we assume that the poloidal field generated by the dynamo is typically of magnitude $`B_{\mathrm{disk}}`$, but of random sign, spread over of order $`mR/H`$ azimuthally distributed patches, then the distribution of net poloidal field from each annulus might be expected to have a mean value of zero, and a standard deviation of $`B_{\mathrm{disk}}/m^{1/2}B_{\mathrm{disk}}(H/R)^{1/2}`$. If this magnitude of poloidal field were reproduced over a range in radius of order $`R`$, then, as discussed above, this could give rise to a jet with $`L_{\mathrm{jet}}L_{\mathrm{acc}}`$. In order to achieve this, we require that $`NR/H`$ neighboring annuli all have a locally generated net poloidal contribution corresponding to of-order-one standard deviation in the same direction. Since the local dynamos vary on time-scales of $`t_{\mathrm{dyn}}`$, the time-scale on which this is likely to occur is comparable to tossing $`NR/H`$ coins once every $`t_{\mathrm{dyn}}`$, and asking how often all the coins come down either all heads or all tails. Thus we estimate the time-scale for establishing a change in the scale of the magnetic field in the disk to be of order

$$t_{\mathrm{jet}}t_{\mathrm{dyn}}2^{R/H}.$$
(4)
If we take $`t_{\mathrm{dyn}}0.01`$ s, corresponding to a couple of orbital periods at three Schwarzschild radii around a 10 M black hole, and we take $`H/R0.1`$ (Belloni et al. 1997a), then we obtain $`t_{\mathrm{jet}}10`$ s, which corresponds to the typical transition time-scale between States A/B and State C (Belloni et al. 1997a,b).

What determines the length of time the source spends in the various states (for example, how long the source remains in the low/hard states before switching) is much harder to estimate without a physical model. The observational data clearly require that both the high and the low states are reasonably stable for long periods of time, and so must correspond to equilibria of the system. Inspection of the various modes of variability categorized by Belloni et al. (2000) indicates that State B, that is the state in which there is no jet and in which all the dissipation occurs locally, is not very stable. All the light curves displayed in Belloni et al. (2000) show that whenever the source is predominantly in State B it switches to and from State C in a sometimes erratic and sometimes quasi-periodic fashion. In terms of our model, this implies that there is always a tendency for the disk to produce a jet. This tendency appear to be a common one among systems which contain accretion disks (Livio 2000). On the other hand, the system spends more than half of the time in State C and is able to remain in that state for prolonged periods of time. Thus State C must be a relatively stable equilibrium. For this to occur we require two conditions. First, once the dynamo process is able to establish a global poloidal field over a significant region of the disk, it is able to maintain that field for a period of time. Numerical simulations of the dynamo process do show significantly different behavior depending on whether or not a weak external poloidal field is applied, so there is at least the hope that this condition can be met. Demonstration of the effect will require global numerical simulations. Second, the poloidal field must be established in the center of the disk. This is a restatement of the effect, noted previously elsewhere (Pringle, 1993; Livio 1997), that the jet launching process from accretion disks cannot be scale-free and that jets are in general launched from the central disk regions.

In this regard we draw attention to two physical effects. The first is that a poloidal magnetic field can diffuse radially through a turbulent accretion disk on a time-scale $`H/R`$ times shorter than the radial flow velocity due to the magnetic viscosity (Lubow et al. 1994a). This is the reason (mentioned above) why it is generally not possible to drag an external ambient field inwards in an accretion disk. The second effect works in the opposite way: if the poloidal field gives rise to a jet/wind which carries away angular momentum at a rate comparable to or greater than that carried by the magnetic viscosity, then the inflow speed in the disk is significantly enhanced. Thus, for example, if in the central disk regions where the jet is launched the inflow speed caused by angular momentum loss in the outflow is increased over the speed due to the magnetic viscosity by a factor of $`R/H`$, then it becomes possible to trap the poloidal flux generated by the dynamo within those regions. This would enable to jet to be launched predominantly from the center of the disk. However, this same effect, also implies (Lovelace, Romanova, & Newman 1994; Lubow et al. 1994b) that if a poloidal field is established over a significant region in the outer part of the disk, to the extent that the outflow generated by it causes significant angular momentum loss, then an instability develops whereby those parts of the disk fall rapidly inwards. This probably implies that it is not possible to establish a steady jet from the outer parts of the disk. However, there is no reason why the dynamo process should not be able to establish a significant global poloidal field in the outer parts of the disk. This would occur on a longer time-scale than that in the inner regions, because the dynamical time-scales are longer at larger radii, and might be the cause of some of the longer term erratic variability. For example, if an outer disk region generates and advects into the center a poloidal field of opposite polarity (to that already in the center and generating a jet), then this might be a mechanism for a sudden switch-off of the jet and a sudden switch from State C to State B.

The switch in magnetic viscosity from radiative dissipation to driving outflow discussed here offers a possible explanation for the very long recurrence times of soft X-ray transients (SXTs). SXT outbursts themselves are convincingly explained as thermal-viscous disk instabilities, greatly prolonged by irradiation by the central X-rays, which trap the disk in the hot, high-viscosity state until most of its mass has been accreted (King & Ritter 1998; in the case of GRS 1915+105 the limit cycle may operate between the standard branch and the slim-disk branches, e.g. Honma, Matsumoto & Kato 1991). But it is currently unclear why the outbursts recur only at intervals of $`t_{\text{rec}}1`$–50 yr or more. In a standard picture where viscosity produces dissipation, the central disk regions would refill on a short viscous timescale, and would thus reach the critical surface density for triggering an outburst on much shorter timescales $`t_{\text{rec}}`$ weeks to months, like those of dwarf novae with comparable disk sizes (and thus similar binary periods). However, if instead, viscosity drives significant outflow, the central disk regions may remain much more depleted than in a dwarf nova, giving much longer recurrence times. A direct consequence of this idea is that quiescent SXTs may show significant radio emission. It appears that this is indeed the case for V404 Cygni (observations by R. Hjellming: M. Rupen, private communication). It would evidently be worthwhile checking other quiescent SXTs for radio emission. Very interestingly, Sokoloski & Kenyon (2002) point out that the symbiotic binary CH Cyg shows a radio jet in an optical low state, and indeed behavior closely analogous to that described above for GRS 1915+105. This supports the view (Livio 1997, 2002; Price, Pringle, & King, 2002) that there may be a common, accretion-powered formation mechanism for all jets. We also note that recent observations of the radio galaxy IC 4296 (Pellegrini et al. 2003), and of M 87 (Di Matteo et al. 2003), find that a sizable fraction of the accretion energy in these sources goes into powering the jet, in support of the scenario presented in this paper.