Sustaining Star Formation in the Galactic Star Cluster M 36?

The Gaia Data Release 2 (Gaia Collaboration et al., 2018, Gaia DR2) furnishes precise astrometry on five parameters, i.e., celestial coordinates, trigonometric parallaxes, and proper motions for more than 1.3 billion objects from the observations and analysis of the European Space Agency Gaia satellite (Gaia Collaboration et al., 2016). In addition to the homogeneous astrometry over the whole sky, Gaia DR2 is further enriched with photometry of three broad-band magnitudes in $G$ (330–1050 nm), $G_{BP}$ (330–680 nm), and $G_{RP}$ (630–1050 nm) with unprecedented accuracy. The $G$ band photometry ranges between 3 mag and 21 mag, although stars with $G\lesssim 6$ mag generally have inferior astrometry due to calibration issues. The median uncertainty in parallax is about 0.04 mas for bright ($G\lesssim 14$ mag) sources, 0.1 mas at $G=17$ mag, and 0.7 mas at $G=20$ mag (Luri et al., 2018). The corresponding uncertainties in proper motion components are 0.05, 0.2, and 1.2 mas yr^-1, respectively (Lindegren et al., 2018). The Gaia DR2 also provides radial velocity measurements for sources with $G\lesssim 13$ mag, and a photometric catalogue of about 0.5 million variable stars (Gaia Collaboration et al., 2018).

The near-infrared $J$-, $H$-, and $K$-band photometric data were obtained from the UKIDSS DR10PLUS Galactic Plane Survey (GPS; Lawrence et al. 2007) and the 2MASS Point Source Catalog (PSC; Skrutskie et al. 2006). UKIDSS has a finer angular resolution ($0\farcs 8$ FWHM, $0\farcs 2$ pixels at $K$) compared to 2MASS (1$\arcsec$ FWHM, 2$\arcsec$ pixels), as well as fainter limits. The UKIDSS GPS has median 5$\sigma$ detection limits at $J=19.77$ mag, $H=19.00$ mag, and $K=18.05$ mag (Lucas et al., 2008), whereas 2MASS is limited to $J=15.8$, $H=15.1$, and $K_{S}=14.3$ mag (S/N = 10; Skrutskie et al. 2006). Selection of reliable sources from the UKIDSS catalog was carried out by utilizing the Structured Query Language (SQL¹¹1http://wsa.roe.ac.uk/sqlcookbook.html) query and adopted the criteria published by Lucas et al. (2008), which eliminates saturated, non-stellar, and unreliable sources near the sensitivity limits. The UKIDSS GPS shows saturation limits at $J=13.25$ mag, $H=12.75$ mag, and $K=12.00$ mag (Lucas et al., 2008). To supplement the UKIDSS saturated sources with 2MASS photometry, we set 0.5 mag fainter than quoted values (Alexander et al., 2013; Dutta et al., 2015). For the near-infrared bands, photometric uncertainty $\lesssim 0.1$ mag was considered as quality criteria, which give S/N $\gtrsim 10$.

The mid-infrared photometric data for point sources toward the M 36 cluster was obtained from the Spitzer Warm Mission (Hora et al., 2012) survey. Magnitudes for Infrared Array Camera (IRAC; Fazio et al. 2004) [3.6] µm and [4.5] µm bands were downloaded from the Glimpse360²²2http://www.astro.wisc.edu/sirtf/glimpse360/ catalog (Program Name/Id: WHITNEY_GLIMPSE360_2/61070) with a pixel scale of 1$\farcs$2 pixel^-1. We restricted the sources with uncertainty $<$ 0.2 mag for all the IRAC bands to achieve good quality photometric catalog.

The Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST, also called the Guoshoujing Telescope) is a reflecting Schmidt telescope (Wang et al., 1996) with an effective aperture of 3.6–4.9 m, a focal length of 20 m, and a field of view of $5\arcdeg$ (Cui et al., 2012; Zhao et al., 2012). After a five-year regular survey³³3http://www.lamost.org/public/node/311?locale=en conducted between 2012 to 2017, more than eight million spectra were obtained, with a spectral resolution of R $\sim 1800$ over the wavelength range of 370–900 nm. The raw 2D spectra are processed uniformly with the LAMOST 2D reduction pipeline (Luo et al., 2015) to generate the 1D spectra, from which stellar parameters, including radial velocity (RV), effective temperature (T_eff), surface gravity ($\log g$), and metallicity ([Fe/H]), are then derived. Both the 1D spectra and the stellar parameters are publicly available via the LAMOST data releases⁴⁴4http://www.lamost.org/ (Luo et al., 2012, 2015). The stellar parameters are collected from LAMOST DR5⁵⁵5http://dr5.lamost.org/ for this work.

The CO observations toward M 36 were conducted on 2 March 2014 using the 13.7 m millimetre-wavelength telescope of the Purple Mountain Observatory (PMO) in Delingha, China, as a part of the Milky Way Imaging Scroll Painting (MWISP) project dedicated to map the molecular gas along the northern Galactic plane (Su et al., 2019). The nine-beam Superconducting Spectroscopic Array Receiver system was used as the front end, and each Fast Fourier transform spectrometer with a bandwidth of 1 GHz provided 16,384 channels and a spectral resolution of 61 kHz (see the details in Shan et al., 2012). The molecular lines of ¹²CO ($J$=1–0) in the upper sideband, and ¹³CO ($J$=1–0) together with C¹⁸O ($J$=1–0) in the lower sideband were observed simultaneously. Typical system temperature is 150–200 K in the lower sideband and 250–300 K in the upper sideband. The on-the-fly mode was applied with typical sample steps $10\arcsec$–$15\arcsec$ and scanned along both the Galactic longitude and latitude, in order to reduce the fluctuation of noise. The on-the-fly raw data were then resampled into $30\arcsec\times 30\arcsec$ grids and mosaicked to a FITS cube using the GILDAS software package⁶⁶6http://www.iram.fr/IRAMFR/GILDAS.

For the data reported here, the typical sensitivity was about 0.3 K (T_mb) in ¹³CO and C¹⁸O at the resolution of 0.17 km s^-1, and 0.5 K (T_mb) in ¹²CO at the resolution of 0.16 km s^-1. The beam widths were about $55\arcsec$ and $52\arcsec$ at 110 GHz and 115 GHz, respectively. The pointing of the telescope has an rms accuracy of about $5\arcsec$. It should be noted that any results presented in the figures and tables are on the brightness temperature scale (T${}_{\rm R}^{*}$), corrected with the beam efficiencies (from the status report of the telescope⁷⁷7http://www.radioast.nsdc.cn/mwisp.php) using $T_{\rm mb}$ = $T_{\rm A}^{\ast}/\eta_{\rm mb}$, with a calibration accuracy within 10%.

M 36 has been known to have an elongated shape, with the aspect ratio of 0.2–0.3, tilted some $20\arcdeg$ away from the plane (Chen et al., 2004; Kharchenko et al., 2009). Using King profiles, Piskunov et al. (2007) estimated an empirical angular radius $16\farcm 2$, and a total mass $\sim 200$ M_☉ within the tidal radius of 9 pc.

In order to determine the cluster extension, we utilized the UKIDSS $K$-band data. The central coordinates ($\alpha_{2000}$ = $05^{\rm{h}}36^{\rm{m}}18^{\rm{s}}$, $\delta_{2000}$ = $+34^{\rm{d}}08^{\rm{m}}24^{\rm{s}}$) as informed in the SIMBAD database⁸⁸8http://simbad.u-strasbg.fr/simbad/sim-fid are adopted to construct the radial density profile with respect to the background. The radial profile is generated by counting the number of stars inside, and divided by the area of, each concentric annulus of 1.5 arcmin width, up to a radius of 35 arcmin. The projected radial stellar density distribution toward M 36 is shown in Figure 1. The peak radial density is estimated to be $\sim 26.4\pm 2.5$ stars arcmin^-2, in contrast to the mean background stellar density of $\sim 14.8\pm 0.3$ stars arcmin^-2. The density profile starts to blend with that of the field population at radius $\sim 13\farcm 5\pm 0\farcm 4$, which we adopt as the cluster boundary. Our estimation of the cluster radius is consistent with the earlier work of Sharma et al. (2006), determined as $14\arcmin$ using optical photometry. The fluctuation around radius $6\arcmin$ arises from the possible boundary of the associated molecular cloud (Section 5.1 and 5.2).

The velocity dispersion of the cluster members is typically much lower than that of the field stars, so the members can be generally distinguished by their uniform kinematics (Kraus & Hillenbrand, 2007). The kinematic membership is tested on the basis of the Gaia DR2 proper motion and parallax data. The proper motion distribution of all the sources within a radius of $30\arcmin$ from the center ($\alpha_{2000}$ = $05^{\rm{h}}36^{\rm{m}}18^{\rm{s}}$, $\delta_{2000}$ = $+34^{\rm{d}}08^{\rm{m}}24^{\rm{s}}$) of M 36 is shown in Figure 2(a). There are two distinct enhancements, one for the cluster members ($\mu_{\alpha}\cos\delta\approx-0.2$ mas yr^-1, $\mu_{\delta}\approx-3.4$ mas yr^-1), and the other for field stars ($\mu_{\alpha}\cos\delta\approx 0.48$ mas yr^-1, $\mu_{\delta}\approx-1.76$ mas yr^-1). Conceivably, around the peak for the cluster, the distribution is dominated by members, whereas away from peak the contamination by field stars becomes prominent. A clear separation between the cluster and field motion enables us to select member stars of M 36 relatively reliably.

The proper motion distribution around the cluster peak is analysed for a square region of width 1.5 mas yr^-1. The zoomed-in distributions for all the sources within this box region for a radius of $30\arcmin$ and of $13\farcm 5$ from the cluster center are shown in Figure 2(b). We estimated the relative contribution of members versus field stars by projecting the distribution along $\mu_{\alpha}~{}cos\delta$ and $\mu_{\delta}$. Each background subtracted rescaled histogram was then fitted with a Gaussian function to quantify the proper motion centroid. The expectation values ($\mu$) of the Gaussian fits, assuming a radius of $30\arcmin$, give $\mu_{\alpha}\cos\delta=-0.17\pm 0.01$ mas yr^-1 and $\mu_{\delta}=-3.34\pm 0.02$ mas yr^-1 with standard deviations ($\sigma$) of $0.12\pm 0.01$ and $0.14\pm 0.02$ mas yr^-1, respectively. Likewise, should a radius of $13.5\arcmin$ be adopted, we derived $\mu_{\alpha}~{}\cos\delta=-0.13\pm 0.01$ mas yr^-1 and $\mu_{\delta}=-3.36\pm 0.02$ mas yr^-1 with standard deviations of $0.16\pm 0.02$ and $0.16\pm 0.02$ mas yr^-1, respectively. We adopted the proper motion center of the cluster as the mean of the above two values, i.e., $\mu_{\alpha}\cos\delta=-0.15\pm 0.01$ mas yr^-1 and $\mu_{\delta}=-3.35\pm 0.02$ mas yr^-1. The optimal range of proper motions is selected within a radius of 0.5 mas yr^-1 ($\simeq 3\sigma$) from this center. In Figure 2(e) and 2(f), the normalization of the histograms permits us to estimate the integrated number of members along $\mu_{\alpha}~{}cos\delta$ as 208, and along $\mu_{\delta}$ as 257, enabling an estimate of the number of missing members (incompleteness) after the parallax constraint is introduced.

We further constrained the membership by including Gaia DR2 parallax measurements. Figure 3 compares the parallaxes for stars in the cluster region (radius $\simeq 13\farcm 5$), with those in a nearby control field, centered around R.A. $=05^{\rm h}36^{\rm m}14^{\rm s}$, Decl. $=+33\arcdeg 21\arcmin 22\arcsec$ (J2000) ($\ell=175\fdg 1877$; $b=0\fdg 6383$), roughly $47\arcmin$ south, with the same sky area as the cluster region. Overplotted are the stars that are located within the cluster region, and satisfying the proper motion criteria (average proper motion $<0.5$ mas yr^-1 from the proper motion center). It is assuring to see a sharp rise in the parallax range between $0.70\pm 0.11$ mas and $0.90\pm 0.08$ mas. The parallax range 0.7–0.9 mas is thus adopted to further refine the membership selection. The number of stars satisfying the positional (inside the cluster region), kinematic (inside the proper motion range), and parallax criteria are 207, which follows to a variation in the incompleteness factor of $\sim$1–19%.

We have identified a total of 200 member candidates by imposing the astrometric and kinematic criteria within the cluster region, that is, within a radius of $\simeq 13\farcm 5$, with proper motion less than 0.5 mas yr^-1 from the cluster proper motion center ($\mu_{\alpha}$ $cos\delta=-0.15$ mas yr^-1, $\mu_{\delta}=-3.35$ mas yr^-1), and parallax between 0.7 mas to 0.9 mas. The median of the parallax for the members is obtained as $0.82\pm 0.07$ mas, corresponding to a distance of $1.20\pm 0.13$ kpc. The photometric and astrometric parameters of all the member candidates are listed in Table 1. The spatial distributions of the members are shown in Figure 4. Among a total of 7948 sources within the cluster region, 757 are found to satisfy the parallax criteria alone and 402 sources satisfy the proper motion criteria alone, with 200 sources satisfying all the above criteria. It is to be noted that, among the members list, 94% sources have parallax errors less than 20%. Moreover, for those sources that satisfy the membership criteria, but have higher parallax errors (20% or above), we compared their distribution with the theoretical isochrones in both color-magnitude diagrams (Figure 5(b) and 5(e)), to be presented in the next Section. Additionally 12 sources are thus included as members. In total, among the 200 members, 182 sources have reliable photometry and 18 sources have higher photometric uncertainty in any of the $J$, $H$, $K$, $[3.6]$, or $[4.5]$ bands.

Our selection does not completely reject field stars; some fraction of the candidates may well be field stars. To estimate the level of contamination, we choose the same control field as in the previous section, with the same sky area, and exercised the same set of selection criteria as in the cluster region. We obtained 16 sources in the control field that share similar kinematics. This gives a false positive rate of 8% in the candidate sample, indicating statistically 184 of the 200 sources to be true members of the cluster.

Figure 5 compares the color-magnitude diagram for stars seen toward the cluster region with that toward the control field region, with photometric data collected from Gaia DR2 and UKIDSS, and converted to absolute magnitudes by adopting the distance values computed by Bailer-Jones et al. (2018) from the Gaia parallax, which provides empirically better estimate of distances for all the stars with parallaxes published in the Gaia DR2, using a probabilistic inference approach, which takes into account for the nonlinearity of the transformation and the positivity constraint of distance.

The average age of the members is estimated by comparing with theoretical PARSEC isochrones (Bressan et al., 2012; Marigo et al., 2017; Pastorelli et al., 2019). The photometric sensitivity curves from Maíz Apellániz & Weiler (2018) are assumed for the Gaia sources, as they give empirically the best fit to our data. The majority of the members are consistent with an age between 5 to 30 Myr, with a best fit age of 15 Myr as seen in Figures 5(b) and 5(e). The estimated age is slightly younger than the values reported in literature (Section 1), where the age consistently varied from 20 Myr to 30 Myr. The candidate sample in this analysis is limited down to $G\sim 20.3$ mag, $G_{BP}\sim 21.1$ mag, $G_{RP}\sim 19.0$ mag, $J\sim 17.6$ mag, and $K\sim 16.7$ mag, equivalent to member masses $\gtrsim 0.6$ M_☉.

The luminosity function of the cluster is derived by statistically subtracting the number of stars in the control field from that toward the cluster region, for each magnitude bin. Field stars were chosen on the basis of satisfying the same proper motion and parallax criteria as for the member candidates. The decontaminated distribution of the luminosity function, generated with both the $G$ and $J$ bands is depicted in Figure 6(a). The luminosity function peaks around 14–15 mag for $G$ band and 13–14 mag for $J$ band, after which it declines swiftly owing to the data incompleteness.

The mass function of the cluster was derived on the basis of luminosity function and is presented in Figure 6(b). Member masses were estimated according to the PARSEC isochrones of an age of 15 Myr (Section 3.4), corrected for a distance of 1.22 kpc (Section 3.3) and reddening $E(B-V)=0.25$ mag (Sanner et al., 2000), with a polynomial interpolation of discrete data points. The slope of a power-law least-squares fit to the mass function, expressed as $\Gamma=d\log N(\log m)/d\log m$, where $N(\log m)$ is the number of stars per unit logarithmic mass, up to the photometric sensitivity, gives a slope of $\Gamma=-1.37\pm 0.18$ for the $G$ band, and a slightly steeper $\Gamma=-1.55\pm 0.14$ for the $J$ band. The mass range used to derive the mass function is $1.38\lesssim M/M_{\sun}\lesssim 7.11$ for the $G$ band and $1.28\lesssim M/M_{\sun}\lesssim 7.54$ for the $J$ band.

In comparison to literature works, Sanner et al. (2000) derived the cluster mass function with a slope of $\Gamma=-1.23\pm 0.17$ for the mass interval $0.72\lesssim M/M_{\sun}\lesssim 9.40$, on the basis of proper motion membership and statistical field star subtraction. Using the $V$ band photometry and a statistical field star subtraction approach, Sharma et al. (2008) reported a slope of $\Gamma=-1.80\pm 0.14$ for the mass range $1.01\lesssim M/M_{\sun}\lesssim 6.82$. Recently, Joshi et al. (2020) estimated $\Gamma=-1.26\pm 0.19$ for the mass range $0.72\lesssim M/M_{\sun}\lesssim 7.32$, by considering members only. Our results, based on reliable membership above $\sim 1$ solar mass selected with Gaia parallax and proper motions plus contamination subtraction, fall within the values published in the literature, and are consistent with the canonical value of $\Gamma=-1.35$ in the solar neighborhood (Salpeter, 1955).

Given that optical extinction decreases with increasing wavelength, observations made at longer wavelengths can detect more background stars through a cloud, and probe deeper cloud depths (Straw & Hyland, 1989; Dickman & Herbst, 1990). We utilized the $H$ and $K$ band photometry from the UKIDSS catalog, and constructed a number density image by defining a grid over the target area, following the method outlined in Gutermuth et al. (2005). Briefly, the region was divided into uniform grids of size $5\arcsec\times 5\arcsec$. The 20 nearest-neighbor sources from the center of each grid were selected to calculate the mean and standard deviation of the ($H-K$) color for each grid, excluding the sources for which the ($H-K$) values deviated $\gtrsim 3\sigma$ from the mean value. The mean ($H-K$) color for each grid then was converted to $A_{K}$, using the reddening law $A_{K}=1.82\times[(H-K)_{\rm obs}-(H-K)_{\rm int}]$, the difference between the observed and the intrinsic color (Flaherty et al., 2007).

From a nearby comparison field with little extinction, the average intrinsic color $(H-K)_{\rm int}$ was estimated to be $\sim 0.2$ mag. To restrict our analysis as far as possible to background objects, we have selected only sources with little extinction ($J-H<1.0$ mag, $H-K<0.6$ mag) for analysis (Panja et al., 2020). The resulting extinction map is displayed in Figure 7. The derived extinction values range from $A_{V}\simeq 1.33$–22.82 mag, or $A_{K}\simeq 0.12$–2.05 mag. Inspection of the extinction map reveals a compact region ($\sim 1\farcm 9\times 1\farcm 2$, centering around $\alpha_{2000}=84\fdg 023$, $\delta_{2000}=+34\fdg 103$) with excessive extinction ($A_{V}\sim 22.8$ mag). The rest of the cluster has otherwise relatively low extinction, with $A_{V}\lesssim 4.56$ mag.

An extinction map thus produced is limited in angular resolution by the availability of detectable background stars. Moreover, the extinction value along a particular line of sight is estimated in a statistical manner (Lada et al., 1994). Empirically, we found a $\sim 5\arcsec$ grid size, and $\sim 20$ nearest neighbor stars to be optimal choices, as a compromise between sensitivity and resolution. Our extinction map with an angular resolution of $5\arcsec$ and sensitivity down to $A_{V}\sim 27$ mag serves to guide the identification of heavily embedded sources such as protostellar objects, as will be discussed in the next section.

A certain combination of infrared colors can be used to distinguish YSOs at different evolutionary stages, as the amount of excessive infrared emission arising from the young circumstellar disks and infalling envelopes diminishes with age. We use such a color-color diagram to diagnose the possible nature of our targets. Lacking longer wavelength IRAC [5.8] and [8.0] $\mu$m detection, we utilized IRAC [3.6] and [4.5] $\mu$m along with UKIDSS $H$ and $K$ band photometry to characterize the YSOs (Dutta et al., 2018; Panja et al., 2020) toward M 36.

It has been shown that [3.6]$-$[4.5] is a useful YSO class discriminant color (Allen et al., 2004; Megeath et al., 2004), provided that the colors are dereddened. The dereddened ([[3.6]$-$[4.5]]₀ versus [$K-$[3.6]]₀) color-color diagram is displayed in Figure 8. To deredden a source, we estimated its line-of-sight extinction from the extinction map (Section 4.1) and exercised the reddening laws from Flaherty et al. (2007) to compute the dereddened colors. We used the set of equations developed by Gutermuth et al. (2009) to identify the YSOs based on different color criteria. To minimise contamination from extragalactic sources, an additional brightness cut was applied on the dereddened [3.6] $\mu$m magnitude by requiring a Class II object to have [3.6]${}_{0}<14.5$ mag and a Class I object to satisfy [3.6]${}_{0}<15$ mag (Gutermuth et al., 2009). While the majority of the field stars and main-sequence stars have dereddened colors about 0.0 mag to 0.2 mag, three Class I and one Class II candidates stand out, with three of them in close physical association with, signifying ongoing star formation in, the dusty clump depicted in Figure 7. The photometric parameters of these YSOs are listed in Table 2.

To reveal the underlying stellar and circumstellar properties of the young objects, we fitted their spectral energy distributions (SEDs) using theoretical models (Zhang et al., 2015). We used the radiative transfer models of Robitaille et al. (2006, 2007) to derive the physical parameters of the YSOs. The grids of models cover a wide range of parameters, consisting of pre-main-sequence stars plus a combination of axisymmetric circumstellar disks, infalling flattened envelopes, and outflow cavities. These models also span a large range of possible evolutionary stages (from deeply embedded protostars to stars surrounded only by optically thin disks) and stellar masses (from 0.1 to 50 $M_{\sun}$). The radiative transfer solution is subject to parameter degeneracy, but fluxes covering a sufficiently wide range of wavelengths can reduce this degeneracy.

The photometric fluxes of the YSOs are collected from NOMAD ($V$; Zacharias et al. 2004), UKIDSS ($J$, $H$, and $K$; Lawrence et al. 2007), IRAC ([3.6] and [4.5] $\mu$m; Fazio et al. 2004), MSX (8.28, 14.65, and 21.34 $\mu$m; Sjouwerman et al. 2009), and IRAS (12, 25, 60, and 100 $\mu$m; Neugebauer et al. 1984), wherever available, detailed in Table 3. We required at least a minimum of five data points for each source to construct the SEDs, which then are fitted with the distance and visual extinction as input parameters to constrain the SEDs. We varied the distance as $1.20\pm 0.13$ kpc (Section 3.2) and extinction $A_{V}\sim 1.3$ mag to 23 mag (Section 4.1). For each source, the “best-fit” model is defined by constraining $\chi^{2}-\chi^{2}_{\rm best}<3$ (per data point), where $\chi^{2}_{\rm best}$ quantifies the goodness-of-fit parameter. The SEDs and their best-fit models of the four YSOs are shown in Figure 9. Expecting that the fitted SEDs are moderately degenerate, we adopted the weighted mean and standard deviations of each physical parameter, and they are given in Table 4. The degeneracies obtained for sources S1, S2, S3, and S4 are 36, 6, 2, and 9, respectively.

The SEDs of the YSOs suggest stellar infancy, with ages less than 0.2 Myr. The stellar masses range from $\sim 0.1M_{\sun}$ to $6.5M_{\sun}$, radii from $\sim 2.5R_{\sun}$ to $24.6R_{\sun}$, and temperatures from $\sim 2830$ K to 6260 K. The envelope and disk parameters serve to probe the evolutionary phase. Class I objects typify earlier stages of star formation compared to Class II objects. At the initial phases, the envelope and disk accretion rates are usually high, and then decrease with time. It is expected that envelopes with accretion rates $\lesssim 10^{-6}M_{\sun}$ yr^-1 do not contribute significantly to the SED (Robitaille et al., 2007). Among the four objects, S4 has the lowest envelope accretion rate ($\sim 10^{-6}M_{\sun}$ yr^-1), the lowest disk mass ($\sim 10^{-4}M_{\sun}$), and the lowest disk accretion rate ($\sim 10^{-9}M_{\sun}$ yr^-1), all consistent with our earlier results of S4 being a Class II object. The contribution to the SED from the disk related to that from the envelope is difficult to disentangle, because the envelope may not be largely dispersed yet. The YSO S3 is associated with IRAS 05327$+$3404 (nickname “Holoea”⁹⁹9In Hawaiian for “flowing gas”), first detected by Magnier et al. (1996) and classified as being transitional between Class I and Class II, and in the process of becoming optically exposed (Magnier et al., 1999a, b). This object, with strong far-infrared fluxes, has a very prominent circumstellar disk, extending from an inner radius of $\sim 25$ AU to an outer radius of $\sim 690$ AU.

Despite being the most abundant molecular species in cold interstellar media, H₂ molecules are difficult to detect due to a lack of a permanent electric dipole moment. CO lines are therefore often used to trace molecular clouds. Using $J$=1–0 ¹²CO and ¹³CO data, we map the molecular cloud structures and estimate their kinematic and physical properties. A compact and massive core is detected around the position $\alpha_{2000}$ = $84\fdg 0308$, $\delta_{2000}$ = $+34\fdg 1094$ for which we investigated the intensity, excitation temperature, H₂ column density, and velocity distribution. The complete procedures and expressions used to derive these parameters are detailed in Sun et al. (2020).

Of the 200 member candidates of M 36, 16 (8%) were observed as LAMOST DR5 targets. RV measurements were not included in our membership selection, because of the limited amounts of data available, and because of possible variations due to binary orbital motion. The RV and metallicity measurements are presented in Figure 13. The RVs of the majority of members range between $-5$ km s^-1 and $-20$ km s^-1. This is consistent with that reported by Frinchaboy & Majewski (2008, RV$=-17.83\pm 0.99$ km s^-1), and matches approximately with the RV ($\approx-21$ km s^-1) of the molecular cloud.

The distance to the cloud is uncertain. Adopting RV$=-21$ km s^-1, the nominal Milky Way rotation (Reid et al., 2019) suggests a kinematic distance of $1.68\pm 0.05$ kpc, placing the cloud somewhat in the background but physical association of the cloud with the cluster ($\sim 1.2$ kpc) cannot be ruled out, in evidence notably of the Gaia distance of $1.395_{-0.26}^{+0.40}$ kpc (Table 2) for the young object S2 in the cloud. While not unambiguously established at this time, the proximate angular separation, radial velocity, and distance provides tantalizing evidence of a physical connection between the cloud and the cluster.

The metallicity distribution of member candidates peaks around $-0.1$ to $-0.2$, with an average of $-0.15\pm 0.10$. Though the dispersion is relatively large, our result, based on a relatively reliable list of members, suggests a subsolar metallicity for M 36.

Star formation activity at the earliest stages (1–2 Myr) is difficult to trace due to heavy dust obscuration. Star formation is believed to be regulated by dense gas in molecular clouds (Gao & Solomon, 2004; Wu et al., 2005; Dutta et al., 2018), and the star formation rate is known to be correlated with the molecular mass (Wong & Blitz, 2002; Lada et al., 2010, 2012). Surveys of star-forming regions have demonstrated that approximately 75% of the stars are formed in groups or clusters (Carpenter, 2000; Allen et al., 2007; Gutermuth et al., 2009), whereas about 80% of all young stars are located in embedded clusters (Lada & Lada, 2003; Porras et al., 2003). Several factors are involved in the formation and early evolution of a cluster, such as the structure of the parental molecular cloud (Samal et al., 2015), fragmentation in the parental cloud due to turbulent motions and/or gravity, dynamical motions of the young stars, and the feedback from young stars (Gutermuth et al., 2009). While infrared surveys are an eminent tool to parameterize young embedded clusters, data at much longer millimeter wavelengths are also powerful to trace the gaseous contents and the distribution of dust throughout a molecular cloud.

A color composite image of the M 36 cluster, made from optical ($B$ band), near-infrared ($K$ band), and far-infrared (WISE/W3 12 µm) data, is presented as Figure 14(a). The WISE 12 $\mu$m emission, sensitive to warm dust (Dunham et al., 2008), coincides well with the extinction complex and CO gas. The ¹²CO and ¹³CO distributions reveal a compact cloud with a linear extent of $\sim 10\farcm 4$ (3.6 pc), extending from north-east to south-west plus a core located at the south-west corner where the three YSOs discussed earlier reside near the densest part, within $\sim 1\arcmin$, of the cloud core. This is where active star formation is taking place.

The timescale of star formation plays a significant role in the evolution of a star cluster. It depends not only on the initial cloud conditions, such as the density and temperature profiles, but also on how the new-born stars interplay with the remnant clouds, so as to induce the birth of next-generation stars, e.g., by the dynamical compression of an H II front, and change in self-gravity (Fukuda & Hanawa, 2000). Alternative to clustering, star formation may proceed in an isolated or distributed mode (Koenig et al., 2008; Evans et al., 2009), leading to multiple stellar populations within a cloud (Gratton et al., 2012; Bouy et al., 2015; Bekki & Tsujimoto, 2017).

The four YSOs identified by this work are in the proximity of the cloud core, with angular separations of $3\farcm 69$ (for S1), $0\farcm 79$ (for S2), $0\farcm 46$ (for S3), and $0\farcm 15$ (for S4). The source S2 has a measured Gaia DR2 distance of $1.395_{-0.26}^{+0.40}$ kpc (Table 2), matching well with the cluster distance. In particular, S3 coincides spatially with an IRAS source (IRAS 05327$+$3404), which drives a powerful ionizing outflow with a velocity of $\sim 650$ km s^-1, that aligns with a CO outflow (Magnier et al., 1996). Its SED indicates stellar infancy with ample circumstellar material. S3 is clearly a part of the continuing stellar formation episode in the dense core, and its outflows may have a profound effect on the surrounding environments to prompt future formation of stars.

Star formation processes in an open cluster may span about 5 Myr to 20 Myr (Lim et al., 2016), with the low-mass members spending up to some $\sim 10^{8}$ years in the pre-main sequence phase (Herbig, 1962). Star clusters older than $\sim 5$ Myr are typified with little parental molecular gas (Leisawitz, Bash & Thaddeus, 1989). M 36, with an age of $\sim 15$ Myr (Section 3.4), is therefore unusual in connection with a nearby voluminous molecular cloud harboring a protostellar population with ages $<0.2$ Myr. More evidence is clearly needed to confirm or to discredit the physical connection between the star-forming cloud and the cluster. If the cloud and the cluster share the same formation scenario, this presents a rare case of sustaining star formation in a cloud complex. Even if the cloud happens to be in the intermediate background of the cluster, our discovery of a cloud active in star formation offers added information of the star formation history in this part of the Galactic disk.