text:
# Modes of Multiple Star Formation

## 1 INTRODUCTION

The past two decades have witnessed the development of a working paradigm for the formation of isolated single stars (e.g., see various reviews from Shu, Adams, & Lizano 1987 to Mannings, Boss, & Russell 2000 \[PPIV\]). In practice, however, stars tend to form in groups and clusters, and many authors have suggested that clusters form the majority of stars in our galaxy (e.g., Evans 1999; Elmegreen 1985). If most star formation takes place within sufficiently dense environments, then the current theory of star formation could require substantial modification, or perhaps even a new paradigm. A vital issue that must be addressed is to determine the distributions of sizes and stellar densities for star forming systems and to estimate how neighboring stars and the background environment affect the formation of a given star. We can view star formation environments – groups and clusters – in two different ways:  We can consider these stellar systems as astronomical objects in their own right and study their birth, evolution, and eventual demise.  We can consider the effects of groups and clusters on the star formation process. These two viewpoints are intrinsically coupled in that the evolution of these systems determines, in part, their effects on star formation. In this work, we necessarily consider both of these points of view.

For stars born within regions of sufficiently high stellar density, many possible effects can influence their formation and subsequent development. If stars form in a stellar system that lives for many dynamical times during the formation stage, then processes such as mergers can affect stellar formation (e.g., Murray & Lin 1996; Bonnell et al. 1997, 1998). If newly formed stars continue to live within a dense cluster-like environment, their final characteristics are influenced by additional effects, including binary capture through disk dissipation (Heller 1993; Ostriker 1994), scattering of planets into more eccentric orbits (Laughlin & Adams 1998), and other related dynamical events (see also Price & Podsiadlowski 1995; Allen & Bastien 1995). Sufficiently large clusters are likely to produce massive stars, which can affect star formation through their intense ultraviolet radiation fields (Störzer & Hollenbach 1999; Hester et al. 1999) and supernova blast waves (e.g., Cameron & Truran 1977).

To assess the importance of these processes, we need to identify the types of stellar systems that form stars and determine the effects of these environments on star formation. Stars form within stellar systems containing $`N`$ members, where $`N`$ is known to vary from 1 to $`10^4`$. The isolated mode of star formation refers to the limit $`N=1`$ (or 2) and allows no influence of neighboring stars. The clustered mode of star formation refers to the limit $`N1`$ and allows for strong influences of neighboring stars. Most stars form in systems with $`N>1`$, but we still seek answers to the following questions: What is the typical size $`N`$ and the distribution of system sizes $`P(N)`$? How does a system of size $`N`$ affect its constituent stars? Is the typical star forming system more nearly isolated or clustered?

In this paper, we take a modest step toward these goals. We first estimate the fraction of star formation that takes place in clusters versus smaller stellar systems that we denote as groups. Along the way, we make a clean distinction between systems with large numbers $`N`$ of stars (clusters) and those with smaller $`N`$ (groups). Specifically, we estimate the critical number of stars $`N_{}100`$ that defines the boundary between these two types of star forming systems. Some distinction between groups and clusters has long been known from observations, especially from infrared imaging of embedded regions. For example, in Orion B (L1630), stars are organized into clusters (Lada et al. 1991) with relatively few groups or distributed (isolated) young stellar objects (Li, Evans & Lada 1997); in contrast, the southern part of Orion A (L1641) contains many groups and distributed stars (Strom et al. 1993, 1994; Chen & Tokunaga 1994).

After making a distinction between clusters and smaller groups, we must assess the effects of these different environments on star formation. Stars born within clusters have a chance to experience disruptive close encounters with other stars, whereas stars born in smaller groups have much smaller odds of such interactions. In this context, close encounters are those that lead to substantial destructive or constructive effects on nascent solar systems (such as binary capture or planetary ejection); all stellar groups evolve through scattering encounters that lead to their eventual dispersal. Because of statistical considerations and the form of the stellar initial mass function (IMF), clusters are large enough to have an appreciable chance of containing massive stars; smaller groups, on average, contain only low mass stars. Stars forming within clusters are thus influenced by nearby supernova explosions and intense ultraviolet radiation fields; stars forming within smaller groups are relatively unaffected by massive stars and their destructive effects.

This paper is organized as follows. Using observed surveys of galactic open clusters, we estimate (in §2) the fraction of stars that are born within clusters. In §3, we find the minimum number $`N_{}`$ of stars required for a stellar system to have its dynamical relaxation time longer than its formation time; only those systems above the threshold $`N>N_{}`$ live long enough ($`100`$ Myr) to be considered as clusters. In §4, we discuss how groups and clusters affect the star formation process through stellar scattering events, supernovae, and background radiation fields. We conclude, in §5, with a summary and discussion of our results.

## 2 FRACTION OF STAR FORMATION IN CLUSTERS

In this section, we estimate the fraction of stars that are initially born within open clusters. Many authors have studied the distribution of open clusters in our galaxy and this fraction is relatively well known ($`f_C0.1`$). In this present discussion, we estimate $`f_C`$ using the observational survey of Battinelli & Capuzzo-Dolcetta (1991; hereafter BC91), who selected a collection of 100 open clusters in the solar neighborhood \[see also the analyses of van den Bergh (1981); Elmegreen & Clemens (1985); Pandey & Mahra 1986; Janes et al. (1988)\]. The principal result of the BC91 survey is the cluster formation rate $``$,

$$=0.45\pm 0.04\mathrm{Myr}^1\mathrm{kpc}^2.$$
(1)
The cluster formation rates from all of the observational surveys are roughly consistent with each other (at the factor of two level) and the BC91 result is among the highest.

To proceed further, we specify the time dependence of the star formation rate, which we take to have an exponentially decreasing form, $``$ e-qt, over the age of the galactic disk. The total number of clusters, per kpc2, produced over the age of the disk is thus

$$N_C=_0^{t_{\mathrm{disk}}}\mathrm{e}^{q(t_{\mathrm{disk}}t)}𝑑t=\frac{}{q}\left\{\mathrm{e}^{qt_{\mathrm{disk}}}1\right\},$$
(2)
where $`t_{\mathrm{disk}}910`$ Gyr (Wood 1992) and where $``$ is the present day cluster formation rate (given by eq. ). The value of $`q`$ can be estimated from models of the chemical evolution of the galaxy (Rana 1991) or from the white dwarf luminosity function (Wood 1992; see also Adams & Laughlin 1996). The white dwarf luminosity function can be fit using a constant star formation rate. With some uncertainty, the chemical evolution models indicate a slowly decreasing star formation rate with $`qt_{\mathrm{disk}}`$ $``$ $`\mathrm{ln}4`$. Since we want to determine the largest fraction of star formation that can take place within clusters, we use the decreasing form with $`q\mathrm{ln}4/t_{\mathrm{disk}}`$ 0.14 Gyr-1 as our standard case.

If the average cluster mass is $`M_C`$, and if all clusters are eventually slated for destruction, then clusters contribute a fixed amount $`\mathrm{\Delta }\mathrm{\Sigma }`$ to the surface density of the galactic disk,

$$\mathrm{\Delta }\mathrm{\Sigma }=t_{\mathrm{disk}}M_C\frac{1}{qt_{\mathrm{disk}}}\left\{\mathrm{e}^{qt_{\mathrm{disk}}}1\right\}.$$
(3)
For comparison, the observed surface density of the galactic disk is 26.4 $`M_{}`$ pc-2 in visible stars, with an additional 18.2 $`M_{}`$ pc-2 in stellar remnants (Binney & Tremaine 1987; hereafter BT87). Correcting for mass loss in the transformation between progenitor stars and stellar remnants, we obtain a total stellar surface density of $`\mathrm{\Sigma }_{}63`$ $`M_{}`$ pc-2. The fraction $`f_C`$ of the stellar disk component contributed by open clusters is thus given by

$$f_C=\frac{\mathrm{\Delta }\mathrm{\Sigma }}{\mathrm{\Sigma }_{}}7\times 10^5\left(\frac{M_C}{1M_{}}\right)\frac{1}{qt_{\mathrm{disk}}}\left\{\mathrm{e}^{qt_{\mathrm{disk}}}1\right\}.$$
(4)
In the BC91 sample, the observed cluster luminosity function implies a typical cluster mass of $`M_C`$ $``$ 500 $`M_{}`$ (which is subject to some uncertainty due to the required transformation from luminosity to mass). Using this result and $`qt_{\mathrm{disk}}`$ = $`\mathrm{ln}4`$, we find the fraction of stars that form in clusters: $`f_C0.077`$. Thus, open clusters contribute $`8\%`$ of the field stars. Although this fraction is substantial, open clusters are not the birth places for the majority of stars (as sometimes claimed). If we consider the limiting case of a constant star formation rate, $`q0`$, equation  implies a smaller cluster fraction $`f_C`$ 0.035.

The above analysis uses results integrated over the age of the galactic disk. We can obtain a consistency check by using the present day values. The cluster formation rate  and $`M_C500M_{}`$ jointly imply a star formation rate in clusters $`(SFR)_C`$ 225 $`M_{}`$ kpc-2 Myr-1, whereas the current star formation rate in the solar neighborhood is substantially larger, $`(SFR)_T`$ $``$ 3000 – 5000 $`M_{}`$ kpc-2 Myr-1 (Rana 1991). These present day values thus indicate that the fraction $`f_C`$ of star formation that takes place in clusters lies in the range 0.045 $`<f_C<`$ 0.075, consistent with the time integrated estimates found previously.

This result for the fraction of star formation in open clusters ($`f_C<0.1`$) is consistent with previous results. For example, Roberts (1957) estimated that 10% of stars born in the galaxy are formed within exposed clusters. Elmegreen & Clemens (1985) suggest that 10% of low mass clouds form bound clusters, with the remaining clouds producing a distributed stellar population. Similarly, Lada (1999) suggests that most embedded clusters emerge from molecular clouds as unbound systems; although most stars may form in embedded “clusters”, the majority of these stellar groups evolve to become unbound associations rather than bound open clusters.

In this discussion, the implied definition of a cluster requires that the system is big and bright enough to be included in the BC91 analysis. The survey has a limiting absolute V magnitude of –4.5, which means that only young clusters ($`\mathrm{log}_{10}[t/\mathrm{yr}]`$ $`6.5`$) more massive than $`500M_{}`$ are directly included. However, BC91 use a mass function for the clusters that extends down to $`100M_{}`$ to derive their cluster formation rate. If this extrapolation were completely accurate, then the cluster formation rate given by equation  would include all clusters larger than $``$100 $`M_{}`$. Nonetheless, selection effects could lead to an underestimate of the number of the smallest clusters. The value of $`f_C`$ derived here should thus be considered as the fraction of “large clusters” and may not include all clusters with only a couple hundred members. An effective boundary (in $`N`$) thus separates the larger clusters included in the observational surveys from smaller systems which are not included; for the BC91 survey, this boundary is at $`N300`$.

As we show next, sufficiently small stellar systems have dynamical relaxation times that are shorter than their formation times. As a result, such small groups of stars are not clusters in a practical sense. These groups differ from clusters not only in their numbers $`N`$ of members but also in their physical properties: Clusters evolve slowly and experience both a collisionless and an interactive phase; groups evolve quickly and have no collisionless phase. Young groups and clusters thus represent different types of physical systems (see §3) and have different effects on star formation (see §4).

## 3 MINIMUM NUMBER OF CLUSTER STARS

In this section, we estimate the minimum number $`N_{}`$ of stars required for a stellar system to behave as a cluster. As a general rule, a system of stars cannot evolve as a cluster unless its relaxation time $`t_{\mathrm{relax}}`$ is sufficiently long. For the sake of definiteness, we require the relaxation time to be longer than the cluster formation time $`t_{\mathrm{form}}`$, i.e.,

$$t_{\mathrm{relax}}\alpha t_{\mathrm{form}},$$
(5)
where $`\alpha `$ is a dimensionless number greater than unity. We thus want to find the minimum number $`N_{}`$ of stars necessary to satisfy inequality .

The dynamical relaxation time, which depends on the cluster size $`N`$, is the time required for a cluster member to change its velocity by a relative amount of order unity. For cluster ages shorter than the relaxation time, $`t<t_{\mathrm{relax}}`$, the stars in a cluster do not interact appreciably. For longer times, $`t>t_{\mathrm{relax}}`$, the effects of interactions add up and the cluster alters its structure. On longer times, $`tt_{\mathrm{relax}}`$, severe structural changes are forced upon the cluster due to stellar loss through evaporation, ejection, and core collapse. The relaxation time varies over the structure of the cluster. One usually adopts the median relaxation time – that evaluated near the median radius – as a characteristic time scale for the system. As a reference point, the typical evaporation time is of order 100 median relaxation times (BT87). The relaxation time can be written in the form

$$t_{\mathrm{relax}}=Q_{\mathrm{relax}}t_{\mathrm{cross}},$$
(6)
where the crossing time is $`t_{\mathrm{cross}}`$ = $`R/v`$ and where $`Q_{\mathrm{relax}}`$ is the number of crossings required to make the velocity of a star change by a relative amount comparable to unity. The velocity $`v`$ is related to the cluster size $`R`$ through the depth of the cluster potential well, i.e.,

$$v^2=\frac{GM}{R},$$
(7)
where $`M`$ is the total mass of cluster, including both stars and any remaining gas. For a purely stellar system, $`Q_{\mathrm{relax}}`$ $`N/10\mathrm{ln}N`$ (BT87). For the present application, however, we generalize this result to include the presence of cluster gas. We define $`ϵ`$ to be the star formation efficiency of the cluster, i.e., $`ϵ`$ = $`Nm_{}/M`$, where $`m_{}`$ is the mean stellar mass and where $`N`$ is the number of cluster stars. The number of crossings per relaxation time is then given by

$$Q_{\mathrm{relax}}\frac{Nϵ^2}{10\mathrm{ln}[N/ϵ]}.$$
(8)
We want to find a lower limit $`N_{}`$ for the number of stars required for the system to behave as a cluster. For a purely stellar system with $`ϵ=1`$, we see immediately that for sufficiently small numbers $`N`$ of stars, $`Q_{\mathrm{relax}}`$ $`<1`$, and the relaxation time is less than the crossing time. This critical number $`N_1`$ of stars (that required for $`N/10\mathrm{ln}N`$ = 1) is $`N_136`$. A firm lower bound on the minimum number $`N_{}`$ of stars required for a system to evolve as a cluster is thus given by $`N_{}>N_136`$.

We can find a more interesting limit by considering the formation of the cluster and hence by invoking the constraint of equation . Unfortunately, we do not yet have a well developed theory of cluster formation. To make a start, we assume that a cluster forms out of the collapse of a molecular cloud. The formation time must be somewhat longer than the sound crossing time of the system (Elmegreen 2000), so we write

$$t_{\mathrm{form}}\beta \frac{r_{\mathrm{}}}{a},$$
(9)
where $`a`$ is the effective sound speed, $`r_{\mathrm{}}`$ is the initial size of the cloud fragment that forms the cluster, and $`\beta `$ is a dimensionless number larger than (but of order) unity. To obtain a simple model of the initial cloud, we assume that it takes the form of an isothermal sphere. In this case, the initial cloud size $`r_{\mathrm{}}`$ is related to the cluster mass $`M`$ and sound speed through the relation $`r_{\mathrm{}}=GM/2a^2`$ (Shu 1977). Notice that any departures from the isothermal model can be incorporated into the parameter $`\beta `$. As an alternate model, e.g., we could specify the cluster formation time using the free fall collapse time of a uniform density gaseous sphere. With the proper choice of $`\beta `$, however, we recover exactly the same mathematical form for our derived constraint \[$`\beta `$(uniform sphere) = $`(4/\pi )\gamma ^{3/2}`$ $`\beta `$ (isothermal sphere), where $`\gamma `$ is defined in equation  below\]. As a rough estimate, the sound crossing time of the cluster is $`10^6`$ yr, an order of magnitude longer than the time scale ($`10^5`$ yr) for individual star formation events (Myers & Fuller 1993; Adams & Fatuzzo 1996). Keep in mind that $`r_{\mathrm{}}`$ is the size of the initial mass distribution and hence is larger than the size $`R`$ of the newly formed cluster. We thus define an additional dimensionless parameter

$$\gamma \frac{r_{\mathrm{}}}{R},$$
(10)
where we expect $`\gamma `$ to be larger than (but of order) unity.

Putting all of the above conditions together (eqs. ), we derive the constraint

$$\frac{N/ϵ}{\mathrm{ln}[N/ϵ]}10\sqrt{2}\alpha \beta \gamma ^{3/2}ϵ.$$
(11)
This constraint (eq. ) grows weaker as the star formation efficiency $`ϵ`$ decreases; for sufficiently low values of $`ϵ`$, the limit becomes too weak to provide a meaningful bound. In order for the cluster to remain gravitationally bound, however, the star formation efficiency cannot become too small. Before analyzing this compromise in detail (see below), we obtain a rough estimate using $`ϵ=1/2`$ as a typical value. We also adopt $`\alpha =1`$ which corresponds to the relaxation time and the cluster formation time being equal. If we had a definitive theory of cluster formation, the remaining parameters $`\beta `$ and $`\gamma `$ (or their distributions of allowed values) would be unambiguously specified. In the absence of a complete theory, we must rely on estimates and hence we adopt $`\beta `$ = 2 = $`\gamma `$. For this case, the limit becomes $`N108`$. Thus, the minimum number $`N_{}`$ of stars required for the cluster formation time to be shorter than the dynamical relaxation time is $`N_{}100`$.

We now derive a more rigorous constraint on $`N_{}`$ by considering the whole range of star formation efficiencies. Virial arguments suggest that clusters remain bound if $`ϵ>1/2`$ and become unbound if $`ϵ<1/2`$ (for rapid gas removal; see Hills 1980, Mathieu 1983, Elmegreen 1983). In practice, however, the stars in a forming cluster have a distribution of velocities. The low velocity stars in the tail of the distribution survive as a gravitationally bound entity even if $`ϵ<1/2`$; the high velocity stars in the opposite tail escape even if $`ϵ>1/2`$. A more accurate description is that a cluster formed with $`N`$ stars (before gas removal) eventually produces a bound cluster with $`N_f=_{}(ϵ)N`$ stars after gas removal. The function $`_{}(ϵ)`$ varies smoothly with star formation efficiency rather than exhibiting step function behavior. The shape of $`_{}(ϵ)`$ depends on the shape of the distribution function for the cluster stars, the rate of gas removal, and the density distributions of the stars and gas (e.g., Adams 2000; Kroupa, Aarseth, & Hurley 2000; Geyer & Burkert 2000; Lada, Margulis, & Dearborn 1984). For clusters with isotropic velocity distributions, for example, the function $`_{}(2ϵϵ^2)`$ provides a good approximation over the expected range of cluster models (see Fig. 3 of Adams 2000).

To incorporate limits on the star formation efficiency into our analysis, we require that the final bound cluster (which contains $`N_f`$ stars after gas removal) have its relaxation time longer than its crossing time, i.e., $`N_fN_1`$, where $`N_136`$ is defined above. We thus impose the additional constraint

$$_{}(ϵ)NN_1.$$
(12)
For a given function $`_{}(ϵ)`$, the coupled constraints  and  define a well posed optimization problem.

The solution is straightforward. The first constraint  says that $`N>f_1(ϵ)`$ = $`\lambda ϵ^2\mathrm{ln}[f_1/ϵ]`$, where $`\lambda 10\sqrt{2}\alpha \beta \gamma ^{3/2}`$ and where the function $`f_1`$ is defined implicitly (notice that $`f_1`$ is not defined for extremely small values $`ϵ<\mathrm{e}/\lambda 0.03`$). This function $`f_1(ϵ)`$ is a monotonically increasing function of the variable $`ϵ`$. Similarly, the second constraint  says that $`N>f_2(ϵ)=N_1/_{}(ϵ)`$, where the function $`f_2(ϵ)`$ is a monotonically decreasing function of $`ϵ`$. Since $`N`$ must be greater than both $`f_1`$ and $`f_2`$ for all values of $`ϵ`$, a lower bound on $`N`$ occurs at the crossover point where $`f_1(ϵ)=f_2(ϵ)`$. This lower bound on $`N`$ is a solution to the equation

$$N=\lambda \left[_{}^1(N_1/N)\right]^2\mathrm{ln}\left[\frac{N}{_{}^1(N_1/N)}\right].$$
(13)
This bound holds for all values of the star formation efficiency. To evaluate this bound, we only need to specify the parameter $`\lambda `$ (which encapsulates our uncertainties regarding cluster formation) and the function $`_{}`$ (which is determined by the escape of stars from the cluster during gas removal). For a representative case of $`\lambda =80`$ ($`\alpha =1`$; $`\beta =2=\gamma `$) and $`_{}=(2ϵϵ^2)`$, we plot the resulting curves $`f_1(ϵ)`$ and $`f_2(ϵ)`$ in Figure 1. The intersection point determines the constraint $`N>N_{}58`$.

Possible uncertainties in our lower bound $`N_{}`$ arise from the form of the function $`_{}(ϵ)`$ and the value of $`\lambda `$. Fortunately, the result is relatively insensitive to these choices. If we use alternate fits for the function $`_{}`$ (Adams 2000), we obtain essentially the same bound. For example, the cruder approximation $`_{}`$ = $`\sqrt{ϵ}`$ changes the crossover point only by $`\mathrm{\Delta }N_{}`$ $``$ 1. Our result also depends only weakly on the value of $`\lambda `$; in the limit $`N_1N`$, the bound obeys the scaling law $`N_{}\lambda ^{1/3}`$ (up to a logarithmic correction), so the constraint is not overly sensitive to $`\lambda `$. Figure 1 illustrates this property by plotting alternate curves for $`f_1(ϵ)`$ using $`\lambda `$ = 40 and 160. The intersections occur at $`N`$ = 47 and 73, respectively, which are close to the values predicted by the $`\lambda ^{1/3}`$ scaling law. Given these uncertainties and the sharpness of the minimum shown in Figure 1, we adopt $`N_{}`$ = 100 as the effective lower limit on the number of stars required for a system to be a cluster.

This analysis defines a critical value $`ϵ_C`$ of the star formation efficiency, i.e., the value corresponding to the crossover point $`f_1(ϵ)`$ = $`f_2(ϵ)`$. For the typical case defined above, this critical value $`ϵ_C0.4`$. For systems with high star formation efficiency, $`ϵ>ϵ_C`$, the constraint of equation  dominates and the minimum number of cluster stars is determined by making the dynamical relaxation time sufficiently long. For systems with low star formation efficiency, $`ϵ<ϵ_C`$, cluster survival depends on having enough stars remaining gravitationally bound after gas removal, as enforced by equation . Notice also that the constraints shown in Figure 1 define a relatively sharp minimum value of $`N`$. For values of the star formation efficiency $`ϵ`$ that depart from the critical value $`ϵ_C`$, the constraints on $`N`$ are considerably more restrictive.

The constraints on $`N_{}`$ derived here depend on both cluster formation parameters and the function $`_{}(ϵ)`$, which, in turn, depends on the gas structure of the cluster and the distribution of stellar velocities. We have used basic considerations of cluster formation and typical parameter values to estimate $`N_{}`$. However, all steps of this calculation can be improved. Although the value of $`N_{}`$ is thus subject to some uncertainty, the existence of a limiting value $`N_{}`$ is not in question – stellar systems with too few members will not behave like clusters in a dynamical sense. Such groups lack a collisionless phase and quickly evolve toward evaporation and dispersal.

The minimum number of stars $`N_{}`$ $``$ 100 is approximately the same as that of the smallest clusters considered in the analysis of the BC91 survey. This result also makes sense: Systems with $`NN_{}`$ evolve so rapidly and are so dim that they would have little chance of being included in an observational survey of this type. As a result, such observational surveys come close to providing a realistic distribution of the stellar systems that can rightfully be considered as open clusters (according to this constraint). We can use the typical age of open clusters ($``$100 Myr; BC91, BT87) to find another consistency check on this argument: Because a typical cluster lasts for $``$100 relaxation times before significant evaporation (BT87), the cluster must have a relaxation time of $``$1 Myr. This requirement, in turn, implies a lower bound of $`N160`$. These considerations thus suggest that stellar systems with $`N<N_{}100`$ are not true clusters in that they evolve and disperse much more quickly than open clusters. These groups should be considered as a different type of astronomical system.

## 4 EFFECTS OF GROUPS AND CLUSTERS ON STAR FORMATION

In the previous sections, we made a dynamical distinction between clusters (large systems with $`N>N_{}`$) and smaller groups ($`N<N_{}`$). We now discuss different ways that groups and clusters affect the star formation process. In particular, we consider scattering interactions involving cluster/group members, supernova explosions, and the background ultraviolet radiation field provided by the cluster (these results are summarized by Fig. 2).

### 4.1 Dynamics and Stellar Scattering Interactions

To illustrate the different dynamical effects that groups and clusters exert upon star formation, we first consider the early evolution of a stellar system near our estimated boundary at $`N`$ = $`N_{}`$ = 100 members. If the cluster forms out of the collapse of a large molecular cloud core with size $`R_0`$ = 1 pc and effective sound speed $`a`$ = 1 km/s, the formation time is a few million years. With $`N`$ = 100, the initial value of $`Q`$ = 7.5 if we assume $`ϵ=1/2`$ so that the cluster is half gas and half stars. With $`v`$ 1 km/s and the typical stellar mass $`m_{}`$ $``$ 0.5 $`M_{}`$, the total mass is initially 100 $`M_{}`$ and the virial size (given by eq. ) is $`R=GM/v^2`$ = 0.43 pc (about $`R_0`$/2). The crossing time in this state is 0.43 Myr and the relaxation time is 3.2 Myr. Because the cluster takes 1–2 Myr to form and another 1–2 Myr to disperse its gas content, it experiences only about one relaxation time while it remains embedded. After gas removal, the cluster retains $`N`$ 75 stars (using $`_{}`$ = $`2ϵϵ^2`$) and the relaxation time drops to 0.25 Myr. After 25 Myr of additional evolution, 100 times this starting relaxation time, the group loses most of its members through evaporation and becomes highly compromised as a stellar system. It would be impossible to observationally identify as a cluster. Over its $``$25 Myr of evolution, the group has an average stellar density less than 50 pc-3.

Given the above evolutionary picture of our transition-sized cluster, we consider the possible effects that the system has on its constituent forming stars and young solar systems. During the formation stage of the cluster, the typical separation between forming stars is $`\mathrm{\Delta }r(4\pi /3N)^{1/3}R`$ 0.35 pc. If we assume that individual stars form through the inside-out collapse of a centrally condensed structure (as in Shu 1977; Adams, Lada, & Shu 1987), the region containing a typical stellar mass ($`0.5M_{}`$) extends over $`r_{\mathrm{}}`$ = 0.027 pc, $`13`$ times smaller than the mean separation (where we assume $`a=0.20`$ km/s for the individual infall region). The mean separation between forming stars is much larger than the size of their protostellar infall regions and interaction effects are minimal. Similarly, the tidal radius due to the tidal forces exerted on an individual star forming site by the gravitational potential of the background cluster is given by $`r_T=\eta (M_{}/M_{\mathrm{clust}})^{1/3}R_0`$, where the constant $`\eta `$ depends on the geometry of the region. For $`R_0`$ = 1 pc, $`M_{}=0.5M_{}`$, $`M_{\mathrm{clust}}`$ = 100 $`M_{}`$, and $`\eta `$ = 1, we thus obtain $`r_T`$ $``$ 0.17 pc, which is much larger than the size $`r_{\mathrm{}}`$ of a protostellar infall region. Tidal effects are thus small in these transition-sized groups.

An important channel for clusters to affect star formation is through stellar encounters within the cluster. Such encounters could lead to binary capture, disk disruption, or changes in planetary orbits. These effects require young solar systems to experience disruptive close encounters; keep in mind that all solar systems experience more distant encounters that lead to dispersal of the cluster. Let $`\sigma _{200}`$ denote the cross section for a close encounter in units of (200 AU)2 or $`9\times 10^7`$ pc2; this is a typical cross section required for an encounter to force binary capture or to strongly disrupt a young solar system (e.g., Heller 1993; Ostriker 1994; Laughlin & Adams 1998, 2000; Kroupa, Petr, & McCaughrean 1999). For our own solar system, for example, a cross section of (200 AU)2 is the value required to eject Neptune, give Uranus an orbital eccentricity $`e>0.75`$, and/or randomize the orbital inclination angles of the giant planets (Adams & Laughlin 2001). In our transition-sized cluster, the probability $`P_D`$ for a disruptive close encounter is given by $`P_Dn\sigma _{200}v(\mathrm{\Delta }t)`$, where the mean density $`n`$ 50 pc-3, $`v1`$ km/s, and $`\mathrm{\Delta }t25`$ Myr. The probability is thus $`P_D`$ $`10^3`$ $`\sigma _{200}`$; the corresponding odds of a disruptive encounter taking place within the expected lifetime of the cluster is only $``$1 in 1000 (for $`\sigma _{200}`$ = 1). Since the cluster contains only 100 stars, the chances of any disruptive close encounters occurring are about 1 out of 10. As a rule, more distant encounters disperse the cluster before disruptive close encounters can greatly alter the constituent solar systems.

The above considerations suggest that stars forming within transition-sized clusters ($`N=N_{}=100`$) experience minimal dynamical effects from their cluster environment. Interactions between protostellar infall regions are rare and tidal influences are small. After solar systems are made, the odds of binary capture or severe disruption are low, only about 1 part in 1000 per star (1 out of 10 per cluster). Stars forming within larger clusters ($`NN_{}`$) experience the aforementioned dynamical effects with high probability, whereas stars forming within smaller groups are less likely to experience such effects. As a result, our estimated boundary between groups and clusters (at $`N_{}`$ $``$ 100) also represents an effective boundary between stellar systems that have a dynamical impact on forming solar systems (clusters) and smaller systems that do not (groups).

We can illustrate the transition from groups to clusters by deriving a rough scaling law for $`P_D`$. The probability $`P_D`$ for a close encounter depends on the mean stellar density $`n`$ and the total lifetime $`\mathrm{\Delta }t`$ of the system. Using equations , we find the scaling laws $`n`$ $``$ $`N/R^3`$ and $`\mathrm{\Delta }t`$ $``$ $`(R/v)N/\mathrm{ln}N`$. Putting these results together, we find that the probability for disruptive encounters takes the form $`P_DN^\mu `$, where the index $`\mu 2`$. According to this relation, stars living in smaller systems with $`N=30`$ are less likely to experience disruptive encounters by an order of magnitude ($`P_D10^4`$). In these small systems, chances are good (about 1 out 300) that no stars ever experience disruptive encounters. For larger systems with $`N=300`$, the probability (per star) of disruptive close encounters increases to $`P_D10^2`$ (1%). For still larger systems with $`N=1000`$, the probability of a disruptive encounter becomes significant, $`P_D0.1`$. In these latter systems, perhaps 100 out of the 1000 cluster members could experience significant disruption through a close encounter (see Fig. 2).

This discussion implicitly assumes that the stellar systems under consideration are large enough (in $`N`$) so that we can make statistical arguments for the evolution (in phase space) of the individual stars. In sufficiently small stellar systems, however, the dynamics of a given star is dominated by a few close encounters rather than many distant encounters and hence the scaling laws used here no longer apply. The criterion for most of the scattering to be due to weak encounters is $`Rb_{\mathrm{min}}=Gm_{}/v^2`$, where $`R`$ is the cluster size and $`m_{}`$ is the typical stellar mass. Using equation , this requirement reduces to $`N1`$, which is met by most systems. Using standard formalism (BT87), we can make this requirement more quantitative: The typical velocity perturbation for a stellar encounter at impact parameter $`b`$ is given by $`\delta v/v`$ $``$ $`b_{\mathrm{min}}/b`$. Suppose we want to calculate the group size $`N`$ such that the velocity perturbations are small, specifically $`\delta v/v<\delta _0`$, for at least half of the encounters (for a given choice of $`\delta _0`$). Half of the encounters have impact parameters $`b<b_{1/2}`$ where $`b_{1/2}`$ is given by $`\mathrm{ln}(R/b_{1/2})=0.5\mathrm{ln}(R/b_{\mathrm{min}})`$, i.e., $`b_{1/2}=(Rb_{\mathrm{min}})^{1/2}`$. We thus require $`b_{\mathrm{min}}/b_{1/2}`$ $`<`$ $`\delta _0`$, which implies $`b_{\mathrm{min}}<\delta _0^2R`$. Using the definition of $`b_{\mathrm{min}}`$ and equation , this constraint can be rewritten in terms of cluster size: $`N\delta _0^2>1`$. For example, if we want at least half of the stellar encounters to have velocity perturbations $`\delta v/v`$ $`<`$ 0.20 = $`\delta _0`$, we need a cluster size of at least $`N>\delta _0^2`$ = 25; similarly, if we require half the encounters to have $`\delta v/v`$ $`<`$ 0.30 = $`\delta _0`$, we would need $`N>\delta _0^211`$. If $`N`$ is too small, the evolution of the stellar aggregate is dominated by a few strong encounters and a statistical description (based on many weak encounters) breaks down. This result implies an effective lower boundary for groups; although the boundary is not sharp, systems with $`N<10`$ are dominated by a few hard collisions and hence their behavior depends sensitively on the initial values of the phase space variables (see Retterer 1979, BT87; for dynamical simulations of small ($`N<10`$) systems, see Sterzik & Durisen 1995, Bonnell et al. 1997).

### 4.2 The Probability of Supernovae

Another way for a star forming environment to affect its constituent stars is through supernova explosions. These energetic events can disrupt star forming regions and remove gas from young clusters; supernovae have also been invoked as a way to trigger star formation (e.g., Cameron & Truran 1977; Boss & Foster 1998). In the present context, we argue that groups have little chance of experiencing a supernova, whereas clusters will often be subjected to their destructive effects.

To support the above claim, we find the probability $`P_{SN}`$ that a stellar system will be subjected to a supernova explosion, as a function of the number $`N`$ of stars in the system. Only stars more massive than $`M_{SN}`$ 8 $`M_{}`$ explode at the end of their nuclear burning lives. The fraction $`f_{SN}`$ of stars that are massive enough to explode ($`M_{}>M_{SN}`$) depends on the stellar IMF and is $`f_{SN}0.004`$ (see Binney & Merrifield 1998 for a discussion of values for $`M_{SN}`$ and $`f_{SN}`$). Next we assume that the IMF is independent of the size $`N`$ of the stellar system. To calculate $`P_{SN}`$, we imagine picking $`N`$ stars at random from the stellar mass distribution. The probability that a given star will not be massive enough to explode is $`\stackrel{~}{p}=1f_{SN}0.996`$. The probability that a system of $`N`$ stars will not contain an exploding star is thus $`(\stackrel{~}{p})^N`$. Finally, the probability that a stellar aggregate (with $`N`$ members) does contain a progenitor star massive enough to explode is given by

$$P_{SN}(N)=1\stackrel{~}{p}^N=1\left[1f_{SN}\right]^N.$$
(14)
Equation  gives the likelihood for supernovae to occur in a stellar system of size $`N`$. For stellar groups, as defined in §3 with $`N<N_{}=100`$, the probability of a supernova is low: $`P_{SN}<0.33`$. A natural break-even point between systems with supernovae and those without occurs where $`P_{SN}(N_{SN})`$ = 0.5; the critical number of stars is $`N_{SN}170`$. Larger clusters thus have an appreciable chance of containing stars large enough to explode as supernovae. Keep in mind that the boundary is not perfectly sharp – stellar aggregates follow a smooth probability distribution (given by eq. ; see Fig. 2).

For supernovae to affect star or planet formation, the stellar system must live long enough for massive stars to develop iron cores and then explode; we thus need a cluster lifetime $`\mathrm{\Delta }t>10`$ Myr. As discussed in §4.1, a transition-sized cluster with $`N=N_{}=100`$ is expected to live for $`\mathrm{\Delta }t25`$ Myr. As a result, clusters that are large enough (in $`N`$) to contain massive stars with high probability are also sufficiently long-lived for their massive stars to evolve and die while the cluster remains intact.

To summarize, small stellar systems are unlikely to have stars large enough to explode as supernovae. The boundary between small systems with no supernovae and larger systems with supernova explosions is about $`N_{SN}170`$. This boundary is roughly coincident with the boundary between groups and clusters ($`N_{}100`$) as defined by dynamical considerations (§3). As a result, stellar groups will not generally contain supernovae, and will not be subjected to their disruptive effects nor the possibility of supernova triggers for star formation. Larger clusters often contain supernovae and can experience both their destructive and (possibly) constructive effects. Throughout this discussion, we assume that the masses of forming stars obey an IMF that is independent of the system size; statistics alone then imply that clusters generally contain massive stars whereas groups generally do not. This distinction becomes even sharper if massive stars form preferentially within larger clusters as some authors have conjectured (e.g., Testi, Palla, & Natta 1999).

### 4.3 UV Radiation Fields Provided by Clusters vs Groups

External radiation fields from the background environment (the group or cluster) can have a substantial impact on the star formation process. For example, radiation fields can remove gas from circumstellar disks and thereby suppress both disk accretion and planet formation (Shu, Johnstone, & Hollenbach 1993; Hollenbach et al. 1994; Störzer & Hollenbach 1999). External radiation fields can also play a role in ending the protostellar infall phase (e.g., Hester et al. 1999). These processes are driven mostly by the ultraviolet (UV) portion of the radiation field, which is dominated by the most massive stars in the system. As in the case of supernovae, the shape of the stellar IMF dictates that massive stars are rare except in sufficiently large systems. As a result, solar systems forming within large stellar aggregates (clusters) receive an appreciable contribution of UV radiation from their background cluster; solar systems forming within small stellar groups receive little UV radiation from the background.

To substantiate this claim, we estimate the UV radiation field provided by a stellar aggregate, as a function of the number $`N`$ of stars in the system. We follow a previous calculation of the UV field for the expected conditions experienced by our own solar system during its planet formation epoch (Adams & Laughlin 2001). This calculation finds the expectation value for the ionizing ultraviolet flux from a background stellar system, i.e.,

$$F_{uv}=1.6\times 10^{12}\mathrm{cm}^2\mathrm{sec}^1\left(\frac{N}{2000}\right)\left(\frac{R}{1\mathrm{p}\mathrm{c}}\right)^2,$$
(15)
where $`N`$ is the number of stars and $`R`$ is the cluster size. This expression was obtained by making two averages: We first integrate over a typical stellar orbit through the cluster to find the mean flux impinging upon a given solar system due to the massive stars, which provide the UV flux and are assumed to reside at the center. We also integrate over the stellar IMF, weighted by the UV luminosity as a function of stellar mass, to find the ionizing UV flux as a function of the number $`N`$ of stars in the system (Adams & Laughlin 2001). For comparison, the ionizing UV luminosity of a 1 $`M_{}`$ star during its pre-main-sequence phase cannot be larger than about $`L_{uv}10^{41}`$ sec-1 (Gahm et al. 1979). The corresponding UV flux from the star is given by

$$F_{uv}=3.5\times 10^{13}\mathrm{cm}^2\mathrm{sec}^1\left(\frac{r_{}}{1\mathrm{A}\mathrm{U}}\right)^2,$$
(16)
where $`r_{}`$ is the radial coordinate centered on the star.

One measure of the importance of the local UV radiation background is the total number of photons intercepted by circumstellar disks. These disks actively form planets during their first $`10`$ Myr (Lissauer 1993). Ultraviolet radiation acts to drive a wind from the disk surface, remove gas from the disk, and eventually compromise planet formation. For a typical disk, we want to compare the number of ionizing UV photons provided by its background cluster (eq. ) with the UV radiation intercepted from its central star (eq. ). The disk is embedded in the background UV radiation field and both sides are exposed; the disk thus receives UV photons from the cluster at a rate

$$\mathrm{\Phi }_{uv}=2\pi R_\mathrm{d}^2F_{uv}10^{41}\mathrm{sec}^1\left(\frac{N}{100}\right)\left(\frac{R}{1\mathrm{p}\mathrm{c}}\right)^2,$$
(17)
where $`R_\mathrm{d}`$ $``$ 30 AU is the radial size of the disk. The nominal value of $`\mathrm{\Phi }_{uv}`$ (which has been rescaled to $`N=100`$) is equal to the total production rate of ionizing UV photons from a 1 $`M_{}`$ star ($`10^{41}`$ sec-1). Sufficiently large clusters (those with $`N>100`$) produce enough ionizing UV radiation to dominate the UV flux intercepted by circumstellar disks; smaller groups ($`N<100)`$ have smaller UV backgrounds and circumstellar disks are primarily irradiated by their central stars (see Fig. 2).

We now derive the break-even point between stellar systems that are large enough to dominate the ionizing UV field experienced by a circumstellar disk and smaller systems in which individual stars provide most of the UV to their disks. In general, disks intercept only a fraction of the UV photons generated by their central stars. Taking into account both disk flaring and scattered (diffuse) photons, Shu et al. (1993) estimate that nearly 50$`\%`$ of the UV photons are intercepted by the disk. The rate of intercepted (intrinsic) stellar UV photons is thus $`\mathrm{\Phi }_{}5\times 10^{40}`$ sec-1 for a 1 $`M_{}`$ star. The scaling law  implies that stellar systems with $`N>50`$ provide more ionizing UV photons to a circumstellar disk than its central star. According to this criterion of ionizing UV radiation, the boundary between large stellar systems (clusters) that affect circumstellar disks and smaller systems (groups) that do not is $`N50`$.

This result was derived using an expectation value for the mean UV flux from a stellar aggregate. With cluster/group sizes as small as $`N=50100`$, the IMF is not completely sampled by any given stellar system and individual groups/clusters will experience sizable fluctuations about the expectation value (eq. ). In addition, equation  was derived without taking into account attenuation of the UV flux by gas and dust in the cluster. Since the gas removal time is relatively short (a few Myr) compared to the expected disk lifetimes and the planet formation times ($`10`$ Myr), an attenuation correction would not appreciably change the estimated boundary at $`N=50100`$. Finally, we have only considered the ionizing radiation in this discussion. Photodissociation can be as important as photoionization and should be studied in a more complete treatment (see Diaz-Miller, Franco, & Shore 1998; Störzer & Hollenbach 1999).

## 5 CONCLUSION

### 5.1 Summary of Results

In this paper, we have obtained three principal results:

 Using the estimated cluster formation rate, we find that less than 10$`\%`$ of stars are formed in systems that become open clusters (§2). In this context, only relatively large systems (with a few hundred stars or more) are considered to be bona fide clusters. Of the nearly 90% of star formation that takes place in other environments, perhaps the majority takes place in smaller systems – denoted here as groups – with $`N=10100`$ stars.

 Viable clusters – those which survive to be observable as open clusters – must have a minimum number $`N_{}`$ of stars. We estimate this limiting value by requiring that the cluster formation time is shorter than the dynamical relaxation time (§3) and thereby find the minimum number of cluster stars: $`N_{}100`$. This critical value $`N_{}`$ provides an effective boundary between clusters ($`N>N_{}`$) and smaller stellar groups ($`N<N_{}`$).

 Small stellar groups and larger clusters affect star formation in different ways (§4). The small fraction ($``$1/10) of stars that form within cluster systems are likely to be affected by disruptive close encounters (§4.1), such as stellar mergers, binary capture via disks, disruption of planetary orbits, strong interactions between protostars, and competitive accretion. Stars forming within clusters are also influenced by nearby massive stars, through both supernova explosions (§4.2) and background UV radiation fields (§4.3). For the majority of stars ($``$9/10), however, none of the aforementioned effects operate with high probability and the background environment is sufficiently diffuse to allow individual stars to form in relative isolation.

### 5.2 Modes of Star Formation

The results of this paper suggest that we must consider star formation to take place in more modes than has been historically recognized. In particular, we need to move beyond the traditional dichotomy between isolated single stars and clusters. The following classes provide a starting point:

\[I\] Isolated single stars and multiple systems, including binaries, triples, and other few body systems. For this class of systems, the number $`N`$ of stellar members is less than about 10, although more typically $`N`$ = 2 or 3. We could subdivide this class further into multiple systems and single stars (thus defining a Class ). Although single stars may well be in the minority, this zeroeth class would include our solar system.

\[II\] Groups, consisting of intermediate numbers of stars with $`10<N<N_{}100`$. The upper limit $`N_{}`$ is determined by the criterion that the dynamical relaxation time must be longer than the formation time of the cluster; the value $`N_{}`$ $``$ 100 is thus approximate and depends on manner in which clusters form. The lower limit at $`N10`$ is also approximate and marks the boundary between few-body systems and larger groups in which stellar dynamics can be described statistically. This mode of star formation may be dominant in that most stars may form in such groups. This work suggests that stars forming within small groups are largely decoupled from their immediate environment.

\[III\] Clusters, consisting of large numbers of stars with $`N>N_{}100`$. This regime corresponds to robust clusters that live for long times, i.e., systems that can be observationally identified as open clusters. Stars forming within these systems are subjected to dynamical effects such as core mergers, binary capture, and planetary scattering. Clusters have large enough stellar membership to contain massive stars with high probability; stars forming in clusters are thus exposed to intense UV radiation and supernova blast waves.

These three modes of star formation represent the types of stellar aggregates that a forming star might find itself within. These modes, or classes, differ from each other in two important respects:  They represent different kinds of stellar systems that exhibit different dynamical behavior and evolution (cf. §2 and §3).  They represent different types of star formation environments that affect the star formation process in different ways (§4). In particular, cluster environments subject their stellar constituents to a host of disruptive (and constructive) influences, whereas smaller groups have relatively little effect on the star formation process.

### 5.3 Discussion

Perhaps the most important result of this paper is to make a clearer distinction between groups and clusters. We have made this distinction in five different ways and thereby obtain five estimates for the boundary $`N_{}`$ between groups ($`N<N_{}`$) and clusters ($`N>N_{}`$):  Clusters are big and bright enough, and live long enough, to be included in observational surveys (§2; $`N_{}300`$).  Clusters have dynamical relaxation times that are longer than their formation times (§3; $`N_{}60100`$).  Clusters are sufficiently dense and long-lived so that disruptive scattering encounters can affect circumstellar disks and their planetary progeny (§4.1; $`N_{}100`$).  Clusters have enough stellar members and live long enough so that supernovae can affect forming stars (§4.2; $`N_{}170`$).  Clusters have enough massive stars so that the ionizing UV radiation field impinging upon forming planetary systems is dominated by the background cluster rather than the central star (§4.3; $`N_{}50100`$). These five determinations are roughly coincident and imply $`N_{}100`$. However, the boundary between groups and clusters is not perfectly sharp: stellar systems exhibit a continuous distribution of properties as a function of stellar membership $`N`$ (see Fig. 2).

These criteria for distinguishing groups and clusters depend on many different (but sometimes coupled) physical processes: the longevity of stellar systems, the stellar IMF, the UV mass-luminosity relationship for stars, the minimum progenitor mass for a supernova, the brightness of stellar populations, the formation time for clusters, and the scattering cross sections for solar system disruption. The complicated interplay between stellar dynamics and stellar physics thus leads to a relatively clean distinction between groups and clusters, both as stellar systems (they behave differently as astronomical entities) and as star formation environments (they affect star formation differently).

Some ambiguity remains in the relative portions of stars that form in the three classes defined here. The observational cluster surveys (e.g., BC91) show that only 10% of star formation takes place within the clusters included in the sample. Although BC91 use a mass function for clusters that extends down to $`m_1`$, where $`40M_{}m_1120M_{}`$ for their various models, the accounting for the lowest mass clusters is not necessarily complete. The observational surveys could thus be missing small clusters in the range 100 $`N`$ 300. One common explanation for why stars are not seen in open clusters is that they form within cluster-sized units, but the units disperse quickly after gas removal. If this scenario is true, then the cluster-sized units must be small enough to avoid detection in the surveys (which include young clusters with $`t`$3 Myr and all bright clusters). It thus remains possible for a substantial fraction of star formation to take place in systems with $`N`$ in the range 100 $`N`$ 300. However, such small stellar systems disperse rapidly and are thus more like groups than clusters; in addition, these small systems have relatively little effect on star formation (see §4 and Fig. 2).

The relative amount of star formation that takes place in groups versus the isolated mode must also be better specified. Using the current data base of dense cores mapped in ammonia (Jijina, Myers, & Adams 1999), we find that most of the dense gas is contained in cores associated with stellar aggregates with $`N>30`$ (safely in the regime of groups as defined here). Although the data base clearly shows that more star formation takes place in groups than in the isolated mode, the heterogeneous nature of the data makes further quantitative determinations difficult. This issue thus requires further work.

When star formation takes place within groups, the most common result is complete dispersal of the system on a relatively short time scale, due to gas removal (a few Myr) and dynamical evolution (10–20 Myr). For example, Orion A (L1641) contains both groups and more widely distributed stars. It has been suggested that the distributed population in L1641 may have been born in larger aggregates like those seen now, but the groups have already dispersed (Strom et al. 1994). On the other hand, groups forming in close proximity could merge and thereby build larger clusters after their initial formative stage. The formation and evolution of intermediate sized stellar systems ($`N100`$) thus constitutes an important area for future work.

Another goal for the future is to construct the distribution $`P(N)`$ which gives the probability that a star will form within a system of $`N`$ members. This work indicates that the probability distribution contains three principle components: isolated and few-body systems with $`N<10`$, small groups with $`10<N<N_{}`$ $`100`$ (perhaps containing the majority of forming stars), and a tail of larger stellar systems representing the clusters with $`N>N_{}`$ (comprising about 10% of the population). Although the available data are not adequate to define a reliable distribution at this time, the construction of $`P(N)`$ should become feasible in the near future.

We would like to thank Gus Evrard, Suzanne Hawley, Charlie Lada, Greg Laughlin, and Frank Shu for useful discussions. We also thank an anonymous referee and Scientific Editor Steve Shore for many helpful suggestions that greatly improved the paper. This work was supported by funding from The University of Michigan and by bridging support from NASA Grant No. 5-2869.


format: text

text:
# Benchmarking Optimization Software with Performance Profiles

## 1 Introduction

The benchmarking of optimization software has recently gained considerable visibility. Hans Mittlemann’s Benchmarks work on a variety of optimization software has frequently uncovered deficiencies in the software and has generally led to software improvements. Although Mittelmann’s efforts have gained the most notice, other researchers have been concerned with the evaluation and performance of optimization codes. As recent examples, we cite HYB00; SCB95a; ASB98; BCGT97; CGT96; HM99; RJV99.

The interpretation and analysis of the data generated by the benchmarking process are the main technical issues addressed in this paper. Most benchmarking efforts involve tables displaying the performance of each solver on each problem for a set of metrics such as CPU time, number of function evaluations, or iteration counts for algorithms where an iteration implies a comparable amount of work. Failure to display such tables for a small test set would be a gross omission, but they tend to be overwhelming for large test sets. In all cases, the interpretation of the results from these tables is often a source of disagreement.

The quantities of data that result from benchmarking with large test sets have spurred researchers to try various tools for analyzing the data. The solver’s average or cumulative total for each performance metric over all problems is sometimes used to evaluate performance HYB00; BCGT97; CGT96. As a result, a small number of the most difficult problems can tend to dominate these results, and researchers must take pains to give additional information. Another drawback is that computing averages or totals for a performance metric necessitates discarding problems for which any solver failed, effectively biasing the results against the most robust solvers. As an alternative to disregarding some of the problems, a penalty value can be assigned for failed solver attempts, but this requires a subjective choice for the penalty. Most researchers choose to report the number of failures only in a footnote or separate table.

To address the shortcomings of the previous approach, some researchers rank the solvers BCGT97; CGT96; SGN91; RJV99. In other words, they count the number of times that a solver comes in $`k`$th place, usually for $`k=1,2,3`$. Ranking the solvers’ performance for each problem helps prevent a minority of the problems from unduly influencing the results. Information on the size of the improvement, however, is lost.

Comparing the medians and quartiles of some performance metric (for example, the difference between solver times BCGT97) appears to be a viable way of ensuring that a minority of the problems do not dominate the results, but in our testing we have witnessed large leaps in quartile values of a performance metric, rather than gradual trends. If only quartile data is used, then information on trends occurring between one quartile and the next is lost; and we must assume that the journey from one point to another proceeds at a moderate pace. Also, in the specific case of contrasting the differences between solver times, the comparison fails to provide any information on the relative size of the improvement. A final drawback is that if results are mixed, interpreting quartile data may be no easier than using the raw data; and dealing with comparisons of more than two solvers might become unwieldy.

The idea of comparing solvers by the ratio of one solver’s runtime to the best runtime appears in SCB95a, with solvers rated by the percentage of problems for which a solver’s time is termed very competitive or competitive. The ratio approach avoids most of the difficulties that we have discussed, providing information on the percent improvement and eliminating the negative effects of allowing a small portion of the problems to dominate the conclusions. The main disadvantage of this approach lies in the author’s arbitrary choice of limits to define the borders of very competitive and competitive.

In Section 2, we introduce performance profiles as a tool for evaluating and comparing the performance of optimization software. The performance profile for a solver is the (cumulative) distribution function for a performance metric. In this paper we use the ratio of the computing time of the solver versus the best time of all of the solvers as the performance metric. Section 3 provides an analysis of the test set and solvers used in the benchmark results of Sections 4 and 5. This analysis is necessary to understand the limitations of the benchmarking process.

Sections 4 and 5 demonstrate the use of performance profiles with results EDD00 obtained with version 2.0 of the COPS cops-home test set. We show that performance profiles eliminate the influence of a small number of problems on the benchmarking process and the sensitivity of results associated with the ranking of solvers. Performance profiles provide a means of visualizing the expected performance difference among many solvers, while avoiding arbitrary parameter choices and the need to discard solver failures from the performance data.

We conclude in Section 6 by showing how performance profiles apply to the data Benchmarks of Mittelmann for linear programming solvers. This section provides another case study of the use of performance profiles and also shows that performance profiles can be applied to a wide range of performance data.

## 2 Performance Evaluation

Benchmark results are generated by running a solver on a set $`𝒫`$ of problems and recording information of interest such as the number of function evaluations and the computing time. In this section we introduce the notion of a performance profile as a means to evaluate and compare the performance of the set of solvers $`𝒮`$ on a test set $`𝒫`$.

We assume that we have $`n_s`$ solvers and $`n_p`$ problems. We are interested in using computing time as a performance measure; although, the ideas below can be used with other measures. For each problem $`p`$ and solver $`s`$, we define

$$t_{p,s}=\text{computing time required to solve problem }p\text{ by solver }s\text{.}$$

If, for example, the number of function evaluations is the performance measure of interest, set $`t_{p,s}`$ accordingly.

We require a baseline for comparisons. We compare the performance on problem $`p`$ by solver $`s`$ with the best performance by any solver on this problem; that is, we use the performance ratio

$$r_{p,s}=\frac{t_{p,s}}{\mathrm{min}\{t_{p,s}:s𝒮\}}.$$

We assume that a parameter $`r_Mr_{p,s}`$ for all $`p,s`$ is chosen, and $`r_{p,s}=r_M`$ if and only if solver $`s`$ does not solve problem $`p`$. We will show that the choice of $`r_M`$ does not affect the performance evaluation.

The performance of solver $`s`$ on any given problem may be of interest, but we would like to obtain an overall assessment of the performance of the solver. If we define

$$\rho _s(\tau )=\frac{1}{n_p}\text{size}\{p𝒫:r_{p,s}\tau \},$$

then $`\rho _s(\tau )`$ is the probability for solver $`s𝒮`$ that a performance ratio $`r_{p,s}`$ is within a factor $`\tau `$ of the best possible ratio. The function $`\rho _s`$ is the (cumulative) distribution function for the performance ratio.

We use the term performance profile for the distribution function of a performance metric. Our claim is that a plot of the performance profile reveals all of the major performance characteristics. In particular, if the set of problems $`𝒫`$ is suitably large and representative of problems that are likely to occur in applications, then solvers with large probability $`\rho _s(\tau )`$ are to be preferred.

The term performance profile has also been used for a plot of some performance metric versus a problem parameter. For example, Higham (NJH96, pages 296–297) plots the ratio $`\gamma /A_1`$, where $`\gamma `$ is the estimate for the $`l_1`$ norm of a matrix $`A`$ produced by the LAPACK condition number estimator. Note that in Higham’s use of the term performance profile there is no attempt at determining a distribution function.

The performance profile $`\rho _s:[0,1]`$ for a solver is a nondecreasing, piecewise constant function, continuous from the right at each breakpoint. The value of $`\rho _s(1)`$ is the probability that the solver will win over the rest of the solvers. Thus, if we are interested only in the number of wins, we need only to compare the values of $`\rho _s(1)`$ for all of the solvers.

The definition of the performance profile for large values requires some care. We assume that $`r_{p,s}[1,r_M]`$ and that $`r_{p,s}=r_M`$ only when problem $`p`$ is not solved by solver $`s`$. As a result of this convention, $`\rho _s(r_M)=1`$, and

$$\rho _s^{}\underset{\tau r_M^{}}{lim}\rho _s(\tau )$$

is the probability that the solver solves a problem. Thus, if we are interested only in solvers with a high probability of success, then we need to compare the values of $`\rho _s^{}`$ for all solvers and choose the solvers with the largest value. The value of $`\rho _s^{}`$ can be readily seen in a performance profile because $`\rho _s`$ flatlines for large values of $`\tau `$; that is, $`\rho _s(\tau )=\rho _s^{}`$ for $`\tau [r_S,r_M)`$ for some $`r_S<r_M`$.

An important property of performance profiles is that they are insensitive to the results on a small number of problems. This claim is based on the observation that if $`\rho _s`$ and $`\widehat{\rho }_s`$ are defined, respectively, by the observed time sets $`t_{p,s}`$ and $`\widehat{t}_{p,s}`$, where

$$\widehat{t}_{p,s}=t_{p,s},p𝒫\{q\},$$

for some problem $`q𝒫`$, then $`\widehat{r}_{p,s}=r_{p,s}`$ for $`p𝒫\{q\}`$. Since only the ratio $`\widehat{r}_{q,s}`$ changes for any $`s𝒮`$,

$$|\rho _s(\tau )\widehat{\rho }_s(\tau )|\frac{1}{n_p},\tau ,$$

for $`s𝒮`$. Moreover, $`\widehat{\rho }_s(\tau )=\rho _s(\tau )`$ for $`\tau <\mathrm{min}\{r_{q,s},\widehat{r}_{q,s}\}`$ or $`\tau \mathrm{max}\{r_{q,s},\widehat{r}_{q,s}\}`$. Thus, if $`n_p`$ is reasonably large, then the result on a particular problem $`q`$ does not greatly affect the performance profiles $`\rho _s`$.

Not only are performance profiles relatively insensitive to changes in results on a small number of problems, they are also largely unaffected by small changes in results over many problems. We demonstrate this property by showing that small changes from $`r_{p,s}`$ to $`\widehat{r}_{p,s}`$ result in a correspondingly small $`L_1`$ error between $`\rho _s`$ and $`\widehat{\rho }_s`$.

###### Theorem 2.1

Let $`r_i`$ and $`\widehat{r}_i`$ for $`1in_p`$ be performance ratios for some solver, and let $`\rho `$ and $`\widehat{\rho }`$ be, respectively, the performance profiles defined by these ratios. If

$$|r_i\widehat{r}_i|ϵ,1in_p$$
(2.1)
for some $`ϵ>0`$, then

$$_1^{\mathrm{}}\left|\rho (t)\widehat{\rho }(t)\right|𝑑tϵ$$

###### Proof

Since performance profiles do not depend on the ordering of the data, we can assume that $`\{r_i\}`$ is monotonically increasing. We can reorder the sequence $`\{\widehat{r}_i\}`$ so that it is also monotonically increasing, and (2.1) still holds. These reorderings guarantee that $`\rho (t)=i/n_p`$ for $`t[r_i,r_{i+1})`$, with a similar result for $`\widehat{\rho }`$. We now show that for any integer $`k`$ with $`1kn_p`$,

$$_1^{s_k}\left|\rho (t)\widehat{\rho }(t)\right|𝑑tk\left(\frac{ϵ}{n_p}\right),$$
(2.2)
where $`s_k=\mathrm{max}(r_k,\widehat{r}_k)`$, and

$$|r_i\widehat{r}_i|ϵ,1ik,\text{and}\widehat{r}_i=r_i,k<in_p.$$
(2.3)
The proof is completed when $`k=n_p`$.

The case $`k=1`$ follows directly from the definition of a performance profile, so assume that (2.2) holds for any performance data such that (2.3) holds. We now prove that (2.2) holds for $`k+1`$ by proving that

$$_{s_k}^{s_{k+1}}\left|\rho (t)\widehat{\rho }(t)\right|𝑑t\left(\frac{ϵ}{n_p}\right),k<n_p.$$
(2.4)
Together, (2.2) and (2.4) show that (2.2) holds for $`k+1`$.

We present the proof for the case when $`\widehat{r}_kr_k`$. A similar argument can be made for $`r_k\widehat{r}_k`$. If $`\widehat{r}_kr_k`$ then $`s_k=r_k`$ and $`\widehat{r}_kr_kr_{k+1}`$. The argument depends on the position of $`\widehat{r}_{k+1}`$ and makes repeated use of the fact that $`\rho (t)=k/n_p`$ for $`t[r_k,r_{k+1})`$, with a similar result for $`\widehat{\rho }`$.

If $`r_{k+1}\widehat{r}_{k+1}`$ then $`\rho (t)=\widehat{\rho }(t)`$ in $`[r_k,r_{k+1})`$. Also note that $`\left|\rho (t)\widehat{\rho }(t)\right|=1/n_p`$ in $`[r_{k+1},\widehat{r}_{k+1})`$. Hence, (2.4) holds with $`s_{k+1}=\widehat{r}_{k+1}`$.

The case where $`\widehat{r}_{k+1}r_{k+1}`$ makes use of the observation that $`\widehat{r}_i=r_ir_{k+1}`$ for $`i>k+1`$. If $`r_k\widehat{r}_{k+1}r_{k+1}`$, then $`\rho (t)=\widehat{\rho }(t)`$ in $`[r_k,\widehat{r}_{k+1})`$, and $`\left|\rho (t)\widehat{\rho }(t)\right|=1/n_p`$ in $`[\widehat{r}_{k+1},r_{k+1})`$. Hence, (2.4) holds. On the other hand, if $`\widehat{r}_k\widehat{r}_{k+1}r_k`$, then we only need to note that $`\left|\rho (t)\widehat{\rho }(t)\right|=1/n_p`$ in $`[r_k,r_{k+1})`$ in order to conclude that (2.4) holds.

We have shown that (2.2) holds for all integers $`k`$ with $`1kn_p`$. In particular, the case $`k=n_p`$ yields our result since $`\rho (t)=\widehat{\rho }(t)`$ for $`t[s_{n_p},\mathrm{})`$. ∎

## 3 Benchmarking Data

The timing data used to compute the performance profiles in Sections 4 and 5 are generated with the COPS test set, which currently consists of seventeen different applications, all models in the AMPL AMPL modeling language. The choice of the test problem set $`𝒫`$ is always a source of disagreement because there is no consensus on how to choose problems. The COPS problems are selected to be interesting and difficult, but these criteria are subjective. For each of the applications in the COPS set we use four instances of the application obtained by varying a parameter in the application, for example, the number of grid points in a discretization. Application descriptions and complete absolute timing results for the full test set are given in EDD00.

Section 4 focuses on only the subset of the eleven optimal control and parameter estimation applications in the COPS set, while the discussion in Section 5 covers the complete performance results. Accordingly, we provide here information specific to this subset of the COPS problems as well as an analysis of the test set as a whole. Table 3.1 gives the quartiles for three problem parameters: the number of variables $`n`$, the number of constraints, and the ratio $`(nn_e)/n`$, where $`n_e`$ is the number of equality constraints. In the optimization literature, $`nn_e`$ is often called the degrees of freedom of the problem, since it is an upper bound on the number of variables that are free at the solution.

The data in Table 3.1 is fairly representative of the distribution of these parameters throughout the test set and shows that at least three-fourths of the problems have the number of variables $`n`$ in the interval $`[400,5000]`$. Our aim was to avoid problems where $`n`$ was in the range $`[1,50]`$ because other benchmarking problem sets tend to have a preponderance of problems with $`n`$ in this range. The main difference between the full COPS set and the COPS subset is that the COPS subset is more constrained with $`n_en/2`$ for all the problems. Another feature of the COPS subset is that the equality constraints are the result of either difference or collocation approximations to differential equations.

We ran our final complete runs with the same options for all models. The options involve setting the output level for each solver so that we can gather the data we need, increasing the iteration limits as much as allowed, and increasing the super-basics limits for MINOS and SNOPT to 5000. None of the failures we record in the final trials include any solver error messages about having violated these limits. While we relieved restrictions unnecessary for our testing, all other parameters were left at the solvers’ default settings.

The script for generating the timing data sends a problem to each solver successively, so as to minimize the effect of fluctuation in the machine load. The script tracks the wall-clock time from the start of the solve, killing any process that runs 3,600 seconds, which we declare unsuccessful, and beginning the next solve. We cycle through all the problems, recording the wall-clock time as well as the combination of AMPL system time (to interpret the model and compute varying amounts of derivative information required by each solver) and AMPL solver time for each model variation. We repeat the cycle for any model for which one of the solvers’ AMPL time and the wall-clock time differ by more than 10 percent. To further ensure consistency, we have verified that the AMPL time results we present could be reproduced to within 10 percent accuracy. All computations were done on a SparcULTRA2 running Solaris 7.

We have ignored the effects of the stopping criteria and the memory requirements of the solvers. Ideally we would have used the same stopping criteria, but this is not possible in the AMPL environment. In any case, differences in computing time due to the stopping criteria are not likely to change computing times by more than a factor of two. Memory requirements can also play an important role. In particular, solvers that use direct linear equation solvers are often more efficient in terms of computing time provided there is enough memory.

The solvers that we benchmark have different requirements. MINOS and SNOPT use only first-order information, while LANCELOT and LOQO need second-order information. The use of second-order information can reduce the number of iterations, but the cost per iteration usually increases. In addition, obtaining second-order information is more costly and may not even be possible. MINOS and SNOPT are specifically designed for problems with a modest number of degrees of freedom, while this is not the case for LANCELOT and LOQO. As an example of comparing solvers with similar requirements, Section 6 shows the performance of linear programming solvers.

## 4 Case Study: Optimal Control and Parameter Estimation Problems

We now examine the performance of LANCELOT ARC92, MINOS BAM95, SNOPT PEG97, and LOQO RJV97b on the subset of the optimal control and parameter estimation problems in the COPS cops-home test set. Figures 4.1 and 4.2 show the performance profiles in different ranges to reflect various areas of interest. Our purpose is to show how performance profiles provide objective information for analysis of a large test set.

Figure 4.1 shows the performance profiles of the four solvers for small values of $`\tau `$. By showing the ratios of solver times, we eliminate any weight of importance that taking straight time differences might give to the problems that require a long run time of every solver. We find no need to eliminate any test problems from discussion. For this reason, solvers receive their due credit for completing problems for which one or more of the other solvers fails. In particular, $`1\rho _s(\tau )`$ is the fraction of problems that the solver cannot solve within a factor $`\tau `$ of the best solver, including problems for which the solver in question fails.

From this figure it is clear that LOQO has the most wins (has the highest probability of being the optimal solver) and that the probability that LOQO is the winner on a given problem is about $`.61`$. If we choose being within a factor of $`4`$ of the best solver as the scope of our interest, then either LOQO or MINOS would suffice; but the performance profile shows that the probability that these two solvers can solve a job within a factor $`4`$ of the best solver is only about $`70\%`$. SNOPT has a lower number of wins than either LOQO or MINOS, but its performance becomes much more competitive if we extend our $`\tau `$ of interest to $`7`$.

Figure 4.2 shows the performance profiles for all the solvers in the interval $`[1,100]`$. If we are interested in the solver that can solve $`75\%`$ of the problems with the greatest efficiency, then MINOS stands out. If we hold to more stringent probabilities of completing a solve successfully, then SNOPT captures our attention with its ability to solve over $`90\%`$ of this COPS subset, as displayed by the height of its performance profile for $`\tau >40`$. This graph displays the potential for large discrepancies in the performance ratios on a substantial percentage of the problems. Another point of interest is that LOQO, MINOS, and SNOPT each have the best probability $`\rho _s(\tau )`$ for $`\tau `$ in some interval, with similar performance in the interval $`[15,40]`$.

An observation that emerges from these figures is the lack of consistency in quartile values of time ratios. The top three solvers share a minimum ratio of $`1`$, and LOQO and MINOS also share first quartile values of $`1`$. In other words, these two solvers are the best solvers on at least 25% of the problems. LOQO bests MINOS’s median value with $`1`$ compared with $`2.4`$, but MINOS comes back with a third quartile ratio of $`4.3`$ versus $`13.9`$ for LOQO, with SNOPT mixing results further by also beating LOQO with $`12.6`$. By looking at Figures 4.1 and 4.2, we see that progress between quartiles does not necessarily proceed linearly; hence, we really lose information if we do not provide the full data. Also, the maximum ratio would be $`r_M`$ for our testing, and no obvious alternative value exists. As an alternative to providing only quartile values, however, the performance profile yields much more information about a solver’s strengths and weaknesses.

We have seen that at least two graphs may be needed to examine the performance of the solvers. Even extending $`\tau `$ to $`100`$, we fail to capture the complete performance data for LANCELOT and LOQO. As a final option, we display a log scale of the performance profiles. In this way, we can show all activity that takes place with $`\tau <r_M`$ and grasp the full implications of our test data regarding the solvers’ probability of successfully handling a problem. Since we are also interested in the behavior for $`\tau `$ close to unity, we use a base of 2 for the scale. In other words, we plot

$$\tau \frac{1}{n_p}\text{size}\{p𝒫:\mathrm{log}_2\left(r_{p,s}\right)\tau \}$$

in Figure 4.3. This graph reveals all the features of the previous two graphs and thus captures the performance of all the solvers. The disadvantage is that the interpretation of the graph is not as intuitive, since we are using a log scale.

Figures 4.1 and 4.2 are mapped into a new scale to reflect all data, requiring at least the interval $`[0,\mathrm{log}_2(1043)]`$ in Figure 4.3 to include the largest $`r_{p,s}<r_M`$. We extend the range slightly to show the flatlining of all solvers. The new figure contains all the information of the other two figures and, in addition, shows that each of the solvers fails on at least $`8\%`$ of the problems. This is not an unreasonable performance for the COPS test set because these problems were generally chosen to be difficult.

## 5 Case Study: The Full COPS

We now expand our analysis of the data to include all the problems in version 2.0 of the COPS cops-home test set. We present in Figure 5.1 a log<sub>2</sub> scaled view of the performance profiles for the solvers on that test set.

Figure 5.1 gives a clear indication of the relative performance of each solver. As in the performance profiles in Section 4, this figure shows that performance profiles eliminate the undue influence of a small number of problems on the benchmarking process and the sensitivity of the results associated with the ranking of solvers. In addition, performance profiles provide an estimate of the expected performance difference between solvers.

The most significant aspect of Figure 5.1, as compared with Figure 4.3, is that on this test set LOQO dominates all other solvers: the performance profile for LOQO lies above all others for all performance ratios. The interpretation of the results in Figure 5.1 is important. In particular, these results do not imply that LOQO is faster on every problem. They indicate only that, for any $`\tau 1`$, LOQO solves more problems within a factor of $`\tau `$ of any other solver time. Moreover, by examining $`\rho _s(1)`$ and $`\rho _s(r_M)`$, we can also say that LOQO is the fastest solver on approximately 58% of the problems, and that LOQO solves the most problems (about 87%) to optimality.

The difference between the results in Section 4 and these results is due to a number of factors. First of all, as can be seen in Table 3.1, the degrees of freedom for the full COPS test set is much larger than for the restricted subset of optimal control and parameter estimation problems. Since, as noted in Section 3, MINOS and SNOPT are designed for problems with a modest number of degrees of freedom, we should expect the performance of MINOS and SNOPT to deteriorate on the full COPS set. This deterioration can be seen by comparing Figure 5.1 with Figure 4.3 and noting that the performance profiles of MINOS and SNOPT are similar but generally lower in Figure 5.1.

Another reason for the difference between the results in Section 4 and these results is that MINOS and SNOPT use only first-order information, while LOQO uses second-order information. The benefit of using second-order information usually increases as the number of variables increases, so this is another factor that benefits LOQO.

The results in this section underscore our observation that performance profiles provide a convenient tool for comparing and evaluating the performance of optimization solvers, but, like all tools, performance profiles must be used with care. A performance profile reflects the performance only on the data being used, and thus it is important to understand the test set and the solvers used in the benchmark.

## 6 Case Study: Linear Programming

Performance profiles can be used to compare the performance of two solvers, but performance profiles are most useful in comparing several solvers. Because large amounts of data are generated in these situations, trends in performance are often difficult to see. As a case study, we use data obtained by Mittelmann Benchmarks. Figure 6.1 shows a plot of the performance profile for the time ratios in the data Benchmark of LP solvers on a Linux-PC (5-25-2000), which includes results for COPL\_LP (1.0), PCx (1.1), SOPLEX (1.1), LPABO (5.6), MOSEK (1.0b), BPMPD (2.11), and BPMPD (2.14).

In keeping with our graphing practices with the COPS set, we designate as failures those solves that are marked in the original table as stopping close to the final solution without convergence under the solver’s stopping criteria. One feature we see in the graph of Mittelmann’s results that does not appear in the COPS graphs is the visual display of solvers that never flatline. In other words, the solvers that climb off the graph are those that solve all of the test problems successfully. As with Figure 4.3, all of the events in the data fit into this log-scaled representation. While this data set cannot be universally representative of benchmarking results by any means, it does show that our reporting technique is applicable beyond our own results.

As in the case studies in Sections 4 and 5, the results in Figure 6.1 give an indication of the performance of LP solvers only on the data set used to generate these results. In particular, the test set used to generate Figure 6.1 includes only problems selected by Mittelmann for his benchmark. The advantage of these results is that, unlike the solvers in Sections 4 and 5, all solvers in Figure 6.1 have the same requirements.

## 7 Conclusions

We have shown that performance profiles combine the best features of other tools for benchmarking and comparing optimization solvers. Clearly, the use of performance profiles is not restricted to optimization solvers and can be used to compare solvers in other areas.

We have not addressed the issue of how to select a collection of test problems to justify performance claims. Instead, we have provided a tool – performance profiles – for evaluating the performance of two or more solvers on a given set of test problems. If the data is obtained by following careful guidelines Crowder:1979:RCE; JacBNP91, then performance profiles can be used to justify performance claims. We emphasize that claims on the relative performance of the solvers on problems not in the test set should be made with care.

The Perl script perf.pl on the COPS site cops-home generates performance profiles formatted as Matlab commands to produce a composite graph as in Figures 4.1 and 4.2. The script accepts data for several solvers and plots the performance profile on an interval calculated to show the full area of activity. The area displayed and scale of the graph can then be adjusted within Matlab to reflect particular benchmarking interests.

###### Acknowledgements.

This work was supported by the Mathematical, Information, and Computational Sciences Division subprogram of the Office of Advanced Scientific Computing, U.S. Department of Energy, under Contract W-31-109-Eng-38, and by the National Science Foundation (Challenges in Computational Science) grant CDA-9726385 and (Information Technology Research) grant CCR-0082807. Version 1.0 of the COPS problems was developed by Alexander Bondarenko, David Bortz, Liz Dolan, and Michael Merritt. Their contributions were essential because, in many cases, version 2.0 of the problems are closely related to the original version. Alex Bondarenko, Nick Gould, Sven Leyffer and Bob Vanderbei contributed interesting and spirited discussions on problem formulation, while Bob Fourer and David Gay generously shared their AMPL expertise with us. Hans Mittelmann deserves special note for paving the way for COPS with his benchmarking work. We also thank the referees and the associate editor, Jorge Nocedal, for their comments. One of the referees provided a thoughtful report that led to additional details on the sensitivity of performance profiles.