<figure><img src="image_1.png" /></figure>

HAL
open science

# Study of morphology and stellar content of the Galactic HII region IRAS 16148-5011

K. K. Mallick, D. K. Ojha, M. Tamura, H. Linz, M. R. Samal, S. K. Ghosh

► To cite this version:

K. K. Mallick, D. K. Ojha, M. Tamura, H. Linz, M. R. Samal, et al. Study of morphology and stellar content of the Galactic HII region IRAS 16148-5011. Monthly Notices of the Royal Astronomical Society, 2015, 447 (3), pp.2307-2321. 10.1093/mnras/stu2584. hal-01439999

HAL Id: hal-01439999

https://hal.science/hal-01439999

Submitted on 15 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Study of morphology and stellar content of the Galactic H II region
IRAS 16148–5011

K. K. Mallick,¹★ D. K. Ojha,¹ M. Tamura,² H. Linz,³ M. R. Samal⁴
and S. K. Ghosh¹,⁵

¹Department of Astronomy and Astrophysics, Tata Institute of Fundamental Research, Homi Bhabha Road, Colaba, Mumbai 400 005, India

²National Astronomical Observatory of Japan, Mitaka, Tokyo 181-8588, Japan

³Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany

⁴Aix Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, F-13388 Marseille, France

⁵National Centre for Radio Astrophysics, Tata Institute of Fundamental Research, Pune 411 007, India

Accepted 2014 December 4. Received 2014 December 4; in original form 2014 September 18

ABSTRACT

An investigation of the IRAS 16148-5011 region – a cluster at a distance of 3.6 kpc – is presented here, carried out using multiwavelength data in near-infrared (NIR) from the 1.4 m Infrared Survey Facility telescope, mid-infrared (MIR) from the archival *Spitzer* GLIMPSE (Galactic Legacy Infrared Midplane Survey Extraordinaire) survey, far-infrared (FIR) from the *Herschel* archive, and low-frequency radio continuum observations at 1280 and 843 MHz from the Giant Metrewave Radio Telescope and Molonglo Survey archive, respectively. A combination of NIR and MIR data is used to identify 7 Class I and 133 Class II sources in the region. Spectral energy distribution (SED) analysis of selected sources reveals a 9.6 $M_{\odot}$ high-mass source embedded in nebulosity. However, Lyman continuum luminosity calculation using radio emission – which shows a compact H II region – indicates the spectral type of the ionizing source to be earlier than B0-O9.5. Free-free emission SED modelling yields the electron density as 138 $cm^{-3}$, and thus the mass of the ionized hydrogen as ~16.4 $M_{\odot}$. Thermal dust emission modelling, using the FIR data from *Herschel* and performing modified blackbody fits, helped us construct the temperature and column density maps of the region, which show peak values of 30 K and 3.3 $\times$ 10$^{22}~cm^{-2}$, respectively. The column density maps reveal an $A_V > 20$ mag extinction associated with the nebular emission, and weak filamentary structures connecting dense clumps. The clump associated with this *IRAS* object is found to have dimensions of ~1.1 pc $\times$ 0.8 pc, and a mass of 1023 $M_{\odot}$.

**Key words:** stars: formation–dust, extinction–H II regions–ISM: individual objects: IRAS 16148–5011–infrared: ISM–radio continuum: ISM.

# 1 INTRODUCTION

Most star formation activity is known to take place in clusters (Lada & Lada 2003), and as such, observational studies of young embedded stellar cluster regions are imperative, because they serve as a template to further investigate various associated processes and their signatures. Due to the youth of such regions and the fact that their natal medium has still not been dispersed, these clusters can be used to scrutinize various theories related to star formation, stellar cluster dynamics, as well as stellar and cloud evolution. Of even more importance are the cluster regions which harbour high-mass stars, partly because such regions are few and far between, and

more so as high-mass star formation is not very well understood (Zinnecker & Yorke 2007).

IRAS 16148–5011 is an infrared nebula (Fig. 1) in the south-
ern sky ($\alpha_{2000} = 16^h 18^m 35.2$, $\delta_{2000} = -50^\circ 18'53''$) associated with
which is an IR cluster found using the Two Micron All-Sky Survey
(2MASS) data by Dutra et al. (2003). It is located at the Galactic
plane ($l \sim 333.047, b \sim +0.037$) and is in the vicinity of the well-
known star-forming region RCW 106. Though other star-forming
regions are present nearby, IRAS 16148–5011 appeared to be a
relatively isolated region in past mappings (Karnik et al. 2001;
Mookerjea et al. 2004). Kinematic distance estimates to this re-
gion vary from ~3.3 to 11.9 kpc (near- and far-distance estimates;
Molinari et al. 2008), and we adopt the distance of 3.6 kpc from
Lumsden et al. (2013, based on the spectrophotometric distance of a
prominent source – ‘G333.0494+00.0324B’ in their nomenclature

<footnote>* E-mail: kshitiz@tifr.res.in</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Colour composite image of the IRAS 16148–5011 region, made using the *K*<sub>s</sub> (red), *H* (green), and *J* (blue) band images. The area of analysis in this paper has been marked by the dashed rectangle. Plus symbol marks the *IRAS* catalogue position of IRAS 16148–5011. The scale bar shows 1 pc extent at a distance of 3.6 kpc.</figcaption></figure>

– associated with the central nebula<fnref n="1" />) for our work. The compilation of Lumsden et al. (2013) also reveals other nearby regions at this distance, providing preliminary indications that this could be a part of larger complex. In the early analyses of Haynes, Caswell & Simons (1979, also see Chan, Henning & Schreyer 1996), radio continuum emission at 6 cm was detected, with peak in the neighbourhood (≥1 arcmin, positional accuracy ~30 arcsec) of this region, and it was predicted to be harbouring massive young stellar objects (YSOs). The *IRAS* colour analysis by MacLeod et al. (1998) was also found to be consistent with that for an H II region, and Molinari et al. (2008) – using the *IRAS* [25–12]’ colour value – have suggested that this region is among the younger *IRAS*-detected regions. A high-mass stellar source was detected by Grave & Kumar (2009) with the help of spectral energy distribution (SED) fitting using near-infrared (NIR) to millimetre (mm) data. Analyses of this region at 1.2-mm dust continuum emission and at molecular lines have revealed the presence of dense gas, with a large column density, as well as massive clumps. The total luminosity estimate (~4.4 × 10⁴ L<sub>☉</sub>, by integrating *IRAS* flux densities) has also been found to be well in the regime of high-mass stellar objects (Fontani et al. 2005; Beltrán et al. 2006). Therefore, taking into account these characteristics, this region makes a good candidate to carry out an investigation for understanding the morphology and stellar population. However, since this is a possible H II region with an embedded cluster, multiwavelength observations are required to fully discern this region’s various constituents and how they relate to each other. In this paper, we have tried to accomplish this using deep NIR observations, archival mid-infrared (MIR) *Spitzer* data, archival far-infrared (FIR) *Herschel* data, and low-frequency radio continuum observations.

In Section 2, we detail the various observations and the corresponding data analysis procedures, followed by an examination of the stellar population in the region in Section 3. The morphology of the region is discussed in Section 4. Discussion and conclusions are presented in Sections 5 and 6, respectively.

## 2 OBSERVATIONS AND DATA REDUCTION

### 2.1 NIR observations

NIR photometric observations in *J* (1.25 µm), *H* (1.63 µm), and *K*<sub>s</sub> (2.14 µm) bands (centred on α<sub>2000</sub> ~ 16<sup>h</sup>18<sup>m</sup>31<sup>s</sup>, δ<sub>2000</sub> ~ -50°17'32'' ) were carried out on 2004 July 29 using the 1.4 m Infrared Survey Facility (IRSF) telescope, South Africa. The observations were taken with the help of the Simultaneous InfraRed Imager for Unbiased Survey (SIRIUS) instrument, a three colour simultaneous camera mounted at the f/10 Cassegrain focus of the telescope. SIRIUS is equipped with three 1k × 1k HgCdTe arrays, each of which, with a pixel scale of 0.45 arcsec, provide a field of view (FoV) of 7.8 arcmin × 7.8 arcmin. Further details can be obtained from Nagashima et al. (1999) and Nagayama et al. (2003). Five sets of frames, with each set containing observations at ten dithered positions (exposure time of 5 s at each dither position), were obtained (i.e. total exposure time = 5 × 10 × 5 s = 250 s, in each band). The sky conditions were photometric, with a seeing size of ~1.35 arcsec. In addition to the target field, a sky region (~10 arcmin to the north of target region) and the standard star P9172 (Persson et al. 1998) were also observed.

Following a standard data reduction procedure, which involved bad pixel masking, dark subtraction, flat-field correction, sky subtraction, combining dithered frames, and astrometric calibration, point spread function (psf) photometry was carried out using the ALLSTAR algorithm of the *DAOPHOT* package in IRAF. About 11–13 sources were used to construct the psf for each band. Finally, the instrumental magnitudes were calibrated using the standard star P9172. The astrometric calibration rms obtained was <0.05 arcsec, and the median photometric error <0.05 mag. On comparing our catalogues with the 2MASS catalogues, sources with *K*<sub>s</sub> magnitude ≤ 9.5 were found to be saturated, and hence had their *J*, *H*, and *K*<sub>s</sub> magnitudes replaced by the corresponding 2MASS magnitudes. The final NIR catalogue is used to identify the YSOs in the region. The area of study in this paper is about 5.5 arcmin × 5.5 arcmin encompassing the nebular cloud region, centred on α<sub>2000</sub> ~ 16<sup>h</sup>18<sup>m</sup>35<sup>s</sup>, δ<sub>2000</sub> ~ -50°19'18''.

Completeness limits were calculated for all three bands by carrying out artificial star experiments using the *ADDSTAR* package in IRAF. A fixed number of stars were added in each 0.5 mag bin and analysis carried out to see how many stars are detected. The ratio of the number of detected stars to the number of added stars gives us the completeness fraction as a function of magnitude. The 90 per cent completeness limit was thus calculated as 16.6, 15.8, and 15.5 for the *J*, *H*, and *K*<sub>s</sub> bands, respectively.

### 2.2 Radio continuum observations

Radio continuum observations at 1280 MHz were obtained on 2012 November 09 using the Giant Metrewave Radio Telescope (GMRT) array. The GMRT array consists of 30 antennas arranged in an approximate Y-shaped configuration, with each antenna having a diameter of 45 m. This translates to a primary beam-size of 26.2 arcmin at 1280 MHz. A central region of ~1 km × 1 km contains 12 randomly distributed antennas, while the remaining 18 are

<footnote n="1">¹ see http://rms.leeds.ac.uk/.</footnote>
along the three radial arms (six along each arm) which extend up to ~14 km. Details about the GMRT can be found in Swarup et al. (1991).

For our observations, the Very Large Array phase and flux calibrators ‘1626-298’ and ‘3C286’, respectively, were used. The total observation time (including the calibrators) was about 3.5 h, limited by the low declination of the source. Data reduction was carried out using the AIPS software. Initial steps involved flagging the bad data (carried out using a combination of ‘VPLOT-UVFLG’ and ‘TVFLG’ tasks) and calibration (carried out using ‘CALIB-GETJY-CLCAL’). After a few iterations of flagging and calibration, the source data were ‘SPLIT’ from the whole, and were used for imaging using the task ‘IMAGR’. A few rounds of (phase) self-calibration were also carried out using the task ‘CALIB’ to remove any ionospheric phase distortion effects. To check that the flux calibration was done correctly, the image of the flux calibrator was constructed, and its flux determined and checked against literature values. The final (target) source images were rescaled to take into account the system temperature corrections for the GMRT. These corrections are required as, at the Galactic plane, the large amount of radiation at metre wavelengths increases the effective temperature of the antennas. This was done as follows. Using the sky temperature map of Haslam et al. (1982) at 408 MHz, and the spectral index of -2.6 given therein, the temperature towards this region was obtained, i.e. $T_{\text{frequency}} = T_{408} \times (\text{frequency}/408 \text{ MHz})^{-2.6}$, where $T_{408}$ is the temperature at 408 MHz and $T_{\text{frequency}}$ is the temperature at required frequency (1280 MHz here). The images were subsequently rescaled by a factor given by $(T_{\text{frequency}} + T_{\text{sys}})/T_{\text{sys}}$, i.e. the ratio of temperature towards the target to that towards the flux calibrator. $T_{\text{sys}}$, the system temperature, was obtained from the GMRT manual.<fnref n="2" />

In addition to the GMRT observations, we also obtained the archival first epoch Molongolo Galactic Plane Survey (MGPS) data at 843 MHz<fnref n="3" /> (Green et al. 1999). The MGPS observations were carried out using the Molongolo Observatory Synthesis Telescope (MOST)<fnref n="4" /> with a resolution of 43 arcsec × 43 arcsec cosec [Declination]. For our analysis purpose, we retrieved the original processed image for this region’s Galactic coordinates (resolution ~ 55.84 arcsec × 43 arcsec) from the website. The radio images are used to obtain the physical parameters and examine the morphology of the ionized gas in the region (see Section 4.2).

## 2.3 Other archival data sets

### 2.3.1 Spitzer MIR observations

Archival MIR observations of this region, obtained using the Spitzer Space Telescope under the Galactic Legacy Infrared Midplane Survey Extraordinaire (GLIMPSE) programme (Benjamin et al. 2003; Churchwell et al. 2009), were retrieved with the help of the InfraRed Science Archive<fnref n="5" /> (the Spring ’07 Archive more complete catalogue). GLIMPSE observations were taken using the InfraRed

Array Camera (IRAC) in 3.6, 4.5, 5.8, and 8.0 µm bands. The final image cutouts (pixel scale of 0.6 arcsec) as well as the final catalogues were downloaded. While the images are further used to analyse features in the region, the photometric catalogue, in conjunction with the NIR catalogue (Section 2.1) was used to identify the YSOs in the region.

### 2.3.2 Herschel FIR observations

This region has been observed at FIR wavelengths, in the range 70–500 µm, using the instruments Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) and Spectral and Photometric Imaging Receiver (SPIRE; Griffin et al. 2010) on the 3.5-m Herschel Space Observatory (Pilbratt et al. 2010), as a part of the Proposal ID ‘KPOT_smolinar_1’ (Molinari et al. 2010). For our analyses, we obtained the PACS 70- and 160-µm level2_5 MADmap images (Cantalupo et al. 2010); and the SPIRE 250-, 350-, and 500-µm level2_5 extended source (‘extdPxW’) map products using the Herschel Science Archive.<fnref n="6" /> The 70, 160, 250, 350, and 500 µm images have pixel scales of 3.2, 6.4, 6, 10, and 14 arcsec, respectively, with resolutions varying from ~5.5–36 arcsec. While the PACS image obtained had the surface brightness unit of Jy pixel⁻¹, the SPIRE images were in the units of MJy sr⁻¹. Detailed information about the data products is provided on the Herschel site.<fnref n="7" /> We use the FIR data to examine the physical conditions of the region.

# 3 STELLAR POPULATION IN THE REGION

## 3.1 Identification of YSOs

The NIR and MIR photometric catalogues from IRSF (Section 2.1) and GLIMPSE (Section 2.3.1), respectively, were cross-matched within 0.6 arcsec matching radius and collated to obtain a combined photometric catalogue. Thereafter, the following set of steps were followed for the identification of the YSOs (similar to Mallick et al. 2013).

(i) First, the YSOs were identified using their MIR magnitudes. The sources with detections in all four IRAC bands with errors ≤ 0.15 mag were used here. Using simple linear regression, the IRAC spectral index ($\alpha_{IRAC} = \log(\lambda F_{\lambda})/\log(\lambda)$; Lada 1987) was calculated for each source (in the wavelength range 3.6–8.0 µm), followed by their classification into Class I and Class II categories using the limits from Chavarría et al. (2008, $\alpha_{IRAC} > 0$ for Class I, and $-2 \le \alpha_{IRAC} \le 0$ for Class II; see Fig. 2a).

(ii) All sources need not have detections at 5.8 and 8.0 µm, but might have good quality detections in NIR bands. To identify the YSOs from amongst such sources, we use a combination of H, $K_s$, 3.6 µm, and 4.5 µm bands following the procedure of Gutermuth et al. (2009). Again, only those sources whose 3.6- and 4.5-µm magnitude errors are $\le 0.15$ are used here. The YSOs were identified from their location in the dereddened – using the colour excess ratios from Flaherty et al. (2007) – $K_s$-[3.6]' versus '[3.6]-[4.5]' colour-colour diagram (CCD), as is shown in Fig. 2(b).

(iii) Additional YSOs were identified using the $J - H/H - K$ CCD (Fig. 2c) according to the following procedure. In Fig. 2(c), the red solid curve marks the dwarf locus from Bessell & Brett (1988) and the blue solid line marks the Classical T Tauri Stars

<footnote n="2">² http://gmrt.ncra.tifr.res.in/gmrt_hpage/Users/Help/help.html</footnote>

<footnote n="3">³ http://www.physics.usyd.edu.au/sifa/Main/MGPS1</footnote>

<footnote n="4">⁴ The MOST is operated by the University of Sydney with support from the Australian Research Council and the Science Foundation for Physics within the University of Sydney.</footnote>

<footnote n="5">⁵ This research has made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.</footnote>

<footnote n="6">⁶ http://www.cosmos.esa.int/web/herschel/science-archive</footnote>

<footnote n="7">⁷ http://www.cosmos.esa.int/web/herschel/data-products-overview</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** (a) Histogram of the IRAC spectral indices. The limits for Class I and Class II sources have been marked by dashed grey vertical lines. (b) NIR-MIR CCD using the procedure of Gutermuth et al. (2009). The regions of Class I and Class II sources have been labelled. The identified YSOs have been shown with black solid circles. (c) NIR CCD in CIT system. The red curve and the blue line show the dwarf locus (Bessell & Brett 1988) and the CTTS locus (Meyer et al. 1997), respectively. The grey dot-dashed line is the extension of CTTS locus. Three parallel slanted dashed lines mark the reddening vectors, drawn using extinction laws from Cohen et al. (1981). Three separate regions, ‘F’, ‘T’, ‘P’ have been labelled on the plot. The source marked with a cross in ‘P’ region is a high-mass YSO (see text).</figcaption></figure>

(CTTS) locus from Meyer, Calvet & Hillenbrand (1997). All the loci curves as well as sources' colours were converted to CIT (California Institute of Technology) photometric system (using Carpenter 2001) for this analysis. The slanted dashed lines are the reddening vectors, drawn using the reddening laws of Cohen et al. (1981) for the CIT photometric system. In this CCD, three separate regions have been marked, similar to Ojha et al. (2004a,b). The sources in 'T' and 'P' regions are taken to be Class II sources (Lada & Adams 1992), as they exhibit IR-excess emission. In the 'P' region, since there could be a slight overlap between Herbig Ae/Be stars and Class II sources (Hillenbrand et al. 1992), we conservatively took only those sources which were above the CTTS locus extended into this region (similar to Mallick et al. 2014). It should be noted that sources in the 'T' region could also contain a few Class III sources with small IR-excess.

Finally, we merged the sources identified in each step. An overlapping source might have different identifications in different steps. Thus, in the final list, the class of a YSO was taken as that in which it was identified as first in the above order of steps. A total of 7 Class I and 133 Class II sources were obtained (for a FoV of ~5.5 arcmin × 5.5 arcmin encompassing the molecular cloud, as marked in Fig. 1) in the final YSO catalogue, which is given in Table 1.

## 3.2 SED of YSOs

SED modelling was carried out for (a subset of) YSOs to get estimates of their physical parameters. The grid of YSO models from Robitaille et al. (2006) – implemented in the online SED fitting tool of Robitaille et al. (2007) – was used for this purpose. The basic model consists of a pre-main sequence (PMS) star surrounded by
<table><caption>**Table 1.** YSOs identified using NIR and MIR data.</caption><thead><tr><th>RA (J2000)</th><th>Dec. (J2000)</th><th>J (mag)</th><th>H (mag)</th><th>K<sub>s</sub> (mag)</th><th>[3.6] (mag)</th><th>[4.5] (mag)</th><th>[5.8] (mag)</th><th>[8.0] (mag)</th><th>YSO classification</th></tr></thead><tbody><tr><td>244.575 165</td><td>-50.365 025</td><td>17.861 ± 0.130</td><td>16.631 ± 0.001</td><td>15.509 ± 0.047</td><td>14.311 ± 0.145</td><td>14.132 ± 0.203</td><td>-</td><td>-</td><td>Class2</td></tr><tr><td>244.575 241</td><td>-50.367 393</td><td>14.580 ± 0.021</td><td>13.876 ± 0.010</td><td>13.643 ± 0.032</td><td>12.983 ± 0.099</td><td>12.872 ± 0.122</td><td>-</td><td>-</td><td>Class2</td></tr><tr><td>244.576 645</td><td>-50.293 938</td><td>16.137 ± 0.016</td><td>14.443 ± 0.010</td><td>13.570 ± 0.020</td><td>12.780 ± 0.102</td><td>12.372 ± 0.132</td><td>-</td><td>-</td><td>Class2</td></tr><tr><td>244.577 545</td><td>-50.314 171</td><td>15.557 ± 0.010</td><td>14.358 ± 0.010</td><td>13.818 ± 0.028</td><td>13.001 ± 0.108</td><td>12.820 ± 0.143</td><td>-</td><td>-</td><td>Class2</td></tr><tr><td>244.579 193</td><td>-50.321 396</td><td>17.491 ± 0.045</td><td>15.921 ± 0.020</td><td>14.881 ± 0.038</td><td>13.412 ± 0.161</td><td>13.265 ± 0.132</td><td>-</td><td>-</td><td>Class2</td></tr></tbody></table>

Note. Table 1 is available in its entirety in a machine-readable form in the online journal. A portion is shown here for guidance regarding its form and content.

a flared accretion disc having a rotationally flattened envelope with cavities carved out by a bipolar outflow. A total of 200 000 SED models are computed in a 14D parameter space (covering properties of the central source, the infalling envelope, and the disc), using the radiation transfer code of Whitney et al. (2003a,b). The online fitting tool attempts to fit the available SED models to the data, characterizing each fitting by a $\chi^2$ parameter. The distance range and the interstellar visual extinction ($A_V$) are free parameters whose range has to be specified by the user. Accounting for the uncertainties, we adopted a large distance range of 3.4 to 3.8 kpc for our sources. From extinction calculations (see Section 4.1), we find that almost all non-YSO sources had $A_V \le 20$ mag, and thus we specified the range of interstellar visual extinction $A_V$ as 1–20 mag.

Since the number of SED models is very large, spanning a wide range of parameter space, the models fitting each source can only be constrained by increasing the number of data points used, and having data points to sufficiently cover the entire wavelength range of fitting. For this reason, we choose the sources with photometry in at least all four IRAC bands for SED modelling. However, in addition to those selected using this criterion, we also carried out SED analysis for the prominent sources associated with the central nebula (even though two of them lacked 8.0 µm photometry). We note that photometry from Herschel images is not useful here as none of the YSOs have counterparts at those wavelengths, which in part is due to poor resolution.

For each source, the SED fitting tool – besides giving the best-fitting model – also gives a set of well-fitted models ranked by their $\chi^2$ values as a measure of their relative ‘goodness-of-fit’. Following a method similar to Robitaille et al. (2007), we consider only those models for further calculation which satisfied the following criterion:

$$ \chi^2 - \chi_{\min}^2 < 3 \text{ (per data point).} \quad (1) $$

As elucidated in Robitaille et al. (2007), though this criterion is based on visual examination of SED plots and has no rigid mathematical backing, a stricter criterion might lead to overinterpretation. For each parameter, the weighted mean value and standard deviation were calculated using the models which satisfied equation (1). The inverse of the respective $\chi^2$ was taken as the weight for each model (similar to Grave & Kumar 2009). Table 2 gives our SED modelling results, listing the physical parameters: age of the central source ($T_*$), mass of the central source, disc mass ($M_{\text{disc}}$), disc accretion rate ($\dot{M}_{\text{disc}}$), envelope mass ($M_{\text{env}}$), temperature of the central source ($T_*$), total system luminosity ($L_{\text{total}}$), interstellar visual extinction ($A_V$), and the $\chi^2_{\min}$ per data point. Errors for some parameters tend to be large as we are dealing with a large parameter space while we have very few data points to constrain the number of models. If a source was simply fit better as a star with high interstellar extinction, it has not been included in the table. Thus finally, SED results are given for a total of three Class I sources,

24 Class II sources (including two central sources), and one extra central source.

Even though the statistics for the YSOs (i.e. Class I and Class II sources) is not very significant, we can still use them to get an idea of the physical parameters of the stellar sources forming in this region. As can be seen from Table 2, there appears to be a considerable age dispersion, with the ages ranging from ~0.05 to 0.5 Myr for most of the sources. Five YSOs even have ages >1 Myr. This is suggestive of ongoing star formation. All the YSOs analysed here yield masses >2 $M_{\odot}$, and five of them >6 $M_{\odot}$. One of the YSOs (#16) appears to be a high-mass star of ~9.6 $M_{\odot}$, and is embedded in the nebular emission associated with this region (see Section 4). The SED plot for this source is shown in Fig. 3. Grave & Kumar (2009) had also done SED analysis for what they mention as an embedded point source in this IR nebula with JHK and four IRAC bands. Since among the possible embedded sources in the central part of this IR nebula, only this source (#16) had all seven magnitudes, most likely their ‘16148-5011mms2near’ refers to this source itself. Grave & Kumar (2009) calculated the age and mass of this source as $4.2 \pm 0.3$ (log $t_*$) and $11 \pm 1$ $M_{\odot}$, respectively. Though the mass estimate appears consistent (within error limits) with our results, their age estimate is much lower. It is probable that this could be because they use different distance estimates and 1.2 mm fluxes from literature. This source is also classified as a YSO by Lumsden et al. (2013, called ‘G333.0494+00.0324B’). Another source (#28) associated with the nebula appears to be a high-mass stellar object, with age > 1 Myr, though it is not classified as a YSO in our analysis.

It should be noted that the SED results are only representative of the actual values, as an empirically consistent fit might not be the correct fit. A deliberate sampling bias of the huge parameter space, to reduce computational time, could give rise to pseudo-trends in the results. The results are also contingent upon the validity of assumed evolutionary tracks from literature. Most importantly, the models are for individual sources, and could be misleading in cases of stellar multiplicity. The caveats are dealt in detail in Robitaille (2008).

## 3.3 Luminosity function

The slope of the $K_s$-band luminosity function (KLF) can serve as an indicator of the age of a stellar cluster (Zinnecker, McCaughrean & Wilking 1993; Lada & Lada 1995; Vig et al. 2014). If we were to assume that the mass function and the mass–luminosity relation for a (coeval) stellar cluster are power laws, i.e. they are of the form $dN(\log m_*) \propto m_*^{-\gamma} d\log m_*$ and $L_K \propto m_*^B$, then it can be shown that the slope of the KLF will be of the form:

$$ \alpha = \gamma/2.5B \quad (2) $$
<table><caption>**Table 2.** Main parameters from SED analysis.</caption><thead><tr><th>S.No.</th><th>RA (J2000)</th><th>Dec. (J2000)</th><th>log T<sub>*</sub> (yr)</th><th>Mass (M<sub>☉</sub>)</th><th>log M<sub>disc</sub> (M<sub>☉</sub>)</th><th>log ˙M<sub>disc</sub> (M<sub>☉</sub> yr<sup>-1</sup>)</th><th>log M<sub>env</sub> (M<sub>☉</sub>)</th><th>log T<sub>*</sub> (K)</th><th>log L<sub>total</sub> (L<sub>☉</sub>)</th><th>A<sub>V</sub> (mag)</th><th>χ<sup>2</sup><sub>min</sub> (per data point)</th></tr></thead><tbody><tr><td colspan="12">Class I sources</td></tr><tr><td>1</td><td>244.591</td><td>080</td><td>-50.285</td><td>355</td><td>6.33 ± 0.74</td><td>2.42 ± 0.79</td><td>-3.38 ± 1.64</td><td>-9.29 ± 1.73</td><td>-4.43 ± 2.73</td><td>3.85 ± 0.15</td><td>1.41 ± 0.30</td><td>4.81 ± 2.35</td><td>0.29</td></tr><tr><td>2</td><td>244.603</td><td>455</td><td>-50.292</td><td>828</td><td>4.71 ± 0.40</td><td>4.20 ± 1.42</td><td>-1.39 ± 0.70</td><td>-6.70 ± 0.99</td><td>0.68 ± 1.06</td><td>3.66 ± 0.07</td><td>2.11 ± 0.28</td><td>15.95 ± 4.98</td><td>0.01</td></tr><tr><td>3</td><td>244.671</td><td>539</td><td>-50.338</td><td>348</td><td>5.30 ± 1.02</td><td>2.58 ± 1.61</td><td>-2.90 ± 1.82</td><td>-8.08 ± 2.02</td><td>-2.70 ± 2.90</td><td>3.73 ± 0.20</td><td>1.56 ± 0.56</td><td>7.45 ± 4.16</td><td>0.10</td></tr><tr><td colspan="12">Class II sources</td></tr><tr><td>4</td><td>244.584</td><td>152</td><td>-50.353</td><td>664</td><td>4.65 ± 0.37</td><td>6.41 ± 0.54</td><td>-2.12 ± 0.75</td><td>-6.94 ± 0.80</td><td>1.65 ± 0.40</td><td>3.66 ± 0.01</td><td>2.49 ± 0.06</td><td>12.17 ± 5.81</td><td>9.90</td></tr><tr><td>5</td><td>244.599</td><td>808</td><td>-50.286</td><td>560</td><td>4.94 ± 0.37</td><td>3.80 ± 1.81</td><td>-1.99 ± 0.84</td><td>-7.55 ± 1.30</td><td>0.54 ± 0.66</td><td>3.66 ± 0.05</td><td>1.93 ± 0.41</td><td>11.51 ± 5.27</td><td>1.93</td></tr><tr><td>6</td><td>244.606</td><td>796</td><td>-50.307</td><td>327</td><td>5.57 ± 0.87</td><td>7.50 ± 3.12</td><td>-4.67 ± 2.72</td><td>-9.18 ± 2.08</td><td>-3.10 ± 4.21</td><td>4.02 ± 0.35</td><td>2.95 ± 0.68</td><td>13.72 ± 5.07</td><td>5.34</td></tr><tr><td>7</td><td>244.609</td><td>207</td><td>-50.351</td><td>631</td><td>4.97 ± 0.27</td><td>4.68 ± 1.40</td><td>-2.11 ± 0.93</td><td>-7.66 ± 1.31</td><td>1.24 ± 0.53</td><td>3.67 ± 0.05</td><td>2.12 ± 0.32</td><td>17.37 ± 3.72</td><td>10.85</td></tr><tr><td>8</td><td>244.610</td><td>397</td><td>-50.279</td><td>411</td><td>5.14 ± 0.76</td><td>2.99 ± 1.59</td><td>-2.25 ± 0.95</td><td>-7.76 ± 1.28</td><td>-0.70 ± 2.22</td><td>3.70 ± 0.15</td><td>1.73 ± 0.36</td><td>13.43 ± 5.74</td><td>0.29</td></tr><tr><td>9</td><td>244.612</td><td>305</td><td>-50.365</td><td>067</td><td>4.84 ± 0.45</td><td>5.60 ± 2.59</td><td>-2.04 ± 1.05</td><td>-7.61 ± 1.55</td><td>1.09 ± 0.87</td><td>3.70 ± 0.07</td><td>2.31 ± 0.66</td><td>4.51 ± 2.81</td><td>0.34</td></tr><tr><td>10</td><td>244.614</td><td>349</td><td>-50.327</td><td>496</td><td>6.28 ± 0.66</td><td>2.50 ± 0.83</td><td>-3.40 ± 1.71</td><td>-9.21 ± 1.89</td><td>-4.37 ± 2.56</td><td>3.85 ± 0.17</td><td>1.44 ± 0.35</td><td>6.20 ± 2.30</td><td>0.13</td></tr><tr><td>11</td><td>244.633</td><td>514</td><td>-50.326</td><td>122</td><td>5.73 ± 0.70</td><td>2.58 ± 0.99</td><td>-3.20 ± 1.77</td><td>-8.84 ± 1.80</td><td>-2.83 ± 2.80</td><td>3.77 ± 0.19</td><td>1.54 ± 0.35</td><td>8.83 ± 4.50</td><td>0.28</td></tr><tr><td>12</td><td>244.639</td><td>954</td><td>-50.283</td><td>981</td><td>5.10 ± 0.07</td><td>4.80 ± 0.22</td><td>-2.36 ± 0.66</td><td>-7.50 ± 0.20</td><td>1.14 ± 0.56</td><td>3.66 ± 0.01</td><td>2.02 ± 0.07</td><td>9.75 ± 5.38</td><td>2.16</td></tr><tr><td>13<sup>a,b</sup></td><td>244.644</td><td>714</td><td>-50.314</td><td>804</td><td>4.18 ± 1.01</td><td>6.17 ± 1.74</td><td>-</td><td>-</td><td>-0.18 ± 1.92</td><td>3.76 ± 0.26</td><td>2.82 ± 0.31</td><td>6.18 ± 3.86</td><td>0.03</td></tr><tr><td>14</td><td>244.649</td><td>155</td><td>-50.346</td><td>909</td><td>5.65 ± 0.90</td><td>2.44 ± 1.12</td><td>-2.92 ± 1.44</td><td>-8.40 ± 1.68</td><td>-2.87 ± 3.05</td><td>3.79 ± 0.20</td><td>1.58 ± 0.35</td><td>7.63 ± 4.29</td><td>0.03</td></tr><tr><td>15</td><td>244.651</td><td>047</td><td>-50.328</td><td>751</td><td>6.29 ± 0.68</td><td>3.46 ± 0.86</td><td>-3.43 ± 1.48</td><td>-8.96 ± 1.60</td><td>-4.83 ± 2.75</td><td>3.99 ± 0.18</td><td>2.01 ± 0.27</td><td>6.70 ± 2.41</td><td>0.13</td></tr><tr><td>16<sup>a,c</sup></td><td>244.652</td><td>878</td><td>-50.316</td><td>849</td><td>5.44 ± 0.61</td><td>9.60 ± 1.44</td><td>-2.15 ± 1.24</td><td>-7.00 ± 1.35</td><td>0.00 ± 3.42</td><td>4.29 ± 0.20</td><td>3.69 ± 0.16</td><td>15.97 ± 4.36</td><td>0.89</td></tr><tr><td>17</td><td>244.663</td><td>803</td><td>-50.352</td><td>959</td><td>6.02 ± 0.78</td><td>2.93 ± 0.82</td><td>-3.23 ± 1.35</td><td>-8.92 ± 1.53</td><td>-3.75 ± 3.08</td><td>3.86 ± 0.19</td><td>1.74 ± 0.24</td><td>4.75 ± 2.79</td><td>0.07</td></tr><tr><td>18</td><td>244.671</td><td>738</td><td>-50.278</td><td>400</td><td>5.30 ± 0.28</td><td>3.02 ± 1.19</td><td>-2.97 ± 1.18</td><td>-8.96 ± 1.72</td><td>-0.74 ± 0.91</td><td>3.67 ± 0.05</td><td>1.62 ± 0.27</td><td>10.13 ± 3.37</td><td>3.23</td></tr><tr><td>19</td><td>244.679</td><td>535</td><td>-50.293</td><td>499</td><td>4.97 ± 0.36</td><td>4.10 ± 1.47</td><td>-2.07 ± 0.86</td><td>-7.70 ± 1.19</td><td>0.45 ± 0.86</td><td>3.67 ± 0.06</td><td>2.01 ± 0.30</td><td>9.06 ± 5.71</td><td>2.25</td></tr><tr><td>20</td><td>244.685</td><td>440</td><td>-50.357</td><td>010</td><td>6.29 ± 0.72</td><td>3.69 ± 0.99</td><td>-3.83 ± 1.67</td><td>-9.31 ± 1.77</td><td>-5.06 ± 2.77</td><td>4.03 ± 0.19</td><td>2.14 ± 0.32</td><td>9.15 ± 3.12</td><td>0.17</td></tr><tr><td>21</td><td>244.689</td><td>316</td><td>-50.343</td><td>143</td><td>4.99 ± 0.39</td><td>3.85 ± 1.38</td><td>-2.13 ± 0.83</td><td>-7.74 ± 1.20</td><td>0.28 ± 1.05</td><td>3.67 ± 0.06</td><td>1.94 ± 0.26</td><td>7.45 ± 5.58</td><td>2.23</td></tr><tr><td>22</td><td>244.689</td><td>453</td><td>-50.366</td><td>440</td><td>4.90 ± 0.48</td><td>4.06 ± 1.33</td><td>-1.97 ± 0.87</td><td>-7.40 ± 1.17</td><td>0.37 ± 1.39</td><td>3.67 ± 0.09</td><td>2.04 ± 0.26</td><td>13.02 ± 5.40</td><td>0.17</td></tr><tr><td>23</td><td>244.693</td><td>146</td><td>-50.279</td><td>778</td><td>5.50 ± 0.67</td><td>2.95 ± 0.96</td><td>-2.39 ± 0.94</td><td>-8.01 ± 1.08</td><td>-1.47 ± 2.53</td><td>3.73 ± 0.16</td><td>1.68 ± 0.19</td><td>7.19 ± 3.64</td><td>4.84</td></tr><tr><td>24</td><td>244.704</td><td>742</td><td>-50.351</td><td>349</td><td>4.78 ± 0.39</td><td>5.70 ± 0.86</td><td>-2.26 ± 0.30</td><td>-7.61 ± 0.94</td><td>1.65 ± 0.30</td><td>3.66 ± 0.02</td><td>2.37 ± 0.09</td><td>17.61 ± 2.63</td><td>14.26</td></tr><tr><td>25</td><td>244.708</td><td>603</td><td>-50.286</td><td>385</td><td>4.89 ± 0.35</td><td>4.59 ± 1.48</td><td>-1.92 ± 0.84</td><td>-7.52 ± 1.14</td><td>0.89 ± 0.71</td><td>3.67 ± 0.06</td><td>2.14 ± 0.34</td><td>12.96 ± 5.63</td><td>3.55</td></tr><tr><td>26</td><td>244.683</td><td>786</td><td>-50.363</td><td>702</td><td>4.75 ± 0.39</td><td>6.24 ± 2.25</td><td>-1.56 ± 0.82</td><td>-6.91 ± 1.09</td><td>1.25 ± 0.72</td><td>3.71 ± 0.11</td><td>2.52 ± 0.51</td><td>6.59 ± 5.53</td><td>1.57</td></tr><tr><td>27</td><td>244.687</td><td>016</td><td>-50.353</td><td>573</td><td>5.46 ± 0.73</td><td>2.87 ± 1.67</td><td>-2.60 ± 1.34</td><td>-8.20 ± 1.54</td><td>-1.73 ± 2.66</td><td>3.74 ± 0.16</td><td>1.62 ± 0.50</td><td>12.36 ± 6.40</td><td>0.01</td></tr><tr><td colspan="12">Extra central source</td></tr><tr><td>28<sup>a,b</sup></td><td>244.650</td><td>894</td><td>-50.314</td><td>579</td><td>6.17 ± 0.35</td><td>8.48 ± 1.17</td><td>-4.25 ± 1.30</td><td>-9.31 ± 1.32</td><td>-4.20 ± 3.65</td><td>4.36 ± 0.06</td><td>3.52 ± 0.21</td><td>18.38 ± 2.26</td><td>3.09</td></tr></tbody></table>

Notes. $^a$Prominent sources associated with the central nebula; $^b$Did not have 8.0-µm photometry; $^c$Identified as 'G333.0494+00.0324B' (and also classified as a YSO) in Lumsden et al. (2013).

<figure><img src="image_6.png" /><figcaption>**Figure 3.** SED fitting using the online tool of Robitaille et al. (2007) for source #16 (a Class II source) in Table 2. The black dots mark the data points. The solid black curve is the best-fitting model, while the grey curves denote the subsequent good fits for $\chi^2 - \chi^2_{min}$ (perdatapoint) < 3. The dashed curve is the photosphere (in the presence of interstellar extinction, but absence of circumstellar dust) of the central source for the best-fitting model.</figcaption></figure>

(Lada, Young & Greene 1993; Megeath 1996). First of all, we try to estimate the KLF slope. We only consider the YSOs here, as opposed to all the observed sources, as they will be much less affected by any field star contamination. Our K-band 100 per cent completeness limit is up to 14 mag, and thus completeness correction was implemented in the (0.5 mag sized) bins after this limit. Thereafter, the (cumulative) KLF of the YSOs was constructed, and is shown in Fig. 4. The fit to the histogram is also shown, and its slope (dlog N/dmK) in [12,15.5] mag range is calculated to be (α =) 0.35 ± 0.04 (the slope will be same for cumulative and differential KLFs here, see Lada et al. 1993). As for the slope of mass-luminosity relation, if we adopt the value of B = 2 (which can be shown to be the approximate value for O-F main-sequence stars; Lada et al. 1993), then the mass function slope comes out to be (γ =) 1.75 ± 0.20 (a value slightly steeper than the Salpeter slope of 1.35).

Alternatively, we could get an estimate of the mass function slope using the SED fitting results. Since most of our sources are of intermediate mass in 2–6 M☉ range (see Table 2), we consider only this range to estimate the slope γ. Fig. 5 shows the mass histogram (for the YSOs). Assuming that the star formation is strictly coeval, the mass function slope is given by $γ = -(dlog N/dlog m_*)$ (Massey 1998), and thus we calculate $γ \sim 1.59 \pm 0.70$ for our case. The grey curve in Fig. 5 shows the fitted function for 2–6 M☉. This is
<figure><img src="image_7.png" /><figcaption>**Figure 4.** The grey dashed histogram shows the cumulative KLF for the YSOs. The black straight line is the fit in [12,15.5] mag range, whose slope is given by $\alpha$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Histogram of stellar masses obtained from the SED fitting (Table 2). The grey curve shows the power-law fit in the intermediate mass range.</figcaption></figure>

consistent with $\gamma$ obtained above. The major source of uncertainty here is the sparse statistics, and thus possible incompleteness in the mass bins.

A similar set of values (both $\alpha$ and $\gamma$) was found by Balog et al. (2004) for the 1-Myr-old NGC 7538 cluster (distance ~ 2.65 kpc; Mallick et al. 2014). $\alpha$ also appears consistent with that derived for the Orion molecular clouds (~0.37–0.38, age ~ 1 Myr; Lada et al. 1991; Lada & Lada 1995). However, other embedded clusters in literature, such as the NGC 1893 cluster (~0.34 ± 0.07, distance ~ 3.25 kpc; Sharma et al. 2007), the Tr 14-16 clusters in Carina nebula (~0.30–0.37, distance ~ 2.5 kpc; Sanchawala et al. 2007), and the IRAS 06055+2039 cluster (~0.43 ± 0.09, distance ~ 2.6 kpc; Tej et al. 2006) do exhibit slightly larger cluster ages, up to ~4 Myr. Hence, it appears that 1 Myr should be the lower age limit of the cluster.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** *J/*J* − *H* CMD for the YSOs with at least *J*- and *H*-band detections. The dashed and dotted black curves plot the 1- and 2-Myr PMS isochrones, respectively, from Siess et al. (2000). The reddening vectors (parallel grey dashed lines) for the 1 Myr isochrone are drawn at 0.1, 1, 2, and 4 M⊙. The blue star symbols are the sources for which SED analysis was done.</figcaption></figure>

## 3.4 Mass spectrum

In addition to the SED, we also use the *J/*J* − *H* colour–magnitude diagram (CMD) to get an estimate of the mass range of the YSOs. We avoid using a CMD involving the *K*-band magnitude as this band is the most affected by the NIR excess flux arising from circumstellar material, which in turn can lead to brightening, and thus erroneous mass estimate, of the sources. Fig. 6 shows the *J/*J* − *H* CMD for the YSOs with at least *J*- and *H*-band detections. The 1-Myr PMS isochrone, along with the 2 Myr isochrone for reference, from Siess, Dufour & Forestini (2000) has been shown on the image. For the 1-Myr PMS isochrone, reddening vectors are shown for 0.1, 1, 2, and 4 M⊙. As can be seen, all but one of the YSOs for which the SED analysis has been done (marked with blue stars) lie in the mass range $\gtrsim$ 2 M⊙, matching well with our SED results. In general, the sources lie in the mass range ~0.1–4 M⊙, or, put differently, the observations probe stellar objects up to the 0.1 M⊙ limit. The (blue star) source at the far-right end of the diagram is the 9.6 M⊙ high-mass source from SED analysis. The wide variation in the colours of YSOs is probably an indication of variable extinction as well as different evolutionary stages of the sources. Major causes of uncertainty here are: uncertainty in distance estimate which is used to obtain the apparent magnitude (for the isochrones), and unresolved (especially since the source distance of ~3.6 kpc is much larger than most of the previously studied clusters) binarity. Additionally, it should be kept in mind that, in general, a particular PMS model will introduce its own systematic error.

# 4 MORPHOLOGY OF THE REGION

## 4.1 Cluster analysis and extinction mapping

Fig. 7 shows the Spitzer 8.0 µm image with overlaid surface density contours (in cyan) as well as visual extinction contours (in blue). The high-mass source in this region (see Section 3.2), as well as the mm
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Spitzer 8.0 µm image of the region with overlaid surface density contours (at 5, 5.75, 6, 7, 7.5, and 8 YSOs pc⁻²) in cyan, and visual extinction contours (at A_V = 4, 4.5, 5, 6, and 6.5 mag) in blue. Green plus symbol marks the IRAS catalogue position of IRAS 16148-5011, green cross the high-mass source (#16 from Table 2), while the diamond and box symbols mark the mm and MSX peaks, respectively, from Molinari et al. (2008).</figcaption></figure>

and MSX peaks from Molinari et al. (2008, 16148-5011 MM 1 and 16148-5011 A from their tables 2 and 3, respectively) have also been marked on the image. The marked mm peak is almost coincident with the (mm) peak from Beltrán et al. (2006, 16148-5011 Clump 2 from their table 2). The surface density and extinction contours were calculated as follows. The surface density analysis was carried out using the nearest-neighbour (NN) method (Casertano & Hut 1985) to discern the YSO clusterings in the region. We chose 20 NN, similar to Schmeja, Kumar & Ferreira (2008), Schmeja (2011), and Mallick et al. (2013). The extinction in the region was estimated by constructing the extinction map with the help of the NIR photometric data. We use the NIR CCD (see Fig. 2c) for this purpose. In this NIR CCD, the ‘F’ region mostly contains the main-sequence field sources along with a few probable Class III sources. Since these are sources which have almost lost their circumstellar material, any extinction they exhibit will come from interstellar – rather than circumstellar – material. The sources from this ‘F’ region, which were not identified as a YSO in Section 3.1, were therefore selected. Subsequently, we used a method similar to the Near-IR-Colour-Excess method of Lada et al. (1994), but here the (H-K) colour excess was estimated by dereddening the sources – along the reddening vector – to the low-mass end (turnover onwards) of the dwarf locus (which was approximated as a straight line). A_V was then calculated using the reddening laws of Cohen et al. (1981). After we obtained the visual extinction A_V for each source, an extinction map of the region was made by using the NN method, where the extinction at each grid point on the map is the median (because median rejects outliers) of A_V of 20 NN. It is possible that the extinction could be slightly underestimated because, due to the large distance to this region, a significant fraction of the detected sources could be foreground sources, especially towards the centre of the nebula.

In Fig. 7, the surface density contours are drawn at 5, 5.75, 6, 7, 7.5, and 8 YSOs pc⁻², while the extinction contours have been drawn at A_V = 4, 4.5, 5, 6, and 6.5 mag levels. As can be seen on

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Spitzer 8.0 µm image of the region with overlaid 1280-MHz contours (at 3σ, 4σ, 5σ, 7σ, 10σ, 12σ, 15σ, 20σ, and 22σ, where σ ~ 8.85 mJy; resolution ~56 arcsec × 35 arcsec) in cyan and 843 MHz contours (at 5σ, 7σ, 10σ, 15σ, 20σ, 30σ, 40σ, and 50σ, where σ ~ 5.57 mJy; resolution ~56 arcsec × 43 arcsec) in blue. The symbols are same as for Fig. 7.</figcaption></figure>

the image, both the surface density as well as extinction contours are in the southern portion of the nebular emission and appear to be along the sharp boundary of the nebula. The highest surface density contour levels coincide with the highest extinction levels. Though we would have expected to see high interstellar extinction along the main body of the nebula, this is not so possibly because the NIR observations are not deep enough to detect stars from behind the nebula and thus the extinction of the nebula should be higher than the highest extinction contour level here (it is later estimated to be A_V > 20 mag; see Section 4.4). It should also be noted that the cluster detected here is in the southern part of the nebula also coincident with extended radio continuum emission (see Section 4.2).

## 4.2 Radio morphology

Fig. 8 shows a Spitzer 8.0 µm image of the region with overlaid 843 MHz contours from MGPS and 1280 MHz contours from GMRT. Rest of the objects marked are same as in Fig. 7. The central core appears to be a compact H II region, near whose peak lie the mm peak and the high-mass source. The 843 MHz contours show extended emission in addition to the central compact region. The extended emission is only in the southern part and none in the northern part of the compact H II region, indicating the presence of dense molecular cloud in the northern part which is not ionized to the same extent as the southern region. The background 8.0 µm image also shows the diffuse nebular emission in the north. Using the AIPS task JMFIT, the compact cores at both the frequencies were fitted with a Gaussian model to determine the source sizes and the fluxes. The obtained results are given in Table 3. The beam-deconvolved source size for MGPS 843 MHz is found to be much larger in area than that for 1280 MHz. The larger integrated flux density for 843 MHz could be because of this. Fig. 9 shows the maximum resolution image of the region which could be obtained at 1280 MHz (~7 arcsec × 2 arcsec), overlaid on the Herschel 70 µm image. This contours show multiple peaks which is probably indicative of
<table><caption>**Table 3.** Fitting results for compact H II region.</caption><thead><tr><th></th><th>1280 MHz</th><th>843 MHz</th></tr></thead><tbody><tr><td>2D Gaussian fit</td><td></td><td></td></tr><tr><td>size</td><td>60.30 arcsec × 39.71 arcsec</td><td>69.41 arcsec × 76.17 arcsec</td></tr><tr><td>Beam-deconvolved</td><td></td><td></td></tr><tr><td>source size</td><td>23.35 arcsec × 18.31 arcsec</td><td>51.25 arcsec × 54.97 arcsec</td></tr><tr><td>Position Angle (deg)</td><td>160.1</td><td>128.2</td></tr><tr><td>Peak flux</td><td></td><td></td></tr><tr><td>Density (mJy beam<sup>-1</sup>)</td><td>204 ± 9</td><td>258 ± 5</td></tr><tr><td>Integrated flux</td><td></td><td></td></tr><tr><td>density (mJy)</td><td>252 ± 18</td><td>552 ± 16</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** The colour map Herschel 70 µm image with overlaid 1280 MHz contours at 3σ, 4σ, 5σ, 7σ, 11σ, 13σ, 15σ, 17σ, 19σ, 20σ, and 21 σ (where σ ~ 0.4 mJy). The contours are from the maximum resolution (7 arcsec × 2 arcsec) radio continuum image that could be constructed at 1280 MHz. The symbols are same as Fig. 7.</figcaption></figure>

the clumpy nature of the ionized matter in the region. These peaks are mostly coincident with large 70 µm emission, which is to be expected as 70 µm also traces the thermal dust emission.

We tried to estimate the physical parameters of the region, using the lower resolution images, as follows. The Lyman continuum luminosity (in photons s⁻¹) required to generate the observed flux density was determined using the following formula (adapted from Kurtz, Churchwell & Wood 1994, see their equations 1 and 3):

$$S_* \geq \left( \frac{7.59 \times 10^{48}}{a(\nu, T_e)} \right) \left( \frac{S_v}{Jy} \right) \left( \frac{T_e}{K} \right)^{-0.5} \left( \frac{D}{\text{kpc}} \right)^2 \left( \frac{\nu}{\text{GHz}} \right)^{0.1}, \quad (3)$$

where $S_v$ is the integrated flux density in Jy, $D$ is the distance in kiloparsec, $T_e$ is the electron temperature, $a(\nu, T_e)$ is the correction factor, and $\nu$ is the frequency in GHz at which the luminosity is to

be calculated. The dynamical age of the H II region ($t$) can be solved for by using the following equation from Spitzer (1978):

$$R(t) = R_s \left( 1 + \frac{7c_{II}t}{4R_s} \right)^{4/7}, \quad (4)$$

where $R(t)$ is the radius of the H II region at time $t$, $c_{II}$ is the speed of sound in H II region (11 × 10⁵ cm s⁻¹; Stahler & Palla 2005), and $R_s$ is the Strömgren radius. The Strömgren radius ($R_s$, in cm) is given by (Strömgren 1939)

$$R_s = \left( \frac{3S_*}{4\pi n_o^2 \beta_2} \right)^{1/3}, \quad (5)$$

where $n_o$ is the initial ambient density (in cm⁻³), and $\beta_2$ is the total recombination coefficient to the first excited state of hydrogen. For our calculation, we assume a typical value of 10 000 K for $T_e$ (which will imply a value of 0.99 for the correction factor ‘a’; see table 6 of Mezger & Henderson 1967), and the corresponding $\beta_2$
<figure><img src="image_13.png" /><figcaption>Figure 10. The fitted free-free emission model for the H II region. The data points at 843 and 1280 MHz have been marked with solid circles.</figcaption></figure>

of $2.6 \times 10^{-13}$ cm$^3$ s$^{-1}$ (Stahler & Palla 2005). For $n_o$, we use the value of $4.8 \times 10^4$ cm$^{-3}$ from Beltrán et al. (2006). $R(t)$ is taken as the geometric mean of fitted Gaussian source sizes from Table 3. Using these formulae and the 1280 MHz data, we calculated $S_*$ and $t$ as $\sim 10^{47.41}$ photons s$^{-1}$ and $\sim 0.3$ Myr, respectively. If we were to use the 843-MHz flux density data point, then $S_*$ and $t$ come out to be $\sim 10^{47.73}$ photons s$^{-1}$ and $\sim 0.5$ Myr, respectively. Assuming Zero Age Main Sequence (ZAMS), a comparison of log $S_*$ with the tabulated values from Panagia (1973) shows that a spectral type of B0-O9.5 corresponds to this luminosity. Recent calibrations, like Martins, Scherer & Hillier (2005), also suggest a (luminosity class V) spectral type of $\sim$O9.5. So, it seems that the spectral type of the source ionizing the region (assuming a single source) has to be earlier than B0-O9.5 for the ionization in the nebula to be sustained, since there can be absorption of ionizing photons by dust in the region which is often significant (Arthur et al. 2004). The dynamical age is approximately of the order of a few tenths of Myr.

We also try to fit to our data the free-free emission model of Mezger & Henderson (1967), according to which (adapted from Mezger, Schraml & Terzian 1967)

$$S_v/\Omega = 3.07 \times 10^{-2} T_e v^2 (1 - e^{-\tau(v)}) \quad (6)$$

$$\tau(v) = 1.643 \times 10^5 a(v, T_e) T_e^{-1.35} v^{-2.1} n_e^2 l, \quad (7)$$

where, $S_v$ is the integrated flux density (in Jy), $v$ is the frequency (in MHz), $n_e$ is the electron density in cm$^{-3}$, $l$ is the extent of the ionized region in pc, $\tau$ is the optical depth, and $\Omega$ is the solid angle subtended by the source (in steradians). $n_e^2 l$ measures the optical depth in the medium (in cm$^{-6}$ pc), and is called the emission measure. Taking $\Omega$ as $(1.133 \times \theta_{\text{major}} \times \theta_{\text{minor}})$, and the two data points, we fit the above equation using non-linear regression keeping $n_e^2 l$ – the emission measure – as a free parameter. The fit is shown in Fig. 10. The emission measure is thus determined to be $\sim 4.00 \pm 0.09 \times 10^4$ cm$^{-6}$ pc. Since most of the central radio emission is confined within a circle of 120 arcsec diameter, we can assume it to be the extent of the H II region ($\sim 2.1$ pc), and thus the electron density ($n_e$) turns out to be $\sim 138$ cm$^{-3}$. Further, using the formula from Mezger et al. (1967, equation A.5), we calculate the total mass of ionized hydrogen ($M_{H II}$) to be $\sim 16.4$ M$_\odot$. These low values of $n_e$ and $M_{H II}$, as well as the extent, would suggest

that this region might be slightly more evolved than a compact H II region (Kurtz & Franco 2002).

### 4.3 Central region

Fig. 11 shows the central 1 arcmin $\times$ 1 arcmin region of IRAS 16148–5011 in the NIR J and $K_s$ bands, IRAC 3.6, 4.5, and 8.0 µm bands, and the MIPS 24 µm band. The IRAC 3.6 and 8.0 µm bands, besides the continuum emission, also encompass the weak polycyclic aromatic hydrocarbon (PAH) emission feature at 3.3 µm, and strong PAH features at 7.7 and 8.6 µm. The 4.5 µm band does not contain any PAH features, but does contain shocked molecular gas emission from $H_2(v=0-0)$ $S(9, 10, 11)$ and $CO(v=1-0)$. 24-µm emission is mainly the thermal continuum emission from the hot dust (Watson et al. 2008; Churchwell et al. 2009).

The $K_s$ band image has been marked with, among others, the high-mass source detected (green cross, #16 in Table 2) and the MSX peak from Molinari et al. (2008, ‘IRAS 16148–5011 A’) (green box). The marked high-mass source is of $\sim$B2–B3 spectral type (as per tabulated values from Stahler & Palla 2005, see their table 1.1), while the radio flux suggests an ionizing star of type $\ge$ B0-O9.5. In addition, Molinari et al. (2008), using SED modelling (assuming an embedded ZAMS source) and the MSX flux values, have estimated the spectral type of ‘IRAS 16148–5011 A’ source as O8. Therefore, it appears that there could be further embedded high-mass source(s) in this region, in addition to those seen in NIR and MIR, for the values to be consistent. It is possible that the cores seen by radio contours in Fig. 9 could be hosting such high-mass source(s) and contributing to the radio luminosity.

The other sources associated with this central part of the nebula have been marked on the 3.6 µm image (#13 and #28 from Table 2, Section 3.2). The source #13 appears to be a highly embedded intermediate-mass young source ($\sim 6.17$ M$_\odot$, $\sim$a few $10^4$ yr; see Section 3.2), even visible at 24 µm. Source #28 appears to be a much older and evolved source ($\sim 8.48$ M$_\odot$, $\sim$a few $10^6$ yr), possibly shrouded in the surrounding nebulosity (interstellar extinction of $A_V \sim 18.38$ mag), leading to a lack of detection in the J band (though it is the second brightest source after the 9.6 M$_\odot$ – #16 from Table 2 – source in the central region in $K_s$ band). The radio continuum emission contribution from this source will be much lower, hardly making a difference even when taken in combination with the emission from the high-mass source, and thus will not affect the inferences regarding radio emission above.

For H II regions, the PAH emissions serve as useful diagnostic tools as they are characteristic of photodissociation regions (PDRs). As a high-mass star ionizes its natal medium, the UV radiation destroys the PAH in its surroundings. However, at the boundary of the H II region produced, the UV intensity falls off, and a PDR is formed where the PAH’s are merely highly excited, leading to strong emission (Povich et al. 2007). The PDR is supposed to be the transition region between the ionized and neutral matter. This usually results in a ‘bubble’ morphology, where H II regions surrounded by ring-like PAH emission are seen (Churchwell et al. 2006). A similar morphology can be seen in the 3.6 and the 8.0 µm image in Fig. 11. A ‘semiring’ in the western half of the image is seen (marked by cyan arrows on the 8.0-µm image). We note that this feature is also faintly seen in the 4.5 µm image, though this could just be continuum emission. That this ring-like feature does not appear symmetric (i.e. there is no clear eastern ‘semi-ring’, though faint indications are seen) could possibly be due to eastern denser and/or non-homogeneous molecular cloud, or projection
<figure><img src="image_14.png" /><figcaption>**Figure 11.** 1 arcmin × 1 arcmin central region of IRAS 16148-5011. The images (left-to-right, top-to-bottom) are J, Ks, 3.6, 4.5, 8.0, and 24 µm. North is up and east is to the left. The symbols on the Ks band image are same as Fig. 7 (green cross is #16 from Table 2). The numbered sources on the 3.6 µm image are from Table 2. The cyan arrows on the 8.0 µm image mark the 'ring-like' morphology. The white patch on the 24 µm image is the saturated region.</figcaption></figure>

effects. Volk et al. (1991), on the basis of low-resolution spectra, had classified IRAS 16148-5011 as a PAH source. Bulk of the radio emission (see Fig. 8), as well as the 24 µm emission from hot dust, is confined within this ring-like feature as expected (Watson et al. 2008). Some of the emission from this ring-like structure could also be due to swept-up material by the expanding ionization front of the H II region.

## 4.4 Herschel results

Thermal emission from cold dust lies in the FIR wavelength range, and thus its analysis can be used to obtain the physical parameters like dust temperature and column density of a region (Battersby et al. 2011; Launhardt et al. 2013). This was carried out by the SED modelling of the thermal dust emission, whose Rayleigh-Jeans regime is covered by the Herschel FIR bands (160–500 µm), using the following order of steps.

First of all, the surface brightness unit for all the images was converted to Jy pixel⁻¹. Since the PACS image is already in Jy pixel⁻¹, this step was only required for the SPIRE images (whose units are in MJy sr⁻¹), and was carried out using the pixel scales for the respective SPIRE bands. Next, the 160–350 µm images were convolved to the resolution of the 500 µm image (~36 arcsec, lowest among all images) using the convolution kernels of Aniano et al. (2011), and regridded to a pixel scale of 14 arcsec (same as 500 µm image). Using these final reworked images with same resolution and pixel scale, a background flux level, $I_{bg}$, was determined from a 'smooth' (i.e. no abrupt clumpy regions) and relatively 'dark' patch of the sky. The distribution of individual pixel values in the dark patch of the sky, for each of the bands, was fitted with a Gaussian iteratively,

rejecting the pixel values outside ± 2σ in each iteration, till the fit converged. The same patch of the sky was used for each band. The background flux level, $I_{bg}$, was thus determined as 0.29, 2.67, 1.27, and 0.45 Jy pixel⁻¹ for the 160, 250, 350, and 500 µm images, respectively. We note that the 70 µm image, though available, was not used here as the optically thin assumption might not hold true at this wavelength (Lombardi et al. 2014).

Modified blackbody fitting was subsequently carried out on a pixel-by-pixel basis using the following formulation (Battersby et al. 2011; Sadavoy et al. 2012; Nielbock et al. 2012; Launhardt et al. 2013):

$$S_v(v) - I_{bg}(v) = B_v(v, T_d)\Omega(1 - e^{-\tau(v)}), \quad (8)$$

with

$$\tau(v) = \mu_{H_2} m_{H} \kappa_{v} N(H_2), \quad (9)$$

where, v is the frequency, $S_v(v)$ is the observed flux density, $I_{bg}(v)$ is the background flux in that particular band (estimated above), $B_v(v, T_d)$ is the Planck function, $T_d$ is the dust temperature, Ω is the solid angle (in steradians) from where the flux is obtained (just the solid angle subtended by a 14 arcsec × 14 arcsec pixel here), τ(v) is the optical depth, $\mu_{H_2}$ is the mean molecular weight (adopted as 2.8 here), $m_H$ is the mass of hydrogen, $\kappa_v$ is the dust opacity, and $N(H_2)$ is the column density. For opacity, we adopt a functional form of $\kappa_v = 0.1$ (v/1000 GHz)β cm² g⁻¹, with β = 2 (see Hildebrand 1983, Beckwith et al. 1990, André et al. 2010). For each pixel, equation (8) was fitted using the four data points, keeping $T_d$ and $N(H_2)$ as free parameters. Pixels for which the fit did not converge, or the error was larger than 10 per cent, had their values taken as the median of eight immediate-neighbour pixels. The final obtained temperature
<figure><img src="image_15.png" /><figcaption>**Figure 12.** (a) Dust temperature map, and (b) column density map of the region around IRAS 16148-5011, derived using the SED fitting to the thermal dust emission.</figcaption></figure>

and column density maps of the wider region (to clearly discern the morphological features) surrounding IRAS 16148-5011 are shown in Fig. 12.

From the temperature (Fig. 12a) and column density (Fig. 12b) maps, we obtain the peak values as ~30 K and $3.3 \times 10^{22} \text{ cm}^{-2}$,

respectively, for IRAS 16148-5011. Fontani et al. (2005) provide a temperature of 38 K (using greybody fit with 60 µm, 100 µm, and 1.2 mm data), and a (beam-averaged) column density of $3 \times 10^{23} \text{ cm}^{-2}$ (using C$^{17}$O molecular line observations) for IRAS 16148-5011. Here our peak temperature is about 20 per cent
lower, and the column density estimate an order of magnitude smaller. The difference in temperature seems to be due to their fitting limitations owing to sparse data points, as well as the fact that the 60 µm emission might not be optically thin. It should be noted, however, that the dust temperature distribution from Fontani et al. (2005, for their analysed catalogue of IRAS sources) peaks at ~30 K. Column density disparity can probably be explained by the dependence of Fontani et al. (2005) on the molecular abundance of the rare C¹⁷O, whose value can have wide variations (Redman et al. 2002; Walsh et al. 2010).

The column density map (Fig. 12b) displays three peaks – the central IRAS 16148–5011 object, a peak to its north-east, and a peak to its south. These north-east and south peaks appear to be associated with cold clumps, as is evident from the temperature map which shows T<sub>d</sub> < 20 K at their positions. The column density map shows that all the peaks appear to be connected by weak filamentary features, similar to previous results for myriad regions (see André 2013, and references therein). For the central IRAS 16148–5011 object, we estimated the associated clump dimensions, using the ‘CLUMPFIND’ software (Williams, de Geus & Blitz 1994), to be ~62 arcsec × 46 arcsec (i.e. 1.1 pc × 0.8 pc at 3.6 kpc). Due to the low resolution of the column density image generated here, it does not appear possible to resolve further sub-clumps. The mass of the clump can be estimated by

$$M_{\text{clump}} = \mu_{\text{H}_2} m_{\text{H}} \Sigma N(\text{H}_2) \text{Area}_{\text{pixel}} \quad (10)$$

$$= \mu_{\text{H}_2} m_{\text{H}} \text{Area}_{\text{pixel}} \Sigma N(\text{H}_2), \quad (11)$$

i.e. calculating the mass in each pixel and then summing over all the pixels which constitute the clump. Using the Σ N(H₂) returned by the ‘CLUMPFIND’ software, the total clump mass was calculated to be ~1023 M<sub>☉</sub>.

Another noticeable feature in this map is that the immediate northern vicinity of the central region exhibits higher column density than the immediate southern portion. This suggests dense nebula in the northern part, affirming the inference also drawn from extinction and cluster analysis in Section 4.1 (also see Fig. 7). This is also consonant with the fact that the dense nebular emission is in the northern part, with no extended radio emission seen there (unlike in the southern part). To get an estimate of the visual extinction in this northern part, we use the relation (N(H₂)/A<sub>V</sub>) = 0.94 × 10<sup>21</sup> molecules cm<sup>-1</sup> mag<sup>-1</sup> (adapted from Bohlin, Savage & Drake 1978, assuming a total-to-selective extinction ratio R<sub>V</sub> = 3.1, and that the gas is in molecular form). Now, typical N(H₂) seen here is ≥2 × 10<sup>22</sup> cm<sup>-3</sup>, which implies an A<sub>V</sub> > 20 mag. This high value of A<sub>V</sub> is probably the reason why very few YSOs are associated with this northern nebular emission (see Fig. 7). If the value of R<sub>V</sub> were to be larger, as has been conjectured for dense environments, then it will lead to an increase in A<sub>V</sub>.

## 5 DISCUSSION

The IRAS 16148–5011 region appears to be hosting an infrared cluster containing high-mass star(s) embedded in the nebula, and could serve as a future template to study the high-mass star formation process. Age estimates for any cluster are usually plagued with uncertainties (see Lada & Lada 2003, for a discussion), and various proxies are often used to ascertain the evolutionary stage of a cluster. One such is the ratio of Class II to Class I sources (as Class II sources are older; Schmeja, Klessen & Froebrich 2005; Gutermuth et al. 2009; Beerer et al. 2010). We have detected 133 Class II and 7 Class I sources in this region. However, here it seems

that the large ratio (19) could be because Class I sources are deeply embedded and their detection is affected by large nebular extinction, though varying ratio values have been estimated for other clusters (Gutermuth et al. 2009).

The two brightest (NIR) sources in the central region appear to be high-mass (see Section 4.3) with indications of further embedded high-mass sources. Also in literature, Beltrán et al. (2006) report two clumps within 90 arcsec of the IRAS 16148–5011 IRAS catalogue position – ‘16148-5011 Clump 2’ (206 M<sub>☉</sub>) and ‘16148-5011 Clump 3’ (42 M<sub>☉</sub>) – detected using 1.2 mm emission and a dust temperature of 30 K. While ‘Clump 2’, is almost coincident with the ‘16148-5011 MM 1’ peak from Molinari et al. (2008, marked with a diamond symbol on Fig. 11), ‘Clump 3’ coordinates are same as the IRAS 16148–5011’s IRAS catalogue coordinates (green plus in Fig. 11). The presence of these massive clumps could also point towards ongoing high-mass star formation in the region.

Though this part of the sky also harbours other IRAS sources (Karnik et al. 2001), IRAS 16148–5011 does not appear to be a part of any larger star-forming complex as such, but it seems to be connected to other clumps via filaments (similar to the hub-filament morphology of Myers 2009). The high-mass stellar cluster formation in the associated molecular cloud appears to be spontaneous. However, along the rather sharp southern boundary of the molecular cloud, the results from cluster analysis show a stellar sub-cluster forming (Fig. 7), which could partly be due to the triggering by the expanding ionization front of the H II region.

## 6 CONCLUSIONS

The main conclusions of this paper, resulting from a multiwavelength study involving NIR, MIR, FIR, and radio continuum data, are as follows.

(i) 7 Class I and 133 Class II YSOs are identified using a combination of NIR and MIR data. A 9.6-M<sub>☉</sub> high-mass source is found to be associated with the central nebula.

(ii) Low-frequency radio emission reveals a compact H II region, with extended emission in the southern part. Lyman continuum photon luminosity calculation gives B0-O9.5 as the lower limit for the spectral type of the ionizing source (assuming single source). The dynamical age of the H II region is in 0.3–0.5 Myr range. SED modelling of the free-free emission yields an electron density of 138 cm<sup>-3</sup>. The mass of the ionized hydrogen is calculated to be ~16.4 M<sub>☉</sub>. The high-resolution 7 arcsec × 2 arcsec contour map shows clumpy ionized gas in the region.

(iii) The central nebular region shows a ring-like PAH emission feature near the borders of the compact H II region, tracing the PDR. There appear to be three NIR- and MIR-visible central sources, of masses ~6.2, 9.60, and 8.5 M<sub>☉</sub>. Based on the incongruence between the total radio flux from these sources and the flux obtained from radio observations, and literature estimates of an early-type star using MIR SED fitting, it is possible that there could be high-mass embedded source(s) present in this region.

(iv) Dust temperature and column density maps are obtained using SED modelling of the thermal dust emission. The peak temperature and column density values are 30 K and 3.3 × 10<sup>22</sup> cm<sup>-2</sup>, respectively, for IRAS 16148–5011. The column density map reveals that the immediate northern vicinity of IRAS 16148–5011, which contains the nebular emission seen prominently at MIR wavelengths, has a large extinction of A<sub>V</sub> > 20 mag. This map also shows that weak filamentary structures join IRAS 16148–5011 to nearby cold clumps. The size and mass of the clump associated with
IRAS 16148–5011 is estimated to be ~1.1 pc × 0.8 pc (at 3.6 kpc) and ~1023 M⊙, respectively.

Future observations of individual objects in spectral lines, deeper IR data to get a full stellar census up to below brown-dwarf limit, further molecular line observations which probe high column densities will help put the star formation scenario on a firm footing and help study high-mass star formation.

## ACKNOWLEDGEMENTS

We thank the anonymous referee for a thorough reading of the manuscript, and for the useful comments and suggestions which helped improve its scientific content. The authors thank the staff of IRSF in South Africa, a joint partnership between S.A.A.O and Nagoya University of Japan; and GMRT managed by National Center for Radio Astrophysics of the Tata Institute of Fundamental Research (TIFR) for their assistance and support during observations. DKO was supported by the National Astronomical Observatory of Japan (NAOJ), Mitaka, through a fellowship, during which a part of this work was done.

## REFERENCES

André P., 2013, Invited Discourse at the IAU XXVIII General Assembly (Beijing, Aug. 2012), preprint (arXiv:1309.7762)

André P. et al., 2010, A&A, 518, L102

Aniano G., Draine B. T., Gordon K. D., Sandstrom K., 2011, PASP, 123, 1218

Arthur S. J., Kurtz S. E., Franco J., Albarrán M. Y., 2004, ApJ, 608, 282

Balog Z., Kenyon S. J., Lada E. A., Barsony M., Vinkó J., Gáspaf A., 2004, AJ, 128, 2942

Battersby C. et al., 2011, A&A, 535, A128

Beckwith S. V. W., Sargent A. I., Chini R. S., Guesten R., 1990, AJ, 99, 924

Beerer I. M. et al., 2010, ApJ, 720, 679

Beltrán M. T., Brand J., Cesaroni R., Fontani F., Pezzuto S., Testi L., Molinari S., 2006, A&A, 447, 221

Benjamin R. A. et al., 2003, PASP, 115, 953

Bessell M. S., Brett J. M., 1988, PASP, 100, 1134

Bohlin R. C., Savage B. D., Drake J. F., 1978, ApJ, 224, 132

Cantalupo C. M., Borrill J. D., Jaffe A. H., Kisner T. S., Stompor R., 2010, ApJS, 187, 212

Carpenter J. M., 2001, AJ, 121, 2851

Casertano S., Hut P., ApJ, 298, 80

Chan S. J., Henning T., Schreyer K., 1996, A&AS, 115, 285

Chavarría L. A., Allen L. E., Hora J. L., Brunt C. M., Fazio G. G., 2008, ApJ, 682, 445

Churchwell E. et al., 2006, ApJ, 649, 759

Churchwell E. et al., 2009, PASP, 121, 213

Cohen J. G., Persson S. E., Elias J. H., Frogel J. A., 1981, ApJ, 249, 481

Dutra C. M., Bica E., Soares J., Barbuy B., 2003, A&A, 400, 533

Flaherty K. M., Pipher J. L., Megeath S. T., Winston E. M., Gutermuth R. A., Muzerolle J., Allen L. E., Fazio G. G., 2007, ApJ, 663, 1069

Fontani F., Beltrán M. T., Brand J., Cesaroni R., Testi L., Molinari S., Walmsley C. M., 2005, A&A, 432, 921

Grave J. M. C., Kumar M. S. N., 2009, A&A, 498, 147

Green A. J., Cram L. E., Large M. I., Ye T., 1999, ApJS, 122, 207

Griffin M. J. et al., 2010, A&A, 518, L3

Gutermuth R. A., Megeath S. T., Myers P. C., Allen L. E., Pipher J. L., Fazio G. G., 2009, ApJS, 184, 18; erratum, 2010, ApJS, 189, 352

Haslam C. G. T., Salter C. J., Stoffel H., Wilson W. E., 1982, A&AS, 47, 1

Haynes R. F., Caswell J. L., Simons L. W. J., 1979, Aust. J. Phys., 48, 1

Hildebrand R. H., 1983, QJRAS, 24, 267

Hillenbrand L. A., Strom S. E., Vrba F. J., Keene J., 1992, ApJ, 397, 613

Karnik A. D., Ghosh S. K., Rengarajan T. N., Verma R. P., 2001, MNRAS, 326, 293

Kurtz S., Franco J., 2002, Rev. Mex. Astron. Astrofis., 12, 16

Kurtz S., Churchwell E., Wood D. O. S., 1994, ApJS, 91, 659

Lada C. J., 1987, in Peimbert M., Jugaku J., eds, Proc. IAU Symp. 115, Star Forming Regions. Reidel, Dordrecht, p. 1

Lada C. J., Adams F. C., 1992, ApJ, 393, 278

Lada E. A., Lada C. J., 1995, AJ, 109, 1682

Lada C. J., Lada E. A., 2003, ARA&A, 41, 57

Lada E. A., Depoy D. L., Evans N. J., II, Gatley I., 1991, ApJ, 371, 171

Lada C. J., Young E. T., Greene T. P., 1993, ApJ, 408, 471

Lada C. J., Lada E. A., Clemens D. P., Bally J., 1994, ApJ, 429, 694

Launhardt R. et al., 2013, A&A, 551, A98

Lombardi M., Bouy H., Alves J., Lada C. J., 2014, A&A, 566, A45; erratum, 2014, A&A, 568, 1

Lumsden S. L., Hoare M. G., Urquhart J. S., Oudmaijer R. D., Davies B., Mottram J. C., Cooper H. D. B., Moore T. J. T., 2013, ApJS, 208, 11

MacLeod G. C., van der Walt D. J., North A., Gaylard M. J., Galt J. A., Moriarty-Schieven G. H., 1998, AJ, 116, 2936

Mallick K. K., Kumar M. S. N., Ojha D. K., Bachiller R., Samal M. R., Pirogov L., 2013, ApJ, 779, 113

Mallick K. K. et al., 2014, MNRAS, 443, 3218

Martins F., Schaeerer D., Hillier D. J., 2005, A&A, 436, 1049

Massey P., 1998, in Gilmore G., Howell D., eds, ASP Conf. Ser. Vol. 142, The Stellar Initial Mass Function. Astron. Soc. Pac., San Francisco, p. 17

Megeath S. T., 1996, A&A, 311, 135

Meyer M. R., Calvet N., Hillenbrand L. A., 1997, AJ, 114, 288

Mezger P. G., Henderson A. P., 1967, ApJ, 147, 471

Mezger P. G., Schraml J., Terzian Y., 1967, ApJ, 150, 807

Molinari S., Pezzuto S., Cesaroni R., Brand J., Faustini F., Testi L., 2008, A&A, 481, 345

Molinari S. et al., 2010, PASP, 122, 314

Mookerjea B., Kramer C., Nielbock M., Nyman L.-Å., 2004, A&A, 426, 119

Myers P. C., 2009, ApJ, 700, 1609

Nagashima C. et al., 1999, in Nakamoto T., ed., Proc. Star Formation 1999, Star Formation 1999. Nobeyama Radio Observatory, Japan, p. 397

Nagayama T. et al., 2003, in Masanori I., Moorwood A. F. M., eds, Proc. SPIE Conf. Ser. Vol. 4841, Instrument Design and Performance for Optical/Infrared Ground-based Telescopes. SPIE, Bellingham, p. 459

Nielbock M. et al., 2012, A&A, 547, A11

Ojha D. K. et al., 2004a, ApJ, 608, 797

Ojha D. K. et al., 2004b, ApJ, 616, 1042

Panagia N., 1973, AJ, 78, 929

Persson S. E., Murphy D. C., Krzeminski W., Roth M., Rieke M. J., 1998, AJ, 116, 2475

Pilbratt G. L. et al., 2010, A&A, 518, L1

Poglitsch A. et al., 2010, A&A, 518, L2

Povich M. S. et al., 2007, ApJ, 660, 346

Redman M. P., Rawlings J. M. C., Nutter D. J., Ward-Thompson D., Williams D. A., 2002, MNRAS, 337, L17

Robitaille T. P., 2008, in Beuther H., Linz H., Henning T., eds, ASP Conf. Ser. Vol. 387, Massive Star Formation: Observations Confront Theory. Astron. Soc. Pac., San Francisco, p. 290

Robitaille T. P., Whitney B. A., Indebetouw R., Wood K., Denzmore P., 2006, ApJS, 167, 256

Robitaille T. P., Whitney B. A., Indebetouw R., Wood K., 2007, ApJS, 169, 328

Sadavoy S. I. et al., 2012, A&A, 540, A10

Sanchawala K. et al., 2007, ApJ, 667, 963

Schmeja S., 2011, Astron. Nachr., 332, 172

Schmeja S., Klessen R. S., Froebrich D., 2005, A&A, 437, 911

Schmeja S., Kumar M. S. N., Ferreira B., 2008, MNRAS, 389, 1209

Sharma S., Pandey A. K., Ojha D. K., Chen W. P., Ghosh S. K., Bhatt B. C., Maheswar G., Sagar R., 2007, MNRAS, 380, 1141

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

Spitzer L., 1978, Physical Processes in the Interstellar Medium. Wiley, New York
Stahler S. W., Palla F., 2005, The Formation of Stars. Wiley-VCH, New York

Strömgren B., 1939, ApJ, 89, 526

Swarup G., Ananthakrishnan S., Kapahi V. K., Rao A. P., Subrahmanya C. R., Kulkarni V. K., 1991, Curr. Sci, 60, 95

Tej A., Ojha D. K., Ghosh S. K., Kulkarni V. K., Verma R. P., Vig S., Prabhu T P., 2006, A&A, 452, 203

Vig S., Ghosh S. K., Ojha D. K., Verma R. P., Tamura M., 2014, MNRAS, 440, 3078

Volk K., Kwok S., Stencel R. E., Brugel E., 1991, ApJS, 77, 607

Walsh A. J., Thorwirth S., Beuther H., Burton M. G., 2010, MNRAS, 404, 1396

Watson C. et al., 2008, ApJ, 681, 1341

Whitney B. A., Wood K., Bjorkman J. E., Wolff M. J., 2003a, ApJ, 591, 1049

Whitney B. A., Wood K., Bjorkman J. E., Cohen M., 2003b, ApJ, 598, 1079

Williams J. P., de Geus E. J., Blitz L., 1994, ApJ, 428, 693

Zinnecker H., Yorke H. W., 2007, ARA&A, 45, 481

Zinnecker H., McCaughrean M. J., Wilking B. A., 1993, Protostars and Planets III. Univ. Arizona Press, Tucson, AZ, p. 429

## SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this article:

<caption>**Table 1.** YSOs identified using NIR and MIR data (http://mnras.oxfordjournals.org/lookup/suppl/doi:10.1093/mnras/stu2584/-/DC1).</caption>

Please note: Oxford University Press are not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

This paper has been typeset from a TeX/LATEX file prepared by the author.