---

A&A 507, 1485-1502 (2009)

Young stars and protostellar cores near NGC 2023

B. Mookerjea¹ - G. Sandell² - T. H. Jarrett³ - J. P. McMullin^4,5

1 - Department of Astronomy & Astrophysics, Tata Institute of Fundamental Research, Homi Bhabha Road, Mumbai 400005, India
2 - SOFIA-USRA, NASA Ames Research Center, MS 211-3, Moffett Field, CA 94035, USA
3 - Spitzer Science Center, California Institute of Technology, Pasadena, CA 91125, USA
4 - Joint ALMA Observatory, Av Apoquindo 3650 Piso 18, Las Condes, Santiago, Chile
5 - The National Radio Astronomy Observatory, 520 Edgemont Road Charlottesville, VA 22903-2475, USA

Received 21 May 2009 / Accepted 14 September 2009

Abstract
Context. We investigate the young (proto)stellar population in NGC 2023 and the L 1630 molecular cloud bordering the H II region IC 434, using Spitzer IRAC and MIPS archive data, JCMT SCUBA imaging and spectroscopy as well as targeted BIMA observations of one of the Class 0 protostars, NGC 2023 MM 1.
Aims. We study the distribution of gas, dust and young stars in this region to see where stars are forming, whether the expansion of the H II region has triggered star formation, and whether dense cold cores have already formed stars.
Methods. We have performed photometry of all IRAC and MIPS images, and used color-color diagrams to identify and classify all young stars seen within a field along the boundary between IC 434 and L 1630. For some stars, which have sufficient optical, IR, and/or sub-millimeter data we have also used the online SED fitting tool for a large 2D archive of axisymmetric radiative transfer models to perform more detailed modeling of the observed SEDs. We identify 5 sub-millimeter cores in our 850 and 450 m SCUBA images, two of which have embedded class 0 or I protostars. Observations with BIMA are used to refine the position and characteristics of the Class 0 source NGC 2023 MM 1. These observations show that it is embedded in a very cold cloud core, which is strongly enhanced in NH₂D.
Results. We find that HD 37903 is the most massive member of a cluster with 20-30 PMS stars. We also find smaller groups of PMS stars formed from the Horsehead nebula and another elephant trunk structure to the north of the Horsehead. Star formation is also occurring in the dark lane seen in IRAC images and in the sub-millimeter continuum. We refine the spectral classification of HD 37903 to B2 Ve. We find that the star has a clear IR excess, and therefore it is a young Herbig Be star.
Conclusions. Our study shows that the expansion of the IC 434 H II region has triggered star formation in some of the dense elephant trunk structures and compressed gas inside the L 1630 molecular cloud. This pre-shock region is seen as a sub-millimeter ridge in which stars have already formed. The cluster associated with NGC 2023 is very young, and has a large fraction of Class I sources.

Key words: dust, extinction - ISM: clouds - infrared: ISM

1 Introduction

The southern part of the Orion B giant molecular cloud complex, L 1630, borders the large H II region IC 434, which is expanding into the molecular cloud and possibly triggering star formation. The interface between the molecular cloud and the H II region is seen as a bright north-south ridge of glowing gas with the Horsehead nebula (B33) and several smaller ``elephant trunks'' seen in silhouette against the bright nebulosity. B33 points directly towards the binary system,  Ori, which is ionizing the H II region IC 434 (Pound Reipurth & Bally; Abergel et al. 2002). Since the Horsehead is seen in absorption against the nebula, it must be at the same distance or closer to us than  Ori, which has a distance of 330-385 pc (Caballero 2008; Brown et al. 1994). Although 350 pc is a more probable distance, in this paper we adopt the distance 400 pc, which allows easy comparison to other studies of this region. The derived luminosities may therefore be overestimated by 20 to 30%, and physical sizes by 10-20%.

The first systematic search for dense cores and embedded star clusters in the L 1630 cloud complex was carried out by Lada et al. (1991a,b), who surveyed 3.6 square degrees in CS (98 GHz) and in K-band (2.2 m) down to 14 mag. They identified 4 embedded clusters, the richest being associated with NGC 2024 with over 300 embedded stars, and the poorest associated with NGC 2023, a reflection nebula just north of the Horsehead with only 21 embedded sources. A more detailed study by Depoy et al. (1990) in J, H, and K with a limiting K magnitude of 15 reduced the number of young stars in the NGC 2023 cluster to 16. Launhardt et al. (1996) surveyed the CS cores for dust emission at 1.3 mm, and found most of them to be associated with compact dust emission.

In the vicinity of NGC 2023 Launhardt et al. (1996) found faint extended emission associated with the southern part of the reflection nebula; to the west of the nebula they found a relatively bright double source LBS 34 SM 1 and SM 2 embedded in fainter extended emission. Sandell et al. (1999) serendipitously re-discovered the latter two sources while carrying out a CO survey of the Horsehead and L 1630 and found an extreme high velocity outflow associated with and driven by LBS 34 SM 1, which they called NGC 2023 MM 1. Sandell et al. (1999) carried out extensive observations of NGC 2023 MM 1 and MM 2 and found that MM 1 had all the characteristics of a Class 0 protostar with a luminosity of 10  . There are other signposts of star formation in the vicinity of NGC 2023 such as Herbig-Haro objects (Malin et al. 1987) and free-free emission sources (Reipurth et al. 2004) without any obvious optical or near-IR counterparts. These observations suggest that there may be a more embedded population of protostars hidden in the dense molecular cloud cores surrounding NGC 2023.

Bowler et al. (2009) have recently derived an infrared census of star formation in a 77 region around the Horsehead nebula. They used deep near-infrared (IRSF/SIRIUS ) and mid-IR (Spitzer IRAC) observations. Their findings support triggered star formation due to radiation driven implosion in the Horsehead, although they see no evidence for sequential star formation in the immediate vicinity of the cloud/H II region interface. They identify three bona fide and five candidate young stellar objects (YSOs) in the Horsehead area.

In this paper we report on millimeter aperture synthesis observations of NGC 2023 MM 1 and MM 2 with BIMA to better characterize the protostars and to see if we could resolve its accretion disk. We also carried out more extensive continuum mapping of the NGC 2023 field with SCUBA on JCMT to look for young deeply embedded stars, which can be too heavily obscured to be detected in the near infrared (e.g. NGC 2023 MM 1). We have also retrieved all IRAC and MIPS images from the archives of Spitzer Space Observatory; these data go much deeper and to longer wavelengths than those available from any ground based observatory. The Bowler et al. (2009) study of the Horsehead area uses the same IRAC data as we analyze in this study, but we examine a much larger field. Since we started this project, two other SCUBA studies of the Orion B South region have been published (Johnstone et al. 2006; Nutter & Ward-Thompson 2007). Both these studies are in part based on the same SCUBA observations. Although we do not go as deep as Johnstone et al. (2006), we have performed a very careful data (re-) reduction, to achieve better astrometry, sharper images and well calibrated data. This has enabled us to identify infrared, radio (and X-ray) counterparts to the sources detected in our SCUBA observations.

2 Observation

2.1 BIMA observations

The observations of NGC 2023 MM 1/MM 2 were made on Dec. 13, 1999 and April 13, 2000 using one frequency setting in the C-array configuration of the Berkeley-Illinois-Maryland Association (BIMA) interferometric array. The weather conditions were marginal in December and somewhat better in April, but both runs had poor phase stability. The correlator was split into four 25 MHz bands giving a velocity resolution of 0.37 km s^-1. Since the single dish observations by Sandell et al. (1999) indicated that NGC 2023 MM 1 has rather strong emission in deuterated lines, we decided to look for NH₂D 1 (at 85.9262703 GHz), which is typically quite abundant in cold cloud cores. We also included SO 2 , which can probe both the outflow and the cold molecular gas in the core. We kept the upper sideband free of any molecular transitions that could be expected to be excited in a cold low mass cloud core. 3C 273 was observed for 10 mn for bandpass calibration and the quasar 0609-157 for phase calibration at intervals of 30 mn. The observations of 0609-157 were done using a bandwidth of 800 MHz for 2.5 mn. Since there were no planets available during the observations we used our phase and bandpass calibrators for flux calibration. At the time of the observations 3C 273 had a flux density of 8.8 Jy, while the flux density of 0609-157 was 5.0 Jy, both with an uncertainty of 10% determined from planet observations within a month from our observing runs.

The data were reduced and imaged in a standard way using MIRIAD software (Sault et al. 1995). The uncertainty in the absolute amplitude scale is 20%. The data were imaged with weights inversely proportional to the variance in order to obtain the best signal to noise ratio. We created a continuum map by adding together the four windows in the upper sideband, which were found to be free of line emission. The synthesized beam for the continuum image is 12 7 7 0 at a position angle (pa) = 4, and the continuum image has an RMS noise of 3.6 mJy beam^-1. We easily detected NGC 2023 MM 1 with a flux density of 110 mJy. However, MM 1 was found to be more extended than what we determined from the JCMT SCUBA images, which appears implausible. We therefore also imaged our phase calibrator, 0609-157, which should be a true point source in the BIMA beam, and found it to have a Full Width Half Maximum (FWHM) of 6 3 , i.e. the beam is broadened by 1 5-2 . Adding 2 5 in quadrature to the theoretical synthesized beam, we find MM 1 to be unresolved. Since 0609-157 was observed only every thirty minutes, it is impossible to judge what the true synthesized beam size was for our observations. It is very likely that if we re-observed MM 1 in perfect sky conditions with the same spatial resolution, we would find it at most marginally resolved. In order to check whether the phase errors were so large that we would lose photons out of the main beam, we also determined the integrated flux over the central 1 1 area centered on MM 1. This gave the same flux density as we obtained from fitting a Gaussian to MM 1. We therefore believe that we retrieved all the flux, i.e. the error in flux density is dominated by the uncertainty of the assumed flux density for 3C 273 and 0609-157, both of which are variable.

We easily detected NH₂D 1 , but not SO, nor any other molecular transition in our data set. The rms in an individual channel, i.e. for a velocity resolution of 0.34 km s^-1, is 140 mJy beam^-1.

2.2 SCUBA observations

2.2.1 SCUBA scan maps

We obtained several large scan maps (8 8) of the NGC 2023 region with the Submillimeter Common User Bolometer Array (SCUBA) (Holland et al. 1999) on JCMT during SCUBA commissioning time on four nights in December, 1997. SCUBA has 37 bolometers in the long (850 m) and 91 in the short (450 m) wavelength array separated by approximately two beam widths and arranged in a hexagonal pattern. Both arrays can be used simultaneously.

For these observations we used the 850 m and 450 m filter combination in traditional scan-map mode, i.e. with continuous scanning in azimuth with a scan rate of 32  s^-1 while chopping in azimuth with a chop throw of 60 . This mode gives much flatter baselines than the more commonly used ``Emerson2'' scan maps. All four nights (Dec. 11, 17, 18 & 20) had excellent submillimeter sky conditions with a precipitable water vapor level, PWV  1 mm. We did a pointing check on HL Tau or CRL 618 before and after each map, and corrected for any pointing drift in the post processing of the data. During the first night, Dec. 11, we had large pointing drifts (>4 ) during the scan map, which resulted in smearing of point-like sources in the map, and we decided we could not use it. In total we acquired 5 maps (12 integrations), half of which were scanned roughly in right ascension (i.e. with the source near transit) and the rest of the maps were obtained with a roughly orthogonal scan direction, which was achieved by observing the source while it was rising. HL Tau was used as our primary calibrator, although calibration was additionally checked using Uranus, CRL 618 and IRC+10216. The HPBW (half power beam width) of SCUBA in scan map mode with linear regridding to a 1 grid was measured to be 8 0 at 450 m and 14 8 at 850 m from scan maps of Uranus.

The maps were reduced in a standard way using the SURF reduction package (Jenness & Lightfoot 1999; Sandell et al. 2001). After reducing each 850 m map separately, the maps were coadded and regridded onto equatorial coordinates. This map was used to determine the position for NGC 2023 MM 1. We then used this position and corrected the individual maps for residual pointing errors (shift and add) and performed a new coadd, where we also weighted each map accounting for the rms-noise in each map. The final positional accuracy is estimated to be 1 5 based on comparison of the positions of MM 1 (this paper) and MM 3 (Wyrowski et al. 2000) deduced from BIMA observations. The rms noise level in the final maps is 45 mJy beam^-1 for the 850 m map, and 0.45 Jy beam^-1 for the 450 m map,

2.2.2 SCUBA jiggle maps

We also took several long integration jiggle maps of the NGC 2023 MM 1/MM 2 field on Dec. 11 1997. These maps were reduced separately and added to earlier jiggle maps discussed in Sandell et al. (1999). A subset of the 850 m maps (only those which had pointing observations immediately before or after the map) were used to determine the position of MM 1. All maps were then corrected to this position and coadded by accounting for the rms-noise in each individual map. The jiggle maps are likely to have higher astrometric accuracy than the scan maps, due to the limited number of pointing checks (Table 1) we did for the scan maps. The rms noise level in the jiggle maps is 30 mJy beam^-1 at 850 m and 0.23 Jy beam^-1 at 450 m. The calibration accuracy is estimated to be better than 10% at 850 m and about 15% at 450 m for both the jiggle and the scan maps. The flux densities deduced for MM 1 and MM 2, which were observed in both modes, agree exceptionally well, see Table 1.

Table 1:   Results for submillimeter sources mapped with SCUBA.

All our maps were converted to FITS-files and exported to MIRIAD (Sault et al. 1995) for further analysis. In order to correct for the error lobe contribution, especially at 450 m, we have deconvolved all the maps using CLEAN and a circular model beam deduced from observations of Uranus. For the scan maps we did the deconvolution for only the three sub-fields where we have strong dust emission: the MM 1/MM 2 region, NGC 2023 and MM 3 and MM 4, and the MM 5 region. The results of this analysis are discussed later in the paper.

2.3 Ancillary molecular line data

We also have access to large ¹²CO and ¹³CO scan maps obtained as bad weather backup programme on JCMT between 1995 and 1998. The maps were obtained by scanning in RA or Dec with sampling every 5 in the scan direction and stepping by 10 in the orthogonal direction. Since the maps were obtained in wet weather conditions and over a long time period, each sub-map was made to overlap with previous maps to ensure that the calibration is uniform over the whole map. Most maps were repeated several times, resulting in integration times per spectral data points between 5-30 s. The ¹²CO and ¹³CO observations cover an area from south of the Horsehead nebula to about the southern part of NGC 2023, and target the L 1630 cloud boundary towards the IC 434 ionization front. To the east these observations stop just short of NGC 2023. The ¹²CO map is 11 14 ( ) and contains 12 000 spectra with a velocity resolution of 0.10 km s^-1. The noise level of the map is on the average 0.15 K km s^-1, although there are areas where the noise level can be five times higher. The ¹³CO map covers a slightly smaller area, 9-10 14, has 9600 spectra and a noise level similar to the ¹²CO map.

We have also retrieved, reduced, and analyzed the ¹³CO and C¹⁸O  maps from the JCMT raw data archive discussed in Johnstone et al. (2006), as well as a smaller C¹⁷O map of the MM 1/MM 2 cloud core.

Any comprehensive analysis of these data sets are well beyond the scope of this paper, but the large CO and ¹³CO maps have been very valuable to look for outflow activity or temperature enhancements in the vicinity of the sub-millimeter sources. Additional details of these data set are therefore given, at those places in the paper where we make use of them in this study.

2.4 Spitzer IRAC and MIPS (24 m) images

We extracted observations from the Spitzer Space Observatory archive (Program ID 43: an IRAC survey of the L1630 and L1641 (Orion) molecular clouds - Fazio et al.). The IRAC data were obtained in the 12 s high dynamic range (HDR) mode using four different AORs (astronomical observation requests), all with overlap. The HDR mode takes pairs of images with 0.6 and 12 s frame times at each position in all four bands. Each star in the field is therefore covered with at least eight image frames. We have processed both short (0.6 s) and long (12 s) integration basic calibrated data (BCD) frames in each channel using the Artifact mitigation software developed by Sean Carey and created mosaics using MOPEX. We used APEX to identify and extract photometry of all point sources in both the short and long integration images of the 3.6 and 4.5 m bands, while we only used the short integration images for the 5.8 and 8.0 m band.

We also extracted MIPS 24 m observations for this region from the Spitzer Space Observatory archive (Program ID 47: A MIPS survey of the Orion L1641 and L1630 molecular clouds - Fazio et al.). The MIPS 70 m image was severely saturated and could not be used, but the MIPS 24 m data were of good quality with only two stars in the southern part of NGC 2023 partially saturated. The MIPS dataset was flatfielded, overlapped and mosaiced using MOPEX.

3 Large scale overview of the NGC2023

A false-color image of the NGC 2023 region is shown in Fig. 1 (blue: 4.5 m; green: 8.0 m; red: 24.0 m). Figure 2 shows the 850 m SCUBA image overlaid as a contour map on the Spitzer IRAC 8 m short integration image. Also marked in Fig. 2 are the locations of all stars identified as young stellar objects based on color-color diagrams incorporating 2MASS, IRAC and MIPS photometry and the free-free emission sources VLA 1-3 identified by Reipurth et al. (2004). The mid-infrared and sub-millimeter images provide a broad spectrum overview of the stellar content and dust emission from the region.

The emission in the mid-infrared is dominated by the reflection nebula NGC 2023, which is illuminated by the early B-type star HD 37903. The images show an infrared dark dust lane at a PA of -20, lying 7 inside the photon dominated region (PDR) boundary between IC 434 and the L1630 dark cloud. To the south-west lies the famous Horsehead nebula, which still can be recognized in the mid-IR, although here we only see the surface layers (PDRs) outlining the nebula. There are a few strong mid infrared sources, but overall the images are characterized by extended emission from the PDR. This makes identification of discrete sources and photometry in the mid-infrared difficult. The physical appearance of the region in the mid-infrared differs quite significantly from the structures seen in low resolution IRAS and far-infrared data (Mookerjea et al. 2000), the moderate resolution 1.3 mm data (Launhardt et al. 1996) and the sub-mm SCUBA data presented here. In the far-infrared wavelengths the emission is from the HD 37903 and the bright sources directly to the south of it. Owing to their poor angular resolution, the far-infrared observations detect only a single blob of emission midway between the dust emission peak MM 3 and HD 37903. At still longer wavelengths, dust emission from the dark clouds dominates, with only the PDR, illuminated by HD 37903, being marginally visible at 850 m.

In the SCUBA images (Fig. 2) we detect emission from the infrared-dark dust ridge. This dust ridge is visible in the more extended (and deeper) 850 m SCUBA image presented by Johnstone et al. (2006), which shows that the dust ridge extends past the Horsehead in the south. (In the IRAC images the dust lane is much narrower with very narrow ``wiggly'' filaments.) We also see strong emission from the molecular cloud and PDR ridge to the south and south east of NGC 2023 and a fainter dust lane connected to the more prominent north-south dust ridge. We securely identify five sub-millimeter continuum sources named MM 1 to MM 5, detected at both 450 and 850 m, all of which appear to be in an early stage of evolution. Four of the five sources detected by us were also detected by Johnstone et al. (2006) and Nutter & Ward-Thompson (2007) using SCUBA observations. The sub-millimeter sources MM 1 and MM 2 are located in the northern part of dust ridge. The only sub-millimeter source detected in the mid-infrared is MM 3, the bright source directly to the south of the NGC 2023. V615 Ori, a heavily reddened T Tauri star lying to the north of the Horsehead nebula, though most likely surrounded by an accretion disk, is not detected in the sub-millimeter. However we detect an outflow in the ¹²CO map mostly to the south-east of V615 Ori with relatively low velocities.

4 Results

4.1 Submillimeter sources in the NGC2023 region

We have identified a total of five sub-millimeter continuum sources in the region mapped with SCUBA. These are the three sources MM 1, MM 2 and MM 5, located in the ridge (extended along north-south) lying to the west of NGC 2023 and the sources MM 3 and MM 4 located directly to the south and the south-east of the NGC 2023 nebula. In order to derive the position and sizes of the sub-millimeter sources we fitted a two component elliptical Gaussian using the task IMFIT in MIRIAD, one for the sub-millimeter source, and the other for the surrounding cloud. The fit to the broader component is mainly to provide a good subtraction of the extended emission, and is not to estimate the flux density of the surrounding cloud. For all the sub-millimeter sources the flux densities at 450 and 850 m along with positions derived from the 850 m image and source sizes are presented in Table 1.

The sub-millimeter emission ridge located to the west of the mapped region extends in the north-south direction. There is no mid-IR emission associated with this dust ridge. Detailed study of the northern part of this ridge using SCUBA images at 850 and 450 m images and CO observation resulted in the identification of MM 1 as a class 0 protostar (Sandell et al. 1999). Johnstone et al. (2006) also incorporated the SCUBA data presented here together with additional data from the SCUBA archive. The latter data have also been published by Nutter & Ward-Thompson (2007). We have completely reprocessed the SCUBA data and have better calibration and astrometric accuracy compared with previously published results (Johnstone et al. 2006; Sandell et al. 1999; Nutter & Ward-Thompson 2007). Launhardt et al. (1996) detected only the sources MM 1, MM 2 and MM 3 at 1.2 mm with the Swedish-ESO Submillimeter Telescope (SEST; beamsize of 23 ) and the IRAM 30 m telescope (beamsize = 12 ). MM 3, located directly to the south of the nebula is the brightest mid-infrared (MIR) source in the region. Surprisingly the coordinates given by Launhardt et al. (1996) for this source are significantly different from those of the 850 and 450 m sources we detect. With the exception of MM 1/ MM 2 all the other sub-millimeter sources were detected at 850 m by Johnstone et al. (2006) and Nutter & Ward-Thompson (2007), but their positions, sizes and deduced cloud masses differ considerably from our results. Both of these previous studies identified only one sub-millimeter source in the MM 1/MM 2 cloud core.

Johnstone et al. (2006) and Nutter & Ward-Thompson (2007) have identified 14 and 12 sources, respectively, in the region we have mapped using SCUBA. We can identify a few of their cores in our 850 m map, but since they are too faint to allow secure identification at 450 m, we do not list them in Table 1. In Sect. 10.2 we discuss the reasons behind the differences between our work and these previous studies.

4.2 Mid-infrared photometry using IRAC and MIPS on Spitzer

We carried out multiframe PSF photometry using the SSC-developed tool APEX on all the Spitzer IRAC images and on the MIPS 24 m images. While all long integration IRAC images are considerably affected by saturation effects, the saturation of the 8 m image is particularly severe and thus could not be corrected for by the available tools. At 8 m we therefore used only the short integration images to identify sources and derive photometry. The nebulous nature of the region and strong emission from the associated PDRs makes it extremely difficult to disentangle the sources from the surrounding clouds and to derive accurate photometry. Therefore a combination of automated routines and eye-inspection was used to extract the photometry on the IRAC and MIPS images. For sources, which APEX failed to detect at one or several wavelengths, we used the APEX user list option to supply the coordinates for the source to successfully derive a PSF fit.

In all we detected over 1000 stars in the long integration 3.6 and 4.5  m IRAC images, while we found only 95 sources in the short integration 8.0 m image. We obtained photometry of 44 sources in the MIPS 24 m image. The two brightest stars, L1630MIR-63 (MM 3), and L1630MIR-73 (star # 219 in Witt et al. 1984) in the southern part of NGC 2023, are heavily saturated in the MIPS image and their flux cannot be measured by conventional photometry. For these we used the tool developed by Jarrett to rectify saturated stars; this tool uses an extended PSF to fit the unsaturated wings of the stellar light distribution. In this way we were able to recover the integrated flux to an accuracy of 3-5% based on the repeatability of several measurements from the MIPS BCDs that contain the two saturated stars. Since we are mainly interested in young stellar objects, we took the 95 sources detected in the short integration IRAC 8 m images as the primary list. We cross-correlated this list with sources detected in other IRAC and MIPS bands, and also with 2MASS point sources. We used the following association radii : 1 for the IRAC images, 2 5 for the MIPS 24 m image and 2 for 2MASS data. Table 2 presents the coordinates of the 95 MIR sources together with the 2MASS magnitudes, Spitzer-IRAC and MIPS flux densities and a preliminary classification based on selected color-color plots.

Table 3:   Coordinates and flux densities of 2MASS & IRAC 3.6/4.5 m PMS sources in NGC 2023/L 1630.

In order to identify the stellar and pre main sequence (PMS) stars we used photometry extracted from the IRAC 3.6 and 4.5 m long integration images, which go much deeper, have a cleaner PSF, and appear to be less affected by nebular emission than the 8 m image. We identify 638 additional sources, detected in both the 3.6 and 4.5 m wavebands. Of these sources 443 are found to have 2MASS point sources associated with them. Based on a classification scheme involving the -[3.6] and [3.6]-[4.5] colors described in Sect. 4.2.1 we identify 31 of the NIR-MIR sources to be type II PMS stars (Table 3).

We find that, with the exception of L1630MIR-35 and L1630MIR-54, all the MIR sources detected at 8.0 m were detected in the other 3 IRAC bands as well. However only 34 MIR sources were detected in the MIPS 24 m image and 86 MIR sources were found to have 2MASS counterparts. Out of the 43 sources detected at 24 m, 10 were detected neither in the IRAC bands nor in the sub-millimeter images (Table 4). Six of the sources detected only in the MIPS band only are in or near NGC 2023. Of the MIPS-only sources (Table 4) L1630MIPS-4 is close to MM 2, and L1630MIPS-5 is in the vicinity of MM 1 and (probably) the young stars associated with the cold MM 1 cloud. Several other MIPS sources lie close to faint sub-millimeter peaks that we do not identify as sub-millimeter sources in this paper (Fig. 2). There are several interesting sources (VLA3, Sellgren's star C, HD 37903, and V615 Ori), which are bright in the MIR, but which were not detected sub-millimeter; see Sect. 9 for further details.

4.2.1 YSO classification based on NIR-MIR colors

We classified the sources detected in the NGC 2023/L 1630 region using the NIR and MIR photometric data presented above. Figure 3 presents the color-color diagrams derived from the 2MASS, IRAC and MIPS 24 m magnitudes of sources detected in NGC 2023 along with several criteria (shown as dashed lines in Fig. 3) used to classify them.

Based on the combined IRAC+MIPS colors ([3.6]-[5.8] vs. [8]-[24]) (see Fig. 3 (left)) we find that almost all the sources have significant infrared excess. Most of the sources that have been identified in the four IRAC and MIPS bands belong to Class II. At 24 m the photospheric colors should be close to zero for all spectral types, hence the [8]-[24] color is very sensitive to excesses.

For the sources detected in all four IRAC bands we find that both the [5.8]-[8.0] and the [3.6]-[4.5] colors for these sources extend over 1.5 mag. (Fig. 3 middle). Sources with the colors of stellar photospheres are centered at ([3.6]-[4.5],[5.8]-[8.0]) = (0,0) and include foreground and background stars as well as diskless (Class III) pre-main sequence stars. The box outlined in Fig. 3 (middle), defines the location of Class II objects in the plane (Megeath et al. 2004; Allen et al. 2004). These sources have colors that can be explained by young, low-mass stars surrounded by disks. Hartmann et al. (2005) have shown from their observations of Taurus-Auriga young stars that the color criteria [3.6]-[4.5] > 0.7 and [4.5]-[5.8] > 0.7 isolate Class 0/I protostars from Class II and Class III T Tauri stars.

Table 4:   Sources detected only in the MIPS 24 m band.

Hartmann et al. (2005) have also shown that the color criteria -[3.6] > 1.6 and [3.6]-[4.5] >0.7 discriminate well between Class II and Class 0/I systems. We have thus used the photometry for the sources presented in Table 2 which have 2MASS counterparts to generate the -[3.6] versus [3.6]-[4.5] color-color plot (Fig. 3 right). In the same color-color plot we also show the 31 (out of the 443 sources detected in the 3.6 and 4.5 m IRAC images and with 2MASS associations) which are identified to be PMS in nature based on the color criteria (Table 3). We do not plot the entire sample of 443 sources in order to avoid crowding the figure. Owing to severe contamination due to PAH emission at 3.6 m L1630NIR-23 and L1630NIR-24 do not show colors typical of YSOs, but since they are known to be H emission stars, we count them as YSOs and included them in Table 3.

A summary of the classification of sources is presented in Table 2. For most sources it was possible to arrive at a unique class based on all three color-color diagrams. For all sources the class derived based on at least two color-color diagrams match. For sources indicating two different classes from the three color-color diagrams, we have indicated both classes. Out of the 95 sources detected in the IRAC 8 m band 32 stars are found to be PMS stars (Classes I, II and I/II) based on their colors. Two sources detected only in the IRAC 5.8 and 8.0 m (L1630MIR-35 and 54) are most likely PMS, while an additional 14 stars marked as type II/III may or may not be PMS stars. The 10 MIPS-only sources (Table 4) are likely to be PMS stars, although we cannot exclude the possibility of some of these stars being background stars. Use of the more sensitive long integration 3.6 and 4.5 m images in combination with 2MASS data yields an additional 31 PMS objects (Table 3). Since the long integration 3.6 and 4.5 m images go deeper and are less affected by extinction than the 2MASS survey, we also checked how many of these IRAC sources are likely to be PMS stars. Since there are no strong PAH (polycyclic aromatic hydrocarbon) features in the 4.5 m band, and the 3.3 and 3.4 m PAH bands will be in emission rather than absorption, the [3.6]-[4.5] color should be sensitive to stars with a very red continuum, i.e. stars with infrared excess. If we use the color criterion, [3.6]-[4.5]  0.2, for sources detected only in 3.6 and 4.5 m which have no 2MASS counterparts, we find another 38 sources which are likely to be PMS stars. In Table 2 we denote the sources which exclusively occupy the region characteristic of the stellar photospheres in the color-color diagram as ``S'', although some of them may indeed be young stars.

Figure 2 shows the PMS objects identified in the region with vertical white crosses. Thus in the NGC 2023 region we identify 73 out of a total of 739 sources (111, if we include the very red sources seen only in the IRAC 3.6 and 4.5 m images) to be PMS in nature and of these five are bona fide Class I sources. Nine are identified as Class I/II. We discuss the distribution of these YSOs in Sect. 10.

5 The NGC2023 MM 1/MM 2 region

As mentioned earlier the SCUBA data for the NGC 2023 MM 1/MM 2 region is of improved quality compared to the data presented by Sandell et al. (1999). We determine the position of MM 1 more accurately with BIMA and we resolve MM 2 in both SCUBA bands. We do not detect MM 1 with either 2MASS, IRAC or MIPS. The flux densities for MM 1 and MM 2 at 450 and 850 m (Table 1) are within 10% of the flux densities measured by Sandell et al. (1999).

Sandell et al. (1999) showed that the sub-millimeter source MM 1 satisfies all the criteria for a low mass class 0 protostar. It is not detected in the near-infrared. It powers a jet-like outflow with extreme velocity (  km s^-1) detected in the CO (Sandell et al. 1999). Its outflow efficiency, defined as the ratio of the momentum flux due to the mechanical force to the radiative momentum flux ( / ), is >1000, a value typically seen in class 0 sources. A large fraction of its luminosity is emitted in the sub-millimeter with   . They found the protostellar disk to be extended in the sub-millimeter with a size 5 , i.e., diameter 2000 AU and a mass of 2-4  . Sandell et al. found the dust emission to be significantly optically thick at 850 m.

5.1 BIMA results

Figure 5 shows the contours of the 3.5 mm continuum and intensity of NH₂D integrated over the main component observed with BIMA, overlaid on the 450 m SCUBA image. We also overplot the contours of integrated intensities of C¹⁷O and C¹⁸O in Fig. 5 for reference. The sub-millimeter source MM 1 is readily detected at 3.5 mm, but not MM 2. NH₂D emission shows a distribution similar to the continuum emission and is detected from both MM 1 and MM 2. Strong emission from the deuterated molecule could be a signature of the youth of the system. The peak of the NH₂D is offset to the north by 1 5 in declination with respect to the continuum peak at MM 1 in the BIMA image. Since the relative positional uncertainty between the spectroscopic and continuum observations done with BIMA simultaneously is negligible, the detected offset is significant. Roueff et al. (2005) also found that the NH₂D emission does not seem to peak at the positions of the embedded protostars, but instead at offset positions, possibly due to interaction between the ambient cloud and the outflow. A plausible explanation could be that NH₂D gets destroyed once disks form, so that only in the pre-stellar phase the emission is peaked on the protostar. For somewhat more evolved protostars it no longer peaks at the center of the source. Further, at moderate velocities the outflow could sputter NH₂D molecules off dust grains, resulting in an enhancement of NH₂D. This of course requires that the outflow velocities are not high enough to destroy the grains completely, but are moderate as found in the shearing layer between the outflow and the ambient cloud.

Although MM 1 appears extended in our BIMA observations (Sect. 2.1), we do not know whether we resolved the continuum emission or not, because of the marginal weather conditions during our observations. The poor phase stability may also affect the position of the continuum emission; (2000.0) = 054124 929, (2000.0) = -0218 07 2, with an uncertainty of 0 8, estimated from dividing up the data into subsets and determining the position of MM 1 in each data set. This position agrees within errors with the faint free-free emission source NGC 2023 VLA 1 (Reipurth et al. 2004), offset by -0 0,+0 5 in RA and Dec, respectively. Since the free-free emission is only 0.14 mJy at 3.6 cm, the contribution from free-free emission is completely negligible at 3.5 mm, i.e. all the continuum emission at 3.5 mm is due to dust.

5.2 Is NH₂D tracing the protostellar disk of MM 1?

Based on CO 2 observations Sandell et al. (1999) concluded that NGC 2023 MM 1 drives a strong outflow at a position angle of 145  (measured from north to east). To explore the presence or detection of a disk from our NH₂D data we assume that the disk is oriented at right angles to the outflow with PA = 55. Figure 6 shows the position-velocity (PV) diagram of the NH₂D emission along the assumed disk as well as the outflow. We see a distinct velocity gradient along the disk, suggesting that the protostellar disk or the circumstellar material is rotating. The emission is slightly redshifted towards northeast and blueshifted to the southwest with marginal asymmetry relative to the position of the continuum source which is at an offset of (+1 2, +1 8) relative to the phase center for the BIMA maps.

Figure 7 shows spectra at selected positions along the disk with MM 1 at the center and at the position of MM 2. The spectra were analyzed using CLASS method HFS to fit the six hyperfine lines. The HFS method of CLASS is a part of the GILDAS software package which calculates the total optical depth , the average linewidth of lines and the radiation temperature of the multiplet, assuming LTE. Frequencies and line strengths were adopted from Tiné et al. (2000). The derived velocities in the LSR frame, , for the main component, the optical depths and linewidths resulting from the fit are presented in Table 5. Based on the fits we find that there is distinct shift in along the disk, the shift is of the same magnitude as the velocity resolution for the BIMA observations. Although the velocity resolution of the BIMA observations is 0.334 km s^-1, based on the assumption that the lines are Gaussian it is possible to determine the center of the line to at least an accuracy of better than a tenth of this value owing to the improved accuracy of the fits due to the hyperfine lines. We find that the velocity changes by 0.13 km s^-1 over a lengthscale of 3 . Most significantly the velocity change per unit length, closer to the center of the disk is similar throughout the length of the disk. If the velocity gradient was due to Keplerian motion, a smaller gradient would have been observed closer to the star. This implies that the rotating disk which probably lies hidden inside is not seen in NH₂D, the emission appears to be primarily due to the large scale infalling envelope rotating approximately as a rigid body.

The narrow line widths seen in NH₂D also support this non-detection. If NH₂D were probing a rotating disk, the emission from the disk in our relatively broad beam, should by itself broaden the line. NH₂D is quickly destroyed in shocks and in warm gas, so that it may be absent in the accretion disk and what we see is the pristine cold gas in the surrounding, slowly rotating infalling envelope. We note that although NH₂D could exist in the dense cold mid-plane of the disk, the size of the disk and the low temperatures of the mid-plane would make detection of the emission from this region extremely difficult.

We find that in comparison to the direction along the disk there is hardly any velocity gradient along the direction of the outflow (Fig. 7). This implies that the spread in velocity seen in the P-V diagram (Fig. 6) is primarily due to the broadening of the spectral lines by the outflow and the difference in velocities of MM 1 and MM 2 (Table 5).

5.3 Mass of MM 1 and MM 2 from greybody SED fits

Figure 8 shows the greybody models fitted to the observed SEDs of NGC 2023-MM 1 and MM 2. Based on the observed sizes ( ) of the sources, for the greybody fitting we have assumed to be 5 and 10 for MM 1 and MM 2 respectively. The dust temperatures (), dust emissivity indices () and masses calculated for the best fit model are presented in Fig. 8. The total dust and gas mass, has been derived from the fitted and using

where d is the distance in kpc, is the total flux (in Jy) at frequency and is given in [ ]. In this equation we have adopted a gas-to-dust ratio of 100 and standard ``Hildebrand'' mass opacity, , defined at 250 m (1200 GHz), i.e.  cm² g^-1(Hildebrand 1983).

The -index derived for MM 1, 0.7, is much smaller than what was found by Sandell et al. (1999). However, since their fit was dominated by single beam sub-millimeter photometry, which includes emission from the surrounding cloud, their -index was too high, and their fit underestimated the SED at long wavelengths. We subtract out the emission from the surrounding cloud, and since we now also have a flux density measurement at 3.5 mm, we obtain a very good fit. A of 0.7, though apparently low for a Class 0 protostar (Froebrich 2005), is not unusual. HH 24 MMS, another Class 0 source in the L 1630 dark cloud complex, has a very similar -index (Kang et al. 2008; Chandler et al. 2005). We estimate the total mass of MM 1 to be 1.04  , lower than the mass estimate of Sandell et al. (1999), primarily because we derive a much lower dust emissivity.

Given that MM 1 drives an energetic outflow the derived seems quite realistic, while the value of derived is quite low. As discussed by Mookerjea et al. (2007) both theoretically and observationally there is a large scatter in the values of . Based on current observations a low value of typically corresponds to evolution of dust grains in circumstellar disks around low mass protostars. Since the evolution of grains can change their shape, size and chemical composition, it is difficult to exactly say what causes the lower values of observed (Beckwith et al. 2000). Miyake & Nakagawa (1993) had theoretically estimated that grainsizes of 1 mm correspond to . Models of grain growth by collisional coagulation predict that while grain growth could occur on very short timescales within 10⁵ years these grains are removed very efficiently within a short timescale. This however is inconsistent with observed SEDs of classical T Tauri stars. Hence models also suggest that small grains probably get replenished, via aggregate fragmentation through high-speed collision (Birnstiel et al. 2009; Dullemond & Dominik 2005). We propose that the cause for the measured low -index is a partially optically thick inner disk, which would mimic the appearance of a low value. This is plausible since many T Tauri stars have optically thick inner disks (Andrews & Williams 2007; Beckwith et al. 1990; Andrews & Williams 2005). Hence the young disk around MM 1, with is far more massive than a Class II disk, could well be partially optically thick at mm-wavelengths.

We note that NGC 2023 MM 2 has a much steeper SED than MM 1. Our isothermal fit gives a dust emissivity, , and gives = 11.4 K, suggesting that it is a cold pre-stellar core. It has a mass of 7.4  . Although there is a faint 24 m source, L1630MIPS-4 (Table 4) 5 from the peak of MM 2, there is no sign of activity in the MM 2 core. The 24 m source is most likely unrelated to MM 2, although it could be a low-mass pre-main-sequence object embedded in the MM 1/MM 2 cloud core.

We have estimated the mass of the cloud harboring MM 1 and MM 2 using both the jiggle and scan maps observed with JCMT. For this purpose we have smoothed the maps to a resolution of 30 and then fitted a single Gaussian to estimate the total emission from the cloud and MM 1 and MM 2 at 450 and 850 m. The flux densities of the cloud from the scan maps and the jiggle maps were determined by subtracting the flux densities of MM 1 and MM 2 in the two types of maps (Table 1). Finally the 450 and 850 m flux densities of only the cloud were obtained by averaging the flux densities obtained from the scan and the jiggle maps. From the SCUBA observations we estimate S₄₅₀ and S₈₅₀ to be 58.2 Jy and 7.25 Jy respectively. Launhardt et al. (1996) using 1.3 mm observations at a resolution of 30 obtained S₁₃₀₀ to be 1.87 Jy. Using single temperature greybody models we derive a temperature of 16.5 K and an emissivity index . Using these parameters we estimate the mass of the cloud in which MM 1 and MM 2 are embedded to be 18  .

Thus, from the present SCUBA observations we derive a total mass of 26  for the entire cloud including MM 1 and MM 2. In contrast, Johnstone et al. (2006); Nutter & Ward-Thompson (2007), who both identified only one core in the MM 1/MM 2 cloud complex, using SCUBA studies, derived a total mass of 14.7 and 12.7  , respectively, for the entire region.

Table 5:   CLASS HFS fitting of selected NH₂D spectra. MM 1 is at offset (1 2,1 8).

5.4 Depletion and Deuteration scenario in MM 1 and MM 2

In order to derive some of the relevant physical parameters, we analyzed the NH₂D spectral profiles using the HFS method in CLASS. This method derives the total optical depth  and average width of the lines in an LTE approximation along with the radiation temperature of the multiplet, . Following Tiné et al. (2000) we have estimated the excitation temperature for the multiplet using the following equation

where is the radiation temperature of a blackbody at temperature T and = 2.7 K (the temperature of the cosmic background radiation).

The total column density for NH₂D is derived using the optical depth, excitation temperature and average FWHM as follows

where is the frequency of the unsplit transition. The values of the partition function, , used for the column density calculations, are 4.47 and 9.37 at 5 K and 10 K respectively (Roueff et al. 2005).

Table 5 also presents the estimated and the column densities calculated for = 5 K and 10 K. The method HFS assumes a beam filling factor of unity, this possibly explains the rather low excitation temperatures that we derive at all positions. The NGC 2023 MM 1 and MM 2 cloud core is extremely rich in deuterium, since we observe optically thick emission at many positions (Table 5). The optical depth and the column density of NH₂D along the disk shows a centrally peaked profile, with the peak being offset from MM 1, while along the outflow the NH₂D column density increases towards the southeast. However keeping in mind the assumptions and uncertainties in the estimation of these column densities variation by factors of 2 to 3 are most probably not significant. Overall, the column densities are similar in magnitude (10¹⁴ cm^-2) to the numbers found in dark cores by Roueff et al. (2005).

The observed high abundance of singly deuterated ammonia in MM 1 is suggestive of freezeout of molecules like CO, since depletion of the heavy molecules results in an increase in the [H₂D⁺]/[H₃⁺] and molecular D/H ratios. Using C¹⁸O and C¹⁷O data from the SCUBA archives we have investigated the depletion scenario in MM 1 and MM 2. The C¹⁸O data (dotted contours in Fig. 5) does not show any enhancement of emission at the positions of MM 1 and MM 2. The C¹⁷O emission appears to peak at a position in between MM 1 and MM 2.

We estimate the molecular hydrogen column densities, N(H₂), for both MM 1 and MM 2 using the C¹⁸O and C¹⁷O observations. For this we have assumed that the LTE approximation holds and have assumed temperatures of 20 K and 11.5 K for MM 1 and MM 2 respectively. We note that the estimated column density changes only by a factor of 2 for a change in temperature from 10 to 100 K. The choice of gas temperatures is based on the dust temperatures estimated from the greybody fits (Sect. 5.3). Using C¹⁸O (C¹⁷O) we estimate N(H₂) to be  cm^-2 (  cm^-2) in a 22 beam for MM 1. For MM 2 we estimate the N(H₂) to be  cm^-2 (  cm^-2) from C¹⁸O (C¹⁷O.) The discrepancy in the N(H₂) calculated from C¹⁸O and C¹⁷O column densities is due to C¹⁸O being optically thick. The relative abundances of NH₂D at the positions of MM 1 and MM 2 are 10^-9 (Table 5).

Alternatively N(H₂) can also be derived from the 450 and 850 m fluxes observed at MM 1 and MM 2. Using the observed flux densities and dust temperatures from the greybody fits we get N(H₂) in the 22 beam to be  cm^-2 and  cm^-2 for MM 1 and MM 2. From dust emission within a beam of 9 4, comparable to the beam of NH₂D, we estimate N(H₂) to be  cm^-2 and  cm^-2 for MM 1 and MM 2 respectively. Thus N(H₂) estimated from dust continuum emission is larger by a factor of 10 than the N(H₂) estimated from the CO isotopomers. This is a clear evidence of depletion of CO and its isotopomers by a factor of at least 8-10. This is similar to the depletion factors observed in prestellar cores, which lie between 5.5 and 15.5 (Bacmann et al. 2002).

6 The molecular cloud south of NGC 2023

The reflection nebula NGC 2023 is bounded to the south and southeast by a dense, extended molecular cloud coinciding with the ridge of strong PAH emission seen in the IRAC and MIPS images. At 850 m the cloud has a size of 98 78 PA 34 measured at the 3- contour level with an integrated flux of 16.7 Jy. Some faint, diffuse dust emission is also seen toward the reflection nebula. The strongest dust emission is just outside or overlapping with the bright PDR ridge seen in IRAC and MIPS images and the ionized carbon (C⁺) ridge. The latter was mapped by Wyrowski et al. (2000) with the VLA in the C91recombination line. At high resolution the dust emission breaks up into several cores (Fig. 10). Both Johnstone et al. (2006) and Nutter & Ward-Thompson (2007) identified four cores in this cloud. We can securely identify only two at both 450 and 850 m. One of them (MM 3) is compact and coincides with a bright mid-IR source, which has the colors of a Class I/0 source (Sect. 4.2). The other core, MM 4, is much more extended and has no near- or mid-IR counterpart. Even though we see 850 m emission at the approximate positions of the other two cores seen by Johnstone et al. (2006) and Nutter & Ward-Thompson (2007), they have no counterparts at 450 m, and we do not consider them as cores. Since the Johnstone et al. (2006) map goes deeper, our failure to see these cores could be due to our limited sensitivity. It is, however, possible that what Johnstone et al. and Nutter & Ward-Thompson identify as cores, are structures that break up into smaller clumps at the higher angular resolution at 450 m, especially since neither core can be identified in any molecular line map.

The molecular gas in the PDR region is hot. The ¹³CO map shown in the Johnstone et al. (2006) paper and re-reduced by us from the JCMT archive show ¹³CO brightness temperatures of 50 K towards MM 4, which means that the gas temperature is 50 K. Such high gas temperature agrees well with the color temperatures of 40-50 K derived from far infrared observations at 50 and 100 m (Harvey et al. 1980) one arcminute south and west of HD 37903. The color temperatures determined between 160 and 100 m are somewhat lower, 30-40 K. The latter color temperatures agree well with Mookerjea et al. (2000), who determined dust temperatures of 30 K from 205 and 138 m over most of the cloud south of NGC 2023. The hot gas is probably restricted to a narrow zone between the reflection nebula and the molecular cloud. If we assume a temperature of 30 K for the cloud as a whole, we find a cloud mass of 30  from C¹⁸O , while the mass seen in 850 m dust emission, which only probes the densest parts of the cloud core, is only 14.3  . Both mass estimates are for the continuum and line fluxes from a region extending over 98 78 . The total continuum emission over this region is 16.7 Jy, from which 1.3 and 1.6 Jy due to the emission from MM 3 and MM 4 are subtracted to get a flux of 13.8 Jy for the cloud only. The average dust temperature is assumed to be 30 K, i.e. similar to what we assumed for C¹⁸O. In this case the cloud mass derived from C¹⁸O is about twice as large than what we derive from dust emission. If C¹⁸O is undepleted, we would expect to derive a larger mass from the C¹⁸O emission than from dust emission, since C¹⁸O is also excited in low density gas, while the dust emission only traces the very dense gas in the cloud core. If the gas temperature is even higher, which could well be the case, it would make the difference in mass even larger. This strongly suggests that in contrast to the MM 1/MM 2 region, C¹⁸O is not depleted in the MM 3/MM 4 region of NGC 2023.

7 The NGC2023 MM 3/MM 4 region

7.1 MM 3

MM 3, coincides with an extremely red IR source, originally discovered by Depoy et al. (1990), #3 in their list of near-IR sources associated with NGC 2023. The source is also in the 2MASS (although not listed by SIMBAD) and the MSX point source catalogues. We detected it in all IRAC bands and in the MIPS 24 m image. Wyrowski et al. (2000) detected MM 3 at 3 mm with a flux density of  mJy, although they interpreted the dust emission to originate in an extended cloud core surrounding the near-IR source detected by Depoy et al. (1990). We find the dust emission to be more compact; at 450 m the emission appears unresolved. The 850 m size quoted in Table 1 is therefore an upper limit. The excellent spatial agreement between our sub-millimeter position for MM 3 and the 2MASS and IRAC position (0 5) therefore supports the classification from the IRAC colors, which suggest that MM 3 is a Class I/0 protostar.

Even though MM 3 is not associated with free-free emission (Reipurth et al. 2004), it may excite the HH objects HH 5, HH 4, and HH 1 A-C, which were discovered by Malin et al. (1987). Examination of the ¹³CO  map shows faint high-velocity wings near the position of MM 3 and faint red-shifted wings to the east of the source. If the outflow is bipolar the blue-shifted molecular outflow would therefore be to the west, i.e. in the direction towards the HH objects. There is, however, considerable velocity structure in ¹³CO in the NGC 2023 PDR region, and this is by no means a clear identification. L1630MIR-62 (Sellgren's star C, see Sect. 9) is somewhat better aligned with the HH objects. Another possibility is a very red near IR source, L1630NIR-15, which is within 1 of HH 1 A.

Because MM 3 lies within the boundary of the reflection nebula, which is very bright in the far infrared, there are no flux density estimates or even stringent upper limits at 100 m, which could help us constrain the temperature of the dust emission. However, MM 3 is a visible IR object, presumably warmer than 20 K (MM 1) but probably not too hot. We therefore assumed a dust temperature of 25 K and estimated the dust emissivity and dust mass with a single temperature greybody fit using the observed sub-millimeter and millimeter flux densities. We obtain a dust emissivity index of 0.8 and a dust mass of 0.33  . Although the dust emissivity is on the low side, most protostellar objects show dust emissivities in the range of 0.7-1.5, and often similar ``disk'' masses as well (Dent, Matthews and Ward-Thompson 1998; Froebrich 2005), re-inforcing our hypothesis that the dust emission originates in a disk/envelope around a young protostar.

Below we analyze in more detail what we can learn about MM 3 with the assumption that the dust emission originates in a disk or compact envelope surrounding the protostar.

7.2 Nature of MM 3 from SED fitting

For a more detailed model which explains the emission from MM 3 over a wide range of wavelengths we have combined the SCUBA data with the 2MASS, Spitzer (IRAC, MIPS, and IRS) and MSX data and measurements at 3 mm (BIMA). Robitaille et al. (2007) have created an archive of two-dimensional (2D) axisymmetric radiative transfer models of protostars calculated for a large range of protostellar masses, accretion rates, disk masses and disk orientations. This archive also provides a linear regression tool which can select all model SEDs that fit the observed SED better than a specified . Each SED is characterized by a set of model parameters, such as stellar mass, temperature, and age, envelope accretion rate, disk mass, and envelope inner radius. We have used this online tool to generate models which fit the observed SED for NGC 2023 MM 3. The models assume that stars of all masses form via free-fall rotational collapse to a disk and accretion through the disk to the star; thus the envelope accretion rate sets the envelope mass. Since the distance to the source is somewhat uncertain, we have explored a range of distance between 0.3 and 0.5 kpc and between 0 and 100 mag. Aperture radii used for the fits are 30 for the MSX data, 5 for the IRAC and IRS data, 9 for the MIPS 24 m data, 14 for the SCUBA data and 13 for the 3 mm BIMA data. Left panel of Fig. 9 shows the SED of the best-fitting model for MM 3 and one more model which fits the data with - .

We have obtained the average values and variances of the fitted parameters by taking a weighted average of the two models for which - , in a manner similar to Simon et al. (2007). Both these models correspond to a distance D = 350 pc. Here we present the parameters obtained from the two best-fitting models, viz., the stellar mass (), temperature () and luminosity (), the envelope mass , disk mass , and the envelope accretion rate . The numerical values are :   ,   ,  K,   ,   , and    yr^-1. The envelope mass matches reasonably well with the mass estimated from the simple greybody fit. The SED fit is very well constrained in the MIR wavelengths owing to the IRS datapoints. As a result we are able to arrive at only two models (as opposed to a family of models), which reproduce most of the observed features reasonably well, with the best fit model fitting the longer wavelength data significantly better. We note that there is a substantial discrepancy in the flux densities measured in the MSX E-band and MIPS 24 m, although both bands have very similar wavelengths. The best-fit model presented here underestimates the MIPS 24 m flux. We have checked that the quality of the fit does not change significantly by not considering the MSX fluxes, since in the MIR the fit is primarily constrained by the IRS fluxes. Thus, only more datapoints in the far-infrared region will enable us to arrive at a better constrained fit particularly at the longer wavelengths of MIR to FIR.

7.3 MM 4

Based on ¹³CO 1 brightness temperatures 50 K, MM 4 appears warmer than MM 3. However in the absence of far-infrared fluxes it is not possible to derive the dust temperature by performing an isothermal fit to the 450 and 850 m fluxes unless the value of is accurately known. We therefore considered two isothermal fits constraining to 30 K and then to 50 K. In both cases we obtained a fitted value of 1.5 for . The masses corresponding to temperatures of 30 and 50 K are 0.87  and 0.46  respectively. Irrespective of the temperature, the value of for MM 4 is lower than , which is typically seen in colder pre-stellar cores. This could be due to heating of the core from the outside, which would suggest that the core is not well described by a single temperature. Another possibility is that the dust has undergone processing and is different in the PDR region. Which explanation is the right one, or whether it is a combination of both, will have to wait until we can observe MM 4 with higher spatial resolution at several wavelengths in sub-millimeter and far infrared.

8 NGC 2023 MM 5 region

MM 5 is a compact object in the southern part of an arc or a shell-like structure in the southern part of the dust ridge (Fig. 11). There are no optical, infrared or radio sources associated with MM 5 region. Examination of our ¹²CO and ¹³CO  maps show no distinctive features in this area and the ¹²CO brightness temperature is 15 K, suggesting that the dust emission is associated with a cold prestellar core. An isothermal greybody fit to the integrated flux densities suggest a dust temperature of 12.2 K, , and a mass of 4.3  . This core was also picked up by Johnstone et al. (2006) and Nutter & Ward-Thompson (2007), who find approximately the same mass. Both identified the dust emission in the northern part of the arc as a separate cloud core (Fig. 11). Johnstone et al. (2006) identified an additional three sub-millimeter ``clumps'' in this region. It is impossible to judge whether all of these are true cores, without observing the region in a molecular line tracer sensitive to dense, cold gas.

9 Other sources of interest in the NGC 2023 region

Comparison of high resolution mid-infrared and sub-millimeter data with medium resolution far-infrared data suggests that the emission from L1630, the molecular cloud associated with NGC 2023 is primarily extended and diffuse emission. However we also find a number of embedded Class I and II objects (see Sect. 4.2), some of which we would have expected to detect in the sub-millimeter, since they have clear IR excess and at least some are expected to be surrounded by a circumstellar disk and/or envelope. Yet we detect only one MIR source in the sub-millimeter. Below we highlight some of the stars for which our IRAC and MIPS photometry provide new information about their nature, viz., HD 37903, V615 Ori, Sellgren's star C, and NGC 2023-VLA3.

HD 37903 or L1630MIR-66 illuminates the reflection nebula NGC 2023 located in the northern part of the L 1630 molecular cloud. SIMBAD lists the spectral type as B1.5 V. This star has also been classified as B2 V (Crawford 1958; Stebbins et al. 1940). Warren & Hesser (1977) assign it a spectral class between B 1.5V to B2 V. and note that it may be variable in V, which is confirmed by examining the measurements in the the SIMBAD data base. Stephenson (1986) identified HD 37903 as a weak H emission line star in an objective prism survey of the northern sky. Yet Hernández et al. (2005) did not classify it as a Herbig Ae/Be star in their study of early type stars in nearby OB associations. Since HD 37903 illuminates a reflection nebula, is an early B-star, is photometrically variable and as we show below, has an infrared excess, it appears to have all the characteristics of a Herbig Ae/Be star (Thé, de Winter & Pérez 1994). Although it is a very bright star (V = 7.8^m) and extremely well studied, almost all spectral classifications of the star are based on photometry. There are, however, several long integration UV spectra of excellent quality available in the archives of the International Ultraviolet Explorer (IUE), which enabled us to perform an improved spectral classification. We therefore looked at the short wavelength high dispersion (SWP) spectra in the IUE archive and used the spectral classification of B stars developed by Rountree & Sonneborn (1991) for ultraviolet spectra to derive a more accurate spectral classification. This comparison shows that the star is definitely later than B1.5 and earlier than B2.5, i.e. a B2 V star is a very good match. We therefore conclude that HD 37903 is a Herbig Ae/Be star of spectral class B2 Ve.

To verify that the star has an infrared excess we compared the observed infrared fluxes to the Kurucz stellar atmosphere model corresponding to solar metalicity, [g] = 3.0 and an effective temperature of 21 000 K (Fig. 9 middle). We find that while the optical and NIR fluxes are satisfactorily reproduced by the Kurucz model, the IRAC fluxes from 4.5 m show clear infrared excess compared to the model and a very large excess at 24 m possibly indicating the presence of a protostellar disk. Although the star is clearly seen at 24 m, the measured flux density is very uncertain due to the intense PDR ridge just south of it. Since HD 37903 illuminates a faint H II region (Wyrowski et al. 2000), some of the 24 m emission also originate in the H II region rather than from a disk. Although HD 37903 has an IR excess, we see no trace of it in the sub-millimeter continuum. This is not surprising, since disks around early B-stars are likely to be short lived. We do not know of any early B-star, that illuminates a large reflection nebula and still would be associated with a cold protostellar disk (see e.g. Fuente et al. 2001,1998).

V615 Ori or L1630MIR-52 is a heavily reddened T Tauri star with spectral type of G7 III-IV (Malin et al. 1987). It is one of the few PMS stars in the field that is a genuine IRAS source (IRAS 05388-0224) and not a false detection, with flux densities of 1.45 Jy at 12 m and 1.67 Jy at 25 m. It is listed in the 2MASS point source catalogue and we detect it in all IRAC bands and in the MIPS 24 m image. Malin et al. (1987) quote an observed B-V color of 2 (suggesting E_B-V = 1.2). In the General Catalogue of Variable stars (Kukarnin et al. 1971) the listed V magnitudes vary between 13.8 and 15.9. These values agree well with more recent data of Droege et al. (2007), who give and . The star is also an X-ray source, detected by the XMM-Newton telescope (XMM-Newton 2nd Incremental Source Catalogue). Using the optical and infrared flux densities and the online SED fitter we find 7 models for which - . However 6 out of the 7 models have stellar temperatures >6000 K, which is inconsistent with the known spectral type of V615 Ori. We thus present here the only model which corresponds to a stellar temperature agreeing with previous observations. The parameters for this model are: Distance = 417 pc,   ,  K,   ,   and   . If we use the power-law dependence between disk mass and 850 m flux density derived from 153 T Tauri stars in Taurus-Auriga (Andrews & Williams 2005), we find that the modeled disk and envelope mass would result in an 850 m flux density of 4 mJy, which is below our detection limit. It is therefore not surprising that we did not detect the star in our SCUBA observations nor would it have been detected in the deeper 850 m image by Johnstone et al. (2006).

Malin et al. (1987) found two HH objects (NGC 2023-HH 2, and HH 3) southeast of V615 Ori. The proper motions of HH 2, however, suggest that the exciting source is south or southeast of HH 2 rather than in the direction of V615 Ori. There is a very red MIR source, L1630MIR-51 (Table 2), south-east of HH 2, which could be the exciting star for HH 2, and possibly HH 3. L1630MIR-51 has very strong infrared excess and colors of a Class II/I object, which indicate that it is very young and deeply embedded. Examination of the large CO map to which we have access, indicates that V615 Ori may drive a weak, extended CO outflow. A position-velocity plot through the star and HH 3 shows faint blue-shifted emission to the south-east, i.e. in the direction of HH 3, and red-shifted emission to the northeast, with a reversal of the line wings approximately at the position of V615 Ori. Nothing stands out in the CO map at the position of HH 2 or L1630MIR-51. We therefore believe that HH 3 is excited by V615 Ori, although we caution that the identification is tentative and has to be confirmed with deeper observations with better spatial resolution. HH 2 may be part of the V615 Ori outflow, although we cannot exclude that this HH complex is excited by L1630MIR-56, which is 30 -40 south east of HH 2 A,B.

Sellgren's star C (L1630MIR-62, HBC 500) was also found to be a T Tauri star by Malin et al. (1987), but with much stronger H emission than V615 Ori. However, since the star is heavily obscured, they were able to only conclude that the star has a spectral type of late G or early K and a luminosity class V, based on weak Na I and Mg I absorption. Weaver & Babcock (2004) also found it on a very deep CTIO objective prism plate from 1976. In their catalogue it is number 54, and they characterize the star as having medium strength H emission with weak continuum. The star is rather bright in IRAC and MIPS images, L1630MIR-62, (Table 2), and has MIR colors of a Class II object.

NGC 2023-VLA 3 or L1630MIR-46 was the strongest radio source (50 mJy) in the 3.6 cm VLA survey by Reipurth et al. (2004). It has no optical or sub-millimeter counterpart, but is detected in all IRAC bands and also in the 24 m MIPS image (Table 2) with colors of a Class II or borderline Class I source. The flux density is too high to be free-free emission from an ionized jet. It is almost certainly gyrosynchrotron emission, which is most often found in weak-line (naked) T Tauri stars, but sometimes in young stellar objects with strong magnetic fields (André et al. 1992). Since the star is an X-ray source (XMM-Newton 2nd Incremental Source Catalogue), it is likely to have magnetic activity.

Additional sources: Weaver & Babcock (2004) identified seven additional H emission line stars within our field of view, WB 53 - 59. WB 53 is most likely the Herbig-Haro object HH 5 all the others have MIR counterparts and all but two, WB 57 and WB 59, have clear IR-excess. VLA 2 has no near- or mid-IR counterpart and is most likely an extragalactic source.

10 Discussion

10.1 The YSO population

Extracting photometry of ``faint'' stars in a bright knotty reflection nebula like NGC 2023 is rather involved. The Spitzer photometry package APEX recovered far more sources than what we present in Table 2. However, many of the ``detections'' had blue [5.8]-[8.0] colors and very red [4.5]-[5.8] colors, i.e. they appeared non-stellar. Visual inspection of each source showed that they all appeared nebulous and fuzzy. It is quite likely that some of them may contain real stars. The PSF photometry, however, was definitely dominated by nebular emission and/or absorption present in the broad IRAC filters, and we therefore rejected them. The final list may still contain some spurious sources, but if there are, they are very few.

With a total of 73 YSOs identified in the NGC 2023/L 1630 the sample is sufficiently large to visualize the distribution of the YSOs (Fig. 2). We find that the majority of them, 23 stars, are located in the southern part of NGC 2023. There is another 7-9 young stars in the sub-millimeter ridge west of NGC 2023, in which the Class 0 protostar, MM 1, is embedded. If we consider these to belong to the NGC 2023 cluster, we therefore have a total of 32 YSOs. We also find a string of YSOs (all Class II) along the whole submillimeter ridge, which runs parallel to the bright PDR rim IC 434. A slightly larger fraction of the YSOs lie on the western side of this dust ridge, i.e. the side facing the large IC 434 H II region. The rest of the YSOs lie in the bright rim created by IC 434 or slightly inside IC 434. Here we observe two distinct groups. One group is associated with the Horsehead nebula, where we see 7 sources just outside the head region and another two down the neck, where the Horsehead nebula merges into the L 1630 molecular cloud. Another small group is seen in the northwestern part of our field. The latter group is centered around an eroded elephant trunk structure, which is now much smaller than the Horsehead. A few PMS stars (including VLA 3) are scattered all over the L 1630 dark cloud, suggesting that occasional star formation can occur in any dense cold cloud core without any external triggering. Of the 5 bona fide Class I/0 stars from Table 2 and the 10 very young MIPS-only sources, a total of 9 are concentrated in the southern part of NGC 2023 and near or in the MM 1 cloud core. The remaining 6 are distributed in a scattered fashion close to the Horsehead nebula and the T Tauri star V615 Ori. The strong concentration of Class I/0 sources near NGC 2023 suggests that this is where the youngest stars are and where star formation may still be on-going.

Our YSO sample includes all the H emission line objects from the Weaver & Babcock (2004) survey. However the two emission line objects, WB 57 and WB 59, which are inside the NGC 2023 reflection nebula and are to faint to be detected at 8.0 m, but are readily detected at 3.6 and 4.5 m. Since both lie in areas of intense PAH emission, it is possible that the neutral or slightly blue [3.6]-[4.5] color is due to intense PAH emission in the 3.6 m band. As we have no way of confirming this at the present, we have counted them as members of our candidate YSO population.

Our study re-enforces the view that a large H II region like IC 434 expanding into dense surrounding cloud complexes triggers the formation of new stars. The densest regions in the surrounding molecular clouds, like the Horsehead nebula, can even form small clusters. The Horsehead may even form additional stars, because we found one Class I/0 source near the base of the Horsehead. There is also a bright, extended pre-stellar core in the neck of the Horsehead, which is seen as an infrared dark cloud in IRAC and MIPS images. The pre-shock from the H II region can compress molecular gas even further into the cloud and trigger another wave of star formation. This is demonstrated by all the young stars we find in submillimeter ridge, which runs parallel to the IC 434 boundary, but almost 0.8 pc from the edge of IC 434.

A comparison of our results with the recent study to derive an infrared census of YSOs within a 7  7 region around the Horsehead nebula (Bowler et al. 2009), shows a very good agreement with the results we get from IRAC photometry. We have identified almost all the sources with an IR excess, which Bowler et al. (2009) present in their Table 1. The sources which we do not detect are B33-13, B33-20, B33-27, B33-29, and B33-30 all of which were selected on the basis of their NIR colors, a criterion which does not match with our MIR-based approach. The classification of these sources derived by Bowler et al. (2009) also match with the classes derived in this paper. Although Bowler et al. (2009) report a significant 8 m excess for B33-31 and classify it as a bona-fide YSO, we do not detect the source at all at 8.0 m. We do retain it as a candidate YSO. Since we have also analyzed MIPS data for the region, we do identify a larger number of candidate YSOs than seen by Bowler et al. (2009).

The present Spitzer-based study has enhanced the number of YSOs in the young star cluster associated with NGC 2023, which has HD 37903 as its most massive member. Depoy et al. (1990) found 16 young stars in this cluster based on deep near-IR imaging. We find that two of them are heavily reddened foreground stars, and the true number is therefore 14. In this study we find 32 or possibly as few as 22 YSOs in the star cluster associated with NGC 2023. Even with our larger estimate of the number of YSOs in the NGC 2023 cluster, we are not even close to the number of stars, 50, Depoy et al. (1990) had estimated to be surrounding HD 37903. Their estimate was based on the assumption that even a small cluster follows a Salpeter or Miller & Scalo IMF, interstellar mass function (Miller & Scalo 1979). Several investigations, especially from groups that study clusters in star forming galaxies, show that there can be substantial stochastic fluctuations when one deals with small number statistics (see e.g., Cerviño & Valls-Gabaud 2003; Carigi & Hernandez 2008). However, we do know that we have missed PMS stars at the faint end, i.e. very low mass stars. If we add stars, which were not seen by 2MASS, but which are very red (0.3) based on their [3.6]-[4.5] colors, we find another 8 stars south of the reflection nebula, one to the north west of HD 37903 inside the reflection nebula, and another 7 in the sub-millimeter ridge. Most of these are likely to be low mass PMS stars. The true, but still conservative number for the NGC 2023 cluster is therefore 32 stars. If we account for stochastic fluctuations and also take into consideration that this a very young cluster, where star formation is still going on, there is nothing anomalous with the mass distribution of these young stars. Some of the low-mass stars may still be in the process of forming or be deeply embedded in the molecular cloud core. They may therefore not even be detectable by deep mid-IR observations, like the IRAC observations that we have employed in this study.

10.2 The sub-millimeter and protostellar candidates

SCUBA maps at 450 and 850 m reliably detect only 5 sub-millimeter sources, MM 1-5. Of these sub-millimeter cores, MM 1 has been shown to be a Class 0 object by Sandell et al. (1999). Based on isothermal grey-body modelling we find that MM 1 has a rather low value of . Since MM 1 is known to drive an energetic outflow, it must have an actively accreting disk. It therefore seems that the inferred low value is unlikely to be due to the formation of large grains. We argue that the low  arises due to the inner disk being partially optically thick. We have imaged the protostellar disk associated with MM 1 using interferometric millimeter observations of the NH₂D emission. However based on the observed velocity structure we conclude that the NH₂D emission appears to trace rotation in the surrounding envelope rather than the protostellar disk. The MM 1/MM 2 region coincides with an IR dark cloud and the presence of two faint MIPS sources suggests that stars have already formed in the cloud core. MM 2 and the surrounding cloud core is extremely dense and cold and is not evident in single dish C¹⁸O and C¹⁷O, while NH₂D shows a very similar morphology to that of the dust continuum, especially for the cold, dense, dark cloud core surrounding MM 1 and MM 2. Analysis of column densities derived from C¹⁸O and C¹⁷O compared to column densities estimated from the dust emission suggest that CO is strongly depleted in the cold core, presumably because a large fraction of the CO is frozen onto dust grains. NH₂D, on the other hand, appears to be strongly enhanced in the cold cloud core. The mass of the gas of MM 2 and the cloud core surrounding MM 1 and MM 2, 25  , is certainly sufficient to form additional low-mass stars in this cloud core.

The sub-millimeter sources MM 3 and MM 4 lie next to the NGC 2023 reflection nebula, which is where we see the highest concentration of young stars. Here the surrounding molecular cloud is much warmer and we detect no CO depletion. MM 3 is a compact sub-millimeter source, which coincides with a bright mid-infrared source, L1630MIR-63, which has a very strong infrared excess and which we classify as an extreme Class I source based on color-color diagrams. In the millimeter and sub-millimeter is looks almost identical to MM 1, although it is not clear whether it drives an outflow and it does not have any observable free-free emission. The dust emissivity index, , deduced for MM 3 is almost identical to MM 1 and although the disk/envelope mass is only about a third of the disk around MM 1, both stars have very massive disk. Yet MM 3 is detected even in the near-IR, while MM 1 is not even seen at 24 m, suggesting that MM 3 is a more evolved star than MM 1. It is, however, possible that the main reason why we do not detect MM 1 in the mid-IR is because we see the disk nearly edge on. If the disk surrounding MM 3 is more edge on, and therefore not blocking the light from the central protostar, both could be in a very similar stage of their evolution. MM 4 is the most extended dust core we identified in our sample. It differs from typical pre-stellar cores, because it is relatively warm and it also has a slightly anomalous -index. MM 5 is a cold prestellar core with both temperature and characteristic of most pre-stellar cores.

As discussed briefly in Sects. 5 and 7 we find substantially different masses for almost all of the sub-millimeter cores that overlap with the sample of Johnstone et al. (2006) and Nutter & Ward-Thompson (2007), except MM 5, (SMM J054149-02212 (Johnstone et al. 2006), and OrionBS-541299-22114 (Nutter & Ward-Thompson 2007)). For this sub-millimeter core all three studies find about the same mass, 3-4  . Here we review why there are such large difference in masses estimated for the cores MM 1-4.

Johnstone et al. (2006) and Nutter & Ward-Thompson (2007) estimate very similar masses for the MM 1/MM 2 region, 15   (SMM J054141-02180 Johnstone et al. 2006) and 13   (OrionBS-541257-21809 Nutter & Ward-Thompson 2007), while we show that this ``core'' breaks up into two sub-millimeter cores embedded in a very dense, cold cloud. MM 1 coincides with a Class 0 protostar, while MM 2 appears to be a pre-stellar core (Sect. 5). They differ in temperature and have very different SEDs. MM 2 has a , while MM 1 has a very flat SED with . If we add the masses for MM 1 and MM 2 we get 8  , only about half of the mass found by treating MM 1 and MM 2 as a single core. For MM 3, which we also identify as a Class 0/I protostar, with an SED similar to MM 1, the differences in mass are large. We find 0.3  , while Johnstone et al. (2006) determine a mass for of 16   (core SMM J054162-02172), and Nutter & Ward-Thompson (2007) find 8   (core OrionBS-541375-21733). Here the large difference is in part due to MM 3 having a very flat SED, but also because both groups deduce a much larger size for MM 3 than we do. We find it point-like and unresolved (Table 1). For MM 4 Johnstone et al. (2006) (SMM J054167-02171) and Nutter & Ward-Thompson (2007) (OrionBS-541404-21702) obtained masses of 8 and 11.0 , respectively, which is about 10 times higher than what we find. In this case the discrepancy primarily arises due to the differences in the derived source size and assuming normal dust properties, while we find that the dust emissivity differs from that of a typical pre-stellar core (Sect. 7.3).

We attribute the primary source of the different mass estimates to the manner in which a core is defined. Automatic clump finding algorithms like clfind appear to find large cores, if the core is embedded in extended emission. Furthermore the previous studies which are both at a single wavelength also identify ``clumps'', which appear spurious, like some of the cores in the molecular cloud south of NGC 2023. A protostellar core should have a correspondence in some dense gas tracers, like HCN, HCO⁺, CS, etc. and stand out with a single velocity. It can have a gradient due to collapse, rotation or expansion. All the cores do not necessarily have the same temperature or dust properties. Cores with embedded protostars have rather different dust properties (-index) and temperature structure. Resolution of the apparent disagreement between different sets of observations will be possible only by using far-infrared observations having similar spatial resolution, so that it is possible to determine the dust emissivity and temperatures of these cores reliably. This will soon be possible with missions like Herschel and SOFIA, both of which have far-infrared imagers providing spatial resolution at a 100 m of 8-10 , very similar to what can be achieved on large ground based single dish telescopes. Only then can we reliably assess how closely the pre-stellar core mass function resembles the interstellar mass function.

11 Summary

We have investigated the distribution of gas, dust and young stars in the vicinity of NGC 2023 and the L 1630 molecular cloud using a combination of the best-available datasets at multiple wavelengths.

We reliably identify 5 sub-millimeter cores MM 1-5 at 450 and 850 m using SCUBA observations. Of these cores MM 1 and MM 3 are bona fide Class I/0 sources. Single temperature greybody models have been fitted to the SEDs of the sub-millimeter cores to derive dust temperatures and emissivity indices (). We find that adjacent cores MM 1 and MM 2 show vastly different values of 0.7 and 2.1 respectively. Large differences in values are obtained for the sources MM 3 (0.8) and MM 4 (1.5) as well. The lower (<1) values of compared to those found in PMS objects could be explained by the presence of partially optically thick disks and/or grain growth in disks. In either case it is a clear signature of the presence of protostellar disks. We find that the values of correlate with the dust temperatures and hence the evolutionary stage of the source. The derived significant variation of values highlights a potential problem in the commonly made assumption of uniform values of when calculating masses of cores distributed over a region.

Using NIR (2MASS) and MIR (Spitzer IRAC and MIPS) photometry we have generated color-color diagrams to estimate the ages of the detected sources. We identify a total of 73 PMS objects in the entire 22 26region, out of which a possible 22-32 are associated with the very young cluster associated with NGC 2023 of which HD 37903 is the most massive member. We categorize 5 YSOs as Class 0/I objects and 10 sources as Class I/II objects. For the 10 sources which are detected only in the MIPS 24 m band the protostellar nature needs to be confirmed. Except for MM 3 none of the young MIR sources has been detected in the sub-millimeter. This, together with the fact that the available low (compared to the present datasets) resolution FIR maps look significantly different from the data presented here, emphasizes the need for follow-up FIR observations with comparable resolution, in order to conclusively determine the evolutionary stages of these objects. The MIR datasets also reveal several sources like HD37903, V615 Ori, Sellgren's star C etc., which are bright in the optical. The UV datasets have been used to refine the spectral classification of HD37903 and the MIR datasets further suggest that it as a Herbig Be star with significant infrared excess. For V615 Ori, which most likely drives a weak extended CO outflow, the infrared excess can be explained by a low-mass (0.06  ) circumstellar disk and envelope, which is too faint to be detected in the sub-millimeter.

We have used high resolution NH₂D observations to probe the protostellar disk associated with the Class 0 source MM 1. Although we detect enhanced NH₂D emission, the velocities suggest that the emission arises primarily from the surrounding ambient cloud which is also rotating. Based on C¹⁸O observations we confirm CO depletion in the MM 1/MM 2 region, consistent with the enhanced NH₂D emission. In contrast MM 3, also identified to be a Class 0/I star and estimated to have a dust temperature, somewhat higher than the temperature of MM 1/MM 2 shows no signs of CO depletion.

Based on the derived distribution of YSOs in NGC 2023/L 1630 region we conclude that the expansion of the H II region IC 434 has triggered the formation of stars in the compressed gas inside the L 1630 molecular cloud, as well as in some of the dense ``elephant trunk'' structures. We have been able to derive reliable age estimates of most of the YSOs in the mapped region based on NIR and MIR photometry. However far-infrared observations at comparable angular resolutions are absolutely essential to determine the evolutionary stages of the extremely young objects which are clearly detected only at 24 m and also to confirm the protostellar nature of three of the sub-millimeter clumps (MM 2, MM 4 and MM 5) for which no other signatures are so far available.

Acknowledgements

We thank Tim Jenness for obtaining some of the SCUBA data discussed in this paper. William Vacca helped us with spectral classification and editing. Brian Fleming provided us with a calibrated Spitzer IRS spectrum of MM 3. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency. This paper has made extensive use of the SIMBAD database, operated at CDS, Strasbourg, France. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

References

• Abergel, A., Bernard, J. P., Boulanger, F., et al. 2002, A&A, 389, 239 [NASA ADS] [CrossRef] [EDP Sciences]
• Allen, L. E., Calvet, N., D'Alessio, P., et al. 2004, ApJS, 154, 363 [NASA ADS] [CrossRef]
• André, P., Deeney, B. D., Phillips, R. B., et al. 1992, ApJ, 401, 667 [NASA ADS] [CrossRef]
• Andrews, S. M., & Williams, J. P. 2005, ApJ, 631, 1134 [NASA ADS] [CrossRef]
• Andrews, S. M., & Williams, J. P. 2007, ApJ, 671, 1800 [NASA ADS] [CrossRef]
• Bacmann, A., Lefloch, B., Ceccarelli, C., et al. 2002, A&A, 389, L6 [NASA ADS] [CrossRef] [EDP Sciences]
• Beckwith, S. V. W., Sargent, A. I., Chini, R. S., et al. 1990, AJ, 99, 924 [NASA ADS] [CrossRef]
• Beckwith, S. V. W., Henning, T., & Nakagawa, Y. 2000, Protostars and Planets IV, 533, ed. V. Mannings, A. P. Boss, & S. S. Russell
• Birnstiel, T., Dullemond, C. P., & Brauer, F. 2009, A&A, 503, L5 [NASA ADS] [CrossRef] [EDP Sciences]
• Bowler, B. P., Waller, W. H., Megeath, S. T., Patten, B. M., & Tamura, M. 2009, AJ, 137, 3685 [NASA ADS] [CrossRef]
• Brown, A. G. A., de Geus, E. J., & de Zeeuw, P. T. 1994, A&A, 289, 101 [NASA ADS]
• Carigi, L., & Hernandez, X. 2008, MNRAS, 390. 582
• Cerviño, M., & Valls-Gabaud, D. 2003, MNRAS, 338, 481 [NASA ADS] [CrossRef]
• Caballero, J. A. 1998, MNRAS, 383,750
• Chandler, C., Koerner, D. W., Sargent, A. I., et al. 1995, ApJ, 449, L139 [NASA ADS] [CrossRef]
• Crawford, D. L. 1958, ApJ, 128, 185 [NASA ADS] [CrossRef]
• Dent, W. R. F., Matthews, H. E., & Ward-Thompson, D. 1998, MNRAS, 301, 1049 [NASA ADS] [CrossRef]
• Depoy, D. L., Lada, E. A., Gatley, I., et al. 1990, ApJ, 356, L55 [NASA ADS] [CrossRef]
• Droege, T. F., Richmond, M. W., & Sallman, M. 2007, VizieR Online Data Catalog, 2271, 0
• Dullemond, C. P., & Dominik, C. 2005, A&A, 434, 971 [NASA ADS] [CrossRef] [EDP Sciences]
• Flaherty, K. M., Pipher, J. L., Megeath, S. T., et al. 2007, ApJ, 663, 1069 [NASA ADS] [CrossRef]
• Froebrich, D. 2005, ApJS, 156, 169 [NASA ADS] [CrossRef]
• Fuente, A., Martín-Pintado, J., Bachiller, R., et al. 1998, A&A, 334, 253 [NASA ADS]
• Fuente, A., Neri, R., Martín-Pintado, J., Bachiller, R., & Palla, F. 2001, A&A, 336, 873 [NASA ADS] [CrossRef]
• Hartmann, L., Megeath, S. T., Allen, L., et al. 2005, ApJ, 629, 881 [NASA ADS] [CrossRef]
• Harvey, P. M., Thronson, H. A., Jr., & Gatley, I. 1980, ApJ, 235, 894 [NASA ADS] [CrossRef]
• Hernández, J., Calvet, N., Hartmann, L., et al. 2005, AJ, 129, 856 [NASA ADS] [CrossRef]
• Hildebrand, R. H. 1983, Royal Astron. Soc. Quart. J., 24, 267 [NASA ADS]
• Holland, W. S., Robson, E. I., Gear, W. K., et al. 1999, MNRAS, 303, 659 [NASA ADS] [CrossRef]
• Jenness, T., & Lightfoot, J. F. 1999, Starlink User Note 216.6, Rutherford Appleton Laboratory, Particle Physics & Astronomy Research Council
• Johnstone, D., Matthews, H., & Mitchell, G. F. 2006, ApJ, 639, 259 [NASA ADS] [CrossRef]
• Kang, M., Choi, M., Ho, P. T. P., et al. 2008, ApJ, 683, 267 [NASA ADS] [CrossRef]
• Kukarkin, B. V., Kholopov, P. N., Pskvsky, Y. P., et al. 1971, in General Catalogue of Variable Stars, 3rd edn. Moskva: Astronomical Council of the Academy of Sciences in the USSR, Sternberg State Astronomical Institute of the Moscow State University
• Lada, E. A., Bally, J., & Stark, A. A. 1991a, ApJ, 368, 432 [NASA ADS] [CrossRef]
• Lada, E. A., DePoy, D. L., Evans, II, N. J., et al. 1991b, ApJ, 371,171
• Launhardt, R., Mezger, P. G., Haslam, C. G. T., et al. 1996, A&A, 312, 569 [NASA ADS]
• Malin, D. F., Ogura, K., & Walsh, J. R. 1987, MNRAS, 227, 361 [NASA ADS]
• Megeath, S. T., Allen, L. E., Gutermuth, R. A., et al. 2004, ApJS, 154, 367 [NASA ADS] [CrossRef]
• Miller, E., & Scalo, J. M. 1979, ApJS, 41, 513 [NASA ADS] [CrossRef]
• Miyake, K., & Nakagawa, Y. 1993, Icarus, 106, 20 [NASA ADS] [CrossRef]
• Mookerjea, B., Ghosh, S. K., Rengarajan, T. N., et al. 2000, AJ, 120, 1954 [NASA ADS] [CrossRef]
• Mookerjea, B., Sandell, G., Stutzki, J., et al. 2007, A&A, 473, 485 [NASA ADS] [CrossRef] [EDP Sciences]
• Muzerolle, J., Megeath, S. T., Gutermuth, R. A., et al. 2004, ApJS, 154, 379 [NASA ADS] [CrossRef]
• Nutter, D., & Ward-Thompson, D. 2007, MNRAS, 374, 413 [NASA ADS] [CrossRef]
• Pound, M. W., Reipurth, B., & Bally, J. 2003, AJ, 125, 2108 [NASA ADS] [CrossRef]
• Reipurth, B., Rodríguez, L. F., Anglada, L., et al. 2004, AJ, 127, 1736 [NASA ADS] [CrossRef]
• Robitaille, T. P., Whitney, B. A., Indebetouw, R., et al. 2007, ApJS, 169, 328 [NASA ADS] [CrossRef]
• Roueff, E., Lis, D. C., van der Tak, F. F. S., Gerin, M., & Goldsmith, P. F. 2005, A&A, 438, 585 [NASA ADS] [CrossRef] [EDP Sciences]
• Rountree, J., & Sonneborn, G. 1991, ApJ, 369, 515 [NASA ADS] [CrossRef]
• Sandell, G., Avery, L. W., Baas, F., et al. 1999, ApJ, 519, 236 [NASA ADS] [CrossRef]
• Sandell, G., Jessop, N., & Jenness, T. 2001, The SCUBA map reduction Cookbook, Starlink Cookbook 11.2 ,Rutherford Appleton Laboratory, Particle Physics & Astronomy Research Council
• Sault, R. J., Teuben, P. J., & Wright, M. C. H. 1995, in Astronomical Data Analysis Software and Systems IV, ed. R. A. Shaw, H. E. Payne, & J. J. E. Hayes (ASP: San Francisco), 433
• Sellgren, K. 1983, AJ, 88, 985 [NASA ADS] [CrossRef]
• Simon, J. D., Bolatto, A. D., Whitney, B. A., et al. 2007, ApJ, 669, 327 [NASA ADS] [CrossRef]
• Stebbins, J., Huffer, C. M., & Whitford, A. E. 1940, ApJ, 91, 20 [NASA ADS] [CrossRef]
• Stephenson, C. B. 1986, ApJ, 300, 779 [NASA ADS] [CrossRef]
• Thé, P. S., de Winter, D., & Pérez, M. R. 1994, A&AS, 104, 315 [NASA ADS]
• Tiné, S., Roueff, E., Falgarone, E., Gerin, M., & Pineau des Forêts, G. 2000, A&A, 356, 1039 [NASA ADS]
• Warren, W. H., Jr, & Hesser, J. E. 1977, ApJS, 34, 115 [NASA ADS] [CrossRef]
• Weaver, Wm. B., & Babcock, A. 2004, PASP, 116, 1035 [NASA ADS] [CrossRef]
• Witt, A. N., Schild, R. E., & Kraiman, J. B. 1984, ApJ, 281, 708 [NASA ADS] [CrossRef]
• Wyrowski, F., Walmsley, C. M., Goss, W. M., et al. 2000, ApJ, 543, 245 [NASA ADS] [CrossRef]

Online Material

Table 2:   Coordinates and flux densities of MIR sources detected in the vicinity of NGC 2023.

Footnotes

... NGC 2023

Table 2 is only available in electronic form at http://www.aanda.org

... BIMA

The BIMA array was operated by the Universities of California (Berkeley), Illinois, and Maryland with support from the National Science Foundation.

... JCMT

James Clerk Maxwell Telescope (JCMT) is operated by the Joint Astronomy Centre on behalf of the Science and Technology Facilities Council of the United Kingdom, the Netherlands Organisation for Scientific Research, and the National Research Council of Canada.

... archive

The JCMT Archive project is a collaboration between the canadian astronomy data centre (CADC), Victoria and the James Clerk Maxwell Telescope (JCMT), Hilo.

... GILDAS

http://www.iram.fr/IRAMFR/GILDAS

... detection

IRAS 05384-0229 and IRAS 05386-0229 in the Horsehead nebula; IRAS 05385-0222, IRAS 05384-0219, and IRAS 05385-0212 in the IC434 ionization front; IRAS 05392-0214, IRAS 05390-0214, and IRAS 05391-0217 (near HD 37903) in NGC2023 are all spurious detections, since they are not seen in any IRAC bands nor in the MIPS 24 m image.

... HBC 500

Both Sellgren's star C (Sellgren 1983) and HBC 500 have incorrect coordinate entries in the SIMBAD data base; neither corresponds to a visible star or IR source.

All Tables

Table 1:   Results for submillimeter sources mapped with SCUBA.

Table 3:   Coordinates and flux densities of 2MASS & IRAC 3.6/4.5 m PMS sources in NGC 2023/L 1630.

Table 4:   Sources detected only in the MIPS 24 m band.

Table 5:   CLASS HFS fitting of selected NH₂D spectra. MM 1 is at offset (1 2,1 8).

Table 2:   Coordinates and flux densities of MIR sources detected in the vicinity of NGC 2023.

All Figures

---

Copyright ESO 2009