text:
# The luminosity functions and stellar masses of galactic disks and spheroids

## 1 Introduction

It has long been known that galaxies in the nearby Universe display a range of morphological characteristics that distinguish between disk-dominated (e.g. spiral) and spheroid-dominated (e.g. elliptical) galaxies, with many fine subdivisions within each class (e.g. Hubble 1926; de Vaucouleurs 1959). The origins of these different types of galaxy and the evolutionary connections, if any, between them are still unclear, although there is a wealth of observational data and several proposed theories.

Traditionally, morphology has been assigned by human classifiers directly from galaxy images, a process that is accurate only to within about two T-types (Naim et al., 1995). More recently, computer-based algorithms have been developed which show a reasonable correlation with “eyeball” estimates, at least for bright galaxies (e.g. Abraham et al. 1996). Unfortunately, all these classifications tend to give considerable weight to detailed morphological features, such as spiral arms or asymmetries in the image, and are difficult to compare with current theoretical predictions which focus on simpler quantities such as the total stellar mass or luminosity in the disk and spheroidal components. Here, we consider a morphological quantifier that is more easily related to theoretical models.

It is widely believed that disks form by slow accretion of gas which acquired angular momentum through tidal torques (Hoyle, 1949; Peebles, 1969), although whether this picture works in detail remains an open question (Navarro, Frenk & White, 1995; Navarro & Steinmetz, 1997, 2000; van den Bosch, Burkert & Swaters, 2001). Spheroids, on the other hand, are thought to form either by a “monolithic collapse” (Eggen, Lyden-Bell & Sandage, 1962; Jimenez et al., 1999), or as a result of mergers of pre-existing galaxies (Toomre 1977; Barnes & Hernquist 1992 and references therein). Detailed theoretical predictions for the statistical morphological properties of galaxies and their evolution have been calculated for the hierarchical merging formation mechanism appropriate to cold dark matter cosmologies (Kauffmann, White & Guiderdoni, 1993; Kauffmann, 1995, 1996; Baugh, Cole & Frenk, 1996a,b; Somerville, Primack & Faber, 2001). Among other things, these models give the relative luminosities and stellar masses of the spheroids and disks of galaxies.

In this work, we measure the I-band bulge-to-total light ratio, B/T, for a large sample of galaxies with spectroscopic redshifts by fitting two-dimensional models to the observed galaxy images. We use this information to estimate the spheroid and disk luminosity functions, as well as the total stellar mass that resides in disks and spheroids. An earlier attempt to estimate the relative contributions to the luminosity from spheroids and disks was carried out by Schechter & Dressler (1987), based on “eyeball” estimates of the bulge-to-total ratio.

The remainder of this paper is arranged as follows. In §2 we describe our basic dataset. In §3 we introduce the method used to fit model images to the data and thereby extract B/T (along with other interesting parameters) and describe how we estimate errors. In §4 we examine the accuracy of our technique and determine how well the B/T ratio can be measured as a function of the apparent magnitude of a galaxy. We then compute luminosity functions and total stellar masses for spheroids and disks. Finally, in §5 we present our conclusions.

## 2 Data

The galaxy sample used in this work is that of Gardner et al. (1997). The reader is referred to that work for a full description of the data. Here we summarize the most important features of the dataset.

Imaging of two fields of total area $`10`$ deg2 was carried out in the B, V and I bands using the T2KA camera on the Kitt Peak National Observatory (KPNO) 0.9-m telescope, resulting in images with $`0.68`$-arcsec per pixel. Exposures of 300s reached $`5\sigma `$ detection depths of B=21.1, V=20.9 and I=19.6 in 10-arcsecond circular apertures. Imaging was also carried out in the K-band using the IRIM camera on the KPNO 1.3-m telescope resulting in 1.96-arcsecond per pixel images, and a $`5\sigma `$ detection depth of K=15.6 in a 10-arcsecond circular aperture. The positions of the fields were chosen randomly (the field centers are RA 14h15m, Dec. $`+00^{}`$ and RA 18h0m, Dec. $`+66^{}`$). The I-band images, which we will use in this work, were bias-subtracted, flattened using twilight flats and with median sky flats. Objects were identified with the sextractor program (Bertin & Arnouts 1996), using a $`3\sigma `$ threshold. The seeing in the optical images varied in the range $`1.3<\mathrm{FWHM}<2.0`$ arcsec. One field contained a nearby rich galaxy cluster.

Spectroscopic follow-up was obtained for a K-selected sample in sub-regions of total area $`4.4`$ deg2, using the Autofib-2 fiber positioner and WYFFOS spectrograph on the 4.2m William Herschel Telescope on La Palma. Spectra were obtained for 567 galaxies with K$`<15`$, which allowed redshifts to be measured for 510 galaxies (a redshift completeness of 90%). Although the spectroscopic sample is K-selected, this does not introduce any incompleteness in the I$`<16`$ sample used extensively in this work (i.e. there are no galaxies with I$``$K$`<1`$ in the sample). We also briefly consider an I$`<18`$ sample, for which the spectroscopic completeness falls to around 50% because of the K-band selection (which also introduces a bias in this fainter sample against objects that are blue in I$``$K). For the I$`<16`$ and I$`<18`$ samples the median redshift is $`z=0.08`$ and $`0.14`$ respectively.

## 3 Method

Wadadekar, Robbason & Kembhavi (1999) have proposed a two-dimensional galaxy decomposition technique that can efficiently recover B/T ratios (and other parameters) of model galaxy images with high accuracy (see also Byun & Freeman 1995; de Jong 1996). They present a detailed study of the effects of uncertainties in the point spread function (PSF), the presence of nearby stars and the stability of the B/T estimates as a function of signal-to-noise. Our approach is similar to theirs, but we apply the technique to a large photometric sample of real galaxies111Wadadekar, Robbason & Kembhavi (1999) applied their technique to three galaxies for which previous estimates of B/T were available and found reasonable agreement.. We consider the “real-world” problems of automated masking of nearby galaxies and stars and make a thorough assessment of the errors in the measured parameters. We also present, in an Appendix, estimators for the luminosity functions of disks and spheroids.

For this analysis, we have used the sample of 636 $`\mathrm{I}<18`$ galaxies of Gardner et al. (1997) imaged in B, V, I and K, as described in §3. Postage stamp images of $`33\times 33`$ pixels ($`0.68^{^{\prime \prime }}`$/pix) around each galaxy were extracted from the I-band data. This was the best observed band by Gardner et al. (1997) and is particularly well-suited for our purposes because it minimizes the effects of young blue stars. This image size is large enough to include the entire region of the galaxy for which reasonable signal-to-noise is achieved (and in most cases extends well beyond it.) To determine the bulge-to-total ratio, B/T, we fit the two-dimensional surface brightness profile of each galaxy using a combination of an exponential disk,

$$\mathrm{\Sigma }_\mathrm{d}(\theta )=\mathrm{\Sigma }_{\mathrm{d},0}\mathrm{exp}(\theta /\theta _\mathrm{d}),$$
(1)
and an $`r^{1/4}`$-law spheroid,

$$\mathrm{\Sigma }_\mathrm{s}(\theta )=\mathrm{\Sigma }_{\mathrm{s},\mathrm{e}}\mathrm{exp}(7.67[(\theta /\theta _\mathrm{e})^{1/4}1]),$$
(2)
where $`\theta `$ is the angular distance from the galaxy center. We will also consider a more general $`r^{1/n}`$ spheroid profile, as Wadadekar, Robbason & Kembhavi (1999) did. The disk is allowed to be inclined, and to have arbitrary position angle. The spheroid is allowed an ellipticity (defined as the ratio of semi-major to semi-minor axes) in the range 1 to 6, and can also have an arbitrary position angle. To mimic seeing, we construct mock images using these profiles which we then smooth with a Gaussian filter (integrated over each pixel to account for the variation of the PSF across the pixel),

$$p(\theta )=\mathrm{exp}[(\theta /\sigma )^2/2]/2\pi \sigma ^2.$$
(3)
The width of the Gaussian is treated as a free parameter, to account for variations in seeing between the images. (We will examine in §4.1 the effect of using a more realistic PSF.) The postage stamp images were centered on the galaxy of interest, but we allow the position of the image center to vary since in many cases the resulting sub-pixel variations lead to lower values of $`\chi ^2`$. We also allow a small contribution from a faint, constant surface brightness background in order to take into account small inaccuracies in sky subtraction. The best fitting parameters for each galaxy were then obtained by minimizing $`\chi ^2`$ using Powell’s algorithm (Brent, 1973). There are a total of 12 fit parameters (13 if we include $`n`$ when fitting $`r^{1/n}`$ spheroid profiles) summarized in Table 1. In §4.1 we will consider how accurately this procedure recovers the B/T ratio of the galaxies.

A significant fraction of the postage stamp images was contaminated by a secondary galaxy (and occasionally by more than one). We use a simple algorithm to identify such contaminants and mask them from the image. The aim is to remove objects which are physically distinct from the galaxy of interest without masking any pixels of the galaxy itself. We first rank the pixels in the image by surface brightness and then proceed to find groups of bright pixels. The brightest pixel is assigned to the first group. Successive pixels are assigned to a pre-existing group if they touch it (i.e. if they are adjacent either horizontally, vertically or diagonally), or else are assigned to a new group. In the case where a pixel touches more than one group, the two groups are merged. This process is continued until pixels of a fixed signal-to-noise are reached (specifically, we consider only pixels more than $`3\sigma `$ above the sky background). If more than one group exists at this point, the group at the center is deemed to be the galaxy of interest and the pixels of all secondary groups are marked as being contaminated and are not included in the $`\chi ^2`$ sum. This simple algorithm works well in the majority of cases, but it fails in a few (28 out of 636 galaxies), either by not removing a contaminating galaxy, or by removing a significant fraction of the primary galaxy. Rather than attempting to use a more complex algorithm in these cases, we resorted to cleaning the image by hand (i.e. we view the image and manually mark the contaminated pixels).

The majority of galaxies (363 out of 626) are fit reasonably well by this procedure. We regard a galaxy as being reasonably well fit when $`Q`$, the probability that the measured value of $`\chi ^2`$ is exceeded by random fluctuations, is greater than 5%. Not surprisingly though, many galaxies are not well fit. These typically show signs of strong morphological disturbance (perhaps due to a recent or imminent merger) or other inhomogeneities. Examples of well-fit and poorly-fit galaxies are given in Figure 1. Poorly-fit galaxies are easily identified by their large $`\chi ^2`$ values and so may be excluded from further analysis if desired. It should be noted that a poor-fit does not indicate a failure of our fitting procedure *per se*, rather it signals that the galaxy is not well described by a combination of a spheroid and a disk. We choose to show results computed using the entire sample, regardless of how well a galaxy was fit, but we will comment on how our results change if badly fit galaxies are excluded from the analysis.

Errors on the fitted parameters could, in principle, be determined using a $`\mathrm{\Delta }\chi ^2`$ approach, but this would require mapping $`\chi ^2`$ in the 12 dimensional parameter space of the fit — an exceedingly time consuming exercise — and, in any case, the errors are unlikely to be normally distributed given that the model is highly nonlinear in the parameters. We therefore adopt a Monte Carlo approach to error estimation. Using the best fitting model for each galaxy, we generate 30 realizations of that model, add random noise at the same level as in the real image, and mask out any pixels which were masked out in the original. We then find the best fitting parameters for each realization and take their distribution as indicative of the uncertainties in the actual fit. It should be noted that this is only a valid procedure if the original image is well fit by the model. In cases where this is not the case, there is no reason to expect the Monte Carlo distributions to give an estimate of the true errors. The B/T distributions for the three galaxies illustrated in Figure 1 are shown in that figure.

## 4 Results

### 4.1 Accuracy Checks

We begin by assessing the reliability of our procedure for recovering the true B/T ratio of a galaxy (assuming, of course, that real galaxies are well described by our model). Our Monte Carlo procedure for error estimation allows a determination of the accuracy of our technique. For each galaxy, the value of B/T input into the Monte Carlo simulations may be compared to the mean and standard deviation of the distribution of 30 recovered B/T values. Figure 2 gives the results of these accuracy tests. The left-hand panel shows the standard deviation of the recovered B/T ratio, as a function of the I-band apparent magnitude of the mock image. For bright galaxies ($`m_\mathrm{I}<16`$), $`\sigma _{\mathrm{MC}}`$ is fairly small, typically less than about 0.1. However, for fainter galaxies, $`\sigma _{\mathrm{MC}}`$ increases very rapidly, resulting in rather poorly constrained B/T values. (In reality, $`\sigma _{\mathrm{MC}}`$ depends also upon the other parameters that describe the mock image, but the correlation with apparent magnitude is the most important.)

In the right-hand panel of Figure 2, we plot the mean value of B/T recovered from the Monte Carlo simulations against the true value for the mock image. The large solid circles indicate those images for which $`\sigma _{\mathrm{MC}}0.1`$. Evidently, for these galaxies the value of B/T is recovered accurately and without any strong systematic bias. The small dots show the results for all other galaxies. Now the scatter is much larger and, more importantly, there are systematic biases in the mean recovered B/T, such that very low and very high values are avoided. This effect is not surprising: the values of B/T for these faint galaxies are almost entirely unconstrained. (Note that the standard deviation for a completely uniform distribution of B/T is approximately 0.3.) As a result, the distribution of B/T from the Monte Carlo simulations becomes close to uniform, with the mean tending towards 0.5 as $`\sigma _{\mathrm{MC}}`$ increases. For the 90 galaxies in our sample with $`m_\mathrm{I}16`$, there is a tight correlation between $`(B/T)_{\mathrm{true}}`$ and the mean value recovered from the Monte Carlo simulations and it is this subsample that we will use below to compute luminosity functions. Unfortunately, its relatively small size limits the statistical accuracy of our estimates quite considerably.

To find the best fit solution we must choose initial values for the parameters to be fit and then use Powell’s algorithm to search for values producing a better fit. For parameters such as the position angles, disk inclination and spheroid/disk sizes we make initial guesses based on the image being fitted. Other parameters are initially assigned “typical” values. We have checked the effect of altering these initial values. For galaxies with $`\sigma _{\mathrm{MC}}0.1`$, the choice of initial value makes almost no difference to the recovered values of the parameters, indicating that our technique is finding the true minimum $`\chi ^2`$. As $`\sigma _{\mathrm{MC}}`$ becomes larger, however, the recovered parameters begin to depend strongly upon the initial values chosen. For these images, the $`\chi ^2`$ surface in the 12 dimensional parameter space does not possess an obvious minimum (i.e. it is very noisy). This is just another way of saying that the values of the fitted parameters for these faint images are highly uncertain.

Andredakis, Peletier & Balcells (1995) have demonstrated that the bulges of spiral galaxies are more accurately fit by an $`r^{1/n}`$, rather than by the more usual $`r^{1/4}`$, surface brightness profile, with values of $`n`$ ranging from around 1 to 6 (similar variations in $`n`$ are seen for elliptical galaxies; Bingelli & Cameron 1991; Caon, Capaccioli & D’Onofrio 1993). They show that the value of $`n`$ is strongly correlated with morphological type. Although our present data provide only rather poor constraints on the value of $`n`$ (due to limited angular resolution and signal-to-noise), we have nevertheless repeated the fitting procedure using $`r^{1/n}`$ profiles for the spheroids, treating $`n`$ as a free parameter. For galaxies where the B/T ratio is well determined, we find that there is a strong correlation between the B/T values obtained with $`r^{1/4}`$ and $`r^{1/n}`$ profiles, although inevitably some scatter is present. The disk and spheroid luminosity functions computed using B/T ratios from $`r^{1/n}`$ fits show no statistically significant difference from those using $`r^{1/4}`$ fits.

Finally, we remind the reader that our analysis makes use of a Gaussian PSF to mimic the effects of seeing in the data. A Gaussian accurately describes the core of the PSF measured from bright stars in the images. However, a profile consisting of a Gaussian core plus power-law wings provides a better match to many of the stellar profiles. (The variation of Gaussian and core components from night to night in the imaging data is not so well characterized, however, and this is why we make use of a simple Gaussian for our main analysis.) Fitting the images using such a profile (keeping the relative proportions of Gaussian and power-law wings fixed, but allowing the overall radial scale of the PSF to be a free parameter) results in small changes in the B/T ratio, typically significantly smaller than the error in the best fit value. Thus, the luminosity functions presented below are unaffected by the exact choice of PSF. However, it is clear that a good characterization of the PSF and its variation will be crucial to obtain accurate disk and spheroid luminosity functions from larger, higher quality datasets.

### 4.2 Luminosity Functions

Using the I$`<16`$ sample of approximately 90 galaxies for which we have good estimates of the B/T ratio we now proceed to estimate the disk and spheroid luminosity functions. Our aim here is to develop the techniques required for this measurement and demonstrate them using a particular dataset. Given the small size of the dataset we must expect that both statistical (due to the small number of galaxies) and systematic (due, for example, to the lack of rich clusters in the dataset) errors will be present. These issues are considered further in §5.

To determine the present-day luminosity functions, we need to apply k$`+`$e corrections to the galaxy luminosities. We use the type-dependent k$`+`$e corrections obtained by Gardner et al. (1997). Briefly, a set of model galaxy colors was computed using an updated version of the Bruzual & Charlot (1993) stellar population models with a range of star formation histories. The observed colors of each galaxy were matched to one of the models and that particular model was then used to extrapolate the observed galaxy luminosity to $`z=0`$. Note that our type-dependent k$`+`$e corrections are based on the *total* (i.e. disk plus spheroid) color of each galaxy. In principle, k$`+`$e corrections could be applied to each component separately if spheroid/disk decompositions were carried out in several bands. Given the uncertainties in our present estimates of B/T, we refrain from this degree of complexity in this analysis.

We use the stepwise maximum likelihood (SWML) estimator proposed by (Efstathiou, Ellis & Peterson, 1988, hereafter EEP) and also the parametric maximum likelihood method proposed by (Sandage, Tammann & Yahil, 1979, hereafter STY) to compute disk and spheroid luminosity functions. The detectability of a spheroid depends on both its apparent magnitude and the B/T ratio, and we must account for this in constructing the likelihood function. This leads us to define a two-dimensional function, $`\mathrm{\Phi }(M,B)`$, such that $`\mathrm{\Phi }(M,B)\mathrm{d}M\mathrm{d}B`$ is the number of galaxies per unit volume with B/T ratio in the range $`B`$ to $`B+\mathrm{d}B`$ and *spheroid* absolute magnitude in the range $`M`$ to $`M+\mathrm{d}M`$ (with an equivalent definition for disks). The application of the maximum likelihood estimator to this function is discussed in detail in Appendix A. The normal luminosity function of spheroids is readily derived using $`\varphi (M)=_0^1\mathrm{\Phi }(M,B)dB`$ (and similarly for disks). For the STY method we must assume some parametric form for the luminosity function. We have tried fitting the disk and spheroid luminosity functions with a “Schechter$``$exponential” form, namely $`\mathrm{\Phi }(M,B)=\varphi (M)\mathrm{exp}(\beta B)`$, where $`\varphi (M)`$ is the normal Schechter function and $`\beta `$ is a parameter to be fit, motivated by the shape of the SWML estimate of these luminosity functions.

Figure 3 shows the resulting I-band luminosity functions with distances computed assuming $`(\mathrm{\Omega }_0,\mathrm{\Lambda }_0)=(0.3,0.7)`$222Assuming $`(\mathrm{\Omega }_0,\mathrm{\Lambda }_0)=(1,0)`$ instead changes our results only slightly, shifting the luminosity function faintwards due to the smaller luminosity distance in this model. and $`H_0=100h`$ km/s/Mpc. Triangles show the total luminosity function; circles and stars show the spheroid and disk luminosity functions separately. These SWML luminosity functions are normalized to the I-band number counts in the Sloan Digital Sky Survey using the procedure described in the Appendix. For the total luminosity function, we plot the standard SWML errorbars (obtained from the covariance matrix of the luminosity function as described by EEP), but for the spheroid and disk luminosity functions, the errorbars are the sum in quadrature of the standard SWML errors and the variance in the luminosity function estimated from the 30 Monte Carlo realizations of the spheroid/disk decomposition process. The errors from each source are of comparable magnitude (although the variance from the Monte Carlo realizations is the smaller of the two). The very heavy solid line shows the best-fitting Schechter function to the total luminosity function (determined using the STY method); the inset shows the values of $`\alpha `$ and $`M_{}`$ for this fit, together with their 1 and $`2\sigma `$ error contours. (The small sample size is reflected in rather large and correlated uncertainties in $`M_{}`$ and $`\alpha `$.) Heavy and thin solid lines show the best fit STY “Schechter$``$exponential” luminosity function fits to the spheroid and disk luminosity functions respectively (with the corresponding confidence ellipses for $`\alpha `$ and $`M_{}`$ shown in the inset, and the values of $`\beta `$ given in the figure caption). A likelihood ratio test (Efstathiou, Ellis & Peterson, 1988) shows that the Schechter$``$exponential luminosity function is not a particularly good fit to the data. With the present small dataset we have been unable to find a better functional form. This situation will be rectified with a larger dataset (assuming that some suitable functional form does actually exist).

The I-band luminosity functions of disks and spheroids are remarkably similar. The only significant difference is that the spheroid luminosity function is somewhat lower at faint magnitudes. However, given the small size of the present sample, this difference may not be robust. The luminosity densities in disks and spheroids obtained by integrating the SWML luminosity functions over the range of absolute magnitudes shown in Fig. 3 are $`5.8\pm 0.8`$ and $`4.7\pm 0.7\times 10^7hL_{}/\text{Mpc}^3`$ respectively (where we have taken $`M_{}=4.14`$ in the I-band; Cox 2000). In principle, we can use our Schechter$``$exponential fits to estimate the total luminosity density, extrapolating to include the contribution from arbitrarily faint spheroids and disks. Doing so yields results which agree with the SWML estimates within the quoted errors, suggesting that our determination may have suitably converged. However, it must be kept in mind that the Schechter$``$exponential form is not a particularly good fit to the current datasets.

In Fig. 4, we show slices through the bi-variate luminosity function, $`\mathrm{\Phi }(M,B)`$, at constant $`M`$ for several values of $`M`$ (as indicated in the figure, and in bins of width $`\mathrm{\Delta }M=0.48`$). It is evident that $`\mathrm{\Phi }(M,B)`$ is *not* independent of $`B`$ and that, in fact, it may not be separable into a simpler form $`\mathrm{\Phi }(M,B)=\varphi (M)g(B)`$. This is particularly noticeable for the spheroid luminosity function. Figure 4 shows that spheroid-dominated systems (i.e. B/T$`>2/3`$) are found only in the brightest spheroids, while disk-dominated systems (i.e. D/T$`>2/3`$) have disks with a much broader range of luminosities. This point is made more clearly in Fig. 5 where we show the spheroid luminosity function of spheroid-dominated systems and the disk luminosity function of disk-dominated systems. We find spheroid-dominated systems in abundance only brightwards of $`M_\mathrm{I}5\mathrm{log}h21`$, but disk-dominated systems across the whole range of luminosities.

### 4.3 The stellar mass in disks and spheroids

The stellar mass associated with the luminosity of each galaxy is easily obtained by a similar procedure to that employed to calculate the k-corrections (namely fitting their BVIK colors to a set of template galaxies). We use these estimates to construct the stellar mass functions of spheroids and disks in the local Universe and thereby estimate the total mass content in each component. The method that we adopt is the same that Cole et al. (2001) used to measure the total stellar mass density in the Universe from a combination of 2dFGRS and 2MASS data. Cole et al. (2001) found the stellar mass density333Specifically, Cole et al. (2001) estimated the mass locked up in stars and stellar remnants which differs from the time integral of the star formation rate due to recycling of material by massive stars. We adopt the same definition of stellar mass here. in units of the critical density to be $`\mathrm{\Omega }_{\mathrm{stars}}=0.0016\pm 0.00024h^1`$ or $`\mathrm{\Omega }_{\mathrm{stars}}=0.0029\pm 0.00043h^1`$, depending on whether a Kennicutt (1983) or a Salpeter (1955) stellar initial mass function (IMF) was assumed. (Both estimates include the effects of dust on galaxy luminosities as described by Cole et al. 2001.) We normalize our stellar mass functions by requiring them to produce the same number density of galaxies more massive than $`M_{}`$ as the Cole et al. (2001) stellar mass function. The total stellar mass density inferred from our $`\mathrm{I}<18`$ sample is then $`\mathrm{\Omega }_{\mathrm{stars}}=0.0009\pm 0.00007h^1`$ or $`\mathrm{\Omega }_{\mathrm{stars}}=0.0017\pm 0.00012h^1`$ for the same two IMFs respectively and the same prescription for dust-extinction. Errors on the stellar mass density were found by summing in quadrature the error from each individual bin in the SWML mass function, together with the error in the overall normalization. Our estimates include contributions from galaxies with stellar masses greater than $`10^9h^2M_{}`$ below which the SWML stellar mass function is not well determined. We can check our result using the STY stellar mass function. For this, we fit a Schechter function convolved with a Gaussian of width $`0.1`$ in $`\mathrm{log}_{10}M_{\mathrm{stars}}`$ to account for the scatter in the relation between stellar mass and I-band absolute magnitude. We find that our estimates of $`\mathrm{\Omega }_{\mathrm{stars}}`$ using the SWML and STY mass functions agree within the errors, suggesting that the result has already converged to sufficient accuracy. Our estimates, however, are lower than those of Cole et al. (2001), a reflection of the flatter faint end slope of our mass functions which may well be due to the small size of our sample.

For our $`\mathrm{I}<16`$ sample, for which the bulge-to-total ratio is well measured, we find a slightly higher total stellar mass density of $`\mathrm{\Omega }_{\mathrm{stars}}=0.0012\pm 0.00014h^1`$ for the Kennicutt (1983) IMF. Splitting into spheroidal and disk components, we find $`\mathrm{\Omega }_{\mathrm{stars},\mathrm{spheroids}}=0.00039\pm 0.00006h^1`$ and $`\mathrm{\Omega }_{\mathrm{stars},\mathrm{disks}}=0.00051\pm 0.00008h^1`$ for this same IMF. (Note that with the SWML method the stellar mass densities of disks and spheroids are not guaranteed to sum to give the total stellar mass density.) If, instead, we assume a Salpeter IMF, the ratio of disk to spheroid stellar mass densities increases slightly from 1.31 to 1.37, but this change is negligible given the current errors in these quantities. Although the small size of our sample is clearly a significant limitation, this initial result suggests that spheroids and disks contribute about equally to the stellar mass density of the Universe. The techniques developed in this paper, when applied to a much larger dataset, should allow their contributions to be more accurately determined.

## 5 Discussion

We have presented a detailed method to determine the bulge-to-total ratios of galaxies by fitting to two-dimensional photometry, and have applied this technique to determine the I-band bulge-to-total luminosity ratios of a sample of approximately 600 galaxies brighter than I$`=18`$ with spectroscopic redshifts. Our approach is designed to work with realistic galaxy images, dealing automatically with contamination by nearby objects, a varying PSF and small changes in the background from image to image. A crucial part of the fitting procedure is a Monte Carlo determination of the errors on the fitted parameters, an approach which is favored since it is fast and automatically accounts for the highly non-linear nature of the model parameters. For the current sample of galaxies, around 60% are well fit by a combination of an exponential disk and an $`r^{1/4}`$-law spheroid. Those that are not well fit frequently show signs of morphological disturbance. We find that bulge-to-total ratios are determined accurately (i.e. with errors of around $`10\%`$) only for galaxies brighter than I$`16`$.

For the 90 galaxies brighter than I$`=16`$ in this sample we measure the B/T ratio with reasonable accuracy. We have used the resulting disk/spheroid decomposition of these bright galaxies to construct separate luminosity functions for disks and spheroids. We find no significant differences between them when considered purely in terms of luminosity, although the statistical uncertainties associated with the small sample size make the detection of any differences difficult. However, when we consider the bi-variate distributions of luminosity and bulge-to-total or disk-to-total light, we find that spheroid dominated systems (B/T$`>2/3`$) only occur for the brightest spheroids, while disk-dominated systems (D/T$`>2/3`$) occur for a much broader range of disk luminosities.

The relative contributions of disks and spheroids to the total stellar mass density in the Universe is a very important constraint on theories of galaxy formation which attempt to describe the assembly of galaxies as a function of time. We find, perhaps surprisingly, that the disks and spheroids in our sample contribute almost equally to the stellar mass density today (in a ratio of $`1.3\pm 0.2`$). Since the stellar populations in disks are generally younger than those in spheroids, it is an interesting coincidence that the total stellar mass in the two kinds of structural components should be so similar at the present time.

Schechter & Dressler (1987) reached a similar conclusion to ours using a photometric comparison technique to estimate the B-band bulge-to-disk ratios of galaxies. While this technique may not be as accurate as our own on a galaxy-by-galaxy basis, it should provide a good estimate of the total contribution of each component to the stellar mass density. It is therefore reassuring that our results agree well with those of Schechter & Dressler (1987). A different result was obtained by Fukugita, Hogan & Peebles (1998), who derived a ratio of disk to spheroid stellar mass density of $`0.33\pm 0.23`$. Although they found comparable B-band luminosity density in spheroids and disks, they adopted a spheroid mass-to-light ratio around four times greater than that for disks, resulting in spheroids making a significantly greater contribution to the stellar mass density. While we use a more accurate technique for converting from luminosity to stellar mass (a technique which could be improved further if B/T ratios were measured for each galaxy in several bands), the small size of our sample limits the accuracy of our results. In particular, our sample may not contain enough rich clusters which are known to contain higher fractions of spheroid dominated galaxies than the field (e.g. Dressler 1980) and this could introduce a small bias in our results.

Clearly the greatest limitation of this work is the small size of the sample of galaxies for which accurate disk/spheroid decompositions can be performed. Fortunately, this problem should be remedied in the near future with the advent of high quality, large area photometric surveys such as that being carried out by the SDSS project.

## Acknowledgments

We thank Jon Gardner and Carlton Baugh for supplying data used in this work and for valuable discussions, and the referee, Alan Dressler, for valuable suggestions. We also thank Istvan Szapudi for his assistance in the early stages of this work.

## Appendix A. Estimators for Spheroid and Disk Luminosity Functions

The traditional $`1/V_{\mathrm{max}}`$ estimator is trivially adapted to the case of disk and spheroid luminosity functions. The estimator is applied just as in the case of the standard luminosity function, except that the *total* luminosity of the galaxy (i.e. disk plus spheroid luminosity) is used to compute $`V_{\mathrm{max}}`$, since it is this total luminosity that determines the volume within which the galaxy could have been detected.

The maximum likelihood estimator of Efstathiou, Ellis & Peterson (1988, hereafter EEP) is also easily generalized to the case of spheroid and disk luminosity functions. Consider the case of the spheroid luminosity function (the same arguments apply to disks). As noted earlier, the detectability of a spheroid depends upon both its absolute magnitude, $`M`$, and on the bulge-to-total ratio which we denote by $`B`$ in this Appendix. We begin therefore by defining a two-dimensional function, $`\mathrm{\Phi }(M,B)`$, such that $`\mathrm{\Phi }(M,B)\mathrm{d}M\mathrm{d}B`$ is the number of galaxies with bulge-to-total ratio in the range $`B`$ to $`B+\mathrm{d}B`$ and *spheroid* absolute magnitude $`M`$ to $`M+\mathrm{d}M`$ per unit volume. The normal luminosity function of spheroids is easily recovered using $`\varphi (M)=_0^1\mathrm{\Phi }(M,B)dB`$. The probability that galaxy $`i`$ with spheroid magnitude $`M_i`$ and bulge-to-total ratio $`B_i`$ is seen in a magnitude limited survey is

$$p_i\mathrm{\Phi }(M_i,B_i)/_0^1_{\mathrm{}}^{M_{\mathrm{lim}}^{}(z_i,B)}\mathrm{\Phi }(M,B)dMdB,$$
(4)
where $`M_{\mathrm{lim}}^{}(z_i,B)=M_{\mathrm{lim}}(z_i)2.5\mathrm{log}_{10}B`$ and $`M_{\mathrm{lim}}(z_i)`$ is the limiting absolute magnitude of the survey at redshift $`z_i`$. The use of $`M_{\mathrm{lim}}^{}`$ is necessary since arbitrarily faint spheroids will make it into the survey provided that they have a sufficiently low bulge-to-total ratio (corresponding to sufficiently bright disks).

From this definition we can construct the usual likelihood function

$$\mathrm{ln}=\underset{i=1}{\overset{N}{}}\mathrm{ln}\mathrm{\Phi }(M_i,B_i)\underset{i=1}{\overset{N}{}}\mathrm{ln}\left\{_0^1_{\mathrm{}}^{M_{\mathrm{lim}}^{}(z_i,B)}\mathrm{\Phi }(M,B)dMdB\right\}+\text{const},$$
(5)
where $`N`$ is the total number of galaxies. There are now two ways to proceed. In the first we assume a simple parametric form for $`\mathrm{\Phi }(M,B)`$ and maximize the likelihood with respect to the parameters. This is analogous to fitting a Schechter function (Schechter, 1976) to the normal luminosity function (e.g. Sandage, Tammann & Yahil 1979). A simple parametric form which we have tried to fit our data is

$$\mathrm{\Phi }(M,B)=\varphi (M)\mathrm{exp}(\beta B),$$
(6)
where $`\varphi (M)`$ is the usual Schechter function and $`\beta `$ is a parameter to be estimated from the fit. With this method, the likelihood function of eqn. (5) can be evaluated for each value of the three parameters $`\alpha `$, $`\beta `$ and $`M_{}`$, and hence the parameter values which maximize the likelihood are readily obtained.

The second approach involves splitting $`\mathrm{\Phi }(M,B)`$ into bins in $`M`$ and $`B`$ and treating each as a parameter. This is equivalent to the SWML method of EEP for estimating the standard luminosity function.

We represent $`\mathrm{\Phi }(M,B)`$ as follows:

$$\mathrm{\Phi }(M,B)=\mathrm{\Phi }_{k,h},\{\begin{array}{cc}M_k\mathrm{\Delta }M/2<M<M_k+\mathrm{\Delta }M/2,\hfill & k=1,\mathrm{},N_p\hfill \\ B_h\mathrm{\Delta }B/2<B<B_h+\mathrm{\Delta }B/2,\hfill & h=1,\mathrm{},N_q.\hfill \end{array}$$
(7)
The likelihood function may then be written as

$$\mathrm{ln}=\underset{i=1}{\overset{N}{}}W(M_iM_k,B_iB_h)\mathrm{ln}\mathrm{\Phi }_{k,h}\underset{i=1}{\overset{N}{}}\mathrm{ln}\left\{\underset{h=1}{\overset{N_q}{}}\underset{k=1}{\overset{N_p}{}}\mathrm{\Phi }_{k,h}\mathrm{\Delta }M\mathrm{\Delta }BH[M_k,B_h,M_{\mathrm{lim}}(z_i)]\right\}+\text{const},$$
(8)
where

$$W_{k,h}(M_i,B_i)=\{\begin{array}{cc}1\hfill & \text{if }M_k\mathrm{\Delta }M/2<M_i<M_k+\mathrm{\Delta }M/2,\text{ and }B_h\mathrm{\Delta }B/2<B_i<B_h+\mathrm{\Delta }B/2,\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}$$
(9)
and

$$H[M_k,B_h,M_{\mathrm{lim}}(z_i)]=\frac{1}{\mathrm{\Delta }M\mathrm{\Delta }B}_{B_h\mathrm{\Delta }B/2}^{B_h+\mathrm{\Delta }B/2}_{M_k\mathrm{\Delta }M/2}^{M_k+\mathrm{\Delta }M/2}Q(M,B)dMdB,$$
(10)
where $`Q(M,B)=0`$ if $`M>M_{\mathrm{lim}}(z_i)2.5\mathrm{log}_{10}B`$ and $`Q(M,B)=1`$ otherwise. Since only the shape of the luminosity function is constrained by the above likelihood function, we introduce an additional constraint, $`g=_k_h\mathrm{\Phi }_{k,h}(L_{k,h}/L_\mathrm{f})^\beta \mathrm{\Delta }M\mathrm{\Delta }B1=0`$, where $`L_{k,h}`$ is the total luminosity of a galaxy with spheroid magnitude $`M_k`$ and bulge-to-total ratio $`B_h`$ and $`L_\mathrm{f}`$ is a fiducial luminosity (which we will take to be that corresponding to $`M_\mathrm{I}5\mathrm{log}h=20.5`$), using a Lagrangian multiplier $`\lambda `$ as did EEP. Maximizing $`\mathrm{ln}^{}=\mathrm{ln}+\lambda g`$ then yields

$$\mathrm{\Phi }_{k,h}=\frac{\underset{i=1}{\overset{N}{}}W(M_iM_k,B_iB_k)}{_{i=1}^NH[M_k,B_h,M_{\mathrm{lim}}(z_i)]/_{l=1}^{N_p}_{m=1}^{N_q}\mathrm{\Phi }_{l,m}H[M_l,B_m,M_{\mathrm{lim}}(z_i)]}$$
(11)
which are easily solved with an iterative procedure. The covariance matrix for the parameters is obtained in a manner entirely analogous to that outlined by EEP.

Normalization of the maximum likelihood luminosity function can be achieved using the actual redshift data as described by Loveday et al. (1992), but using $`M_{\mathrm{lim}}^{}`$ in the selection function to account for the effects of the bulge-to-total ratio. A better approach is to normalize by performing a least squares fit to the number counts of galaxies from a wide area survey. The cumulative number count to apparent magnitude $`m`$ is given by

$$n(m)=_0^{\mathrm{}}_0^1_{\mathrm{}}^{M_{\mathrm{lim}}^{}}\mathrm{\Phi }(mD(z)K(z)2.5\mathrm{log}_{10}B,B)\frac{\mathrm{d}V}{\mathrm{d}z}dMdBdz,$$
(12)
where $`D(z)`$ and $`K(z)`$ are the distance modulus and k$`+`$e correction respectively at redshift $`z`$, which we can compute from the SWML estimate of $`\mathrm{\Phi }`$.


format: text

text:
# Extended Near-Infrared Emission from Candidate Protostars in the Taurus-Auriga Molecular Cloud

## 1 INTRODUCTION

Standard models for the formation of a single, low-mass star begin with a slowly rotating cloud of gas and dust with a centrally concentrated density distribution (e.g., Larson 1969, 1972; Stahler, Shu, & Taam 1980, 1981; Shu, Adams, & Lizano, 1987; Shu et al., 1993; Safier, McKee, & Stahler, 1997; Ciolek & Basu, 2000; Whitworth & Ward-Thompson, 2001). When the inward gravitational pressure exceeds the outward thermal and magnetic pressures, the cloud collapses. Material near the axis of rotation with low specific angular momentum contracts into a central star-like core. Once the radius of the centrifugal barrier exceeds the radius of the central stellar core, material away from the rotation axis with large specific angular momentum follows a curved infall trajectory, falls onto the equatorial plane, and forms a circumstellar disk (Cassen & Moosman, 1981; Terebey et al., 1984; Lin & Pringle, 1990; Yorke et al., 1993; Stahler et al., 1994). As the collapse proceeds, the protostellar core becomes more massive and more luminous; the surrounding infalling envelope becomes less massive and less opaque at optical and infrared wavelengths (Lada, 1991). Eventually, a wind from the star-disk system begins to clear gas and dust away from the rotational axis, resulting in the production of collimated jets and a bipolar molecular outflow (Adams & Lin, 1993; Konigl & Ruden, 1993; Konigl & Pudritz, 2000; Shu et al., 2000). Once the infalling envelope disperses or falls into the circumstellar disk, the central star-disk system becomes optically visible as a pre-main-sequence star. As the material in the disk accretes onto the the central star, dissipates, or condenses into planets, the star moves onto the main sequence (Stahler & Walter, 1993; Strom et al., 1993; Cameron, 1995).

This picture of star formation generally accounts for an apparent evolutionary sequence in the observed properties of young stars in nearby molecular clouds (Lada, 1984; Adams et al., 1987; Lada, 1987). Lada (1991) describes three of the classes in this sequence. Class I sources have spectral energy distributions (SEDs) which are broader than a single temperature blackbody function and rise from $`2\mu `$m to $`25100\mu `$m. Class I sources are star-disk systems still surrounded by massive infalling envelopes. The envelope absorbs optical and near-infrared photons from the central protostar and re-emits this radiation at far-IR and mm wavelengths where the envelope is optically thin. Class II sources also have broad SEDs, but the SED falls from $`2\mu `$m to $`25100\mu `$m. These sources contain a pre-main-sequence star surrounded by a circumstellar disk, which absorbs and reradiates emission from the central star. Class III sources have roughly blackbody SEDs that peak at $``$ 1 $`\mu `$m; radiation from a surrounding disk or envelope contributes negligible emission at 10–100 $`\mu `$m (Lada 1991). André et al. (1993) added the class 0 sources to this sequence to include very young objects at the beginning of the accretion process. These objects have narrow SEDs which peak at $``$ 100–300 $`\mu `$m (e.g., Motte & André, 2001, and references therein). In the current picture, the central star of a class 0 source has accumulated $``$ 50% of its final mass; in class I sources, the central star has $``$ 50% of its final mass.

The Taurus-Auriga cloud is one of the best sites for testing models of isolated, low mass star formation (see Kenyon & Hartmann, 1995, and references therein). This collection of filamentary clouds contains good samples of class I, II, and III pre-main sequence stars, a class 0 source, and several examples of interesting transition objects. The clouds are close enough (d = 140 pc; Kenyon, Dobrzycka, & Hartmann, 1994a) to allow high signal-to-noise observations at all wavelengths with moderate to high spectral and spatial resolution. Although the clouds are poor examples of clustered star formation, pre-main sequence stars in Taurus-Auriga provide interesting constraints on the formation of small close groups of pre-main sequence stars (Gómez et al., 1993; Armitage & Clarke, 1997; Vanhala & Cameron, 1998; Scalo & Chappell, 1999).

Here, we consider a comparison between the millimeter (mm) and near-infrared (near-IR) properties of class I sources in Taurus-Auriga. Motte & André (2001) demonstrate that moderate spatial resolution mm continuum observations divide these objects into unresolved, point-like sources with negligible infalling envelopes, extended sources with significant infalling envelopes, and several transition objects with modest infalling envelopes. In this division, the unresolved Taurus-Auriga sources might be edge-on disks similar to those discovered in HH30 IRS and HK Tau (Burrows et al., 1996; Stapelfeldt et al., 1998). Because the derived lifetimes in each evolutionary phase depend on the relative numbers of objects in each phase (e.g. Kenyon et al., 1990; Hogerheijde, 2001), the mm continuum data also provide better limits on the class I lifetime in Taurus-Auriga. Our goal is to test whether near-IR images yield a similar division of sources into extended and point-like objects. This study extends previous imaging surveys, where near-IR imaging data constrains radiative transfer models for the class I sample (Kenyon et al., 1993b; Whitney et al., 1997).

We describe the observations in §2, comment on source variability in §3, and consider the image morphology in §4. The paper concludes with a brief summary.

## 2 OBSERVATIONS

We acquired JHK images of candidate protostars in Taurus-Auriga with Stelircam at the Fred L. Whipple Observatory 1.2-m telescope on 5–12 October 1998 and 18–20 November 1999. Stelircam is a dual-beam IR imager equipped with two 256 $`\times `$ 256 InSb SBRC arrays. The field-of-view for these images is $``$ 300′′ $`\times `$ 300′′. We took nine 6 $`\times `$ 10 sec exposures of each field at HK and later at JK and offset the telescope by 45′′–60′′ after each integration. Poor weather prevented a complete set of J band observations in photometric conditions; the 1998 HK observations were photometric. We repeated observations of some targets to improve the signal-to-noise of faint nebulosity.

We produce calibrated object frames in each filter using standard techniques in the Interactive Data Language (IDL)111IDL is distributed by Research Systems, Inc. of Boulder, Colorado.. After correcting each image for non-linearity, we construct a nightly flat-field for each filter from the median of dark-subtracted image frames. The sky frames for each object image are a median-filtered set of flattened frames scaled to have the same background as the object image. We offset and coadd the nine frames in each filter to make a final deep field image which reveals fainter extended emission without saturating bright stars.

Table 1 lists HK photometry for each source. We use the NOAO IRAF222IRAF is distributed by the National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc. under contract to the National Science Foundation. routine PHOT to extract magnitudes in three apertures. We analyze the nine HK frames of each source – instead of the coadded frames – to obtain error estimates for the photometry directly from the data. Once we remove bad pixels, the typical standard deviation of nine measurements for one source is $`\pm `$0.01–0.02 mag at K and $`\pm `$0.02–0.03 mag at H. Our nightly photometric calibration is based on 15–25 observations of Elias et al. (1982) standard stars. The uncertainty in the calibration is $`\pm `$0.03 mag in each filter on each night. This uncertainty is comparable to the range in the zero-points and airmass corrections for the HK filters on each observing run. Adding the calibration error in quadrature with the standard deviations of the source data yields a typical photometric error of less than 0.05 mag in each filter.

To measure morphological information from the individual and coadded frames, we use Sextractor (Bertin & Arnouts, 1996). Sextractor is a modern software package which automatically detects, measures, and classifies sources from digital images. For each image, we set up the package to measure the FWHM $`f`$, elongation $`e`$, central surface brightness $`\mu `$, and the magnitude in 7′′ and 14′′ apertures for each candidate protostar. Table 2 lists the results. The FWHM is the full width at half maximum of the best-fitting one dimensional gaussian to the intensity profile of the source; the elongation is the ratio of the major axis to the minor axis of the best-fitting two dimensional gaussian to the intensity profile. The fits for both parameters use intensity levels 3$`\sigma `$ or more above the local noise background. We use the central surface brightness and aperture magnitudes to verify that fits for each source reach the same limiting surface brightness and to check the PHOT photometry. The Sextractor photometry agrees with PHOT to $`\pm `$ 0.02 mag or better. The limiting surface brightness in the coadded frames is $``$ 19.5 mag arcsec-2.

Figure 1 gives examples of contour maps derived from coadded K band images for Taurus-Auriga protostars (see also Tamura et al., 1991; Kenyon et al., 1993b; Whitney et al., 1997; Lucas & Roche, 1997). The images have a variety of morphologies, ranging from apparent point sources (04108+2803 and 04295+2251) to monopolar (04287+1801 and 04288+1802) and bipolar (04325+2402 and 04361+2547) nebulae with central point sources. Two objects are faint patches with no obvious point source (04166+2706 and 04368+2557). Whitney et al. (1997) measure large near-IR polarizations for sources with extended nebulae, confirming models where dust grains within the nebulosity scatter radiation from the central star into our line-of-sight (see also Lucas & Roche, 1998a). Gómez et al. (1997) identify optical and near-IR jets on \[S II\] and H2 images of many Taurus-Auriga protostars (see also Lucas & Roche, 1998b); Kenyon et al. (1998) note \[S II\] emission on optical spectra of several other sources. Jet emission and bipolar outflows are characteristic of embedded protostars (e.g., Reipurth, 1991; Moriarty-Schieven et al., 1992; Bontemps et al., 1996; Tamura et al., 1996; Reipurth & Bally, 2001).

The IRAS source 04166+2706 deserves special mention. Kenyon et al. (1990) discovered a near-IR counterpart with raster scans on the single channel MMT photometer. From near-IR imaging data, Kenyon et al. (1993b) identified a point source with very red colors near the IRAS position. This object is $``$ 2 mag brighter than the Kenyon et al. (1990) source and lies near the reddening band in the J–K,H–K color-color diagram. Our images recover the Kenyon et al. (1990) source as a patch of low surface brightness nebulosity similar to the nebula of 04368+2557 (see Whitney et al., 1997). This object is 52′′ N and 44′′ W of the Kenyon et al. (1993b) point source, with J2000 coordinates of RA = 04h19m43.2s and Dec = +271339′′. The error in these coordinates is roughly $`\pm `$ 2′′. We suggest that this object is the IRAS counterpart.

The triplet of embedded sources associated with IRAS sources 04181+2654 and 04181+2655 also merit a short comment. Kenyon et al. (1993a, 1993b) associate the northern component with IRAS 04181+2655 and the southern pair with IRAS 04181+2654. In this convention, the southern embedded source is 04181+2654A; the middle component of the three is 04181+2654B. Figure 1 labels these three sources. For completeness, the J2000 coordinates are RA = 04h21m11.6s and Dec = +270110′′ for 04181+2654A, RA = 04h21m10.5s and Dec = +270138′′ for 04181+2654B, and RA = 04h21m08.1s and Dec = +270222′′ for 04181+2655.

## 3 Variability

Most pre-main sequence stars vary at optical and infrared wavelengths (e.g., Joy, 1945; Bouvier et al., 1993; Choi & Herbst, 1996; Stassun et al., 1999). Many systems have regular variations with amplitudes of 0.1–1.0 mag and periods of 2–10 days (Eaton, Herbst, & Hillenbrand, 1995; Herbst et al., 1994, 2000). Other objects vary irregularly with similar amplitudes on similar or longer timescales. The periodic variations yield information on stellar rotational periods and the sizes, temperatures, and lifetimes of dark and bright spots on the stellar photosphere (Bouvier et al., 1993; Stassun et al., 1999). These data also provide tests of models for the inner disks and magnetospheres of low mass pre-main sequence stars (e.g., Mahdavi & Kenyon, 1998; Wood & Whitney, 1998; Wood et al., 2000).

To test whether any of the candidate protostars are variable, we compare with previous measurements obtained with single channel photometers (Kenyon et al., 1990) and near-IR imagers (Kenyon et al., 1993b, 1994b; Whitney et al., 1997). For sources where we have only one or two of our own measurements, we supplement the photometry with data from the compilation of Kenyon & Hartmann (1995). We derive the average brightness and the standard deviation for the candidate protostar and for another stellar source in the same field using 8′′ and 13′′–14′′ apertures. Although there are too few measurements for each source to construct a robust average, the standard deviation provides a reasonable first estimate of source variability.

The observations in this study use different instruments; we use the comparison stars to verify that the photometry from each instrument is on the same magnitude scale. We could not perform this test with single channel photometry; previous comparisons show that the differences between photometric systems are small compared to the source variations in our sample (Kenyon & Hartmann, 1995). As a final test, we check whether the mean variation of a protostar from one observing run to the next is zero. Using data from Kenyon et al. (1993b), Whitney et al. (1997), and this paper, the average change in brightness of the complete sample is 0.11$`\pm `$0.37 from 1991 to 1993–94 and 0.05$`\pm `$0.31 from 1993–94 to 1998–99.

Our analysis indicates that candidate protostars in Taurus-Auriga are variable. Histograms of standard deviations $`\sigma `$ in K and in H–K (Figure 2) demonstrate that candidate protostars have larger variations than comparison stars in the same field. The variation of a typical source has an amplitude of 0.2–0.3 mag at K. These fluctuations are similar to variations observed in optically visible T Tauri stars (e.g., Kenyon & Hartmann, 1995; Choi & Herbst, 1996; Stassun et al., 1999). Many sources have variable colors. Roughly half of the sample varied by 0.1–0.3 mag in H–K. Although large color variations are observed in some T Tauri stars, the variations are less than 0.1–0.2 mag in more than 70% of pre-main sequence stars in Taurus-Auriga (Kenyon & Hartmann, 1995, see also Carpenter, Hillenbrand, & Skrutskie 2001).

Without repeat observations on short timescales, we are unable to analyze the variability of Taurus-Auriga protostars in detail. The large brightness and color changes are consistent with the behavior expected from bright spots, dark spots, and obscuration models (see also Carpenter et al., 2001, and references therein). There is no correlation between brightness and color changes; some sources become bluer as they brighten while others become redder. Most models predict strong correlations between color and brightness evolution, but our observations are too few in number to discriminate between possible models.

The variations of several objects deserve special mention.

04018+2803B – This source is the more embedded of a pair of IRAS sources separated by 17$`^{\prime \prime }.`$5 on the sky. The large variation is based on a single observation in October 1998.

04181+2655 – The brightest of a triplet of IRAS sources, this object has varied consistently from one observing run to the next. The color variation of this source is among the largest in our sample.

04239+2436 – This source has a large nebulosity at J which is nearly invisible at K. It has one of the largest amplitudes in our sample, with a variation of more than 0.7 mag at K during the past decade.

Six other sources – 04016+2610, 04181+2654A, 04287+1801 (L1551 IRS5), 04361+2547, 04365+2535, and 04381+2540 – vary by 20% or more in brightness and in color. All other sources change by less than 0.1–0.2 mag (Table 1; columns 10–11). The historical variation of L1551 IRS5 is difficult to quantify due to large differences in observing techniques, but published data are consistent with the variation in our data.

## 4 Morphology

Image morphology provides a good constraint on models for embedded protostars. Whitney et al. (1997) demonstrate that the size of synthetic near-IR images correlates with the centrifugal radius $`R_c`$ and the infall rate $`\dot{M}`$ in the Terebey et al. (1984) infalling envelope model. The elongation correlates with the inclination of the source and the size of the outflow cavity (see also Whitney & Hartmann, 1993). Both parameters vary with the grain albedo. Grains with large albedo produce larger and more elongated near-IR images of protostars. Whitney et al. (1997) emphasize that the size and elongation provide better constraints on models when combined with polarization measurements.

The image morphology also provides a good comparison with mm continuum observations. Motte & André (2001) describe moderate spatial resolution mm-wave observations of the Taurus-Auriga sample. They divide the sources into two groups based on the ratio of the flux $`F`$ in 11′′ and 60′′ beams. Extended sources have $`F_{11}/F_{60}`$ 0.5–0.7; point sources have $`F_{11}/F_{60}`$ 1. Roughly half of the Taurus-Auriga sample is extended by this criterion.

To compare the near-IR morphology with mm continuum morphology, we use two measures to divide objects into extended sources and point-like sources. We follow Motte & André (2001) and define extended sources as objects where the flux in a 26′′ aperture exceeds the fluxes in 13′′ and 8′′ apertures (Table 1). Larger apertures yield poor quality photometry on our images. In a sample of $``$ 2000 field stars with moderate signal-to-noise photometry, the mean aperture corrections for our data are $`m_8`$$`m_{13}`$ = 0.09$`\pm `$0.02 and $`m_{13}`$$`m_{26}`$ = 0.11$`\pm `$0.02. The aperture corrections are independent of the filter. With this result, extended sources have $`m_8`$$`m_{13}>`$ 0.15 and $`m_{13}`$$`m_{26}>`$ 0.17. The size (FWHM) $`f`$ and elongation $`e`$ also provide good measures of source extension. In a sample of $``$ 3000 moderate signal-to-noise detections from our K-band images, the median size is $`f_{med}`$ = 1$`^{\prime \prime }.`$75$`\pm `$0$`^{\prime \prime }.`$30 where the error bar refers to the inter-quartile range (Figure 3). The median elongation in the K-band images is $`e_{med}`$ = 1.10$`\pm `$0.04. Point sources are well-defined in our images. We thus define extended sources to have $`f_K>`$ 2$`^{\prime \prime }.`$5 and $`e_K>`$ 1.2.

Our analysis indicates a good correlation between various measures for extended emission in the Taurus-Auriga class I sample (Figure 4). For sources with $`f_K`$ 5′′, there is a one-to-one relation between $`f_K`$ and $`K_8`$$`K_{13}`$ (Figure 4; lower left panel). Although this relation is weaker at larger sizes, the Spearman rank order test (Press et al., 1988) still indicates a significant correlation with a probability, $`p_s`$ 3 $`\times 10^4`$, of being drawn from a random sample. The elongation and size derived from Sextractor show a stronger correlation (Figure 4, lower right panel), with $`p_s`$ 5 $`\times 10^5`$.

The good correlation between size and elongation confirms predictions of radiative transfer models for protostars with infalling envelopes. In the Terebey et al. (1984) infall models with a bipolar outflow cavity, the size and the elongation of near-IR images are related to infall rate and the geometry of the outflow cavity. Protostars with denser envelopes, larger infall rates, and curved or dust-filled outflow cavities are larger and more elongated than protostars with less dense envelopes, smaller infall rates, or streamline cavities (Kenyon et al., 1993b; Whitney et al., 1997). Our data confirm the general correlation predicted by the models. Further comparisons between the data and the models will yield better tests of the models.

The K-band size and elongation are well-correlated with the mm continuum measurements. Motte & André (2001) define the concentration $`c`$ as the ratio of the fluxes in 11′′ and 60′′ apertures expressed as a percentage. Highly concentrated sources have $`c`$ = 100; very extended sources have $`c<`$ 30. Our near-IR measures $`f_K`$ and $`e_K`$ are inversely correlated with $`c`$ (Figure 4, middle panels). The relation between $`c`$ and $`e_K`$ is significant at the 3$`\sigma `$ level, with $`p_s`$ 4 $`\times 10^3`$; the size is more weakly correlated with $`c`$ ($`p_s`$ $`10^2`$). Motte & André (2001) use the concentration and the absolute fluxes to measure the envelope mass, $`M_{env}`$. Our K-band measures are more strongly correlated with $`M_{env}`$ than with $`c`$ (Figure 4, upper panels); $`p_s`$ 2 $`\times 10^3`$ is the Spearman rank probability for the correlation between $`f_K`$ and $`M_{env}`$ and $`p_s`$ $`10^3`$ is the probability for the correlation between $`e_K`$ and $`M_{env}`$.

The near-IR measurements at J and H also correlate with the mm continuum data. Because the short wavelength images have lower signal-to-noise, the correlations at J and H are 10% to 20% weaker than at K. These results exclude the deeply embedded source 04368+2557 (L1527 IRS) from the analysis. Including this source with a nominal size of $`f_K`$ 25′′ and an elongation of $`e_K`$ 2 yields smaller Spearman rank probabilities by a factor of 2–3.

These correlations provide additional confirmation of infall models for protostars. In the standard infall picture (Terebey et al., 1984; Adams et al., 1987), an optically thick envelope absorbs radiation from the central object and re-radiates this energy at mm wavelengths. Protostars with massive infalling envelopes are therefore more extended at mm continuum wavelengths than class II or class III sources with little or no infalling envelope. Once the central object generates a bipolar outflow, optical and near-IR radiation from the central object can escape through the optically thin cavity (e.g. Whitney & Hartmann, 1993; Kenyon et al., 1993b; Whitney et al., 1997). Some of this radiation scatters off material in the cavity and in the outer envelope, producing an extended, elongated, and highly polarized optical or near-IR image (Whitney et al., 1997). As material from the envelope continues to fall onto the disk, the envelope opacity declines and more optical and IR radiation from the central star can reach the observer directly. This process reduces the amount of mm continuum emission and scattered optical and near-IR emission. These objects thus appear more concentrated and more symmetrical at all wavelengths than more highly embedded objects. Our data confirm this general trend: class I sources that are extended at mm continuum wavelengths are extended and more elongated on near-IR images than their point-like counterparts.

These correlations may also offer a test of infalling envelope models. To the best of our knowledge, there are few quantitative predictions of the relationship between the optical or near-IR properties and the mm continuum properties of protostars. Published radiative transfer models have successfully modeled independently the spectral energy distributions, density distributions, and near-IR images of individual protostars. However, we are not aware of models that predict the observed correlations between the concentration or envelope mass derived from mm data and the size or elongation derived from optical or near-IR data. Our results suggest that observations can now test these predictions.

The near-IR data divide Taurus-Auriga class I sources into nearly the same two groups as the mm continuum data. We base this conclusion on Table 3, which compares our results with those of Motte & André (2001). Nine sources are extended at mm and near-IR wavelengths; eight sources are not extended at these wavelengths. Two sources – 04016+2610 and 04302+2247 – have large near-IR nebulae with little evidence for extended millimeter continuum emission from a massive envelope33304016+2610 is resolved at 2 cm and 3.5 cm (Lucas, Blundell, & Roche, 2000); this emission probably traces radiation from an ionized jet or inner disk as in L1551 IRS5 (Rodriguez et al., 1986). (Hogerheijde, 2001; Motte & André, 2001). Three sources – 04181+2654A, 04181+2654B, 04365+2535, and 04381+2540 – are extended at mm continuum wavelengths but are not extended at near-IR wavelengths based on our measurements of size and elongation. Three of these – 04181+2654A, 04365+2535, and 04381+2540 – show some evidence for an extended nebula at K from $`K_8`$$`K_{13}`$ and a visual inspection of the coadded images (see Figure 1). Higher signal-to-noise K-band images would provide definitive tests for extended emission in these sources.

Our data confirm the Motte & André (2001) conclusion that most ($``$ 60%) of the Taurus-Auriga class I sources are embedded protostars with massive infalling envelopes (see also Hogerheijde et al., 1998; Hogerheijde & Sandell, 2000; Chandler & Richer, 2000, and references therein). Because these objects are definitely extended in the mm continuum and probably extended at K, we include 04181+2654A, 04181+2654B, and 04381+2540 in this group (see Table 2). The envelope masses inferred from the millimeter data range from $``$ 0.09 $`M_{}`$ for 04239+2436 to $``$ 0.9 $`M_{}`$ for L1551 IRS5. The density distributions for these sources, $`\rho (r)r^n`$ with $`n`$ = 1.5–2, agree with the predictions of infall models on scales of $``$ 10,000 AU (Terebey et al., 1984; Chandler, Barsony, & Moore, 1998; Hogerheijde et al., 1998; Hogerheijde & Sandell, 2000; Motte & André, 2001). The SEDs and the near-IR morphologies and polarizations of these sources also agree with model predictions (Kenyon et al. 1993a,b; Whitney et al., 1997; Motte & André, 2001). All of these sources have molecular outflows (e.g., Bontemps et al., 1996) and 66% have emission from optical or near-IR jets (e.g., Gómez et al., 1997).

At least two Taurus-Auriga class I sources without extended mm emission might be transition objects between class I protostars and class II T Tauri stars. In addition to their large near-IR nebulae, 04016+2610 and 04302+2247 have molecular outflows, optical or near-IR emission from a jet, and large near-IR polarization as in true protostars (see also Heyer et al., 1990; Tamura et al., 1991; Kenyon et al., 1993b; Gómez et al., 1997; Whitney et al., 1997). However, both sources fail to meet the “true protostar” criterion of Motte & André (2001), $`M_{env}/M_{}>`$ 0.1. To estimate these ratios in 04016+2610 and 04302+2247, we assume the central star in each source lies on the birthline of Palla & Stahler (1999). The observed bolometric luminosities then yield approximate stellar masses of $`M_{}`$ 0.6 $`M_{}`$ for 04016+2610 and $`M_{}<`$ 0.1 $`M_{}`$ for 04302+2247. Motte & André (2001) estimate envelope masses $`M_{env}<`$ 0.04 $`M_{}`$ for 04016+2610 and $`M_{env}<`$ 0.01 $`M_{}`$ for 04302+2247. Thus, both sources have $`M_{env}/M_{}<`$ 0.1 and are not true protostars.

Based on mm aperture synthesis observations, Hogerheijde (2001) previously proposed that 04016+2610 is a transition object between true protostars and class II sources. The mm data are consistent with a Keplerian disk with a radius of $``$ 2000 AU surrounding a 0.65 $`M_{}`$ central star (see also Saito et al., 2001). Our mass estimate favors this interpretation over the proposal of sub-Keplerian motion around a more massive star. If 04016+2610 is a transition object, several common near-IR properties suggest that 04302+2247 might also be a transition object444Padgett et al. (1999) detect an edge-on disk on HST near-IR images of 04302+2247, but our ground-based near-IR images of 04302+2247 show more near-IR nebulosity than class II sources with edge-on disks such as HH30 IRS and HK Tau.. Aperture synthesis observations at mm wavelengths would test this proposal.

Other Taurus-Auriga class I sources are peculiar in having class I SEDs but little or no extended near-IR or mm continuum emission. The properties of this group differ significantly from the true protostars. Most do not have molecular outflows, optical/near-IR jet emission, or large near-IR polarization. As a group, these objects thus appear older, with smaller extended envelopes and less accretion activity than the protostars (see also Reipurth & Bally, 2001). There are several exceptions to this statement. Four sources – 04158+2805, 04260+2642, 04264+2433, and 04489+3042 – have \[S II\] emission; 04264+2433 has large near-IR polarization; near-IR spectra of 04295+2251 and 04489+3042 more closely resemble spectra of class I sources than T Tauri stars Greene & Lada (1996). This diversity of properties suggests that the peculiar class I sources do not belong to a single class of objects.

Table 2 lists our best estimates for the nature of these peculiar class I sources. We consider 04016+2610, 04260+2612, 04264+2433, 04302+2247 as transition objects between class I and class II. Extended, highly polarized near-IR nebulae demonstrate that 04016+2610 and 04302+2247 are not class II sources; 04181+2655, 04260+2612 and 04264+2433 have similar spectral energy distributions but no obvious near-IR nebulosity. The class I sources 04108+2803B, 04295+2251, and 04489+3042 have larger optical to far-IR flux ratios than 04016+2610 and other ‘transition candidates.’ These sources could be transition objects viewed along the rotational axis of the disk. The strong \[S II\] emission in 04295+2251 recalls similar emission in the edge-on disk systems HK Tau and HH30 IRS. Finally, 04108+2803A and 04158+2805 are more similar to class II sources and are probably heavily reddened T Tauri stars.

## 5 SUMMARY

Our near-IR observations provide a first detection of near-IR variability in a sample of class I protostars. The variations are similar in amplitude to those observed in T Tauri stars (Kenyon & Hartmann, 1995; Carpenter et al., 2001). The poor time coverage of the imaging data for this sample does not provide a useful measure of the timescale for these variations. If the variations are similar to those in T Tauri stars, observations on successive nights over weeks or months can measure the timescales. As an example, Hillenbrand et al. (2001) derive an 8 day period for the variation in the Becklin-Neugebauer object, an embedded protostar in Orion (see also Wood et al., 2000).

The near-IR morphologies of Taurus-Auriga class I sources are consistent with the mm continuum measurements of Motte & André (2001). Most ($``$ 60%) of these objects are true protostars, with massive envelopes surrounding a young star with $``$ 50% of its final mass. Some ($``$ 25%) of the class I sources are probably transition objects between true protostars and optically visible T Tauri stars; one or two of these might have edge-on disks similar to those observed in HH30 IRS and HK Tau. Two other class I sources are probably heavily reddened T Tauri stars;.

Making accurate classifications for the “peculiar” class I sources requires additional data. Radio observations of some class I sources using higher spatial resolution might reveal compact nebulae or large disks as in 04016+2610 (Hogerheijde, 2001). Deeper near-IR images may also reveal extended nebulosity. We have marginal detections of faint nebulosity on coadded near-IR images of 04181+2655 and 04264+2433; higher signal-to-noise observations might provide a better indication of the extended nature of these and other point-like class I sources. Because molecular outflows and optical/near-IR jets are more common among class I sources and transition objects than among T Tauri stars with class II SEDs (Gómez et al., 1997; Kenyon et al., 1998; Reipurth & Bally, 2001), deeper optical and near-IR surveys for jet emission and more extensive outflow surveys can also help to classify these objects. These new data would provide a better understanding of the evolution from class I to class II among low mass protostars.

Finally, our results offer a new test of infall models for class I protostars. In the standard picture, a massive envelope produces mm and far-IR continuum emission; near-IR and optical radiation from the central object escape through the optically thin outflow cavity. As the infalling envelope becomes less massive and less opaque, the angular size of the outflow cavity should shrink. Because optical and near-IR radiation can escape the envelope more easily, protostars should become less elongated as the optical depth through the envelope decreases. Our near-IR morphological data confirm this expectation: $`f_K`$ and $`e_K`$ of class I sources correlate with $`M_{env}`$. The standard model also predicts a general anti-correlation between the near-IR measurements $`f_K`$ and $`e_K`$ and the mm measurement $`c`$ of Motte & André (2001). As the envelope mass declines, the unresolved central star and the disk produce more of the mm continuum flux; sources thus become more concentrated at mm wavelengths as they become less elongated and less extended at optical and near-IR wavelengths. To the best of our knowledge, however, theoretical radiative transfer models of protostars do not make specific predictions for the variation of $`f_K`$ or $`e_K`$ as a function of $`c`$ or $`M_{env}`$. Our results demonstrate that near-IR and mm continuum data can test these predictions.

We thank C. Lada for advice and comments on this project. M. Geller and an anonymous referee made many useful comments on the manuscript. W. Wyatt and J. Grimes provided technical assistance. S. P. thanks the National Science Foundation for their support of the REU program at the Smithsonian Astrophysical Observatory where most of this work was carried out. S. K. thanks Dr. S. Starrfield and Arizona State University for their hospitality during the completion of this project.


text:
# Classical Computation of Elliptic Flow at Large Transverse Momentum

## I Introduction

The azimuthal anisotropy of the particles produced in heavy ion collisions is quantified by the Elliptic flow parameter $`v2`$. It is defined to be the second Fourier moment of the azimuthal distribution,

$`v_2=\mathrm{cos}(2\varphi ){\displaystyle \frac{_\pi ^\pi 𝑑\varphi \mathrm{cos}(2\varphi )p_t𝑑p_t\frac{d^3N}{dyp_tdp_td\varphi }}{_\pi ^\pi 𝑑\varphi p_t𝑑p_t\frac{d^3N}{dyp_tdp_td\varphi }}}.`$ (1)
The RHIC (Relativistic Heavy Ion Collider) experiments have measured the centrality dependence of $`v_2`$ and find a value for $`v_2`$ Expt1  that is significantly larger for non-central collisions than previous measurements at lower energies CERN . Hydrodynamic models do a reasonable job of describing the centrality and mass dependence of $`v_2`$ Hydro1 ; Hydro2 . The RHIC experiments have also measured $`v_2(p_t)`$, which is defined to be

$`v_2(p_t)={\displaystyle \frac{_\pi ^\pi 𝑑\varphi \mathrm{cos}(2\varphi )\frac{d^3N}{dyp_tdp_td\varphi }}{_\pi ^\pi 𝑑\varphi \frac{d^3N}{dyp_tdp_td\varphi }}}.`$ (2)
The data Expt1 ; Expt2  show that, for a wide range of centralities, $`v_2(p_t)`$ rises up to $`1.5`$ GeV, flattens and then remains flat up to the highest transverse momentum measured $`p_t6`$ GeV. For peripheral collisions (specifically 34%–85% central collisions), $`v_2(p_t)`$ can be as large as $`25\%`$ at $`p_t3`$ GeV.

Hydrodynamic models describe the initial rise in $`v_2(p_t)`$ but overpredict the data at large $`p_t`$. The models predict a continued rise except perhaps in a Blast wave parametrization with extreme assumptions Blast . Partonic descriptions which account for the quenching of high $`p_t`$ “jets” as they traverse the hot QCD matter do not describe the data either GVW -unless unusually large partonic cross-sections are assumed Molnar .

In this note, we discuss the computation of $`v_2`$ in the Color Glass Condensate (CGC) approach to nuclear collisions MV ; KMW ; KV . At very high energies, the density of partons per unit area in the colliding nuclei becomes very large and saturates at a scale $`Q_s\mathrm{\Lambda }_{QCD}`$ GLR . For RHIC, estimates give $`Q_s12`$ GeV GyulassyMcLerran ; Mueller ; KV ; KN . The typical occupation number of partons is $`1/\alpha _s(Q_s)>1`$. Thus classical methods can be used to describe nuclear collisions at high energies. The CGC approach has been used to describe successfully various aspects of the bulk properties of charged hadrons in heavy ion collisions, including the multiplicity KV , the centrality dependence of the multiplicity KN , the energy and rapidity dependence KL  and the $`p_t`$ dependence of inclusive hadron spectra Juergen .

Elliptic flow, as defined in Eq. 1, is dominated by contributions from small momenta. In the classical approach, this involves computing diagrams to all orders in $`Q_s/p_t`$. This computation is very difficult to perform analytically. It can be computed numerically Yasushi . However, when the transverse momenta is large, $`p_tQ_s`$, the classical equations can be linearized and the contribution to $`v_2(p_t)`$ can then be computed analytically. This computation is performed in this paper. We find that the contribution to $`v_2(p_t)`$ from classical gluon fields at large momenta is small and differs significantly from the measured distribution at high $`p_t`$. The origin of the empirically observed behavior of $`v_2(p_t)`$, at large $`p_t`$, must arise from a source outside the classical perturbative expansion. Interestingly, even though the contribution from high momenta is small, this contribution is generated at very early times.

The paper is organized as follows. In the following section, we set up the problem of computing the classical gluon distribution at large transverse momentum. In section 3, we perform the computation. A brief final section summarizes our results. Some details of the computation are relegated to an appendix.

## II The Classical Gluon Distribution

Within the framework of the McLerran-Venugopalan model, Kovner, McLerran, and Weigert (KMW) KMW  (see also Refs. KR  and GyulassyMcLerran ) computed the classical spectrum of hard gluons (with $`p_t>Q_s`$) radiated from a collision of two infinitely large nuclei with uniform color charge distributions. In order to calculate the azimuthal distribution of the radiated gluons, we generalize this work to a finite nucleus with a non-uniform color charge distribution. Gradients in the charge distribution then give rise to space-momentum correlations of the radiated gluons. The typical momentum scale of the produced gluons is $`Q_s(x_t)`$. A typical spatial gradient is $`1/R`$. Therefore, the typical space-momentum correlation is of order $`(Q_sR)^1`$. The final expression for $`v_2`$ is actually of order $`(Q_sR)^2`$ since elliptic flow responds to the quadrupole moment of the color charge distribution.

We first review the perturbative solutions to the classical Yang Mills equations in order to establish the notation. In the McLerran-Venugopalan model, two nuclei (A and B) collide at very high energies. The valence partons generate classical fields at very small Bjorken $`x`$. The subsequent evolution of these classical fields is determined by the equations of motion. Therefore, the solution to the equations of motion is a functional of the valence color charge density (per unit area per unit rapidity) in nuclei A and B – $`\rho _A(x_t)`$ and $`\rho _B(x_t)`$, respectively. The expectation of any quantity $`O`$ is first expressed in terms of the classical solution and then averaged over all possible valence charge distributions with a weight function. In the McLerran-Venugopalan model, a Gaussian ansatz is taken for the weight function. Specifically, for a functional O($`\rho _A`$,$`\rho _B`$) we have the following average:

$$O_\rho \left[D\rho _A\right]\left[D\rho _B\right]O(\rho _A,\rho _B)\mathrm{exp}\left(d^2x_t\frac{\text{Tr}\left(\rho _A^2(x_t)\right)}{\mu _A^2(x_t)}d^2y_t\frac{\text{Tr}\left(\rho _B^2(y_t)\right)}{\mu _B^2(y_t)}\right).$$
(3)
At large transverse momentum, the density of charges is small and the field strength is weak. Therefore the classical Yang-Mills equations can be linearized and the solution can be expressed as a functional of the sources in closed form. The final distribution of gluons can then be related to the averages of the gluon field. The analysis of KMW  shows that the distribution of radiated gluons is given by

$`{\displaystyle \frac{dN}{dyd^2k_t}}`$ $`{\displaystyle d^2x_td^2y_te^{ik_tx_t}e^{ik_ty_t}\text{Tr}[^i\varphi _A,^k\varphi _B](x_t)[^j\varphi _A,^l\varphi _B]^{}(y_t)_\rho }`$ (4)
$`\times (\delta ^{ij}\delta ^{kl}+ϵ^{ij}ϵ^{kl}),`$
where $`\varphi _A`$ ($`\varphi _B`$) is the 2-d Coulomb potential associated with the color charge distribution $`\rho _A`$ ($`\rho _B`$),

$$_T^2\varphi _A(x_T)=g\rho _A.$$
(5)
In what follows, we evaluate the average over $`\rho `$ and and perform the necessary Fourier transforms in order to determine the gluon distribution to order $`(Q_sR)^2`$. With this gluon distribution, we determine the azimuthal anisotropy, $`v_2(p_t)`$.

## III Computation of $`v_2(k)`$ for $`k>Q_s`$

The averages over the color sources in Eq. 4 can be performed most simply in Fourier space. Below, we adopt the convention that repeated spatial indices are integrated over and repeated Fourier indices are integrated over with the appropriate two dimensional measure, $`\frac{d^2k}{(2\pi )^2}`$. In Fourier space, the average of $`O(\rho _A,\rho _B)`$ is then

$$O_\rho =\left[D\stackrel{~}{\rho _A}\right]\left[D\stackrel{~}{\rho _B}\right]O(\rho _A,\rho _B)\mathrm{exp}\left(\text{Tr}(\stackrel{~}{\rho }_A)_k^{}(C_A)_{kk^{}}(\stackrel{~}{\rho }_A)_k^{}\text{Tr}(\stackrel{~}{\rho }_B)_k^{}(C_B)_{kk^{}}(\stackrel{~}{\rho }_B)_k^{}\right)$$
(6)
where $`(\stackrel{~}{\rho }_A)_k=e^{ikx}\rho _A(x)`$ is the Fourier transform of $`\rho _A(x)`$ and $`(C_A)_{kk^{}}`$ is given by

$$(C_A)_{kk^{}}=e^{ikx}\frac{\delta _{xx^{}}}{\mu _A^2(x)}e^{ik^{}x^{}}.$$
(7)
Then pairwise correlations for nucleus A are given by

$$(\stackrel{~}{\rho _A})_k^{}(\stackrel{~}{\rho _A})_l=(C_A^1)_{kl}=e^{i(kl)x}\mu _A^2(x)\stackrel{~}{\mu ^2}_A(kl).$$
(8)
Fourier transforming each of the $`\varphi `$’s, the gluon distribution of Eq. 4 becomes

$$\frac{dN}{dyd^2k}e^{ikx_t}e^{il_1x_t}e^{il_2x_t}e^{iky_t}e^{il_3y_t}e^{il_4y_t}\frac{l_1}{l_1^2}(\stackrel{~}{\rho }_A)_{l_1}\frac{l_2}{l_2^2}(\stackrel{~}{\rho }_B)_{l_2}\frac{l_3}{l_3^2}(\stackrel{~}{\rho }_A)_{l_3}^{}\frac{l_4}{l_4^2}(\stackrel{~}{\rho }_B)_{l_4}^{}_\rho .$$
(9)
Contracting the $`\rho _A`$ variables and the $`\rho _B`$ variables we find

$$\frac{dN}{dyd^2k}\frac{\delta ^{ij}\delta ^{kl}+ϵ^{ij}ϵ^{kl}}{l_1^2l_3^2(l_1+k)^2(l_3+k)^2}l_1^il_3^k(l_1+k)^j(l_3+k)^l(C_A^1)_{l_1l_3}^{}(C_B^1)_{(l_1+k)(l_3+k)}.$$
(10)
Specializing now to the case where the target and projectile are identical, restoring the implicit integration, and using $`(C_A^1)_{kl}=(C_B^1)_{kl}=\stackrel{~}{\mu ^2}(kl)`$ we obtain

$$\frac{dN}{dyd^2k}\frac{d^2l_1}{(2\pi )^2}\frac{d^2l_3}{(2\pi )^2}\frac{\delta ^{ij}\delta ^{kl}+ϵ^{ij}ϵ^{kl}}{l_1^2l_3^2(l_1+k)^2(l_3+k)^2}l_1^il_3^k(l_1+k)^j(l_3+k)^l\left|\stackrel{~}{\mu ^2}(l_1l_3)\right|^2.$$
(11)
For the purposes of illustration, we consider a Gaussian model charge distribution comment

$`\mu ^2(\stackrel{}{x})`$ $`=`$ $`{\displaystyle \frac{Q^2}{2\pi R_xR_y}}e^{\frac{x^2}{2R_x^2}\frac{y^2}{2R_y^2}}`$ (12)
$`\stackrel{~}{\mu ^2}(\stackrel{}{k})`$ $`=`$ $`Q^2e^{\frac{1}{2}k_x^2R_x^2\frac{1}{2}k_y^2R_y^2},`$ (13)
where $`R_x`$ and $`R_y`$ denote the transverse spatial extent (in the $`x`$ and $`y`$ direction) of the overlap region of the two nuclei and $`Q^2`$ is the net color charge squared in the overlap region.

For this model distribution, we see that $`\left|\stackrel{~}{\mu ^2}(l_1l_3)\right|^2`$ is a sharply falling function. Thus $`\left|\stackrel{~}{\mu ^2}(l_1l_3)\right|^2`$ is approximately a delta function. Indeed, when $`\left|\stackrel{~}{\mu ^2}(l_1l_3)\right|^2`$ is a delta function Eq. 11 reduces to the result of KMW  for the gluon distribution. We therefore approximate the integral over $`\mu ^4`$ as

$`{\displaystyle \frac{d^2z}{(2\pi )^2}\left|\stackrel{~}{\mu ^2}(z)\right|^2f(z)}`$ $``$ $`(\stackrel{~}{\mu ^4})^{(0)}f(0)+{\displaystyle \frac{1}{2}}(\stackrel{~}{\mu ^4})_{ij}^{(2)}\left({\displaystyle \frac{^2f}{z_iz_j}}{\displaystyle \frac{\delta ^{ij}}{2}}{\displaystyle \frac{^2f}{_{z_l}_{z_l}}}\right)_{z=0}`$ (14)
$`+{\displaystyle \frac{1}{2}}(\stackrel{~}{\mu ^4})_{ij}^{(2)}\left({\displaystyle \frac{\delta ^{ij}}{2}}{\displaystyle \frac{^2f}{_{z_l}_{z_l}}}\right)_{z=0}.`$
In writing the above equation we have assumed that the charged distribution $`\left|\stackrel{~}{\mu ^2}(z)\right|^2`$ is symmetric in $`z`$ and defined the moments of this distribution as

$`(\stackrel{~}{\mu ^4})^{(0)}`$ $``$ $`{\displaystyle \frac{d^2z}{(2\pi )^2}\left|\stackrel{~}{\mu ^2}(z)\right|^2}`$ (15)
$`(\stackrel{~}{\mu ^4})_{ij}^{(2)}`$ $``$ $`{\displaystyle \frac{d^2z}{(2\pi )^2}z_iz_j\left|\stackrel{~}{\mu ^2}(z)\right|^2}.`$ (16)
For the problem at hand, we have

$$\frac{dN}{dyd^2k}\frac{d^2z}{(2\pi )^2}\left|\stackrel{~}{\mu ^2}(z)\right|^2f(k,z)$$
(17)
where $`f(k,z)`$ is given by

$$f(k,z)\frac{d^2p}{(2\pi )^2}\frac{\delta ^{ij}\delta ^{kl}+ϵ^{ij}ϵ^{kl}}{p^2(pz)^2(p+k)^2(pz+k)^2}p^i(pz)^k(p+k)^j(pz+k)^l.$$
(18)
We may write down the general form of the second derivatives of $`f(k,z)`$

$`\left({\displaystyle \frac{^2f(k,z)}{z_iz_j}}{\displaystyle \frac{\delta ^{ij}}{2}}{\displaystyle \frac{^2f(k,z)}{_{z_l}_{z_l}}}\right)_{z=0}`$ $``$ $`f^{(2)}(k^2)({\displaystyle \frac{k_ik_j}{k^2}}{\displaystyle \frac{\delta ^{ij}}{2}})`$ (19)
$`\left({\displaystyle \frac{^2f(k,z)}{_{z_l}_{z_l}}}\right)_{z=0}`$ $``$ $`g^{(2)}(k^2).`$ (20)
Now the elliptic flow parameter $`v_2`$ is given by

$$v_2(y,k)\frac{\frac{d\varphi }{2\pi }\mathrm{cos}(2\varphi )\frac{dN}{dykdkd\varphi }}{\frac{d\varphi }{2\pi }\frac{dN}{dykdkd\varphi }}.$$
(21)
Substituting Eq. 14 and Eq. 19 into this definition of $`v_2`$, we obtain to leading order in $`(Q_sR)^1`$

$$v_2(y,k)=\frac{(\mu ^4)_{xx}^{(2)}(\mu ^4)_{yy}^{(2)}}{8(\mu ^4)^{(0)}}\frac{f^{(2)}(k^2)}{f(k,z=0)}.$$
(22)
$`f(k,z=0)`$ and $`f^{(2)}(k^2)`$ are calculated in appendix A. The results are logarithmically divergent. As discussed in KMW , this is a consequence of the weak field expansion; higher order non-linear corrections cut off the behavior in the infra-red at a scale $`\alpha _s\mu `$. To leading log accuracy we find

$`f(k,z=0)`$ $``$ $`{\displaystyle \frac{1}{(2\pi )}}{\displaystyle \frac{1}{k^2}}\mathrm{log}\left({\displaystyle \frac{k^2}{(\alpha _s\mu )^2}}\right)`$ (23)
$`f^{(2)}(k^2)`$ $``$ $`{\displaystyle \frac{1}{(2\pi )}}{\displaystyle \frac{4}{k^4}}\mathrm{log}\left({\displaystyle \frac{k^2}{(\alpha _s\mu )^2}}\right).`$ (24)
With these expressions we obtain our final result for $`v_2`$

$$v_2(y,k)\frac{1}{2k^2}\frac{(\mu ^4)_{xx}^{(2)}(\mu ^4)_{yy}^{(2)}}{(\mu ^4)^{(0)}}.$$
(25)
It is remarkable (see the appendix) that there are no power or even logarithmic divergences from the ratios of integrals in our expression for $`v_2`$. For the Gaussian distribution of Eq. 12 the moments of the charge distribution are easily calculated. The resulting elliptic flow is

$$v_2(y,k)\frac{1}{4k^2}(\frac{1}{R_x^2}\frac{1}{R_y^2}).$$
(26)
The contribution of momenta $`k>Q_s`$ to the integrated $`v_2`$ (Eq. 1) is

$$v_2=\frac{_{Q_s}^{\mathrm{}}k𝑑kv_2(k)\frac{dN}{d^2k}}{_{Q_s}^{\mathrm{}}d^2k\frac{dN}{d^2k}}\frac{1}{Q_s^2}\left(\frac{1}{R_x^2}\frac{1}{R_y^2}\right).$$
(27)
For realistic values of $`Q_s1`$ GeV for RHIC and $`R_x3`$ fm for a peripheral collision, this contribution is a fraction of a percent. The significant contribution of classical fields to $`v_2`$ comes from momenta $`k<Q_s`$. Interestingly, our results suggest that $`v_2`$ for $`k>Q_s`$ is generated very early. The spatial and temporal components decouple already at proper time $`\tau =0`$. The temporal contributions are Bessel functions $`J_0(k\tau )`$ and $`J_1(k\tau )`$ which can be linearized into plane waves for times $`\tau <1/Q_s`$ when $`k>>Q_s`$.

## IV Summary

In this paper, we have computed analytically the contribution to the azimuthal anisotropy ($`v_2(k)`$) from classical fields at large transverse momenta, $`k^2>Q_s^2`$. We find that the azimuthal anisotropy has the simple form shown in Eq. 26 for a Gaussian charge distribution. Clearly, this behavior disagrees with the data and another mechanism must be found to explain the large momentum behavior. Non-flow azimuthal correlations may provide a natural explanation Yuri . An interesting feature of our result is that the temporal dependence of the perturbative classical fields decouples already at proper time $`\tau =0`$ KMW . Thus as Eq. 4 suggests, the momentum anisotropy is generated at very early times in the collision. This analysis may thus explain why (in numerical simulations of the Yang-Mills equations Yasushi ) a significant anisotropy is seen at early times after the collision.

## V Acknowledgements

We would like to thank D. Kharzeev, Y. Kovchegov, A. Krasnitz, L. McLerran, Y. Nara and K. Tuchin for useful discussions. R.V and D.T.’s work was supported by DOE Contract No. DE-AC02-98CH10886. R.V would also like to acknowledge RIKEN-BNL for support.

## Appendix A Evaluation of Integrals

To complete the evaluation of $`v_2`$ we expand Eq. 17 as a function of $`z`$ and extract the coefficients of the symmetric traceless tensor structures. We find after some algebra

$`f(k,z=0)`$ $`=`$ $`{\displaystyle \frac{d^2p}{(2\pi )^2}\frac{1}{p^2(k+p)^2}}`$ (28)
$`f^{(2)}(k^2)({\displaystyle \frac{k_ik_j}{k^2}}{\displaystyle \frac{\delta _{ij}}{2}})`$ $`=`$ $`(I^1)_{ij}+(I^2)_{ij}+(I^3)_{ij}`$ (29)
where the integrals are given by

$`(I^1)_{ij}`$ $`=`$ $`{\displaystyle \frac{d^2p}{(2\pi )^2}\frac{2}{p^2(k+p)^6}\left(k_ik_jk^2\frac{\delta _{ij}}{2}\right)}`$ (30)
$`(I^2)_{ij}`$ $`=`$ $`{\displaystyle \frac{d^2p}{(2\pi )^2}\frac{4}{p^2(k+p)^6}\left(\frac{k_ip_j+p_ik_j}{2}pk\frac{\delta _{ij}}{2}\right)}`$ (31)
$`(I^3)_{ij}`$ $`=`$ $`{\displaystyle \frac{d^2p}{(2\pi )^2}\left(\frac{2}{p^2(k+p)^6}+\frac{2}{p^6(p+k)^2}\right)\left(p_ip_jp^2\frac{\delta _{ij}}{2}\right)}.`$ (32)
$`f(k,z=0)`$ is logarithmically divergent at $`p=0`$ and $`\stackrel{}{p}=\stackrel{}{k}`$. It is identical to the expression obtained previously in KMW ; GyulassyMcLerran . As discussed there and shown explicitly in Ref. KV  this infra-red divergence is an artifact of the weak field expansion. To leading log accuracy we find

$$f(k,z=0)\frac{1}{(2\pi )}\frac{1}{k^2}\mathrm{log}\left(\frac{k^2}{(\alpha _s\mu )^2}\right).$$
(33)
Consider now the sum of the integrals in Eqs. 32. While each of these individually have power divergences, taken together, the result is again only logarithmically divergent. Note, as previously, the integrals are divergent at $`p=0`$ and $`\stackrel{}{p}=\stackrel{}{k}`$ and one obtains logarithmic contributions from both regions to the final result. For illustration, consider the contribution to the sum from the divergence around $`p=0`$. $`(I^1)_{ij}`$ is logarithmically divergent, $`(I^2)_{ij}`$ is finite, and $`(I^3)_{ij}`$ appears power divergent. However, $`(I^3)_{ij}`$ is actually only logarithmically divergent upon integration over the angle between $`p`$ and $`k`$. Contracting both sides of Eq. 29 with $`2(\frac{k_ik_j}{k^2}\frac{\delta _{ij}}{2})`$ and keeping only logarithmically divergent terms we obtain the contribution of the $`p=0`$ divergence to $`(I^1)_{ij}`$ and $`(I^3)_{ij}`$ respectively to be

$`2(I_1)_{ij}\left({\displaystyle \frac{k^ik^j}{k^2}}{\displaystyle \frac{\delta ^{ij}}{2}}\right)_{|\mathrm{p}|\alpha _\mathrm{s}\mu }`$ $`=`$ $`k^2{\displaystyle \frac{d^2p}{(2\pi )^2}\frac{2}{p^2(k+p)^6}}{\displaystyle \frac{1}{(2\pi )}}{\displaystyle \frac{1}{k^4}}\mathrm{log}\left({\displaystyle \frac{k^2}{(\alpha _s\mu )^2}}\right)`$
$`2(I_3)_{ij}\left({\displaystyle \frac{k^ik^j}{k^2}}{\displaystyle \frac{\delta ^{ij}}{2}}\right)_{|\mathrm{p}|\alpha _\mathrm{s}\mu }`$ $`=`$ $`{\displaystyle \frac{d^2p}{(2\pi )^2}\frac{2p^2\mathrm{cos}(2\varphi _{pk})}{p^6(k+p)^2}}{\displaystyle \frac{1}{(2\pi )}}{\displaystyle \frac{1}{k^4}}\mathrm{log}\left({\displaystyle \frac{k^2}{(\alpha _s\mu )^2}}\right)`$ (34)
where $`\varphi _{pk}`$ is the angle between $`p`$ and $`k`$. A detailed analysis shows that the contribution to the sum $`(I^1)_{ij}+(I^2)_{ij}+(I^3)_{ij}`$ from the logarithmic divergence around $`\stackrel{}{p}=\stackrel{}{k}`$ is identical to the sum of the divergent pieces in Eqs. 34.

Substituting the results of the above analysis into the expression for $`f^{(2)}(k^2)`$ we find

$$f^{(2)}(k^2)=\frac{1}{(2\pi )}\frac{4}{k^4}\mathrm{log}\left(\frac{k^2}{(\alpha _s\mu )^2}\right).$$
(35)
With these expressions for $`f(k,z=0)`$ and $`f^{(2)}(k^2)`$, the final result for $`v_2`$ given in Eq. 25 follows.


text:
# Strongly correlated wave functions for artificial atoms and molecules

## References