DYNAMICAL EVOLUTION OF YOUNG EMBEDDED CLUSTERS: A PARAMETER SPACE SURVEY

Cluster membership, N. As outlined above, in this study we consider intermediate-sized clusters with stellar memberships ranging from N = 100 to 3000 (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 IMF observed in young stellar clusters is almost universal for clusters with more than N ∼ 100 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 0.5 M_☉ (the average stellar mass is somewhat higher than the median stellar mass, which is roughly ∼0.3 M_☉), consistent with observations of young stellar clusters (Muench et al. 2002; Luhman et al. 2003).

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

$$\begin{equation} R_{c}= R_{\rm sc}\left(\frac{N}{300}\right)^{\alpha }. \end{equation} \tag{ 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 R_sc = 1 pc and α = 1/2. Thus, a cluster with N = 300 stars typically has a radius of 1 pc.

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 ρ_* ∼ r⁻¹. 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 the 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 Q_i. The virial ratio is defined as Q ≡ |K/W|, where K is the total kinetic energy and W is the total potential energy of the cluster. A cluster that is in virial equilibrium has a virial ratio Q = 0.5. Most simulations considered in this study are initialized with a virial ratio Q_i = 0.04, 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 Δt = 1 Myr. The formation time of each star is sampled from a uniform distribution over the range of 0 to 1 Myr, 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 ∼3 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

$$\begin{equation} \Psi = \frac{2 \pi G \rho _0 r_s^2}{1 + \xi }, \quad \rho = \frac{\rho _0}{\xi (1 + \xi)^3}, \quad {\rm and} \quad M = \frac{M_\infty \xi ^2}{(1+\xi)^2}, \end{equation} \tag{ 2 }$$

where ξ ≡ r/r_s. Here the parameter r_s is a scale length, which is chosen to be equal to the cluster radius, i.e., r_s = R_c. In the inner limit, this profile has the form of ρ ∼ r⁻¹, 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 ε_SF in young star-forming regions vary from ∼10% to 50% (Lada & Lada 2003). In our simulated clusters, a standard star formation efficiency of ε_SF = 33% is assumed. This value corresponds to a total stellar mass M_T* in the cluster that is one half of the mass of the embedding gas M_gas (from Equation (2), M_gas = M_∞/4 is the effective gas mass within the cluster radius R_c).

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 ∼5 Myr 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 t_gas. 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 M_* = 20 M_☉ 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 M_*>20 M_☉ will finish their evolution during the 10 Myr window considered here, and only about one 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).

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), 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 parameter(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 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 ε_SF = 0.33. The stellar membership in the simulated clusters ranges from N = 100 to 3000. 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 Q_i = 0.04 and 0.5, respectively. The subvirial clusters thus have initial stellar velocities that are approximately 1/3 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 α = 1/2 is a reasonable fit to the observed data within approximately 2 kpc of the Sun (Lada & Lada 2003; Porras et al. 2003), where this sample contains intermediate-sized clusters with N ≲ 2000. 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 of which are thought to be progenitors of globular clusters, contain as many as N ∼ 10⁶ stars and have sizes on the order of R_c ∼ 10 pc (Mengel et al. 2008). If we extend the cluster membership–radius relation out to stellar memberships as high as N ∼ 10⁶, the choice of α = 1/2 would overestimate the cluster radius by a factor of ∼5. A power-law index of α = 1/4 more closely reproduces the observed data points over the full range of N. In this study, we investigate the evolution of intermediate-sized clusters using both α = 1/2 and α = 1/4 power-law indices in the cluster membership–radius relation. In both cases, we chose R_sc = 1.0 pc so that the power-law passes through the point where N = 300 and R_c = 1.0 pc.

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 n₀ gives the relation

$$\begin{equation} n_0 \sim \frac{N}{R_{c}^3} \sim \frac{N}{N^{3\alpha }} = N^{1-3\alpha }. \end{equation} \tag{ 3 }$$

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

Note that an intermediate value of the index, α = 1/3, implies a constant stellar density. This benchmark density value is n₀ ∼ 100 stars pc⁻³. 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, n_M = 65 pc⁻³.

Zoom In Zoom Out Reset image size

2.2.2. Initial Virial Parameter

2.2.3. Cluster Scaling Radius

2.2.4. Star Formation Efficiency

2.2.5. Gas Removal Timescale

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 approximately given by

$$\begin{equation} t_{\rm relax} \approx \frac{N}{8 {\rm ln} N}\tau _c\approx \frac{N R_{c}}{8 {\rm ln} N \left<v\right>}, \end{equation} \tag{ 4 }$$

where τ_c is the cluster crossing time and <v> 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 N ∼ 300, R_c ∼ 1 pc, <v> ∼ 1 km s⁻¹, and a corresponding relaxation time of roughly t_relax = 6.5 Myr. 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 F_seg as the fraction of the cluster membership, which has been ordered by mass at the center of the cluster, F_seg = N_seg/N. More specifically, at the start of the simulation, the stellar masses are sampled from a standard IMF and stellar positions are sampled from a ρ_* ∼ r⁻¹ 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 F_seg>1/N, additional mass segregation is implemented by rearranging the stellar positions so that the N_seg most massive stars are located at the inner N_seg radial positions. The mass segregation parameter is varied over the range 1/N ⩽ F_seg ⩽ 0.99 in subvirial clusters with N = 300, 1000, and 2000 members.

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 10%) of the stellar population is born within "robust" clusters that are destined to become open clusters (which then live for ∼100 Myr 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, 2003b). 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 f_b of stars that remain gravitationally bound as a function of time. The bound fraction is defined as f_b ≡ N_bound/N, where N is the initial stellar membership, and N_bound 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 t_gas. 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 ∼10⁴ 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 f_b decreases significantly (by as much as 50%) over a short period of time. The fraction f_b then continues to decrease (more slowly) until the end of the simulations. Note that our temporal cutoff is chosen to be 10 Myr, but the clusters will continue to evolve and f_b will continue to decrease on longer timescales.

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

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The value of the bound fraction at t = 10 Myr 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 t = 10 Myr 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 ε_SF 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 f_b at t = 10 Myr depends sensitively on the star formation efficiency ε_SF of the cluster. Panel (a) in Figure 4 displays the cluster bound fraction f_b as a function of star formation efficiency, ε_SF. The data are well fitted by a power law in ε_SF:

$$\begin{equation} f_{b}= 2.23 \left(\varepsilon _{\rm SF}\right)^{1.2}, \qquad {\rm where} \qquad \varepsilon _{\rm SF}\le 0.5 \,. \end{equation} \tag{ 5 }$$

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

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the suite of simulations used to investigate the effects of SFE, the clusters are initially subvirial. After gas removal, subvirial clusters are more tightly bound than virial clusters (APFM). For even relatively high SFEs (ε_SF = 0.3) and small initial virial parameters, conditions that produce the most tightly bound systems, the clusters are significantly disrupted by gas removal and promptly lose ∼40% of their stars. SFEs larger than ∼30% 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 2003a, 2003b).

As mentioned above, clusters with subvirial initial velocities are more tightly bound than clusters with virial initial conditions. Panel (b) of Figure 4 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 Q_i over the range considered, 0.01 ⩽ Q_i ⩽ 0.5. 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 Q_i = 0.025, approximately 40% 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 N, for both the R_c ∼ N^1/2 and the R_c ∼ N^1/4 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 (R_c ∼ N^1/2). The bound fraction at t = 10 Myr is plotted as a function of cluster size N in panels (c) and (d) of Figure 4. The upper curves in these panels indicate the bound fraction at 10 Myr 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 (R_c ∼ N^1/4), the bound fraction decreases as a function of the cluster membership N. The bound fraction decreases roughly as f_b ∼ N^−1/4 (see Figure 4, panel (c), lower curve). This decrease in f_b(N) is due to a combination of effects arising from the relationship between cluster radius and cluster membership defined by Equation (1). In clusters with R_c ∼ N^1/4, the mean velocity and velocity dispersion scale approximately as $\left<v\right> \sim \sigma _v\sim \sqrt{GmN/R_{c}} \sim N^{3/8}$. 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 N results in more stars with velocities high enough to escape from the cluster, $v_{\rm esc} \sim \sqrt{2GmN/R_{c}}$. In addition, the interaction rate between cluster members increases with the stellar density, which increases as n ∼ N^1/4 in these clusters. In virial clusters with R_c ∼ N^1/2 (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 N is not as strong: <v> ∼ σ_v ∼ N^1/4. In addition, the stellar density actually decreases with N, n ∼ N^−1/2; 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 t_gas or the scaling radius R_sc, i.e., the bound fraction f_b is not a monotonic function (see Figure 5, panels (a) and (b)). This result occurs because in subvirial clusters, such as those 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 f_b 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 20% 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 F_seg = 0.2, as shown in panel (c) of Figure 5.

In summary, the results of this part of the study indicate that the star formation efficiency ε_SF 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 and 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, Q_i ≲ 0.2, the bound fraction is not a sensitive function of the initial stellar density (as indicated by the suite of simulations varying N 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 ∼10–20 Myr, 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 ∼10 Myr will lead to even lower bound fractions f_b 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 f_b evaluated at time t = 10 Myr for the simulations presented here. Each table lists the bound fractions f_b as a function of a given input variable, including the stellar membership N (Table 2), initial virial parameter Q_i (Table 3), cluster scaling radius R_sc (Table 4), star formation efficiency ε_SF (Table 5), gas removal time t_gas (Table 6), and the mass segregation parameter F_seg (Table 7).

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 the 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 b = 10⁴ 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 t = 0 − t_gas, the exposed phase t = t_gas − 10 Myr, or the entire 10 Myr interval). Specifically, the interaction rate Γ(b) is defined as the number of close encounters with distance of closest approach r ⩽ b 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 ∼3000 AU. In other words, the interaction rate has the form

$$\begin{equation} \Gamma = \Gamma _0\left(\frac{b}{1000\, \rm {AU}}\right)^{\gamma } \,, \end{equation} \tag{ 6 }$$

where b is the distance of closest approach. The fiducial interaction rate Γ₀ and the power-law index γ are fitted 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 b less than 1000 AU per star per million years. The fiducial interaction rate Γ₀ is displayed as a function of cluster initial conditions in Figures 6 and 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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

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

$$\begin{equation} \delta n = \frac{2N}{R_{c}^2}b \delta b \,. \end{equation} \tag{ 7 }$$

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

$$\begin{equation} \delta \Gamma \approx \frac{2N}{R_{c}^2}\left(\frac{v}{R_{c}}\right)b \delta b = 2 n v b \delta b \,. \end{equation} \tag{ 8 }$$

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

In panel (d) of Figure 6, the Γ₀ values are plotted as a function of initial stellar density n₀ = N/R³_c for all of the simulations with varying cluster membership N and scaling radius R_c (using the data from panels (a)–(c) of Figure 6). 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 ×'s. Straight lines (indicating a power-law with index = 1) are included in the panel; note that the numerically determined data are roughly consistent with the simple scaling relation Γ₀ ∼ n₀.

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 ∼1–2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray} \Gamma _b &\approx& 0.122 \, \, N^{-1/2} \left(\frac{b}{1000\ {\rm AU}}\right)^2\;{\rm Myr}^{-1} \nonumber\\ &&\qquad {\rm for} \qquad R_{c}= 1.0\;{\rm pc} \left(\frac{N}{300}\right)^{1/2}, \end{eqnarray} \tag{ 9 }$$

and

$$\begin{eqnarray} \Gamma _b &\approx& 0.0095 \, \, N^{1/4} \left(\frac{b}{1000\ {\rm AU}}\right)^2\;{\rm Myr}^{-1} \nonumber\\ &&\qquad {\rm for} \qquad R_{c}= 1.0\;{\rm pc} \left(\frac{N}{300}\right)^{1/4}. \end{eqnarray} \tag{ 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 2 (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 eight times larger than the rates for virial clusters of the same starting density (as defined by n₀ ∝ N/R³_c). 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 (Q_i = 0.04) behave as if they have nearly zero-temperature starting states and thus collapse to roughly $\sqrt{2}$ of their initial size. This decrease in radius corresponds to an increase in density by a factor of $2\sqrt{2}$. In addition, subvirial clusters retain more of their members after the gas is removed from the cluster—we find that the bound fraction is ∼ three 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 three 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 ${\sim} 6\sqrt{2} \sim 8$. The results of the parameter survey varying Q_i also indicate that subvirial clusters have higher interaction rates. In Figure 7(a), the interaction rate clearly decreases as a function of initial virial parameter Q_i.

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(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 five 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 N = 300 stars and F_seg = 0.05 have fiducial interaction rates that are about five times larger than those found in subvirial clusters with N = 300 members and minimal mass segregation (see Figure 7(c), top curve).

The average interaction rate also increases as a function of gas removal time t_gas as shown in Figure 7(b). This interaction rate is averaged over the 10 Myr 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 ∼2 Myr, the rate of close encounters during the embedded phase is roughly constant. When the embedded phase lasts less than ∼2 Myr, the clusters have lower encounter rates due to lower densities during the first ∼1 Myr 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 ε_SF (see Figure 7(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 N (Table 8), initial virial parameter Q_i (Table 9), cluster scaling radius R_sc (Table 10), star formation efficiency ε_SF (Table 11), gas removal time t_gas (Table 12), and the mass segregation parameter F_seg (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).

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 v_enc 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 R_c = 1 pc, N = 300 stars, and subvirial initial speeds (Q_i = 0.04). A binning size of 0.25 km s⁻¹ has been used to construct the histogram, and error bars are included on the distribution to indicate the dispersion within each velocity bin.

Zoom In Zoom Out Reset image size

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 v_enc ≲ 0.5 km s⁻¹ (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 <v_enc> 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. These figures demonstrate 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 $\ssty\sqrt{v_{\rm rel}^2} \sim \sqrt{2\left<v^2\right>}$ (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 b ∼ 10³ AU is a typical encounter distance, and that the orbit speed v_orb ∼ (GM_*/b)^1/2 ∼ 1 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 b ≲ 10⁴ AU are included in the interaction velocity distribution, and this subsample is likely to have somewhat larger relative velocities.

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 t can be characterized by the cumulative mass distribution M(r, t)/M_T* or the cumulative number distribution N(r, t)/N, where M_T* and N are the total masses and numbers of the stars that are gravitationally bound to the cluster at time t, respectively. In the simulated clusters, each of these distributions is calculated at intervals of 0.25 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 fitted by simple functions of the form

$$\begin{equation} \frac{M(r)}{M_{T\ast }} = \left(\frac{\xi ^a}{1+\xi ^a}\right)^{p}, \end{equation} \tag{ 11 }$$

$$\begin{equation} \frac{N(r)}{N} = \left(\frac{\xi ^a}{1+\xi ^a}\right)^{p}, \end{equation} \tag{ 12 }$$

where ξ = r/r₀ and the scale radius r₀ and index p are free parameters that are fitted to the distributions observed in the simulated clusters. The parameter a may also be varied to fit the data. We find that the choice a = 2 gives the best fit for the subvirial clusters, and a = 3 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 Q_i is varied, the choice a = 2 works best for simulations for which Q_i < 0.25, and a = 3 works best for those with Q_i ⩾ 0.25. The parameters r₀ and p 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.

These radial profiles of clusters provide insight into the general evolution of a cluster, and perhaps more importantly, into the expected radiation fields that young solar systems in the cluster will experience. Circumstellar disks and forming solar systems residing in a cluster will be subjected to the FUV and EUV radiation fields produced by the cluster population, and these fields are dominated by the large stars in the cluster. If these radiation fields are strong enough, they are capable of photoevaporating circumstellar disks and thereby preventing (or at least limiting) giant planet formation. The massive young stars that produce the majority of the UV radiation are often located near the center of the cluster (Hillenbrand & Hartmann 1998; Bonnell & Davies 1998; Lada & Lada 2003). To first approximation, the UV radiation in a cluster can be considered as a point source at the center. An understanding of the EUV and FUV fields associated with young clusters combined with the average radial distributions of stars in young clusters thus provides a framework with which to predict how effectively cluster radiation can restrict planet formation (Johnstone et al. 1998; Fatuzzo & Adams 2008).

Although one can use the full distributions (see Equations (11) and (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 r_med calculated from the fits to the cumulative radial distributions for the entire parameter space. The scale r_med is defined as the radius at which N(r)/N = 0.5, and thus represents the radius which, on average, contains half of the cluster members. As these plots indicate, the median cluster radius r_med scales with the initial virial parameter Q_i, the gas removal timescale t_gas, the cluster membership N, and the cluster scaling radius R_sc. On the other hand, the radius r_med does not vary strongly with the either star formation efficiency ε_SF or the degree of primordial mass segregation F_seg.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Scaling the median radius r_med by the initial cluster radius R_c removes the dependence on this initial cluster parameter and more readily identifies trends that are distinct from the initial assumptions concerning the cluster size. Figure 15 displays the median radius r_med normalized by the initial cluster radius R_c. Panel (a) of this figure clearly shows that the median radius depends almost linearly on the initial virial parameter Q_i for Q_i ⩽ 0.5. 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 $\sqrt{2}$ of the median radii of virial clusters.

Zoom In Zoom Out Reset image size

The median radius also decreases as a function of the gas removal time t_gas. The data points connected by a solid line in panel (b) of Figure 15 correspond to the time averaged (over 0–10 Myr) normalized median radius r_med of clusters with differing values of t_gas. 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 10 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 (0 − t_gas Myr). Removing the apparent dependence on t_gas that is actually due to the time averaging, we find that in clusters with dispersal times greater than ∼2 Myr, the cluster median radius does not depend sensitively on the gas removal time t_gas.

Note that clusters with early gas dispersal times (t_gas ≲ 2 Myr) have significantly larger median radii than clusters with later gas dispersal times. The average crossing time in a subvirial cluster is ∼1 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 t_gas ≲ 2 Myr are very low (see panel (a) of Figure 5, 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 R_sc used in Equation (1). The normalization of the median cluster radius includes the intrinsic dependence on R_sc, and hence the trend observed in the normalized median radius r_med/R_c must be accounted for by another mechanism. The larger normalized median radius observed in clusters with smaller initial values of R_sc 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 R_sc series of simulations because the interaction rates have the widest dynamical range, varying by three orders of magnitude (see Figure 6, panel (c)).

Zoom In Zoom Out Reset image size

Tables 14–19 list the values of the parameters (p, r₀, a) 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 N (Table 14), initial virial parameter Q_i (Table 15), cluster scaling radius R_sc (Table 16), star formation efficiency ε_SF (Table 17), gas removal time t_gas (Table 18), and the mass segregation parameter F_seg (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 r_med is given by the formula

$$\begin{equation} r_{\rm med} = r_0 \left(2^{1/p} - 1 \right)^{-1/a} \,, \end{equation} \tag{ 13 }$$

where (p, r₀, a) are the parameters that specify the number density profile. As another example, the mass density of the cluster is given by

$$\begin{equation} \rho (r) = \frac{1}{4 \pi r^2} \frac{d M}{dr} = \frac{M_{T\ast }}{4 \pi r^3} \frac{ap \, \xi ^{ap}}{(1 + \xi ^a)^{p+1}} \, = \frac{M (r)}{4 \pi r^3} \frac{ap}{1 + \xi ^a} \,, \end{equation} \tag{ 14 }$$

where the parameters (p, r₀, a) are those of the mass profile. Similarly, the magnitude of the gravitational force is given by |F| = GM(r)/r², and the corresponding potential of the stellar component can be obtained from the integral Ψ_* ∼ ∫Fdr. 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.

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 1/3 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 b_C as the closest encounter that an average star will experience over the course of a given time interval, taken here to be 10 Myr. The value of b_C is readily calculated from the distribution of interaction rates. For the form given in Equation (6), we define the characteristic radius b_C to be

$$\begin{equation} b_C \equiv 1000\ \rm {AU} \left(t_0 \Gamma _0 \right)^{1/\gamma } \,, \end{equation} \tag{ 15 }$$

where t₀ = 10 Myr.

Figures 17 and 18 display the characteristic impact parameter b_C as a function of initial cluster parameter for the clusters considered in this parameter space survey. The horizontal dashed and solid lines (at 190 AU and 120 AU) indicate encounter distances that restrict planet formation in a circumstellar disk to distances less than 55 AU and 30 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 55 AU Kuiper Belt radius, or even within 100 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 R_sc (Figure 17, panel (c)), and those with significant amounts of mass segregation (Figure 18, panel (b)), the characteristic scale b_C is small enough that planet-forming disks may be truncated at radii r_d ≲ 100 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 ⩾ 0.05 (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 ⩾ 0.05 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 (Δi ⩾ 35). 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 b_C is somewhat larger than the 230 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 N = 300 stars, subvirial initial conditions characterized by Q_i = 0.04, and an initial cluster radius R_c = 0.67 pc. In this cluster, b_c ≈ 400 (see Figure 17, panel (c)). The rate of interactions with an encounter distance less than or equal to b can be determined from Equation (6) using the values Γ₀ = 0.333 and γ = 1.47. The rate of close encounters with b ⩽ 230 AU is thus about 0.038 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., N = 300, Q_i = 0.04, and R_c 1 pc), this number drops to ∼50 planetary systems (Γ₀ = 0.1290, γ = 1.39).

This example illustrates the types 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 b_C are typically much larger (∼1000 AU). For comparison, for solar type stars, the peak of the binary period distribution occurs at a period P ≈ 10⁵ days (Duquennoy & Mayor 1991), which corresponds to a distance of ∼40 AU. As a result, the closest encounter experienced by 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 N, 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 N 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 N (see their Figure 6). For each cluster membership size N, the distribution of possible FUV luminosities is quite wide due to the undersampling of the stellar IMF; for clusters with N ≲ 700, 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 N≲ 700).

Each individual cluster member will experience an 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 L_FUV(N) 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 L_FUV(N) and median radius r_med, the corresponding median FUV flux ${\cal F}_{\rm FUV}$ experienced by the protoplanetary disks in the cluster is defined to be

$$\begin{equation} {\cal F}_{\rm FUV}\equiv \frac{L_{\rm FUV}(N)}{4\pi r_{\rm med}^2} \,, \end{equation} \tag{ 16 }$$

where r_med is defined by the radial probability distribution of the cluster (Equation (12)) and L_FUV(N) 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 ∼100 AU due to the central star (Bergin et al. 2004). For three T Tauri stars, these authors determine FUV fluxes of G₀ = 240, 340, and 1500, where G₀ = 1 corresponds to a benchmark value of 1.6 × 10⁻³ 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 G_host ≈ 500. 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 ${\cal F}_{\rm FUV}$ (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 G₀ = 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 10 Myr, an FUV flux of G₀ = 3000 is capable of truncating a circumstellar disk down to a radius r_d given by

$$\begin{equation} r_d \approx 36\, {\rm AU}\left(\frac{M_\ast }{M_{\odot }}\right) \,, \end{equation} \tag{ 17 }$$

where the result depends on the stellar mass M_*. 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 G₀ ≳ 3000, with a weaker requirement for smaller stars.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The FUV luminosity of a cluster is an increasing function of the cluster membership N, and this trend is visible in the FUV fluxes shown in Figures 19 and 20. In subvirial clusters with more than N ∼ 1000 members, circumstellar disks are exposed to median FUV fluxes that are of order G₀ ≈ 3000 (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 ∼35 AU for solar type stars. For the more numerous red dwarf stars, planet formation is limited to the inner ∼12 AU. In clusters with smaller memberships N, 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 timescales (t_gas⩽ 3 Myr), the FUV fluxes are lower due to early expansion of the cluster and larger values of r_med (see panel (b) of Figure 19).