A New Method for Obtaining the Star Formation Law in Galaxies
Abstract
We present a new observational method to evaluate the exponent of the star formation law as initially formulated by Schmidt, i.e. the powerlaw expression assumed to relate the rate of star formation in a volume of space to the local total gas volume density present there.
Total volume densities in the gas clouds surrounding an OB association are determined with a simple model which considers the atomic hydrogen as a photodissociation product on the cloud surfaces. The photodissociating photon flux incident on the cloud is computed from the farUV luminosity of the OB association and the geometry. As an example, we have applied this “PDR Method” to a sample of starforming regions in M33 using VLA 21cm data for the Hi and GALEX imagery in the farUV. With these two observables, our approach provides an estimate of the total volume density of hydrogen (atomic + molecular) in the gas clouds surrounding the young star cluster.
A graph in logarithmic coordinates of the cluster UV luminosity versus the total density in the surrounding gas provides a direct measure of the exponent of the star formation law. However, we show that this plot is severely affected by observational selection, which renders large areas of the diagram inaccessible to the data. An ordinary leastsquares regression fit to a straight line therefore gives a strongly biased result. In the present case, the slope of such a fit primarily reflects the boundary defined when the 21cm line becomes optically thick and is no longer a reliable measure of the Hi column density. We use a maximumlikelihood statistical approach which can deal with truncated and skewed data, and also takes account of the large uncertainties in the total gas densities which we derive. The exponent we obtain for the Schmidt law in M33 is 1.4 0.2.
Subject headings:
galaxies: individual (M33) — ultraviolet: galaxies — galaxies: ISM — ISM: clouds — ISM: molecules — ISM: atoms1. Introduction
Some 50 years ago, Maarten Schmidt published a pair of seminal papers titled “The Rate of Star Formation” that discussed the relationship between the interstellar gas and the present and past rates of star formation in the Galaxy (Schmidt, 1959, 1963). This was a topic of considerable interest at the time, since radio telescopes had begun to map the distribution of atomic hydrogen over the Galaxy in some detail. Schmidt postulated a simple powerlaw relationship between the volume density of the interstellar gas in a region of the Galaxy and the number of stars per unit volume and time formed there. He wrote this postulate as:
(1) 
where SFR is the local star formation rate e.g. in solar masses per cubic parsec per year, and is the total volume density of the interstellar gas. For the former, Schmidt used the Zdistribution and counts of young stars. For the latter, he used estimates of the density and Zdistribution of atomic hydrogen, while at the same time acknowledging that these values might be wrong owing to the unknown amount and Zdistribution of molecular hydrogen which may also be present. With these simplifying assumptions, Schmidt found that the exponent in equation 1 was approximately 2 in the local Galactic neighborhood.
The difficulty of obtaining estimates of total gas volume density in the interstellar medium (ISM) has prevented further attempts to obtain estimates of the star formation law as initially formulated by Schmidt. Instead, the suggestion was made to relate a star formation rate per surface area to the surface density of atomic hydrogen as determined from the surface brightness in the 21cm Hi line (e.g. Buat et al., 1989). While this approach necessarily ignores differences in the lineofsight distributions of the gas and young stars, it is straighforward to implement. For the last years this approach has been the mainstay of observational studies of the star formation law in disk galaxies at low inclinations, and the formulation is now called the “KennicuttSchmidt Law” in recognition of the extensive work in this area by R. Kennicutt. These studies are summarized in Kennicutt (1998), where it is concluded that the exponent equals 1.4 in an expression for surface densities of the form . It is then generally assumed that the determined from surface densities is the same as the for volume densities in our equation 1. In more recent studies the 21cm Hi data has been augmented with estimates of the molecular content using the surface brightness of the CO(10) line as a tracer, but the general approach has otherwise remained the same. Dobbs & Pringle (2009) have recently summarized the subject and shown that a simple physical model of star formation in the ISM is capable of reproducing a powerlaw form for the KennicuttSchmidt law.
The existence of a powerlaw form for the star formation has been confirmed in many galaxies, although a wide range of values has been reported (Kennicutt, 1997). The recent addition of CO as a tracer for the molecular gas has not reduced the scatter in the results. For instance in M33, Heyer et al. (2004) find a slope of 1.38 for the molecular gas (based on CO measurements) but 3.3 for the total gas relation (with a monotonicallyincreasing atomic gas fraction at increasing galactocentric distances) and Verley et al. (2010) report values varying from 1.1 to 2.9. A number of concerns can be expressed about the meaning of these results given the input data. For instance, what are the consequences of averaging along the line of sight when we know from Galactic studies that the spatial distributions of the tracers are not the same on scales below a few hundred parsec? And how reliable are the tracers we are using for the gas (the beamsmoothed 21cm and CO(10) line surface brightnesses) in providing column densities independent of local physical conditions and radiative transfer effects?
In view of such concerns and of the wide range of values reported for the exponent in the KennicuttSchmidt Law, it seems useful to develop other methods for examining the quantitative relationship between gas and young stars in galaxies. In this paper, we present a new method based on treating the Hi found near OB associations as arising in photodissociation regions (PDRs) which develop on the surfaces of the parent giant molecular clouds (GMCs) under the action of the farUV radiation produced by nearby young stars. This approach, which we call the “PDR method”, provides us with a way of estimating the total gas volume density (atomic + molecular) in the parent GMCs. The farUV luminosity of the nearby OB association is a measure of the star formation rate at that location in the galaxy, so that a diagram of these luminosities vs. the estimates of volume densities in the surrounding clouds (of which there may be several) is amenable to interpretation as a star formation law in the sense originally described by Schmidt.
2. The PDR method
Molecular hydrogen is difficult to observe directly in the ISM, owing primarily to its lack of a dipole moment. The most common indirect means of inferring its presence in the ISM is to observe the lowlevel rotational lines of the carbon monoxide (CO) molecule. The alternative method we present here is based on the physics of photodissociation of molecular hydrogen and the discovery by Allen et al. (1986) that this process was responsible for the production of atomic hydrogen in the prominent spiral arms of M83. The method was first featured in Allen et al. (1997), then in Smith et al. (2000), Heiner et al. (2008a), and Heiner et al. (2008b). The method provides estimates of the total hydrogen volume densities (atomic + molecular) in gas clouds located in close proximity to regions of recent star formation. These clouds may be considered representative of the parent GMCs out of which the young stellar association formed. In Heiner et al. (2009) we presented the first results of applying the PDR method to a study of GMCs in M33 at a linear resolution of 20 pc, the highest available in our observations. Here we use the data from a larger survey of clouds (Heiner et al., in preparation) at a linear resolution of pc to specifically address the star formation law in this galaxy. We begin by describing the geometry and reviewing the relevant physics of our model.
2.1. A simple PDR model
Clusters of young OB stars are copious producers of farUV (FUV) photons in galaxies, and the FUV luminosity of the stellar association is a measure of the local star formation rate (see e.g. Kennicutt, 1998). The FUV photons also radiate into the surrounding interstellar medium and dissociate some of the molecular gas found in the remaining GMCs. Owing to the high porosity of the ISM, the dissociating photons penetrate hundreds of parsecs into the interstellar medium (as confirmed by e.g. Heiner et al., 2009), creating “skins” of atomic hydrogen on the surfaces of the GMCs which they encounter; see Figure 1. As we shall show, the Hi column densities so produced on the surfaces of the GMCs can be directly related to the total hydrogen volume density in the cloud; that gas will be mostly atomic on the cloud surfaces, but predominantly molecular deep inside the GMCs. The resulting total gas volume densities obtained for samples of GMCs surrounding a selection of OB associations is then combined together with the local FUVdetermined star formation rate for each association to obtain an estimate for the local star formation law in the galaxy.
A key step in the modeling is to use the physics of H dissociation by FUV photons in order to relate the observed Hi column densities on the cloud surfaces to the total gas densities in the cloud. We use a simple “slab” model described initially by Sternberg (1988) and in more detail by Allen (2004), with improved coefficients provided by Allen et al. (2004).
The Hi column density in a PDR is calculated with the same physics used to determine the excitation of the H nearinfrared fluorescence lines (Sternberg, 1988). We use the formulation in Appendix A of Allen et al. (2004). Briefly, the model is a simple semiinfinite slab geometry in statistical equilibrium with FUV radiation incident on one side, and a Hi H dissociationreformation equilibrium in the slab on the right side. The solution appropriate for our present purposes gives the steady state Hi column density along a line of sight perpendicular to the face of the slab as a function of , the incident UV intensity scaling factor (see Appendix B of Allen et al., 2004, for a definition of ), and the total volume density of H nuclei in the slab. The result is:
(2) 
where is the unattenuated H photodissociation rate in the average interstellar radiation field, is the H formation rate coefficient on grain surfaces, is the effective grain absorption cross section per H nucleus in the FUV continuum, is the incident UV intensity scaling factor, is the Hi column density, and is the volume density of H nuclei in the slab. This equation has been developed using a simplified threelevel model (Sternberg, 1988) for the excitation of the H molecule. It is applicable for lowdensity ( cm), cold (T K), isothermal, and static conditions, and neglects contributions to from ion chemistry and direct dissociation by cosmic rays. The quantity is a dimensionless function of the effective grain absorption cross section , the absorption self–shielding function , and the column density of molecular hydrogen :
The function becomes constant for large values of due to self–shielding (Sternberg, 1988). Using the parameter values in this equation adopted by Allen et al. (2004), and writing explicitly the dependence on the dusttogas ratio (Allen, 2004), we have:
(3) 
where is , the (background subtracted) atomic hydrogen column density (in ), is the dusttogas ratio in the ISM expressed in terms of the value in the neighborhood of the sun, is the (background subtracted) incident FUV flux calculated at the affected Hi patch, and is the total baryonic gas volume density (atomic + molecular). This gas will be mostly atomic on the surface of the GMC and mostly molecular deep inside the cloud. If we can measure , , and , equation 3 can be inverted to solve for ; this is the essence of our approach to obtaining volume densities in interstellar clouds. Although the exponential form of the result provides only a noisy estimate of the volume density, many measurements of can be made on GMCs in the immediate neighborhood of an OB association, thereby improving the overall precision.
Note that the method we have described does not make use of the observed Hi column density to calculate the gas volume density from some estimated (but not measured) thickness. Such an estimate of Hi volume density may, however, provide a useful “lower limit” check on the value of calculated by the PDR method, since the transition region from Hi to H is known to be very abrupt in most PDRs.
To summarize, the steps involved in applying the PDR method are:

Identify the sample of OB star clusters in the galaxy around which candidate PDRs can be sought. The clusters are knots of bright UV emission, and the candidate PDRs are patches of Hi surrounding them.

Determine the FUV flux of the OB clusters, e.g. on GALEX FUV images. The effective wavelength of the GALEX FUV band is nm, and since the spectrum is expected to be quite flat in this wavelength range the observed GALEX FUV is a reasonable proxy for dissociating radiation at nm (see Draine, 1978).

Identify the associated Hi patches on the surfaces of the surrounding GMCs, measure their surface brightnesses on VLA 21cm images, and convert those brightnesses to Hi column densities. Here we need to be cognizant of the fact that, with our relatively high spatial resolution, we may finally be resolving the Hi features, thus becoming sensitive to optical depth effects. Indeed, evidence for the presence of opticallythick Hi features in our M33 data will be presented later.

Measure, and as far as possible deproject, the separation between the central UV source(s) and the Hi patches. This requires that the GALEX and VLAHi images be accurately aligned.

Obtain (or estimate) dusttogas ratios in the gas, ideally local to the candidate PDR.
This procedure allows us then to calculate the total hydrogen volume densities as local spot measurements at the location of candidate PDRs.
3. Data
The Hi observations of M33 used in this paper were provided by David Thilker and Rob Braun (2007, private communication) and are presented and discussed in Heiner et al. (in preparation). We used the PDR method on these data to calculate the total hydrogen volume densities used here. The linear resolution of the Hi data is pc. The local star formation rate at the location of the candidate PDRs is estimated using the farUV luminosity of the parent OB associations from the GALEX data, assuming that the star formation rate is directly proportional to the luminosity at the GALEX farUV wavelength (see Kennicutt, 1998; Madau et al., 1998)). In that case, using the farUV luminosity results in the same power law slope when relating to . In the absence of regionspecific extinction measurements, we have not made any corrections for extinction; a global extinction would merely shift all logluminosities equally, and therefore would not influence the calculated value of the exponent in the star formation law. The volume Schmidt Law correlation we seek is of the form:
(4) 
4. Results
The computed values of and from Heiner et al. (in preparation) are shown in Figure 2 in log coordinates, along with several fits to straight lines. This figure shows a rough correlation, as might be expected from a star formation law of the form of equation 4. The circles correspond to regions inside , whereas the triangles are candidate PDRs outside of . We have made no attempt to fit these regions separately. The straightforward ordinary leastsquares (OLS) fit to the data represented by the dashdotted line in the left panel yields .
4.1. Selection effects
The data plotted in Figure 2 clearly avoid large regions of the graph, as follows:
4.1.1 Low values of
Our survey of star clusters in M33 cuts off below a limit of ergs/sec/Å). Most of these faint clusters (which contain few O stars) are found outside the isophotal radius in regions of low confusion; here the lower limit is set by the sensitivity of the GALEX data. Such faint clusters also exist inside (along with a profusion of widelyscattered FUVproducing B stars), but confusion prevents their accurate tally. However, complete counts are not required for our study, since we subtract away the mean background levels on the FUV and Hi images in order to establish the excess of the cluster and the excess Hi which the cluster’s FUV photons produce on the surfaces of nearby GMCs. In this paper we have therefore not pursued FUVfaint objects inside the main disk of the galaxy^{1}^{1}1A quantitative study of these “backgrounds” is nevertheless highly relevant to the interesting question of what fraction of the total Hi content of the entire galaxy arises through photodissociation, i.e., what fraction of the Hi present is not “primordial”, but rather a product of the star formation process..
4.1.2 High values of
At the bright end of the FUV luminosity we see that there are no OB associations more luminous than NGC 604 at , the equivalent of several dozens of O4 stars (as inferred from Bruhweiler et al., 2003). This limit is apparently set by the physics of “remaining” steps in the star formation process in M33, which we do not consider here (nor even pretend to understand).
4.1.3 Region II
The sloping dashed line (partially) bounding the empty triangular region labelled II in the lower right corner of the left panel of Figure 2 indicates an Hi column density lower limit of , at a typical dusttogas ratio of around 0.33. The available data progressively disappears as we approach this line, which marks the approximate sensitivity limit of our VLAHi data. This is a common (but largely benign) observational selection effect; the computed values of are more noisy here, but they are not strongly biased. A small bias could creep in if the slope of the line reflecting the selection effect were significantly different from the slope of the actual star formation relation, and at the same time there is an abundance of data points near this selection limit. Although the former is likely to be true, the latter is most assuredly not.
4.1.4 Region I
The other empty triangular region labeled I in the upper left corner of the left panel of Figure 2 appears to arise from a more interesting (and also a more sinister) observational selection effect. The sloping dashed boundary line here is defined by an Hi column density upper limit of at a characteristic dusttogas ratio; higher column density estimates are apparently very rare. We suggest that this is a consequence of the 21cm line becoming progressively more optically thick at about this limiting value, thereby underestimating the true Hi column density and hence, following equation 3, overestimating the volume densities . Points can lie somewhat to the left of this boundary owing, for example, to different dusttogas ratios compared to the typical value of 0.25 that we adopted. The different typical dusttogas ratio from the one we used to delineate region II (0.33) is a reflection of the fact that the higher column densities occur in the inner regions of M33, where the metallicity is higher as well. On the other hand, the lower column densities are measured mostly in the outer regions, where the metallicity is systematically lower. This selection effect is an optical depth limit, appearing as a consequence of the improved spatial resolution of the present data set that leads to a (partial) resolution of the M33 GMCs in Hi.
4.2. A simple correction
The appearance of an optical depth limit to the 21cm surface brightness leads to a serious bias in determining the true exponent of the star formation law from the observations. The dashdot OLS fit line in the left panel of Figure 2 clearly shows a slope (of 0.6) which is too low, a consequence of biasing the calculated gas density to values that are artificially high. To explore this selection effect further, we divided the data into six (log) luminosity ranges, and computed the histograms of the values in these six horizontal strips; the results are shown in Figure 3. These distributions appear to be approximately lognormal, but are also clearly biased by the selection effect. As a first attempt to correct for this bias, we identified the value of corresponding to the peak of the histogram in each strip (shaded darker gray in the figure) and fitted an OLS line to these six values. The exponent found by this (oversimplified) method is 1.2 (graphed as the dashed line labeled “Hist” in the right panel of figure 2). This value must still be too small, since we have assumed that the histograms of figure 3 are merely cut off to varying degrees on their left sides, rather than being cut off and skewed towards their right sides in some nonlinear way by the bias. Removing the outer bins does not change the slope.
4.3. A regression analysis
We have shown that selection effects truncate and skew the observed distribution of cloud densities, biasing the value of the exponent obtained by a simple leastsquares fit to the star formation law to values that are too low. These observational limits cause selection effects similar to the wellknown “Malmquist bias”. In order to obtain a more statisticallysound value for the exponent , a regression analysis based on the use of “Monte Carlo Markov Chains” (MCMC) can be used on the cloud density and the UV luminosity data (Equation 3), taking into account both the bias and the (relatively large) uncertainties in our cloud density measurements. This form of Bayesian inference (see e.g. Tolstoy & Saha, 1996), uses a priori assumptions of the data, that are adjusted by the measured data. It provides a statistical formalism to obtain the true exponent based on simulations of unbiased data, where the “simple approach” in our previous section provided an intuitive way to estimate the missing data with similar assumptions. In both cases the assumption is the lognormal distribution of the values of in limited ranges of the UV luminosity.
The full procedure is described in Kelly (2007), where it is shown how this approach (referred to as a maximum likelihood estimate MLE using a Gaussian mixture model) is robust in the presence of censored data and large measurement errors. We used the IDL routines described in that paper to perform our analysis. The result of the MLE is an a posteriori median estimate, which means that the slope (and intersect) resulting from this method are direct measurements of their true values.
The IDL routine constructs a Monte Carlo Markov Chain (Jewell et al., 2009, review the formalism in their Appendices), where new data points are simulated using the “Metropolis  Hastings” algorithm. This particular algorithm is appropriate when errors dominate, and when the variables may not be independent, as is the case here. The regression parameters (generally) converge after several thousand iterations to provide the parameters drawn from the posterior distribution. It is important to realize that the value for the exponent in the star formation law obtained in this way is a random draw from the true distribution of all exponents consistent with the noisy and censored input data. The method can therefore also be used repeatedly to provide an estimate of the stochastic error in the exponent.
To use the MCMC method as presented by Kelly (2007), we need to adopt reasonable estimates for the () errors in and . This is essential for the algorithm to estimate the true distribution of and . The crude but intuitive attempt at correcting for the selection effects described in the previous section shows that the cloud densities (ignoring the truncation and skewing) roughly follow a lognormal distribution. Also, the errors in dominate the correlation between log() and log(). We take absolute uncertainties in the values of log() and log() of 0.5 and 0.2 respectively, conservative values reflecting the relatively large error in and the smaller error in , where is directly based on the measured UV flux^{2}^{2}2Our forthcoming paper on the input data set will have the values of observed UV fluxes.. We then run the regression routines repeatedly. The resulting slopes are shown in Figure 4. The distribution of the values of the slope is nonnormal, but the average value is 1.4 with a range of . Note that the MCMC method is invariant to swapping the coordinate axes. It is also robust against small variations in the estimated uncertainties for and , and setting these uncertainties to nearly 0 brings back the OLS bias. The solid line with slope labeled MCMC in the right panel of Figure 2 reflects the outcome of our Bayesian approach to correct for our selection effects, and provides our best estimate of the exponent in the Schmidt Law of star formation.
5. Discussion and conclusions
The new observational method to estimate the slope of the star formation law that we present here can be further improved in various ways. While some selection effects may have been dealt with acceptably, others remain, and all such effects need to be considered carefully as we have attempted to do here.

The method yields cloud densities with a relatively high uncertainty (we used a conservative estimate of 50%). This uncertainty could be reduced by improving the quality of the data; however, the spread in cloud densities may very well be intrinsic, similar to the scatter in the local metallicities measured at different places in the disks of nearby galaxies (e.g. Rosolowsky & Simon, 2008).

The star formation rate is derived from the UV fluxes of OB associations, and these may be affected by extinction to various degrees which we have ignored here. A detailed extinction map of M33 would be helpful, as only a spatiallyvarying extinction correction would influence our determination of the exponent in the star formation law.

A better calibration of the observed UV flux to the actual star formation rate could be carried out. However, Madau et al. (1998) found that the star formation rate is directly proportional to the FUV luminosity at 150 nm, which is very close to the effective wavelength of the GALEX farUV data, so it is unlikely that further precision here would change our results.

A correction for the effects of high optical depth in the 21cm line could be developed in order to obtain more accurate estimates of the Hi column density. While such a correction would surely help, it is unfortunately not clear how one might determine it.

Lower atomic hydrogen column densities can only be used if they can be distinguished clearly from the general Hi background. More sensitivity and high spatial resolution would be required here, although we have argued that the biasing effect of limited Hi sensitivity on the power law slope is minor.
Unique to our approach is the use of cloud volume densities, not surface densities, and the enhanced sensitivity of the PDR method to lower gas densities (less than a few hundred baryons ) when compared to the higherdensity sensitivity of the CO(10) estimates. To our knowledge, the only exception to this is the work of Abramova & Zasov (2008); these authors specifically aimed to explore a volume density Schmidt Law. However, they could make only rough (CObased) gas volume density estimates, based on column density measurements and an (uncertain) theoretical model of the thickness of galactic disks.
To conclude, we have carried out the first direct determination of the exponent in the star formation law as initially formulated by Schmidt over 50 years ago. The star formation rate has been estimated using the farUV luminosities of OB associations in M33. The total gas volume density (atomic + molecular) has been obtained using a method which regards the atomic hydrogen as the dissociated “skins” of the molecular clouds out of which the young stars have formed. We have used a simple onedimensional slab model for the dissociated layer of Hi on the cloud surfaces. We also show how observational selection leads to a serious bias in the determination of the exponent of the star formation law. We have used a “maximum likelihood” technique to account for this selection bias. Our result for the slope of the star formation law in M33 is . We note that in this paper we have addressed only the question of the value of the exponent in the Schmidt Law, and not the normalization. Application of our method to a larger sample of galaxies will provide additional information on this normalization and its universality.
Finally we point out that the value we have obtained for the exponent in the Schmidt Law for star formation is close to that used in many recent successful numerical simulations of galaxy formation and evolution. This correspondence lends further credence to our view that the Hi in galaxy disks is a product of the star formation process, and that this process can be usefully analysed using the simple physical model described in this paper.
References
 Abramova & Zasov (2008) Abramova, O. V., & Zasov, A. V. 2008, Astronomy Reports, 52, 257
 Allen (2004) Allen, R. J. 2004, in Astrophysics and Space Science Library, Vol. 319, Penetrating Bars Through Masks of Cosmic Dust, ed. D. L. Block, I. Puerari, K. C. Freeman, R. Groess, & E. K. Block, 731
 Allen et al. (1986) Allen, R. J., Atherton, P. D., & Tilanus, R. P. J. 1986, Nature, 319, 296
 Allen et al. (2004) Allen, R. J., Heaton, H. I., & Kaufman, M. J. 2004, ApJ, 608, 314
 Allen et al. (1997) Allen, R. J., Knapen, J. H., Bohlin, R., & Stecher, T. P. 1997, ApJ, 487, 171
 Bruhweiler et al. (2003) Bruhweiler, F. C., Miskey, C. L., & Smith Neubig, M. 2003, AJ, 125, 3082
 Buat et al. (1989) Buat, V., Deharveng, J. M., & Donas, J. 1989, A&A, 223, 42
 Dobbs & Pringle (2009) Dobbs, C. L., & Pringle, J. E. 2009, MNRAS, 396, 1579
 Draine (1978) Draine, B. T. 1978, ApJS, 36, 595
 Heiner et al. (2008a) Heiner, J. S., Allen, R. J., Emonts, B. H. C., & van der Kruit, P. C. 2008a, ApJ, 673, 798
 Heiner et al. (2009) Heiner, J. S., Allen, R. J., & van der Kruit, P. C. 2009, ApJ, 700, 545
 Heiner et al. (2008b) Heiner, J. S., Allen, R. J., Wong, O. I., & van der Kruit, P. C. 2008b, A&A, 489, 533
 Heyer et al. (2004) Heyer, M. H., Corbelli, E., Schneider, S. E., & Young, J. S. 2004, ApJ, 602, 723
 Jewell et al. (2009) Jewell, J. B., Eriksen, H. K., Wandelt, B. D., O’Dwyer, I. J., Huey, G., & Górski, K. M. 2009, ApJ, 697, 258
 Kelly (2007) Kelly, B. C. 2007, ApJ, 665, 1489
 Kennicutt (1997) Kennicutt, R. C. 1997, in Astrophysics and Space Science Library, Vol. 161, Astrophysics and Space Science Library, 171–195
 Kennicutt (1998) Kennicutt, Jr., R. C. 1998, ARA&A, 36, 189
 Madau et al. (1998) Madau, P., Pozzetti, L., & Dickinson, M. 1998, ApJ, 498, 106
 Rosolowsky & Simon (2008) Rosolowsky, E., & Simon, J. D. 2008, ApJ, 675, 1213
 Schmidt (1959) Schmidt, M. 1959, ApJ, 129, 243
 Schmidt (1963) —. 1963, ApJ, 137, 758
 Smith et al. (2000) Smith, D. A., Allen, R. J., Bohlin, R. C., Nicholson, N., & Stecher, T. P. 2000, ApJ, 538, 608
 Sternberg (1988) Sternberg, A. 1988, ApJ, 332, 400
 Tolstoy & Saha (1996) Tolstoy, E., & Saha, A. 1996, ApJ, 462, 672
 Verley et al. (2010) Verley, S., Corbelli, E., Giovanardi, C. & Hunt, L. K. 2010, A&A, 510, 64