text:
# The Rotation of Low-Mass Pre-Main-Sequence Stars

## 1. The Angular Momentum Problem

The specific angular momenta of the interstellar clouds from which stars form are significantly higher that that of the Sun, a fact that has been recognized for many decades. The numbers shown in Table 1 give scale to the “angular momentum problem”, with dense molecular cores having six orders of magnitude higher specific angular momentum than the Sun. In many, and perhaps all, cases the clouds partially solve this problem by depositing much of their angular momentum in the orbital motions of stellar and planetary companions. As shown by Vogel & Kuhi (1981) and Hartmann et al. (1986), most of the angular momentum problem has been solved by the pre-main-sequence (PMS) stage of evolution. Nonetheless, a factor 10-100 reduction in specific angular momentum is still required to derive the solar rotation at 4.5 Gyr, and indeed to produce the rotation of main-sequence solar-mass stars at 100 Myr.

| Table 1. Specific angular momenta (cm2s-1) | |
| --- | --- |
| Dense molecular cores | $`10^{2122}`$ |
| Wide binary stars | $`10^{1920}`$ |
| Pre-main-sequence stars | $`10^{1617}`$ |
| Sun | 1015 |

Thus many recent investigations into the rotation of low-mass stars, both theoretical and observational, have focussed on identifying the torques modifying stellar angular momentum from protostars to the Sun. Figure 1 shows a schematic of those mechanisms presently thought to be in play and the timescales on which they operate. During the PMS phase much recent work has been done within the context of the magnetic disk-locking mechanism developed by Shu et al. (1994). Here the ionization fraction in a protostellar disk is taken to be sufficient to couple the stellar magnetic field to the disk. At a certain radius in the disk, set by the stellar magnetic field strength and the mass accretion rate, the magnetic torque on the disk material balances the viscous torque. The equilibrium configuration is one where the star is in corotation with Keplerian rotation at that radius. Higher mass accretion rates drive that radius toward the star, leading to more rapid stellar rotation and shorter rotation periods. Higher magnetic field strengths move the radius away from the star, leading to slower stellar rotation and longer rotation periods. Importantly, in this picture the rotation period is independent of the stellar radius. Thus, if the magnetic field and mass accretion rate are constant, the star remains ”locked” at a rotation period independent of its PMS contraction.

In this paper we review recent observational studies of PMS stellar rotation, with particular emphasis on their implications for this disk-locking model of early angular momentum evolution.

## 2. Rotation Period Distributions of Pre-Main-Sequence Stars

The basic data for stellar angular momentum studies are rotation period distributions. PMS stars have large star spots which, as the stars rotate, modulate the flux observed at Earth by as much as several percent. Thus periodicities derived from light curves provide a very accurate ($``$ 1%) measure of stellar rotation periods. Indeed, stellar rotation period is the most accurate physical parameter that we know for PMS stars.

A great deal of observational work has focussed on the Orion Nebula Cluster (ONC), since the high surface density of PMS stars allows CCD imaging to yield large numbers of rotation periods (e.g., Choi & Herbst 1996, Stassun et al. 1999, Rebull 2001, Herbst et al. 2002). The total count is over 1000, with $``$ 40% of all ONC stars yielding rotation periods. (See summary of Herbst in Stassun & Terndrup 2003.) The most recent ONC rotation period distributions are shown in Figure 2, as derived by Herbst et al. (2002). Most of the rotation periods are found to be between 1 and 10 days. There has been some concern that the marked decrease in number of rotation periods longer than 10 days is a selection effect resulting from reduced dynamo activity and spot sizes in slow rotators, but vsin i studies tend to argue against this (Rhode, Herbst, & Mathieu 2001).

Bimodality in the ONC period distribution has been the subject of extensive discussion in the recent literature. A bimodal structure is present at a statistically significant level for stars with masses greater than 0.4 M$`_{}`$ in the central region of the ONC (using the mass calibration of Siess, Dufour, & Forestini 2000). However, in the regions flanking the cluster core, Rebull (2001) does not find a bimodal distribution among stars of similar numbers, mass and age (specifically, similar location in the HR diagram), a difference that remains to be resolved. The importance of resolution lies in the attribution of the long-period peak ($``$ 8 days) in the bimodal distribution to slowly rotating, disk-locked stars, while more rapidly rotating stars have been suggested to have spun up as a result of PMS contraction after the disk lock is lost (Bouvier et al. 1993, Choi & Herbst 1996).

Perhaps as significant as the bimodality itself is that the range of rotation periods extends over an order of magnitude. With the discovery in the early 1980’s of rapidly rotating low-mass stars in the Pleiades, a challenge to theorists has been to create a broad distribution in stellar rotation at 100 Myr (see the excellent review of past work by Barnes, Sofia & Pinsonneault 2001 as well as in these proceedings). In fact, the relative width of the period distribution in the ONC at 1 Myr already mimics that in the Pleiades, as shown in Figure 3. Recently Tinker, Pinsonneault, & Terndrup (2001) have modeled the evolution of the ONC distribution to 100 Myr, with some success in recreating the observed Pleiades rotation period distribution. There remains much work to be done, however, as these models use a parameterized approach to the braking mechanisms, and rather long disk-locking lifetimes are required (3-8 x 106 yr).

Importantly, the shape of the PMS rotation period distribution depends upon stellar mass (Figure 2; Herbst et al. 2000, 2002). The period distribution for stars with masses less than 0.4 M$`_{}`$ is substantially more concentrated toward shorter rotation periods than the period distribution for higher mass stars. Interestingly, this mass dependence is also seen in the Pleiades at 100 Myr (Terndrup et al. 2000); once again the basic structure of the older rotation period distribution is established by 1 Myr.

## 3. Pre-Main-Sequence Stellar Rotation and Protostellar Disks

Considered simply, the disk-locking model might be expected to predict a prevalence of protostellar disks associated with the most slowly rotating PMS stars (i.e., the stars that are disk-locked) and a lack of disks associated with more rapidly rotating PMS stars (i.e., the stars no longer locked which have spun up with their contraction). This is in fact what Edwards et al. (1993) found in their initial study of Taurus and ONC PMS stars, using 2 $`\mu `$m excess emission as their diagnostic for the presence of disks.

Thus it was surprising when Stassun et al. (1999) found no correlation between 2 $`\mu `$m excess and stellar rotation period among the ONC PMS stars (Figure 4). Stassun et al. also used H$`\alpha `$ and/or Ca triplet emission line strength (filled circles in Figure 4) as a diagnostic for active mass accretion, again finding no correlation with rotation period. The same conclusion was drawn by Rebull (2001) using ultraviolet excess as a diagnostic for mass accretion. More recently Herbst et al. (2002) have argued for a weak but significant correlation of 2 $`\mu `$m excess with rotation period for ONC stars both above and below 0.4 M$`_{}`$. Even so, in their sample many rapid rotators show accretion disk diagnostics, while many slow rotators do not. A simple association of protostellar disks with slow PMS stellar rotation is not present in the ONC.

Importantly, the predictions of disk-locking theory for emission excesses are not yet secure. For example, disk-locking theory predicts that gaseous disks will not exist within the corotation radius, since accretion of disk material proceeds along the magnetic field lines to the stellar surface. For a stellar magnetic field strength of 1 kG, a mass accretion rate of 10-7 M$`_{}`$/yr, and a stellar mass of 0.5 M$`_{}`$, the corotation radius lies at 14 R$`_{}`$. Using simple disk emission models, Stassun et al. (2001) showed that a disk with such a large central hole may not show a 2 $`\mu `$m excess. Indeed, in their models the peak 2 $`\mu `$m excess derives from a disk with a central hole of radius 6 R$`_{}`$, which corresponds to a rotation period of 2 days for a 0.5 M$`_{}`$ star. Thus, if disk locking is in force, the predicted correlation of 2 $`\mu `$m excess with rotation period would be in the opposite sense of that found by Edwards et al. and Herbst et al.

This simple model for disk structure and emission is meant only to demonstrate that the observational predictions of disk-locking theory for excess emission are not yet clearly established. Furthermore, if we are to continue to use 2 $`\mu `$m excess as our disk diagnostic, we must develop a better understanding of the dynamics and emission properties of the dust. For example, Shu (private communication) has pointed out that while the gas will follow the magnetic field lines, dust may diffuse into the central hole. Alternatively, the emission from dust traveling along stellar field lines is not yet known; presumably this will depend greatly on the surface area of the flux tubes.

Significantly, Stassun et al. (2001) used 10 $`\mu `$m excess to search for disks among PMS stars of all rotation periods showing no 2 $`\mu `$m excess. In only three cases among 32 stars was 10 $`\mu `$m disk emission found in the absence of near-infrared emission. Evidently, slow and rapid rotators without 2 $`\mu `$m excess in fact do not have protostellar disks to brake the stars.

More fundamentally, the disk-locking model does not in itself predict a preferred rotation period; the rotation period of a given star is set by the combination of stellar magnetic field, mass accretion rate, and stellar mass. Thus it is not evident that all disk-locked PMS stars should have similar rotation periods. Neither is it evident that all disk-locked PMS stars should rotate slowly. Indeed, given the certainty that at least the mass accretion rates decrease steadily from the protostar phase to the late-PMS phase, one might expect that, given disk locking, stellar rotation would evolve from rapid to slow (due to decreasing mass accretion rate) to rapid again (after the release of disk locking and subsequent stellar contraction). The findings of Greene & Lada (2002) that protostars are more rapidly rotating than PMS stars are compatible with this conjecture. Of course, the stellar magnetic fields themselves also may vary significantly over the same timescales, and both mass accretion rates and stellar magnetic fields will vary from star to star at any given time. So a strong correlation of stellar rotation with the simple presence of disks would not be expected. Perhaps the weak correlation found by Herbst et al. (2002) might be explained by the one-directional trend of any star to spin up once disk-locking ends.

Finally, it is worth noting that mass accretion rates have been observed to vary stochastically on short timescales, and if FU Ori events are part of every PMS star’s history the accretion rates also have large excursions on long timescales. Indeed such events of high mass accretion rates might act to spin up PMS stars. Critically, the response time of stellar rotation to variations in the disk-locking radius must be explored. Hartmann (2002) suggests that the response times are quite long ($``$ 106-7 yr), perhaps so much so that the effectiveness of any disk-locking during the PMS phase is limited.

## 4. Rapid Stellar Rotation

With a high frequency of observation (1 hr-1), Stassun et al. (1999) were able to identify a number of very rapid rotators with sub-day periods. The shortest PMS rotation period discovered to date is 0.27 days. As shown in Figure 5, some of these rapid rotators are spinning within a factor of two of their critical velocities. Furthermore, ONC stars with masses of 0.2 - 0.5 M$`_{}`$ will contract by factors of several before reaching the ZAMS. Finally many of the most rapidly rotating stars do not show any evidence of associated protostellar disks.

Thus these stars pose an interesting challenge. In the absence of angular momentum loss these stars would exceed their critical velocity as they approach the ZAMS. In addition, they have no evident means for rapid braking over that contraction time, although winds will have some effect. Rebull (2001) has pointed out that the calculation of a star’s critical velocity is dependent upon uncertain mass and radius determinations, and indeed the atmospheres and colors of these rapidly rotating stars may differ from standard models. With that caveat acknowledged, these stars give motivation for theorists to investigate the pre-main-sequence evolution of very rapidly rotating stars.

It may be physically significant that these rapidly rotating stars were discovered in the ONC and not in the earlier Taurus-Auriga studies, where no sub-day rotation periods have been found. Clarke & Bouvier (2000) compared the Taurus-Auriga and ONC period distributions and determined that they differ at a confidence level in excess of 99%, in the sense that the ONC has a larger population of rapid rotators. Their conjectures for the origin of this physical difference include a larger population of tidally locked binaries in the ONC, weaker stellar magnetic fields in the ONC, a higher frequency of stellar encounters in the ONC cluster destroying disks and ending disk-locking at an earlier age, or differing initial conditions. As they discuss, none of these explanations are without difficulties. (As an aside, Stassun and Mathieu have underway a spectroscopic monitoring study of the ONC rapid rotators. As yet a high frequency of spectroscopic binaries have not been found.)

If this environmental difference stands the test of further observation, it may prove seminal for our understanding of the origins of stellar angular momentum. Most critical will be extensive rotation period studies in other Taurus-Auriga- like regions so that the population differences can be put on more secure ground.

## 5.  Pre-Main-Sequence Evolution of Stellar Angular Momentum

As many modeling investigations have found, stellar winds alone do not seem capable of providing the angular momentum loss required by young stars before reaching the main sequence (e.g., Tinker et al. 2002). Thus most recent scenarios for PMS angular momentum evolution have the general form of strong braking due to disk-locking during the protostar phase and early in PMS evolution, followed by weak braking by winds after the loss of these disks. The external torques provided by winds are thought to be sufficiently weak as to permit substantial spin-up of the stars as they contract to the main sequence.

Such a scenario is, in principle, straightforwardly testable by observing the evolution of stellar rotation period distributions with age. This approach has been somewhat stymied by the paucity of rich clusters with ages between the ONC and the Pleiades (roughly 1 Myr to 100 Myr), but nonetheless observational work has begun. Perhaps befitting the early stage of the work, the first studies have not obtained findings consistent with either these scenarios or with each other.

Rebull et al. (2002 and these proceedings) have analyzed both the rotation period and vsin i distributions in the core of the ONC, in the fields flanking the ONC, and in NGC 2264 (an older cluster with an age of roughly 2-4 Myr). They take stellar radius as a diagnostic of age for stars of a given spectral type. Here stellar radii are derived from luminosities and effective temperatures , while spectral type is taken as a diagnostic for stellar mass for stars on the Hayashi portion of their PMS evolution (i.e., evolution at nearly constant effective temperature). They find no evidence up to an age of 3 Myr for stars spinning up as they contract down the Hayashi track.

This observational result is essentially that expected from strong disk locking, where the rotation period is held fixed by coupling to the disk. However, this interpretation is problematic in that most of the stars in the sample lack 2 $`\mu `$m evidence for disks. Conceivably 2 $`\mu `$m excess is not a sufficiently sensitive diagnostic for disks (e.g., Lada et al. 2000), but as noted above Stassun et al. (2001) found very few cases of stars not showing 2 $`\mu `$m excess emission yet having disk emission at 10 $`\mu `$m.

In contrast to this result, Lamm et al. (these proceedings) find that the rotation period distribution in NGC 2264 is significantly different from that in the ONC. Specifically, they find that the peaks in both the higher- and lower-mass rotation period distributions are at periods roughly a factor of 2 shorter than found in the ONC. (Although, importantly, the relative shapes of the higher- and lower-mass distributions are very similar to the ONC, with the lower mass stars NGC 2264 stars again rotating more rapidly.) If between 1 Myr and 2-4 Myr angular momentum is conserved, then the more rapid rotation in NGC 2264 could be attributed simply to the contraction of the stars to smaller radii. Of course, in this case the conclusion would be that disk locking is not significant after 1 Myr, precisely the opposite of the indication from the Rebull et al. result.

Clearly we are at an early stage of this line of investigation, particularly since both Lamm et al. and Rebull et al. base their analyses on (independently derived) NGC 2264 rotation periods. Equally clearly such observations have the potential to define stellar angular momentum evolution from the protostar to the ZAMS, and so must be pursued with vigor.

## 6. Future Directions

While the rich statistical studies have provided important foundations and insights, we now need direct attacks on the physics of angular momentum evolution during the protostar and PMS phases of stellar evolution. Within the context of testing the disk-locking model, future observations should include:

1) Direct observation of the kinematics and spatial distribution of gas in disks. New high-spectral-resolution infrared spectrographs are permitting observation of CO lines in disks. Not only are these gas tracers more sensitive to the presence of disks than dust emission, they provide direct spatial information on the emission sources under the presumption of Keplerian rotation. A recent study of CO fundamental lines by Najita, Carr, & Mathieu (2003) suggested that the inner edges of circumstellar gas disks lie within the corotation radii, but these are only the first observations of a rich direction of study.

2) Extensive measures of mass accretion rates and stellar magnetic fields. Mass accretion rates have been under intense study for over a decade. Recently Johns-Krull & Gafford (2002) have studied relationships predicted by disk-locking among stellar mass, radius, accretion rate and rotation period. Under the assumption that the magnetic field strength does not vary appreciably from star to star, they find support for the Shu et al. (1994) formulation of disk-locking, extended to include non-dipole field topologies. Magnetic field strengths of PMS stars have only recently been measurable (Johns-Krull et al. 2001), and will ultimately allow detailed tests of whether inner disk structures are being molded by magnetic fields.

On the theoretical side, we need:

1) Dynamical studies of stellar rotation evolution in the presence of stellar magnetic fields and accretion disks. Of particular interest will be timescales for rotation evolution given evolving and variable accretion rates. If magnetic coupling to accretion disks occurs, do young stars achieve equilibrium configurations? More generally, what are the evolutionary paths for stellar rotation from the protostar to the birthline to the main sequence? As briefly discussed in Section 3, the paths may be complex, with even the sign of the evolution changing with time. These studies will be a major computational challenge, but they are necessary in order to predict observables, such as stellar rotation distributions.

2) Fully consistent magnetohydrodynamic modeling across the boundary of the stellar surface. In recent studies either the braking mechanisms are parameterized and the stellar interior is carefully done (e.g, Barnes et al. 2001, Tinker et al. 2002) or the physics of braking is carefully done with the interior structures approximated (Collier Cameron & Campbell 1993). Recently a new paradigm for low-mass stellar rotation developed by Barnes (2003) rather naturally explains the evolution of stellar rotation distributions in main-sequence stars, including bimodalities and mass dependences observed in the period distributions. This model provides a framework for connecting internal and external magnetohydrodynamic processes. Ultimately these new ideas must interface PMS period distributions at 1 Myr, and the same magnetohydrodynamic issues addressed for PMS stars and younger.

Every indication is that a strong braking mechanism is needed if the rotational properties of protostars and young PMS stars are to evolve into the rotational properties of main-sequence stars. Magnetic disk-locking is a theoretically attractive mechanism for such strong braking. Observational indications of active disk locking at ages of 1 Myr and earlier are beginning to appear, most notably the detailed tests of Johns-Krull & Gafford (2002), the rapid rotation of protostars compared to PMS stars at 1 Myr, and perhaps the weak correlation of near-infrared excess with longer rotation period in the ONC. A critical outstanding question is the duration of disk locking, and in particular whether it can remain effective for several Myr into the pre-main-sequence phase of evolution.

Acknowledgements - It is a pleasure to acknowledge the important contributions of many colleagues past and present to the results presented here, with particular thanks for long-standing collaborations and associations with Bill Herbst, Luisa Rebull, and Keivan Stassun.

## References

Barnes, S. 2003, ApJ, in press

Barnes, S., Sofia, S. & Pinsonneault, M. 2001, ApJ, 548, 1071

Bouvier, J., Cabrit, S., Fenandez, M., Martin, E.L., & Matthews, J.M. 1993, A&A, 272, 176

Choi, P.I. & Herbst, W. 1996, AJ, 111, 283

Clarke, C.J. & Bouvier, J. 2000, MNRAS, 319, 457

Cameron Collier, A. & Campbell, C.G. 1993, A&A, 274, 309

Edwards, S. et al. 1993, AJ, 106, 372

Greene, T.P. & Lada, C.J. 2002, AJ, 124, 2185

Hartmann, L. 2002, ApJ, 566, L29

Hartmann, L., Hewett, R., Stahler, S., & Mathieu, R.D. 1986, ApJ, 309, 275

Herbst, W., Bailer-Jones, C.A.L., Mundt, R., Meisenheimer, K. & Wackermann, R. 2002, A&A, 396, 513

Herbst, W., Rhode, K.L., Hillenbrand, L.A. & Curan, G. 2000, AJ, 119, 261

Hillenbrand, L.A. 1997, AJ, 113, 1733

Johns-Krull, C. M., Valenti, J. A., Saar, S. H., & Hatzes, A. P. 2001, ASP Conf. Ser. 223: 11th Cambridge Workshop on Cool Stars, Stellar Systems and the Sun, 11, 521

Lada, C.J. et al. 2000, AJ, 120, 3162

Najita, J., Carr, J.S. & Mathieu, R.D. 2003, ApJ, in press

Rebull, L.M. 2001, AJ, 121, 1676

Rebull, L.M., Wolff, S.C., Strom, S.E. & Makidon, R.B. 2002, AJ, 124, 546

Rhode, K.L., Herbst, W. & Mathieu, R.D. 2001, AJ, 122, 3258

Seiss, L., Dufour, E. & Forestini, M. 2000, A&A, 358, 593

Shu, F., Najita, J., Ostriker, E., Wilkin, F., Ruden, S. & Lizano, S. 1994, ApJ, 429, 781

Stassun, K.G., Mathieu, R.D., Mazeh, T., & Vrba, F.J. 1999, AJ, 117, 2941

Stassun, K.G., Mathieu, R.D., Vrba, F.J., Mazeh, T., & Henden, A. 2001, AJ, 121, 1003

Stassun, K.G., & Terndrup, D. 2003, PASP, in press

Terndrup, D.M., Stauffer, J.R., Pinsonneault, M.H., Sills, A., Yuan, Y., Jones, B.F., Fischer, D. & Krishnamurthi, A. 2000, AJ, 119, 1303

Tinker, J., Pinsonneault, M. & Terndrup, D. 2002, ApJ, 564, 877

Vogel, S.N. & Kuhi., L.V. 1981, ApJ, 245, 960


format: text

text:
# Envelope Emission in Young Stellar Systems: A Sub-Arcsecond Survey of Circumstellar Structure

## 1 Introduction

Stars form deep within dense molecular cores which provide the raw material for star formation. Their structure and kinematics determine how the formation proceeds and the eventual outcome: single star, binary, or multiple system. For isolated low mass stars, we have a good overall picture of the process and well developed models for an idealized collapse to form a single star system. At a broad-brush level, these idealized collapse models are able to provide a good understanding of the process, despite shortcomings with respect to stellar multiplicity and core geometry: dense cores are roughly spherical and stellar multiplicity typically occurs on scales much smaller than scales characteristic of cores. This paper utilizes high resolution observations of the dust continuum emission from cores around young stars to study the morphology of the inner cores, or envelopes, around young systems. The primary aim is to compare the inferred density structure to that expected from collapse models.

The solution to the collapse of an isothermal sphere can be thought of as a continuum of solutions (Whitworth & Summers, 1985, but also see Hunter 1986) bracketed by the two standard solutions: the Larson-Penston solution, or LP solution, put forth by Larson (1969) and Penston (1969) and expanded by Hunter (1977) and the Shu solution proposed by Shu (1977). The LP solution is valid for the isothermal spherical collapse of an initially uniform medium. The core density profile is characterized by a uniform density central region surrounded by a $`\rho r^2`$ envelope. As collapse progresses, the uniform density region shrinks until the entire sphere follows $`\rho r^2`$; continued collapse sets up an $`r^{3/2}`$ central density profile around the central forming star (Hunter, 1977). In the other line of heritage, Shu (1977) defined a class of self-similar collapse solutions based on an initial singular isothermal spherical density distribution, $`\rho =\frac{a^2}{2\pi G}r^2`$, where a is the local sound speed in the gas and G is the gravitational constant. The sphere collapses most rapidly at the center, leading to a collapse wave which moves outward in radius at the local sound speed. This is referred to as an “inside-out” collapse; the density profile inside the wave approaches a free-fall density profile, $`\rho r^{3/2}`$; outside the collapse wave, the profile remains $`\rho r^2`$. Both solutions have density profiles of $`\rho r^2`$ in the outer core and both tend toward $`\rho r^{3/2}`$ in the centers some time after the star begins forming.

The biggest differences between the two solutions involve the mass infall rate, infall velocity structure, and the density structure evolution. The Shu solution has a constant mass infall rate at all times; the mass in the central object is nearly equal to the mass within the collapse wave radius at any time; and material outside the collapse wave is nearly static because the collapse starts from a static isothermal sphere. The LP solution has collapse motions wherein all radii acquire a velocity of 3.3 times the sound speed, a mass infall rate that peaks at the beginning then asymptotically approaches the Shu value over time, and a maximum infall velocity that is not in the center, as in the Shu solution, but at a finite radius.

Looney, Mundy, & Welch (2000) (hereafter LMW) presented images of 24 young stellar sources with sensitivity to spatial scales from 0$`\stackrel{}{\mathrm{.}}`$5 to 50$`\mathrm{}`$ in the $`\lambda `$ = 2.7 mm continuum. A discussion of the data acquisition and images can be found in that paper. In this paper, we present modeling of the envelope emission for many of the embedded objects in LMW. Our primary goals in the modeling are to address three issues: (1) What constraints, if any, can we place on the power-law of the density and on the two collapse solutions? Does the power-law index resemble the isothermal sphere ($`\rho r^2`$) or the free-fall profile ($`\rho r^{3/2}`$)? Or something else entirely? For example, Ward-Thompson et al. (1994) found that starless clouds tended to have flat-topped density profiles, more like Bonner-Ebert spheres. (2) What constraints can we place on the circumstellar disks in these systems? With the highest resolution to date at these wavelengths, we will be able to place limits on the size of the embedded disk. Can the data be fit with or without a central circumstellar disk? (3) What constraints can we place on the inner and outer radii of the envelopes?

## 2 The Observational u,v Data

To compare the observations and models, we analyze the data in u,v space where the data are not effected by the CLEAN algorithm or u,v sampling. The interferometric data are binned in u,v annuli around the source location from LMW and averaged vectorially. The resulting u,v amplitudes for the bins are shown in Fig. 1 with $`\mathrm{log}`$(u,v distance) versus $`\mathrm{log}`$(amplitude) plots. The displayed error bars are statistical error bars based on the standard deviation of the mean of the data points in the bin with a minimum of 10%, reflecting the uncertainty in the overall calibration. In the cases where the binary systems were separated by more than 10$`\mathrm{}`$, the companion sources were subtracted out of the u,v data using the large scale images of LMW. The new u,v data were remapped to confirm that the large scale emission from the companion sources was not detected. Although there may exist some residual of small-scale emission from the companion envelope in the u,v data, vector averaging in u,v annuli will minimize its contribution.

For most sources, the circular symmetry of the observed emission and the lack of significant internal structure (see images of the Class 0 sources in LMW) justifies the simplicity of a spherical model. The spherical envelope can be characterized by two radial distributions: $`\rho (r)`$ and $`T(r)`$. These quantities are typically and reasonably assumed to be power-laws in radius (Adams, Shu, & Lada, 1988; Beckwith et al., 1990; Adams, Emerson, & Fuller, 1990; Keene & Masson, 1990; Beckwith and Sargent, 1991; Terebey, Chandler, & André, 1993) $`\rho (r)=\rho _o(r/r_o)^p`$ and $`T(r)=T_o(r/r_o)^q`$.

What shape do we expect for the classic envelope power-law emission model in the u,v plane? One can get a sense from simple consideration of the interdependence of the radial brightness distribution in the image and u,v planes. A power-law radial dependence of the emission in the image plane, $`B(r)`$, transforms to a power-law dependence of the visibility amplitudes on u,v distance, $`V(s)`$. The index of the power-law in u,v distance is related to the index of the power-law in radial distance. Specifically, as long as the emission is optically thin and in the Rayleigh-Jeans limit, and the envelope infinite in size,

$$B(r)r^{(p+q)+1}V(s)s^{(p+q3)}$$

where $`s`$ is u,v distance (cf. Looney, 1998). Note that the slope is only sensitive to the sum of the temperature and density power-law indices. In the simplest terms, the large u,v spacings convey information about the inner part of the envelope and the short u,v spacing reflect the structure of the outer portion of the envelope. If the envelope is truncated at a finite radius, the cutoff in radius is equivalent to the convolution of a modified first-order Bessel function with the power-law in the u,v plane, causing a flattening of the slope at short u,v spacings and a ringing effect in the u,v distance versus amplitude plot. This ringing arises from the sharp edged envelope which is probably not physical, but the short spacing turnover is robust to the slope of the cutoff.

The bottom right panel of Fig. 1, shows two example simple power-law models– a finite envelope and compact circumstellar disk at the distance of Taurus. The envelope has an outer radius of 5000 AU, a density profile of $`\rho r^2`$, and a temperature profile of $`Tr^{0.4}`$. The disk has an outer radius of 100 AU, a surface density profile of $`\mathrm{\Sigma }r^{1.5}`$, and a temperature profile of $`Tr^{0.5}`$. The disk model shows that the flux remains constant with increasing u,v distance until the u,v distance is large enough to resolve the radial structure in the source. For the envelope, its finite size also creates a rounding off of the curve at very small u,v spacings, as the source becomes unresolved. When this happens the flux transitions to follow the power law with index $`p+q3`$. However, since the envelope is larger than a disk, the related power-law of the sum of the temperature and density is more quickly seen at shorter u,v spacings.

Comparing the panels in Fig. 1, it is immediately obvious that the data and models are not well represented by a single power-law extending to all radii. The most significant deviation can be seen at short u,v distances where the flux levels off in most sources, indicating finite sized envelopes. In two cases, IRAS 4 A and IRAS 16293-2422, the shortest u,v data, corresponding to the largest spatial scales, rises abruptly. For these two objects, it is most likely that there is some contamination from the very large scale cloud within which the sources are embedded. In addition, IRAS 2 B and IRAS 4 C do not follow the expected profile of an envelope dominated source and may be either very compact envelopes or circumstellar disks. This was also suggested by their near unity 5 k$`\lambda `$/50 k$`\lambda `$ ratio in LMW. As the very smallest spatial scales (the largest u,v distances) are probed, there are two separate profiles that stand out. Some sources (e.g. IRAS 4A or L1448 B) are consistently dropping in amplitude with increasing u,v distance, suggesting that the power-law is being traced down to very small scales. Others (e.g. SVS 13 A or IRAS 2A) flatten toward a constant flux level at very small scales, perhaps suggesting that there is an unresolved point source at the center of the envelope, arguably a circumstellar disk.

## 3 Sources and u,v Data

### 3.1 L1448 IRS3

L1448 IRS3 (Bachiller & Cernicharo, 1986) is comprised of three distinct sources in the $`\lambda `$ = 2.7 mm continuum (Terebey, Chandler, & André, 1993; Terebey & Padgett, 1997, Fig. 7a LMW). Sources A and B, a proto-binary system (separation $``$ 7$`\mathrm{}`$), share a common large-scale envelope. Only source B provides sufficient signal-to-noise in the u,v plane to be modeled effectively. However, the data for source C in Fig. 1 drops sharply in amplitude with u,v distance implying that the source is extended and dominated by the envelope component.

For the modeling, we assume that the circumbinary envelope is centered on source B. The total luminosity of sources A and B is estimated to be 8 L(Barsony et al., 1998); we assume that the luminosity of source B is 6.8 L, based on the ratio of fluxes of sources A and B at $`\lambda `$ = 2.7 mm (LMW). Fig. 1 displays the u,v data for L1448 IRS3 B, and the data are smooth with a well defined slope of -0.3 within the inner 40 k$`\lambda `$ (corresponding to $`p+q`$ 2.7, see §2). At 40 k$`\lambda `$, the curve transitions into a steeper slope of -1.1 (corresponding to $`p+q`$ 1.9, see §2).

### 3.2 NGC 1333 IRAS 2

The young, binary system NGC 1333 IRAS2 (Jennings et al., 1987; Blake, 1997, LMW) has a separation of $``$ 31$`\mathrm{}`$. In our high resolution image (Figure 8, LMW), the envelope of source A appears to be completely resolved out, and the remaining emission is consistent with a point source. Source B (Fig. 9, LMW) appears to be slightly resolved, and we concluded from Fig. 17 in LMW that source B is a very compact source. Overall, there may be an extended envelope that surrounds both of the sources, but our data can not place useful limits on this structure. The luminosity for sources A and B is estimated to be 40 L(Jennings et al., 1987). The luminosity for source A is assumed to be 30 L, based on the ratio of fluxes of sources A and B at $`\lambda `$ = 2.7 mm (LMW).

Fig. 1 displays the u,v data for NGC 1333 IRAS 2 A, after subtraction of source B from the u,v data. The curve has a well defined slope of -0.4 within the inner 12 k$`\lambda `$ (corresponding to $`p+q`$ 2.6, see §2). Unlike L1448 IRS3 B, the visibility at larger u,v distance flattens to a constant flux, rather than sharply descending. This indicates that the model point source component will play an important role.

As suggested in Fig. 17 of LMW, IRAS2 B is a compact source; the amplitude does not significantly drop off until around 100 k$`\lambda `$. It seems likely that this source has either a small, compact envelope or has a circumstellar disk with very little envelope emission. The slight dip in the visibility at 20 k$`\lambda `$ may be an artifact from the subtraction of source A from the u,v data.

### 3.3 SVS 13

SVS 13 (Strom, Vrba, Strom, 1976) has at least four young stellar systems (LMW). Besides the large-scale envelope that may surround all of the sources, A and B are embedded in individual envelopes, and, in the case of source A, the separate envelope is probably a circumbinary envelope that enshrouds both A1 and A2 (separation of $``$5$`\mathrm{}`$). In the $`\lambda `$ = 2.7 mm high resolution images (Figs. 10 & 11, LMW), SVS 13 A1 and B are both amorphous structures. The total luminosity for sources A, B, and C is 45 L (Jennings et al., 1987); we assume that the luminosity for source A and B are 19 and 22 L, respectively, based on the ratio of fluxes of sources A and B at $`\lambda `$ = 2.7 mm (LMW).

For source A, Fig. 1 shows the u,v data flatten to a constant in the outer u,v spacings. The curve has a well defined slope of -0.2 within the inner 8 k$`\lambda `$ (corresponding to $`p+q`$ 2.8, see §2). At 8 k$`\lambda `$, there is an inflection point, and the curve transitions to a steeper slope of -0.8 (corresponding to $`p+q`$ 2.2, see §2).

In Fig. 1, the u,v data are shown for NGC 1333 SVS 13 B. There is a slight “bump” in the curve at 80 k$`\lambda `$, which corresponds to a fringe of around 3$`\mathrm{}`$. This excess in amplitude may be a beating effect from incomplete subtraction of one of the other sources, or the excess may be from a non-spherical symmetric component of source B on a size scale of $`3\mathrm{}`$. In either case, the 80 k$`\lambda `$ u,v point cannot be fit with a simple-model, increasing the reduced $`\chi ^2`$ and biasing the amplitude upward in the outer parts of the u,v plane. The slope of the u,v data from 5 k$`\lambda `$ to 35 k$`\lambda `$ is -0.6 (corresponding to $`p+q`$ 2.4, see §2).

### 3.4 NGC 1333 IRAS4

NGC 1333 IRAS4 is comprised of at least four distinct young stellar objects (Sandell et al., 1991; Lay, Carlstrom, & Hills, 1995, LMW). Source A, the northern source, is a 1$`\stackrel{}{\mathrm{.}}`$7 binary system that shares a common circumbinary envelope (Lay, Carlstrom, & Hills, 1995, LMW). Source C is a compact source that has a brightness distribution that more resembles an optical/IR source than its IRAS 4 companions (LMW). The total luminosity for sources A and B is assumed to be 21 L (Jennings et al., 1987), and the luminosity for source A and B is assumed to be 16 and 5 L, respectively, based on the ratio of fluxes of sources A and B at $`\lambda `$ = 2.7 mm (LMW).

In Fig. 1, the u,v data are shown for NGC 1333 IRAS 4 A1— the brighter source at $`\lambda `$=2.7mm in the 1$`\stackrel{}{\mathrm{.}}`$7 binary. We did not subtract out the binary companion source A2; running test models with a second fixed point source at the location of source A2 does not significantly alter the fits. The shortest u,v spacing data point in Fig. 1 is excessively high. The best explanation is that the data are beginning to pick-up very large-scale structure, perhaps from the cloud. Since we are not including such a large-scale envelope in this modeling, we only use the shortest u,v spacing point as a maximum constraint in the models. In that case, the data has a well defined slope of -0.3 in the 3 k$`\lambda `$ to 35 k$`\lambda `$ range (corresponding to $`p+q`$ 2.3, see §2) and a slope of -0.7 in the 22 k$`\lambda `$ to 125 k$`\lambda `$ range (corresponding to $`p+q`$ 1.9, see §2).

In Fig. 1, the u,v data are shown for NGC 1333 IRAS 4 B. The data has a well defined slope of -0.2 in the 5 k$`\lambda `$ to 90 k$`\lambda `$ range (corresponding to $`p+q`$ 2.8, see §2). Unlike the data from IRAS 4 A, the slope in the outer u,v spacings is particularly steep. In Fig. 1, the u,v data are shown for NGC 1333 IRAS 4 C, after subtraction of sources A and B from the u,v data. As stated in LMW, the source is mostly compact with nearly constant u,v amplitude from 2 k$`\lambda `$ to 80 k$`\lambda `$. As shown in Fig. 1, IRAS 4 C becomes resolved for u,v distances longer than 80 k$`\lambda `$, but there are not enough data in this region to constrain the density power-law.

### 3.5 VLA 1623

VLA 1623 is the prototype Class 0 source (Andrè, Ward-Thompson, & Barsony, 1993). However, our high resolution $`\lambda `$ = 2.7 mm observations (Figure 15, LMW) suggest that the system is actually a binary system with a separation of 1$`\stackrel{}{\mathrm{.}}`$1. The u,v data were binned around the center of the system, between the two point sources. The first three u,v spacing data points are amplitude biased by large scale emission from the nearby regions of SMM1 and SMM2 (Ward-Thompson et al., 1989; Andrè, Ward-Thompson, & Barsony, 1993). From 7 k$`\lambda `$ to 25 k$`\lambda `$, the circumbinary envelope structure dominates the u,v data, with a shallow slope of -0.1 (corresponding to $`p+q`$ 2.9, see §2). Although, the exact index could vary significantly based upon embedded point source(s), inner radius, and outer radius. At u,v spacings greater than 30 k$`\lambda `$, the two circumstellar regions beat against each other, making it difficult to estimate the density power-laws.

We attempted to fit the u,v data with two point sources of flux equal to the peak flux in Figure 15 of LMW. However, this model overestimated the flux in the outer u,v spacings. A correct model requires two circumstellar disks (or point sources) embedded within two circumstellar envelopes, and perhaps a larger circumbinary envelope. Since this type of model has too many free parameters to be constrained by the current data, VLA 1623 was not modeled further.

### 3.6 IRAS 16293-2422

IRAS 16293-2422 is a deeply embedded binary with two molecular outflows (Walker et al., 1986; Wootten, 1989; Mundy et al., 1992). In high resolution observations at $`\lambda `$ = 2 cm, source A has two peaks, A1 and A2 (Wootten, 1989). In the $`\lambda `$ = 2.7 mm high resolution image (Figure 16, LMW), source A and B are clearly detected. Source A, the most extended object in the survey, appears slightly elongated along the position angle of the $`\lambda `$ = 2 cm sources. The two sources, which may have individual circumstellar envelopes and disks, are embedded within a circumbinary envelope.

The u,v data were binned around source B, since source A was known to be very extended. Like NGC 1333 IRAS 4 A, the first u,v data point is biased high, due to a large scale structure. There is a very steep slope of -1.7 from 8 k$`\lambda `$ to 20 k$`\lambda `$, which indicates a small $`p+q`$ value of 1.3. For the u,v values greater than 20 k$`\lambda `$, the beating effect of the two sources dominates the curve. We attempted to fit the u,v data with two point sources of flux equal to the peak flux in (Figure 16, LMW). As in the case of VLA 1623, this model overestimated the flux in the outer u,v spacings. A more complicated model is required, but these data would not constrain the model; thus, IRAS 16293-242 was not modeled further.

## 4 Modeling of the Emission

The simple power-law discussed in §2 provides the first insights into the behavior of the brightness distribution in the u,v plane. However, detailed physical interpretation of the data requires radiative transfer modeling. Unfortunately, a number of assumptions must be made about the dust properties and source morphology in order to simplify the model and reduce the number of free parameters; they make the model calculation easier but also reflect the limitations of the observational data and of our knowledge of dust and source properties. In the spirit of modest complexity, our basic model is a spherical dust envelope with an embedded emitting point source, representing an unresolved circumstellar disk.

### 4.1 Dust Properties

Dust emissivity is dependent on the dust grain composition, size, and distribution (c.f. Krügel & Siebenmorgen, 1994; Pollack et al., 1994), all of which can vary significantly with environment. Since grain emissivity affects grain thermal equilibrium and the observable emission, the determination of the dust temperature, column density, and optical depth from the observed emission is explicitly dependent on knowledge of the underlying dust properties.

For the broad dust description we adopt the model from Wolfire & Cassinelli (1986) which uses a MRN (Mathis, Rumpl, & Nordsieck, 1977) dust grain distribution (“bare” graphite plus silicate grain distribution) and the Draine & Lee (1984) optical constants to describe grain properties. However, for this grain model the long wavelength dust emissivity is assumed to follow a power-law ($`\kappa _\nu \lambda ^\beta `$) with $`\beta `$ = 2. Studies have suggested that the dust emissivity in circumstellar disks and dense cores has a $`\lambda ^1`$ dependence at submillimeter wavelengths (Weintraub, Sandell, & Duncan, 1989; Beckwith et al., 1990; Beckwith and Sargent, 1991) rather than the $`\lambda ^2`$ dependence found in calculations for bare grains (cf. Draine, 1990) and as observed in the large scale ISM (Reach et al., 1995). Accordingly, we have modified the long wavelength characteristics of our dust model ($`\lambda >100\mu `$m) to assume a $`\lambda ^1`$ emissivity. This was done by adopting a value of $`\kappa _\nu `$(1 mm) = 0.0033 cm<sup>2</sup> g<sup>-1</sup> and following that back to shorter wavelengths where the $`\lambda ^1`$ curve intersected the bare grain emissivity curve. This hybrid model preserves the optical and infrared properties of the MRN dust grain model, while forcing the long wavelength $`\lambda ^1`$ behavior and $`\kappa _\nu `$(2.7 mm) = 0.009 cm<sup>2</sup> g<sup>-1</sup>. Recent works looking at the $`\kappa _\nu `$(1.2 mm) appropriate for dense molecular clouds and circumstellar disks yield reasonable consistent values, $`\kappa _\nu `$(1.2 mm) = 0.004 cm<sup>2</sup> g<sup>-1</sup> with range of a factor of 2 (Kramer et al., 1998). Our mass determinations, all based on $`\lambda `$ = 2.7 mm data, depend roughly linearly on the assumed $`\kappa `$. We assume that dust properties are not a function of radius in the envelope.

### 4.2 Temperature and Density Profiles

As discussed in §2, the young protostellar envelope can be described by two radial power-law distributions for density and temperature: $`\rho (r)=\rho _o(r/r_o)^p`$ and $`T(r)=T_o(r/r_o)^q`$. A power-law is a reasonable first-order approximation for the density profile in many situations, as discussed in §1. Also, in a realistic system the density distribution has an inner and outer cutoff radius, $`R_{in}`$ and $`R_{out}`$. These cutoffs are physically motivated by the presence of the circumstellar disk which truncates the inner envelope and the outer radius which is set by either a companion forming star or merging with the extended cloud.

Although a power-law is a reasonable assumption for the density profile, a radial power-law temperature profile is more problematic for our application. The temperature profile in an optically thin, dust envelope heated by a central star will have a power-law dependence on radius and stellar luminosity ($`L_{}`$) if dust opacity has a power-law dependence on frequency: $`T(r)=T_o(r/r_o)^q(L_{}/L_o)^{\frac{q}{2}},`$ where $`T_o`$ is the dust temperature at the radius $`r_o`$ for a stellar luminosity of $`L_o`$. The index $`q`$ depends on the dust opacity power-law index, $`q=\frac{2}{4+\beta }`$. For values of $`\beta `$ from 0 to 2, the temperature power-law index is between 0.33 and 0.5.

The use of simple power-laws in density and temperature allow for a quick interpretation of an observed emission profile, as shown in §2. However, the temperature profile can diverge from a simple power-law as the envelope becomes optically thick at the primary wavelengths of energy transport. For a centrally peaked density profile, the envelope can become optically thick over the inner 100’s of AU, resulting in a steeper inner temperature profile. To account for the effects of this heating we have calculated self-consistent dust temperature profiles using the code of (Wolfire & Cassinelli, 1986, WC hereafter).

The WC code assumes a central heating source embedded within a spherical dust envelope. The central source is characterized by a luminosity and an effective temperature. For the effective temperature, we used 10,000 K, which is both consistent with the effective temperature from models of zero age low mass stars (Hejlesen, 1980; Girardi et al., 2000) and the temperature derived for T Tauri stars to explain the veiling continuum (Hartigan, Edwards, & Ghandour, 1995). The dust envelope is specified by an outer radius, the power-law density index, the density at the outer radius, and the destruction temperature of the dust, which specifies the inner radius. Given these parameters, the WC code self-consistently calculates the dust temperature profile by conserving luminosity at all radii.

Figs. 2 and 3 display plots of the temperature profile for various values of the density power-law index p and envelope mass, respectively. As can be seen in Fig. 2, the density power-law index can have significant impact upon the temperature profile at radii less than 200 AU. Indeed, the temperature radial distribution can depart significantly from the optically-thin assumption of $`Tr^{0.4}`$, the line at the bottom of the figure. As the density power-law index increases, the inner opacity increases, and the interior material can not radiate away its luminosity efficiently, causing a steep increase in temperature. At the outer radii, the emission escapes easily and the temperature approaches the optically thin slope although the absolute temperature can differ by as much as 50%. A similar effect occurs with increasing mass (Fig. 3). As the mass of the envelope is increased the radius at which the temperature approaches the optically thin limit increases. The temperature at inner radii can be as much as a factor of three times the optically thin solution.

The WC code shows that the optically thin temperature assumption is incorrect for envelope masses $`>`$ 0.1 M and density profiles steeper than p = 1. The largest changes occur at small radii, $`<`$ 200 AU, (0$`\stackrel{}{\mathrm{.}}`$6 at the distance of Perseus), a region beginning to be probed by our observations. Beyond 200 AU all of the temperature profiles tend toward the optically thin $`Tr^{0.4}`$ but with a shift in the absolute temperatures. Finally, increasing luminosity for fixed envelope properties approximately scales the temperature profiles by $`(L/L_{\mathrm{}})^{0.25}`$ as found by Wilner, Welch, & Forster (1995).

## 5 Modeling Details and Results

For each model, the dust temperature distribution is derived self-consistently with the luminosity; the envelope emission is calculated as an image; the central point source flux, attenuated by the envelope, is added; and the model is multiplied by the BIMA primary beam. The resulting image is Fast Fourier Transformed and sampled with the same u,v spacings as the data. The data and model are both vector-averaged into u,v distance bins (annuli in the u,v plane) with the phase center at the source position given in LMW. A reduced $`\chi ^2`$ of the fit is calculated using the model, observed flux amplitudes, and the $`\sigma `$ of the data in each radial bin.

The model has three fixed quantities (luminosity, dust emissivity, and distance), as well as five degrees of freedom: power-law index of the density ($`p`$), total mass of the envelope ($`M`$), inner cutoff radius ($`R_{in}`$), outer cutoff radius ($`R_{out}`$), and point source flux ($`S_p`$). For each model, the $`\chi ^2`$ was minimized with respect to the envelope mass to produce the best fit to the data. We were only able to fully model six of the eleven sources due to either a lack of signal-to-noise (L1448 C), inability to constraint fit (IRAS 2B and IRAS 4C), or dominant binarity (VLA 1623 and IRAS 16293-2422).

In this section, we will first examine the constraints that the data can place on the power-law index of the density, then, in more detail, the interesting $`p=1.5`$ and $`p=2.0`$ power-law indices, and finally, specific fitting of the Larson-Penston and Shu collapse solutions to the data. In the latter case, the envelope parameters are uniquely determined by the sound speed and the time since the formation of a centralized source, t = 0 in the standard collapse scenario (cf. Larson, 1969; Penston, 1969; Shu, 1977; Hunter, 1977).

### 5.1 Results for Power Law Envelopes

The first question to examine is how well does the data constrain the models parameters? The most interesting parameter is the density power-law index since it is of central interest to collapse models. To answer this, we explored parameter space including: $`p`$ from 0.5 to 2.9 in steps of 0.2, $`R_{out}`$ from 2000 AU to 10000 AU in steps of 1000 AU, $`R_{in}`$ of 10 AU and 100 AU, and a central point source flux ($`S_{ps}`$) in 1$`\sigma `$ steps ($`\sigma `$ = RMS noise in the last u,v distance amplitude bin) starting with no point source, and increasing to a maximum of the amplitude in the last u,v distance bin.

We were able to find acceptable fits for a range of parameter values, but the data do not always well constrain $`p`$, even with the highest signal to noise (S/N) sources (Table 1). Fig. 4 shows the best fits for NGC 1333 IRAS2 A, a modest S/N example, for different values of $`p`$; the data can be acceptably fit by $`p`$ in the range 0.5 to 2.3. Fig. 5 presents the range of the acceptable models for the six sources. The fits follow general trends, and there are typically a number of models with good fits within a region of parameter space. Although most of the envelope models are fit with density power-law indices between $`p`$ = 1.5 and 2.0, some sources prefer steep power-law indices. In particular, the three sources with the highest S/N only have acceptable fits for $`p>1.7`$.

In general, but especially in modest S/N cases, a range of fits are possible as a result of (1) the outer radius of the envelope is typically weakly constrained for models with $`p>1.5`$ since the combination of the density and temperature laws ($`p+q`$) cause the total emission to grow slowly with increasing radius, (2) an increasing point source flux can be compensated by increasing the inner radius cutoff of the envelope, and (3) the outer radius and power-law index can conspire to fit the inner u,v region. If one or more of these parameters is constrained, the range in $`p`$ can become much smaller. Nonetheless, the higher S/N systems do constrain the power-law in the envelope, and these data can be used to make model significant statements.

### 5.2 Results for $`p=1.5`$ and $`p=2.0`$

Now that we have an idea of how constrained the $`p`$ values are by our data, we can explore in more detail $`p=1.5`$ and $`p=2.0`$ models. As discussed in §1, we expect a density structure of $`\rho `$ r<sup>-2</sup> at young ages evolving from the inside-out at later times to $`\rho `$ r<sup>-3/2</sup>. The other parameters are explored over a finer grid: inner radii ($`R_{in}`$) of 10, 30, 50 and 100 AU, outer radius ($`R_{out}`$) from 1000 AU to 10000 AU, in steps of 1000 AU, and the same step size for the embedded point source as in the general cases.

Fig. 6 shows the best fit $`p=1.5`$ and $`p=2.0`$ models along with the u,v data for the six modeled sources; the best fit parameters are listed in Table 2. All six of the sources can be fit with a density power-law index of $`p=2.0`$; in half of the systems, those with the highest S/N ratios, $`p=1.5`$ power-law models can be rejected at the $`>95`$% confidence level. What do these fits tell us about the nature of the systems? In order to understand the constraints that these higher S/N systems place on descriptions of the envelope collapse, it is necessary to examine the Larson-Penston (LP) and Shu solutions in more detail. Fig. 7 shows the density profiles for the two solutions at time zero and 10<sup>4</sup> years, assuming a sound speed of 0.2 km/s (e.g. Myers & Benson, 1983). There are two obvious differences: the LP solution has a factor of 4.4 higher density and the inner envelope region evolves more quickly into a $`p=1.5`$ power-law density profile.

In the Shu model, an expansion wave travels from the inner region outwards— the “inside-out” collapse. The head of the expansion wave is roughly where the $`\rho r^2`$ region transitions to the $`\rho r^{3/2}`$ region, the kink near 400 AU radius in Fig. 7. The head of the expansion wave is only dependent upon the time since the formation of the singular core and the local sound speed. Thus, for a sound speed, $`a=0.2`$ km/s, the head of the expansion wave will be located at 400 and 4000 AU for ages of $`10^4`$ and $`10^5`$ years, respectively. Note that in Fig. 7 there is a significant interface region (100 to 400 AU) in the Shu solution where the density slope is shallower than the $`p=1.5`$ power-law, whereas the LP solution gently increases from $`p`$ = 1.5 to 2.

In the LP solution, there is no “inside-out” expansion wave since the initial velocity profile is already collapsing, not at rest as in the Shu solution. Nonetheless, there is an infall scale for the interface between the $`\rho r^2`$ and $`\rho r^{3/2}`$ regions. Since the LP infall velocity is 3.3 times that in the Shu solution, the effective spatial scale is 3.3 times as large for the same sound speed. Thus, the radius at which one expects a p = 1.5 is $``$ 1300 and 13000 AU for ages of $`10^4`$ and $`10^5`$ years, respectively. In addition, the larger mass infall rate creates a higher density envelope as well as a larger mass protostar earlier in the collapse; the mass of the envelope at the onset of the core formation is 4.4 times larger than in the Shu solution.

In order to match the $`\rho r^2`$ power-law fits, the infall regions in these models must not dominate the radial density profiles within a few hundred AU. This implies that these systems must be young, which is consistent with their classification as Class 0 sources. The youthfulness of the systems is also supported by the dynamical age estimates for the outflows. For example, L1448 IRS3 has the largest outflow, 14.5$`\mathrm{}`$ (Bally et al., 1997). For an assumed jet speed of $``$ 200 km/s, the minimum age of L1448 IRS3 is $``$ 6,200 yrs. The corresponding times for the other systems are shorter, but it must be stressed that these are minimum ages as the outflow region should propagate out at a fraction of the basic jet speed. In addition, all current outflow mechanisms require the presence of a stellar-size object which is not present until after the collapse of the first hydrostatic core. This collapse occurs when about 0.05 M has accumulated within 5 AU, effectively extending the age of these systems with respect to the classic self-similar solutions (e.g. Masunaga, Miyama, & Inutsuka, 1998).

A simple interpretation of these arguments for the infall region size would suggest that our fits to the u,v data require an age of $`<`$ 5000 yrs for the Shu model, and an age of $`<`$ 1500 yrs for the LP model (assuming $`a=0.2`$ km/s and ignoring any differences in absolute density). At such young ages, the luminosity of the system is also a constraint since it is determined by the infall rate and the mass of the stellar object. In the above example of L1448 IRS3, which has an estimated luminosity of 6.8 L, the Shu solution underestimates the luminosity; a 75% increase in the sound speed is needed to match the luminosity for an age of 6,200 years. Of course, if the sound speed is increased, the expansion wave moves faster and the radius of the transition region increases; this pushes allowed solutions to even earlier times. The luminosity is not as much of a problem for the LP solution since its mass accretion is nearly 46.9 times larger at early times.

### 5.3 Results with Fits using the LP and Shu Collapse Solutions

As the above reasoning suggests, if one wants to compare the LP and Shu solutions, the age and sound speed must be included with a complete description of the two solutions. We have performed modeling of the data using the LP and Shu solutions for the density structure of the envelope as a function of time and sound speed. Reasonable fits for both models can be achieved for an age around 1000 years, with $`a`$ = 0.46 km/s and 0.22 km/s for the Shu and LP models, respectively. For example in the source IRAS 4A, acceptable fits were found for ages around 1000 years for the Shu model and 1000 - 2000 years for the LP model. These results emphasize that the expected difference in the extent of the infall region of the two solutions at young ages is neutralized by the high sound speed necessary for the Shu solution to fit the physical densities inferred from the observations. We stress that we do not propose that the sources are 10<sup>3</sup> yrs old, but show that the strict fit to the data require the standard models to be that age.

In our modeling, the LP solution provides an overall better match to the data, but neither solution is the ideal fit. The computational simulations of the collapse of isothermal spheres demonstrate that it is incorrect to think of the envelopes of young stars as the realization of either collapse model because the collapse behavior ranges between the LP solution at early times and the Shu solution at later times (cf. Hunter, 1977; Foster & Chevalier, 1993). With that in mind, our modeling of small scale structure in the envelopes of young stars suggests that additional physics is necessary to model the structure of the envelope. Most likely, inclusion of magnetic fields (e.g. Mouschovias, 1976) and/or turbulence (e.g. Gammie & Ostriker, 1996) is important to understand properly the structure in these systems. Turbulent creation of the initial mass concentration which then evolves into the core and forming star (Ostriker, Gammie, & Stone, 1998) is one example of an intriguing alternative model.

### 5.4 Alternative Assumptions in the Model

Is there another way in which a $`p=1.5`$ power-law can mimic a steeper power-law in our fits? Two of our simplifying assumptions can have an impact on the fitted power-law index of the density: the dust opacity and the geometry. We adopted a constant dust opacity in our models. However, dust properties can change with environment (e.g. Gehrz, 1989; Weintraub, Sandell, & Duncan, 1989; Henning, Michel, & Stognienko, 1995). There are several grain alterations that may explain an increased dust opacity in circumstellar regions, such as chemical evolution (Begemann et al. 1994; van Dishoeck & Blake 1988), formation of dirty ice mantles (Draine 1985; Henning, Chini, & Pfau 1991; Preibisch et al. 1993), altering of the grain geometry (long “needle-like” grains; Wright 1982), or grain coagulation into fluffy grains (Wright 1987; Jones 1988; Bazell & Dwek 1990; Ossenkopf 1991; Stognienko, Henning, & Ossenkopf 1995). In any of these scenarios we may expect the dust opacity in the circumstellar region to become a function of radius: the outer regions could have grain properties similar to the interstellar medium while the inner, denser portion of the envelope is likely to have the most processed, perhaps larger grains due to the short timescales for grain alteration.

Since this timescale depends on the density, the dust opacity could have a power-law dependence on radius. In the standard power-law envelope model, the emission is dependent on the optical depth, $`d\tau =\kappa _\nu \rho (r)dl.`$ With grain alteration, the optical depth could become the product of power-laws in density and dust opacity,

$$d\tau =\kappa _0\left(\frac{\nu }{\nu _o}\right)^\beta \left(\frac{r}{r_o}\right)^s\rho _o\left(\frac{r}{r_o}\right)^pdl,$$

where $`s`$ is the radial power-law index for optical depth. Thus, a radial dependency in the dust opacity could be indistinguishable from that in density and erroneously produce steeper density power-law indices in our simple model fits.

Another possible explanation for the steeper density profiles is non-spherical or more complicated geometries. All of the embedded systems are known to be multiple systems or driving large molecular outflows that are evacuating material out of the envelope, both of which may significantly affect the envelope structure. In order to explore the effect of geometry on our models, we stretched one of our envelope models such that the elongation ratio of the semi-major to semi-minor axis is 1.2. Fig. 8 shows that in this simple model, an elongation can make the density power-law index appear steeper than expected. Note however that the observed emission from the sources modeled here are roughly circular symmetric.

## 6 Conclusions

We have presented the first modeling of circumstellar envelope dust emission that is uniformly sensitive to core structure from sub-arcsecond out to 50 arcseconds. Other groups have done modeling of the density structure in the outer envelope (c.f. Chandler & Richer 2000; Shirley et al. 2000), but probing the inner envelope, where the collapse solutions differ, requires high spatial resolution. Even with our high spatial resolution, only high signal-to-noise data can provide significant constraints to the standard collapse models. With the models presented in this paper, we place some of the first constraints on the emission contributions from the envelope, inner envelope, and disk in the youngest protostellar sources in the nearest star forming regions. Morphologically, the u,v data of all the sources typically show extended objects that either appear to follow a power-law down to small scales (large u,v ) or level out at some constant residual value, presumably a disk.

When modeling the inner envelope with high resolution, the standard power-law temperature model is not valid, and a self-consistent, luminosity preserving technique is required. Using a self-consistent temperature model and a dust structure model consisting of a simple power-law density, inner and outer radius cutoffs, and an embedded point source to emulate an embedded circumstellar disk, we have fit the u,v data from our $`\lambda `$ = 2.7 mm sub-arcsecond survey (LMW), placing constraints upon the conditions in the early stages of star formation.

Specifically, we found that the density structure in L1448 IRS3 B, NGC 1333 IRAS 4 A, and NGC 1333 IRAS 4 B are best-fit with a power-law index of p = 2.0. However models of star formation (Larson-Penston and Shu solutions to the isothermal sphere) predict that a $`\rho r^{3/2}`$ density power-law should be established quickly in the collapse process. So neither of the two self-similar collapse solutions exactly fit our data, although the LP solution provides the overall better fit with more likelihood and allowing the observed luminosity. The strict Shu inside-out collapse solution fit to the data requires very young, unphysical ages ($`<`$ 1000 yrs) and very high sound speeds that in some cases are still not enough to match the luminosity. The LP solution easily fits the luminosity constraints, but the LP solution evolves too quickly into a $`\rho =r^{3/2}`$ density power-law. The strict LP solution fit to the data also requires very young ages 1000 to 2000 yrs that are non-physical: all of the modeled systems have central stellar objects and outflows. Age estimates based on the apparent outflow extent suggest typical ages for these systems of one to a few times 10<sup>4</sup> years.

It is clear that there is a fundamental problem with the time scale in these simple self-similar models; either the models are not accurately including all the physics important to evolution at early times or in a real system the t = 0 point is not well determined. Incorporation of magnetic fields and turbulence, on all scales, is likely to be important new pieces. Such an improvement of models will become even more critical in the future with the advent of higher resolution and higher sensitivity submillimeter and millimeter observations become available from the SMA, CARMA, and ALMA.

One of the primary strengths of this study is the ability of the interferometer to separate large scale emission from compact emission, allowing us to probe for a circumstellar disk component embedded within a circumstellar envelope. Our data solidly show that most of the emission, a minimum of 85% up to 100% arises from the circumstellar envelope. But, point source flux values in the range of 0 to 40 mJy are found for the best fits (Table 1). We can make a mass estimate for an embedded disk using the circumstellar disk of HL Tauri as a standard. A HL Tauri type disk (disk mass $``$ 0.05 M; LMW) has a flux of $``$ 100 mJy at the distance of Taurus (140 pc). If placed at the distance of Perseus (350 pc), the flux of HL Tauri would be 16 mJy. So the range of acceptable point source fluxes would represent circumstellar disk masses of 0 to 0.12 M. This is a small fraction of the circumstellar envelope mass, typically $``$ 1 M. In other words, since the disks are not over-massive compared to the Class I/II systems, the larger quantity of infalling material in the Class 0 stage must process through the disk and onto the protostar quickly.

We especially thank Mark Wolfire for all of his help and advice on his self-consistent radiative transfer code without which this paper could never have been realized. We thank Pedro Safier for discussions on cloud collapse and Eve Ostriker, Steve Lubow, Chris Hunter, and Fred Adams for useful discussions. LWL would also like to thank the Infrared Group at MPE and its Director Reinhard Genzel in particular. This work was supported by NSF Grants NSF-FD93-20238 and AST-9314847. LWL and LGM acknowledge support from NASA grant NAGW-3066.


text:
# Neutrino Superbeams and the Magic Baseline

## Abstract

We examine the sensitivity to $`\nu _\mu \nu _e`$ of a conceptual experiment with a neutrino superbeam incident on a Megaton-scale water Cherenkov detector over a ”magic” baseline $``$7300 km. With realistic beam intensity and exposure, the experiment may unambiguously probe $`\mathrm{sin}^22\theta _{13}`$ and the sign of $`\mathrm{\Delta }m_{31}^2`$ down to $`\mathrm{sin}^22\theta _{13}10^3`$.

Detecting the subdominant oscillation $`\nu _\mu \nu _e`$ on the ”atmospheric” scale of $`L/E`$ has emerged as a priority for long-baseline accelerator experiments. This is because the $`\nu _\mu \nu _e`$ and $`\overline{\nu }_\mu \overline{\nu }_e`$ probabilities are sensitive to yet-unknown parameters of neutrino mixing: the mixing angle $`\theta _{13}`$, the sign of the ”atmospheric” mass-squared difference $`\mathrm{\Delta }m_{31}^2`$, and the $`CP`$-violating phase $`\delta _{CP}`$ . However, extracting the values of these parameters from measured probabilities will encounter the problem of degenerate solutions . In particular, the asymmetry between $`P(\nu _\mu \nu _e)`$ and $`P(\overline{\nu }_\mu \overline{\nu }_e)`$ may arise from either the intrinsic $`CP`$ violation and the matter effect that is correlated with the sign of $`\mathrm{\Delta }m_{31}^2`$ . The degeneracies can be resolved by comparing the data taken with a shorter and longer baselines . Selecting the latter as the ”magic” baseline $`L_{\mathrm{magic}}7300`$ km will render this strategy particularly effective: for $`L=L_{\mathrm{magic}}`$, all $`\mathrm{\Delta }m_{21}^2`$-induced effects like $`CP`$ violation are predicted to vanish up to second order of the small parameter $`\mathrm{\Delta }m_{21}^2/\mathrm{\Delta }m_{31}^2`$ . Therefore, selecting $`L=L_{\mathrm{magic}}`$ may allow to uniquely determine $`\mathrm{sin}^22\theta _{13}`$ and the sign of $`\mathrm{\Delta }m_{31}^2`$, but not $`\delta _{CP}`$ which should be probed with a shorter baseline.

In this paper, we discuss a conceptual experiment that involves a neutrino ”superbeam” incident on a water Cherenkov detector over a magic baseline of $`L=7340`$ km111This is chosen to match the distance from Fermilab to Gran Sasso or from CERN to Homestake.. A water Cherenkov target is selected on the merit of good separation and spectrometry of electromagnetic showers , and is assumed to be a megaton-scale detector like UNO or Hyper-Kamiokande . In tuning the energy of the neutrino beam, one must take into account that the $`E_\nu `$-dependence of oscillation probability for $`L=7340`$ km is strongly affected by Earth matter: for $`\mathrm{\Delta }m_{31}^2>0`$, the matter effect  shifts the first maximum of $`P(\nu _\mu \nu _e)`$ down to $`E_\nu /\mathrm{\Delta }m_{31}^22.5\times 10^3`$ GeV/eV2 from the vacuum value of $`5.9\times 10^3`$ GeV/eV2. Assuming $`\mathrm{\Delta }m_{31}^2=0.003`$ eV2, the oscillation maximum is at $`E_\nu 7.5`$ GeV which conveniently matches the peak of $`\nu _\mu `$ flux in the ”Medium-Energy” (or PH2me) beam of Fermilab’s Main Injector, as designed for the NuMI–MINOS program . Therefore, this is selected as the model beam in our simulation. We assume $`1.6\times 10^{21}`$ protons on neutrino target per year, as expected upon the planned upgrade of Main Injector’s intensity . In the absence of oscillations, the beam will produce some 58 $`\nu _\mu `$CC (21 $`\overline{\nu }_\mu `$CC) events per 1 kton$`\times `$yr in the far detector with the $`\nu `$ ($`\overline{\nu }`$) setting of the focusing system.

At neutrino energies below 1 GeV, as in the proposed JHF–Kamioka experiment , $`\nu _e`$ appearance can be efficiently detected in a water Cherenkov apparatus by selecting 1-ring $`e`$-like events of the reaction $`\nu _eNe^{}X`$ that is dominated by quasielastics. (Here and in what follows, $`X`$ denotes a system of hadrons other than the $`\pi ^0`$, in which the momenta of all charged particles are below the Cherenkov threshold in water.) At substantially higher energies considered in this paper, using the 1$`e`$ signature of $`\nu _\mu \nu _e`$ is complicated by more background from the flavor-blind NC reaction $`\nu N\nu \pi ^0X`$: its cross section increases with $`E_\nu `$, and so does the fraction of $`\pi ^0`$ mesons whose $`\gamma \gamma `$ decays produce a single $`e`$-like ring in the water Cherenkov detector222This happens when the opening angle is too small for the two showers to be resolved .. In , we have demonstrated that $`\nu _e`$ appearance can be analyzed with less NC background by detecting the reactions $`\nu _eNe^{}\pi ^+X`$ and $`\nu _eNe^{}\pi ^0X`$ that involve emission of a charged or neutral pion333Here and below, corresponding antineutrino reactions are implicitly included.. We proceed to briefly describe the selections of these CC reactions, as formulated in .

The reaction $`\nu _eNe^{}\pi ^+X`$ is selected by requiring two rings in the detector, of which one is $`e`$-like and the other is non-showering and has a large emission angle of $`\theta _\pi >50^0`$. This is referred to as the ”$`e\pi `$ signature”. The selection $`\theta _\pi >50^0`$ is aimed at suppressing the NC reaction $`\nu p\nu \pi ^0p`$ in which the momentum of the final proton is above the Cherenkov threshold444This reaction may also be rejected by identifying relativistic protons by ring shape, as proposed in .. The residual NC background is largely due to the reaction $`\nu N\nu \pi ^0\pi ^\pm X`$ with two pions in the final state. The $`\nu _\mu `$CC background arises from the reaction $`\nu _\mu N\mu ^{}\pi ^0X`$ in which the muon is emitted at a broad angle. The $`\nu _\tau `$CC background arises from the dominant oscillation $`\nu _\mu \nu _\tau `$ followed by $`\nu _\tau N\tau ^{}\pi ^+X`$ and $`\tau ^{}e^{}\nu \overline{\nu }`$.

The reaction $`\nu _eNe^{}\pi ^0X`$ is selected by requiring either three $`e`$-like rings of which two fit to $`\pi ^0\gamma \gamma `$, or two $`e`$-like rings that would not fit to a $`\pi ^0`$. This is referred to as the ”multi-$`e`$ signature”. The NC background arises from the reaction $`\nu N\nu \pi ^0\pi ^0N`$ in which at least one of the two $`\pi ^0`$ mesons has not been reconstructed. Note that in the latter reaction the two $`\pi ^0`$ mesons are emitted with comparable energies, whereas in $`\nu _eNe^{}\pi ^0X`$ the $`e^{}`$ tends to be the leading particle. This suggests a selection based on the absolute value of asymmetry $`A=(E_1E_2)/(E_1+E_2)`$, where $`E_1`$ and $`E_2`$ are the energies of the two showers for the two-ring signature, and of the reconstructed $`\pi ^0`$ and the ”odd” shower—for the three-ring signature. In this paper, we use the selection $`|A|>0.6`$. The $`\nu _\tau `$CC background is largely due to electronic decays of $`\tau `$ leptons produced in association with a $`\pi ^0`$. The $`\nu _\mu `$CC background originates from CC events with a muon below the Cherenkov threshold and two $`\pi ^0`$ mesons in the final state, and is negligibly small.

In the simulation, the matter effect is accounted for in the approximation of uniform matter density along the neutrino path ($`\rho =4.3`$ g/cm3 for $`L=7340`$ km), which adequately reproduces the results of exact calculations for the actual density profile of the Earth . Relevant neutrino-mixing parameters are assigned the values consistent with the atmospheric and reactor data : $`\mathrm{\Delta }m_{31}^2=\pm 0.003`$ eV2, $`\mathrm{sin}^22\theta _{23}=1`$, and $`\mathrm{sin}^22\theta _{13}=0.01`$ (the latter value is ten times below the upper limit imposed in ). The simulation relies on the neutrino-event generator NEUGEN based on the Soudan-2 Monte Carlo , that takes full account of exclusive channels like quasielastics and excitation of baryon resonances.

The $`E_{\mathrm{vis}}`$ distributions of 1$`e`$-like, $`e\pi `$-like, and multi-$`e`$-like events are illustrated in Fig. 1, assuming $`\mathrm{\Delta }m_{31}^2>0`$ and incident neutrinos. Here, $`E_{\mathrm{vis}}`$ stands for the net energy of all $`e`$-like rings. Total background to the $`\nu _\mu \nu _e`$ signal is seen to be the greatest for 1$`e`$-like events, and therefore we drop these from further analysis. Combined $`E_{\mathrm{vis}}`$ distributions of $`e\pi `$-like and multi-$`e`$-like events are shown in Fig. 2 for either beam setting and either sign of $`\mathrm{\Delta }m_{31}^2`$. With equal $`\nu `$ and $`\overline{\nu }`$ exposures of 1 Mton$`\times `$yr, the oscillation signal reaches some 250 events for $`\mathrm{\Delta }m_{31}^2>0`$ and incident neutrinos, and some 140 events for $`\mathrm{\Delta }m_{31}^2<0`$ and incident antineutrinos.

The experimental strategy we adopt is to share the overall exposure between the $`\nu `$ and $`\overline{\nu }`$ running so as to equalize the expected backgrounds under the $`\nu _\mu \nu _e`$ and $`\overline{\nu }_\mu \overline{\nu }_e`$ signals, and then analyze the difference between the $`E_{\mathrm{vis}}`$ distributions for the $`\nu `$ and $`\overline{\nu }`$ beams. The motivation is that many systematic uncertainties on the background should cancel out in the difference555This is particularly important here, as the large dip angle of the neutrino beam ($`35^0`$) will rule out the construction of a ”near” water Cherenkov detector.. The $`\nu `$ and $`\overline{\nu }`$ backgrounds are approximately equalized by running 1.7–1.8 times longer in the $`\overline{\nu }`$ mode than in the $`\nu `$ mode (see Fig. 2). The difference between the $`E_{\mathrm{vis}}`$ distributions for the $`\nu `$ and $`\overline{\nu }`$ beams, assuming $`\nu `$ and $`\overline{\nu }`$ exposures of 1.0 and 1.8 Mton$`\times `$yr, is illustrated in Fig. 3. Depending on the sign of $`\mathrm{\Delta }m_{31}^2`$, this distribution shows either a bump or a dip at oscillation maximum with respect to the background that corresponds to $`\mathrm{sin}^22\theta _{13}=0`$.

In order to estimate the significance of the oscillation signal in Fig. 3, we vary the $`E_{\mathrm{vis}}`$ interval so as to maximize the ”figure of merit” $`F=(S_\nu S_{\overline{\nu }})/\sqrt{B_\nu +B_{\overline{\nu }}}`$. Here, $`S_\nu `$ and $`S_{\overline{\nu }}`$ are the numbers of $`\nu _\mu \nu _e`$ and $`\overline{\nu }_\mu \overline{\nu }_e`$ events falling within the $`E_{\mathrm{vis}}`$ interval, and $`B_\nu `$ and $`B_{\overline{\nu }}`$ are corresponding numbers of background events. We obtain $`F=+19.6`$ for $`\mathrm{\Delta }m_{31}^2>0`$, and $`F=20.8`$ for $`\mathrm{\Delta }m_{31}^2<0`$. Recalling that these figures refer to $`\mathrm{sin}^22\theta _{13}=0.01`$, we estimate that at 90% CL the sensitivity to either $`\mathrm{sin}^22\theta _{13}`$ and the sign of $`\mathrm{\Delta }m_{31}^2`$ will be maintained down to $`\mathrm{sin}^22\theta _{13}8\times 10^4`$. Still lower values of $`\mathrm{sin}^22\theta _{13}`$ may perhaps be probed with a neutrino factory in combination with a magnetized iron–scintillator detector . Note however that the experimental scheme proposed in this paper is based on proven technology and involves a multi-purpose facility  rather than a dedicated detector.

To summarize, we have examined the physics potential of an experiment with a neutrino superbeam that irradiates a Megaton-scale water Cherenkov detector over the ”magic” baseline $``$7300 km. With realistic beam intensity and exposure, the experiment may probe $`\mathrm{sin}^22\theta _{13}`$ and the sign of $`\mathrm{\Delta }m_{31}^2`$ down to $`\mathrm{sin}^22\theta _{13}`$ values below $`10^3`$. Thus obtained values of these parameters, that are not affected by degeneracies, can then be used as input for extracting $`\delta _{CP}`$ from the data collected with a shorter baseline as in the JHF–Kamioka experiment .


text:
# Search for Lepton Flavor Violation Process 𝐽/𝜓→𝑒𝜇

## 1 Introduction

In the minimal standard model of electroweak theory, the three separate lepton numbers, electron number, muon number, and tau number, are conserved, but this conservation law may be broken in many extensions of the Standard Model, such as grand unified models , supersymmetric standard models , left-right symmetric models  and models where electroweak symmetry is broken dynamically . Recent experimental results from Super-Kamiokande , SNO  and KamLAND  indicate strongly that neutrinos have masses and mix with each other; hence lepton flavor symmetries are broken. There have been many studies both experimentally and theoretically on searching for the lepton flavor violating processes , mainly from $`\mu `$, $`\tau `$ and Z decays . Theoretical studies on the possibility of searching for the lepton flavor violation in decays of charmonium and bottomonium systems are discussed in Refs. -. In this paper we report on a search for lepton flavor violation via the process J/$`\psi e^\pm \mu ^{}`$ using $`5.8\times 10^7J/\psi `$ events at BEPC/BESII .

## 2 BES detector

The BEijing Spectrometer (BES)  is a large general purpose solenoidal detector at the Beijing Electron Positron Collider (BEPC). The beam energy of BEPC ranges from 1.0 to 2.5 GeV, and the peak luminosity at the $`J/\psi `$ is around $`5\times 10^{30}`$ cm-2s-1. The upgrade from BESI to BESII  includes the replacement of the central drift chamber with a straw-tube vertex chamber, composed of 12 tracking layers arranged around a beryllium beam pipe and with a spatial resolution of about 90 $`\mu `$m; a new barrel time-of-flight (TOF) counter with a time resolution of 180 ps; and a new main drift chamber (MDC), which has 40 tracking layers with a $`dE/dx`$ resolution of $`\sigma _{dE/dx}=8.0\%`$ and a momentum resolution of $`\sigma _p/p=1.78\%\sqrt{1+p^2}`$ ($`p`$ in GeV/$`c`$) for charged tracks. These upgrades augment the pre-existing calorimeter and muon tracking systems. The barrel shower counter (BSC), which has a total thickness of 12 radiation lengths and covers 80 $`\%`$ of 4$`\pi `$ solid angle, has an energy resolution of $`\sigma _E/E=0.21/\sqrt{E}`$ ($`E`$ in GeV) and a spatial resolution of 7.9 mrad in $`\varphi `$ and 2.3 cm in z. The $`\mu `$ identification system consists of three double layers of proportional tubes interspersed in the iron flux return of the magnet. They provide coordinate measurements along the muon trajectories with resolutions in the outermost layer of 10 cm and 12 cm in $`r\varphi `$ and $`z`$, respectively. The absorber thicknesses in front of the three layers are 120, 140, and 140 mm, and the solid angle coverage of the layers is 67 %, 67 %, and 63 % of $`4\pi `$, respectively.

## 3 Event selection

The initial selection of $`J/\psi e^\pm \mu ^{}`$ requires that candidate events have two oppositely charged tracks within the polar angle region, $`|\mathrm{cos}\theta |<0.8`$, and fewer than four neutral tracks. A charged track should have a good helix fit; a momentum, $`P`$, with 1.45 GeV/$`c`$ $`<P<1.65`$ GeV/$`c`$; and originate from the interaction region, defined by $`|x|<0.015`$ m, $`|y|<0.015`$ m and $`|z|<0.15`$ m. To reject cosmic rays, the TOF time difference of the two charged tracks must be less than 1.2 ns. The invariant mass of the $`e\mu `$ system, $`M_{e\mu }`$, must satisfy 2.95 GeV/$`c^2<M_{e\mu }<`$ 3.25 GeV/$`c^2`$, and the angle between the two charged tracks, $`\theta _{12}`$, is required to be greater than 178.5. The $`\theta _{12}`$ distribution is shown in Fig. 1.

Isolated photons are defined as those photons having an energy deposit in the BSC greater than 50 MeV, an angle with any charged track greater than $`15^{}`$, and an angle between the direction defined by the first layer hit in the BSC and the developing direction of the cluster in the $`xy`$-plane less than $`18^{}`$. There must be no isolated photon in the selected event.

The above criteria select back-to-back, two prong events, such as $`J/\psi e^+e^{}(\gamma ),\mu ^+\mu ^{}(\gamma ),K^+K^{},\pi ^+\pi ^{}`$ and $`e^+e^{}e^+e^{}(\gamma )`$, etc. To select electrons and muons, the BSC and $`\mu `$ counter information is used. The actual selection criteria are based on distributions determined from data. Fig. 2 shows distributions of $`E/P`$, where $`E`$ is the energy deposited in the shower counter, for candidate electron and muon tracks. To be an electron, a track must have no hits in the muon counter and satisfy $`E/P>0.7`$.

To select muon tracks, the differences, $`\delta _i(i=x,y,z)`$, between the closest muon hit and the projected MDC track in each layer are used. Fig. 3 shows these differences in the third layer of the $`\mu `$ counter. The distributions are Gaussians, and standard deviations, $`\sigma _i(i=x,y,z)`$, are determined for each layer of the $`\mu `$ counter.

To reduce background, we make a tight cuts on the muon. A good hit in the $`\mu `$ counter requires $`|\delta _i|<\sigma _i`$ for $`i`$ = x, y, and z. The total number of layers hit in the $`\mu `$ counter is denoted as $`\mu _{hit}`$ and can range from zero to three. The number of good hits is denoted by $`\mu _{hit}^{good}`$. A track is considered as a muon if $`E/P<0.3`$ and $`\mu _{hit}^{good}`$ is greater than zero if the transverse momentum of the track, $`P_{xy}`$, is less than 0.75 GeV/$`c`$, $`\mu _{hit}^{good}>1`$ if 0.75 GeV/$`c`$ $`<P_{xy}<0.95`$ GeV/$`c`$, or $`\mu _{hit}^{good}=3`$ if $`P_{xy}>0.95`$ GeV/$`c`$.

Using these selection criteria, four $`J/\psi e\mu `$ candidates are found in $`5.8\times 10^7J/\psi `$ events. The characteristics of these events are listed in Table 1.

## 4 Efficiencies and backgrounds

The main backgrounds for $`J/\psi e\mu `$ are nearly back-to-back $`J/\psi e^+e^{}(\gamma )`$, $`\mu ^+\mu ^{}(\gamma )`$, $`K^+K^{}`$, $`\pi ^+\pi ^{}`$, $`e^+e^{}e^+e^{}(\gamma )`$, and $`e^+e^{}\mu ^+\mu ^{}(\gamma )`$ events where one or more tracks is misidentified as a $`\mu `$ or $`e`$. It is therefore important to determine the $`\mu `$ and $`e`$ particle identification (PID) misidentification efficiencies for the background channels, as well as the PID efficiencies for the signal channel. In this study, we determine these efficiencies using information from the BSC and muon counters from data. The overall efficiencies include the PID efficiencies and the geometric efficiencies, $`ϵ_{MC}`$, determined using our Monte Carlo (MC) program SIMBES (SIMulation at BES), which is based on GEANT3.21. For the signal channel $`J/\psi e\mu `$, the total selection efficiency is then

$$ϵ_{J/\psi e\mu }=ϵ_{e\mu MC}\times ϵ_{ee}\times ϵ_{\mu \mu },$$

where $`ϵ_{e\mu MC}`$ is the geometric efficiency, $`ϵ_{ee}`$ is the electron PID efficiency, and $`ϵ_{\mu \mu }`$ is the muon PID efficiency. Using 30000 $`J/\psi e\mu `$ Monte Carlo events with an angular distribution of $`1+\mathrm{cos}^2\theta `$, the geometric efficiency is determined to be $`ϵ_{e\mu MC}=(53.7\pm 0.3)\%`$.

For the background channels $`J/\psi `$ $`(e^+e^{})XX`$, where $`X=e`$, $`\mu `$, $`\pi `$, $`K`$, the efficiency after $`J/\psi e\mu `$ selection is

$$ϵ_{XX}=ϵ_{XXMC}\times ϵ_{Xe}\times ϵ_{X\mu }\times 2,$$

where $`ϵ_{XXMC}`$ is the Monte Carlo geometric acceptance, shown in Table 2, and $`ϵ_{Xe}`$ and $`ϵ_{X\mu }`$ are the particle misidentification efficiencies for X being identified as an electron or a muon, respectively.

The $`e`$ and $`\mu `$ PID efficiencies are determined using $`e`$-pair samples of $`e^+e^{}e^+e^{}(\gamma )`$ (Bhabha, resonance, and continuum production) events and $`\mu `$-pair samples of $`e^+e^{}\mu ^+\mu ^{}(\gamma )`$ (resonance and continuum production) events. The highest momentum track is required to satisfy $`1.45<P<1.65`$ GeV/c and to be identified as an $`e`$ or $`\mu `$ according to the PID selection criteria. The other track in the event is then assumed to be the same type of lepton, and the fraction of these satisfying the lepton selection determines the efficiency. The PID efficiencies are shown in Table 3, and the total efficiency for $`J/\psi e\mu `$ is $`(9.7\pm 0.5)\%`$.

The misidentification efficiencies between electrons and muons cannot be studied from the data because the events selected for this study are the same as our $`J/\psi e\mu `$ candidates. A rough estimate of this background using a Monte Carlo simulation predicts a total of 7 background events from $`e^+e^{}e^+e^{}(\gamma )`$ and $`e^+e^{}\mu ^+\mu ^{}(\gamma )`$. Since the Monte Carlo simulation is not precise enough to determine this background, we conservatively ignore these backgrounds in the estimation of the upper limit on the branching ratio for $`J/\psi e\mu `$.

For the misidentification efficiencies of hadrons, $`\pi `$ and $`K`$ tracks are selected from $`J/\psi \rho ^\pm \pi ^{}`$ and $`K^\pm K^{}`$ decays. The track with the highest momentum is considered as the $`\pi `$ and $`K`$ in the respective samples, and a total of 60551 $`\pi `$ tracks and 9275 $`K`$ tracks are selected. The fraction of the tracks above 1.2 GeV/c passing the electron and muon selection criteria are used to obtain the misidentification efficiencies, that are listed in Table 3. The resultant background rates from particle misidentification of $`J/\psi K^+K^{},\pi ^+\pi ^{}`$ events are listed in Table 4, taking into account their branching ratios from PDG2002 .

There are a number of sources of systematic error, including the error on the total number of $`J/\psi `$ ($`4.7\%`$) events, the errors on the branching ratios from PDG2002 , the uncertainties of the muon and electron identification efficiencies, shown in Table 3, and the errors on the hadronic misidentification efficiencies, also shown in Table 3. Table 5 summarizes the errors in the background from the hadronic channels. The total background, ignoring the contributions from $`J/\psi e^+e^{}(\gamma )`$ and $`\mu ^+\mu ^{}(\gamma )`$, is $`2.3\pm 0.3`$ events. Therefore the observed four events are consistent with background.

## 5 Results and discussion

The four observed $`J/\psi e\mu `$ candidates are consistent with the estimated background, hence an upper limit for $`J/\psi e\mu `$ is determined. Only the backgrounds from the hadronic channels are taken into account in the estimation of the upper limit.

There are several ways to determine the upper limit   . The following formula is used to calculate the upper limit for $`Br(J/\psi e\mu )`$:

$$Br<\lambda (N_{OB},N_{BG})/[N_T\times ϵ_{J/\psi e\mu }],$$

where $`\lambda `$ is the upper bound of a $`90\%`$ C.L. for an unknown Poisson signal mean, for total observed events, $`N_{OB}`$, and known mean background, $`N_{BG}`$. $`N_T`$ is the total number of $`J/\psi `$ events, and $`ϵ`$ is the detection efficiency. $`\lambda (N_{OB},N_{BG})`$ can be calculated using the method of Ref. .

As a conservative estimation, we take the number of background events as the central value of the number of background events reduced by one standard deviation, $`i.e.`$ $`N_{BG}=2.0`$ and $`N_{OB}`$=4, then $`\lambda (N_{OB},N_{BG})`$=5.98. Therefore we obtain:

$$Br(J/\psi e\mu )<1.1\times 10^6.$$

In summary, we searched for $`J/\psi e\mu `$ based on $`5.8\times 10^7J/\psi `$ events and observed four $`J/\psi e\mu `$ candidates passing our selection criteria, which are consistent with the estimated background. The upper limit of $`Br(J/\psi e\mu )`$ is determined to be $`1.1\times 10^6`$ at the 90 $`\%`$ C.L.

## 6 Acknowledgments

The BES collaboration thanks the staff of BEPC for their hard efforts. We also thank Profs. Xinmin Zhang and Jianxiong Wang for helpful discussions. This work is supported in part by the National Natural Science Foundation of China under contracts Nos. 19991480, 10175060 and the Chinese Academy of Sciences under contract No. KJ 95T-03, the 100 Talents Program of CAS under Contract Nos. U-24, U-25, and the Knowledge Innovation Project of CAS under Contract Nos. U-602, U-34 (IHEP); by the National Natural Science Foundation of China under Contract No.10175060(USTC); and by the Department of Energy under Contract Nos. DE-FG03-93ER40788 (Colorado State University), DE-AC03-76SF00515 (SLAC), DE-FG03-94ER40833 (U Hawaii), DE-FG03-95ER40925 (UT Dallas).