format: text

text:
# ISOCAM observations of the 𝜌 Ophiuchi cloud: Luminosity and mass functions of the pre-main sequence embedded cluster This work is based on observations with ISO, an ESA project with instruments funded by ESA Member States (especially the PI countries: France, Germany, the Netherlands, and the United Kingdom) with the participation of ISAS and NASA.,Table 1 is only available in electronic form at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsweb.u-strasbg.fr/

## 1 Introduction

Recent observations suggest that most stars in our Galaxy and other galaxies form in compact clusters. In particular, near-IR imaging surveys of nearby molecular cloud complexes have shown that the star formation activity is typically concentrated within a few rich clusters associated with massive dense cores which constitute only a small fraction of the total gas mass available (e.g. Lada, Strom, & Myers lada93 (1993), Zinnecker, McCaughrean, & Wilking zinnecker (1993)). These embedded clusters comprise various types of young stellar objects (YSOs) – from still collapsing protostars to young main sequence stars – but are usually dominated in number by low-mass pre-main sequence (PMS) stars, i.e., T Tauri stars. Young clusters provide excellent laboratories for studying the formation and early evolution of stars through the observational analysis of large, genetically homogeneous samples of embedded YSOs. Two key characteristics of these young stellar populations are their luminosity distribution and their mass spectrum, which give important observational constraints on the stellar initial mass function (IMF) in the Galaxy. Observations of young embedded clusters can also help us understand possible links between parent cloud properties and the resulting stellar masses. They are however hampered by: (1) dust extinction from the parent cloud which hides most of the newly formed stars at optical wavelengths, (2) the difficulty to recognize the nature of individual sources (e.g. protostars, T Tauri stars, or background sources), (3) the youth (and thus poorly known intrinsic properties) of most cluster members.

Owing to these difficulties, the census of embedded YSOs provided by IRAS and near-IR studies is far from complete even in the nearest clouds (e.g. Wilking, Lada, & Young wly (1989) – hereafter WLY89). Thanks to its high sensitivity and good spatial resolution in the mid-IR, the ISOCAM camera on board ISO (Cesarsky et al. cesarsky (1996); Kessler et al. kessler (1996)) was a powerful tool to achieve more complete surveys for YSOs in all major nearby star-forming regions (see Nordh et al. nordh1 (1996), nordh2 (1998); Olofsson et al. ohk99 (1999); Persi et al. pmo00 (2000)).

The nearby $`\rho `$ Ophiuchi cloud is one of the most actively studied sites of low-mass star formation. Its central region harbors a rich embedded cluster with about 100 members recognized prior to the present work (e.g. WLY89, Casanova et al. cmfa (1995)). While a dispersed population of optically visible young stars, associated with the Upper-Scorpius OB association, has a typical age of several million years (Myr) (e.g. Preibisch & Zinnecker pz99 (1999)), the central embedded cluster is recognized as one of the youngest clusters known with an estimated age on the order of $``$ 0.3–1 Myr (e.g. WLY89, Greene & Meyer gm95 (1995), Luhman & Rieke lr99 (1999)). This young cluster has been extensively studied at wavelengths ranging from the X-ray to the radio band. The satellites Einstein and $`ROSAT`$ have revealed $`70`$ highly variable X-ray sources associated with magnetically-active young stars, including deeply embedded protostellar sources (Montmerle et al. montmerle83 (1983), Casanova et al. cmfa (1995), Grosso et al. grosso (2000)). In the near-IR, the cloud has been deeply surveyed from the ground using large-format arrays (Greene & Young gy (1992), Comerón et al. crbr (1993), Strom et al. sks (1995), Barsony et al. bklt (1997)). Unfortunately, due to difficulties in discriminating between background sources and embedded YSOs without performing time-consuming mid-IR photometry (e.g. Greene et al. gwayl (1994)) or near-IR spectroscopy (e.g. Greene & Lada gl96 (1996); Luhman & Rieke lr99 (1999)), these recent near-IR surveys have only partially increased the number of classified, recognized members. Finally, while only relatively poor angular resolution $`IRAS`$ data are available so far in the far-IR (e.g. WLY89), deep imaging surveys at an angular resolution of $`1015\mathrm{}`$ or better exist at (sub)millimeter wavelengths (e.g., Motte, André, & Neri man (1998) – hereafter MAN98 –; Wilson et al. wilson99 (1999)).

It is thanks to the illuminating example of the $`\rho `$ Ophiuchi embedded cluster that the now widely used empirical classification of YSOs was originally introduced. Three IR classes were initially distinguished based on the shapes of the observed spectral energy distributions (SEDs) between $``$ 2 $`\mu `$m and $``$ 25–100 $`\mu `$m (Lada & Wilking lw84 (1984), WLY89). Objects with rising SEDs in this wavelength range were classified as Class I, sources with SEDs broader than blackbodies but decreasing longward of $``$ 2 $`\mu `$m as Class II, and sources with SEDs consistent with (or only slightly broader than) reddened stellar blackbodies as Class III. These morphological SED classes are interpreted in terms of an evolutionary sequence from (evolved) protostars (Class I), to T Tauri stars with optically thick IR circumstellar disks (Class II), to weak T Tauri stars with at most optically thin disks (Class III) (Lada lada87 (1987), Adams, Lada, & Shu als (1987), André & Montmerle am (1994) – hereafter AM94). A fourth class (Class 0) was subsequently introduced by André, Ward-Thompson, & Barsony (awb (1993)) to accommodate the discovery in the radio range of cold sources with large submillimeter to bolometric luminosity ratios and powerful jet-like outflows, such as VLA 1623 in $`\rho `$ Oph A (e.g. André et al. andre90 (1990); Bontemps et al. batc (1996)). Class 0 objects, which have measured circumstellar envelope masses larger than their inferred central stellar masses, are interpreted as young protostars at the beginning of the main accretion phase (e.g. André, Ward-Thompson, & Barsony awb00 (2000)). The fact that the $`\rho `$ Ophiuchi central region contains at least two Class 0 protostars (André, Ward-Thompson, & Barsony awb (1993))), as well as numerous ($`\stackrel{>}{}60`$) pre-stellar condensations (MAN98), demonstrates that it is still actively forming stars at the present time.

The distance to the $`\rho `$ Ophiuchi cloud is somewhat uncertain. Usually, a value of 160 pc is adopted (e.g. Chini chini81 (1981)). However, recent Hipparcos results on the Upper-Scorpius OB association (de Zeeuw et al. dezeeuw (1999)) provide a reasonably accurate estimate of $`145\pm 2`$ pc for the average distance to the stars of the OB association. The embedded cluster is located at the inner edge of the molecular complex on the outskirts of the OB association (e.g. de Geus degeus (1989)), and not very far, in projection, from the association center (less than $`4\mathrm{°}`$ apart which corresponds to $`10`$ pc; see Fig. 9 of de Zeeuw et al. dezeeuw (1999)). In this paper, we therefore adopt a distance of $`d=140\pm 10`$ pc for the $`\rho `$ Ophiuchi IR cluster which corresponds to a distance modulus $`5\times `$log$`{}_{10}{}^{}(d/10\mathrm{p}\mathrm{c})=5.73\pm 0.15`$.

The layout of the paper is as follows. Sect. 2 gives observational details (§ 2.2) and describes the way the data have been reduced to obtain mid-IR images and extract point-sources (§ 2.3) along with the photometric uncertainties and the sensitivity levels (§ 2.4). In Sect. 3, the identification of detected sources is discussed (§ 3.1) and the selection of a new population of 123 Class II YSOs, as well as 16 Class I and 77 Class III YSOs, is described (§ 3.2-3.5). We then derive luminosities for these YSOs and build the corresponding luminosity functions in Sect. 4. In Sect. 5, we model the luminosity functions for Class II and Class III YSOs in terms of the underlying mass function and star formation history. In Sect. 6, we discuss the resulting constraints on the IMF of the $`\rho `$ Ophiuchi cloud down to $`0.06\text{M}\text{}`$, (§ 6.1-6.2) as well as related implications (§ 6.3-6.5).

## 2 Observations and data reduction

### 2.1 Region surveyed by ISOCAM

The survey encompasses the $`\rho `$ Ophiuchi central region associated with the prominent dark cloud L1688, as well as the two subsidiary sites L1689N and L1689S. The L1688 field is a $`45\mathrm{}\times 45\mathrm{}`$ square, while the L1689N and L1689S fields each cover an area of $`16.5\mathrm{}\times 16\mathrm{}`$ (see Fig. 1). Most previously known members of the $`\rho `$ Ophiuchi cluster lie within these fields. In particular, this is the case for 94 of a total of 113 recognized members from WLY89, AM94, Greene et al. (gwayl (1994)), and MAN98. The known young stars which lie outside the boundaries of our survey are mostly optically visible, weak-line or post T Tauri stars (belonging to Class III) spread over a large area on the outskirts of the molecular complex (e.g. Martín et al. martin98 (1998)).

### 2.2 Observational details

The mapping was performed in the raster mode of ISOCAM in which the mid-IR $`32\times 32`$ pixel array imaged the sky at consecutive positions along a series of scans parallel to the right-ascension axis. The offset between consecutive array positions along each scan ($`\mathrm{\Delta }\alpha `$) was 15 pixels, while the offset between two scans ($`\mathrm{\Delta }\delta `$) was 26 pixels. Each set of scans was then co-added and combined into a single raster image. The final image of the L1688 field (Fig. 1 and Abergel et al. 1996) actually results from the combination of six separate rasters. A pixel field of view of 6$`\mathrm{}`$ was used for four of these rasters, but smaller 3$`\mathrm{}`$ pixels were employed for the other two rasters in order to avoid saturating the array on the brightest sources of the cluster. The L1689N and L1689S fields were imaged with one raster each using 6$`\mathrm{}`$ pixels.

In order to avoid saturation, the individual readout time for the L1688 rasters was set to $`t_{\mathrm{int}}=0.28`$s. About 55 of these readouts (i.e. an integration time of $`15`$s) were performed per sky position. Thanks to the half-frame overlap between subsequent individual images, each sky position was observed twice, yielding an effective total integration time of $`30`$s. For L1689N and L1689S it was possible to use $`t_{\mathrm{int}}=2.1`$s, and about 15 readouts were performed per sky position with the same half-frame overlap, giving an integration time per sky position of $`60`$s. A total of 1104 individual images were necessary to mosaic the L1688 field, and an additional 60 images each were used to map L1689N and L1689S.

All three fields were mapped in two broad–band filters of ISOCAM: LW2 (5–8.5 $`\mu `$m) and LW3 (12–18 $`\mu `$m). These filters are approximately centered on two minima of the interstellar extinction curve and are situated apart from the silicate absorption bands (at roughly 10 and 18 $`\mu `$m). However, they include most of the Unidentified Infrared Bands (UIBs, likely due to PAH-like molecules) which constitute a major source of background emission toward star-forming clouds (e.g. Bernard et al. bernard (1993), Boulanger et al. boulanger (1996)). The ISOCAM central wavelengths adopted here for LW2 and LW3 are 6.7 $`\mu `$m and 14.3 $`\mu `$m respectively.

### 2.3 Image processing, source extraction, and photometry

Each raster consists of a temporal series of individual integration frames (i.e. of $`32\times 32`$ pixel images) which was reduced using the CAM Interactive Analysis software (CIA)111 CIA is a joint development by the ESA Astrophysics Division and the ISOCAM consortium led by the ISOCAM PI, C. Cesarsky, Direction des Sciences de la Matière, C.E.A., France.. We have subtracted the best dark current from the ISOCAM calibration library, and as a second step we improved it with a second order correction using a FFT thresholding method (Starck et al. starck2 (1999)). Cosmic-ray hits were detected and masked using the multi-resolution median transform algorithm (Starck et al. starck1 (1996)). The transients in the time history of each pixel due to detector memory effects were corrected with the inversion method described in Abergel et al. (abergel (1996)). The images were then flat-fielded with a flat image obtained from the observations themselves. Since these various corrections applied to the images are not perfect, the extraction of faint sources from the images is a difficult task. We have developed an interactive IDL point-source detection and photometry program for raster observations which works in the CIA environment. This program helps to discriminate between astronomical sources and remaining low-level glitches or ghosts due to strong transients (see also Nordh et al. nordh1 (1996); Kaas et al. kaas01 (2001)). The fluxes of the detected sources were estimated from the series of flux measurements made in the individual images (usually 2 to 4 individual images cover each source) which were obtained from classical aperture photometry. The emission was integrated in a sky aperture, the background emission subtracted, and finally an appropriate aperture correction was applied based on observed point-spread functions available in the ISOCAM calibration library. In practice, the radius of the aperture used was 9$`\mathrm{}`$ (i.e., 3 and 1.5 pixels for a pixel size of $`3\mathrm{}`$ and 6$`\mathrm{}`$, respectively). For the weakest sources, however, we reduced the aperture radius to 4.5$`\mathrm{}`$ (i.e. 1.5 pixels for a pixel size of $`3\mathrm{}`$), in order to improve the signal-to-noise ratio. Finally, we applied the following conversion factors: $`2.33`$ and $`1.97`$ ADU/gain/s/mJy for LW2 and LW3 respectively (from in-orbit latest calibration - Blommaert blommaert (1998)). These calibration factors are strictly valid only for sources with a flat SED ($`F_\nu \nu ^1`$). Here, a small but significant ($`\stackrel{>}{}1`$ %) color correction needs to be applied to the bluest sources, recognized as Class III YSOs in Sect. 3 below. For these sources, the conversion factors quoted above were divided by 1.05 for LW2 and 1.02 for LW3 to account for the color effect. The 212 ISOCAM sources recognized as cluster members (see Sect. 3 below) are listed in Table 1 (available only in electronic form at http://cdsweb.u-strasbg.fr/) with their J2000 coordinates, their flux densities and associated rms uncertainties (see Sect. 2.4), as well as the corresponding near-IR identifications.

### 2.4 Photometric uncertainties and point-source sensitivity

The uncertainties on the final photometric measurements result from systematic errors due to uncertainties in the absolute calibration and the aperture correction factors, and from random errors associated with the flat-fielding noise, the statistical noise in the raw data, the noise due to remaining low-level glitches, and the imperfect correction for the transient behavior of the detectors. The in-orbit absolute calibration has been verified to be correct to within 5 % (Blommaert blommaert (1998)), and we estimate that the maximum systematic error on the aperture correction is $`10`$% (by comparing theoretical and observed point-spread functions). The maximum systematic error on our photometry is thus $`15`$%. The magnitudes of the random errors were directly estimated from the data by measuring both a “temporal” noise (noise in the temporal sequence of individual integrations) and a “spatial” noise (due to imperfect flat-fielding and/or spatial structures in the local mid-IR background emission) for each source in the automatic detection procedure. The temporal noise was computed as the standard deviation of the individual aperture measurements divided by the square root of the number of measurements. The spatial noise was estimated as the standard deviation around the mean background (linear combination of the median and the mean of the pixels optimized for the source flux estimates) in the immediate vicinity of each source.

The sensitivity limit of the survey was estimated by calculating the average value of the quadratic sum of the temporal and spatial noises measured on the weakest detected sources. The total rms flux uncertainty found in this way, $`\sigma _{tot}=(\sigma _{temp}^2+\sigma _{spat}^2)^{1/2}`$, is $`\sigma _{6.7}=2.2`$ mJy at 6.7 $`\mu `$m and $`\sigma _{14.3}=4.1`$ mJy at 14.3 $`\mu `$m, $``$ 75 % of which is due to the spatial noise component. The large contribution of the spatial noise originates in the highly structured diffuse mid-IR emission from the ambient molecular cloud itself (see Abergel et al. abergel (1996)). Figure 2 displays the distributions of fluxes at 6.7 $`\mu `$m and 14.3 $`\mu `$m for all the detected ISOCAM sources. We used the Wainscoat et al. (wainscoat (1992)) Galactic model of the mid-IR point source sky to estimate the expected number of foreground and background sources up to a distance of 20 kpc. The model predictions are shown by solid and dashed curves in Fig. 2 for cloud extinctions of $`A_V=0`$ and $`A_V=10`$, respectively (see Kaas et al. kaas01 (2001) for more details). It can be seen that the flux histograms of the ISOCAM sources not associated with YSOs (light shading in Fig. 2) are remarkably similar in shape to the model distributions down to $`6`$ mJy at 6.7 $`\mu `$m and $`10`$ mJy at 14.3 $`\mu `$m. These flux densities can be used to estimate the completeness level of our observations which is not uniform over the spatial extent of the survey. The histograms with light (grey) shading in Fig. 2 are dominated by background sources preferentially located in low-noise regions (i.e., outside the crowded central part of L1688), where the total rms flux uncertainty is $`2.0`$ mJy at 6.7 $`\mu `$m and $`3.5`$ mJy at 14.3 $`\mu `$m. The effective completeness level in these regions is thus $`3\sigma _{tot}`$, where $`\sigma _{tot}`$ is the total222 The detection level is more directly related to the temporal noise than to the spatial noise. In terms of $`\sigma _{temp}`$, the effective completeness level is $`12\sigma _{temp}`$. flux uncertainty (see above). However, most of the YSOs are located in regions where the noise is somewhat larger. The largest rms noise is reached in the Oph A core area (see Fig. 1), where $`\sigma _{tot}3.4`$mJy at 6.7 $`\mu `$m and $`\sigma _{tot}5.0`$mJy at 14.3 $`\mu `$m. Therefore, we conservatively estimate the completeness levels of the global ISOCAM survey to be $`10`$ mJy at 6.7 $`\mu `$m and $`15`$ mJy at 14.3 $`\mu `$m.

Finally, we note that the $`A_V=10`$ model curve in Fig. 2b accounts for essentially all the “blue” sources detected at 14.3 $`\mu `$m and not associated with known YSOs. At 6.7 $`\mu `$m, the predictions of the Wainscoat et al. model suggest that there might still be a slight excess of $`30`$ unidentified sources belonging to the cloud (Fig. 2a).

## 3 Identification and nature of the mid-IR sources

### 3.1 Statistics of detections

Within the $`0.7`$ square degree imaged by ISOCAM, a total of 425 sources has been identified, among which 211 are seen at both 6.7 $`\mu `$m and 14.3 $`\mu `$m. The spatial distribution of the sources is shown in Fig. 1, where the “red” sources \[those with log10($`F_\nu ^{14.3}/F_\nu ^{6.7})>0.2`$, see below\] are indicated as filled circles. These “red” sources appear to be clustered into four main groupings: three sub-clusters in L1688, i.e., Oph A (West), Oph B (North-East), Oph EF (South) (see also Strom, Kepner, & Strom sks (1995)), as well as a new sub-cluster in L1689S.

Four bright embedded stars (S1, SR3, WL16, WL22) are spatially resolved by ISOCAM in both filters. Their extended mid-IR emission is most likely due to PAH-like molecules excited by relatively strong far-ultraviolet (FUV) radiation fields333 S1 and SR3 are optically-visible stars of spectral types B3 and B7, respectively, which clearly emit enough UV photons to excite PAHs. The other two objects, WL16 and WL22, show PAH features in their mid-IR spectra (e.g. Moore et al. moore (1998)), and may be young embedded early-type (i.e., early A or late B) stars.. These bright sources are displayed as open star symbols in Fig. 1 (SR3, S1, WL22, WL16 from right to left).

A total of 89 previously classified YSOs lie within the area of the present survey: 2 Class 0, 69 Class I/II, and 18 Class III YSOs. ISOCAM detected 84 of these 89 YSOs (94 %): 97 % of the Class I/IIs, 94 % of the Class IIIs, and none of the Class 0s. The two undetected Class I/IIs (GY256, GY257), and the undetected Class III (IRS50), are located very close ($``$ 10″–15″) to bright sources (WL6 and IRS48, respectively), which may account for their non-detection. The two Class 0 objects (VLA 1623 and IRAS 16293$``$2422) are deeply embedded within massive, cold circumstellar envelopes which are probably opaque at 6.7 $`\mu `$m and 14.3 $`\mu `$m (see André, Ward-Thompson, & Barsony awb (1993))) and too weak to be detected by ISOCAM.

### 3.2 A population of new YSOs with mid-IR excess

As can be seen in Fig. 3, the mid-IR regime is ideal to detect and characterize the excess emission due to circumstellar disks around young stars. In Fig. 3a, the 207 sources detected at 6.7 $`\mu `$m and 14.3 $`\mu `$m (excluding S1, SR3, WL16, WL22) are shown in a diagram which displays the logarithmic flux ratio log10($`F_\nu ^{14.3}/F_\nu ^{6.7}`$) against the mid-IR flux $`F_\nu ^{14.3}`$. On the right axis of the diagram, the flux ratio $`F_\nu ^{14.3}/F_\nu ^{6.7}`$ has been converted into a classical IR spectral index, $`\alpha _{\mathrm{IR}}=`$ d$`\mathrm{log}_{10}(\lambda F_\lambda )/`$d$`\mathrm{log}_{10}(\lambda )`$ (e.g. WLY89), calculated between 6.7 $`\mu `$m and 14.3 $`\mu `$m, i.e., $`\alpha _{\mathrm{IR}}^{714}`$. In this diagram, two groups of sources can clearly be distinguished. The lower ratio group has log10($`F_\nu ^{14.3}/F_\nu ^{6.7})0.7`$, i.e., a spectral index on the order of $`\alpha _{\mathrm{IR}}^{714}3.0`$, which is the value expected for simple photospheric blackbody emission in the Rayleigh-Jeans regime. The dispersion around $`\alpha _{\mathrm{IR}}^{714}3.0`$ is obviously larger for weaker sources. This is mainly due to increasing photometric uncertainty with decreasing flux. The higher ratio group consists of “red” sources defined by log10($`F_\nu ^{14.3}/F_\nu ^{6.7}`$$`>0.2`$ (cf. Nordh et al. nordh1 (1996)) . Most of these actually have $`0.2`$ log10($`F_\nu ^{14.3}/F_\nu ^{6.7})0.3`$, i.e., $`1.6\stackrel{<}{}\alpha _{\mathrm{IR}}^{714}\stackrel{<}{}0`$, which is typical of classical T Tauri stars (hereafter CTTS). This range of $`\alpha _{\mathrm{IR}}^{714}`$ roughly delineates the domain of Class II YSOs (e.g. Adams et al. 1987, AM94, Greene et al. 1994), and is usually interpreted in terms of optically thick circumstellar disk emission (e.g. Lada & Adams la92 (1992)). In Fig. 3, the new ISOCAM sources are distinguished from the previously known cluster members by different symbols. One can see in Fig. 3a that all the previously known Class II sources but one lie above the dividing line for red sources, while all the previously known Class III sources but two lie below it.

The object WL19, classified as a Class I YSO by WLY89 and as a reddened Class II by AM94 (see also Lada & Wilking lw84 (1984)), is here found to be a “blue” source in the mid-IR range ($`\alpha _{\mathrm{IR}}^{714}=1.8\pm 0.2`$). This object may correspond to a luminous Class III star located behind the cloud (see also Comerón et al. crbr (1993)). Although GY12 formally has a “red” mid-IR spectral index ($`\alpha _{\mathrm{IR}}^{714}=0.5\pm 0.5`$) here, we still consider it as a Class III object (cf. Greene et al. gwayl (1994)). (The mid-IR color is highly uncertain since GY12 is only marginaly resolved at 14.3 $`\mu `$m from its bright Class I neighbor GSS30.) Finally the Class III source DoAr21 has a borderline mid-IR spectral index ($`\alpha _{\mathrm{IR}}^{714}=1.5\pm 0.1`$) but is kept as a Class III object since its color between 2.2 $`\mu `$m and 14.3 $`\mu `$m corresponds to $`\alpha _{\mathrm{IR}}^{214}=2.0\pm 0.1`$.

A total of 71 sources are identified for the first time as mid-IR excess objects in Fig. 3a. These new “red” sources are most likely all embedded YSOs, i.e., members of the $`\rho `$ Ophiuchi cluster. Since dust extinction is roughly the same at 6.7 $`\mu `$m and 14.3 $`\mu `$m (Rieke & lebofsky rl (1985); Lutz lutz (1999)), any background source should be intrinsically red in order to contaminate the sample of red YSOs. Based on the Galactic model by Wainscoat et al. (wainscoat (1992)), the vast majority of background objects should appear “blue” (cf. Fig. 2). Seven background giant stars (GY45, GY65, VSSG6, GY232, GY351, GY411, and GY453 – cf. Luhman & Rieke lr99 (1999)), and a known foreground dwarf (HD148352 – Garrison garrison (1967)) are detected, which all have blue mid-IR colors ($`\alpha _{\mathrm{IR}}^{714}=2.7,2.0,3.1,3.3,2.9,3.2,3.1`$, and $`3.1`$, respectively).

The newly identified cloud members nicely extend the previously known Class I/II population toward low IR fluxes. While previous studies could only identify sources with $`F_\nu ^{14.3}\stackrel{>}{}100`$ mJy, the present census is complete for objects down to $`F_\nu ^{14.3}15`$ mJy. Altogether, a sample of 139 Class I/Class II YSOs is identified, of which 71 are new members. The present survey has thus allowed us to more than double the number of recognized YSOs with circumstellar IR excess in the $`\rho `$ Ophiuchi cloud.

The 139 red YSOs are listed in Table 2 (Class I YSOs) and Table 3 (Class II YSOs) by decreasing order of $`\alpha _{\mathrm{IR}}^{214}`$. Twelve of them (last entries of Table 3) are completely new sources with respect to published IR surveys.

### 3.3 Class I versus Class II YSOs

In Fig. 3, the YSOs classified as Class I and Class II by AM94 and Greene et al. (gwayl (1994)) are shown as filled stars and filled circles, respectively. In the log10($`F_\nu ^{14.3}`$/$`F_\nu ^{6.7}`$) versus $`F_\nu ^{14.3}`$ diagram, there is no clear color gap between the two classes of objects, even though the Class I YSOs tend to lie in the upper part of the “red” group (see Fig. 3a). Since extinction has a negligible effect on the mid-IR ratio $`F_\nu ^{14.3}`$/$`F_\nu ^{6.7}`$, this suggests that Class I and Class II objects have fairly similar intrinsic colors between 6.7 $`\mu `$m and 14.3 $`\mu `$m.

The classical IR spectral index calculated from 2 $`\mu `$m to 10 $`\mu `$m (or 25 $`\mu `$m) (e.g. Lada & Wilking lw84 (1984) and WLY89) appears to provide a better way of discriminating between envelope-dominated Class I YSOs and disk-dominated Class II sources (see Fig. 3b). In particular, millimeter continuum mapping of optically thin circumstellar dust emission confirms that, apart from a few important exceptions (e.g., WL22, WL16, WL17, IRS37, IRS47), the $`\rho `$ Ophiuchi objects selected on the basis of $`\alpha _{\mathrm{IR}}^{210}\stackrel{>}{}0`$ are indeed Class I protostars surrounded by spheroidal envelopes (AM94, MAN98).
As expected, the previously known Class I YSOs are well concentrated in the upper-right part of the $`\alpha _{\mathrm{IR}}^{27}`$ versus $`\alpha _{\mathrm{IR}}^{214}`$ diagram of Fig. 3b (where $`\alpha _{\mathrm{IR}}^{27}`$ and $`\alpha _{\mathrm{IR}}^{214}`$ are close to the $`\alpha _{\mathrm{IR}}^{210}`$ index used in previous studies – e.g., AM94 and Greene et al. gwayl (1994)). There is also a hint of two gaps in this diagram at $`\alpha _{\mathrm{IR}}^{214}0.5`$ and $`\alpha _{\mathrm{IR}}^{214}0.0`$, which roughly bracket the regime of flat-spectrum sources as defined by Greene et al. (gwayl (1994)). These may represent a distinct population of transition objects between Class I and Class II (e.g. Calvet et al. calvet94 (1994)).
Here, we thus consider sources with $`\alpha _{\mathrm{IR}}^{214}>0.55`$ as Class I YSOs (Table 2), sources with $`0.05<\alpha _{\mathrm{IR}}^{214}<0.55`$ as candidate flat-spectrum objects (see Table 3), and sources with $`\alpha _{\mathrm{IR}}^{214}<0.55`$ as Class II YSOs (Table 3). These limiting indices are displayed in Fig. 3b. In the following, the candidate flat-spectrum sources will be treated as Class II YSOs.

### 3.4 Mid-IR excess versus near-IR excess

Most of the ISOCAM sources (e.g. 90 % of the Class II sources) were also detected in the near-IR (JHK) survey of Barsony et al. (bklt (1997)). Comparison between Fig. 3b and Fig. 3c illustrates the advantage of mid-IR measurements for selecting sources with intrinsic circumstellar IR excesses. While the red and blue groups of Fig. 3a are well separated in the mid-IR diagram of Fig. 3b, they blend together in the near-IR diagram of Fig. 3c.

Figure 3c also shows that most of the ISOCAM-selected YSOs lie within the reddening band associated with the intrinsic locus of CTTSs as derived by Meyer et al. (meyer97 (1997)) in Taurus. The few exceptions, which lie to the right of the reddening band, correspond to Class I and flat-spectrum YSOs.

### 3.5 The Class III YSO census in $`\rho `$ Ophiuchi

While Class I/II YSOs are easily recognized in the mid-IR range thanks to their strong IR excesses, Class III objects are difficult to identify without deep X-ray and/or radio centimeter continuum observations. We have used the ROSAT X-ray surveys of Casanova et al. (cmfa (1995)) and Grosso et al. (grosso (2000)), along with the VLA radio surveys by, e.g., André, Montmerle, & Feigelson (andre87 (1987)) and Stine et al. (stine88 (1988)) to build up a sample of bona-fide Class III YSOs covered by the present survey. With the additional Class III candidate WL19 (see Sect. 3.2 above), there are 38 such YSOs which are listed in order of decreasing $`\alpha _{\mathrm{IR}}^{214}`$ in Table 4.

This Class III sample is unfortunately not as complete as the Class I and Class II samples discussed above. According to Grosso et al. (grosso (2000)), the number of Class IIIs may be roughly as large as the number of Class IIs: above their typical X-ray detection limit of $`L_X3\times 10^{29}`$ erg s-1 (corresponding to $`L_{}0.3L_{}`$ – see Fig. 7 of Grosso et al.), they found a Class III/Class II number ratio of 19/22 in the $`ROSAT`$-$`HRI`$/$`ISO`$-$`ISOCAM`$ overlapping survey area. If this ratio is representative of the complete population of young stars in $`\rho `$ Ophiuchi, the total number of Class IIs found here (123 objects) suggests that as many as $`106`$ Class IIIs may be present in the cluster down to $`L_{}0.03\text{L}\text{}`$ (our completeness level for Class IIs, see Sect. 4.4). A total of 38 Class IIIs are already known within the ISOCAM survey area, so that $`68`$ unknown Class IIIs may remain to be found. Since it was noted in Sect. 2.4 that $`30`$ sources detected at 6.7 $`\mu `$m might be unidentified cluster members, about half of the missing Class IIIs may have been actually seen by ISOCAM.

We also note that a large proportion (80/123) of the Class II sources are closely associated with the densest part of L1688 (see Fig. 1). Assuming the same proportion applies to Class IIIs, we would expect $`44`$ unknown Class III sources to be located within the CS contours of Fig. 1. A total of 39 unclassified ISOCAM sources (also detected by Barsony et al. bklt (1997)) lie within these CS contours where the number of detected background stars should be small due to high cloud extinction. Most of these 39 sources might thus be yet unidentified Class III YSOs. These candidate Class III sources are listed in order of decreasing $`\alpha _{\mathrm{IR}}^{214}`$ in Table 5.

## 4 Luminosity estimates and luminosity function

Most of the new YSOs identified by ISOCAM are weak IR sources which were not detected by IRAS and were not observed in previous ground-based mid-IR surveys (dedicated to bright near-IR sources). They likely correspond to low-luminosity, low-mass young stars. In Sect. 4.1 below, we derive stellar luminosity estimates for Class II and Class III objects using published near-IR photometry from Barsony et al. (1997). In Sect. 4.2, we provide mid-IR estimates of the disk luminosities, $`L_{\mathrm{disk}}`$, for Class II YSOs. Finally, calorimetric estimates of the bolometric luminosities, $`L_{\mathrm{bol}}`$, for Class I YSOs are calculated in Sect. 4.3. The luminosity function of the $`\rho `$ Ophiuchi embedded cluster is then assembled in Sect. 4.4.

### 4.1 Dereddened J-band stellar luminosities for Class II and Class III sources

The J-band flux provides a good tracer of the stellar luminosity for late-type PMS stars (i.e., T Tauri stars) because the J-band is close to the maximum of the photospheric energy distribution for such cool stars. It is also a good compromise between bands too much affected by interstellar extinction at short wavelengths (very few $`\rho `$ Ophiuchi YSOs have been detected in the V, R, or I bands), and the H, K and mid-IR bands which are contaminated by intrinsic excesses. Greene et al. (gwayl (1994)) showed that there is a good correlation between the dereddened J-band flux and the stellar luminosity derived by other methods. They pointed out that in $`\rho `$ Ophiuchi this correlation is roughly consistent with a more theoretically based correlation expected for 1-Myr old PMS stars following the D’Antona & Mazzitelli (dm94 (1994)) evolutionary tracks. More recently, Strom et al. (sks (1995)) and Kenyon & Hartmann (kh95 (1995)) used the same model PMS tracks to directly convert dereddened J-band fluxes into stellar masses. We adopt a similar approach here.

#### 4.1.1 Extinction estimates from near-IR colors

The main difficulty and source of uncertainty with this method is due to the foreground extinction affecting the J-band fluxes. One must estimate the interstellar extinction toward each source in order to correct the observed J-band fluxes. We have used the observed near-IR colors to estimate the J-band extinction. The $`(JH)`$ color excess is most suitable for this purpose (e.g. Greene et al. gwayl (1994)) since the dispersion in the intrinsic $`(JH)`$ colors of CTTSs is small and observationally well determined (cf. Strom et al. strom89 (1989), Meyer et al. meyer97 (1997)).

The reddening law quoted by Cohen et al. (cohen (1981)), which is determined for the standard CIT system, should be applicable to the JHK photometry of Barsony et al. (bklt (1997)). We have thus used:

$$A_V=9.09\times [(JH)(JH)_0]$$
(1)
when the $`(JH)`$ color is available, and

$$A_V=15.4\times [(HK)(HK)_0]$$
(2)
otherwise, where we have adopted $`(JH)_0=0.85`$, and $`(HK)_0=0.55`$ from the CTTS results of Meyer et al. (meyer97 (1997)). The dereddened $`J_0`$ magnitude is then obtained as $`J_0=J0.265\times A_V`$, or $`H_0=H0.155\times A_V`$ for those stars not detected in J (with $`A_V`$ estimated from H–K). The dereddened $`J_0`$ (or $`H_0`$) can then be converted into an absolute J-band (H-band) magnitude, $`M_J`$ ($`M_H`$), via the distance modulus of the $`\rho `$ Ophiuchi cloud equal to 5.73 (for $`d=140`$ pc).

The uncertainties on $`M_J`$ and $`M_H`$ result from the typical uncertainties on the J, H, K magnitudes and on the intrinsic colors $`(JH)_0`$ and $`(HK)_0`$. With $`\sigma (J)=0.1`$, $`\sigma (H)\sigma (K)=0.05`$ (Barsony et al. bklt (1997)), and $`\sigma `$($`\left(JH\right)_0`$)$`=0.1`$, $`\sigma `$($`\left(HK\right)_0`$)$`=0.2`$ (Meyer et al. meyer97 (1997)), we obtain the following typical uncertainties: $`\sigma `$($`M_J`$$`=0.39`$ mag, and $`\sigma `$($`M_H`$$`=0.60`$ mag. In addition, the uncertainty on the cluster distance ($`140\pm 10`$ pc) induces a maximum systematic error of $`\pm 0.15`$ mag on $`M_J`$ and $`M_H`$.

#### 4.1.2 Relationship between $`M_J`$ and $`L_{}`$ for young T Tauri stars

The absolute J-band magnitude $`M_J`$ can be directly converted into a stellar luminosity $`L_{}`$ if the effective stellar temperature $`T_{\mathrm{eff}}`$ is known: $`\mathrm{log}_{10}(L_{})=1.890.4\times (M_J+BC_J(T_{\mathrm{eff}}))`$, where BCJ is the bolometric correction for the J band depending only on $`T_{\mathrm{eff}}`$. Pre-main sequence objects in the mass range $`0.1\stackrel{<}{}M_{}\stackrel{<}{}2.5M_{}`$ are cool sub-giant stars with typical photospheric temperatures $`T_{\mathrm{eff}}`$ 2500–5500 K (e.g. Greene & Meyer gm95 (1995), D’Antona & Mazzitelli dm94 (1994)). In this temperature range (0.34 dex wide), the photospheric blackbody peaks close to the J band (1.2$`\mu `$m), so that the J–band bolometric correction spans only a limited range, $`1\stackrel{<}{}BC_J\stackrel{<}{}2`$, corresponding to a total shift in luminosity of only 0.4 dex. Therefore, if we use a (geometrical) average value for the effective temperature, $`T_{\mathrm{eff}}^03700`$ K, we should not make an error larger than $`\pm 0.2`$ dex on $`\mathrm{log}_{10}(L_{})`$.

Furthermore, since bright $`M_J`$ sources tend to be of earlier spectral type than faint sources, we can in fact achieve more accurate luminosity estimates. Indeed, PMS stars are predicted to lie within a well-defined strip of the HR diagram. This is illustrated in Fig. 4 which displays model evolutionary tracks and isochrones from D’Antona & Mazzitelli (dm97 (1998)) on a $`M_J`$–log$`{}_{10}{}^{}(T_{\mathrm{eff}})`$ diagram for PMS stars with ages between $`0.1`$Myr and $`3`$Myr. We have used the compilations of Hartigan et al. (hss94 (1994)), Kenyon & Hartmann (kh95 (1995)), and Wilking et al. (wgm (1999)) to derive an approximate linear interpolation for $`BC_J`$: $`BC_J=1.65+3.0\times `$log10($`3700/T_{\mathrm{eff}}`$). Figure 4 shows that, for a given observed value of $`M_J`$, the possible range of log10($`T_{\mathrm{eff}}`$) is reduced to less than $``$ 0.15 dex, inducing a maximum error of $`\pm 23`$% ($`\pm 0.09`$ dex) on $`L_{}`$. Based on Fig. 4, we have adopted a linear relationship between $`M_J`$ and $`T_{\mathrm{eff}}`$: log$`(T_{\mathrm{eff}}/3700)=0.055\times (M_J4.0)`$ (see the heavy dashed line in Fig. 4). This leads to the following $`L_{}M_J`$ conversion:

$$\mathrm{log}_{10}(L_{})=1.490.466\times M_J.$$
(3)
The uncertainty resulting from this conversion should be on the order of 0.045 dex for $`0.01L_{}10.0\text{L}\text{}`$. We therefore estimate that the total uncertainty on $`\mathrm{log}_{10}(L_{})`$ is $`\sigma `$($`\mathrm{log}_{10}L_{}`$) $`=(0.045^2+(0.466\times \sigma `$($`M_J`$)$`)^2)^{1/2}`$ $`=0.19`$ dex.

The $`L_{}M_J`$ conversion described above cannot be applied to sources undetected in the J band. Instead, we use the H-band magnitude, along with the extinction estimate derived from the $`(HK)`$ color, but with the additional complication that the circumstellar (disk) emission cannot be neglected. Meyer et al. (meyer97 (1997)) found that the H-band circumstellar excess is on the order of 20 % of the stellar flux, on average, for the Taurus CTTS sample of Strom et al. (strom89 (1989)) (for a $`\rho `$ Ophiuchi sample, see Greene & Lada gl96 (1996)). This excess, expressed as a veiling index $`r_H=F_\nu ^H{}_{\mathrm{exc}}{}^{}/F_\nu ^H_{}`$ (e.g. Greene & Meyer gm95 (1995)), is equal to $`0.2`$. Accordingly, we have applied a systematic correction $`\delta `$$`=2.5\mathrm{log}_{10}(1+r_H)0.2`$ mag. The $`L_{}M_H`$ relationship obtained in a way similar to the J-band relation is then:

$$\mathrm{log}_{10}(L_{})=1.260.477\times (M_H+0.2)$$
(4)
The total uncertainty on $`\mathrm{log}_{10}(L_{})`$ derived from $`M_H`$ is then $`\sigma `$($`\mathrm{log}_{10}L_{}`$) $`=(0.045^2+(0.477\times \sigma `$($`M_H`$)$`)^2)^{1/2}`$ $`=0.29`$ dex.

This method can also be applied to Class III YSOs, using different values for $`(JH)_0`$ and $`(HK)_0`$. We have derived $`L_{}`$ for all Class III YSOs using the following relationships: $`A_V=9.09\times [(JH)0.6]`$, or $`A_V=15.4\times [(HK)0.15]`$; and log$`(L_{})=1.490.466\times (J0.265\times A_V5.73)`$, or log$`(L_{})=1.260.477\times (H0.155\times A_V5.73)`$ (the H-band IR excess for Class III YSOs is negligible).
The resulting $`M_J,M_H,A_V`$, and $`L_{}`$ estimates are listed in Table 3 for Class II YSOs, and in Table 4 and Table 5 for Class III YSOs.

### 4.2 Disk luminosities for Class II YSOs

Since the SED of an embedded Class II YSO peaks in the mid-IR range, the ISOCAM fluxes should be approximately valid tracers of the total, bolometric luminosities ($`L_{\mathrm{bol}}`$). To estimate $`L_{\mathrm{bol}}`$ for weak Class II sources, an empirical approach thus consists in using this $`F_\nu ^{\mathrm{MIR}}`$$`L_{\mathrm{bol}}`$ relationship after proper calibration on a sub-sample of (brighter) objects for which the luminosity can be derived by a more direct method. This approach has been adopted by, e.g., Olofsson et al. (ohk99 (1999)). Here, we have used the $`L_{}`$ estimates of Sect. 4.1 to check that a correlation is actually present between $`L_{}`$ and the mid-IR fluxes. Figure 5 displays $`F_\nu ^{14.3}`$ (corrected for extinction) as a function of $`L_{}`$ for the 104 Class II sources detected both in the near-IR and in the mid-IR range.

A correlation is found, showing that, despite some scatter, the ISOCAM fluxes can be used to give rough estimates of the stellar luminosities of Class II YSOs. This is useful for the few ISOCAM sources of our sample which have not been detected at near-IR wavelengths.

The mid-IR emission of Class II YSOs is usually interpreted as arising from warm dust in an optically thick circumstellar disk. Using a simplified disk model (e.g. Beckwith et al. beckwith90 (1990)), it is easy to show that any observed monochromatic flux in the optically thick, power-law range of the disk SED is simply proportional to the total disk luminosity $`L_{\mathrm{disk}}`$ divided by the projection factor cos($`i`$), where $`i`$ is the disk inclination angle to the line of sight. In the Beckwith et al. (1990) model, the disk is parameterized by a power-law temperature profile with three free parameters, $`T_0`$, $`r_0`$, and $`q`$, such that: $`T(r)=T_0\times (r/r_0)^q`$. Here, we have adopted $`T_0=1500`$ K, meaning that the disk inner radius is at the dust sublimation temperature, and $`q=2/3`$, corresponding to an IR spectral index $`\alpha _{\mathrm{IR}}=2/q4=1.0`$ typical of CTTS spectra.

We must, however, account for the fact that the stellar emission itself is not completely negligible in the mid-IR bands, especially at 6.7 $`\mu `$m. A simple blackbody emission at $`T_{}=3700`$ K gives

$`F_\nu ^{\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}6.7}}=0.082\times (L_{}/1\text{L}\text{})\mathrm{Jy}`$ (5)
$`F_\nu ^{14.3}=0.021\times (L_{}/1\text{L}\text{})\mathrm{Jy}.`$ (6)
Since $`F_\nu ^{14.3}`$ is significantly smaller than $`F_\nu ^{6.7}`$, only $`F_\nu ^{14.3}`$ is used to estimate $`L_{\mathrm{disk}}`$. Assuming cos$`(i)=0.5`$ ($`i=60\mathrm{°}`$) in all cases, we have computed $`L_{\mathrm{disk}}`$ for each of the Class II YSOs with $`L_{}`$ and $`A_V`$ estimates from Sect. 4.1, based on the following relationship:

$$L_{\mathrm{disk}}/1\text{L}\text{}=(F_\nu ^{14.3}\times 10^{+0.4\times A_{14.3}}F_\nu ^{14.3})/1.80\mathrm{Jy},$$
(7)
where $`A_{14.3}=0.03\times A_V`$ and $`F_\nu ^{14.3}`$ from Eq. 6. The results are given in the last column of Table 3.

According to this model, the mid-IR flux is a direct tracer of the disk luminosity, and the $`F_\nu ^{\mathrm{MIR}}L_{}`$ correlation of Fig. 5 simply expresses that $`L_{\mathrm{disk}}`$ correlates with $`L_{}`$. The origin of $`L_{\mathrm{disk}}`$ is either the release of gravitational energy by accretion in the disk, or the absorption/reprocessing of stellar photons by the dusty disk. In the latter case, $`L_{\mathrm{disk}}`$ is naturally proportional to $`L_{}`$. The fraction $`L_{\mathrm{disk}}/L_{}`$ of stellar luminosity reprocessed by the disk depends on the spatial distribution of dust. In the ideal case of an infinite, spatially flat disk, this fraction is 0.25 (Adams & Shu as86 (1986)). If the disk is flared, $`L_{\mathrm{disk}}/L_{}`$ is larger, while it is smaller if the disk has a inner hole. The theoretical $`F_\nu ^{\mathrm{MIR}}L_{}`$ correlations plotted in Fig. 5 correspond to $`L_{\mathrm{disk}}=0.25\times L_{}`$ and to three representative inclination angles. The fact that this simple model accounts for the observed correlation quite well, suggests that the disks of most $`\rho `$ Ophiuchi Class II YSOs are passive disks dominated by reprocessing. This is consistent with recent estimates of the disk accretion level in Taurus CTTSs (e.g. Gullbring et al. gullbring (1998)). The typical disk accretion rate of a CTTS is estimated to be $`10^8\text{M}\text{}/`$yr, corresponding to an accretion luminosity $`L_{\mathrm{acc}}=0.025\text{L}\text{}`$ for $`R_{}=3\text{R}\text{}`$ and $`M_{}=0.25\text{M}\text{}`$ (i.e., $`L_{}0.25\text{L}\text{}`$ at 1 Myr). In this case, the luminosity due to reprocessing is $`5`$ times larger than the accretion luminosity in the disk.

On the other hand, 37 sources (among a total of 104) are located above the passive disk model lines in Fig. 5. These are good candidates for having an active disk with an accretion rate typically larger than $`5\times (L_{}/0.25\text{L}\text{})\times 10^8\text{M}\text{}/`$yr.

Overall, we find that the median $`L_{\mathrm{disk}}/L_{}`$ ratio is 0.41 for the 93 Class II sources detected in the near-IR and with $`F_\nu ^{14.3}>\mathrm{\hspace{0.17em}15}`$mJy (i.e. the completeness level derived in Sect. 2.4). Using this ratio, we have derived rough estimates of the stellar luminosities of the 15 Class II YSOs which have no near-IR photometry (see Table 3) as follows:

$$L_{}(\text{L}\text{})L_{\mathrm{disk}}/0.41(0.97\times 1.6/0.41)\times (F_\nu ^{14.3}/1.8\mathrm{Jy}).$$
(8)
\[The factor 0.97 corresponds to a typical stellar contribution of $`3`$% at 14.3 $`\mu `$m, while the factor 1.6 is an average extinction correction at 14.3 $`\mu `$m (corresponding to $`<A_V>=17`$ mag).\]

### 4.3 Calorimetric luminosities for Class I sources

The most direct method of estimating the total luminosities $`L_{\mathrm{bol}}`$ of embedded YSOs consists in integrating the observed SEDs (cf. WLY89). However, since most of the $`\rho `$ Ophiuchi Class II and Class III YSOs are deeply embedded within the cloud ($`A_V\stackrel{>}{}10`$), only a negligible fraction of their bolometric luminosity can be recovered by finite-beam IR observations (e.g. Comerón et al. crbr (1993)). We thus do not attempt to derive calorimetric estimates of $`L_{\mathrm{bol}}`$ for these sources. In contrast, the calorimetric method is believed to be appropriate for Class I YSOs since these are self-embedded in substantial amounts of circumstellar material which re-radiate locally the absorbed luminosity (cf. WLY89 and AM94). Using our new mid-IR measurements, we have evaluated the calorimetric luminosities ($`L_{\mathrm{cal}}`$) of the 16 Class I YSOs observed in our survey. Only 7 of them have reliable IRAS fluxes up to 60 or 100 $`\mu `$m (IRS54, IRS44, GSS30, IRS43, EL29, IRS48, IRS51). For these, the median of the ratio of $`L_{\mathrm{cal}}`$(6.7–14.3$`\mu `$m) to $`L_{\mathrm{bol}}`$ is found to be 9.8, suggesting that the typical fraction of a Class I source’s luminosity radiated between 6.7 and 14.3 $`\mu `$m is $`10`$ % . Assuming that this ratio is representative of all Class I YSOs, we have derived estimates of $`L_{\mathrm{bol}}`$ for the remaining 9 weaker Class I sources (i.e., CRBR85, LFAM26, LFAM1, WL12, IRS46, CRBR12, IRS67, CRBR42, WL6). These luminosities are listed in Table 2.

### 4.4 Luminosity function of the $`\rho `$ Ophiuchi embedded cluster

Combining the $`L_{}`$ luminosities determined in Sect. 4.1 for the sources detected in the near-IR with the $`L_{}`$ estimates from $`L_{\mathrm{disk}}`$ for the sources without near-IR measurements (Sect. 4.2), we have built a luminosity function for Class II YSOs which represents a major improvement over previous studies (see Fig. 6a). In terms of $`L_{}`$, the completeness level for this population can be estimated from the $`F_\nu ^{14.3}`$ completeness limit derived in Sect. 2.4 ($`F_\nu ^{14.3}15`$mJy) using Eq. 8: $`L_{}^{comp}(\mathrm{ClassII})=0.032\text{L}\text{}`$. While the luminosity function previously published by Greene et al. (gwayl (1994)) included only 33 (bright) Class II sources and suffered from severe incompleteness below $`L_{}12\text{L}\text{}`$, our present completeness level is a factor $``$ 30–50 lower. The new luminosity function shows a marked flattening in logarithmic units at $`L_{}2\text{L}\text{}`$, well above our completeness limit. This important new feature is discussed in Sect. 5 below.

Based on the $`L_{\mathrm{cal}}`$ estimates of Sect. 4.3, a new bolometric luminosity function for the 16 Class I YSOs of $`\rho `$ Ophiuchi is displayed in Fig. 6b. The associated completeness level is derived from $`F_\nu ^{6.7}=10`$mJy and $`F_\nu ^{14.3}=15`$mJy (Sect. 2.4) using $`L_{\mathrm{bol}}/L_{\mathrm{cal}}`$(6.7–14.3$`\mu `$m)$`9.8`$ (Sect. 4.3): $`L_{}^{\mathrm{comp}}(\mathrm{ClassI})=0.017\text{L}\text{}`$. The median $`L_{\mathrm{bol}}`$ for Class I YSOs is $`1.6\text{L}\text{}`$, which is $`8`$ times larger than the median $`L_{}`$ of Class II YSOs ($`0.20\text{L}\text{}`$). The luminosities of Class I YSOs span a range of two orders of magnitude between $`0.1\text{L}\text{}`$ and $`10\text{L}\text{}`$, which is roughly as wide as the luminosity range spanned by Class II YSOs. The comparatively large value of $`L_{\mathrm{bol}}`$ for Class I YSOs is probably due to a dominant contribution of accretion luminosity as expected in the case of protostars.

In Fig. 6c, we plot the luminosity function of the 55 Class III sources that are located within the CS contours of Fig. 1 and for which we have enough near-IR data to derive $`L_{}`$ according to the procedure described in Sect. 4.1. This sample comprises 19 confirmed Class IIIs from Table 4 together with 36 candidate Class III sources from Table 5. It might be contaminated by a few background/foreground sources but is characterized by a relatively well defined completeness luminosity. From the completeness level $`F_\nu ^{6.7}=10`$mJy and using Eq. 5 with an average extinction correction corresponding to $`A_V=17`$ mag, we get $`L_{}^{comp}(\mathrm{ClassIII})=0.20\text{L}\text{}`$, which is $`6`$ times higher than $`L_{}^{comp}(\mathrm{ClassII})`$. Deep X-ray observations with $`XMM`$ should improve the completeness luminosity for Class III YSOs by an order of magnitude in the near future (cf. discussion by Grosso et al. grosso (2000) in a companion paper).

## 5 Modeling of the luminosity function

If the mass-luminosity relationship of PMS stars were a simple, fixed power-law function, the shape of the luminosity function would directly reflect the underlying mass function. Unfortunately, the mass-luminosity relation is a complex function with inflexion points and is strongly age-dependent. Deriving the mass function from the luminosity function thus requires knowledge of the stellar age distribution. The $`\rho `$ Ophiuchi embedded cluster is believed to be younger than most known star clusters (e.g. Luhman & Rieke lr99 (1999)), with a typical age on the order of, or less than, 1 Myr (e.g. WLY89, Greene & Meyer gm95 (1995)).

The 123 Class II sources identified with ISOCAM are presently the most complete sample of young stars available in the cluster. The corresponding mass function is estimated down to $`0.055\text{M}\text{}`$ in Sect. 5.1. We investigate the effect of including Class III YSOs by modeling the luminosity function of the 135 Class II and Class III objects located inside the CS contours of Fig. 1 (Sect. 5.3).

### 5.1 Constraints on the mass function of Class II YSOs

We have modeled the observed luminosity function using a two-segment power-law mass function with two free parameters, $`M_{\mathrm{flat}}`$ and $`\alpha _1`$: d$`N/`$d log$`M_{}M_{}^{\alpha _1}`$ for $`M_{}<M_{\mathrm{flat}}`$, and d$`N/`$d log$`M_{}M_{}^{1.7}`$ for $`M_{}>M_{\mathrm{flat}}`$. The high-mass range (for $`M_{}>M_{\mathrm{flat}}`$) is simply taken to be a power-law with index $`1.7`$ (whereas the Salpeter index would be $`1.35`$), in agreement with the IMFs favored by Kroupa et al. (kroupa (1993)) and Scalo (scalo (1998)). We assume that the cluster forms stars with a fixed IMF, independent of time. The evolutionary tracks of D’Antona & Mazzitelli (dm97 (1998))444We here adopt the tracks of D’Antona & Mazzitelli (dm97 (1998)) as they provide the most extensive set of PMS models available to date, including a tight grid of masses and ages. In the future, we plan to check the details of our results using other tracks such as those of Baraffe et al. (b98 (1998)) or Palla & Stahler (ps99 (1999)). are used to derive a set of mass-luminosity relations adapted to the $`\rho `$ Ophiuchi PMS population, under the following three hypotheses about the star formation history: Constant star formation rate over a) 0.5 Myr (Fig. 7a), b) 1 Myr (Fig. 7b), and c) 2 Myr (Fig. 7c). Cases b) and c) correspond to the simplest scenarios consistent with current observational constraints (e.g. WLY89). Case a) mimics a recent “burst” of star formation suggested by recent studies (e.g. Greene & Meyer gm95 (1995) and Luhman & Rieke lr99 (1999)). These three scenarios are roughly representative of our current (imperfect) knowledge of the stellar age distribution in the cluster.

The fitting analysis is performed on the 12 (logarithmic) luminosity bins above the completeness level of 0.032 L and below $`10\text{L}\text{}`$ by varying $`\alpha _1`$ and $`M_{\mathrm{flat}}`$. Due to the age ambiguity for a particular ($`M_{}`$, $`L_{}`$) in the PMS tracks for $`M_{}\stackrel{>}{}2\text{M}\text{}`$ and $`age<3`$Myr as a result of the transition from convective to radiative interiors (e.g. Fig. 4), our simplified modeling is not valid for $`L_{}\stackrel{>}{}10\text{L}\text{}`$. The best fit is shown as a heavy curve for each of the three assumed age distributions (Fig. 7a, 7b, and 7c). The flattening of the luminosity function below $`2\text{L}\text{}`$ is well reproduced in all three cases. The best-fit model in the 1 Myr case (Fig. 7b) is however not quite as good as the two others since it predicts a peak at $`L_{}0.40.5\text{L}\text{}`$ where the data show a dip. This is reflected by the largest value $`1.41`$ for the reduced $`\chi _\nu ^2\frac{\chi ^2}{Nn}`$ (with $`N=12`$ fitted data points and $`n=2`$ free parameters). The model of Fig. 7c reproduces the data somewhat better ($`\chi _\nu ^2=1.10`$) as the predicted peak luminosity moves down to $`L_{}0.20.3\text{L}\text{}`$. A marginally better fit is obtained with the 0.5 Myr model (Fig. 7a, $`\chi _\nu ^2=0.66`$). In this case, the model luminosity function has a local maximum at $`L_{}\stackrel{>}{}1\text{L}\text{}`$ as observed and the flattening for $`L_{}2\text{L}\text{}`$ is particulary well reproduced.

The best-fit values of $`M_{\mathrm{flat}}`$ along with the formal $`1\sigma `$ errors resulting from the fitting procedure are as follows: $`0.33\pm 0.04\text{M}\text{}`$, $`0.55\pm 0.09\text{M}\text{}`$, and $`0.74\pm 0.13\text{M}\text{}`$ for Fig. 7a, 7b, and 7c, respectively. These $`M_{\mathrm{flat}}`$ values follow the expected trend that the older the stars, the higher the derived masses for the same luminosities. The final error on the determination of $`M_{\mathrm{flat}}`$ is clearly dominated by the uncertainties on the stellar age distribution. We adopt the following average value: $`M_{\mathrm{flat}}=0.55\pm 0.25\text{M}\text{}`$. The best-fit values for $`\alpha _1`$ are $`0.02\pm 0.15`$, $`0.27\pm 0.13`$, and $`0.19\pm 0.13`$, respectively. There is no correlation with age in this case. We conservatively adopt $`\alpha _1=0.15\pm 0.2`$. These constraints on the mass function are valid for $`0.032\text{L}\text{}<L_{}<10\text{L}\text{}`$, which approximately corresponds to $`0.055\text{M}\text{}\stackrel{<}{}M_{}\stackrel{<}{}2\text{M}\text{}`$.

Under the same three assumptions about the stellar age distribution as above, we have also derived a range of stellar masses directly from the $`M_J`$ magnitude for each of the 123 Class II YSOs. We display the results as a mass function in Fig. 8, where the vertical error bars reflect the uncertainties induced by the various age assumptions rather than the statistical errors. The best two-component power-law mass function derived above is superposed as a heavy curve. The error bars are particulary large for the two mass bins close to 1 M, i.e., just above $`M_{\mathrm{flat}}`$. This is due to a shift of a large number of stars beyond $`M_{\mathrm{flat}}`$ when the stellar ages are varied upward. The global shape of the mass function is however fairly well determined: It is basically flat (in logarithmic units) at low masses down to 0.055 M and shows a steep decline beyond $`0.5\text{M}\text{}`$.

### 5.2 Possible deuterium feature in the luminosity function of Class II YSOs

Some of the features (peaks and dips) apparent in the observed luminosity function might be real as they are also predicted by the models (see Fig. 7). The feature expected in the luminosity function of young stars as a result of deuterium burning during PMS contraction is particularly interesting since its location is age-dependent (e.g. Zinnecker et al. zinnecker (1993)). As already noted, the model fit obtained under the 0.5 Myr scenario is particularly good as it predicts a peak at $`L_{}1.5\text{L}\text{}`$, consistent with the observations (see Fig. 7a). In the model, this $`1.5\text{L}\text{}`$ peak corresponds to an inflexion point in the mass-luminosity relation due to deuterium burning in $`0.4`$-Myr old stars with $`M_{}0.5\text{M}\text{}`$. If this peak in the luminosity function is confirmed, it will provide strong evidence that the population of Class II YSOs is particularly young in $`\rho `$ Ophiuchi.

### 5.3 Effect of Class III YSOs

There is most probably a significant population of Class III objects embedded in the cloud (e.g. Sect. 3.5). To evaluate how the presence of this population may affect our conclusions on the IMF of the cluster (see Sect. 6 below), we here consider the luminosity function of a combined sample of Class II and Class III sources. The sample comprises the 80 Class II and 55 Class III YSOs located inside the CS contours of Fig. 1. It should be complete down to $`0.2\text{L}\text{}`$. (Note, however, that only 19 of the 55 Class III sources of this sample are confirmed YSOs at this stage – see Sect. 4.4.) The associated luminosity function, displayed in Fig. 7d, is very similar to the luminosity function of Class II YSOs. In particular, it exhibits a flattening below $`12\text{L}\text{}`$.

The results of a modeling similar to that of Sect. 5.1, but performed on only 8 luminosity bins between 0.2 and 10 L, are shown in Fig. 7d. The best fits found under the same three star formation scenarios as in Sect. 5.1 are good (with $`\chi _\nu ^2=0.64`$, 0.94, and 1.15, respectively). The best-fit values for the parameters $`M_{\mathrm{flat}}`$ and $`\alpha _1`$ of the mass function are as follows: $`M_{\mathrm{flat}}=0.32\pm 0.04\text{M}\text{}`$, $`0.48\pm 0.08\text{M}\text{}`$, $`0.66\pm 0.10\text{M}\text{}`$, and $`\alpha _1=0.08\pm 0.45`$, $`0.00\pm 0.35`$, $`+0.05\pm 0.19`$ for star formation durations of 0.5, 1, and 2 Myr, respectively. Combining these results, we obtain $`M_{\mathrm{flat}}0.5\pm 0.2\text{M}\text{}`$ and $`\alpha _10.0\pm 0.4`$. Within the uncertainties these values are essentially identical to those found in Sect. 5.1 for the sample of Class II YSOs, although they apply to a more limited mass range, approximately $`0.17\text{M}\text{}\stackrel{<}{}M_{}\stackrel{<}{}2\text{M}\text{}`$.

In contrast to the Class II case, the luminosity function of the Class II$``$Class III combined sample does not show a peak at $`L_{}=1.5\text{L}\text{}`$ (cf. Fig. 7d). This may be understood if the average stellar age of the combined population is somewhat larger than that of the population of Class II sources (as is expected).

## 6 Discussion

### 6.1 The initial mass function in $`\rho `$ Ophiuchi

Since the $`\rho `$ Ophiuchi molecular cloud is still actively forming stars, the final IMF of the cluster cannot be directly measured. The masses of the stars already formed in the cluster provide only a “snapshot” of the local IMF (cf. Meyer et al. meyer00 (2000)). However, assuming that the mass distribution of formed stars does not change significantly with time during the cluster’s history (an assumption made in the models of Sect. 5.1), any snapshot of the mass distribution taken on a large, complete population of PMS objects should accurately reflect the end-product IMF.

Compared to previous investigations of the IMF in the $`\rho `$ Ophiuchi cloud based exclusively on near-IR data (e.g. Comerón et al. crbr (1993), Strom et al. sks (1995), Williams et al. williams (1995), Luhman & Rieke lr99 (1999)), the use of mid-IR photometry has allowed us to consider a much larger sample of young stars (see, e.g., Fig. 6). Since our sub-sample of Class II YSOs is fairly large (123 objects) and complete down to low luminosities, it provides an excellent opportunity to derive improved constraints on the $`\rho `$ Ophiuchi IMF down to low masses.

In a statistical sense at least, Class II YSOs are believed to represent a specific phase of PMS evolution which follows the (Class 0 and Class I) protostellar phases, and precedes the Class III phase (see Sect. 1). Due to their short lifetime ($`\stackrel{<}{}10^5`$yr), protostars make up only a small fraction of a young cluster’s population (cf. Fletcher & Stahler 1994), and can be neglected in the global mass function. Furthermore, in contrast to PMS stars, protostars have not yet reached their final stellar masses.

Class III objects are more numerous and should thus contribute significantly to the global mass distribution. Furthermore, it has been suggested that some YSOs evolve quickly to the Class III stage, perhaps as early as the “birthline” for PMS stars (e.g. Stahler & Walter 1993), and spend virtually no time in the Class II phase (cf. André et al. 1992, Greene & Meyer 1995). Since such objects cannot be identified through IR observations, their exact number and mass distribution will not be known until the results of deep X-ray (and follow-up) surveys are available (see Sect. 3.5). However, providing the (range of) evolutionary timescale(s) from Class II to Class III is independent of mass, both classes of PMS objects should have identical mass functions. The results of Sect. 5.3 do seem to support this view, as they suggest that the mass functions of the Class II and Class III samples do not differ down to $`M_{}0.17\text{M}\text{}`$.

We therefore conclude (and will assume in the following) that the mass distribution of Class II YSOs determined in Sect. 5.1 and shown in Fig. 8 currently represents our best estimate of the IMF in the $`\rho `$ Ophiuchi embedded cluster down to $`M_{}0.055\text{M}\text{}`$. As discussed in Sect. 6.2 below, this mass function applies to stellar systems rather than single stars.

### 6.2 Influence of binary stars

A large proportion ($`\stackrel{>}{}50\%`$) of field stars are in fact multiple (e.g. binary) systems (e.g. Duquennoy & Mayor dm91 (1991)). This is also true for young PMS stars, and there is a growing body of evidence that the binary fraction is even larger in a young PMS cluster like $`\rho `$ Ophiuchi than for main sequence stars in the field (e.g. Leinert et al. leinert93 (1993); Simon et al. simon95 (1995)). Since the present study is based on ISOCAM/near-IR observations which do not have enough angular resolution to separate most of the expected binaries, a significant population of low-mass stellar companions are presumably missing from the mass function derived above. These low-mass companions are hidden by the corresponding primaries.

To estimate the magnitude of this binary effect, we show three simple models in Fig. 8 which assume a population of hidden secondaries corresponding to a binary fraction $`f`$ of 50%, 75%, and 100 %, respectively. In each case, we start from a primary mass function with the two-segment power law form derived in Sect. 5.1. We then add a population of secondaries with masses distributed uniformly (in logarithmic units) between a minimum mass of $`0.02\text{M}\text{}`$ and the mass of the primary. We thus assume that the component masses are uncorrelated and drawn from essentially the same mass function (cf. Kroupa et al. kroupa (1993)).

The global mass functions resulting from addition of companions to the primary mass function are shown in Fig. 8. It can be seen that the global mass functions are similar in form to the original primary mass function. Neither the position of the break point ($`M_{\mathrm{flat}}`$) nor the slope in the high-mass range are affected by the addition of companions. However, the slope in the low-mass range ($`\alpha _1`$) steepens as the binary fraction increases. Indeed, the power-law index between $`0.055\text{M}\text{}`$ and $`M_{\mathrm{flat}}=0.55\text{M}\text{}`$ changes from $`\alpha _1=0.15`$ to $`\alpha _1=0.31`$ for $`f=50\%`$, $`\alpha _1=0.37`$ for $`f=75\%`$, and $`\alpha _1=0.42`$ for $`f=100\%`$.

In summary, if we account for uncertainties in the binary fraction ($`50\%\stackrel{<}{}f\stackrel{<}{}100\%`$), our best estimate of the single-star mass function in $`\rho `$ Ophiuchi is well described by a two-segment power-law with a low-mass index $`\alpha =0.35\pm 0.25`$ down to $`0.055\text{M}\text{}`$, a high-mass index $`\alpha _2=1.7`$, and a break (flattening) occurring at $`M_{\mathrm{flat}}0.55\pm 0.25\text{M}\text{}`$.

### 6.3 Luminosity/mass functions in individual cloud cores

Our sample of Class II YSOs is large enough that we can study the properties of the four sub-clusters associated with the dense cores Oph A, Oph B, Oph EF, and L1689S (see Fig. 1). All 123 Class II sources but 5 belong to these 4 sub-clusters. The luminosity functions of Class II objects in the individual sub-clusters are displayed in Fig. 9 along with the ‘best-fit’ model of Sect. 5.1 (solid curve). They all agree reasonably well in shape with both the total luminosity function and the model: all four luminosity functions are essentially flat over two orders of magnitude in luminosity and appear to have a peak at $`L_{}1.5\text{L}\text{}`$. The agreement is particularly good for the sub-cluster with the largest number of stars, Oph A (see Fig. 9a), but even the smallest sub-cluster, L1689S, tends to reproduce the shape of the global luminosity function on a smaller scale (Fig. 9d). In contrast, for instance, the luminosity function derived for the Chamaeleon I cloud based on ISOCAM data (Persi et al. pmo00 (2000)) differs markedly from the $`\rho `$ Ophiuchi luminosity functions. It does not exhibit any peak at 1.5 L and is consistent with an older ($`3`$ Myr) PMS population having a similar underlying mass function (see Kaas & Bontemps kb01 (2001)).

The similarity of the individual luminosity functions suggests similar distributions of stellar ages and stellar masses in each of the four sub-clusters.

### 6.4 Global aspects of star formation in $`\rho `$ Ophiuchi

After correction for unresolved binaries (assuming a binary fraction $`f=75\%`$), the total number of Class II sources (including companions) down to 0.055 M is $``$ 145. Assuming a Class III to Class II number ratio of 19/22 as found by X-ray surveys (Grosso et al. grosso (2000) – see Sect. 3.5), we infer the presence of $`125`$ Class III stars (including associated companions) in the same mass range. The typical number ratio of Class Is (plus Class 0s) to Class IIs is 18/123 suggesting an additional $`21`$ embedded YSOs. Altogether, we therefore estimate that there are currently $`N_{}291`$ YSOs down to $`0.055\text{M}\text{}`$ including $`19`$ % of brown dwarfs. The average and median masses of these objects are $`0.35\text{M}\text{}`$ and $`0.20\text{M}\text{}`$ respectively. The total mass of condensed objects (including brown dwarfs) in the cluster is thus estimated to be $`M_{}^{\mathrm{clust}}=291\times 0.35102\text{M}\text{}`$. (The brown dwarfs with $`0.055\text{M}\text{}<M_{}<0.08\text{M}\text{}`$ contribute only $`4\%`$ of this mass.) Restricting ourselves to L1688 (thus subtracting the $`10\%`$ contribution from L1689) whose average radius is approximately 0.4 pc (cf. CS contours in Fig. 1), we find $`N_{}^{\mathrm{L1688}}262`$, $`M_{}^{\mathrm{L1688}}92\text{M}\text{}`$, $`n_{}^{\mathrm{L1688}}980`$ stars/pc3, and $`\rho _{}^{\mathrm{L1688}}340\text{M}\text{}`$/pc3, where $`n_{}`$ and $`\rho _{}`$ are the stellar number density and stellar mass (volume) density of the cluster, respectively.

Adopting a conservative value of 2 Myr for the cluster age, the total mass of $`M_{}^{\mathrm{clust}}=102\text{M}\text{}`$ translates into an average star formation rate of $`5.1\times 10^5\text{M}\text{}/`$yr, corresponding to one new YSO (of typical mass 0.20 M) every $`4000`$ yr.

Lastly, we can derive the star formation efficiency (SFE) in L1688, defined as SFE$`=M_{\mathrm{star}}/(M_{\mathrm{star}}+M_{\mathrm{gas}})`$. The total molecular gas mass, $`M_{\mathrm{gas}}`$, of L1688 has been estimated to range between $`550\text{M}\text{}`$ (from C18O measurements – Wilking & Lada wl83 (1983)) and $`1500\text{M}\text{}`$ (from CS(2-1) data – Liseau et al. liseau95 (1995)). Using $`M_{}^{\mathrm{L1688}}=92\text{M}\text{}`$, we thus get SFE$`{}_{}{}^{\mathrm{L1688}}614`$ %, which is somewhat lower than previous estimates ($`22\%`$ – WLY89). Note, however, that active star formation in L1688 appears to be limited to the three sub-clusters/dense cores Oph A, Oph EF, and Oph B (see Fig. 1 and Loren et al. loren90 (1990)), where the local star formation efficiency is significantly higher: $`SFE31\%`$, using a total core mass of $`200\text{M}\text{}`$ (Loren et al. loren90 (1990)).

### 6.5 Comparison with the $`\rho `$ Ophiuchi protocluster condensations

In an extensive 1.3 mm dust continuum imaging survey of L1688 with the IRAM 30 m telescope (11′′ resolution), MAN98 could identify 58 compact starless condensations. Molecular line observations (e.g. Belloche, André, & Motte bam (2001)) indicate that the condensations are gravitationally bound and thus likely pre-stellar in nature. MAN98 noted a remarkable similarity between the mass spectrum of these pre-stellar condensations and the IMF of Miller & Scalo (ms79 (1979)).

In Fig. 10, we compare the pre-stellar mass spectrum determined by MAN98 with the mass distribution of Class II YSOs derived in Sect. 5.1. (As such, both distributions are uncorrected for the presence of close binary systems.) It can be seen that there is a good agreement in shape between the two mass spectra. A Kolmogorov-Smirnov test performed on the corresponding cumulative distributions confirms that they are statistically indistinguishable at the 95% confidence level. This supports the suggestion of MAN98 that the IMF of embedded clusters is primarily determined by cloud fragmentation at the pre-stellar stage of star formation. The fact that both the pre-stellar and the YSO spectrum of Fig. 10 present a break at roughly the same mass $`0.5\text{M}\text{}`$ is quite remarkable. A small, global shift of the masses by only $``$ 30% upward or downward in one of the spectra would make them differ at the 2$`\sigma `$ statistical level. Although in absolute terms, both sets of masses are probably uncertain by a factor of $`2`$ (due to uncertainties in the 1.3 mm dust opacity and in the cluster age, respectively), this suggests that the protocluster condensations identified at 1.3 mm may form stars/systems with an efficiency larger than $``$ 50–70%.

## 7 Conclusions

We have used ISOCAM to survey the $`\rho `$ Ophiuchi main cloud for embedded YSOs down to a completeness level of $`1015`$mJy at 6.7$`\mu `$m and 14.3$`\mu `$m. Our main findings are as follows:

1. A total of 425 point sources are detected in the $`0.7`$ deg2 field covered by the survey, of which 211 are seen at both 6.7$`\mu `$m and 14.3$`\mu `$m. The observed distribution of the flux density ratio $`F_\nu ^{14.3}/F_\nu ^{6.7}`$ is clearly bimodal with a gap separating two distinct groups of sources: a “red” group corresponding to Class I and Class II YSOs with optically thick mid-IR excesses, and a “blue” group consisting of Class III YSOs and background stars (cf. Fig. 3).
2. The red group consists of 139 cluster members, 71 of them being newly identified YSOs. Based on their mid-IR colors from 2$`\mu `$m to $`6.7\mu `$m and $`14.3\mu `$m, essentially all of these new YSOs are low-luminosity Class II objects. This brings the total number of Class II members to 123, a factor of 2 larger than previously known. Only 50 % of these Class II sources have a near-IR excess large enough to be recognizable in a near-IR color-color diagram.
3. Combining near-IR data from Barsony et al. (bklt (1997)) with our mid-IR photometry, we derive stellar luminosities for 123 Class II and 74 Class III YSOs. We also estimate bolometric luminosities for 16 Class I objects. The corresponding luminosity functions (Fig. 6) are complete down to $`0.02\text{L}\text{}`$, $`0.03\text{L}\text{}`$, and $`0.2\text{L}\text{}`$ for Class I, Class II, and Class III objects, respectively. The luminosity function of Class II YSOs is essentially flat (in logarithmic units) below $`2\text{L}\text{}`$, with a possible peak at $`1.5\text{L}\text{}`$ and dip at $`0.5\text{L}\text{}`$ seen in the four sub-regions Oph A, Oph B, Oph EF, and L1689S (Fig. 9).
4. The large proportion of Class II objects observed above the completeness level for Class III sources ($`M_{}0.2\text{M}\text{}`$) suggests that more than 50% of the $`\rho `$ Ophiuchi embedded YSO population have optically thick circumstellar disks at mid-IR wavelengths. In a majority of cases, the luminosities of these disks are consistent with pure reprocessing of stellar light. Only $`35`$% of the Class II objects have excess mid-IR luminosities suggestive of substantial disk accretion rates.
5. The luminosity function of Class II objects is well modeled by a population of PMS stars with ages in the range $``$ 0.3$``$2 Myr and a roughly flat mass distribution below $`M_{\mathrm{flat}}0.5\text{M}\text{}`$, i.e., d$`N/`$d log$`M_{}M_{}^{0.15}`$ down to $`0.06\text{M}\text{}`$.
6. Pending the results of a deep X-ray census of Class III objects, we argue that the mass distribution of Class II YSOs is representative of the emergent IMF of the embedded cluster. If we account for the presence of unresolved binaries, this emergent mass function is well described by a two-segment power law with a low-mass index $`\alpha _1=0.35\pm 0.25`$, a high-mass index $`\alpha _2=1.7`$, and a transition mass $`M_{\mathrm{flat}}=0.55\pm 0.25\text{M}\text{}`$ (where d$`N/`$d log$`M_{}M_{}^\alpha `$). We find no evidence for a sharp turnover at low masses down to at least $`0.06\text{M}\text{}`$.
7. The shape of the mass function for Class II systems is statistically indistinguishable from the mass spectrum determined at 1.3 mm by Motte, André, & Neri (man (1998)) for the pre-stellar condensations of the protocluster. This supports the conclusion of these authors that the IMF may be primarily determined by fragmentation at the pre-stellar stage of star formation. It also suggests that the 1.3 mm protocluster condensations should form stars/systems with an efficiency larger than $``$ 50–70%.

###### Acknowledgements.

S.B. was supported by an ESA Research Fellowship during his stay at the Stockholm Observatory. The authors thank the referee, Andrea Moneti, for constructive criticisms.


text:
# Retarded Learning: Rigorous Results from Statistical Mechanics

## Abstract

We study learning of probability distributions characterized by an unknown symmetry direction. Based on an entropic performance measure and the variational method of statistical mechanics we develop exact upper and lower bounds on the scaled critical number of examples below which learning of the direction is impossible. The asymptotic tightness of the bounds suggests an asymptotically optimal method for learning nonsmooth distributions.

PACS numbers: 87.10.+e, 05.20.+m, 02.50.-r

In recent years, methods of Statistical Physics have contributed important insights into the theory of learning with neural networks and other learning machines (see e.g. ). Among the most prominent discoveries of statistical mechanics in this field is the occurence of phase transitions in the progress of learning, when the number of presented example data is gradually increased and the dimensionality of data and model parameter is large.

Besides the ubiquity of phase transitions in discrete parameter models, they are typically observed when the learning problem contains symmetries which are spontaneously broken when the scaled number of examples increases beyond a critical value . Although phase transitions in neural networks have been analysed extensively by the method of replicas  it is usually hard to present a rigorous analysis (for an exception see e.g.  and the recent attempts of Talagrand ). Hence, this often precludes a digestion of the interesting results by researchers outside the community of statistical physicists working on disordered systems. Unfortunately, also other standard techniques based on asymptotic expansions  will not apply in these cases. They are only valid when the number of data is much larger than the number of parameters.

In this letter we will present a rigorous and simple approach to these problems. We combine information theoretic bounds for the performance of statistical estimators (see e.g. ) with an elementary variational principle of statistical physics . This will allow us to compute rigorous upper and lower bounds for the critical number of examples at which a transition occurs.

We will explain our method for the case of retarded unsupervised learning which has been analysed before using the replica framework (see e.g. ). The goal of unsupervised learning is to find a nontrivial structure in a set of data which reflects the properties of the underlying data generating mechanism rather than being an artefact of statistical fluctuations. The phenomenon of retarded learning describes the fact that for some high dimensional probability distributions, it is impossible at all to predict the underlying structure (usually a symmetry axis) if the (scaled) number of data is below a certain critical value. Only above this value, estimation of the structure can start.

We adopt a probabilistic, Bayesian formulation of unsupervised learning following . We model a situation where the probability distribution of the data is characterized by a single unknown rotational-symmetry direction $`𝐰^{}`$. More specifically, we assume that a set of $`t`$ data $`𝒙^t=𝒙_1,\mathrm{},𝒙_t`$, has been generated independently by $`t`$ samplings from a distribution of the form

$$P(𝒙|𝐰^{})=P_0(𝒙)\mathrm{exp}\left(V(\lambda )\right)$$
(1)
where $`P_0(𝒙)=\frac{1}{(2\pi )^{N/2}}\mathrm{exp}(𝒙^2/2)`$ is a spherical Gaussian distribution and $`\lambda 𝐰^{}𝒙`$ is the projection of the N-dimensional data vector $`𝒙`$ on the direction defined by $`𝐰^{}`$. The distribution of the projection is given by $`p(\lambda )=\frac{1}{\sqrt{2\pi }}\mathrm{exp}(\lambda ^2/2V(\lambda ))`$. In the following, averaged quantities with respect to $`p(\lambda )`$ will be denoted with an overline $`\overline{(..)}=d\lambda p(\lambda )(..)`$. Based on the set of data $`𝒙^t`$, the goal of a learner is to produce an estimate $`\widehat{P}_t(𝒙|𝒙^t)`$ for the true distribution $`P(𝒙|𝐰^{})`$. $`\widehat{P}_t`$ will not necessarily belong to the given parametric class (1). Our approach relies very much on the choice of a specific measure for the quality of the estimation. Rather than computing the overlap between an estimated direction and the true $`𝐰^{}`$, we choose a quantity which is directly measuring our ability to compress the data based on the information we have gained on the structure of $`P`$. This is related to the averaged relative entropy (Kullback Leibler (KL) divergence) between the true distribution and the estimate

$$L(\widehat{P}_t,𝐰^{})𝑑𝒙^tP(𝒙^t|𝐰^{})𝑑𝒙P(𝒙|𝐰^{})\mathrm{ln}\frac{P(𝒙|𝐰^{})}{\widehat{P}_t(𝒙|𝒙^t)}$$
(2)
where $`P(𝒙^t|𝐰^{})`$ is shorthand for the product distribution $`_{i=1}^tP(𝒙_i|𝐰^{})`$ and $`d𝒙^t=_{i=1}^td𝒙_i`$. We will further adopt a Bayesian approach where we assume that ”nature” draws the true parameter $`𝐰^{}`$ at random from a (noninformative) prior distribution $`p(𝐰^{})`$ and associate measure $`d\mu (𝐰^{})=p(𝐰^{})d𝐰^{}`$ given by the uniform distribution on the sphere with radius $`𝐰^{}^2=1`$. The case of a discrete prior will be discussed later.

The progress of the learning will be measured by the cumulative risk defined by

$$R_m(\widehat{P})=\underset{t=0}{\overset{m1}{}}𝑑\mu (𝐰^{})L(\widehat{P}_t,𝐰^{})$$
(3)
This measure of loss has a variety of important applications in information theory, game theory and mathematical finance (see e.g. ). E.g., it is proportional to the expected extra number of bits (assuming a reasonable quantization of the $`𝒙_i`$) we have to suffer in compressing the observed data when their distribution is unknown and a sequential estimate is used instead . As we will see in a moment it has also an important meaning in statistical physics.

A first attempt to study retarded learning by using bounds on (3) was undertaken by . However the bounds were too weak to give a nonzero bound on the critical number of examples below which learning is impossible.

An elementary calculation shows that the posterior probability $`\widehat{P}_t^{Bayes}(𝒙|𝒙^t)=\frac{{\scriptscriptstyle 𝑑\mu (𝐰)P(𝒙^t|𝐰)P(𝒙|𝐰)}}{{\scriptscriptstyle 𝑑\mu (𝐰^{})P(𝒙^t|𝐰^{})}}`$ achieves the minimum risk $`R_m^{Bayes}=R_m(\widehat{P}^{Bayes})`$ over all choices of estimators. Inserting this estimator into (3) and using (1) we get

$`R_m^{Bayes}`$ $`=`$ $`{\displaystyle 𝑑\mu (𝐰^{})𝑑𝒙^mP(𝒙^m|𝐰^{})\left[\mathrm{ln}𝑑\mu (𝐰)e^{_i\left\{V(𝐰𝒙_i)V(𝐰^{}𝒙_i)\right\}}\right]}`$ (4)
The last line looks very much like an averaged free energy in statistical mechanics for a system with hamiltonian $`_i\left\{V(𝐰𝒙_i)V(𝐰^{}𝒙_i)\right\}`$. Hence we can expect that useful bounds for this quantity can be derived using the standard variational principle of statistical mechanics  for the free energy

$`\mathrm{ln}{\displaystyle 𝑑\mu (𝐰)e^{H(𝐰)}}\mathrm{ln}{\displaystyle 𝑑\mu (𝐰)e^{H_0(𝐰)}}+HH_0_0`$ (5)
which bounds the free energy of a system with hamiltonian $`H`$ in terms of the free energy of a trial hamiltonian $`H_0`$ plus a correction term. The brackets $`(..)_0=d\mu (𝐰)e^{H_0(𝐰)}(..)/d\mu (𝐰^{})e^{H_0(𝐰^{})}`$ denote an average with respect to the Gibbs distribution defined by $`H_0`$. Using appropriate choices for $`H`$ and $`H_0`$, we will get both upper and lower bounds on (4).

We begin with the lower bound. We set $`H_0=_i\left\{V(𝐰𝒙_i)V(𝐰^{}𝒙_i)\right\}`$ and $`H=_i\left\{\lambda V(𝐰𝒙_i)\gamma V(𝐰^{}𝒙_i)\right\}`$ where $`\lambda ,\gamma >0`$ are variational parameters. Averaging both sides of (5) over $`P(𝒙^m|𝐰^{})`$ and $`p(𝐰^{})`$ using Jensen’s inequality in the second line, we derive the lower bound on $`R_m^{Bayes}`$

$`R_m^{Bayes}`$ $``$ $`{\displaystyle 𝑑\mu (𝐰^{})𝑑𝒙^mP(𝒙^m|𝐰^{})\left[\mathrm{ln}𝑑\mu (𝐰)e^H\right]}+m(\gamma \lambda )\overline{V}`$
$``$ $`{\displaystyle 𝑑\mu (𝐰^{})\mathrm{ln}𝑑\mu (𝐰)\left[𝑑𝒙P_0(𝒙)\frac{e^{\lambda V(𝒙𝐰)}}{e^{(1\gamma )V(𝒙𝐰^{})}}\right]^m}+m(\gamma \lambda )\overline{V}`$
$`=`$ $`\mathrm{ln}\left\{{\displaystyle _1^1}𝑑qW_N(q)\left[F_{\lambda \gamma }(q)\right]^m\right\}+m(\gamma \lambda )\overline{V}`$
where $`W_N(q)=d\mu (𝐰)\delta (q𝐰.𝐰^{})(1q^2)^{\frac{N3}{2}}`$ and

$$F_{\lambda \gamma }(q)=DxDye^{\lambda V(x)(1\gamma )V(xq+y\sqrt{1q^2})},$$

with the gaussian measure $`Dx=e^{\frac{1}{2}x^2}dx/\sqrt{2\pi }`$. This bound holds for every $`N`$ and every $`m`$.

To show the phenomenon of retarded learning we will compare the cumulative risk of the Bayes estimator to the risk of a trivial estimator which assumes that there is no specific structure in the data and always predicts with the spherical distribution $`\widehat{P}_t^{triv}(𝒙)=P_0(𝒙)`$ thereby achieving the trivial total risk $`R_m^{triv}=R_m(\widehat{P}^{triv})=m\overline{V}`$. Note, that $`\overline{V}`$ is a nonpositive quantity. We are interested in the difference $`\mathrm{\Delta }R_m=R_m^{triv}R_m^{Bayes}`$ between the trivial risk and the Bayes risk. Taking the thermodynamic limit $`N\mathrm{}`$ with $`\alpha =m/N`$ fixed, we can evaluate the integral in (Retarded Learning: Rigorous Results from Statistical Mechanics) by Laplace’s method which gives the asymptotic upper bound for $`\mathrm{\Delta }R_m`$

$`\underset{N\mathrm{}}{lim\; sup}\mathrm{\Delta }R_{\alpha N}/N`$ $``$ $`\underset{q}{\mathrm{min}}\left\{{\displaystyle \frac{1}{2}}\mathrm{ln}(1q^2)+\alpha \mathrm{ln}F_{\gamma \lambda }(q)\right\}+\alpha (\lambda \gamma 1)\overline{V}`$ (7)
For sufficiently small $`\alpha `$, the bound (7) is optimized for $`\gamma =0`$ and $`\lambda =1`$. For any potential $`V`$ having the property $`\overline{\lambda }=0`$ (ie. when the problem is not trivially learnable by computing the mean of the data) there is a critical value $`\alpha _{lb}=(1\overline{\lambda ^2})^2`$ such that as long as $`\alpha \alpha _{lb}`$ the minimizer is $`q=0`$ and $`lim_N\mathrm{}\mathrm{\Delta }R_{\alpha N}/N0`$ (see Fig. 1). Since the Bayes risk is minimal, we have $`\mathrm{\Delta }R_m0`$ and we conclude that $`lim_N\mathrm{}\mathrm{\Delta }R_{\alpha N}/N=0`$ at least for $`\alpha \alpha _{lb}`$. This proves the existence of a region of retarded learning, where even the risk of the optimal Bayes estimator is to leading order in $`N`$ as large as the risk of a trivial estimator which assumes that there is no any spatial structure at all. The bound $`\alpha _{lb}`$ agrees with the critical $`\alpha `$ obtained in the replica analysis of .

We next derive a lower bound on $`\mathrm{\Delta }R_m`$. Using the fact that $`R_m^{Bayes}`$ is the minimimum cumulative risk over any choice of estimators, for any distribution $`Q(𝒙|𝐰)`$ and estimator $`\widehat{Q}(𝒙|𝒙^t)=\frac{{\scriptscriptstyle 𝑑\mu (𝐰)Q(𝒙^t|𝐰)Q(𝒙|𝐰)}}{{\scriptscriptstyle 𝑑\mu (𝐰^{})Q(𝒙^t|𝐰^{})}}`$ we have $`R_m^{Bayes}R_m(\widehat{Q})`$. Setting $`H=\mathrm{ln}\frac{Q(𝒙^m|𝐰)}{P(𝒙^m|𝐰^{})}`$ in (5) and restricting ourselves to the class of trial Hamiltonians $`H_0`$ which do not depend on $`𝒙^m`$, it can be shown that the optimal choice is the data average $`H_0=𝑑𝒙^mP(𝒙^m|𝐰^{})H`$, for which, on average, the correction term in (5) vanishes. This yields

$`R_m^{Bayes}`$ $``$ $`R_m(\widehat{Q})`$
$`=`$ $`{\displaystyle 𝑑\mu (𝐰^{})𝑑𝒙^mP(𝒙^m|𝐰^{})\mathrm{ln}\frac{𝑑\mu (𝐰)Q(𝒙^m|𝐰)}{P(𝒙^m|𝐰^{})}}`$
$``$ $`{\displaystyle 𝑑\mu (𝐰^{})\mathrm{ln}𝑑\mu (𝐰)\mathrm{exp}\left[m𝑑𝒙P(𝒙|𝐰^{})\mathrm{ln}\frac{Q(𝒙|𝐰)}{P(𝒙|𝐰^{})}\right]}`$
We now have to find a good choice for $`Q(𝒙|𝐰)`$. For $`Q(𝒙|𝐰)`$ with a structure of the form (1), we set $`Q(𝒙|𝐰)=P_0(𝒙)\mathrm{exp}\left(U(𝐰𝒙)\right)`$. In the thermodynamic limit we get the lower bound for $`\mathrm{\Delta }R_m`$

$`\underset{N\mathrm{}}{lim\; inf}\mathrm{\Delta }R_{\alpha N}/N\underset{q}{\mathrm{max}}\left\{{\displaystyle \frac{1}{2}}\mathrm{ln}(1q^2)\alpha {\displaystyle Dx\mathrm{exp}(U_q(x))U(x)}\right\}`$ (9)
with

$$U_q(x)=\mathrm{ln}Dy\mathrm{exp}\left[V(xq+y\sqrt{1q^2})\right].$$
(10)
It is easy to see that for any $`q`$, the expression in the curly brackets of (9) is maximised for $`U(x)=U_q(x)`$. With this choice for $`U`$, we find that there exists an $`\alpha _{ub}`$ such that for $`\alpha >\alpha _{ub}`$, we have $`lim_N\mathrm{}\mathrm{\Delta }R_{\alpha N}/N>0`$ which means that now the performance of the Bayes risk is better than the trivial risk and a nontrivial estimation of the direction $`𝐰^{}`$ is possible (see Fig. 1). $`\alpha _{ub}`$ gives an upper bound on the region of retarded learning but has no simple analytical expression.

Our approach is also easily applied to a discrete prior distribution, e.g. a uniform distribution on the hypercube . Again, for small $`\alpha `$ we find a region of retarded learning $`\mathrm{\Delta }R_m^{Bayes}/N=0`$ at least for $`\alpha \alpha _{lb}`$ where $`\alpha _{lb}`$ is exactly the same as for the spherical prior.

For illustration, we apply the bounds to the simple case of a Gaussian distribution for which the integrals can be done analytically. Other distributions will be discussed in . We set $`P(𝒙|𝐰^{})=\frac{1}{(2\pi )^{N/2}(1+A)}\mathrm{exp}\left(\frac{𝒙^2}{2}+\frac{A}{2(1+A)}(𝒙𝐰^{})^2\right).`$ The data are normally distributed with unit variance in all directions perpendicular to $`𝐰^{}`$ and with variance $`1+A`$ in the direction $`𝐰^{}`$. The upper and lower bounds (7) and (9) (optmized w.r.t. $`\lambda `$ and $`\gamma `$) are shown in Figure 1 for $`A=0.5`$ for which we obtain $`\alpha _{lb}=1/A^2=4`$. We have compared the bounds with numerical simulations. Since it is hard to compute the Bayes optimal estimation algorithmically, we have used the following (suboptimal) algorithm instead. We have computed the direction $`\widehat{𝐰}(𝒙^t)`$ for $`𝐰^{}`$ which maximises the posterior probability of the data. The estimate of the distribution is given by the plugin estimator $`\widehat{P}_t^p(𝒙|𝒙^t)P(𝒙|\widehat{𝐰}(𝒙^t))`$ which has the KL divergence

$`{\displaystyle 𝑑𝒙P(𝒙|𝐰^{})\mathrm{ln}\frac{P(𝒙|𝐰^{})}{\widehat{P}_t^p(𝒙|𝒙^t)}}={\displaystyle \frac{A^2}{1+A}}(1(𝐰^{}\widehat{𝐰}(𝒙^t))^2)`$
and the cumulative entropic risk $`R_m(\widehat{P}^p)`$ can be easily approximated numerically by averaging over a large number of data sets. Figure 1 shows the difference $`R_m(P^{triv})R_m(\widehat{P}^p)`$. Since the Bayes risk is minimal, the upper bound on $`\mathrm{\Delta }R_m`$ is also an upper bound on every estimator while the lower bound is only a lower bound on the Bayes risk. We see that until $`\alpha 4`$, $`\mathrm{\Delta }R_m`$ is negative and decreases linearly. Since the slope of the dash-dotted curve is proportional to the (relative) instantaneous loss (i.e. the non cumulative risk), the plugin estimator is in the retarded learning regime, its instantaneous loss is even bigger than that of the trivial estimator.

This is due to the fact that the plugin estimator has to keep the elliptic form of the distribution $`P(𝒙|𝐰)`$ which is always different from spherical when $`\widehat{𝐰}(𝒙^t)0`$. Between $`\alpha 4`$ and $`\alpha 7.5`$, the performance of the plug-in estimator improves, the slope increases but is still negative. For $`\alpha >7.5`$, the performance of the plugin estimator is better than the trivial one and the curve starts to increase. The Bayes estimator does not have this kind of disadvantage and can take a form closer to the spherical distribution in the retarded learning region by smoothing over the parameters $`𝐰^{}`$.

Both bounds in figure 1 show the same type of asymptotic growth for $`\alpha \mathrm{}`$. The asymptotics can be found analytically by expanding both bounds (7) and (9) for $`q1`$. For a smooth potential $`V`$ both bounds give asymptotically the same logarithmic growth $`R_m^{Bayes}/N1/2\mathrm{ln}\alpha `$ which can also be obtained by well known asymptotic expansions involving the Fisher information matrix . On the other hand, our bounds can also be used when these standard asymptotic expansions do not apply, e.g. when the potential exhibits a discontinuity  of the form $`V(\lambda )=\mathrm{ln}\mathrm{\hspace{0.25em}2}\mathrm{\Theta }(\lambda )+U(\lambda )`$ with corresponding projection distribution $`p(\lambda )=\mathrm{\Theta }(\lambda )\frac{2}{\sqrt{2\pi }}\mathrm{exp}(\lambda ^2/2U(\lambda ))`$ where $`\mathrm{\Theta }(\lambda )`$ is the Heaviside function and $`U`$ is smooth, Our bounds yield the asymptotic scaling $`R_m^{Bayes}/N\mathrm{ln}\alpha `$.

The asymptotic matching of our bounds has an important consequence for computing asymptotically good approximations to Bayesian predictions. Such approximations are easily derived for smooth potentials by local expansions of the posterior distribution around its maximum (). However, this technique will obviously fail for nonsmooth potentials. On the other hand, our results show that the estimate $`\widehat{Q}`$ in (Retarded Learning: Rigorous Results from Statistical Mechanics) which uses the smooth optimizing potential $`U_q`$ (10), has the same asymptotic performance as the Bayes optimal estimate $`R_m^{Bayes}`$. The smoothness of $`U_q`$ will again enable local expansions. For example, following (10), the case $`V(\lambda )=\mathrm{ln}\mathrm{\hspace{0.25em}2}\mathrm{\Theta }(\lambda )`$ can be well estimated using $`U_q(\lambda )=\mathrm{ln}\mathrm{\hspace{0.25em}2}H(q\lambda /\sqrt{1q^2})`$ with $`H(x)=_x^{\mathrm{}}Dx`$ and where the maximizer $`q`$ of the right hand side of (9) is a function of $`\alpha `$.

In this letter, we have put the phenomenon of retarded learning first established by the replica method on a rigorous footing. Our method relies on a general information theoretic performance measure for learning probability distributions which is related to the free energy of statistical physics. A variational principle yields a controlled approximation to this quantity by providing exact upper and lower bounds which are valid for arbitrary dimensionality of the problem. We expect that this framework is flexible enough to be applicable to more complex and realistic probabilistic models. It may also be useful for constructing criteria that help to decide if structures estimated from a dataset in a high dimensional space reflect a real feature of the underlying data generating mechanism or if the result is expected to be a spurious effect of random fluctuations.

We are grateful to J.-P. Nadal for helpful discussions.