Dynamical Evolution of Young Clusters

Dynamical Evolution of Young Embedded Clusters: A Parameter Space Survey

Eva-Marie Proszkow Michigan Center for Theoretical Physics, Physics Department, University of Michigan, Ann Arbor, MI 48109; emproszkow@gmail.com Fred C. Adams Michigan Center for Theoretical Physics, Physics Department, University of Michigan, Ann Arbor, MI 48109; fca@umich.edu
Astronomy Department, University of Michigan, Ann Arbor, MI 48109
Abstract

This paper investigates the dynamical evolution of embedded stellar clusters from the protocluster stage, through the embedded star-forming phase, and out to ages of 10 Myr — after the gas has been removed from the cluster. The relevant dynamical properties of young stellar clusters are explored over a wide range of possible star formation environments using -body simulations. Many realizations of equivalent initial conditions are used to produce robust statistical descriptions of cluster evolution including the cluster bound fraction, radial probability distributions, as well as the distributions of close encounter distances and velocities. These cluster properties are presented as a function of parameters describing the initial configuration of the cluster, including the initial cluster membership , initial stellar velocities, cluster radii, star formation efficiency, embedding gas dispersal time, and the degree of primordial mass segregation. The results of this parameter space survey, which includes simulations, provide a statistical description of cluster evolution as a function of the initial conditions. We also present a compilation of the FUV radiation fields provided by these same cluster environments. The output distributions from this study can be combined with other calculations, such as disk photoevaporation models and planetary scattering cross sections, to ascertain the effects of the cluster environment on the processes involved in planet formation.

open clusters and associations: general –â stellar dynamics — stars: formation â– planets: formation

1 Introduction

The formation of stars and planets constitutes a fundamental problem in astrophysics. Although a working theory of star formation has been constructed over the past two decades (e.g., Shu et al., 1987), much of the theoretical development applies specifically to the formation of isolated stars. In contrast, recent observational work underscores the fact that most star formation takes place in embedded stellar groups and clusters (e.g., Lada & Lada, 2003; Porras et al., 2003; Megeath et al., 2004; Allen et al., 2007). Given that most stars form in clusters, we are faced with the following overarching question: If stars form in clusters, how does the background cluster environment affect star formation and the accompanying process of planet formation? The fundamental goal of this paper is to help address the second portion of this question.

In rough terms, the cluster environment can influence star and planet formation through two channels: through direct dynamical effects (e.g., scattering events between cluster members) and through the background radiation fields produced by massive stars in the cluster. For the latter issue, the massive stars tend to reside in the cluster center, so that radiation exposure depends on the radial locations of the constituent solar systems, and these locations are determined dynamically. As a result, both scattering interactions and radiation exposure depend on the dynamical evolution of the cluster. To investigate these issues, this paper presents an extensive parameter space study of the dynamics of young embedded clusters spanning a wide range of initial conditions and other properties.

Star-forming clusters can be viewed in two ways. One can consider the cluster itself as an astrophysical object, and study its properties as a function of time. For example, we can track the fraction of stars retained, the half-mass radius, the virial parameter, and other variables as the cluster is born, lives, and dies. On the other hand, we can focus on the effects of the cluster environment on its constituent solar systems. Although this paper presents results relevant to both points of view, we concentrate on the latter.

This paper will focus on clusters with membership sizes in the range = 100 – 3000. The current observational surveys in the solar neighborhood indicate that the majority of star formation takes place within clusters within this size range (Lada & Lada, 2003; Allen et al., 2007). Large clusters () are relatively rare and their dynamics are well-studied (e.g., Portegies Zwart et al. 1999; see also Heggie & Hut 2003). Small systems () do not have a large impact on star/planet formation (Adams & Myers, 2001), except for few-body effects that have already been considered (Sterzik & Durisen, 1998). This paper thus works in an intermediate regime of parameter space, between the extremes where clusters are highly disruptive (e.g., Bate et al., 2003) and more isolated cases where individual collapse events take place unimpeded.

This parameter space survey will perform dynamical calculations spanning the first 10 Myr of cluster evolution. This time scale is comparable to the typical lifetime for gaseous disks (Hernández et al., 2007), the time required to form gas giant planets (Lissauer & Stevenson, 2007), and the lifetime of massive stars. In addition, the gaseous component of embedded young clusters is observed to disperse in only 3 – 5 Myr (e.g., Allen et al., 2007), so that the clusters expand appreciably by the time they reach an age of 10 Myr (Bastian et al. 2008), and this expansion leads to reduced interaction rates (Adams et al., 2006, hereafter APFM). As a result, embedded clusters exert their greatest effects on forming solar systems during their first 10 Myr. Notice also that protostellar collapse takes place over a much shorter time scale, only about 0.1 Myr, so that the cluster environment has relatively little time to affect the star formation process per se. The most important effects of the background environment act on circumstellar disks and planet formation, i.e., forming solar systems.

We note that a great deal of dynamical work on young clusters has been done previously (e.g., Lada et al., 1984; Rasio et al., 1995; Kroupa, 1995; Boily & Kroupa, 2003a, b, and many others). In addition to expanding the parameter space under consideration, this work differs from previous dynamical studies of young clusters in several ways. The clusters considered here, with = 100 – 3000, are highly chaotic so that different, but equivalent, realizations of the system can produce different dynamical results (e.g., different numbers of close approaches). In addition, the sampling of the stellar initial mass function is not complete for these intermediate-sized clusters. To address these issues, we must adopt a statistical approach. In order to characterize the dynamics of these systems, we run multiple equivalent realizations of the simulations to build up robust distributions of the output measures, e.g., the distributions of closest approaches and distributions of radial positions (which determine radiation exposure). Previous work (APFM, Proszkow et al. 2009; Proszkow 2009; see also Malmberg et al. 2007) indicates that for these moderate-sized clusters and intermediate time scales one needs realizations for each set of initial conditions to obtain robust statistics. As a result, this survey of parameter space includes the results from numerical simulations.

The initial conditions for the -body simulations of this paper also differ from most previous studies. Past studies often assume that the initial phase space variables of the stars are close to virial equilibrium. In many regions, however, pre-stellar clumps are observed to move subsonically before the clumps form stars (Peretto et al. 2006; Walsh et al. 2004; André 2002); these data imply that newly formed stars begin their dynamical evolution with subvirial speeds. Recent results from the Spitzer Space Telescope show that the class I objects (protostars) are segregated from the class II objects (star/disk systems) in embedded clusters, but their positions are highly correlated with each other and with the gas ridges (T. Megeath, private communication); this result suggests that the protostars are not moving in a fully dynamical manner. As shown previously (APFM), and reinforced in this study, subvirial starting states have a significant impact on the resulting cluster properties.

This paper is organized as follows. In Section 2, we outline the parameter space that the -body simulations explore, including membership size , stellar profiles, initial velocities, radial sizes, star formation efficiency, gas removal timescale, and the degree of primordial mass segregation. The results from these simulations are then presented in Section 3, including both the parameters that describe cluster properties (e.g., fraction of bound stars, mass distributions, and number density distributions) and parameters that affect solar system disruption (e.g., rate of close encounters). The effects of these cluster environments on the solar systems forming within them is addressed in Section 4, including a discussion of the radiation fields produced by these cluster systems. We conclude in Section 5 with a summary of our results and a discussion of their implications.

2 Numerical Simulations of Young Embedded Clusters

Throughout this work, a modified version of the NBODY2 code developed by Aarseth (2001) is employed to numerically calculate the dynamics in young stellar clusters from the embedded stage out to ages of . The modifications made to this code allow for the cluster’s initial conditions to be more like those observed in young stellar clusters. More specifically, these modifications allow us to specify the form and time evolution of the embedding gas, provide differing degrees of initial mass segregation, and define the geometry and velocity structure of the stellar distribution. Additional modifications are implemented to produce the output parameters of interest in this study, including closest approach distributions, velocity distributions, and mass profiles.

In addition to focusing on clusters with initial conditions similar to those found in nearby young clusters, this work is distinguished by its statistical character. The -body problem is by nature chaotic, so that clusters with the similar initial configurations can produce dissimilar results. In order to produce robust statistical descriptions of cluster evolution, we must carry out realizations of the cluster simulations for each initial condition configuration considered in this study. Specifically, for a given set of initial cluster conditions (i.e., cluster membership, radius, velocity distribution, etc.), simulations are completed using a different random number seed to sample the relevant distributions (e.g., stellar positions, IMF). The resulting output parameters are then averaged over the set of realizations to provide a statistical description of how a similar cluster is likely to evolve. In this context, a “similar cluster” is one produced by an independent realization of the initial conditions, which are sampled from the same distribution of values.

We note that the version of the -body code used for these simulations (NBODY2) does not include binaries in the initial conditions (see Aarseth 1999, 2001); we also use softening parameter = 0.001 (which under-resolves binaries that may form during the course of cluster evolution). With this value of , the code resolves stellar encounters down to AU, the smallest approach distance considered here. Since this integration package is relatively fast, it can produce many realizations for each set of initial conditions, as required here to obtain good statistics. In dense and/or long-lived clusters, binaries can affect the energetics of the system by absorbing and storing energy in their orbits. In the systems of interest here, however, the densities are low and the evolutionary times are short, so that binaries have relatively little impact on cluster evolution. In other words, interactions are sufficiently rare and sufficiently distant, so that binarity has a relatively small effect on overall energy budget (for further discussion, see Kroupa 1995, Kroupa et al. 2003). In our previous work (APFM), we checked this approximation for consistency by running a parallel set of simulations including binaries (using the NBODY6 integration package; Aarseth 1999) and found that the results were essentially unchanged. This same test case implies that our choice of softening parameter does not greatly affect the global results. As an additional check on binary effects, we can use the distributions of closest approaches resulting from our ensemble of simulations; these results also indicate that that binary interactions are not energetically important.

In this section we outline the standard initial conditions used in the simulated clusters. Specifically, we discuss the qualities most commonly observed in nearby young embedded clusters and identify these qualities as the center of our parameter space. These initial conditions thus define the typical cluster, and the parameter space survey is conducted by varying one or more of the initial conditions at a time. In the following discussion, in Section 2.1, we enumerate the variables of interest. The particular range of parameter space investigated in this survey is then outlined in Section 2.2. The results of the cluster simulations are presented in Section 3, where we also discuss the implications for planet formation within these clusters.

2.1 Specification of Input Variables

Cluster Membership, . As outlined above, in this study we consider intermediate-sized clusters with stellar memberships ranging from to (Lada & Lada, 2003; Allen et al., 2007). This range corresponds to that observed in the solar neighborhood, where we have the most direct window into the star formation process.

Initial Mass Function. The shape of the stellar initial mass function (IMF) observed in young stellar clusters is almost universal for clusters with more than members (Lada & Lada, 2003). In our simulations, stellar masses are sampled from the log-normal analytic fit to the standard IMF of Miller & Scalo (1979) presented by Adams & Fatuzzo (1996). With this choice, the average stellar mass in a cluster is  M(the average stellar mass is somewhat higher than the median stellar mass which is roughly  M), consistent with observations of young stellar clusters (Muench et al., 2002; Luhman et al., 2003).

Cluster Radius, . Stars are initially distributed within the cluster radius according to the density profiles described below. The cluster radius is taken to be a function of the cluster membership , the scaling radius , and the power-law index , so that

 Rc=Rsc(N300)α. (1)

This membership-radius relationship is observed in young clusters in the solar neighborhood (Lada & Lada, 2003; Porras et al., 2003), and typical values of the parameters are = and = . Thus, a cluster with stars typically has a radius of .

Initial Stellar Profile. Many young embedded clusters display degrees of central concentration (see Lada & Lada 2003, and references therein). The simulated clusters in this study are correspondingly centrally condensed and have initial stellar density distributions of the power-law form . Density profiles of this form are consistent with the observed density profiles of the embedding gas in cluster-forming cores (Jijina et al. 1999; see also below). For simplicity, we take the initial profiles to be spherically symmetric. However, we note that irregular initial distributions (Goodwin & Whitworth 2004) can have interesting effects, e.g., accelerated mass segregation.

Mass Segregation. Young stellar clusters exhibit varying degrees of mass segregation, even though the clusters themselves are not old enough to have undergone dynamical mass segregation (their age is generally less than a relaxation time, e.g., Bonnell & Davies 1998). However, recent work shows that mass segregation can occur more rapidly when the clusters have subvirial initial conditions and clumpy substructure (McMillan et al. 2007, Allison et al. 2009, Moeckel & Bonnell 2009). For most of this survey, the simulated clusters contain minimal mass segregation implemented by a straightforward algorithm: At the start of the simulation, the most massive star in the cluster is relocated to the center of the cluster, while the positions of the rest of stars are not correlated with mass.

Initial Stellar Velocities. Kinematic observations of the youngest stellar objects (Class I sources) and of starless dense cores indicate that these objects are moving at speeds that are a fraction of the virial speed in many young clusters (Walsh et al., 2004; André, 2002; Peretto et al., 2006; Kirk et al., 2006). In the cluster simulations, initial stellar velocities are sampled from a spatially isotropic distribution and then scaled by the initial virial ratio of the cluster . The virial ratio is defined as , where is the total kinetic energy and is the total potential energy of the cluster. A cluster that is in virial equilibrium has a virial ratio . Most simulations considered in this study are initialized with a virial ratio which results in stellar velocities that are approximately one-third of the virial velocity of the cluster and is consistent with the kinematic observations of young embedded clusters.

Star Formation History. The stars in the simulated clusters have a spread in formation times of . The formation time of each star is sampled from a uniform distribution over the range from to , independent of position within the cluster and independent of stellar mass. We assume that the forming stars are tied to their formation site (they do not move dynamically) until the collapse phase is complete, i.e., until the star is formed. The stars are included in the simulations as static point masses until their formation time. After they form, the stars are allowed to move through the gravitational potential of the cluster with an initial velocity sampled from the distribution described above.

Embedding Gas Profile. Extremely young stellar clusters (with ages less than Myr) are almost always associated with a molecular cloud core (Leisawitz et al., 1989). These cores are often centrally concentrated (Larson, 1985; Myers & Fuller, 1993; Jijina et al., 1999). In the simulated clusters this embedding gas is represented as a static gravitational potential with a Hernquist profile (Hernquist, 1990), with potential, density, and mass profiles of the form

 Ψ=2πGρ0r2s1+ξ,ρ=ρ0ξ(1+ξ)3,andM=M∞ξ2(1+ξ)2, (2)

where . Here the parameter is a scale length, which is chosen to be equal to the cluster radius, i.e., . In the inner limit, this profile has the form of , and outside of the cluster the density profile matches onto a force-free background. We are thus neglecting external forces on the cluster (e.g., galactic tides).

Star Formation Efficiency. Estimates of the star formation efficiency in young star-forming regions vary from to (Lada & Lada, 2003). In our simulated clusters, a standard star formation efficiency of is assumed. This value corresponds to a total stellar mass in the cluster that is one half of the mass of the embedding gas (from equation [2], is the effective gas mass within the cluster radius ).

Gas Removal History. Although young stellar clusters are associated with embedding molecular gas, the gas is quickly dispersed from the cluster by a collection of processes, including stellar winds from young stars, ionizing radiation from massive stars, and other processes. Clusters with ages greater than are rarely associated with molecular gas. In the cluster simulations, the depth of the potential well associated with the embedding gas is reset to zero instantaneously at a time denoted as . The gas removal mechanism is thus assumed to rapidly disperse gas from the vicinity of the cluster at this epoch.

Stellar Evolution. This parameter space survey considers only the first 10 Myr of cluster evolution and does not include stellar evolution effects. We expect these corrections to be relatively small for the cases considered here: Note that a star with initial mass is expected to burn hydrogen for 8.13 Myr and then burn helium for another 1.17 Myr (Woosley et al. 2002). The total time before exploding as a supernova is thus about 9.3 Myr. As a result, only those stars with masses will finish their evolution during the 10 Myr window considered here. And only about 1 star in 800 (or 1000, depending on the IMF) has a starting mass this large. Stellar evolution thus has a small effect on the mass budget of the clusters. Although supernovae would have a large impact on the gas content of clusters, the gas is removed earlier, as indicated by observations (and presumably driven by other mechanisms).

2.2 Parameter Space Overview

As described in Section 2.1, a large number of initial parameters must be specified in order to characterize a cluster at the start of a simulation. As a result, the parameter space available for studying the evolution of embedded stellar clusters is extremely large. In this current work, we target our parameter space survey on embedded cluster environments similar to those observed in our solar neighborhood (Lada & Lada, 2003; Allen et al., 2007; Megeath et al., 2004, Gutermuth et al. 2009, in preparation), with an extrapolation to somewhat larger clusters. In this section we identify the range of parameter space for which our survey is conducted. As described below, we perform many different series of simulations, where each series explores the effects of one (or more) specific parameters(s). It is important to note that this range, while motivated by observations of nearby clusters, does not necessarily encompass all of the possible initial conditions spanned by these these cluster environments. The range of parameter space surveyed and the initial conditions assumed in the simulated clusters are summarized in Table 1.

2.2.1 Cluster Membership

We perform a series of simulations to study the effect that stellar membership has on the dynamics of young embedded clusters. We consider spherical clusters embedded in centrally concentrated gas potentials with a star formation efficiency . The stellar membership in the simulated clusters ranges from to . Clusters of this size roughly span the range of young clusters observed in the solar neighborhood (Lada & Lada, 2003; Porras et al., 2003). Motivated by observations of young stellar objects with subvirial velocities, this study considers embedded clusters with both subvirial and virial initial velocity distributions. The subvirial and virial clusters have and respectively. The subvirial clusters thus have initial stellar velocities that are approximately of the virial velocity.

It is possible that the index appearing in the membership-radius relation (equation [1]) takes on different values for different cluster samples. For example, the value is a reasonable fit to the observed data within approximately of the Sun (Lada & Lada, 2003; Porras et al., 2003), where this sample contains intermediate-sized clusters with . In environments with star formation rates much higher than that of the solar neighborhood, a significant amount of activity occurs in clusters more massive than those found in our solar neighborhood (e.g., Chandar et al. 1999, Pfalzner 2009). These extremely massive young clusters, some which are thought to be progenitors of globular clusters, contain as many as stars and have sizes on the order of (Mengel et al., 2008). If we extend the cluster membership-radius relation out to stellar memberships as high as , the choice of would overestimate the cluster radius by a factor of . A power-law index of more closely reproduces the observed data points over the full range of . In this study, we investigate the evolution of intermediate-sized clusters using both and power-law indices in the cluster membership-radius relation. In both cases, we chose so that the power-law passes through the point where and .

These two choices of the index result in clusters whose average number density varies differently as a function of cluster membership. Specifically, substituting the membership-radius relation into the equation for average number density gives the relation

 n0∼NR3c∼NN3α=N1−3α. (3)

For the choice the average stellar density decreases as a function of , whereas for the stellar density is an increasing function of . In the results summarized in Section 3, many of the trends observed as a function of cluster membership are more fundamentally trends in average stellar density as a function of .

Note that an intermediate value of the index, , implies a constant stellar density. This benchmark density value is stars . Figure 1 displays the average number densities found in clusters in the solar neighborhood. The data are taken from the cluster catalogs of Lada & Lada (2003) (diamonds) and Carpenter (2000) (triangles). Number densities are calculated assuming spherical symmetry in the stellar clusters. Nearby young clusters may have higher densities in the cluster cores (Hillenbrand & Hartmann, 1998; Gutermuth et al., 2005; Teixeira et al., 2006), but their average stellar densities are relatively constant. The horizontal line in the figure shows a constant density reference with the median value of the data set, = 65 pc.

2.2.2 Initial Virial Parameter

As discussed in Section 2.1, recent observations of young embedded clusters indicate that stars are formed with initial velocities lower than the virial speed of the cluster. During the early evolution of a subvirial cluster, the average stellar velocities increase as individual stars fall through the global potential well of the cluster. Stars with initially subvirial velocities thus trade potential energy for kinetic energy during the early phases of cluster evolution. On somewhat longer timescales, interactions share energy among stellar orbits and the cluster approaches virial equilibrium. Here we present a series of numerical experiments designed to investigate the effect that the initial virial ratio has on the evolution of clusters. While the simulations completed as a part of the cluster membership parameter study make a rough comparison between subvirial and virial clusters, this set of simulations samples a much wider range of virial ratios with much higher resolution. Simulations are completed for clusters with initial membership , and that have a membership-radius relation characterized by , similar to that observed in the solar neighborhood.

One question that this study attempts to address is: How small must the initial virial parameter be in order for cluster evolution to differ significantly from that of a cluster in virial equilibrium? Our results indicate that even moderately subvirial clusters display characteristics significantly different from virial clusters (see Section 3). For instance, the bound fraction in a cluster with initial virial parameter is almost larger than in a virial cluster with . A value of corresponds to an average stellar velocity that is approximately of the virial velocity of the cluster. As a result, initial stellar velocities can be an appreciable fraction of the virial value and still lead to significant differences.

The cluster membership-radius relation presented in equation (1) depends on both the power-law index and the fiducial scaling parameter that sets the radius of a cluster with members. Although this scaling relationship is robust, the radius of observed clusters still contains significant scatter as a function of stellar membership (see Figure 1 of APFM). Some of this scatter results from the observational difficulty of determining the outer radius of a cluster as the surface density approaches that of the background sky. In addition, it is difficult to determine cluster radii for non-spherical clusters and clusters with small memberships (Gutermuth et al., 2005; Allen et al., 2007).

A series of simulations are completed to investigate how cluster evolution varies with different values of scaling parameter . Specifically, cluster simulations are completed for scaling radii in the range , where we use power-law index in the membership-radius relation. The clusters are assumed to have stellar memberships , and and subvirial initial velocities with . Changing the scaling radius effectively changes the average stellar density in a cluster, and many of the trends observed in the cluster evolution as a function of scaling radius are linked to this change in density.

2.2.4 Star Formation Efficiency

The star formation efficiency (SFE) of a region is defined as , where and are the total stellar and gaseous mass contained in the region, respectively. The stellar mass thus corresponds to the value after the interval , i.e., after star formation has been completed. Estimates of the SFE for young embedded clusters in the solar neighborhood range between and (Lada & Lada, 2003; Allen et al., 2007). These efficiencies are significantly higher than the SFEs of entire giant molecular clouds, which have typical SFEs less than 0.05 (Duerr et al., 1982; Evans & Lada, 1991). As a part of this parameter space study, we complete a suite of cluster simulations in which the SFE is varied over the range . The clusters are assumed to be in an initially subvirial state, and have stellar memberships of = 300 and 1000. The value of the star formation efficiency parameter is attained by varying the mass of the gas in the simulated cluster (for fixed ). The SFE of a cluster is a major factor in determining the fraction of the stars in a cluster that remain bound after the embedding gas is dispersed (due to outflows from young stars or ionizing radiation from the most massive star; see Section 3).

2.2.5 Gas Removal Timescale

Although the youngest star-forming clusters are deeply embedded in their natal molecular clouds, clusters with ages greater than Myr are associated with relatively little molecular gas (Leisawitz et al., 1989). In this series of simulations, we consider the evolution of young clusters with gas removal times in the range from to . The gas is assumed to be removed instantaneously at . The study considers subvirial clusters with stellar memberships , and . A significant fraction of the stars in a given cluster become gravitationally unbound at the time of gas dispersal and the cluster begins to expand radially outward. As the cluster expands, the average density decreases and close interactions between stellar members become less frequent. As a result, in addition to affecting how much of the cluster remains gravitationally bound, the gas removal time places important constraints on the close encounter rates in young clusters.

2.2.6 Mass Segregation

Observations indicate that massive stars are preferentially found near the center of both evolved open clusters and young embedded clusters. Mass segregation in the evolved clusters can be explained by dynamical theory: high mass stars lose energy to low mass stars through two-body interactions and subsequently sink toward the cluster’s center. This process takes place on timescales comparable to a cluster’s dynamical relaxation time, which is given by

 trelax≈N8 ln Nτc≈NRc8 ln N⟨v⟩, (4)

where is the cluster crossing time and is the average stellar velocity (Binney & Tremaine, 1987). Open clusters have typical ages 20 – 500 Myr and thus are old enough for dynamical mass segregation to have occurred. However, observations of mass segregation in young embedded clusters are more difficult to explain (Bonnell & Davies, 1998). A (logarithmically) average embedded cluster in the solar neighborhood has , , , and a corresponding relaxation time of roughly . As a result, dynamical evolution is unlikely to be responsible for the mass segregation observed in young clusters such as the Trapezium, NGC 2071, or NGC 2074 (Hillenbrand & Hartmann, 1998; Lada et al., 1991; Bonnell & Davies, 1998). This finding suggests that the mass segregation is due to a primordial tendency to form massive stars near the center of clusters. We note that primordial mass segregation is naturally produced in embedded clusters through some proposed massive star formation scenarios. For example, competitive accretion preferentially forms massive stars in the deepest part of the cluster potential well, near the center of the cluster (Bonnell et al., 2001; Beuther et al., 2007).

One experiment in the parameter space survey explores the evolution of clusters with varying amounts of primordial mass segregation. We define the primordial mass segregation parameter as the fraction of the cluster membership which has been ordered by mass at the center of the cluster, . More specifically, at the start of the simulation, the stellar masses are sampled from a standard IMF and stellar positions are sampled from a density profile (regardless of mass). After this initial sampling, the most massive star in the cluster is moved from its initial (randomly assigned) position to the center of the cluster. This resulting state represents a cluster with minimal mass segregation. For values of , additional mass segregation is implemented by rearranging the stellar positions so that the most massive stars are located at the inner radial positions. The mass segregation parameter is varied over the range in subvirial clusters with = 300, 1000, and 2000 members.

3 Summary of Results

This parameter space survey includes a large number (25,000) of -body simulations. In this section, we package the results of these numerical calculations by extracting relevant output variables. The fraction of stars that remain bound to the cluster is discussed in Section 3.1. Next we determine the interaction rates between cluster members (Section 3.2) and the distributions of interaction speeds (Section 3.3). Finally, we construct mass profiles and number density profiles in Section 3.4.

3.1 Bound Fraction

Observational studies that compare the formation rates of embedded clusters and open clusters have shown that the embedded cluster formation rate is significantly higher than that of open clusters (Elmegreen & Clemens, 1985; Battinelli & Capuzzo-Dolcetta, 1991; Piskunov et al., 2006). This discrepancy in the formation rates leads to the conclusion that while most star formation occurs within clusters, only a fraction (about ) of the stellar population is born within “robust” clusters that are destined to become open clusters (which then live for or longer). This result suggests that very few embedded clusters remain gravitationally bound after their embedding molecular gas is removed. The process by which gas removal leads to the unbinding of a cluster has been denoted as “infant mortality” for embedded clusters, and has been addressed via both analytical studies (Hills, 1980; Elmegreen, 1983; Verschueren & David, 1989; Adams, 2000) and numerical simulations (Lada et al., 1984; Geyer & Burkert, 2001; Boily & Kroupa, 2003a, b). Evidence that this process is occurring in extragalactic young massive clusters has been presented by Bastian & Goodwin (2006).

An important output parameter explored in our simulations is the fraction of stars that remain gravitationally bound as a function of time. The bound fraction is defined as where is the initial stellar membership, and is the number of stars that have total energy (kinetic plus potential) less than zero. Throughout the embedded phase of cluster evolution, the bound fraction remains equal to unity (note that the timescale of the embedded phase is shorter than the relaxation time). The embedding gas potential is removed from the simulated clusters instantaneously at time . This event significantly reduces the depth of the potential well in which the cluster members reside. Rapid gas removal is an appropriate approximation for gas expulsion due to high mass star formation, which removes the embedding gas over timescales as short as years (Whitworth, 1979).

When the gravitational potential of the gas is removed from the system, the high-velocity stars become gravitationally unbound while the low-velocity stars remain bound to the cluster’s gravitational potential. As a result, the bound fraction decreases significantly (by as much as ) over a short period of time. The fraction then continues to decrease (more slowly) until the end of the simulations. Note that our temporal cutoff is chosen to be , but the clusters will continue to evolve and will continue to decrease on longer timescales.

Figures 2 and 3 display as a function of time for the range of parameter space surveyed in this study. Each panel illustrates the temporal evolution of for a specific cluster parameter, where the individual curves correspond to different values of that cluster parameter. For example, in Figure 2 panel (a), the evolution of is plotted for clusters with different star formation efficiencies . The top curve corresponds to a cluster with and the bottom curve corresponds to a cluster with . For each of the curves in Figures 2 and 3, the other initial conditions (i.e., cluster membership , initial cluster radius , degree of mass segregation , etc.) are held constant. All of the panels show a rapid decrease in at (except for Figure 3 panel (a), which represents simulations in which the gas removal timescale is varied). As expected, the downward jump in corresponds to the time at which the gas is removed from the cluster.

The value of the bound fraction at provides one measure of how tightly bound a cluster remains after the embedding gas is removed. Figures 4 and 5 display the value of the bound fraction at as a function of the initial cluster parameter values for the range of parameters considered in this survey. Previous theoretical and numerical work has identified a cluster’s star formation efficiency as the most important parameter in determining whether or not a cluster will remain gravitationally bound (Hills, 1980; Elmegreen, 1983; Lada et al., 1984). In clusters with high SFEs, a large proportion of the total cluster mass remains behind (in the form of stars) after the embedding gas is removed. Clusters with high SFEs remain more tightly bound after gas dispersal than clusters with low SFEs. In our cluster parameter survey, we also find that the bound fraction at depends sensitively on the star formation efficiency of the cluster. Figure 4 panel (a) displays the cluster bound fraction as a function of star formation efficiency, . The data is well fit by a power-law in :

 fb=2.23(εSF)1.2whereεSF≤0.5. (5)

This fit is shown as the solid curve (line) in the figure.

In the suite of simulations used to investigate the effects of star formation efficiency, the clusters are initially subvirial. After gas removal, subvirial clusters are more tightly bound than virial clusters (APFM). For even relatively high star formation efficiencies () and small initial virial parameters, conditions which produce the most tightly bound systems, the clusters are significantly disrupted by gas removal and promptly lose of their stars. Star formation efficiencies larger than are rarely observed (Lada & Lada, 2003, and references therein), and are difficult to attain theoretically (Matzner & McKee, 2000). Our results are thus in agreement with previous studies indicating that a significant fraction of the stellar population is lost from a cluster during the gas removal phase (e.g., Lada et al. 1984, Adams 2000, Boily & Kroupa 2003ab).

As mentioned above, clusters with subvirial initial velocities are more tightly bound than clusters with virial initial conditions. Figure 4 panel (b) demonstrates this trend by plotting the bound fraction as a function of the initial virial ratio. The bound fraction decreases steadily with the initial virial ratio over the range considered, . Gas removal has a weaker effect on spherical clusters with subvirial initial conditions because as a subvirial cluster collapses, more of the stars spend more of their time inside of the embedding gas (which is assumed to be static, i.e., not in a state of global collapse). When the embedding gas is removed from the cluster, many of the cluster members are interior to the gas and are less affected by the change in potential. The results from this suite of simulations again indicate that a significant fraction of cluster members are lost due to the change in the gravitational potential that occurs during the dispersal of the natal gas. In the most tightly bound subvirial clusters with , approximately of stars become unbound due to dispersal of the gas.

In clusters with subvirial initial conditions, the bound fraction remains constant as a function of cluster membership , for both the and the cluster membership-radius scaling relations. This finding also holds true for virial clusters that have cluster membership-size relations similar to those observed in the solar neighborhood (). The bound fraction at is plotted as a function of cluster size in panels (c) and (d) of Figure 4. The upper curves in these panels indicate the bound fraction at in the more tightly bound subvirial clusters, whereas the lower curves correspond to the virial clusters.

In virial clusters with a lower power-law index in the cluster membership-size relation (), the bound fraction decreases as a function of the cluster membership . The bound fraction decreases roughly as (see Figure 4, panel (c), lower curve). This decrease in is due to a combination of effects arising from the relationship between cluster radius and cluster membership defined by equation (1). In clusters with , the mean velocity and velocity dispersion scale approximately as . The velocity distributions in the clusters are nearly Gaussian during the embedded phase (rather than perfectly Maxwellian as would be expected in a collisionless isothermal sphere of stars). As a result, the increased velocity dispersion in clusters with larger results in more stars with velocities high enough to escape from the cluster . In addition, the interaction rate between cluster members increases with the stellar density, which increases as in these clusters. In virial clusters with (Figure 4, panel (d), lower curve), the bound fraction is roughly constant. This trend occurs because although the average velocity and velocity dispersions increase as a function of cluster membership, the dependence on is not as strong: . In addition, the stellar density actually decreases with , ; the competing effects of increased velocities and lower interaction rates are comparable and act to cancel each other out.

The bound fraction does not appear to be simply related to either the gas removal timescale or the scaling radius , i.e., the bound fraction is not a monotonic function (see Figure 5, panels (a) and (b)). This result occurs because in subvirial clusters, such as the ones considered in these parameter space surveys, changing either the scaling radius or the gas removal time affects the relationship between the gas removal time and the initial collapse and relaxation time. The resulting is sensitive to the particular dynamical state of the system at the time of gas removal. For example, if the cluster is re-expanding (after its initial collapse) when the gas is removed, many stars will have trajectories that are directed radially outward and are thus more likely to become gravitationally unbound.

Next, we consider the effects of primordial mass segregation. In general, the effects of mass segregation saturate when more than approximately of the stars are segregated by mass. Mass segregation only slightly affects the bound fraction, and clusters with minimal mass segregation (where the largest star is located at the cluster center) have slightly lower bound fractions than clusters with , as shown in panel (c) of Figure 5.

In summary, the results of this part of the study indicate that the star formation efficiency is the parameter that most significantly affects the fraction of a cluster that remains bound after the embedding gas is removed from the system. In addition, the initial virial state of the cluster, as well as the specific dynamical state at the time of gas dispersal, are important parameters in determining how many members remain bound to the cluster. We find that in sufficiently subvirial clusters, , the bound fraction is not a sensitive function of the initial stellar density (as indicated by the suite of simulations varying and the cluster membership-radius relations), but rather is dominated by the fact that the initial global collapse produces a cluster whose members reside interior to the bulk of the embedding gas and thus are not strongly affected by the gas removal.

Two caveats should be included in this discussion. First, current observations of young emerging clusters cannot determine whether a cluster member is gravitationally bound or unbound from its host cluster. Over the first , bound and unbound clusters are visibly similar, and the results of simulations such as those presented here are not easily compared directly to observations. Second, this parameter space study focuses on the early evolution of embedded clusters. Additional dynamical evolution of the clusters on timescales greater than will lead to even lower bound fractions at later times. As a result, the bound fractions presented in this work should be considered as upper limits on the expected bound fractions for clusters with older ages.

Finally we note that Tables 2 – 7 provide compilations of the bound fractions evaluated at time = 10 Myr for the simulations presented here. Each table lists the bound fractions as a function of a given input variable, including the stellar membership (Table 2), initial virial parameter (Table 3), cluster scaling radius (Table 4), star formation efficiency (Table 5), gas removal time (Table 6), and the mass segregation parameter (Table 7).

3.2 Stellar Interaction Rates

A significant consequence of living in high density environments, such as those found in young embedded clusters, is that close encounters with other cluster members may be relatively frequent. If these interactions are sufficiently close, they can have important ramifications for planet formation in circumstellar disks and for solar system survival. During early stages of solar system formation, encounters can disrupt protoplanetary disks and limit the mass reservoir for planet formation (Ostriker, 1994; Heller, 1993, 1995; Kobayashi & Ida, 2001). At later times, close encounters can disrupt planetary systems themselves, by significantly altering the eccentricities of planets and, in sufficiently close encounters, ejecting planets from the solar system entirely (Adams & Laughlin 2001, APFM).

Throughout the cluster simulations, close encounters with distance of closest approach less than  AU are recorded. Note that the distance of closest approach is somewhat smaller than the impact parameter of the encounter due to gravitational focusing. A cumulative distribution of the close encounters is then constructed, and the interaction rate is calculated by averaging the encounter distributions over the time span of interest (here we use the embedded phase , the exposed phase , or the entire interval). Specifically, the interaction rate is defined as the number of close encounters with distance of closest approach per star per million years. We find that the interaction rates have the form of power-laws for encounters with closest approach distances less than  AU. In other words, the interaction rate has the form

 Γ=Γ0(b1000AU)γ, (6)

where is the distance of closest approach. The fiducial interaction rate and the power-law index are fit to the cumulative closest approach distribution for each set of cluster simulations. The fiducial interaction rate corresponds to the number of encounters with impact parameter less than  AU per star per million years. The fiducial interaction rate is displayed as a function of cluster initial conditions in Figures 6 and 7.

For the parameter space considered in this study, the interaction rate depends most sensitively on a single parameter, the stellar number density . The trends observed as a function of stellar membership and cluster scaling radius (see Figure 6 panels (a) – (c)) are more fundamentally trends indicating how varies as a function of the average stellar density (keep in mind that the stellar density varies with and ).

We can understand this behavior in simple terms as follows: Consider a cluster with stars and radius . For simplicity, we ignore the difference between impact parameters and distances of closest approach . A star passing through the cluster will experience, on average, a number of close encounters with impact parameters within the range to (Binney & Tremaine, 1987), where is given by

 δn=2NR2cbδb. (7)

The crossing time in a cluster is given by , where is the average stellar velocity. As a result, the star will experience close encounters in the given annulus at the rate

 δΓ≈2NR2c(vRc)b δb=2nvb δb. (8)

The interaction rate for all encounters with impact parameter less than is then given by . Since the mean dynamical speeds are slowly varying over the regime of cluster parameter space considered in this paper, the interaction rate as claimed.

In Figure 6 panel (d), the values are plotted as a function of initial stellar density for all of the simulations with varying cluster membership and scaling radius (using the data from Figure 6 panels (a) – (c)). The plus symbols indicate the interaction rate for clusters with initially subvirial velocities, whereas the interaction rates for clusters with initially virial velocities are indicated by the x’s. Straight lines (indicating a power-law with index = ) are included in the panel; note that the numerically determined data are roughly consistent with the simple scaling relation .

Figures 8 and 9 present the values of the index that provide the best fit to the close encounter distributions as a function of the initial cluster parameters. The value of does not vary strongly as a function of the initial conditions, but rather remains in the range .

Using the simple argument constructed above, the total rate of close encounters with impact parameter less than in a cluster of membership size is given by

 Γb≈0.122N−1/2(b1000 AU)2Myr−1forRc=1.0 pc(N300)1/2, (9)

and

 Γb≈0.0095N1/4(b1000 AU)2Myr−1forRc=1.0 pc(N300)1/4. (10)

These estimates are similar to the interaction rates found for the virial clusters, although the fitted value of the index is slightly lower than (due to gravitational focusing) in the numerically determined distributions (see Figures 8 and 9), and the fiducial interaction rate is somewhat higher.

The subvirial clusters have interaction rates that are about 8 times larger than the rates for virial clusters of the same starting density (as defined by ). This trend is due to a combination of the smaller effective cluster radius that a subvirial cluster attains after initially collapsing and its higher bound fraction after gas dispersal. During the embedded phase, subvirial clusters () behave as if they have nearly zero-temperature starting states and thus collapse to roughly of their initial size. This decrease in radius corresponds to an increase in density by a factor of . In addition, subvirial clusters retain more of their members after the gas is removed from the cluster — we find that the bound fraction is times higher for the subvirial starting states. Since the close encounter profiles are averaged over the initial number of stars in the cluster, the interaction rates in subvirial clusters will be roughly 3 times higher than in virial clusters over the exposed phase of cluster evolution (due to increased stellar retention). Note that combining these two factors increases the interaction rates by . The results of the parameter survey varying also indicate that subvirial clusters have higher interaction rates. In Figure 7 panel (a), the interaction rate clearly decreases as a function of initial virial parameter .

We also find that the interaction rates are somewhat higher in clusters that have more of their massive stars residing near the cluster center (see Figure 7, panel (c)). This finding is consistent with the modeling results of the observed cluster NGC 1333 presented in APFM. The fiducial interaction rate in the simulated NGC 1333 cluster was approximately times higher than the calculated in equivalent subvirial clusters with minimal mass segregation. We suggested that the increased interaction rate was due to the primordial mass segregation observed in NGC 1333 (see Figure 15 of APFM). For comparison, subvirial clusters with = 300 stars and have fiducial interaction rates that are about times larger than those found in subvirial clusters with members and minimal mass segregation (see Figure 7, panel (c), top curve).

The average interaction rate also increases as a function of gas removal time as shown in Figure 7 panel (b). This interaction rate is averaged over the simulation time. However, the majority of close encounters occur during the embedded phase, and hence the average interaction rate increases as the length of the embedded phase increases. For clusters with embedded phases lasting more than , the rate of close encounters during the embedded phase is roughly constant. When the embedded phase lasts less than , the clusters have lower encounter rates due to lower densities during the first while the subvirial cluster is still contracting. For completeness, we note that the interaction rate does not display strong trends with varying star-formation efficiency (see Figure 7, panel (d)).

For all of the simulations discussed in this section, Tables 8 – 13 provide listings of the parameters (, ) that specify the interaction rates through equation (6). Each table provides the values of and as a function of a given input variable, including the stellar membership (Table 8), initial virial parameter (Table 9), cluster scaling radius (Table 10), star formation efficiency (Table 11), gas removal time (Table 12), and the mass segregation parameter (Table 13). These interaction rates, as determined by the values of , are one of the primary products of this investigation. They can be used to calculate interaction rates as a function of closest approach distance for a wide variety of cluster environments (see Section 4.1 below for one such application).

3.3 Interaction Velocities

In addition to constructing the distribution of closest approach distances associated with close encounters in the simulated clusters, we also determine the distribution of encounter velocities. The distribution of encounter velocities provides additional information regarding the effect that close encounters may have on the constituent solar systems. For example, the interaction cross sections depend on the encounter speeds.

We define the encounter velocity as the magnitude of the relative velocities of the stars at the moment of closest approach. We then create distributions of the frequency of encounter velocities throughout the simulations. Figure 10 presents the resulting distribution of encounter velocities for a cluster with radius , stars, and subvirial initial speeds (). A binning size of has been used to construct the histogram, and error bars are included on the distribution to indicate the dispersion within each velocity bin.

We find that the encounter velocity distribution can be approximated reasonably well by a Gaussian curve where the mean and the width (as measured by the variance ) of the Gaussian are varied to fit the encounter velocity distribution for each particular set of initial conditions. The Gaussian fits to the velocity distributions in Figure 10 are indicated by the dashed curves. We note that the Gaussian fit slightly overestimates the number of low velocity encounters with (and formally even predicts a few encounters with negative velocities). However, the general shape and width of the distribution are well represented by these gaussian forms.

In Figures 11 and 12, the mean encounter velocity is plotted as a function of the initial cluster parameter. The encounter velocity has been normalized by the mean velocity within the cluster’s half-mass radius (the regime where most of the interactions occur within the cluster). The error bars indicate the normalized width (FWHM) of the Gaussian that best fits the velocity distribution. This figure demonstrates that the normalized encounter velocity distributions do not vary strongly as a function of the initial conditions, but rather are a robust function of mean cluster velocity.

The encounter velocities are about twice the average velocity in the interactive region of the cluster. This result is roughly consistent with an analytic estimate of the relative velocities of cluster members whose velocities are sampled from a Maxwellian distribution, so that (Binney & Tremaine, 1987). The encounter velocities are somewhat larger than those predicted by this estimate (about twice the mean velocity) due in part to gravitational focusing. Note that AU is a typical encounter distance, and that the orbit speed km s, so that the gravity of the stars taking part in the encounter does matter. The numerically calculated distribution of interaction velocities includes only a subsample of the relative velocities, because only encounters with impact parameter  AU are included in the interaction velocity distribution, and this subsample is likely to have somewhat larger relative velocities.

3.4 Mass and Number Distributions

As a cluster evolves, interactions between stars, and between the stars and the background gas potential, produce a distribution of stellar positions and velocities. The distribution of stars within a cluster at a given time can be characterized by the cumulative mass distribution or the cumulative number distribution , where and are the total masses and numbers of the stars that are gravitationally bound to the cluster at time , respectively. In the simulated clusters, each of these distributions is calculated at intervals of Myr. The profiles are then averaged over the cluster lifetime and over 100 realizations of the cluster used to produce a statistical description of the mass and number profiles. We find that both of these distributions may be fit by simple functions of the form

 M(r)MT∗=(ξa1+ξa)p, (11)
 N(r)N=(ξa1+ξa)p, (12)

where and the scale radius and index are free parameters that are fit to the distributions observed in the simulated clusters. The parameter may also be varied to fit the data. We find that the choice gives the best fit for the subvirial clusters and gives the best fit for the initially virial clusters (this finding is consistent with the results found in APFM). In the series of simulations where the initial virial parameter is varied, the choice works best for simulations for which , and works best for those with . The parameters and that provide the best fit for the radial distributions (equation [12]) and the mass profiles (equation [11]) are similar, although not identical, because stars of different masses have somewhat different radial profiles.

Although one can use the full distributions (see equations [11, 12]), it is sometimes useful to characterize the distributions in terms of a single parameter. Toward this end, Figures 13 and 14 present the median cluster radius calculated from the fits to the cumulative radial distributions for the entire parameter space. The scale is defined as the radius at which and thus represents the radius which, on average, contains half of the cluster members. As these plots indicate, the median cluster radius scales with the initial virial parameter , the gas removal timescale , the cluster membership , and the cluster scaling radius . On the other hand, the radius does not vary strongly with the either star formation efficiency or the degree of primordial mass segregation .

Scaling the median radius by the initial cluster radius removes the dependency on this initial cluster parameter and more readily identifies trends that are distinct from the initial assumptions concerning cluster size. Figure 15 displays the median radius normalized by the initial cluster radius . Panel (a) of this figure clearly shows that the median radius depends almost linearly on the initial virial parameter for . This result is consistent with what we expect from the initial collapse associated with the evolution of a cluster with subvirial velocities: The quasi-equilibrium radius (that obtained after the initial dynamical adjustment) scales linearly with the initial virial parameter. Clusters with completely subvirial starting states have median radii that are approximately of the median radii of virial clusters.

The median radius also decreases as a function of the gas removal time . The data points connected by a solid line in Figure 15 panel (b) correspond to the time averaged (over Myr) normalized median radius of clusters with differing values of . During the embedded phase, these (initially subvirial) clusters remain bound and do not expand. As a result, clusters that become unbound early in their history have larger median radii simply due to time averaging (over the Myr time interval of the simulations). The data points connected by the dashed curve correspond to the normalized median radius averaged over the embedded stage of the cluster evolution ( Myr). Removing the apparent dependence on that is actually due to the time averaging, we find that in clusters with dispersal times greater than Myr, the cluster median radius does not depend sensitively on the gas removal time .

Note that clusters with early gas dispersal times ( Myr) have significantly larger median radii than clusters with later gas dispersal times. The average crossing time in a subvirial cluster is Myr, and thus gas removal within the first couple crossing times prevents the cluster from approaching a state of virial equilibrium. This behavior also explains why the bound fractions in clusters with Myr are very low (see Figure 5 panel (a), and the discussion in Section 3.1). In other words, the process of gas dispersal in a cluster that is not in virial equilibrium is more destructive to the cluster than if gas removal occurs after the system approaches an equilibrium state (Goodwin & Bastian, 2006).

Panel (a) in Figure 16 displays the normalized median cluster radius as a function of the initial scaling radius used in equation (1). The normalization of the median cluster radius includes the intrinsic dependence on , and hence the trend observed in the normalized median radius must be accounted for by another mechanism. The larger normalized median radius observed in clusters with smaller initial values of can be understood in terms of the higher interaction rates observed in these clusters; the interactions keep the cluster cores slightly inflated. This trend should be present, to some extent, in all clusters with high interaction rates; however, it is easiest to observe in the series of simulations because the interaction rates have the widest dynamical range, varying by three orders of magnitude (see Figure 6, panel (c)).

Tables 14 – 19 list the values of the parameters that specify the mass and number profiles for the clusters considered herein, where the profiles have the form given by equations (11) and (12). Each table lists the fitted parameters as a function of a given input variable, including the stellar membership (Table 14), initial virial parameter (Table 15), cluster scaling radius (Table 16), star formation efficiency (Table 17), gas removal time (Table 18), and the mass segregation parameter (Table 19). These tabulated values, in conjunction with equations (11) and (12), provide analytic descriptions of the mass and number density distributions for a wide variety of clusters. These analytic forms, in turn, can be used to calculate related physical quantities. For example, the cluster median radius is given by the formula

 rmed=r0(21/p−1)−1/a, (13)

where are the parameters that specify the number density profile. As another example, the mass density of the cluster is given by

 ρ(r)=14πr2dMdr=MT∗4πr3apξap(1+ξa)p+1=M(r)4πr3ap1+ξa, (14)

where the parameters are those of the mass profile. Similarly, the magnitude of the gravitational force is given by , and the corresponding potential of the stellar component can be obtained from the integral . Note that the force can be expressed in terms of elementary functions, but the potential integral leads to hypergeometric functions for general values of the indices. These analytic forms can be used in a wide variety of applications to help determine the effects of cluster environments on forming stars and young solar systems.

4 Effects of Cluster Environment on Planetary Formation

In this section, we consider two mechanisms through which young embedded clusters affect their constituent members. Note that for systems in the regime of parameter space considered here, the background cluster environment has more effect on circumstellar disks and planet formation than on the star formation process itself. In Section 4.1, we use the interaction rates calculated from our ensemble of -body simulations (Section 3.2) to determine the effectiveness of close encounters for disrupting circumstellar disks and newly formed planetary systems. In Section 4.2, we use the output number density profiles (from Section 3.4) to determine how much UV radiation the cluster provides to the circumstellar disks surrounding its members.

4.1 Close Encounters

One way that planet formation may be compromised in stellar clusters is through close encounters between planet-forming circumstellar disks and other cluster members. A close encounter will truncate a circumstellar disk down to a radius that is roughly the distance of closest approach (Kobayashi & Ida, 2001). We thus need a metric to determine the importance of close encounters during the 10 Myr embedded phase of cluster evolution (note that this timescale is comparable to the observed lifetimes of circumstellar disks and the expected time required for giant planet formation). Using data from our cluster simulations, we define the characteristic impact parameter as the closest encounter that an average star will experience over the course of a given time interval, taken here to be . The value of is readily calculated from the distribution of interaction rates. For the form given in equation (6), and we define the characteristic radius to be

 bC≡1000AU(t0Γ0)1/γ, (15)

where = 10 Myr.

Figures 17 and 18 display the characteristic impact parameter as a function of initial cluster parameter for the clusters considered in this parameter space survey. The horizontal dashed and solid lines (at  AU and  AU) indicate encounter distances that restrict planet formation in a circumstellar disk to distances less than  AU and  AU, respectively (Kobayashi & Ida, 2001). These latter distances correspond to the outer edge of the Kuiper Belt and the orbit of the planet Neptune in our solar system.

These figures indicate that most of the clusters in this parameter space survey do not have interaction rates that are high enough to seriously compromise planet formation within the  AU Kuiper Belt radius, or even within  AU. Equivalently, disruptive encounters at these close distances are predicted to be rare. In the most interactive clusters, those with the highest number densities due to small values of (Figure 17, panel (c)), and those with significant amounts of mass segregation (Figure 18, panel (b)), the characteristic scale is small enough that planet-forming disks may be truncated at radii  AU and hence giant planet formation in this class of clusters can be partially inhibited. Nonetheless, the cluster environment provides only moderate constraints on the planet forming process.

Interactions remain important even after planet formation in the circumstellar disk has taken place. Close encounters with other cluster members are capable of disrupting newly formed planetary systems. In this context, the primary channel of “disruption” is to increase the orbital eccentricity and/or the inclination angles of the planetary orbits. Sufficiently close encounters can eject planets from their orbits entirely, and produce free floating planets in the clusters. The stellar interaction rates presented in this paper may be combined with scattering calculations to investigate the rate of solar system disruption in young clusters.

One way to characterize the possible effects of clusters on newly formed planetary systems is to pose the following question: Under the assumption that planet formation naturally produces systems similar to our solar system (with planets in relatively circular orbits with small inclination angles), how many planetary systems in a young cluster will be noticeably different due to an encounter with another cluster member? Previous work provides cross sections for disrupting solar systems in a variety of ways (e.g., Adams & Laughlin 2001, APFM); as expected, scattering interactions are most effective at altering the orbits of the outermost planets. For the sake of definiteness, we consider a collection of solar systems that have the same architecture as our own (same masses and semimajor axes for the outer planets) except that the orbits are circular. We then designate a solar system to be “noticeably different” (due to an encounter) when the analog of Neptune has its orbital eccentricity increased (from zero) to values (about twice the value observed in our solar system). We note that nothing especially dramatic happens when orbital eccentricities are increased by this amount, only that the solar system is changed enough to be noticed. Previous work shows that the cross section for increasing Neptune’s eccentricity to during an encounter with a cluster member is 167,000 AU (Adams & Laughlin, 2001), which corresponds to a closest approach distance of approximately 230 AU. This cross section for increasing Neptune’s eccentricity is about the same as the cross section for increasing the spread in inclination angles beyond that observed in our solar system (). For completeness, we note that in these solar system scattering calculations, the perturber is assumed to be a binary star system, which is reasonable since a large fraction of stars are members of multiple systems.

As shown in Figures 17 and 18, the typical interaction distance is somewhat larger than the  AU encounter distance required to “disrupt” the solar system. However, some fraction of stars in a cluster will experience encounters that are significantly closer. For instance, consider a cluster with = 300 stars, subvirial initial conditions characterized by = 0.04, and an initial cluster radius pc. In this cluster, (see Figure 17, panel (c)). The rate of interactions with an encounter distance less than or equal to can be determined from equation (6) using the values = 0.333 and . The rate of close encounters with  AU is thus about interactions per star per Myr. As a result, in this cluster approximately 115 stars will experience close encounters severe enough to change their planetary systems. In contrast, in a cluster similar in all respects except with a slightly larger initial radius (i.e., , , and pc), this number drops to planetary systems ().

This example illustrates the type of calculations that are possible using the interaction rates and fitting parameters presented in this numerical study. The close encounter distributions may also be combined with results presented by other scattering calculations to determine the efficacy of circumstellar disk and planetary system disruption in young stellar clusters (Heller, 1993, 1995; de La Fuente Marcos & de La Fuente Marcos, 1997, 1999; Kobayashi & Ida, 2001; Adams & Laughlin, 2001; Adams et al., 2006; Pfalzner, 2008).

Before leaving this section, we can use the interaction rates determined here to provide a consistency check. Although this ensemble of simulations does not include binary systems (only single stars), the presence of binaries could affect the energy budget of the system. As outlined above, relatively few stellar interactions take place with distances of closest approach less than 100 AU, and the characteristic radii are typically much larger ( AU). For comparison, for solar type stars, the peak of the binary period distribution occurs at a period days (Duquennoy & Mayor 1991), which corresponds to a distance of AU. As a result, the closest encounter experienced by the the vast majority of stars (during the 10 Myr window of interest) is much wider than the typical binary separation. These results thus vindicate the approach of ignoring binarity for this class of simulations. However, we note that for more extreme regimes of parameters space (e.g., dense stellar systems destined to become globular clusters), neglecting binaries is not a good approximation.

Another mechanism through which cluster environments may affect planet formation is by photoevaporation of protoplanetary disks due to the enhanced FUV radiation fields produced by massive young stars. It is well known that radiation from the central host star can heat and dissipate its surrounding planet-forming disk (Shu et al., 1993; Hollenbach et al., 1994). Most stars are not massive enough to produce large quantities of FUV radiation, the wavelength range that most effectively photoevaporates the disk. In young clusters, however, the most massive cluster stars can provide FUV radiation fields that are strong enough to photoevaporate the disks associated with other cluster members (Johnstone et al., 1998; Adams & Myers, 2001) and dominate over the radiation produced by the host star.

Recent studies have determined the FUV (and EUV) luminosities of clusters as a function of cluster membership , the mass function of the cluster, and various amounts of extinction within the cluster (e.g., Armitage 2000, APFM, Fatuzzo & Adams 2008). In this section, we combine these previous determinations of the typical FUV background luminosities in stellar clusters with the radial profiles presented in Section 3.4. Taken together, these stellar positions and the FUV luminosities determine the expected FUV flux that impinges upon the circumstellar disks in young stellar clusters. In these systems, this FUV flux places limits on both the timescale over which planets may form and the region of the disk that has the potential to form planets.

The total FUV luminosity of a cluster originates primarily from the most massive stars in the cluster. As a result, the cluster FUV luminosity is sensitive to the membership size of the cluster and the mass function of the stars in the cluster. APFM presented detailed calculations of the expected FUV luminosity in a cluster of size (see their Figure 6). For each cluster membership size , the distribution of possible FUV luminosities is quite wide due to under-sampling of the stellar initial mass function; for clusters with , the width of the FUV luminosity distribution is larger than the mean (or median) value. As a result, clusters are predicted to display substantial system to system variation in the radiation fields they produce (especially for 700).

Each individual cluster member will experience a FUV flux that is time-dependent as it orbits through the cluster potential. For purposes of this paper, we estimate the typical flux experienced by a circumstellar disk in a simulated cluster by combining the FUV luminosity with the median radial position of the stars in a simulated cluster. The vast majority of the FUV radiation in a cluster is produced by the most massive stars, which typically reside near the cluster’s center (see Section 2.2.6). The FUV source is thus modeled as a point source at the center of the cluster. For a given cluster with FUV luminosity and median radius , the corresponding median FUV flux experienced by the protoplanetary disks in the cluster is defined to be

 FFUV≡LFUV(N)4πr2med, (16)

where is defined by the radial probability distribution of the cluster (equation [12]) and is taken from the results presented in APFM.

In addition to the FUV flux provided by the background cluster, a circumstellar disk is also subjected to FUV radiation from its host star. Observations of T Tauri stars provide estimates of the FUV flux experienced by their disks at distances of  AU due to the central star (Bergin et al., 2004). For three T Tauri stars, these authors determine FUV fluxes of = 240, 340, and 1500 where corresponds to a benchmark value of  erg s cm (close to the value of the interstellar radiation field at FUV wavelengths). As a rough estimate, the FUV flux associated with the host star can be taken to be . When the background FUV radiation field of the cluster exceeds this benchmark value, the environment can, in principle, affect the evolution of circumstellar disks and planet formation. However, in the outer regions of a cluster, where flux levels from the central massive stars are relatively low, the FUV flux from the host star can be more important than the FUV flux from the background cluster.

Figures 19 and 20 present the median cluster FUV flux (defined by equation [16]) as a function of the initial cluster parameters for the clusters considered in this parameter space survey. The dashed horizontal lines correspond to FUV radiation levels of = 300 and 3000, values for which the effects of photoevaporation on circumstellar disks have been calculated in detail (Adams et al., 2004). The results of these photoevaporation models provide a rough scaling law: Over the course of  Myr, an FUV flux of is capable of truncating a circumstellar disk down to a radius given by

 rd≈36 AU(M∗M⊙), (17)

where the result depends on the stellar mass . This nominal radius is close to the size of our solar system. As a rough rule of thumb, significant evaporation around solar type stars thus requires flux levels , with a weaker requirement for smaller stars.

The FUV luminosity of a cluster is an increasing function of the cluster membership , and this trend is visible in the FUV fluxes shown in Figures 19 and 20. In subvirial clusters with more than members, circumstellar disks are exposed to median FUV fluxes that are of order (see Figure 19, panels (c) and (d)). In this class of clusters, significant photoevaporation of the outer disk is predicted to take place. Many of the circumstellar disks in such cluster systems will be truncated, so that planet formation is limited to the inner  AU for solar type stars. For the more numerous red dwarf stars, planet formation is limited to the inner AU. In clusters with smaller memberships , the radiation fields are more modest and disks will not be significantly photoevaporated by cluster radiation. We note that in this latter class of clusters, the FUV radiation from the host star is comparable to that of the background cluster. Finally, we note that in clusters where gas expulsion takes place on short time scales ( 3 Myr), the FUV fluxes are lower due to early expansion of the cluster and larger values of (see panel (b) of Figure 19).

5 Conclusion

In this paper, we have presented the results from a large ensemble of numerical simulations designed to study the dynamics young embedded clusters over a wide range of conditions. The choice of parameter space was motivated by recent catalogs and surveys of star-forming environments in the solar neighborhood. We consider clusters with a range of stellar memberships , a variety of parameter values and that define the cluster membership-radius relation (equation [1]), a range of gas removal timescales , initial virial states , star formation efficiencies , and amounts of primordial mass segregation . The range of parameter space surveyed is summarized in Table 1.

The results of this survey show how the evolutionary parameters that describe the properties of evolving young clusters vary as a function of the initial cluster conditions. Section 3.1 considers how the cluster’s bound fraction (generally evaluated at time ) varies with the input parameters. The bound fraction depends most sensitively on the star formation efficiency , but also depends quite strongly on the initial virial state of the system. As a general rule, clusters in which the stars are formed with subvirial velocities have higher bound fractions than clusters in which stars are formed with initially virialized velocities.

In Sections 3.2 and 3.3, we considered the distributions of close encounters between cluster members over the first 10 Myr of cluster evolution. We find that the encounter rates scale linearly with the average stellar density, so that as expected, and that the interaction rates increase as the initial virial parameter decreases. Subvirial clusters have interaction rates that are roughly times higher than those found in clusters that start in a virialized state. We also show that the interaction rates are higher in clusters that have larger amounts of primordial mass segregation. In addition, the distributions of interaction velocities are nearly Gaussian (see Figure 10). The interaction velocities themselves do not vary strongly as a function the initial conditions in the cluster; instead, the interaction velocities are always about twice the average stellar velocity in the cluster.

In Section 3.4, we present the results of empirical fits to the cluster radial profiles and . We find that, in general, the median cluster radius scales as the initial cluster radius , so that the ratio is nearly invariant. In some clusters, however, significant interactions between stellar members produce a break in the relationship. More specifically, clusters that are initially subvirial develop smaller median radii than those that are initially in virial equilibrium. In addition, if the embedding gas is removed from a subvirial cluster early in the cluster’s evolution (), the resulting median cluster radius is much larger than it would be if gas removal occurred at a later time. This difference occurs because the cluster is not in (or near) virial equilibrium at the time of gas removal; as a result, gas dispersal is more destructive in a cluster that has not yet approached virial equilibrium.

We note that these general trends in the output parameters found in this study (e.g., the interaction rates are proportional to the mean density ) are not unexpected. However, this work puts these results on a firm, statistically significant footing. We find quantitative results, with accuracy and precision well beyond that obtained from qualitative scaling arguments. In addition, we obtain these results for varying values of the input parameters that characterize the cluster.

In summary, the principle contribution of this work is to provide a collection of output parameters that describe cluster properties and cluster evolution as a function of initial conditions (see Tables 2 – 19). The interaction rates for close encounters between cluster members are given by the power-law form of equation (6), where the parameters are listed in Tables 8 – 13. Taken together, these results provide us with an analytic description of the interaction rates for a wide range of cluster input parameters. Similarly, the radial profiles of mass and stellar number take the forms given by equations (11) and (12), where the parameters are listed in Tables 14 – 19. We thus also have an analytic description of cluster mass and number distributions. From these functions, additional physical quantities can be derived, including mass density, number density, stellar potential, median radius, and many others (see equations [13] and [14]). Furthermore, all of these analytic quantities are specified as a function of stellar membership , initial virial parameter , scale radius , star formation efficiency , gas removal time , and degree of mass segregation . These results can be combined with calculations of the radiation fields in young clusters and cross sections for planetary system disruption to determine the effects of a wide variety of cluster environments on planetary formation. In Section 4, we have provided a few examples of the types of information that can be readily extracted from this data set. These results — along with analogous follow-up studies — will provide a more complete statistical description of cluster evolution and can be used to help understand the impact of the cluster environment on planet formation.

We thank Lori Allen, Rupali Chandar, Rob Gutermuth, Tom Megeath, and Phil Myers for useful discussions. This work was supported in part by the Michigan Center for Theoretical Physics. EMP acknowledges support of the University of Michigan Physics Department through a dissertation fellowship. FCA is supported by the NASA Origins of Solar Systems Program via grant NNX07AP17G, and by the NSF Division of Applied Mathematics via grant DMS-0806756.