text:
# The Minimum Orbital Period in Thermal–Timescale Mass Transfer

## 1 Introduction

The importance of episodes of mass transfer in a semidetached binary on a thermal timescale has recently been emphasized in a number of contexts. Such episodes can arise in either of two ways:

(i) the donor star fills its Roche lobe while already undergoing thermal expansion across the Hertzsprung gap, or

(ii) the donor fills its Roche lobe with a mass ratio $`q=M_2/M_1`$ ($`M_2`$ denotes the donor mass, $`M_1`$ the accretor mass) large enough that the Roche lobe radius $`R_\mathrm{L}`$ shrinks on mass transfer more rapidly than the thermal equilibrium radius $`R_{\mathrm{te}}(M_2)`$ (where this exists) appropriate to the donor’s current mass $`M_2`$.

In case (i) the thermal–timescale episode will end once the donor attains a new thermal equilibrium radius, e.g. at some point after reaching the Hayashi line. In case (ii) the donor is continually trying to expand thermally beyond $`R_\mathrm{L}`$ in order to reach $`R_{\mathrm{te}}`$, so mass is transferred on a thermal timescale while this condition holds. Generally this case involves shrinkage of the orbit to some minimum Roche lobe size, followed by orbital expansion as the mass ratio $`q`$ reverses. The thermal–timescale episode ends only when $`R_{\mathrm{te}}<R_\mathrm{L}`$ (cf. Fig. 1; see also van den Heuvel 1992 and references therein). Thereafter normal mass transfer continues driven either by systemic angular momentum losses (descreasing $`R_\mathrm{L}`$) or nuclear evolution of the donor (increasing $`R_{\mathrm{te}}`$).

Case (i) arises in the formation of intermediate–mass X–ray binaries with black–hole accretors, such as GRO J1655–40 (Kolb et al., 1997; Kolb, 1998). (In similar systems with neutron–star accretors, such as Cyg X–2, the condition for case (ii) may hold simultaneously – King & Ritter, 1999; Podsiadlowski & Rappaport, 2000, Kolb et al., 2000; Tauris et al., 2000). In this paper we are mainly concerned with case (ii). This has received most attention in connection with the supersoft X–ray binaries (van den Heuvel et al., 1992). Thermal–timescale mass transfer from a donor initially on or close to the main sequence, on to a white dwarf accretor, offers a way of driving accretion rates $`\dot{M}_110^7\mathrm{M}_{}\mathrm{yr}^1`$ high enough to allow steady nuclear burning. As well as potentially explaining the observed supersoft systems, this process allows the white dwarf mass $`M_1`$ to grow. If one can arrange that $`M_1`$ reaches the Chandrasekhar mass $`M_C1.44\mathrm{M}_{}`$ this suggests a way of making Type Ia supernovae. However the difficulty of computing mass transfer on these timescales meant that early studies of this process simply used the assumption that mass transfer occurred on a thermal timescale, and were therefore unable to predict the evolution of the binary parameters (masses, period, mass transfer rate). Detailed calculations of these have only recently begun to emerge (Deutschmann, 1998) but are not yet exhaustive.

Discussion of the evolution of supersoft X–ray binaries would be greatly eased if observation provided reliable masses. However this is very difficult for several reasons. For example Greiner et al. (2000) suggest rather low masses $`M_10.6\mathrm{M}_{},M_20.35\mathrm{M}_{}`$ for the short–period system RX J0537.7–7034 we shall discuss extensively in this paper. But these and similar estimates use the assumption that the donor is close to the main sequence. By definition this cannot be true in thermal–timescale mass transfer, since the star is not in thermal equilibrium; although it may be quite close to its main–sequence radius, it can also be considerably smaller than this (cf. Deutschmann, 1998). Moreover Greiner et al.’s mass function uses the HeII emission line, and thus may not reflect the dynamical motion of the accretor.

Since mass information is so hard to come by, a crucial test for the thermal–timescale model for the supersoft X–ray binaries is provided by the discovery of systems with fairly short orbital periods. In particular RX J0537.7–7034 (Greiner et al., 2000) has a period of about 3.5 hr, and 1E 0035.4-7230 (= SMC 13) (Schmidtke et al., 1996) has a period of 4.126 hr. Since the initial mass ratio $`q_i`$ must be $`1`$ for this type of evolution (see below), one expects initial donor star masses $`M_{2i}1\mathrm{M}_{}`$ for a typical white dwarf accretor, and thus an initial orbital period $`P_i10`$ hr (see Section 3 and Fig. 3). Evidently considerable orbital shrinkage would be needed for such systems to reach the periods of RX J0537.7–7034 and 1E 0035.4-7230. This is unlikely if the mass transfer is conservative, as we shall show. However it is probable that much of the transferred mass is not accreted by the white dwarf, but blown away from it as a wind (e.g. Li & van den Heuvel, 1997), allowing greater orbital shrinkage. Mass loss induced in some way by the mass transfer process is probably the only way of significantly increasing the orbital shrinkage, as other angular momentum loss processes such as magnetic braking or gravitational radiation generally take place on timescales far longer than the $`10^7`$ yr characteristic of mass transfer in the supersoft X–ray binaries.

In this paper we consider thermal–timescale mass transfer and investigate how much mass loss from the accretor is required if RX J0537.7–7034 and 1E 0035.4-7230 are products of the standard picture of the supersoft binaries. Our method is to compute the minimum orbital period analytically for specified rates of mass and angular momentum loss from the binary.

## 2 Orbital Evolution

We consider the orbital evolution of a semidetached binary in which a fixed fraction $`1\eta `$ of the mass transferred from the donor (star 2) is lost from the accretor (star 1) with $`\beta `$ times the specific angular momentum of the latter. (The quantity $`1\eta `$ is called $`\alpha `$ by King & Kolb, 1995 in their general treatment of such ‘consequential angular momentum loss’ \[CAML\] mechanisms; the use of $`\eta `$ allows more compact formulae in what follows.) Following the general method of van Teeseling & King (1998), setting their quantities $`\dot{M}_{w2}=\dot{J}_{\mathrm{sys}}=0,\beta _1=\beta `$, we find

$$\dot{M}_{w1}=(1\eta )\dot{M}_2,$$
(1)
so that

$$\dot{M}_1=\dot{M}_2+\dot{M}_{w1}=\eta \dot{M}_2.$$
(2)
These definitions give

$$\frac{\dot{J}}{J}=(1\eta )\dot{M}_2\beta \frac{M_2}{M_1M},$$
(3)
with $`M=M_1+M_2`$. Kepler’s 3rd law links angular momentum $`J`$ and period $`P`$ of the orbit as

$$J^3=\mathrm{G}^2\frac{M_1^3M_2^3}{M}\frac{P}{2\pi }.$$
(4)
The evolution of the orbital period then obeys

$$\frac{\dot{P}}{P}=\frac{3\dot{M}_2}{M_2}+\frac{3\dot{M}_2}{M_1}\frac{2(1\eta )\dot{M}_2}{M}+3(1\eta )(\beta 1)\frac{M_2\dot{M}_2}{MM_1}.$$
(5)
From (2) we have $`M_1=C\eta M_2,M=C+(1\eta )M_2`$, with $`C=M_{1i}+\eta M_{2i}`$, where $`M_{1i},M_{2i}`$ are the initial values of the two masses. Thus (5) integrates to

$$PM_2^3M_1^{3(1\beta )3\beta /\eta }M^{13\beta }$$
(6)
for $`\eta 0,1`$, the corresponding expression for $`\eta =0`$ being given by the limit of this expression as $`\eta 0`$, while we get the well–known result

$$PM_2^3M_1^3,$$
(7)
for conservative mass transfer, i.e. $`\eta =1`$. In terms of the mass ratio $`q=M_2/M_1`$ we can express the relation between $`M_2`$ and $`M_1`$ as

$$qM_1=M_{2i}+\frac{1}{\eta }(M_{1i}M_1),$$
(8)
leading to

$$\frac{M_{1i}}{M_1}=\frac{q\eta +1}{q_i\eta +1}.$$
(9)
Thus with $`M=(1+q)M_1`$, etc, (6) gives finally

$$\frac{P}{P_i}=\left(\frac{q_i}{q}\right)^3\left(\frac{1+q_i}{1+q}\right)^{(3\beta 1)}\left(\frac{1+q\eta }{1+q_i\eta }\right)^{5+3\beta /\eta },$$
(10)
where $`P_i`$ is the initial period.

As a check we note that if all the transferred mass is blown away from star 1, with the specific angular momentum of that star, we have $`\eta =0,\beta =1`$; we take the limit of (10) as $`\eta 0`$ by noting that $`(1+q\eta )^{3/\eta }=\mathrm{exp}[(3/\eta )\mathrm{ln}(1+q\eta )]`$ and using l’Hôpital’s rule on the exponent to give the limit as $`\eta 0`$ as $`e^{3q}`$. We then get

$$\frac{P}{P_i}=\left(\frac{q_i}{q}\right)^3\left(\frac{1+q_i}{1+q}\right)^2e^{3(qq_i)},$$
(11)
or $`PM_2^3M^2\mathrm{exp}[3M_2/M_1]`$, as found by e.g. King & Ritter (1999) for this case.

We are interested in the minimum orbital period attained during the binary evolution. Regarding $`P`$ in (10) as a function of $`q`$, i.e. with $`P_i,q_i`$ fixed, we find

$$\frac{1}{P}\frac{P}{q}=\frac{3}{q}\frac{3\beta 1}{1+q}+\frac{3\beta +5\eta }{1+q\eta },$$
(12)
so that $`P`$ is an extremum at

$$q=q_m=\frac{(1\eta )+\sqrt{(1\eta )^2+9(\eta +\beta \eta \beta )}}{3(\eta +\beta \eta \beta )}$$
(13)
From (12) we can show that

$$\frac{1}{P}\frac{^2P}{q^2}=\frac{3}{q^2(1+q\eta )}+\frac{(3\beta 1)(1\eta )}{(1+q)^2(1+q\eta )}$$
(14)
at $`q=q_m`$, so $`P`$ has a minimum there, assuming $`0\eta 1`$. Hence for given $`P_i,q_i`$, the smallest orbital period which can be attained in evolution with mass loss is given by (10) with $`q=q_m`$. For a given initial period, we still have the freedom to vary the initial mass ratio $`q_i`$. (This corresponds to the fact that the period $`P_i`$ essentially determines the donor’s mean density, almost independently of the primary mass. So a donor of given initial mass $`M_{2i}`$ may initiate mass transfer at various $`P_i`$, depending on its evolutionary state, with different $`q_i`$, determined by $`M_{1i}`$.) For fixed $`P_i`$, the minimum value of $`P`$ clearly has a maximum, regarded as a function of $`q_i`$, at $`q_i=q_m`$ (see also Fig. 2 which shows the special case $`\beta =1`$). In other words, $`\mathrm{d}(P(q=q_m))/(\mathrm{d}q_i)_{[P_i\mathrm{fixed}]}<0`$. Since $`q`$ decreases from its initial value $`q_i`$ through the evolution this means that given an initial period $`P_i`$, the smallest value of the minimum period is given by the largest possible value $`q_{il}`$ of $`q_i`$. Hence given an initial period $`P_i`$, the minimum possible orbital period $`P_m(P_i)`$ is given by

$$\frac{P_m}{P_i}=\left(\frac{q_{il}}{q_m}\right)^3\left(\frac{1+q_{il}}{1+q_m}\right)^{3\beta 1}\left(\frac{1+q_m\eta }{1+q_{il}\eta }\right)^{5+3\beta /\eta },$$
(15)
with $`q_m`$ given by (13), and $`q_{il}`$ the largest possible value of the initial mass ratio $`q_i`$. Although $`\beta `$ is of course irrelevant in the conservative case ($`\eta =1`$), larger angular momentum losses in the wind allow even smaller ratios $`P_m/P_i`$ for the same $`q_i`$.

It is likely that thermal–timescale mass transfer is still going on at the minimum orbital period, since by Roche geometry we have

$$R_\mathrm{L}f(q)M^{1/3}P^{2/3}.$$
(16)
where $`f(q)`$ is a slowly increasing function of $`q`$ given e.g. by the approximation from Eggleton (1983). Thus

$$\frac{\dot{R}_\mathrm{L}}{R_\mathrm{L}}=\frac{\dot{q}}{f}\frac{\mathrm{d}f}{\mathrm{d}q}+\frac{\dot{M}}{3M}+\frac{2\dot{P}}{3P},$$
(17)
which tells us immediately that the period reaches its minimum value at larger mass ratios than $`R_\mathrm{L}`$ as $`\dot{q}<0,\dot{M}0`$. In general $`R_\mathrm{L}`$ initially shrinks much more rapidly than $`R_{\mathrm{te}}`$. Thus the condition $`R=R_{\mathrm{te}}`$ signalling the end of thermal–timescale mass transfer is reached only once $`R_\mathrm{L}`$ is increasing after passing through its minimum value, i.e. thermal timescale mass transfer usually stops only after the minimum orbital period is reached. However it is conceivable that in some cases $`R_{\mathrm{te}}`$ shrinks very rapidly on mass loss, and individual cases must be checked. For the purposes of this paper we note than the true minimum period for thermal–timescale could if anything be longer than the values of $`P_m`$ we find here.

In what follows we will assume that $`\beta =1`$, i.e. that the mass lost from the accretor has the same specific angular momentum as this star. This for example includes any form of mass loss from orbits with circular symmetry about the accretor, as might occur from an accretion disc. In this case (15) simplifies to

$$\frac{P_m}{P_i}=\left(\frac{q_{il}}{q_m}\right)^3\left(\frac{1+q_{il}}{1+q_m}\right)^2\left(\frac{1+q_m\eta }{1+q_{il}\eta }\right)^{5+3/\eta },$$
(18)
with $`1q_m1.387`$ now given by

$$q_m=\frac{(1\eta )+\sqrt{(1\eta )^2+9}}{3}.$$
(19)
and $`q_{il}`$ the largest possible value of the initial mass ratio $`q_i`$. Figure 2 shows the ratio $`P_m/P_i`$ as a function of $`q_i`$ for the two extreme cases $`\eta =1,0`$.

## 3 Minimum Periods

The work of the last Section, especially eqn (15), shows that to compute the ratio $`P_m/P_i`$ of minimum to initial period we need to specify the largest possible value $`q_{il}`$ of the initial mass ratio. For a white–dwarf accretor the short observed periods of RX J0537.7–7034 and 1E 0035.4-7230 then constrain the mass–loss parameter $`\eta `$. Clearly to reach such short periods it is preferable to start from the shortest initial periods $`P_i`$, which in turn are given by assuming that the donor is still very close to the ZAMS when it initiates mass transfer. We can easily show that the initial donor mass $`M_{2i}`$ must be $`1.9\mathrm{M}_{}`$ in order to explain the 3.5 hr period of RX J0537.7–7034, by iterating the formula (18), using the fact that we must clearly have $`q_{il}>q_m`$. For the lowest–mass white dwarf we consider ($`M_{1i}=0.7\mathrm{M}_{}`$) this implies $`M_{2i}>0.7q_m\mathrm{M}_{}`$. Thus in the conservative case $`\eta =1`$ we have $`q_m=1,M_{2i}>0.7\mathrm{M}_{}`$, and Fig. 3 shows that $`P_i6`$ hr for a ZAMS star of this mass. We can now use Fig. 2 with the restriction $`P_m/P_i<3.6/6=0.58`$ to find $`q_{il}>2.36`$ and thus $`M_{2i}>1.66\mathrm{M}_{}`$. Figure 3 now shows that $`P_i10`$ hr, and we may iterate using Fig. 2 to get $`q_{il}>3.38,M_{2i}>2.37\mathrm{M}_{}`$. Further iteration fails, as Figure 3 shows that the new higher estimate for $`M_{2i}`$ does not increase the estimate for $`P_i`$. In a similar way we find $`q_{il}>2.7,M_{2i}>1.89\mathrm{M}_{}`$ for the other extreme case $`\eta =0`$.

This argument essentially fixes the minimum value of $`P_i`$ at about 10 hr (see Figure 3), so the orbital shrinkage required to explain the 3.5 hr period of RX J0537.7–7034 has $`P_m/P_i<0.35`$.

Given an initial white dwarf mass $`M_{1i}`$ and mass ratio $`q_{il}`$ we can predict the period ratio $`P_m/P_i`$ for any given $`\eta `$. The initial mass ratio $`q_{il}`$ must obey two constraints (added as vertical lines in Fig. 2):

(i) the system must avoid the ‘delayed dynamical instability’ (DDI; see Webbink, 1977; Hjellming, 1989), which occurs when sustained thermal–timescale mass transfer exposes inner layers with a flat entropy gradient, and

(ii) the white dwarf mass cannot exceed $`M_C`$ before the minimum period $`P_m`$ is reached.

In practice the first of these constraints requires $`q_{il}<q_{\mathrm{DDI}}3`$. Hjellming (1989) found this value for a donor with $`M_{2i}=3\mathrm{M}_{}`$ near the terminal main sequence (TMS), whereas Kalogera & Webbink (1996) seem to prefer even smaller limiting ratios around $`2.5`$. Kolb et al. (2000) used Mazzitelli’s stellar code to calculate a test sequence with constant primary mass $`0.75\mathrm{M}_{}`$ where mass transfer starts from a $`3\mathrm{M}_{}`$ near–TMS star. A DDI occured at donor mass $`2.6\mathrm{M}_{}`$, in perfect agreement with Hjellming’s prediction. The critical maximum mass ratio $`q_{\mathrm{DDI}}`$ depends on the stellar structure, and therefore on the stellar input physics. In particular, $`q_{\mathrm{DDI}}`$ is probably sensitive to the degree of convective core overshooting, as this determines the size of the convective core. This is highlighted by the fact that Kolb et al. (2000) find $`q_{\mathrm{DDI}}=2.9`$ for early massive case B mass transfer, while Tauris et al. (2000) find $`q_{\mathrm{DDI}}=3.6`$, using an updated Eggleton code.

The second constraint is quite severe for large initial white dwarf masses and for $`\eta 1`$. In particular, for $`M_{1i}=0.7\mathrm{M}_{}`$ it requires $`q_{il}<3.11`$, in mild contradiction with the requirement $`q_{il}>3.38`$ we found above.

Figure 4 shows $`P_m/P_i`$ versus $`\eta `$ for various values of $`q_{il}`$ The stars denote combinations $`M_{1i},q_{il},\eta `$ where $`M_1`$ reaches $`M_C`$ precisely at $`q=q_m`$. To the right of these positions we take $`P_m`$ as the minimum period actually achieved, i.e. the period where the white dwarf reaches $`M_C`$.

Each panel of this Figure shows:

1. The horizontal line $`P_m/P_i=0.35`$, i.e. the upper limit required by RX J0537.7–7034.

2. Only those curves which actually manage to cross this line for $`\eta 0`$. In particular we plot the curve for the limiting value of $`q_i`$ such that the curve just crosses $`P_m/P_i=0.35`$ at $`\eta =0`$. This gives a lower limit on $`q_i`$ which is the same for each mass of the WD.

3. At some of the extreme solutions, i.e. wherever the line $`P_m/P_i=0.35`$ is crossed before a SN occurs, the current values of $`M_1,M_2`$ are given. Note that with growing $`M_{1i}`$ the donor mass has to go up as well for the same $`q_i`$. Furthermore $`M_1`$ at the minimum period comes closer and closer to the Chandrasekhar mass as we increase $`\eta `$.

## 4 Discussion

Figure 4 shows already that the short orbital period of RX J0537.7–7034 poses very severe constraints if this system results from thermal–timescale mass transfer. In particular

(a) Conservative evolution ($`\eta =1`$) is possible only for initial mass ratios $`q_i>3`$ which probably make the system vulnerable to the delayed dynamical instability. Full evolutionary calculations are required to check if there are any evolutionary tracks which avoid it. If such tracks exist, the initial white dwarf mass in RX J0537.7–7034 must have been low ($`0.7\mathrm{M}_{}`$), but the current system masses must be fairly high (e.g. $`q_i=4`$ requires $`M_1=1.17\mathrm{M}_{},M_2=2.33\mathrm{M}_{}`$). These are of course in conflict with the mass estimates of Greiner et al. (2000), but this may not by itself be fatal (see the remarks in the Introduction).

(b) Non–conservative evolution ($`\eta <1`$) does allow tracks with $`q_i3`$ which probably avoid the delayed dynamical instability. However the ranges of $`\eta `$ and $`q_i`$ are still very tightly constrained. Figure 4 shows that $`2.7q_i3`$, $`0<\eta 0.3`$, with $`\eta `$ and $`q_i`$ correlated as in Figure 5. For $`M_{1i}1.0\mathrm{M}_{}`$ $`\eta `$ has to be even lower, and generally upper limits on $`q_i`$ are given for each $`M_{1i}`$ in addition to DDI. Note that for the largest values of $`\eta `$ the white dwarf in RX J0537.7–7034 is predicted to be close to the Chandrasekhar limit unless the DDI limit imposes a more severe restriction than the lowest possible WD mass (as is the case in Fig. 5 for $`q_i3,M_{1i}0.7\mathrm{M}_{}`$).

Even given these tight constraints, the evolutions discussed above require the system to have come into contact with the donor still very close to the ZAMS. As orbital angular momentum losses (e.g. via magnetic braking) are probably negligible for the likely initial masses $`M_{2i}2\mathrm{M}_{}`$, this also requires the initial orbital separation to lie in an extremely narrow range. Evidently to make a system like RX J0537.7–7034 or 1E 0035.4-7230 by the thermal–timescale route is a very rare event. In line with this, Deutschmann (1998) found no orbital periods shorter than about 6 hr in his detailed calculations with solar metallicity. Furthermore these requirements have been derived assuming that RX J0537.7–7034 and 1E 0035.4-7230 are just at the minimum period $`P_m`$, so the conditions might be even harder to meet. The probability of observing such a system is larger near the end of the thermally unstable phase because the mass transfer rate decreases (in both Deutschmann’s and our own full computations).

We have so far neglected the effect of tidal interactions on the orbital evolution of the binary. Tauris & Savonije (2000) show that these can be important in low–mass X–ray binaries by translating spin angular momentum losses (via e.g. magnetic stellar wind braking) into orbital losses when tidal synchronization occurs. However for the binaries we consider here, any effect before the beginning of mass transfer simply allows a shorter $`P_i`$ for a wider range of systems. Once synchronism is achieved, the angular momentum of a lobe–filling donor is less than about 10 per cent of the orbital angular momentum (because the gyration radii of typical donor stars is small, $`r_\mathrm{g}^20.10\mathrm{}0.20`$). Even transferring all of this to the orbit in the course of the evolution would lead to only a marginal shift towards longer minimum periods, leaving our conclusions unchanged.

The work of this paper suggests that while the thermal–timescale mass transfer model for the supersoft X–ray binaries has many desirable features, it may not be possible to use it to describe all of the supersoft binaries, as well as SNe Ia progenitors. For example, if highly non–conservative mass transfer is as common as seems to be required to explain RX J0537.7–7034 or 1E 0035.4-7230, this would make building up the white dwarf mass to produce a Type Ia supernova highly problematical. It may therefore be necessary to consider other possibilities (e.g. van Teeseling & King, 1998; King & van Teeseling, 1998).

## 5 Acknowledgment

ARK gratefully acknowledges the support of a PPARC Senior Fellowship. This work was partially supported by a PPARC Rolling Grant for Theoretical Astrophysics.


format: text

text:
# The Formation of a Bound Star Cluster: From the Orion Nebula Cluster to the Pleiades

## 1 INTRODUCTION

The Orion Nebula Cluster (ONC) is a rather peculiar object. It is extremely young ($`<2.5`$ Myr; Hillenbrand 1997; Palla & Stahler 1999), contains 4000–10000 stars within a region with a diameter of at least 5 pc and has a very high central number density of about $`10^{4.7}`$ stars/pc3 (McCaughrean & Stauffer 1994). It also has a velocity dispersion that is too large for it’s known mass (Jones & Walker 1988). Thus, it ought to be expanding, which may be a result of very recent ($`\begin{array}{c}<\hfill \\ \hfill \end{array}0.5`$ Myr ago) loss of a substantial amount of gas when the central OB stars “ignited”. It has been debated for a long time whether the ONC actually is a proto-Galactic cluster, or if it is in the process of forming an unbound OB association (e.g. Zinnecker, McCaughrean & Wilking 1993). Indeed, a large fraction of the ONC stars are optically visible and thus not embedded in gas (Hillenbrand & Hartmann 1998), and its central region, the Trapezium Cluster, contains less than $`10M_{}`$ in gas (Wilson et al. 1997). The large velocity dispersion implies that an unbound association is likely to result, but a recent investigation of the set of possible dynamical solutions strengthens the possibility that the ONC is close to forming a bound cluster (Kroupa 2000a). There is no discernible sub-structure in the Trapezium so that it is at least a few crossing times ($`t_{\mathrm{cross}}`$) old (Bate, Clarke & McCaughrean 1998).

Under which conditions some embedded clusters survive to become bound entities, such as the Pleiades Cluster, remains something of a mystery (see Elmegreen et al. 2000; Clarke, Bonnell & Hillenbrand 2000 for reviews). It is known though that in the process of forming an embedded cluster a relatively large mass of gas, typically $`\begin{array}{c}>\hfill \\ \hfill \end{array}500M_{}`$, is squeezed into a volume of typically one pc3. The gas dynamics is complex, but once stars form, they remove a substantial fraction of this gas through outflows and winds. Especially in those circumstances when O stars form and start shining, gas removal can be expected to be very rapid and faster than the dynamical crossing time of the embedded cluster (Whitworth 1979). Observations indicate that embedded massive stars have dynamically young ($`\begin{array}{c}<\hfill \\ \hfill \end{array}10^4`$ yr) outflows that contain more mass than the powering star (Churchwell 1997, 1999). In addition to the expansion of the HII region by virtue of it’s overpressure, and the mass loss from the central region containing the massive star through the outflow, a massive star produces a powerful wind that is very effective in expelling mass from the cluster (Garay & Lizano 1999).

Star-formation is typically less than 40 per cent efficient (Lada 1999; Elmegreen et al. 2000; Clarke et al. 2000), and once more than 50 % of the mass of a cluster is removed instantaneously, an unbound association is deemed to result, unless the cluster was beginning to collapse, or was virialised after collapse, at the instant of gas removal. Such results that laid the conceptual framework for this work and which were obtained at a time when infrared observations did not exist to constrain the star-formation efficiency (sfe) or embedded cluster morphology and membership, were arrived at using either a direct $`N`$-body code being limited to $`100`$ stars (Lada, Margulis & Dearborn 1984), or necessarily simplifying analytical arguments (Hills 1980; Elmegreen 1983; Mathieu 1983) with more recent generalisations (Pinto 1987; Verschueren & David 1989; Verschueren 1990). A particularly interesting result that emerged from this early ground-breaking work was the realization that, unless the sfe is larger than 50 per cent, Galactic (i.e. ’open’) clusters can not form containing OB stars, because these drive out the gas faster than the dynamical crossing time leading to complete dispersion of the stellar system. With a state-of-the-art $`N`$-body code and advanced observational constraints, it is now possible to re-visit this problem with the aim of investigating the likely fate of the ONC as a first application.

The ONC can be in one of three dynamical states (Kroupa 2000a): (i) in virial equilibrium, (ii) collapsing, or just after collapse and in or past the associated violent relaxation phase, or (iii) expanding. Case (i) certainly leads to a bound cluster, but is very unlikely because it implies an sfe of essentially 100 % ($`ϵ=1`$), which has never been observed. Case (ii) may occur if stars form out of the molecular cloud with a smaller-than-virial velocity-dispersion, so that the stellar system contracts. Rapid gas removal will then leave a bound cluster, if the gas is removed rapidly at onset of collapse, or after virialisation. If even a relatively small fraction of the gas is removed during the main collapse phase, however, then an unbound association is likely to result (Hills 1980). In this case the observer should be seeing the ONC in the special time just before or after maximum central density with a radial infall or expanding bulk motion. However, such a scenario requires star formation to be finished within less than one proto-cluster crossing time ($`\begin{array}{c}<\hfill \\ \hfill \end{array}0.4`$ Myr for the ONC), and thus a physically implausible degree of synchronisation across the whole cluster-forming cloud.

In reality, stars appear to continue forming in a cluster for a few Myr until the gas is expelled (Palla & Stahler 2000), and once a star forms it decouples from the gas and either falls towards the potential minimum or is already on a “virial equilibrium orbit”. In any case, the accumulating stellar system will be virialising constantly as new stars are added, while the stars formed less than a crossing time ago may be falling towards the cluster centre. Gas-drag may cause the stellar system to contract enhancing the effective sfe (Saiyadpour, Deiss & Kegel 1997). Thus, most of the cluster is likely to be in virial equilibrium when the first massive stars “ignite” with devastating effects for the residual (but still mass-dominating) gas. This scenario is supported by various observations which show star-formation on-going in a clustered environment (e.g. Megeath et al. 1996).

In this paper the one scenario (iii) that is least likely to form a bound cluster and which is the most intuitive, is studied in the realistic situation of very rapid gas removal using a newly completed extension of the code Nbody6 (Aarseth 1999), GasEx. It is a modified version of Nbody6, by incorporating the gas in an embedded cluster as a time-evolving background potential. The essential physics of gas expulsion is thus described in a computationally-efficient way, while retaining the necessary highly accurate treatment of point-mass stellar encounters. GasEx also allows the setup of a realistic primordial binary population according to Kroupa (1995b), and a new five-part power-law IMF-generating routine extends the mass range into the brown dwarf (BD) regime.

Calculations are performed of clusters that initially resemble the ONC and which contain twice as much mass in gas than in stars. After 0.6 Myr the gas is removed on a time-scale given by the velocity of the heated ionised gas (10 km/s, e.g. Hills 1980). Despite the essentially instantaneous gas loss, roughly $`1/3`$ of the cluster remains bound and forms a cluster remarkably similar to the Pleiades. This result differs from previous work, as is discussed in Section 5.

The code is described in Section 2, and the adopted initial conditions are detailed in Section 3. Section 4 contains the results and predictions. A discussion is provided in Section 5, and the paper is concluded with Section 6.

## 2 THE CODE

A few points concerning Nbody6 are stressed, and the inclusion of the gas potential is described in this section.

### 2.1 The code Nbody6

The stellar-dynamical interactions are treated with Nbody6, which is described elsewhere in detail (Aarseth 1999, 2000).

A direct $`N`$-body code must deal efficiently with a range of dynamical time-scales spanning many orders of magnitude. To achieve this, special mathematical techniques are employed to transform the space-time coordinates of closely interacting stars such that the resulting equations of motion of the sub-system are regular (Mikkola & Aarseth 1993), that is, the calculation can proceed accurately through the near time-step singularity that occurs during an encounter. This regularisation extends from the two-body case to a chain containing up to six members. Perturbations acting on the chain from close non-chain-members are incorporated, and allowance is made for approaching stars or binaries to become new members, and old members to be ejected from the usually unstable small-body group.

In addition to the efficient treatment of close encounters and transient compact sub-systems, the code assigns a neighbour radius for each star which typically contains the nearest $`\sqrt{N}`$ stars, where $`N`$ is the total number of stars in the cluster. The irregular forces from these neighbours are computed on the smallest time-step of their central star, and the regular forces from the stars lying outside the neighbour sphere are added at less frequent intervals. The Ahmad-Cohen (Ahmad & Cohen 1973) neighbour-scheme makes the CPU time scale approximately with $`N^{1.6}`$ rather than the usual $`N^2`$ in a “primitive” direct $`N`$-body code. Each star has it’s own regular and irregular time-step, which is adjusted throughout the computation, but is alway commensurate with 2. That is, block-time-steps are used in this version, as opposed to the basic individual time-step scheme in Nbody6’s precursor, Nbody5, which was previously applied to the problem of young clusters (e.g. Kroupa, Petr & McCaughrean 1999, hereinafter KPM). Block-time-steps are used in connection with Hermite integration, and CPU-timing tests show that it is more accurate and at least as fast as Nbody5 for the same number of integration steps.

State-of-the-art stellar evolution is incorporated using analytical fitting functions (Hurley, Pols & Tout 2000). These functions allow the generation of synthetic HR diagrams for a wide range of the metallicity. However, in these simulations we adopt the standard solar value (Z = 0.02). The resulting mass loss is implemented in a nearly continuous manner which facilitates an energy-conserving scheme (Aarseth 1999). Although the code contains a simple scheme for tidal circularisation, this process has not been included up to now, with consequent neglect of Roche-lobe mass transfer. The code assigns kick velocities when the massive stars explode, so that high-velocity neutron stars are readily produced. This, however, is not a topic of the present study.

### 2.2 The gas

To represent the essential physics of the dynamical evolution of an embedded cluster when the gas is removed, a time-varying analytical spherical background potential is added, following the pioneering study of Lada et al. (1984). This avoids the daunting task of a cluster-scale three-dimensional magneto-hydrodynamical computation with stellar feedback, but includes the important dynamical effect of reducing the containing potential as the gas is blown out. That the analytical approximation for the gas leads to physically realistic behaviour is concluded by Geyer & Burkert (2000), who apply collisionless $`N`$-body computations to this problem assuming an analytical gas as well as a more realistic treatment with smoothed-particle hydrodynamics (SPH).

Both the embedded star cluster and the parent cloud are assumed to be Plummer spheres. A star with a position vector $`𝐑`$ relative to the centre of the cluster experiences an acceleration from the background potential,

$$𝐚=\frac{GM_\mathrm{g}(t)}{(R^2+R_{\mathrm{pl},\mathrm{g}}(t)^2)^{3/2}}\frac{𝐑}{R},$$
(1)
where $`M_\mathrm{g}(t)`$ and $`R_{\mathrm{pl},\mathrm{g}}(t)`$ are the time-varying mass and Plummer radius of the gas cloud, and G is the gravitational constant. The background potential is assumed to begin evolving after a delay time $`t_\mathrm{D}`$ according to

$$M_\mathrm{g}(t)=M_\mathrm{g}(0)e^{\left(\frac{(tt_\mathrm{D})}{\tau _\mathrm{M}}\right)},tt_\mathrm{D},$$
(2)
and similarly the radius evolves as

$$R_{\mathrm{pl},\mathrm{g}}(t)=R_{\mathrm{pl},\mathrm{g}}(0)\left[1+\left(\frac{(tt_\mathrm{D})}{\tau _\mathrm{R}}\right)^{1/2}\right],tt_\mathrm{D},$$
(3)
where $`\tau _\mathrm{M}`$ and $`\tau _\mathrm{R}`$ are the time-scales for mass and radius evolution, respectively. In the present application $`R_{\mathrm{pl},\mathrm{g}}(t)`$ remains unchanged.

In general, the changes to $`M_\mathrm{g}`$ and $`R_{\mathrm{pl},\mathrm{g}}`$ are applied at constant time intervals, $`\delta t_\mathrm{a}<<\mathrm{min}(\tau _\mathrm{M},\tau _\mathrm{R})`$. To conserve energy at every adjustment the resulting change in potential energy of the system of $`N`$ stars is obtained by summing over the changes of each star with mass $`m_\mathrm{i}`$, giving $`\delta \mathrm{\Phi }=_{i=1}^Nm_\mathrm{i}\delta \mathrm{\Phi }_\mathrm{i}`$, where

$$\delta \mathrm{\Phi }_\mathrm{i}=\frac{GM_\mathrm{g}(t+\delta t_\mathrm{a})}{\left(R_\mathrm{i}^2+R_{\mathrm{pl},\mathrm{g}}^2(t+\delta t_\mathrm{a})\right)^{1/2}}+\frac{GM_\mathrm{g}(t)}{\left(R_\mathrm{i}^2+R_{\mathrm{pl},\mathrm{g}}^2(t)\right)^{1/2}}.$$
(4)
The contributions by the gas potential are neglected once $`M_\mathrm{g}(t)<10^5M_{}`$.

## 3 INITIAL CONDITIONS

The initial properties of the stellar population and of the embedded cluster models are described.

### 3.1 The initial stellar and binary population

Initial stellar masses are distributed according to the universal three-part power-law IMF (Kroupa 2000c) with lower and upper mass limits $`m_\mathrm{l}=0.01M_{}`$ and $`m_\mathrm{u}=50M_{}`$, respectively:

$$\xi (m)m^{\alpha _i},$$
(5)
where

$$\begin{array}{cc}\alpha _0=+0.3,\hfill & 0.01m/M_{}<0.08,\hfill \\ \alpha _1=+1.3,\hfill & 0.08m/M_{}<0.50,\hfill \\ \alpha _2=+2.3,\hfill & 0.50m/M_{}.\hfill \end{array}$$
(6)
Here, $`\xi (m)dm`$ is the number of stars in the mass interval $`m`$ to $`m+dm`$. An often useful description is achieved via the logarithmic IMF,

$$\xi _\mathrm{L}(m)=\xi (m)m\mathrm{ln10},$$
(7)
where $`\xi _\mathrm{L}`$ is the number of stars in the logarithmic mass interval log$`{}_{10}{}^{}m`$ to log$`{}_{10}{}^{}m+dm`$. This is the form used in later figures of this paper.

The above IMF leads to the following stellar population: 37 % BDs ($`0.010.08M_{}`$) contributing 4.3 % to the stellar mass, 48 % M dwarfs ($`0.080.5M_{}`$) contributing 28 % mass, 8.9 % “K” dwarfs ($`0.51.0M_{}`$) contributing 17 % mass, 5.7 % “intermediate mass (IM) stars” ($`1.08.0M_{}`$) contributing 34 % mass, and 0.37 % “O” stars ($`>8M_{}`$) contributing 17 % mass.

The total binary proportion is

$$f_{\mathrm{tot}}=\frac{N_{\mathrm{bin}}}{(N_{\mathrm{bin}}+N_{\mathrm{sing}})},$$
(8)
where $`N_{\mathrm{bin}}`$ and $`N_{\mathrm{sing}}`$ are the number of binary and single-star systems, respectively.

Initially, all stars are assumed to be in binary systems ($`f_{\mathrm{tot}}=1`$) with component masses $`m_1,m_2`$ chosen randomly from the IMF. This gives approximately a flat distribution of mass-ratios, $`q=m_2/m_11`$, with a maximum at $`q=1`$ as a result of the model adopted by Kroupa (1995b) for system-internal processes, such as star-disk interactions, tidal circularisation and mass ejection that redistribute mass, angular momentum and energy within the binary system while it is still very young ($`<10^5`$ yr). For convenience, this evolution is collectively referred to as ’pre-main sequence eigenevolution’, to differentiate it from the externally-induced perturbations from nearby stars. Random pairing over the entire mass range implies that most massive primaries will initially have M dwarf and BD companions, which may not be consistent with observational constraints. Their numbers are so small though, that they do not significantly alter the overall mass-ratio distribution, $`f_q`$. The fast sinking to the cluster centre leads to rapid (within a few $`t_{\mathrm{cross}}`$) changes in $`f_q`$ for massive primaries as they pair up with massive secondaries near the cluster core, but details await a further study.

The birth period distribution, with periods ranging from about 1 d to $`10^9`$ d, is constructed following Kroupa (1995b, eq.8), and pre-main sequence eigenevolution gives the observed correlation between orbital period and eccentricity for short-period systems with stellar masses in the range $`0.081.1M_{}`$. Only systems with a period $`P\begin{array}{c}<\hfill \\ \hfill \end{array}10^4`$ d are affected, the overall distribution of masses not being changed greatly. As a result, some binaries merge, so that at $`t=0`$ the binary proportion is slightly below unity. The initial binary proportion is reduced further through crowding, that is, immediate disruption because some long-period systems have companion separations larger than the system–system distances in the initial cluster, and because of the cluster tidal field, giving the true initial binary proportion.

### 3.2 Cluster models

The star clusters are assumed to have spherical Plummer number-density profiles initially (e.g. Aarseth, Hénon & Wielen 1974) with half-mass radii $`R_{0.5}`$ and mass $`M_{\mathrm{st}}`$. The true initial conditions of a cluster are not known, and are likely to be very complex with significant sub-structure, which evolves on a crossing time-scale with star formation continuing for some time until the gas is expelled (e.g. Klessen, Burkert & Bate 1998). A more detailed model will not be attempted until the evolution of the simpler spherical case with gas expulsion as set up here is understood.

The number of stars, $`N_{\mathrm{in}}`$, and initial central density, $`\rho _\mathrm{C}`$, are taken from Kroupa (2000a), where possible initial configurations of the ONC are constrained. $`N_{\mathrm{in}}=10^4`$ is an upper limit to the initial number of stars in the ONC, and in expanding models $`\rho _\mathrm{C}`$ must be at least as large as the current ONC value.

Stellar masses and binary properties are assumed to be uncorrelated with their location in the cluster, giving an average stellar and system mass that is constant with $`R`$. Again, this assumption is not likely to be realistic because observational evidence and theoretical arguments suggest that massive stars may form near cluster centres (Bonnell, Bate & Zinnecker 1998). The present assumption is made because it allows measurement of the rapidity with which dynamical mass segregation occurs when calculated with a high-precision $`N`$-body code (touched upon in Fig. 13). This, however, will be the focus of another contribution.

The gas is assumed to have the same spatial distribution as the stars with a half-mass radius $`R_{0.5,\mathrm{g}}=R_{0.5}`$. It is “expelled” by reducing its mass, $`M_\mathrm{g}(t)`$, at time $`tt_\mathrm{D}`$ on a time-scale $`\tau _\mathrm{M}`$. $`M_\mathrm{g}(0)=2M_{\mathrm{st}}`$, i.e. an sfe $`ϵ0.3`$, is assumed.

The Galactic tidal field is instrumental in tidally truncating an expanding cluster. A standard Galactic tidal field is adopted using the linearised approximation (Terlevich 1987) with Oort’s constants $`A=14.4`$ km/s/kpc and $`B=12.0`$ km/s/kpc, and a mass density in the solar neighbourhood of $`0.11M_{}`$/pc3. Stars are removed from the $`N`$-body calculation once they reach a distance of 100 pc from their cluster, but are kept in memory without advancing them on their Galactic orbit, but with stellar evolution included, to facilitate data reduction.

Two models are discussed here. The model parameters are listed in Table 1.

## 4 RESULTS

The evolution of the two models is discussed, with a comparison to some observational constraints available for both the ONC and the Pleiades. These constraints include the structural and kinematic properties, as well as the binary proportion and period distributions. The section ends with a discussion of the stellar and system mass functions.

Throughout this paper, $`R`$ refers to the 3D radial distance from the density centre of the cluster, whereas $`r`$ is the projected 2D distance. The projection onto the “observational plane” is arbitrary, except that the vertical direction is always taken to be perpendicular to the Galactic disc. Pictures of the distribution of stars in the model clusters at $`t=100`$ Myr thus show flattened clusters with shortest extension along the y-direction (c.f. Terlevich 1987; Portegies Zwart et al. 2000).

### 4.1 Cluster properties vs observations

Not much is known about the stellar-dynamical evolution of a cluster for which the velocity dispersion is increased as a result of an additional containing potential. In the cases considered here, the 3D velocity dispersion $`\sigma _{3\mathrm{D}}=\sqrt{3}\sigma _{3\mathrm{D},\mathrm{st}}`$, where $`\sigma _{3\mathrm{D},\mathrm{st}}`$ is the velocity dispersion if no gas were present. The crossing time is thus shortened by the factor $`\sqrt{3}`$, and stellar interactions occur $`\sqrt{3}`$ times as often, because the stellar density remains unchanged. This topic requires more study, and one result of the present calculations is that for example $`f_{\mathrm{tot}}`$ (eq. 8) decays in the same manner as for a cluster without gas but the same density distribution, which is somewhat surprising. Hence the increased number of encounters is compensated by the reduced encounter duration.

Evaluating the mean radial distance of stars in various mass bins as a function of time shows that mass segregation is very weak in model A at the time of gas expulsion, whereas it develops to be quite pronounced in model B, leading to a slight expansion of stars with mass $`m<8M_{}`$, as a result of the associated heating of the cluster. This will be re-visited briefly in Section 4.4 below.

At time $`tt_\mathrm{D}=0.6`$ Myr the gas mass evolves according to eq. 2 with $`\tau _\mathrm{M}`$ given in Table 1. In both models, $`M_\mathrm{g}(t)0`$ for $`t\begin{array}{c}>\hfill \\ \hfill \end{array}1`$ Myr, the background potential being insignificant thereafter, and most of the original cluster expands as an OB association. Fig. 1 depicts this expansion, but also shows that between 20 and 30 % of the stellar population remains within the approximate tidal radius (Binney & Tremaine 1987),

$$R_{\mathrm{tid}}(t)=\left(\frac{M_{\mathrm{st}}(t)}{3M_{\mathrm{gal}}}\right)^{\frac{1}{3}}R_{\mathrm{GC}},$$
(9)
with $`M_{\mathrm{gal}}=5\times 10^{10}M_{}`$ being approximately the Galactic mass enclosed within the distance of the Sun to the Galactic centre, $`R_{\mathrm{GC}}=8.5`$ kpc. To estimate $`R_{\mathrm{tid}}(t)`$, $`M_{\mathrm{st}}(t)`$ is calculated by summing only those stars which have $`R(t)2R_{\mathrm{tid}}(t\delta t_{\mathrm{op}})`$, where the data output time interval $`\delta t_{\mathrm{op}}<<t_{\mathrm{relax}}(t)`$. $`R_{\mathrm{tid}}(t)`$ is shown in Fig. 1 for $`t>0.6`$ Myr merely as an aid to the eye to facilitate a rough estimate of how many stars form the bound cluster by 100 Myr ($`R_{\mathrm{tid}}`$ is not used to remove stars from the $`N`$-body calculation). The bound cluster that forms from the expanding OB association contracts with time as $`R_{\mathrm{tid}}`$ decreases, and the density profiles at 100 Myr are discussed further below (Fig. 6).

The core radius, $`R_\mathrm{C}`$, is a useful quantity because it can be constrained observationally much more easily than the half-mass radius or tidal radius. In the $`N`$-body computations it is approximated by calculating the density-weighted radius $`R_{\mathrm{HA}}`$ (Heggie & Aarseth 1992),

$$R_{\mathrm{HA}}^2=\frac{\underset{i=1}{\overset{N_{20}}{}}R_\mathrm{i}^2\rho _\mathrm{i}^2}{_{\mathrm{i}=1}^{N_{20}}\rho _\mathrm{i}^2},$$
(10)
where $`\rho _\mathrm{i}=3m_{\mathrm{i},5}/(4\pi d_{\mathrm{i},5}^3)`$ is the density around star i estimated within the closest distance $`d_{\mathrm{i},5}`$ from star i containing $`n=5`$ additional stars with combined mass $`m_{\mathrm{i},5}`$, and $`R_\mathrm{i}`$ is the distance of star i from the density centre of the cluster. The summation extends only over the innermost 20 % of all stars in the cluster, $`N_{20}`$, since this is sufficient to ensure convergence. For a comparison with the core radius, $`R_\mathrm{C}`$, obtained from fitting King density profiles, as is often done for observed clusters such as the ONC and the Pleiades,

$$R_\mathrm{C}R_{\mathrm{HA}}/0.8,$$
(11)
(Heggie & Aarseth 1992; Giersz & Spurzem 2000). The evolution of $`R_\mathrm{C}=R_{\mathrm{HA}}/0.8`$ is shown in Fig. 1. It follows the 5 and 10% Lagrange radii rather well.

The formation of a bound cluster despite $`ϵ0.3`$, despite explosive mass loss and additional mass-loss from evolving massive stars, and despite the Galactic tidal field, is evident by a contraction of $`R_\mathrm{C}`$ for $`t\begin{array}{c}>\hfill \\ \hfill \end{array}5`$ Myr (model A) and $`t\begin{array}{c}>\hfill \\ \hfill \end{array}1`$ Myr (model B). Evolving massive stars cause $`R_\mathrm{C}`$ to expand in model B for $`t\begin{array}{c}>\hfill \\ \hfill \end{array}2.0`$ Myr, and in both models $`R_\mathrm{C}1.5`$ pc when $`t=150`$ Myr is reached. The core radius fluctuates by about 0.5 pc as a result of Poisson noise. It is rather remarkable that $`R_\mathrm{C}`$ ends up being so similar for both models despite the significantly different initial density.

Thus, in the presence of a tidal field, models of the same cluster mass converge to a similar structural state independent of the initial central density used. That is, after some time the inner regions begin to feel the effect of the contracting tidal radius and “self-similar” evolution sets in: the core- and half-mass radii then depend only on the total cluster mass even though, for each individual model, the evolution path of the core radius to this point may be quite different. This was evident in the cluster models with low-mass stars computed by Kroupa (1995c), and in more realistic models of M67 by Hurley (2000) and Hurley et al. (2000b).

The core radius is compared with observational constraints for the ONC and the Pleiades in Fig. 2. Both models are in reasonable agreement with the constraints. Model B, however, fits better near $`t=1`$ Myr when Poisson fluctuations are still small. The same is true for the number of systems shown in Fig. 3. Comparison of the available velocity dispersion measurements with the models also leads to agreement (Fig. 4).

The projected radial density profiles are further constraints. Ideally, the observed ONC and Pleiades profiles should be reproduced by one model at $`t1`$ Myr and $`100`$ Myr, if the Pleiades is a later rendition of the ONC. Model A fails with the ONC (Fig. 5) in that the profile is too flat for $`r\begin{array}{c}<\hfill \\ \hfill \end{array}1`$ pc, and too steep at larger radii, when $`t0.9`$ Myr. The model, however, fits the Pleiades rather well (Fig. 6). Conversely, model B fits the ONC reasonably well for $`t1.1`$ Myr (Fig. 5), but contains too many stars compared to the Pleiades at $`t100`$ Myr (Fig. 6).

The deviations and agreements found here are very useful in discerning what density profile and $`N`$ the ONC and Pleiades precursors are likely to have had for a better reproduction of either cluster, and perhaps both simultaneously. However, it would have to be considered an impressive “quirk” of nature if it turns out that the ONC actually can be viewed as a precursor of the Pleiades. The present results have demonstrated that this is indeed the case, to a surprising degree of accuracy. Further work will be needed to study how much the various parameters (radial profile, $`N`$, IMF, binaries, $`ϵ`$) must be varied in order to get one model with improved agreement for both clusters.

### 4.2 Prediction for the bulk motion and velocity dispersion

If the ONC is the precursor of a Pleiades-like open cluster, then it is expanding now. The possible expansion of the Trapezium, which is essentially the core of the ONC (Hillenbrand & Hartmann 1998), was already considered by KPM. In their model, which assumed $`ϵ=0.42`$ and a higher central number density but with only 1600 stars, the projected bulk radial velocity $`<v_\mathrm{r}>=2.9`$ km/s for stars with $`r<0.41`$ pc about 60000 yr after instantaneous gas expulsion.

Models A and B can be evaluated to give the predicted expansion rate of the ONC. The projected bulk radial motion, $`<v_\mathrm{r}>`$, which is the average of the radial component of the projected 2D velocity vectors of all stars with $`rr_\mathrm{o}`$ in the observational plane, is plotted in Fig. 7. Model A gives an expansion velocity of $`<v_\mathrm{r}>=2.4`$ km/s ($`r_\mathrm{o}=2.5`$ pc) and $`<v_\mathrm{r}>=0.6`$ km/s ($`r_\mathrm{o}=0.41`$ pc) at an age $`t=0.9`$ Myr (log$`{}_{10}{}^{}t=0.05`$), whereas model B expands with $`<v_\mathrm{r}>=2.7`$ km/s ($`r_\mathrm{o}=2.5`$ pc) and $`<v_\mathrm{r}>=0.4`$ km/s ($`r_\mathrm{o}=0.41`$ pc) at the same age. The bulk radial velocity drops when the fastest stars leave the field of view.

The line-of-sight velocity dispersion, $`\sigma _{\mathrm{los}}`$, decays at a different rate after gas expulsion than the projected velocity dispersion, because fast moving stars remain in the observational field for a longer time along the line of sight. Thus a difference between the two is expected. This is shown in Fig. 8, where the line-of-sight and the projected 1D velocity dispersions are compared. As expected, $`\sigma _{\mathrm{los}}>\sigma _{1\mathrm{D},\mathrm{c}}`$, where $`\sigma _{1\mathrm{D},\mathrm{c}}`$ is the projected velocity dispersion corrected for the bulk radial motion. This is the velocity dispersion obtained from ground-based proper motion studies, which apply special photographic-plate-reduction techniques to eliminate distortions which remove bulk rotation and bulk radial motions (see e.g. Jones & Walker 1988; this is also discussed in KPM). For model A at $`t=0.9`$ Myr, $`\sigma _{\mathrm{los}}=2.6`$ km/s and $`\sigma _{1\mathrm{D},\mathrm{c}}=2`$ km/s for $`r2.5`$ pc and $`r0.41`$ pc. For model B, $`\sigma _{\mathrm{los}}=2.8`$ km/s and $`\sigma _{1\mathrm{D},\mathrm{c}}=2`$ km/s for $`r2.5`$ pc with smaller values if the measurement is confined to $`r0.41`$ pc.

The planned astrometric satellite missions DIVA (Röser 1999) and GAIA (Lindegren & Perryman 1996; Gilmore et al. 1998) should verify these predictions of $`<v_\mathrm{r}>`$ for different $`r_\mathrm{o}`$.

A further generic prediction of these models is that the expanding population of unbound stars (about $`2/3`$ of $`N`$ here) forms, at the age of the Pleiades, a distinct co-eval group with an extension of roughly 1 kpc ($`\sigma _{3\mathrm{D}}\times 100\mathrm{Myr}\times 2/\sqrt{3}`$, $`\sigma _{3\mathrm{D}}9`$ pc/Myr from Table 1) and co-moving with the Pleiades. This moving group, which we refer to as ’group I’, differs from the classical moving ’group II’ that forms later from the long-term evaporation of Pleiades members through two-body relaxation. The difference is kinematical in nature, in that members of group II diffuse away from their parent cluster mostly through the two Lagrange points with a small relative velocity of typically $`\begin{array}{c}<\hfill \\ \hfill \end{array}1`$ km/s (Terlevich 1987; de La Fuente Marcos 1997; Portegies Zwart et al. 2000). The planned astrometric missions will be able to distinguish these two moving groups associated with the Pleiades, in which case this cluster-formation theory will be ascertained. Group I may in fact have already been observed as the ’Pleiades supercluster’ by Eggen (e.g. Eggen 1998).

### 4.3 The binary population

The distribution of periods of binary systems measured in the Pleiades (and the ONC) are not the birth distributions, since these have been changed through dynamical interactions (KPM). The dynamical models studied here allow this evolution to be quantified. The most obvious hypothesis that one would consider testing is whether the Pleiades and ONC binary populations are consistent with the same populations seen in Taurus–Auriga, i.e. that the birth binary properties do not depend on environment, apart from disruption through crowding and the cluster tidal field. Thus, models A and B assume Kroupa’s (1995b, eq. 8) birth distribution, which is essentially the Taurus–Auriga period distribution of pre-main sequence binaries, and in this section the evolved distributions are compared to the observational constraints available for the ONC and the Pleiades. The stellar-dynamical interaction between a primordial binary population and a young cluster is described in some detail in Kroupa (2000b).

The evolution of the binary proportion of primaries in five mass ranges for stars remaining in the cluster ($`R3.2`$ pc) is shown in Fig. 9. Initially, the binary proportion of the massive stars, $`f_\mathrm{O}`$, is reduced through disruption from crowding in the dense clusters and through interactions, but it increases as the massive stars gain new companions. A topic for future study will be the comparison of the period and mass-ratio distributions of the massive stars as a function of time with observational constraints. For example, Mason et al. (1998) find a very high proportion of binaries among massive stars, as is obtained in the present calculations by $`t>2.5`$ Myr in model A, and $`t>16`$ Myr in model B. For the massive binaries, $`f_\mathrm{O}`$, Preibisch et al. (1999) find that in the ONC the O and B primaries have on average 1.5 companions. A more detailed analysis of the model data will have to be performed to quantify the number of triple and quadruple systems among the massive stars, before a final conclusion concerning the primordial characteristics of massive binaries can be made.

At the other extreme, BD binaries are more easily disrupted at an early stage owing to their weaker binding energy, and by $`t1`$ Myr, $`f_{\mathrm{BD}}0.35`$ (model A) and 0.2 (model B). Overall, the Pleiades and ONC observational constraints for late-type binaries are in acceptable agreement with the models.

The distribution of periods for late-type stellar ($`0.081.5M_{}`$) and BD primaries is shown in Fig. 10 at $`t=0.9`$ Myr and Fig. 11 at $`t=100`$ Myr. The distributions in the former figure are compared with the initial (after pre-main sequence eigenevolution) distributions, and with observational constraints available for the ONC. By $`t=0.9`$ Myr substantial depletion of the period distribution has occurred in both models. However, only model B leads to agreement with the measured proportion of binaries with periods in the range $`P10^510^{6.5}`$ d, model A retaining too many binaries with these periods. Model B is also consistent with the absence of long-period ($`P>10^7`$ d) binaries in the ONC (Scally, Clarke & McCaughrean 1999).

The models are compared with the Pleiades constraints in Fig. 11, where a comparison is also made between the period distributions for stars inside and outside the cluster. Both models are consistent with the Pleiades data, although an intermediate model in terms of central density would lead to improved agreement.

The distributed population of stars formed after gas expulsion (moving group I, Section 4.2), has a higher binary proportion than the stars remaining in the cluster. This is evident in Fig. 11 through the larger number of orbits at $`P>10^5`$ d. In model A there are slightly too many field binaries relative to the Galactic field with periods in the range $`P=10^510^7`$ d, although the disagreement is marginal. Model B produces too few field binaries with $`P10^7`$ d. It is thus somewhat doubtful if either model can produce a Galactic field population, but it is to be noted here that unlike the empirical G-dwarf Galactic-field period distribution, the model period distributions do not measure triple and quadruple systems. The data reduction performed only finds the closest bound pairs, so that a third long-period companion is missed. This will be the subject of a future study.

In both figures, the BD binaries have a significantly depleted binary proportion over all periods, which is a result of their small binding energy.

A comparison of Figs. 10 and 11 shows that the birth distribution of orbital periods in the ONC could have been similar to that in Taurus–Auriga, which also leads to the Galactic field, if the ONC began with a central density between that assumed in models A and B. The same conclusion holds true for the Pleiades. Thus, the observational constraints do not exclude the hypothesis that the star formation conditions in the ONC and the Pleiades led to the same distribution of binary periods as in Taurus–Auriga. That is, the initial period distribution function may be universal.

### 4.4 The mass function

There is much interest in the IMF as it is one of the most important constraints on star-formation theories, as well as being the fundamental quantity entering many astrophysical problems. It is thus useful to study the changes in the stellar and system mass function (MF) as the clusters evolve, to obtain an insight into the sort of changes that make the quasi-observed MF (i.e. the system MF) differ from the IMF.

Pre-main sequence eigenevolution leads to small deviations of the IMF from eq. 6 shown by the solid dots in Fig. 12, because of the mass-gain of the secondary in close binary systems. This is evident in Fig. 12: the stellar IMF (thin dashed histograms) contains fewer BDs and more massive stars. The deviations are small and consequently we refer to both, the thin dashed histogram and eq. 6, as the IMF.

Assuming random pairing and $`f_{\mathrm{tot}}=1`$ gives the system MFs, which are plotted at $`t=0`$ (thin solid histograms) and $`t=0.9`$ Myr (thick solid histograms). At $`t=0`$, the system MF has a broad maximum near $`m=0.35M_{}`$. It falls below the stellar IMF by a factor of about 4 at the hydrogen burning mass limit and has $`\alpha _1=+0.6`$ (instead of $`+1.3`$) between $`0.08`$ and $`0.5M_{}`$. It falls below the IMF by an order of magnitude in the BD regime with $`\alpha _0=0.8`$ (instead of $`+0.3`$). This difference is reduced by $`t=0.9`$ Myr because most of the soft binaries have been disrupted, but significant discrepancies persist (this is studied in greater depth in Kroupa 2000c).

Mass segregation becomes evident through a surplus of massive stars in the central region. That this is not the case in model A by $`t=0.9`$ Myr is shown in the upper panel of Fig. 13. However, by virtue of the significantly shorter initial $`t_{\mathrm{cross}}`$, model B has evolved much more before gas expulsion, and shows significant mass segregation (lower panel in same figure). The binary proportion has also decreased much more (Fig. 9), so that the system and single star MFs are quite similar.

Similarly, Fig. 14 shows the IMF for stars and systems in comparison to the stellar and system MFs at $`t=100`$ Myr for stars with $`R15`$ pc, which is approximately $`R_{\mathrm{tid}}`$ of the model open cluster. Stellar evolution has removed all stars with $`m\begin{array}{c}>\hfill \\ \hfill \end{array}10M_{}`$. The system MF falls significantly below the stellar MF for $`m<0.25M_{}`$ in model A, by virtue of the larger surviving binary proportion. In model B (lower panel) the smaller binary proportion implies that the system MF (thick solid histogram) lies below the stellar MF (thick dashed histogram) only by an approximately constant factor of 2 for $`m<0.4M_{}`$. Table 2 summarises the census of systems with $`R15`$ pc.

The central ($`R2`$ pc) cluster region shows significant mass segregation at $`t=100`$ Myr in both models (Fig. 15), although it is much more pronounced in model B (lower panel). Matching the stellar MFs (dashed histograms) to the IMF (solid dots) near $`1M_{}`$, it can be seen that the central region has been relatively depleted in BDs. Table 3 summarises the census of systems with $`R2`$ pc.

Comparison to observational estimates of the MF for the Pleiades is not the main topic of this paper, but will be a focus in a future study. Kroupa (1995c) showed that the Pleiades LF is consistent with the Galactic field IMF, but those cluster models did not incorporate gas expulsion nor massive stars. At present, it suffices to state that the system MF presented in Fig. 14 has $`\alpha _0=0.0`$ for both models, i.e. the slope is similar to the actual IMF ($`+0.3`$). This is quite consistent with the estimate by Martin et al. (2000) ($`\alpha _00.53`$), and the number of BD candidate systems found (34), may be consistent with the present models (Tables 2 and 3) since their survey area extends to $`r8`$ pc for a Pleiades distance of 120 pc and they are not likely to have found the least massive BDs. The overall shape of the system MF is also similar to the observed MF (Meusinger, Schilbach & Souchay 1996, their fig.10; Hambly et al. 1999, their fig.11).

## 5 DISCUSSION

In the previous section it was seen that models A and B both reproduce the ONC and Pleiades to a remarkable degree, suggesting that the Pleiades may have looked similar to the ONC about 99 Myr ago. Additional calculations with GasEx but without a background gas potential, i.e. assuming $`ϵ=1`$, but with stellar evolution, demonstrate that, for the Pleiades to have its present mass and structure, $`R_\mathrm{C}1`$ pc and $`N3000`$ at birth. Such models are studied in detail by Portegies Zwart et al. (2000). An initial model in virial equilibrium without gas and as concentrated as the ONC can never evolve to have the low concentration (central number density $`\rho _\mathrm{C}11`$ stars/pc3, Pinfield, Jameson & Hodgkin 1998) which the Pleiades has now. It is natural to assume that the large core radius is a result of gas expulsion, and this applies to all known nearby open clusters.

However, before the explanation presented here of how open clusters form can finally be accepted as being correct, the bulk radial motion predicted by these models (Fig. 7) have to be verified through observations. The ONC must be expanding with $`<v_\mathrm{r}>23`$ km/s for stars with $`rr_\mathrm{o}=2.5`$ pc. The expansion velocity should be smaller for a smaller $`r_\mathrm{o}`$. This will be an important goal of the upcoming astrometric space observatories, DIVA (Röser 1999) and GAIA (Lindegren & Perryman 1996; Gilmore et al. 1998). If these satellites can measure such a bulk expansion of the ONC, and detect the Pleiades moving groups I and II (Section 4.2), then the claim that star-cluster formation has been solved can be made.

### 5.1 Cluster formation despite instantaneous expulsion of more than 50 % gas mass

Why does a substantial bound cluster form despite an sfe of $`ϵ0.3`$?

The underlying reason is very simple and has been realised for some years when the first expanding computations were done with Nbody5. In KPM the formation of a bound cluster was noted although the expectation was an expanding OB association.

In any stellar-dynamical system, the velocity distribution of the stars has a tail with radial velocities near 0. In the event of the expulsion of a large part of the mass, those stars which have turned around on their orbits remain bound. Moreover, the richness of the emerging cluster depends on the number of stars in the velocity distribution that have a velocity, $`v`$, smaller than the escape velocity, $`v_{\mathrm{esc}}`$, from the system after the mass (gas$`+`$unbound stars) is expelled. After mass expulsion, the population with initially $`v<v_{\mathrm{esc}}`$ has a non-equilibrium distribution of velocities, and as this distribution relaxes on the new dynamical time-scale, more stars will be lost. The final population of the cluster is therefore expected to depend sensitively on the violent relaxation process, and thus on the accurate treatment of close encounters as well as on the tidal field. For this reason GasEx was developed during the subsequent years 1998–1999 with additional improvements in the underlying Nbody6 code.

If the requirement of detailed predictions is relaxed, important insights into the above cluster-forming process can be gained through analytical work. Adams (2000) considers the stellar distribution function, which is a solution of the collisionless Boltzmann equation (CBE), in a potential that is a combination of stars and gas. Use of the CBE, and thus neglect of near-neighbour interactions, should be a reasonable approximation, as long as the crossing time of the cluster, which sets the time-scale for stellar interactions, is longer than the gas-expulsion time-scale. The distribution function gives the fraction of stars that remains bound after the gas is removed instantaneously. Adams finds that for $`ϵ=0.3`$ (the case assumed here) between $`=35`$ and 70 % of the stars in the initial cluster have velocities smaller than the escape velocity after gas removal. This work shows that the central concentration and a velocity anisotropy that is predominantly radial in the outer initial cluster regions lead to more massive remnant clusters.

However, in Adams’ models, the gas distribution is more extended than that of the stars, so that the effective $`ϵ`$ is larger within the stellar cluster, being as high as 0.9 (Geyer & Burkert 2000). Furthermore, the estimates obtained using this approach do not take into account the additional stars lost during the relaxation process to new dynamical equilibrium, nor the Galactic tidal field, and so constitute upper limits on $``$ without knowledge of the lower limits, which are the relevant ones for the problem.

It may thus not be surprising that Adams’ conclusions deviate substantially from those of Geyer & Burkert (2000), who employ a collisionless $`N`$-body code and adopt the same density distribution for the gas and stars, and show that no cluster forms if $`ϵ0.40`$ for instantaneous gas expulsion. An external Galactic tidal field (not applied by Geyer & Burkert) will further increase the critical $`ϵ`$ needed for cluster formation in the collisionless regime. The computations with GasEx, on the other hand, result in the formation of substantial clusters containing about $`1/3`$ of the initial number of stars, despite $`ϵ=0.3`$, despite applying a Galactic tidal field, and despite evolving stars, which together worsen the survival chance of an expanding system.

The reason for the formation of substantial clusters in our computations are gravitational interactions between neighbours during the radially expanding flow, which are not correctly treated in the collisionless approximation that underlies all hitherto available work. This mechanism leads to the emergence and amplification of a non-radial velocity dispersion even in the event of an initially purely radial outwards flow, and it operates most efficiently for the stars populating the slow tail of the velocity distribution. Kinetic energy is thus redistributed from the radial flow into orbital motions with non-zero angular momentum about the origin of the expanding flow, allowing a substantial part of the stellar system to condense as a bound entity, thereby defining and filling it’s tidal radius. Animated additional computations with GasEx sometimes show multiple clusters emerging from the radial flow, which, if massive enough so that two-body relaxation does not evaporate them too fast, survive and merge with the dominating sub-cluster to form the Galactic cluster. Essentially, the same mechanism operates in the inverse problem of a finite-$`N`$ system that collapses from rest. The increasing non-radial velocity dispersion that results from the grainy potential limits the maximum collapse factor in dependence of $`N`$, as has been studied in detail by Aarseth, Lin & Papaloizou (1988). Clearly, much work remains to be done along similar lines in the case of a radial outwards flow.

### 5.2 Comparison with similar work

That a bound “core” remains despite substantial mass-loss has also been noted by Lada et al. (1984) in their pioneering computations with $`100`$ stars. The models considered here correspond to the case $`\tau _\mathrm{R}<<t_{\mathrm{cross}}`$ ($`\tau _\mathrm{R}`$ being the gas-removal time) in their fig.2, although they varied the radius of the background gas, $`R_{\mathrm{pl},\mathrm{g}}`$, in contrast to the variation of the mass, $`M_\mathrm{g}`$, applied here. The Galactic tidal field is not incorporated in their investigation.

By scaling their results to a Pleiades-like cluster, they find that it’s precursor must have had a central density $`\rho _\mathrm{C}10^{5.2}`$ stars/pc3 assuming a mean stellar mass $`<m>=0.4M_{}`$ (their table 2 for $`\tau _\mathrm{R}=0`$ and $`ϵ=0.4`$). This is close to the actual central density observed in the Trapezium, and is intermediate to that assumed in models A and B here, which, however, assume $`ϵ=0.3`$ (note that Geyer & Burkert 2000, who revisit the Lada et al. problem with a collisionless code, find no cluster formation if $`ϵ=0.3`$).

The high central density of the Trapezium (McCaughrean & Stauffer 1994) was not known at the time of Lada et al. (1984), and they concluded that such high densities are not realistic, because this has not been observed in molecular clouds. They reason therefore that Pleiades-like clusters must form from less dense initial conditions with $`\tau _\mathrm{R}>>t_{\mathrm{cross}}`$. This can only be achieved in the absence of O stars, a suggestion also made by Elmegreen (1983) and Elmegreen & Efremov (1997). In view of this reasoning, a possible precursor to a Pleiades-like cluster appears to be the $`\rho `$ Oph cluster which contains no O stars. However, today we know that the stellar population in the $`\rho `$ Oph cluster is too sparse. It amounts to about 100 stars, and the total mass in gas and stars is about $`500M_{}`$ (Zinnecker et al. 1993; Luhman & Rieke 1999). In contrast, the Pleiades contains a stellar mass of between 500 and about 5000$`M_{}`$ (Pinfield et al. 1998; Raboud & Mermilliod 1998).

The line of argument followed here thus differs from the conclusions arrived at by Lada et al. (1984) and Elmegreen & Efremov (1997), because today we know that the very high central densities needed to form a Pleiades-like cluster from the ONC are available in the presence of O stars, and that O stars form in profusion in rich clusters. The calculations presented here demonstrate that even in the presence of O stars, which lead to explosive gas ejection, a substantial bound cluster remains.

The reaction of young clusters to mass loss is also studied by Goodwin (1997a; 1997b) and Geyer & Burkert (2000), but in the context of the formation of globular clusters. Goodwin uses the Nbody2 code with $`N=1000`$ particles, and applies a time-varying background potential in much the same way as done here, but also studies the effects of gas expulsion via an expanding shell. The gas is removed over a time-scale of a few Myr initially, and finally, after 10 Myr as a result of supernovae. Bound clusters form for $`ϵ0.25`$. Using a softened (“collisionless”) potential is, in principle (but see Section 5.1), a reasonable approximation for this purpose, because the evolution is restricted to being shorter than a relaxation time, but the initially slow gas depletion favours the formation of bound clusters, since the stellar orbits can adjust adiabatically to the varying potential, until the residual gas is blown out through supernovae. The results give useful insights into the likely initial conditions of massive clusters.

However, observational evidence appears to imply that gas expulsion occurs on a shorter time-scale. For example, the massive central R136 cluster in the 30 Doradus Nebula in the Large Magellanic Cloud has an age younger than about 2 Myr (Selman et al. 1999), but is already devoid of gas. In this respect, the collisionless computation by Geyer & Burkert (2000) are interesting, since they consider instantaneous expulsion and find that no cluster forms unless $`ϵ0.4`$. Their more realistic SPH treatment of the gas arrives at the same results as with an analytical background potential, lending credence to our approach.

### 5.3 An ONC peculiarity

One observational finding further supports the present result that the Pleiades may have formed from an ONC-like object, namely the short life-time of the circum-stellar features found around stars in the Trapezium. Henney & O’Dell (1999) measure the mass-loss rates of four of these objects, and find that they could not have been exposed to the UV flux from the central O stars for longer than about $`10^4`$ yr. These particular objects could be crossing the inner cluster region after spending more time at larger radii. However, the large proportion of stars with circum-stellar material does indicate that the destructive irradiation may indeed have turned on relatively recently. If this was the case, then it would appear natural to associate this event with the gas-expulsion time.

## 6 CONCLUDING REMARKS

This contribution presents the first results that were obtained using GasEx, thereby also being the first time that a high-precision $`N`$-body code that treats close stellar encounters accurately is applied to the problem of star-cluster formation. The results differ to those obtained with ’collisionless’ $`N`$-body codes and analytical estimates, in that substantial clusters form despite a low sfe ($`ϵ=0.3`$) and rapid gas-expulsion ($`\tau _\mathrm{M}<<t_{\mathrm{cross}}`$) due to the embedded OB stars, and despite applying a Galactic tidal field and mass-loss from evolving stars, which further limit the survival chances for an emerging star cluster.

Following Kroupa’s (2000a) identification of possible initial states of the ONC, two calculations were performed of a binary-rich embedded cluster with brown dwarfs. The extreme situation of very rapid gas expulsion together with a large number of massive stars is studied in order to investigate the likely fate of the ONC.

The rather startling result is that a Pleiades-like cluster forms containing about $`1/3`$ of the initial number of stars. This result is startling because all previous and contemporary work leads to complete cluster dissolution for such a low sfe, the reason being the neglect of the collisional nature of a finite-$`N`$ system.

A possible (non-causal) connection between the ONC and the Pleiades is thus established. This scenario suggests that Galactic clusters form with relatively large core radii and by filling their tidal radii, as a result of the expansion after gas loss from a compact state similar to the ONC. The scenario also suggests that they form as nuclei of expanding OB associations. The calculations show that the primordial binary population in the ONC and the Pleiades could have been indistinguishable from the Taurus–Auriga population, indicating that the initial period distribution function may be universal. An analytical description of this function is given by eq. 8 in Kroupa (1995b).

The most important prediction of this scenario of cluster formation is that the ONC must now be expanding, for it to be a precursor of a Pleiades-like cluster. The bulk expansion velocity and the line-of-sight and projected (proper motion) velocity dispersions are predicted in Figs. 7 and 8. This expanding population forms a moving group (referred to here as ’group I’), which differs from the classical moving ’group II’ that results from secular loss of cluster stars due to two-body relaxation. Future astrometric satellites will be able to confirm these predictions and differentiate between groups I and II. Either way, a very important lesson about cluster formation will be learned. An additional diagnostic of the dynamical state of a young cluster is the radial profile of the binary proportion, as stressed by KPM.

In this model of cluster formation, roughly $`2/3`$ of the initial stars become field stars immediately after gas expulsion. These stars retain a somewhat higher binary proportion than is seen in the Pleiades, and form an expanding OB association and later moving group I. The presence of wide binaries in the Galactic field poses important constraints on the possible overall contribution to Galactic field stars from such gas-expelled cluster stars (c.f. Scally et al. 1999). The same applies to massive field stars and their multiple properties. Clearly this poses a very rich field for further study, and will improve our knowledge of the Galactic field population.

The initial conditions assumed here (Table 1) are only two from the parameter space of possible solutions for an expanding ONC constrained by Kroupa (2000a). There is some redundancy in the parameters. Notably, a smaller $`N`$ can be compensated by a larger $`ϵ`$ and/or larger $`\alpha _3`$ (fewer massive stars) for the Pleiades, and it will be important to delineate the allowed range of these parameters for the Pleiades, Praesepe and Hyades clusters. Collapsing models will also be studied to obtain predictions.

The findings presented here have important implications for the formation of massive star clusters, some of which evolve to globular clusters. Observations of young massive clusters, for example R136 in the Large Magellanic Cloud, indicate that the gas is expelled extremely rapidly and possibly before the first massive stars explode. The observed haloes of stars (Elson, Fall & Freeman 1987) around some of the clusters suggest that a large proportion of the stars formed in massive embedded clusters become unbound, much in the way investigated here, and as also stressed by Goodwin (1997b). A significant proportion of Galactic spheroid (or ’stellar halo’) stars, that outnumber the stars in all globular clusters by at least 10:1 (Binney & Merrifield 1998), may thus be stars that formed in globular-cluster precursors but became unbound as a consequence of gas-expulsion. The distribution of the dynamical properties (binarity, and mass-ratio and period distributions) of the stars reflect such events, and the stars that got unbound during gas expulsion should have different phase-space properties than the stars lost from globular clusters as a result of their long-term evolution in the Galactic potential, for example by virtue of the general shrinkage of the globular cluster system as a result of dynamical friction.

Finally, a more realistic variation of the present models would include a primordial segregation of masses such that the most massive stars are located in the core of the embedded cluster. In the event of gas expulsion, the less massive stars located at larger radii would leave the cluster preferentially. This scenario implies that the MF in young clusters such as the Pleiades might be depleted in low-mass stars and BDs relative to their IMF, so that the study of the IMF via such young clusters would be compromised. Future work with GasEx will address this issue.

Acknowledgements PK acknowledges support through DFG grant KR1635.

### REFERENCES

Aarseth S.J., 1999, PASP, 111, 1333

Aarseth S.J., 2000, Gravitational N-Body Simulations, in prep.

Aarseth S.J., Hénon M., Wielen R., 1974, A&A, 37, 183

Aarseth S.J., Lin D.N.C., Papaloizou J.C.B., 1988, ApJ, 324, 288

Adams F.C., 2000, ApJ, 542, 964

Ahmad A., Cohen L., 1973, J.Comput.Ph., 12, 389

Bate M.R., Clarke C.J., McCaughrean M.J., 1998, MNRAS, 297, 1163

Binney J., Tremaine S., 1987, Galactic Dynamics, Princeton University Press, Princeton

Binney J., Merrifield M., 1998, Galactic Astronomy, Princeton University Press, Princeton

Bonnell I.A., Bate M.R., Zinnecker H., 1998, MNRAS, 298, 93

Bouvier J., Rigaut F., Nadeau D., 1997, A&A, 323, 139

Churchwell E.D., 1997, ApJ, 479, L59

Churchwell E.D., 1999, in The Origin of Stars and Planetary Systems, ed. C.J. Lada, N.D. Kylafis, NATO Science Series C, Vol 540, (Dordrecht: Kluwer), p.515

Clarke C.J., Bonnell I.A., Hillenbrand L.A., 2000, in Protostars and Planets IV , eds. V.G. Mannings, A.P. Boss, S.S. Russell (Tucson: Univ. Arizona Press), p.151

Duquennoy A., Mayor M., 1991, A&A, 248, 485

Eggen O.J., 1998, AJ, 116, 1810

Elmegreen B.G., 1983, MNRAS, 203, 1011

Elmegreen B.G., Efremov Yu.N., 1997, ApJ, 480, 235

Elmegreen B.G., Efremov Yu.N., Pudritz R.E., Zinnecker H., 2000, in Protostars and Planets IV , eds. V.G. Mannings, A.P. Boss, S.S. Russell (Tucson: Univ. Arizona Press), p.179

Elson R.A.W., Fall S.M., Freeman K.C., 1987, ApJ, 323, 54

de La Fuente Marcos R., 1997, A&A, 322, 764

Garay G., Lizano S., 1999, PASP, 111, 1049

Geyer M.P., Burkert A., 2000, MNRAS, submitted (astro-ph/0007413)

Giersz M., Spurzem R., 2000, MNRAS, 317, 581

Gilmore G., Perryman M., Lindegren L., et al., 1998, in Astronomical Interferometry, ed. R.D. Reasenberg, Proc. SPIE Vol. 3350, p. 541

Goodwin S.P., 1997a, MNRAS, 284, 785

Goodwin S.P., 1997b, MNRAS, 286, 669

Hambly N.C., Hodgkin S.T., Cossburn M.R., Jameson R.F., 1999, MNRAS, 303, 835

Heggie D.C., Aarseth S.J., 1992, MNRAS, 257, 513

Henney W.J., O’Dell C.R., 1999, AJ, 118, 235

Hillenbrand L.A., 1997, AJ, 113, 1733

Hillenbrand L.A., Hartmann L.W., 1998, ApJ, 492, 540

Hills J.G., 1980, ApJ, 225, 986

Hurley J.R., 2000, PhD Thesis, Trinity College, University of Cambridge

Hurley J.R., Pols O.R., Tout C.A., 2000, MNRAS, 315, 543

Hurley J.R., Tout C.A., Aarseth S.J., Pols O.R., 2000b, MNRAS, in prep.

Jones B.F., Walker M.F., 1988, AJ, 95, 1755

Klessen R.S., Burkert A., Bate M.R., 1998, ApJ, 501, L205

Kähler H., 1999, A&A, 346, 67

Köhler R., Leinert Ch., 1998, A&A, 331, 977

Kroupa P., 1995a, MNRAS, 277, 1491

Kroupa P., 1995b, MNRAS, 277, 1507

Kroupa P., 1995c, MNRAS, 277, 1522

Kroupa P., 2000a, NewA, 4, 615

Kroupa P., 2000b, in ASP Conf. Ser. Vol. 211, Massive Stellar Clusters, ed. A. Lancon, C. Boily (San Francisco: ASP), p.233 (astro-ph/0001202)

Kroupa P., 2000c, MNRAS, in press

Kroupa P., Petr M.G., McCaughrean M.J., 1999, NewA, 4, 495 (KPM)

Lada E.A., 1999, in The Origin of Stars and Planetary Systems, ed. C.J. Lada, N.D. Kylafis, NATO Science Series C, Vol 540, (Dordrecht: Kluwer), p.441

Lada C.J., Margulis M., Dearborn D., 1984, ApJ, 285, 141

Lindegren L., Perryman M.A.C., 1996, A&AS, 116, 579

Luhman K.L., Rieke G.H., 1999, ApJ, 525, 440

Martin E.L., Brandner W., Bouvier J., et al., 2000, ApJ, 543, 299

Mason B.D., Gies D.R., Hartkopf W.I., Bagnuolo W.G., et al., 1998, AJ, 115, 821

Mathieu R.D., 1983, ApJ, 267, L97

Mathieu R.D., 1994, ARA&A, 32, 465

McCaughrean M.J., Stauffer J.R., 1994, AJ, 108, 1382

Megeath S.T., Herter T., Beichman C., Gautier N., et al., 1996, A&A, 307, 775

Mermilliod J.-C., Rosvick J.M., Duquennoy A., Mayor M., 1992, A&A, 265, 513

Meusinger H., Schilbach E., Souchay J., 1996, A&A, 312, 833

Mikkola S., Aarseth S.J., 1993, Cel. Mech. Dyn. Astron., 57, 439

Palla F., Stahler S.W., 1999, ApJ, 525, 772

Palla F., Stahler S.W., 2000, ApJ, 540, 255

Petr M.G., 1998, Binary Stars in the Orion Trapezium Cluster: A High angular Resolution Near-Infrared Imaging Study, PhD thesis, University of Heidelberg

Pinfield D.J., Jameson R.F., Hodgkin S.T., 1998, MNRAS, 299, 955

Pinto F., 1987, PASP, 99, 1161

Portegies Zwart S.F., McMillan S.L.W., Hut P., Makino J., 2000, MNRAS, submitted

Preibisch T., Balega Yu., Hofmann K.-H., Weigelt G., Zinnecker H., 1999, NewA, 4, 531

Prosser C.F., Stauffer J.R., Hartmann L.W., Soderblom D.R., Jones B.F., Werner M.W., McCaughrean M.J., 1994, ApJ, 421, 517

Raboud D., Mermilliod J.-C., 1998, A&A, 329, 101

Richichi A., Leinert Ch., Jameson R., Zinnecker H., 1994, A&A, 287, 145

Röser S., 1999, Reviews in Modern Astronomy, 12, 79

Saiyadpour A., Deiss B.M., Kegel W.H., 1997, A&A, 322, 756

Scally A., Clarke C., McCaughrean M.J., 1999, MNRAS, 306, 253

Selman F., Melnick J., Bosch G., Terlevich R., 1999, A&A, 347, 532

Terlevich E., 1987, MNRAS, 224, 193

Verschueren W., 1990, A&A, 234, 156

Verschueren W., David M., 1989, A&A, 219, 105

Whitworth A., 1979, MNRAS, 186, 59

Wilson T.L., Filges L., Codella C., Reich W., Reich P., 1997, A&A, 327, 1177

Zinnecker H., McCaughrean M.J., Wilking B.A. 1993, in Protostars and Planets III, eds. E.H. Levy, J.I. Lunine, (Tucson: Univ. Arizona Press), p.429