{'fraction_non_alphanumeric': 0.033280440955403376, 'fraction_numerical': 0.02341183675741885, 'mean_word_length': 4.556028149408398, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Massive star formation and evolution in starburst galaxies: mid-infrared spectroscopy with ISO-SWS

We present new ISO-SWS data for a sample of 27 starburst galaxies, and with these data examine the issues of formation and evolution of the most massive stars in starburst galaxies. Using starburst models which incorporate time evolution, new stellar atmosphere models for massive stars, and a starburst model geometry derived from observations of the prototypical starburst M82, we model the integrated mid-infrared line ratio [NeIII](15.6 microns)/[NeII](12.8 microns). This line ratio is sensitive to the hardness of the stellar energy distribution and therefore to the most massive stars present. We conclude from our models, with consideration of recent determinations of the stellar census in local, high-mass star forming regions, that the [NeIII]/[NeII] ratios we measure are consistent with the formation of massive (~50-100 solar mass) stars in most starbursts. In this framework, the low nebular excitation inferred from the measured line ratios can be attributed to aging effects. By including estimates of the ratio of infrared-to-Lyman continuum luminosity for the galaxies in our sample, we further find that most starbursts are relatively short-lived (1-10 million years), only a few O-star lifetimes. We discuss a possible cause of such short events: the effectiveness of stellar winds and supernovae in destroying the starburst environment.


Introduction
Starburst galaxies constitute an important class of extragalactic objects. They contribute a significant fraction of the total high-mass star formation in the local universe (e.g. Soifer et al. 1987;Gallego et al. 1995), and about 25% of the highmass star formation within 10 Mpc occurs in four starburst galaxies (Heckman 1998). At intermediate and high redshifts, important populations of galaxies are observed which exhibit properties in-NASA. The SWS is a joint project of SRON and MPE. 1 dicating intense star formation activity (e.g. Colless et al. 1994;Babul & Ferguson 1996;Steidel et al. 1996;Lowenthal et al. 1997). A strong correlation between the most prodigious starbursts and interaction or merger events (e.g. Condon et al. 1982;Kennicutt et al. 1987;Telesco, Wolstencroft & Done 1988;Joseph 1990) further emphasizes the importance of starbursts in galaxy evolution.
Despite the obvious importance of starbursts, a canvas of recent starburst galaxy studies suggests that a full "prescription" for these events has yet to be written. Discrepancies exist even for the most massive stars formed, which represent the largest energy contributions. For instance, it has been suggested by some authors that the production of the highest mass stars may be suppressed in starburst galaxies, based on low measured values of diagnostic line ratios such as He I (2.06µm)/ HI(Brγ) and [Ne II]/[Ar III] detected in well-known starburst sources. These results have been interpreted as possibly indicating severe upper mass cutoffs to the initial mass function (IMF), some as low as ∼25-30M ⊙ (e.g. Puxley et al. 1989;Doyon et al. 1994;Doherty et al. 1995;Achtermann & Lacy 1995;Beck, Kelly, & Lacy 1997).
A lack of massive stars in starbursts is difficult to reconcile with a growing body of evidence for the formation of very massive stars in nearby regions of active star-formation. In recent studies of local high-mass star-forming regions at both high and low metallicity, stars up to at least 100 M ⊙ are observed; these regions include the Galactic Center (e.g. Krabbe et al. 1995;Najarro et al. 1997;Serabyn, Shupe, & Figer 1998;Figer et al. 1998), the Galactic star-forming region NGC 3603 (Drissen et al. 1995;Eisenhauer et al. 1998), and the R136 cluster at the center of 30 Doradus (e.g. Hunter et al. 1995, Massey & Hunter 1998. There is also a growing body of indirect evidence for the presence of very massive stars in starbursts. HST optical imaging has revealed the presence of young compact "super star clusters" in several nearby starburst galaxies, including M 82 (O'Connell et al. 1995), NGC 4038/4039 (Whitmore & Schweizer 1995, NGC 5253 (Meurer et al. 1995), He 2-10 (Conti & Vacca 1994), NGC 1569and NGC 1705Ho & Filippenko 1996;Sternberg 1998), and NGC 1140 . These super star clusters have extreme luminosities, and in some cases bright emission features characteristic of Wolf-Rayet stars (e.g., He II 4696Å emission); such properties are difficult to explain without a large contribution from very massive stars (M≥60M ⊙ ; see e.g., Gonzalez Delgado et al. 1997;Schaerer et al. 1997).
Several authors have recently argued that observed low-excitation nebular line ratios in starbursts may be due to aging effects rather than a severe upper mass cutoff (e.g. Rieke et al. 1993;Genzel, Hollenbach, & Townes 1994;Achtermann & Lacy 1995;Satyapal et al. 1997;Engelbracht et al. 1998; see also Vanzi & Rieke 1997). Studying a sufficiently large sample of starbursts is therefore valuable in breaking the degeneracy between age and upper mass cutoff. As we are not likely to catch all starbursts at a "late" age, and the ultraviolet, optical, and far-infrared luminosities decline with time after the peak star formation activity, we can determine the relative importance of variations in star formation parameters and aging effects by studying a larger sample.
In this work, we present results of a starburst modeling program to interpret mid-infrared (MIR) spectroscopy from the Short Wavelength Spectrometer (SWS; de Graauw et al. 1996) aboard the Infrared Space Observatory (ISO; Kessler et al. 1996). To this end, we have gathered atomic fine structure line fluxes (in particular for the [Ne II] 12.8µm and [Ne III] 15.6µm lines) for 27 starburst galaxies. This project is ideally suited for investigating the properties of the high-mass stellar population in heavily obscured star-forming regions of starburst galaxies, since the extinction at MIR wavelengths is only a few percent of the optical extinction. This investigation thus constitutes an important contribution to the understanding of star formation in starburst systems, a field that has been dominated by studies in the optical and ultraviolet regimes in the intervening years between the IRAS and ISO missions. As recent studies of the cosmic infrared background suggest a significant contribution from dusty starburst systems (e.g., Puget et al. 1996;Hughes et al. 1998;Hauser et al. 1998;Elbaz et al. 1998), the understanding of dusty starbursts is important in defining the picture of the most powerful star-forming events in the universe.
We present the sample data in §2. We outline our modeling procedure and compare our ISO-SWS observations with model predictions in §3.
We examine additional constraints on this starburst modeling study and examine the robustness of the modeling analysis in §4. We discuss the implications of our findings and summarize this work in §5.

Database selection
The database for this project was created by combining data from the MPE guaranteed-time program on bright galactic nuclei (15 objects) and an open-time observing proposal (12 objects), resulting in an ISO-SWS spectroscopic dataset for a total sample of 27 galaxies. The objects from the guaranteed-time program are all well-studied, infrared-bright starburst galaxies, and the opentime sample objects were selected on the basis of their 60µm flux densities (S 60µm >20 Jy) and their observability by ISO in the last few months of the mission.
Within these programs, we have targeted observations at a specific set of spectral lines, particularly the [Ne II] 12.8µm and [Ne III] 15.6µm fine structure lines; in this work we will concentrate on the ratio of these two lines as a tracer of the massive star population, with additional lines introduced to constrain metallicities and gas densities. The neutral Ne atom and the Ne + ion have ionization potentials of 21.56 eV and 40.95 eV respectively, making the [Ne III]/[Ne II] line ratio very sensitive to the spectral shape of the UV radiation field and thus to the properties of the most massive stars present. This information can be used to provide constraints on various star formation parameters, including the upper mass cutoff of the IMF and the age and duration of the starburst (see Kunze et al. 1996 andRigopoulou et al. 1996 for earlier models using this and other MIR line ratios).
We have chosen the neon line ratio as the focus of this larger study, both for its sensitivity to massive, hot stars as well as the minimal uncertainties introduced in constructing the ratio of these two lines from ISO-SWS observations.  Kunze et al. 1996;Rigopoulou et al. 1996). In addition, the [S IV]/[Ne II] ratio has inherent uncertainties due to the determination of metallicities, and all of the alternatives listed above are subject to corrections for aperture size variations and/or systematic uncertainties due to measurements in different ISO-SWS bandpasses. Systematic effects are minimized for the neon lines because they were observed with the same ISO aperture and bandpass. In addition, because the wavelengths of the two MIR neon lines are close, extinction effects can be neglected. For example, for a dust distribution in M82 which is well mixed with the stars and obeys a standard extinction law (e.g. Draine 1989) with A V ∼50 mag (McLeod et al. 1993;Schreiber 1999), the observed line ratio is higher than the intrinsic line ratio by only a few percent. The same is true in the ultraluminous galaxy Arp 220 for a uniform foreground dust screen providing A V ∼50 mag for the starburst emission (Sturm et al. 1996) lines lie in the same spectral region as a broad silicate absorption feature at 9.8µm, further complicating the measurement of line fluxes.

Observations and Data Reduction
The majority of the data described herein was taken in grating line scan mode SWS02, with data from some galaxies taken from full-grating-scan SWS01 observations. Data reduction was completed using the ISO SWS Interactive Analysis (IA) package, which includes interactive tools for dark current subtraction, flat fielding, and removal of fringing caused by interference within the detectors themselves (see e.g., Schaeidt et al. 1996). The reduction completed for this paper used post-mission calibration files (June 1998).
Line fluxes were determined by integration over the observed line profiles, which were generally close to Gaussian in shape. The systematic calibration of line fluxes has an uncertainty of ≤30% (Schaeidt et al. 1996), though the effect of this uncertainty is minimized by comparing two lines from the same bandpass. The uncertainties in the measured fluxes range from ≤10% for strong line sources to 30-40% for weaker sources (generally from the open-time proposal). For comparison with models ( §4), we will assume uncertainties in the line ratios for the entire sample to be ±50%. The aperture size for the bandpass containing the [Ne II] and [Ne III] lines is 14 ′′ x 27 ′′ . Because of the weak sensitivity of this line ratio to extinction, no corrections for extinction were made (though see §4.3 for other effects of dust). Table 1 lists the MIR neon line strengths for the galaxy sample. Figure 1 shows the corresponding neon line spectra, and Figure 2 shows the range of observed [Ne III]/[Ne II] ratios as a function of the infrared luminosity. The infrared luminosity L IR is computed from the IRAS flux densities following the prescription of Sanders & Mirabel (1996, their Table 1): L IR = L 8−1000 µm = 5.6 × 10 5 D 2 Mpc (S 100 + 2.58S 60 + 5.16S 25 + 13.48S 12 ), where S λ is the flux density in Jy in the band centered at wavelength λ. For several of the galaxies, the sizes of the emission regions exceed the SWS field of view, so we have applied a aperture-size correction to L IR based on the spatial distribution of MIR or radio continuum emission, millimetric CO emission, or near-infrared hydrogen recombination line emission from the literature, assuming they all trace the same emission regions as the MIR neon lines. For eight galaxies, the observed L IR was multiplied by a factor less than 0.5 to correct for the smaller SWS field of view, and for two galaxies there was insufficient data to determine an aperture-size correction. For the remaining galaxies, 11 required no correction, and seven galaxies were corrected by factors between 0.5 and 1.0.
In general, all of the sample objects have neon ratios between 0.05 and 1.0; however, the lowmetallicity systems NGC 5253 and IIZw40 have neon ratios of 3.5 and 12, respectively. In §4.3, we discuss the effects of metallicity on model line ratios, in order to address the observed range of emission properties for the entire sample. The sample presented in this paper includes nine objects for which [Ne II] line fluxes were measured by Roche et al. (1991), who took observations from the ground in the 8-13µm atmospheric window. It is difficult to compare the fluxes derived from these two sets of observations due to differences in aperture size and the uncertain contribution of the broad 12.7µm PAH feature at the end of the 8-13µm band. This is particularly true for galaxies such as M82 which exhibit strong PAH features (e.g., Roche et al. 1991;Sturm et al. 2000); in the case of M82, we estimate from the SWS data that the broad 12.7µm feature could contribute a factor of 50% more flux to the apparent [Ne II] line at the resolution observed by Roche et al. . In addition, for some objects in the Roche et al. sample the [Ne II] line is redshifted out of the 8-13µm wavelength range and so cannot be measured from the ground. However, for the 7 objects for which comparisons could be made, the line fluxes are consistent to within these rather large uncertainties. NGC 5253 and IIZw40 were also singled out by Roche et al. (1991) as displaying high-excitation spectra, though the [Ne II] line was not detected. Roche et al. noted that these galaxies had little emission from broad PAH features, which is consistent with high excitation and UV energy density in a low-metallicity environment (see also Vigroux 1997;Thuan et al. 1999;and Madden 2000).
For further comparison, we include in Figure 2 the neon ratios for similar observations of W51, 30 Doradus, and the central parsec of the Galaxy. The neon ratios for W51 and 30 Doradus are higher than for the starburst galaxy sample, while that for the Galactic Center is lower (a further discussion of the Galactic Center can be found in §4.2). Our ISO-SWS data thus generally confirm findings from NIR and MIR spectroscopy that powerful dusty starbursts exhibit low nebular excitation in comparison to Galactic HII regions (e.g., Doherty et al. 1994;Cox et al. 1999). However, we find a range of values, with the most excited regions having neon ratios similar to those of young, Galactic HII regions.

Nebular emission from young stellar clusters
The starburst modeling method presented here predicts MIR spectral line ratios for evolving stellar clusters. Observations of nearby starburst galaxies reveal that starburst regions are in fact composed of many individual units. In addition to the HST detections of super star clusters in starburst galaxies, bright compact sources seen in high-resolution near-infrared broad-band images have been interpreted as young clusters of red supergiants (e.g. Tacconi-Garman, Sternberg & Eckart 1996;Satyapal et al. 1997). Furthermore, important substructure in the ISM on 10-parsecor even parsec-scales is observed or inferred from modeling (e.g. Carral et al. 1994;Shen & Lo 1995;Achtermann & Lacy 1995).
In §3.4, we present a set of models in which the starburst region is composed of a heterogeneous ensemble of recently-formed star clusters. Therefore, we first lay the groundwork of modeling the emission from the HII region around stellar clusters.

The SED produced by a stellar cluster
We have used the evolutionary synthesis code STARS (Sternberg 1998) to model the integrated properties of individual stellar clusters. This code is based on the most recent Geneva stellar evolutionary tracks (Schaller et al. 1992;Schaerer et al. 1993aSchaerer et al. , 1993bCharbonnel et al. 1993;Meynet et al. 1994), and is similar to other models developed by various investigators (e.g. Tinsley 1972;Huchra 1977;Bruzual 1983;Guiderdoni & Rocca-Volmerange 1987;Mas-Hesse & Kunth 1991;Bruzual & Charlot 1993;Leitherer et al. 1999). STARS follows the evolution in the H-R diagram (or, alternatively, in a log T eff −log g diagram) of a stellar population whose composition is specified by a timeindependent IMF. The star formation rate (SFR) is assumed to decline exponentially as e −t b /tsc , where t b is the age of the cluster and t sc is the burst timescale. At a given age, the integrated SED is obtained by summing over the contributions of all stars present.
The SED assigned to each location in the log T eff − log g diagram is taken from a hybrid grid created from two separate libraries, in order to best represent the properties of all stars. For very hot stars, the effects of line blanketing and blocking in rapidly expanding atmospheres are not well represented by LTE models; thus for stars with temperatures above 25,000 K (M ZAMS 13M ⊙ ), we have used new, non-LTE SEDs by Pauldrach et al. (1998). These models include the effects of stellar winds and mass loss, and cover effective temperatures up to 60,000 K (M ZAMS ∼ 120M ⊙ ). The "cornerstone" grid of Pauldrach models con-sists of 8 main sequence models and 6 supergiant models, from which SEDs for intermediate values of Teff and log g are interpolated. Table 2 lists the (T ef f , log g) combinations in the cornerstone grid, and Figure 3 shows the sampling of the interpolated grid of solar-metallicity models. For stars with effective temperatures below 19,000 K (M ZAMS < 8M ⊙ ), we use the Kurucz (1992) library of LTE SEDs. "Hybrid" SEDs for stars with intermediate effective temperatures were created by interpolating between the 19,000 K Kurucz model and the 25,000 K Pauldrach model. Figure  4 compares examples of Kurucz, Pauldrach, and hybrid SEDs over a range of temperatures.
We have created SED grids for two metallicities, solar (Z ⊙ ) and 0.2Z ⊙ . The low-metallicity grid allows us to better represent the stellar populations in NGC 5253 (∼0.24Z ⊙ , Hunter, Gallagher, & Rautenkranz 1982) and IIZw40 (∼0.15Z ⊙ , Masegosa, Moles, & Campos-Aguilar 1994), and provides a benchmark for the effects of metallicity variations on our models. The low-metallicity grid differs from the solar-metallicity grid only in the sampling of the cooler Kurucz models. Starburst models were computed using each grid. Table 3 summarizes the model parameters used in this work. For each metallicity grid, only the upper mass cutoff (m up ), burst age (t b ), and burst timescale (t sc ) were allowed to vary. The IMF was taken to have a Salpeter power-law index, dN/dm ∝ m −2.35 (Salpeter 1955), between a lower mass cutoff (m low ) fixed at 1 M ⊙ and a variable m up . Though some studies suggest that the IMF in starburst galaxies is deficient in stars below a few solar masses (e.g. Rieke et al. 1980Rieke et al. , 1993Shier, Rieke & Rieke 1996;Engelbracht et al. 1998), variations in the low-mass IMF do not affect the ionizing continuum and the results presented here; thus we will not consider variations in m low .

Modeling the nebular emission around clusters
To model the nebular emission excited by a stellar cluster, we have used the photoionization code CLOUDY (version C90.04, Ferland 1996), with the source of ionizing radiation described by SEDs from STARS and the ionized nebula represented by a shell around a central stellar cluster. The nebular emission from young clusters will depend on the properties of the surrounding ISM. The nebular conditions are specified by the gas and dust composition, the hydrogen number density n H of the gas, the distance R between the ionizing source and the illuminated surface of the cloud, and the ionization parameter U . To estimate the appropriate inputs for these models, we have collected measurements from our data and from the literature which can be used to infer metallicities, n H , R, and U for as many sample objects as possible.
Using our neon line measurements and published hydrogen recombination line measurements, we have calculated the abundance of neon relative to hydrogen for the sample using the strong-line method. This method provides abundances good to within factors of 2-3 (Osterbrock 1989), sufficient for our purposes here. Excluding NGC 5253 and IIZw40, the calculated metallicities are consistent with solar metallicity, to within the uncertainties: the average metallicity for the 13 objects for which metallicities could be determined was 1.9±1 Z ⊙ . The metallicities of these objects will be more rigorously derived by Spoon et al. (in prep).
We have computed the photoionization models assuming a solar gas-phase abundance for the majority of the sample, as well as models assuming a 0.2 Z ⊙ gas-phase abundance for NGC 5253 and IIZw40. As changes in gas-phase abundances have relatively little effect on the neon ratio, we have not attempted to further tailor the gas-phase abundances for each galaxy. The comparison between solar and sub-solar metallicity models will be discussed further in §4.3.
Due to uncertainties about the dust properties in the starburst environments we have surveyed, we neglect the effects of dust grains mixed with the ionized gas. The gas density is assumed to be uniform, and we take n H = n e , where n e is the electron density. The electron density was determined from the [S III] 18.7µm/33.5µm line ratio for the galaxies in the guaranteed-time sample, with the line fluxes corrected for the extinction determined in §4.1, assuming the extinction law recommended by Draine (1989). The typical uncertainties on the line ratios are 50%, and are dominated by extinction and aperture-size corrections (the [S III] 18.7µm and 33.5µm lines were measured through 14 ′′ × 27 ′′ and 20 ′′ × 33 ′′ aper-tures, respectively), and the absolute flux calibration. Within the uncertainties, the dereddened line ratios for this subsample lie in the low-density limit, and imply n e 10 3 cm −3 . We adopt a value of 300 cm −3 for all models, inferred from the average of the line ratios for the galaxies in the guaranteed-time sample.
The ionization parameter (U) of a nebula is one of the most important parameters in photoionization models. A measure of U gives the number of ionizing photons impinging at the surface of the nebula per hydrogen atom: where Q Lyc is the rate of production of ionizing photons, R is the radius of the shell, and c is the speed of light. The determination of the ionization parameter U from observed properties depends on the assumptions made about the geometry of the ionizing clusters and of the nebulae. As discussed above, entire starburst regions such as included in the SWS field of view likely do not consist of a gas shell illuminated by a single, centrallyconcentrated cluster. If the ISM and a number of young star clusters are in a mixed distribution, gas clouds will shield each other partially from ionizing radiation from clusters distributed throughout a starburst region of radius R (e.g., Wolfire, Tielens, & Hollenbach 1990). Both representative geometries are shown schematically in Figure 5. In the distributed cluster geometry, we can describe the conditions at the illuminated face of each cloud structure in terms of an effective ionization parameter, U eff which may be written The distributed cluster geometry is analogous to that in "mixed extinction" models, where R/λ is a measure of the optical depth and λ is the mean free path of photons. The parameter λ depends on the number density and sizes of the clouds which reside in the starburst region. Therefore, the determination of U eff , and the corresponding effective radius R eff (from replacing U with U eff in Eqn. 1), requires a knowledge of the relative distributions of gas clouds and ionizing clusters, as well as the properties of the clouds and clusters themselves. Such information is available for the archetypal starburst galaxy M 82: by combining new near-and mid-infrared data with the results of modeling of photodissociation regions by Lord et al. (1996) and observations of the molecular gas by Shen & Lo (1995), Schreiber (1999) reassessed the geometry of clusters and gas clouds in the starburst region of M82. The effective values, from an investigation of the ionization parameter in regions of varying size within the starburst core of M 82, are log U eff = −2.3 and R eff = 25 pc (Schreiber 1999).
The investigation of M 82 by Schreiber (1999) further shows that the effective ionization parameter is essentially independent of the location and the size of the regions studied, from individual burst sites ≈ 20 pc in radius to the entire starburst core extending over 450 pc. This strongly suggests that the regions of intense starburst activity in M 82, which dominate the nebular line emission, have very similar properties over a wide range of physical sizes. On average, the nebular conditions of the gas in the starburst regions of M82 can thus be described locally and globally with the same parameters.
A similar analysis is possible for the entire starburst regions of NGC 3256 and NGC 253. By combining published data on the integrated Lyman continuum luminosity from Rigopoulou et al. (1996) and Engelbracht et al. (1998) with the results from PDR modeling by Carral et al. (1994), we infer log U eff = −2.3 and −2.6, respectively. The value of log U eff we derive for NGC 253 is consistent with that derived by Engelbracht et al. 1998 (log U∼ −2.2 to −2.5). The determination of log U for these three galaxies supports the suggestion by Carral et al. (1994) that the ISM properties in starburst galaxies are independent of the global luminosity and the triggering mechanism, and are endemic to starburst activity itself. We therefore use the log U eff and R eff determined in M 82 as representative values for the entire sample.
In our models, n e and R = R eff are assumed to be time-independent. We first computed the neon line ratio for a very short-lived starburst, with t sc =1 Myr. In this model, U varies proportionally with Q Lyc over time, via Equation (1). To set the absolute flux scale of the cluster SEDs, we choose the condition that when each cluster reaches its maximum Q Lyc , U max = U eff ; the shape and relative flux scale as the clusters evolve are determined from the composite SEDs output by STARS. The computations are stopped when log U becomes smaller than −5.5; at this point, the SFR has dropped down by several orders of magnitude compared to the initial SFR, and the stellar population produced by the starburst is expected to have faded away into the galaxy's background population.
With the nebular parameters determined above as constraints, the SEDs produced by STARS are input into CLOUDY. Approximating a thin shell around the distributed cluster geometry, CLOUDY was run with plane-parallel geometry and an ionization parameter specified by Equation (1). Models for the burst timescales presented in this paper were computed by convolving the neon fluxes for the t sc =1 Myr burst with an exponential SFR with the appropriate value of t sc , where F λ is the flux in the line at wavelength λ. This method inherently assumes that U(t b ) is proportional to the ionizing photon rate produced in the 1 Myr-timescale burst, rather than that indicated by the integrated Q Lyc for the longer burst timescales. Thus at later times, the [Ne III]/[Ne II] ratio will stay constant, as the line emission will be produced only by the most recently formed clusters. It is these recently-formed clusters which will dominate the ionizing radiation field at all ages, and are all assumed to have the same U max .

Results for a homogeneous sample of evolving clusters
The model neon line ratios output by CLOUDY are indicated by the solid lines in the upper panel of Figure 6, as a function of burst age and upper mass cutoff for various t sc . These ratios assume no dust, n e ∼300 cm −3 , and log U max = −2.3. The model predictions for five different upper mass cutoffs (m up =25, 30, 35, 50, and 100 M ⊙ ) and three burst timescales (t sc =1, 5 and 20 Myr) are shown to represent a range of prescriptions for the formation of stars in clusters. There is a slight decline in the computed ratio for t b ∼5 Myr for the longest timescales, due to a build up of evolved clusters that soften the shape of the integrated ionizing radiation field. However, these evolved clusters must still be relatively young to affect the integrated ratio significantly. The asymptotic behaviour of the curves for timescales longer than ≈ 5 Myr can be attributed mainly to our choice of an exponentially decaying SFR, but also reflects the dominance of the youngest generations, for which log U eff is close to the maximum.
The vertical bar to the right of the plot represents the range of observed neon line ratios for all sample galaxies except the two low-metallicity blue dwarf galaxies NGC 5253 and IIZw40, which will be discussed separately in §4.3. From the comparison with models, values of m up < 25 M ⊙ can be ruled out for the entire sample. In addition, the largest neon ratios in the sample are consistent with upper mass cutoffs at least as large as 50 M ⊙ , independent of the burst age and timescale. Figure 6 shows clearly the degeneracy between age and m up : the range of observed neon ratios may either be due to true variations in the population of the most massive stars, or to short bursts of star formation where m up is always large but the SED is softened through aging. However, a similar modeling analysis of the nearby star formation sources plotted in Figure 1 ( Thornley et al. , in prep.) and extensive modeling of M82 (Schreiber 1999) support the need for stars more massive than 50M ⊙ , in agreement with stellar census measurements in a variety of nearby star forming regions.

Modeling infrared emission from a heterogeneous ensemble of evolving clusters
The models described above are appropriate for individual clusters (and surrounding nebulae) or for a homogeneous collection of clusters, enclosed within a larger region. However, the starburst galaxies we are examining may contain evolving clusters with a wide range of total masses, and low-mass clusters may not have enough material to fully sample the IMF. Here, we explore the effect of an assumed cluster mass spectrum on the integrated [Ne III]/[Ne II] ratio. On average, the ionizing spectrum is softer if small clusters contribute a significant fraction of the total ionizing luminosity.
Such an ensemble of clusters can be conveniently described by a luminosity function (LF). Various studies have explored the LF of young clusters and HII regions in a wide range of star-forming environments, including the disk of our own Galaxy (McKee & Williams 1997), nearby spiral and irregular galaxies (e.g. Kennicutt, Edgar & Hodge 1989;Elson & Fall 1985;Gonzalez-Delgado et al. 1995), as well as merger systems like NGC 4038/4039 (Whitmore & Schweizer 1995). These studies show that the number distribution of clusters in optical light and in ionizing luminosity follows a power-law LF, dN /d(log L) ∝ L −β , with remarkably similar indexes in the range β = 0.5 − 1, down to the faintest luminosities observed. The typical completeness limits correspond to an absolute V -band magnitude of ≈ −10 mag or to a Lyman continuum luminosity of ∼ 10 50 s −1 . Hence only bright clusters are directly observed. For comparison, the Orion nebula is powered by stars which emit a total of 10 48.85 Lyman continuum photons s −1 (Kennicutt 1984).
In order to model the integrated emission properties of an ensemble of clusters, we must include the contribution from all clusters in a model starburst event. While this calculation is relatively straightforward in large, luminous clusters, it is more uncertain in lower-mass clusters where statistical fluctuations at the high-mass end of the IMF become important: for instance, the formation of a single 20M ⊙ star instead of a single 15M ⊙ star could change the cluster luminosity by an order of magnitude. To determine the shape of the LF consistently for the brightest clusters as well as clusters fainter than ∼10 50 s −1 , we ran a Monte Carlo simulation to determine the stellar contents of clusters of different masses using Poisson statistics.
Our Monte Carlo simulation was run for a sample of 2x10 5 clusters, with the number of stars per cluster ranging from 100 to 10 7.5 and mass bins populated from m low =0.1M ⊙ upward in accordance with the specified Salpeter IMF. In this way, the upper mass cutoff of an individual cluster is dependent on both the total mass of the cluster and statistical variations in the population of the IMF. Using published values to represent the Q Lyc (m) relation for individual stars (Vacca, Gar-many, & Shull 1996;Panagia 1973), we then determined the average cluster mass, m cl , as a function of ionizing luminosity Q cl Lyc for the range of clusters falling in a given luminosity bin of width 0.1 dex. The resulting relation, m cl (Q cl Lyc ), can be approximated by two power laws, with m cl ∝ Q cl Lyc for clusters with log Q cl Lyc >49.5 and m cl ∝ (Q cl Lyc ) 0.19 for clusters with log Q cl Lyc <49.5. To convert from a function m cl (Q cl Lyc ) to a luminosity function, we assume a simple, power-law cluster With these conditions, we derive a LF covering Q cl Lyc from 10 44 to 10 53 photons s −1 . For the cluster mass function, we assume the simple case α=1, which produces a LF shape consistent with that determined by McKee & Williams (1997) for HII regions in the disk of our own Galaxy. The lower and upper limits for Q cl Lyc in our LF were chosen to represent the smallest associations still capable of ionizing an H II region (i.e. containing only one early-B star), and the most luminous super star clusters detected (e.g., in NGC 4038/4039). The resulting LF is again represented by a broken power law, with β=0.19 at the lower end and β=1.0, with the break occuring at Q cl Lyc ∼ 10 49.5 s −1 as for the m cl (Q cl Lyc ) relation for individual clusters. Figure 7 shows the derived LF, with the McKee & Williams LF overplotted for comparison.
We assume that the LF determined above describes the cluster ensemble at birth, and then follow the evolution of the ensemble of clusters for various star formation histories. The evolution of the LF with burst age for the ensemble of clusters is obtained by following the evolution of Q cl Lyc for the clusters in each luminosity bin. For this purpose, we computed a library of cluster models for a t sc =1 Myr burst using STARS and CLOUDY. The Monte Carlo simulations described above provide the relation between Q cl Lyc and m cl up , the mass of the most massive star in an individual cluster given a fixed IMF and an input stellar count (note that m cl up should be distinguished from an intrinsic, galaxy-wide upper mass cutoff of the IMF, which we have thus far designated as m up ).
Models were generated for t sc =1 Myr and were obtained by interpolation of the library models. At zero age, the clusters with different initial m cl up (and Q cl Lyc ) were distributed in the luminosity bins according to the LF. As the burst age increases, the evolution of each cluster in Q cl Lyc is followed (using the library models), and the new LF at each time step is determined from the distribution of the clusters in the different Q cl Lyc bins after evolution has taken place.
The slope of the high-luminosity end of the LF ( 10 50 s −1 ) changes very little during this evolution; this is mainly due to the fact that these clusters contain stars with masses >50M ⊙ . All such massive stars have similarly short main-sequence lifetimes (e.g. Schaller et al. 1992), which implies that these luminous clusters move into lower Q cl Lyc bins at similar rates. The constant slope of the high-luminosity end of the LF is thus consistent with the observed LF in various sources, which have presumably a range of ages and starburst histories.

Results for a heterogeneous ensemble of evolving clusters
The neon ratio at any given age for the t sc =1 Myr ensemble of clusters was obtained by summing over the [Ne II] and [Ne III] line fluxes of all the model clusters still contributing; the integrated properties for longer burst timescales were again obtained by convolving those for the t sc =1 Myr burst. The neon ratios predicted for an ensemble of clusters distributed according to the derived cluster LF are indicated, for the longest and shortest timescales we explored, by the dashed lines in the upper panel of Figure 6. As expected, the main effect of accounting for a cluster size distribution is that the smaller, less luminous clusters soften the integrated ultraviolet radiation field. In particular, the predicted [Ne III]/[Ne II] for an ensemble of clusters is lower than for a single cluster, assuming the same galaxy-wide m up . The effect of the LF is larger for higher m up , with a reduction in predicted ratios by a factor of ∼2 for a galaxywide m up =100 M ⊙ , and ∼1.2 for m up =25 M ⊙ . This dependence on m up reflects primarily the stellar properties themselves: the hardness of the ionizing spectrum decreases steeply as the stellar mass decreases. The magnitude of the effect of incorporating a cluster LF depends mainly on the power-law index for the LF and on the upper limit in Q cl Lyc . The steep power-law we have adopted maximizes the differences between single and ensemble cluster models, in comparison with other plausible LFs with indices above the break of β = 0.5 − 1. The choice of a high upper limit in Q cl Lyc is justified by observations in some starburst galaxies, notably NGC 4038/4039 and M 82 (Whitmore & Schweizer 1995;O'Connell et al. 1995).
Although the "down-weighting" effects of the smaller clusters, which produce softer ionizing radiation, are measurable, they do not require a significantly different conclusion than that reached by comparison of our data with models of homogeneous clusters. Accounting for a LF of the ionizing clusters which excite the observed nebular emission lines, we confirm that very highmass stars can form in starburst galaxies, allowing m up = 50−100 M ⊙ even more easily in the sources with the highest measured ratios. If we posit that all starbursts form stars in the manner that nearby, massive star forming regions do, then the most plausible explanation for low observed line ratios is that starbursts are events of short duration ( Figure 6 suggests t sc ∼1 Myr) which produce very massive stars, but whose aging rapidly softens diagnostic ratios such as the ones we use here. Our models show that as long as very massive stars are formed, even in small numbers, they strongly dominate the ionizing radiation field and thus maintain high neon line ratios. The ratios can only decrease to the observed range long enough after the exhaustion of starburst activity for the most massive stars to have evolved off the main sequence. For ratios near 0.1, corresponding to stars with initial masses near 30 M ⊙ , this will take about 5 Myr. In principle, the neon ratios will start declining very rapidly after the last massive stars have formed in the burst, so that the neon ratio alone does not fully constrain the timescale. However, by combining with other measurements which are indicative of starburst properties, we can explore the robustness of our conclusion that starburst timescales are generally quite short.
The ratio of the bolometric to Lyman continuum luminosities (L bol /L Lyc ) is a useful contrasting probe of the properties of massive stars: since this ratio is sensitive to a somewhat lower mass range than that to which [Ne III] 15.6 µm/[Ne II] 12.8 µm is sensitive, it varies significantly with time even for longer burst timescales. This reflects the buildup of a population of stars in a lower mass range which contribute more to L bol , and which have longer main-sequence lifetimes. Therefore, the combination of both diagnostics can help address the degeneracy between aging effects and variations in the upper mass cutoff.
As we are characterizing the properties of massive star forming regions, we assume L bol = L IR (as defined in §2.2). This is usually a good approximation since a large fraction of the energy output of OB stars is absorbed by the surrounding interstellar dust, present in large amounts in starbursts, and re-emitted in the thermal infrared. We derived L Lyc using two different methods. First, we use measurements of hydrogen recombination line fluxes and thermal radio continuum emission, when such data were available from observations with the SWS and in the literature. This allowed us to constrain simultaneously the extinction toward the sources from the comparison of observed and theoretical relative line fluxes. Second, we used our own neon line fluxes corrected for extinction and assuming all neon atoms are either singlyor doubly-ionized.
We assumed case B recombination coefficients and line emissivities (Hummer & Storey 1987), with an electron density of n e = 100 cm −3 and temperature of T e = 5000 K, except when individual determinations were available in the literature. For the extinction correction, we adopted a composite extinction law made up of the Rieke & Lebofsky (1985) law for λ ≤ 0.9µm and the Draine (1989) law for 0.9µm< λ <40µm. The effects of obscuration were neglected for λ > 40µm. If sufficient data were available, we constrained the geometry of the emission sources and obscuring dust as well. Two models were considered: a uniform foreground screen (UFS) of dust and a homogeneous mixture of dust and sources (MIX). Otherwise, we considered only the UFS model.
For the sources for which the determination from recombination lines was possible, the estimates from the hydrogen and neon lines agree to within a factor of three or better, and we adopted the average as the final L Lyc . The cases where the extinction could not be reliably constrained yield lower limits on L Lyc , and thus upper limits on L IR /L Lyc . The lower panel of Figure 6 shows the models for L IR /L Lyc obtained with STARS. As in the upper panel, the ratios expected for a homogeneous cluster population are shown as solid lines and those expected for a cluster ensemble defined by our LF are shown as dashed lines. The effect of including a cluster LF is of similar magnitude for L IR /L Lyc as for the neon ratio. It is clear that the average neon and L IR /L Lyc ratios are consistent with conditions where the clusters have high upper mass cutoffs, m up 50−100 M ⊙ , and short burst timescales of a few million to ∼10 7 years. In fact, by plotting the neon ratios against the L IR /L Lyc ratios in Figure 8, we see that the models suggest very short timescales, so short as to be difficult to produce even with the t sc =1 Myr burst, the burst with the shortest timescale considered in this study. Though this discrepancy is model-dependent, it is generally clear that short timescales are needed to reproduce the range of ratios observed in our starburst sample.

The Galactic Center as a template for a short starburst
The short timescales and ages inferred for the star forming activity in starburst galaxies are reminiscent of the ones in the young stellar clusters in the central parsec of our Galaxy, and can be checked there on the basis of the existing stellar census. The observed [Ne III]/[Ne II] ratio is even lower than in most starbursts, despite the high ionization parameter log U = −1 derived for the spatially resolved ionized region ( 1 pc, Lutz et al. 1996). The picture of a short but aged burst is supported by the direct stellar census (Genzel et al. 1994;Krabbe et al. 1995;Najarro et al. 1997) which suggests a star formation event of age ∼7 Myr and duration approximately 4 Myr. This event is most directly indicated by the presence of both cool red supergiants like IRS7 and massive, moderately hot (20,000 to 30,000 K) blue helium-rich supergiants or Wolf-Rayet (WR) stars.
We have computed models optimized for this Galactic Center region, assuming a single cluster with a Salpeter IMF, t sc =1 Myr, m up =100, solar metallicity, and fixed log U=−1, n=3000 cm −3 as derived by Lutz et al. (1996). The burst timescale is somewhat shorter than that suggested by the stellar census and was chosen as a conservative assumption, since a shorter burst will produce softer radiation fields at late ages and thus minimize any need to invoke other effects for explaining the observed soft radiation. The decay of line ratios with time is similar to Figure 6, but the ratios at any given time are higher due to the higher log U in this smaller region. At the age of 7 Myr preferred by the stellar census, the [Ne III]/[Ne II] ratio is still about 1 to 2, well above the value of 0.05 observed by Lutz et al. (1996). The model reaches the observed value only after more than 13 Myr, which is difficult to reconcile with the presence of evolved massive stars approaching 100 M ⊙ (Najarro et al. 1997;Ott, Eckart, & Genzel 1999). While aging effects push the Galactic Center neon ratio more closely into agreement with models, they are not enough to fully account for the low observed neon ratio.
The stellar census for the Galactic Center gives an independent view through direct analysis of the contributions of different stellar types to the ionizing continuum. Figure 9 presents an Hertzsprung-Russell diagram weighted by Lyman continuum luminosity, derived from the same starburst model for an age of 7 Myr. It is evident that most of the ionizing luminosity still originates in stars close to the main sequence at log T= 4.5 to 4.6, with smaller contributions by stars somewhat evolved towards lower temperatures and a population of hot WR stars to the left of the main sequence. This is in stark contrast to the finding by Najarro et al. (1997) that 7 of their stars, found at log T < 4.5 and log L > 5.75, contribute half of the ionizing luminosity of the Galactic Center. The same region holds less than 1% of the ionizing luminosity for the model in Figure 9. At all other model ages, this fraction does not exceed about 1%.
The fact that the discrepancy between 'hard' models and 'soft' observations is seen both in the photoionization modelling and the stellar census leads us to the conclusion that this discrepancy is not primarily due to uncertainties in the hot star SEDs used for the photoionization models. Such uncertainties affect mainly the shape of the ionizing flux and thus the fine structure line ratios, but much less the total ionizing flux of the star and the census. Another possibility is that the adopted Geneva stellar evolutionary tracks (Schaller et al. 1992) do not predict the large number of 20,000 K to 30,000K supergiants seen in the Galactic Center, or that there exists a mismatch between the stellar effective temperatures from evolutionary tracks and those from stellar atmosphere models. The corresponding region is populated by these tracks but only for short intervals, possibly because massive, post-mainsequence stars rapidly move off to hotter parts of the tracks. The disagreement between tracks and observations may be related to the difficulty of defining mass loss, atmospheres, and effective temperatures for late stages of massive star evolution governed by strong winds. Our conclusions are unchanged when using evolutionary tracks with twice solar metallicity (Schaller et al. 1992) and higher mass loss (Meynet et al. 1994), which might be adequate for the somewhat higher metallicity in the Galactic Center. The latter tracks in fact increase the disagreement by adding more very hot WR stars. In fact, WR stars may contribute to the uncertainties in the evolutionary tracks, as our knowledge of the ionizing spectra of WR stars is uncertain. Indeed, there is both observational and theoretical evidence (Crowther et al. 1999;Hillier & Miller 1998) which suggests that the ionizing spectra of WR stars may be much softer than commonly assumed.
The case of the Galactic Center has obvious implications for aging starbursts with similar populations. If the current Geneva tracks, interpreted in the fashion described above, indeed mispredict the post-main-sequence evolution of massive stars by postulating very hot stars where the Galactic Center shows us a cooler luminous population, then starburst models will tend to predict too high values for ratios like [Ne III]/[Ne II]. The resultant "need" to invoke low upper mass cutoffs (or aging effects) to explain soft radiation fields may be due to this effect. Corrections for this effect would likely lower the neon ratios, and thus relax the stringent conditions of short burst time scales, leading to more plausible timescales of 10 7 years. While a careful study is needed to test the reliability of current evolutionary tracks, this analysis of the Galactic Center supports our general conclusion that starbursts are "normal" star-forming environments, in the sense that they produce very massive stars just as local star forming regions do.

The effects of parameter variations
Starburst modeling is, by the very nature of the star formation process, a multi-parameter problem. To assess the robustness of the results we have presented, we now discuss the effects of varying the possible input parameters of our model, considering variations in starburst SEDs and nebular parameters. The most significant changes in a starburst SED will arise because of our choice of metallicity, model SED libraries for individual stars, or the shape of the IMF. Changes in nebular parameters have effects of similar magnitude as those of changes in SEDs. The most significant effects will be due to variations in the ionization parameter, though variations in the dust population as well as the density of the gas in these regions also have minor effects. For simplicity, these variations are illustrated for single-cluster models with m up =100M ⊙ and t sc =1 Myr in Figure 10.

Metallicity
Of the above effects, those of metallicity are the most significant. Metallicity effects are twofold: at sub-solar metallicities, the stellar SEDs are harder due to reduced line blanketing and blocking (e.g., Pauldrach et al. 1998), and the evolutionary track each star follows through the H-R diagram changes, producing a "hotter" main sequence. Figure 10a compares the neon ratios for two models, the hybrid SED grid at solar metallicity, and the corresponding hybrid grid with stars and gas at 0.2 Z ⊙ , which is more appropriate for the two low-metallicity dwarf galaxies in our sample. Lowering the metallicity results in considerably increased neon ratios: for Z=0.2Z ⊙ , the predicted neon ratios increase by factors of ∼4-10. The increase in neon ratios at low metallicity is dominated by changes in the evolutionary tracks and SEDs of model stars, with the corresponding changes to gas-phase abundances playing a relatively minor role. NGC 5253 and IIZw40 have neon ratios 3 and 13 times higher, respectively, than in any of the other starbursts in our sample; the observed ratio values are also plotted as horizontal grey lines in Figure 10a. Therefore, even accounting for low metallicity, the neon ratios measured for these systems are consistent with a stellar population with m up >50-100M ⊙ . In contrast with the example of Figure 10a, models with Z> Z ⊙ will have correspondingly lower predicted neon ratios; this factor is worth considering if more accurate determinations of abundances become available.

SED libraries
If we were to choose Kurucz instead of Pauldrach SEDs to represent the most massive stars, the contribution of the softer Kurucz spectra for high-mass stars would produce lower neon ratios (cf. Figure 4). The significance of this effect can be seen in Figure 10b, which shows the neon ratio for two input stellar grids: the interpolated, hybrid grid we are using, and a standard Kurucz model grid. At a burst age of ∼3 Myr, the predicted neon ratio using Kurucz SEDs is a factor of ∼2 lower than that predicted by our hybrid grid. Comparison with our ISO-SWS data would suggest m up >50M ⊙ for more than half the sample, even at zero-age.

IMF slope
We have taken a Salpeter IMF to be the most representative form of the initial mass function; however, the universality of the IMF is still debated (see, e.g., Massey 1998 andScalo 1998 for contrasting views). Some authors have suggested that the IMF in active star-forming environments is well-represented by a Salpeter form (Hunter 1995;Massey & Hunter 1998). While the alternative model favored by Scalo (1998) exhibits a steeper function at intermediate masses, it is similar to Salpeter at the high-mass end (M≥10M ⊙ ), where the neon ratio is most sensitive. To illustrate the effects of changing the IMF, we show the resulting neon ratios for a Miller & Scalo (1979) IMF in Figure 10b, where this IMF is represented as a power law with index −1.4 for M=1-10M ⊙ and −2.5 for M=10-100M ⊙ . The Miller-Scalo IMF is generally steeper than Salpeter, and thus there are fewer massive stars formed relative to a given number of low-mass stars. The net result is a softer composite SED and a prediction of lower neon line ratios. The decrease in the predicted neon ratio caused by changing the IMF in this way is smaller than that caused by changing SED libraries. Figure 10c shows the predicted neon ratios for two alternate values of the ionization parameter: log U max = −1.5, the upper limit for U assuming that all clusters lie at the center of the starburst nebular emission region of M82 (Schreiber 1999), and log U max = −3.5, a value more similar to that derived by Wang, Heckman, & Lehnert (1997) for the diffuse ionized medium (DIM). An ionization parameter of log U max = −1.5 would result in predicted neon line ratios 3 times higher out to ∼7 Myr for t sc =1 Myr. Such an increase would imply upper mass cutoffs between 25 and 50 M ⊙ for any bursts with t sc 5 Myr, without accounting for any aging effects. This value of U max is a factor of ∼3 greater than the highest value which is consistent with the M82 analysis. The effect of reducing the ionization parameter to log U max = −3.5 is larger, causing a decrease in the predicted neon ratios by an order of magnitude. The highest observed ratios in our starburst sample would not be reproduceable, even for the most massive stars for which we have models, if log U max = −3.5. Previous starburst modeling studies along these same lines (see, e.g., Kunze et al. 1996, Rigopoulou et al. 1996 have used log U=−2.5, similar to the M82-based value used in this study.

Dust population and gas density
The presence of dust within HII regions can affect the efficiency of nebular photoionization, but introducing a dust component has a relatively small effect on our models. Figure 10d shows the effect of adding dust grains similar to those in Orion (we use the Orion dust population that is incorporated in CLOUDY, from Baldwin et al. 1991). Adding such a dust component causes a ∼15% increase in the predicted line ratios. This variation is smaller than the uncertainties in the measured neon line ratios, and negligible compared to the other parameter variations we have explored. We conclude therefore that we have introduced no significant uncertainties by neglecting dust in earlier sections. Note, however, that the uncertainties are larger for determinations of L bol /L Lyc , where dust could have a much more significant effect.
Due to the large uncertainties in inferring the gas density from the SIII (18.7/33.5µm) ratio, we also show in Figure 10d the effects of increasing the gas density to 10 3 cm −3 . This is the highest density consistent with the range of measured sulfur ratios and their uncertainties. For the model shown, the neon line ratio increases by less than five percent over the entire age range, indicating that variations in gas densities have a negligible effect on the output neon ratios. However, we note that the ionization parameter changes inversely with gas density, such that this increase in gas density would also imply a drop in the ionization parameter to approximately log Usim − 2.8, thus solidifying the case for very massive stars being present in the starbursts in our sample.

Extra-starburst contributions?
The models we have presented thus far assume that the only contribution to the MIR neon line fluxes comes from direct photoionization by stars in the starburst region itself. In this section we examine constraints, from ISO spectroscopy, on the possible contributions of two other processes: excitation by active galactic nuclei (AGNs), and contributions from a diffuse ionized medium (DIM) in the surrounding galaxy (e.g., Lehnert & Heckman 1994;Wang et al. 1997).

AGN contributions
The MIR spectra of AGNs are distinctive in their display of strong emission lines from highly ionized species, such as [Ne V] and [O IV], which require higher excitation that can be produced even by the hottest stars. The absence of these lines, or their weakness relative to lower excitation lines, has been used to demonstrate the dominant contribution of star formation to the power produced in ultraluminous IR galaxies (Lutz et al. 1996;Genzel et al. 1998). Though some starbursts show very weak [O IV] emission, the most plausible explanation is ionizing shocks from supernovae or superwinds (Lutz et al. 1998 emission. we therefore consider the measured neon ratios for NGC 6240 and NGC 7469 to be upper limits to the emission arising from the starburst region. With these exceptions in mind, we consider the contribution from an AGN to be unlikely across the sample, and any contribution must have a neglible effect on the results presented here.

DIM contributions
Studies of the Milky Way and other nearby galaxies suggest the presence of a "diffuse ionized medium" (DIM), a gas component with a relatively low ionization state and large scale height which permeates the galaxy disk (e.g., Reynolds 1990;Dettmar 1992). It may be difficult to exclude a DIM component as a contributor to the neon line fluxes measured for this sample. Several studies have shown the existence of diffuse, ionized emission in the disks of galaxies, which is generally not associated with individual star forming regions. For the more distant objects in our sample, the aperture covers a large physical area (the long axis of the SWS aperture corresponds to linear diameters of 0.3-14 kpc for the galaxies observed), which may encompass non-starburst emitting regions elsewhere in the galaxy. Thus, we must consider the possibility that some form of DIM emission may influence the integrated line fluxes. Lehnert & Heckman (1994) and Ferguson et al. (1996) showed evidence that the DIM may be produced by ionizing starlight escaping from HII regions in the disks of galaxies; Wang et al. (1997) measured its effects on large scale measurements of optical excitation ratios such as [NII/Hα] and [SII/Hα]. These studies suggest that the DIM may contribute as much as 50% of the global, integrated Hα flux in spiral galaxies. If we assume that the DIM consists of ionizing radiation escaping from young star clusters, with the same average parameters as derived by , n e ∼1 cm −3 ), then the model neon ratios could drop by a factor of two, making the effect of a DIM contribution similar in magnitude to that of variations in the ionization parameter that were discussed in §4.3.
We cannot exclude the contribution from a lowdensity component, as the [S III] ratio generally provides only an upper limit to the density. How-ever, we have a qualitative constraint supplied by the range of neon line ratios observed in our sample: the observed neon ratios show no correlation with distance, as seen in Figure 11. If the DIM were a significant contributor, we might expect the neon ratio to decrease with increasing source distance, as the aperture encloses an ever larger physical area. In the nearest galaxy in our sample, IC 342, the aperture encloses an area ∼200 pc in diameter, and for galaxies at a distance of 30 Mpc the aperture still covers a region less than 5 kpc in diameter. Thus, unless we are observing objects in which the size of the starburst area grows in proportion to its distance, we conclude that any DIM-like component in the sample galaxies does not make a significant contribution. Any DIM contributions will result in lower predicted neon ratios; thus, our conclusion that the most massive stars are generally formed in all starburst environments is not affected.

Discussion and Summary
The cluster models that we have presented in this paper support the formation of very massive stars (50-100M ⊙ ) in starburst galaxies. While the quantitative estimate of m up for each galaxy is model-dependent, it is clear that the formation of very massive stars is necessary to explain the ionized line diagnostics observed in this starburst sample. This result suggests that while starbursts produce prodigious amounts of energy and stars, the high-mass stellar populations in starburst galaxies are not radically different than those in high-mass star-forming regions observed locally.

Short timescales for starburst activity
As illustrated in Figures 6 and 8, the combination of the neon line ratio with the L IR /L Lyc ratio strongly favors the scenario for starburst activity where very massive stars form, as in local smaller-scale starburst templates, and where the burst last typically a few million to ∼10 7 years. Such timescales are shorter than previously thought (10 7 − 10 8 yr; e.g., Thronson & Telesco 1986;Heckman 1998). It is clear that detailed modeling is required to secure this result, including additional constraints (e.g. K-band luminosity, the rate of supernova explosions, the depth of the near-infrared CO bandheads) and spatially resolved information. It is nonetheless in agreement with other recent detailed studies of a few starburst galaxies, some of which are also included in our sample (e.g. M 82, Schreiber 1999;NGC 253, Engelbracht et al. 1998). As a result of instrumental progress, there is now growing evidence that starburst activity occurs in individual burst sites on physical scales of a few tens of parsecs or less. Short timescales are therefore naturally understandable locally. Our data, in conjunction with the other studies cited above, provide evidence for short timescales on much larger scales, suggesting that starburst activity also occurs globally on short timescales, presumably as a result of one brief gas compression event, or of successive episodes of such events separated by more than one typical timescale. Short burst durations thus imply strong negative feedback effects of starburst activity, globally as well as locally.
A simple argument can be invoked to explain the physical arguments behind this result. We can compare the cumulative kinetic energy injected in the ISM by the supernova explosions over time (E kin ) with the binding energy of the gas (E grav ), and assume the starburst activity will stop when E kin just balances E grav . This is a simplistic way of expressing the conditions for a starburst wind to break out of the galaxy (e.g. Heckman, Armus & Miley 1990), but it is sufficient for orderof-magnitude estimates.
In order to relate E kin to observed quantities, we have considered the relationship between the rate of supernova explosions ν SN and the bolometric luminosity L bol . Model predictions obtained with STARS for a variety of star formation histories and upper mass cutoffs of the IMF indicate that as soon as the massive stars start to explode as supernovae, and as long as substantial star formation takes place. It thus seems reasonable to assume that Equation (4) holds for the bulk of the sample, likely having a range in age and timescale but all exhibiting signs of significant, recent starburst activity. For simplicity, we here assume a spherical geometry for the starbursts, with uniform mass distribution. In addition, we assume that each su-pernova explosion liberates E SN mech = 10 51 erg of mechanical energy, transferred as kinetic energy to the ISM with an efficiency η. The timescale τ for our condition above satisfies : where dE SN mech /dt is the rate of mechanical energy injection from the supernovae, G is the gravitational constant, M dyn is the dynamical mass of the system, and R is the radius of the starburst region. Equation (5) can be re-written as where we have substituted L IR for L bol , appropriate for our sample galaxies. Application of Equation (6) to M82 (M dyn = 8 × 10 8 M ⊙ , R = 0.25 kpc, L IR = 4 × 10 10 L ⊙ ), yields τ ≃ 5 η −1 Myr, so for efficiencies 10%, the estimated timescales are ∼ 10 6 − 10 7 yr.
Our argument above is based on "gas-disruption timescale" estimates. This differs from the conventional "gas-consumption" arguments, which compare the star formation rates with the mass of the gas reservoir. In such estimates, the star formation rates are based on comparison of absolute fluxes (e.g. Hα fluxes, L IR ) with predictions from starburst models. The estimates are thus very sensitive to the assumed age and history of the starburst. Our estimates of the gasdisruption timescales are also model-dependent, but have the advantage of being based on a quantity (10 12 ν SN /L bol ) which varies by smaller factors. Neither point of view accounts for further fueling processes, or dynamical evolution of the systems (e.g. starbursts in barred galaxies, interacting/merging systems, etc.). However, the discussion presented here gives an alternative perspective to the issue of global burst timescales, and provides a plausible explanation for our results.

Summary
Starburst models predicting the [Ne III]/[Ne II] ratio from ISO-SWS spectra of 27 starburst galaxies show that the observed data are consistent with the formation of very massive stars in starbursts, thus precluding the need for the restrictive upper mass cutoffs suggested by some earlier studies (m up ∼25-30M ⊙ ). Combining the neon line ratios with starburst modeling and the consideration of the stellar content measured in local star forming regions, we find that starburst events may be generally described as short bursts of star formation which produce very massive stars, and which exhibit relatively soft integrated line ratios as a result of aging the stellar population.
In particular, our modeling of neon and L IR /L Lyc ratios, together with results on local high-mass star-forming regions, suggest: • very massive stars (m up 50M ⊙ ) form in typical starbursts.
These results suggest strong negative feedback from starburst activity; the galactic superwinds frequently observed in starburst galaxies are particularly striking examples of the consequences of such feedback.
In our analysis, we have examined the degeneracy between aging effects and model parameter variations in the assessment of upper mass cutoffs to the IMF. There is still room for significant improvements in modeling the properties of starbursts: determination of metallicities and the radiation environment (e.g., for measurements of the ionization parameter U) compete with the characterization of stellar properties (SEDs, evolutionary tracks) as the largest contributors to uncertainty in the modeling of star formation properties such as the upper mass cutoff. Other datasets, such as additional MIR line ratios (e.g., Kunze et al. 1996;Rigopoulou et al. 1996;Engelbracht et al. 1998) or K-band luminosities and near-infrared spectroscopy (e.g., Forbes et al. 1993;van der Werf et al. 1993;Genzel et al. 1995;Tacconi-Garman et al. 1996;Böker, Förster-Schreiber, & Genzel 1997;Engelbracht et al. 1998;Schreiber 1999) would be very useful in further constraining the properties of starbursts. However, it will be important to compile such additional data for a large sample in order to proscribe further, general constraints on the way in which starbursts form stars. Observations with higher spatial resolution would better isolate regions of active star formation, making it possible to confirm whether highand low-excitation lines arise from the same region; the spectroscopic capabilities of SIRTF will be well-suited to addressing this issue. By accounting for a reasonable range of uncertainties which constrain the present observations, we find that the observed MIR neon ratios are generally consistent with the formation of very massive stars in starburst events; we offer this hypothesis up to future datasets for increasingly rigorous testing.
We would like to thank A. Pauldrach and R.-P. Kudritzki for providing model atmospheres, Tal Alexander for assistance in introducing lowmetallicity SEDs into STARS, and Eckhard Sturm for providing the ISO-SWS spectra of Arp 220. We would also like to thank Jack Gallimore and Dan Tran for interesting discussions. MDT would like to thank the Alexander von Humboldt-Stiftung and the NRAO 5 for support during the production of this work. This research also received support from the German-Israeli Foundation under grant I-0551-186.07/97. SWS and the ISO Spectrometer Data Center at MPE are supported by DARA under grants 50 QI 8610 8 and 50 QI 9402 3.   Squares indicate Kurucz models, triangles indicate Pauldrach models, and stars indicate "hybrid" models interpolated between Kurucz and Pauldrach models (see text). Original ("cornerstone") models are shown as filled symbols, and interpolated and extrapolated models are shown as open symbols. Models for log g>4 are identical to log g=4 models, and models with T ef f >60,000 K are extrapolated from the highest temperature Pauldrach model.     Fig. 9.-Contributions of different stellar types to the total Lyman continuum luminosity of a short starburst (t sc =1Myr) at a burst age of 7 Myr. Contours are at 10 −6 , 10 −5 , ..., 10 −2 contribution of a ∆ log L = 0.25, ∆ log T=0.0125 bin to the total luminosity. The greyscale coded regions are above the 10 −3 level and contribute more than 97% of the total flux for this model. The box enclosed by a dotted line indicates the region which hosts stars contributing about 50% to the ionizing luminosity of our Galactic Center. The measured neon ratios for IIZw40 and NGC 5253 are indicated by the upper and lower grey lines, respectively. (b) Neon ratios for a burst composed of SEDs from our hybrid grid, and from a grid composed solely of Kurucz models. In addition, the neon ratios for the hybrid grid using a Miller-Scalo IMF are also shown. (c) Neon ratios for log U=-1.5,-2.3, and -2.5. (d) Neon ratios determined in environments with gas densities of 300 and 1000 cm −3 ; for the case of 300 cm −3 , we also show the effect of adding an Orion-like dust component to the CLOUDY modeling.  a Infrared (λ = 8 -1000 µm) luminosities computed from the point source flux densities listed in the IRAS Faint Source Catalog, version 2.0 (1990), using the prescription from Sanders & Mirabel (1996). The infrared luminosities were further scaled to match the SWS aperture (see §2.2). In some cases, no information was available to provide a reliable scaling between the emission seen by IRAS and that expected in the SWS aperture; for these cases, the total L IR computed from the IRAS Faint Source Catalog flux densities is indicated as an upper limit.
b From open-time sample.
j L Lyc used in this ratio is the average of the values determined from HI recombination lines and from MIR neon lines. See §4.1.