Molecular Bubble and Outflow in S Mon Revealed by Multiband Data Sets

We identify a molecular bubble, and study the star formation and its feedback in the S Mon region, using multiple molecular lines, young stellar objects (YSOs), and infrared data. We revisit the distance to S Mon, ∼722 ± 9 pc, using Gaia Data Release 3 parallaxes of the associated Class II YSOs. The bubble may be mainly driven by a massive binary system (namely 15 Mon), the primary of which is an O7V-type star. An outflow is detected in the shell of the bubble, suggesting ongoing star formation activities in the vicinity of the bubble. The total wind energy of the massive binary star is 3 orders of magnitude higher than the sum of the observed turbulent energy in the molecular gas and the kinetic energy of the bubble, indicating that stellar winds help to maintain the turbulence in the S Mon region and drive the bubble. We conclude that the stellar winds of massive stars have an impact on their surrounding environment.


INTRODUCTION
The mass-loss phase is a very common phenomenon during the early evolutionary stages of a star.Strong winds from young stars inject momentum and energy into the surrounding environment, thus affecting the dynamics and structure of their parent clouds (e.g., Lada 1985;Arce et al. 2011;Frank et al. 2014;Bally 2016;Li et al. 2020).High-mass stars evolve very quickly and they may still be embedded in their parent cloud when they reach the main-sequence stage.When a high-mass star forms in the molecular cloud, its ultraviolet (UV) radiation can ionize and heat the surrounding gas to create an H II region.The expanding H II region can reshape the surrounding molecular gas and drive a bubble (e.g., Churchwell et al. 2006;Zhang et al. 2015).During the main-sequence stage, high-mass stars can drive spherical winds that blow away the gas around the star and then form bubbles (e.g., Castor et al. 1975;Arce et al. 2011).Bubbles and outflows are common stellar-feedback phenomena in the star-formation process and provide information about the physical properties of their surroundings (e.g., Lada 1985;Arce et al. 2010Arce et al. , 2011;;Li et al. 2018Li et al. , 2020)).
There have been quite a few studies concerned with bubbles, outflows and their feedback (e.g., see Frank et al. 2014;Dale 2015; Bally 2016, for reviews).Since directly measuring the proper motion of extended structures, like molecular clouds, is very difficult, many previous studies only analyzed the line-of-sight dynamics in a region.Additionally, recent studies started to investigate the 3D structure of molecular clouds in the solar neighborhood (e.g., Green et al. 2019;Lallement et al. 2019;Guo et al. 2021), mainly enabled by the Gaia space mission (Gaia Collaboration et al. 2016).Young stellar objects (YSOs) that only recently formed in star-forming molecular clouds are still close to their birth sites; hence, they still, on average, share the distance and kinematic properties of their parent clouds (e.g., Fűrész et al. 2008;Tobin et al. 2009;Hacar et al. 2016;Großschedl et al. 2021).On this basis, cloud distances can be estimated by using the distances of YSOs located within the clouds as proxies (e.g., Großschedl et al. 2018Großschedl et al. , 2021;;Zhang 2023).The recent release of astrometric data from Gaia Data Release 3 (DR3, Gaia Collaboration et al. 2023) allows one to obtain a large number of high-quality parallaxes and proper motions of nearby YSOs, which also makes it possible to explore further the 3D kinematics of molecular clouds that harbor these YSOs (e.g., Großschedl et al. 2021;Flaccomio et al. 2023).
NGC 2264 is a nearby high-mass star-forming complex in Monoceros with a distance of about 700-800 pc (e.g., Zucker et al. 2020;Flaccomio et al. 2023), and it is well studied at multiple wavelengths (for a review see Dahm 2008).Studies of the spatial and dynamical structure of NGC 2264 revealed that it is composed of two complexes: S Mon in the northern part and the Cone Nebula in the southern part (e.g., Sung et al. 2008;Tobin et al. 2015;Venuti et al. 2018).The gas structure of S Mon is dominated by massive stars at its center.For instance, 15 Mon, which consists of an O7V primary star with a mass of 35 M ⊙ and an O9.5V secondary star with a mass of 24 M ⊙ (Dahm 2008).15 Mon exhibits variations in its UV resonance line profile as well as fluctuations in its soft X-ray flux, which can be interpreted as being induced by variations in its mass-loss rate (Snow et al. 1981;Grady et al. 1984;Dahm 2008).Sung & Bessell (2010) and Venuti et al. (2018) studied the star formation history of S Mon by analyzing the properties of young clusters, and found that stars are still forming actively in the region.Buckle et al. (2012) found that the CO emission surrounding S Mon shows many filaments and arcs, and these structures are unlikely due to protostellar outflow activity.Based on 13 CO emission, Tobin et al. (2015) speculated that there exists a bubble in the S Mon region which may be driven by 15 Mon.Since studies of the S Mon region either used only young stars or molecular gas, or did not focus on the details of the bubble, the properties of this region, especially the bubble, outflows, and their feedback remain unclear.
In this work, we carried out molecular line observations and collected YSO information in the S Mon region, aiming to provide insights into the star formation activities and their feedback.This paper is organized as follows.In Section 2, we present the data used in this work.The correlation between YSOs and molecular gas, the distance of the molecular gas, and details of the bubble and outflow are presented in Section 3. We discuss the driving source of the bubble and the energy cascade in S Mon, as well as the kinematic features of the bubble based on the YSOs in Section 4. In Section 5, we summarize the main results.

Molecular Gas
Molecular line observations of S Mon were carried out between 2022 April and November using the Purple Mountain Observatory Delingha (PMODLH) 14 m millimeter-wavelength telescope.Four molecular lines were observed: 12 CO (J = 1 -0), 13 CO (J = 1 -0), C 18 O (J = 1 -0), and HCO + (J = 1 -0).The three CO lines were observed simultaneously.All lines were observed with the nine-beam Superconducting Spectroscopic Array Receiver system in the sideband separation mode (Shan et al. 2012), with the 13 CO and C 18 O lines in the lower sideband and the 12 CO and HCO + lines in the upper sideband.The on-the-fly (OTF; Sun et al. 2018) mode was used during the observations, and the OTF raw data were gridded in fits cube with a pixel size of 30 ′′ using the GILDAS software package (Pety 2005). 1 The typical integration time per position was approximately 3 minutes.Each fast Fourier transform spectrometer with a bandwidth of 1 GHz provides 16 384 channels, yielding a spectral resolution of 61 kHz. 2 All results presented in this work are expressed as brightness temperatures, T * R = T * A /η mb , where T * A is the antenna temperature and η mb is the main beam efficiency.Table 1 lists the observational parameters of the molecular lines.

Young Stellar Object Sample
The classification of YSOs was initially proposed by Lada & Wilking (1984) to delineate their evolutionary stages, and is defined as follows.Class I YSOs are protostars that are still embedded in an envelope and surrounded by a circumstellar disk, causing significant infrared-excess and a rising slope of the spectral energy distribution (SED) in the mid-to far-infrared range; Class II YSOs are pre-main-sequence (PMS) stars with circumstellar disks, which show a negative SED slope, that is significantly above the slopes of main-sequence stars (see Lada et al. 2006).Class II YSOs   with inner holes in their accretion disks are named transition disks and Class III YSOs are PMS stars surrounded by thin disk remnants (anaemic disks, ADs) or they could be already disk free.
To study the young cluster NGC 2264, a variety of multiwavelength datasets and criteria have been employed to identify the evolutionary stages of YSOs.For example, Sung et al. (2008Sung et al. ( , 2009) ) classified YSOs based on SED slopes of infrared data combined with X-ray data from the Chandra X-ray Observatory (Weisskopf et al. 2002).Broos et al. (2013) adopted a Naive Bayes approach to classify YSOs, using data from the Massive Young star-forming Complex Study in Infrared and X-rays (MYStIX) project (Feigelson et al. 2013).Cody et al. (2014) classified YSOs based on the slope of SEDs at near-and mid-infrared wavelengths.Rapson et al. (2014) utilized a color-based classification scheme to classify YSOs with Spitzer (Fazio et al. 2004) data and Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006) photometry.Venuti et al. (2018) classified YSOs using spectroscopic and photometric data from the Gaia-ESO Survey (GES, Gilmore et al. 2012;Randich et al. 2013) and the Coordinated Synoptic Investigation of NGC 2264 Survey (CSI 2264, Cody et al. 2014).
We describe in detail the YSO data collection, classification criteria, and astrometric parameters acquisition in Appendix A, and briefly introduce these procedures here.We collected YSO candidates in the S Mon region from Sung et al. (2009), Broos et al. (2013), Cody et al. (2014) and Rapson et al. (2014), within an area of 202.8 • < l < 203.1 • , and 2.0 • < b < 2.3 • .Duplicates were eliminated through coordinate cross-matching with a radius of 1 ′′ , yielding in an initial sample of 1 023 YSO candidates.As these candidates were classified by different methods, we reclassified them through a homogeneous classification scheme based on infrared data.We adopted the SED slope (α = d log λF λ /d log λ, Lada 1987) and the classification scheme of Greene et al. (1994) to classify the evolutionary stage of the YSOs.Flatspectrum sources are classified as Class II YSOs in this study.In line with Sung et al. (2009), we use α ≥ 0.3 to classify objects as Class I and −1.8 to 0.3 for Class II.When using selection criteria based on infrared photometry, only the part of Class III YSOs with infrared-excess (i.e., ADs) can be classified.Therefore, X-ray data collected by Flaccomio et al. (2023) were also employed to classify additional Class III YSOs.The Class III YSOs in this work are comprised of sources with ADs and disk-less PMS stars.For the former, the lower limit for a source with infrared-excess is given by Lada et al. (2006) with α > -2.56.This slope is typical for stellar photospheres when there is no additional infrared-excess from a disk present, and when using the Spitzer IRAC bands to calculate the infrared SED (see e.g., Lada et al. 2006;Teixeira et al. 2012).Given the uncertainties in photometry, however, we took -2.3 as a rough limit to exclude possible contaminants from field stars.For disk-less PMSs, they have a SED slope of below -2.3 and X-ray emission.The SED fitting was performed without reddening corrections following Cody et al. (2014).In this case, we gain a sample of seven Class I, 160 Class II, and 245 Class III YSOs (101 ADs + 144 disk-less PMSs).The parallaxes and proper motions of these YSOs were taken from Gaia DR3 (Gaia Collaboration et al. 2023).The radial velocities (RVs) of the YSOs were collected from the Sloan Digital Sky Survey (SDSS) DR17 APOGEE-2 (Abdurro'uf et al. 2022) and the catalog of Tobin et al. (2015).
To summarize, we end up with a sample of one Class I, 145 Class II, and 231 Class III YSOs that had matched with the Gaia DR3 catalog, and 47 Class II and 91 Class III YSOs with available RV information.
To compare the gas and stellar RVs, the heliocentric RVs (V Helio ) of the YSOs were converted to the local standard of rest velocity (V LSR ).We followed the conversion method described by Reid et al. (2009) via: where l and b are the Galactic longitude and latitude, respectively, and U ⊙ , V ⊙ , and W ⊙ denote the solar motion components in the direction of the Galactic Center, along Galactic rotation, and toward the north Galactic pole, respectively.The used values for U ⊙ , V ⊙ , and W ⊙ are 10.3, 15.3, and 7.7 km s −1 , respectively (Kerr & Lynden-Bell 1986).

RESULTS
In the following we present the relationship between the molecular gas and YSOs in Section 3.1.Then, the distance to the S Mon region is estimated in Section 3.2 using YSOs that are closely correlated to the gas.In Section 3.3, we present the bubble observed in the S Mon region.An outflow is detected in one of clumps, as described in Section 3.4.

Relationship between the Molecular Gas and YSOs
Figure 1 shows a position-velocity (P-V) diagram (a) and the normalized V LSR distribution (b) for YSOs and molecular gas traced by 12 CO in S Mon.The 12 CO emission shows two velocity components in this region, with one ranging from 0 to 7 km s −1 and the other from 7 to 15 km s −1 .The presence of a gap or minimum in the CO emission at V LSR ∼ 7 km s −1 is in agreement with the results of Tobin et al. (2015).In panel (b) of Figure 1, the V LSR of the majority (62%) of the Class II YSOs peak at about 10 km s −1 , which is consistent with the 12 CO emission component ranging from 7 to 15 km s −1 .Most (71%) Class III YSOs have LSR velocities between 0 and 7 km s −1 , with a peak at 5 km s −1 .
We investigate the two velocity components of molecular gas in S Mon separately.Figure 2 shows integrated maps of 12 CO and HCO + within the velocity range from 0 to 7 km s −1 (panels (a) and (b)) and from 7 to 15 km s −1 (panels (c) and (d)).The integrated map of molecular gas is superposed on a Herschel (Pilbratt et al. 2010) H 2 column density map of this region, with a resolution of 18 ′′ (Nony et al. 2021).The 12 CO and HCO + emission between 0 and 7 km s −1 appears to be diffuse and is not consistent with the Herschel H 2 column density map.This component seems to dominate in the easternmost area of the region.In particular in the 12 CO map (Figure 2 (a)) there appears to be a somewhat stronger component starting from the left of the image, beyond the shown sky cutout.Our investigation does not include the easternmost region.The 12 CO and HCO + emission between 7 and 15 km s −1 shows an arc-like structure, possibly indicative of a "bubble" (centered at (0, 0) offset, see more details in Section 3.3).We note that the distribution of the molecular gas within this velocity range is in agreement with the distribution of the Herschel H 2 column density map.
Figure 3 depicts the projected distribution of the YSOs in combination with 12 CO emission and Herschel H 2 column density map.Both the Class II and Class III YSOs, which have LSR velocities between 7 and 15 km s −1 (yellow pluses), generally coincide with the denser regions of the molecular cloud.More members associated with the molecular gas are found belong to the Class II YSOs.
Compared with the Class III YSOs, the Class II YSOs are better correlated with the molecular gas in both spatial distribution and kinematic properties.Considering the association between the molecular gas and the Class II YSOs, we are able to use the Class II YSOs to characterize the properties of the molecular gas, such as its distance and kinematics.The proper motions of the YSOs are used to derive the kinematics of the molecular gas in the plane of the sky, since they are the missing part, and cannot be measured directly for clouds.

Distance to S Mon
Accurate distances to molecular clouds are essential to study their physical properties.We employ the Gaia DR3 parallaxes of the Class II YSOs to estimate the distance to S Mon.In this work, only Class II YSOs with parallax accuracies better than 10% are selected for distance estimation.First, outliers of parallaxes beyond 3σ from the median value are excluded.Class II YSOs, which do not coincide with the gas distribution (CO at 7-15 km s −1 , magenta pluses in Figure 4 (a)), are also excluded from the distance estimation.Finally, a total 35 YSOs that meet these criteria are selected to estimate the distances shown (yellow pluses in Figure 4 (a)).
The distance is determined from Gaia parallaxes with d(pc) = 1 000/ϖ, where ϖ is the parallax in units of mas.Lindegren et al. (2021) provided a method to correct the parallax zero points, however for red sources with m GBP − m GRP > 1.6, the correction could lead to a large parallax bias of about 40 µas.Since the photometric data of the used Class II YSOs have generally redder colors, we do not apply the parallax zero-point correction in this work.
We estimate the distance to the cloud directly using the YSOs associated with the cloud and their Gaia DR3 parallaxes.The cumulative distribution function (CDF) of Gaia distances is shown in panel (b) of Figure 4, where the median distance estimate is of 722 ± 9 pc.The uncertainty is derived by σ/ √ N , where σ represents the dispersion of the distance, and N represents the number of stars.
Flaccomio et al. ( 2023) conducted a study on the NGC 2264 star-forming region using X-ray data obtained with XMM-Newton telescope.The S Mon core region in their study appears to be close to the region we studied, although we note that they did not give the specific coordinates of the region they investigated.Although only half of the samples we used for distance estimation aligned with their sample counterparts, our distance estimation is similar to their results (719 ± 9 pc).However, the distance we derive is smaller than the distance estimated (759 +10 −26 ± 37 pc) by Zucker et al. (2020), who inferred the distance by fitting a line-of-sight dust model with Gaia DR2 data.We conducted a parallax difference test between Gaia DR2 and DR3 using a selected set of Class II YSOs in our region.These YSOs have LSR velocities ranging from 0 to 15 km s −1 , and both DR2 and DR3 exhibit parallax accuracies better than 10%.The difference between these two groups is approximately 25 pc, meaning the two estimates agree within their given uncertainties.Although we did not test the method of Zucker et al. (2020) explicitly, variations in methods could also introduce some differences in distance estimations.Therefore, the difference in distance between this work and that of Zucker et al. (2020) could be caused by the difference in both method and data.In the subsequent analysis, we adopt 722 ± 9 pc as the distance to the S Mon cloud, which has a velocity ranging between 7 and 15 km s −1 .

Bubble
We identified the bubble in S Mon based on four criteria summarized from previous studies (e.g., Arce et al. 2011;Li et al. 2015;Feddersen et al. 2018): (1) the shell shows a circular or arc-like structure in integrated CO or HCO + maps; (2) the P-V diagram of the shell shows an expansion signature (i.e., a "U" or "V" shape); (3) the CO or HCO + shell matches the morphology of infrared data in at least one band; and (4) the shell contains a candidate driving source.An expanding bubble model requires three parameters: bubble radius (R bub ), expansion velocity (V exp ), and central velocity (V cnt ).We followed the method of Arce et al. (2011) to estimate these parameters.The radius is determined by a Gaussian fit to the azimuthally averaged profile of the integrated intensity map.The central and expansion velocities can be roughly read from a P-V diagram, where the upper part of the "U" or "V" structure corresponds to the central velocity and the difference between the central velocity and the lowest point of "U" or "V" structure indicates the expansion velocity.Then, the central and expansion velocities are corrected by visual fitting based on the channel maps.As indicated by Arce et al. (2011), the uncertainty in the central velocity of a bubble is about ±0.5 km s −1 , and since the bubble is not detected over its entire velocity range, the expansion velocity is obtained as a lower limit.
An arc-like structure is clearly visible in the channel maps of 12 CO and HCO + (see the details in Figures 5-6), spanning a broad V LSR range from 7.5 to 12.5 km s −1 .This structure is also evident in the integrated maps of 12 CO and HCO + from 7.5 to 12.5 km s −1 and the Herschel H 2 column density map (see panels (a) and (d) of Figure 7, Nony et al. 2021).The P-V diagrams of 12 CO and HCO + (panels (b) and (e) of Figure 7) both show an upside down "U" or "V" shape.Based on the P-V diagram of 12 CO (HCO + ), V exp and V cnt are derived as ∼ 3.1 km s −1 (∼ 3.4 km s −1 ) and ∼ 8.2 km s −1 (∼ 8.4 km s −1 ), respectively.The V exp values estimated from 12 CO and HCO + are similar, thus an average of 3.3 km s −1 is adopted as a lower limit for the expansion velocity.The center of the bubble is set to be (0, 0).The radius and thickness of the bubble traced by the 12 CO emission (see panel (c) of Figure 7) are ∼ 5.4 and ∼ 4.2 arcmin, respectively, corresponding to 1.13 and 0.88 pc at a distance of 722 pc.These values are compared with the bubble traced by HCO + emission (see panel (f) of Figure 7), i.e., a radius of ∼ 5.3 arcmin (∼ 1.11 pc) and a thickness of ∼ 4.4 arcmin (∼ 0.92 pc).Based on both the 12 CO and HCO + emission and considering the uncertainties from the distance and the Gaussian fits, we derive the radius and thickness of the bubble to be 1.1 ± 0.1 pc and 0.9 ± 0.1 pc, respectively.The mean column density of the bubble is ∼ 2.5 × 10 20 cm −2 and the mass (M shell ) is ∼ 1 400 ± 400 M ⊙ (see the calculation and estimation of uncertainties in Appendix B).
Another arc-like structure appears to be located in the southwest region of the Herschel H 2 column density map (see the boxed region of Figure 8 (a)).The low-sensitivity CO data, observed using PMODLH, were employed to investigate the structure.The rms noise of the 12 CO, 13 CO and C 18 O spectra is 0.37, 0.20, and 0.20 K, respectively.A small arc-like structure is only evident in the 13 CO integrated map (see panel (b) of Figure 8), so we use the 13 CO data to produce a P-V diagram for the arc-like structure.However, there is no clear "U" or "V" shape in the P-V diagram (see panel (c) of Figure 8).In addition, we do not find a potential driving source for this arc-like structure.Moreover, using high-spatial-resolution CO maps (20 ′′ ), Tauber et al. (1993) found that the properties of the northern region of the arc-like structure (∆b > −7 ′ ) are consistent with the ionizing radiation from 15 Mon, while those of the southern region are not.In conclusion, the arc-like structure shown is unlikely to be a bubble.Given the wind-blown appearance, this feature might be caused by feedback from the northeast side, coincident with the general location of the massive O-and B-type stars in the region.
There are eight clumps, labeled as A-H in Figure 7, identified from the integrated map of C 18 O emission.Figure 9 plots the four molecular lines we observed toward the C 18 O peak position of each clump.Table 2 lists the fits for the C 18 O peak emission for each clump.There are three main components among these eight clumps, i.e., 3 km s −1 < V LSR < 7 km s −1 , 7 km s −1 < V LSR < 10 km s −1 , and 10 km s −1 < V LSR < 13 km s −1 .As discussed in Section 3.1, the bubble is not associated with the component of 3 km s −1 < V LSR < 7 km s −1 , but instead with the other two components.Clumps A to D correlate with the velocity component of 7 km s −1 < V LSR < 10 km s −1 and centered at ∼ 8.5 km s −1 , and clumps E to H correlate with the velocity component of 10 km s −1 < V LSR < 13 km s −1 and centered at ∼10.5 km s −1 .Figure 10 shows integrated maps of C 18 O with different velocity components, which are in agreement with the parameters given in Table 2.These two groups of clumps correspond to different lobes of the bubble, so the different kinematic characteristics may be caused by the expansion of the bubble.Furthermore, the velocity dispersion around the bubble is larger (see Figure 11), indicating possible interaction between the bubble and surrounding gas.Table 3 lists the H 2 column density, mass, and turbulent energy of each clump dervied using the optically thin 13 CO line at the peak of the clump (see the calculations in Appendix B).Note-(1) Source name.( 2), ( 5), ( 8), (11) Brightness temperature(s) of 12 CO, 13 CO, C 18 O, and HCO + , respectively, at the emission peak of C 18 O.(3), ( 6), ( 9), (12) FWHM of 12 CO, 13 CO, C 18 O, and HCO + , respectively.( 4), ( 7), ( 10), (13) Central velocity of 12 CO, 13 CO, C 18 O, and HCO + , respectively.

The outflow of clump E
In this section we investigate clump E in more detail, since the velocity profile of this clump alone displays typical outflow characteristics and because the spectra of the other clumps are too complex to identify reliable outflows.Outflows are a direct signature of ongoing star formation, and they can be identified by examining their line profiles, integrated intensity maps, and P-V diagrams.The central velocity and position of the driving source are estimated through the C 18 O line.The emission peaks of C 18 O can be seen in Figure 7.The initial velocity ranges of the blue and red wings are determined from the line wings of the spectral line where the C 18 O emission reaches the 1σ noise level.According to morphologies of the P-V diagram and the integrated map of the line wings of the 13 CO and HCO + emission, then visually adjusted the velocity range of the line wings.To improve the signal-to-noise ratio, we smooth the CO and HCO + lines to a velocity resolution of 0.25 km s −1 .
Figures 12 and 13 show the identified outflow in clump E based on 13 CO and HCO + emission, respectively.The spatial distributions of the 13 CO and HCO + outflows are consistent with each other and are aligned along the east-west direction.The blue and red lobes are symmetrically distributed with respect to the C 18 O emission peak, which are similar in size.We note that the C 18 O and HCO + emission at V LSR ∼ 8.6 km s −1 , and the offset of ∼ 1.5 arcmin in the P-V diagram, are not associated with the blue lobe of the outflow, thus we did not cut the edge of the blue wing to the 1σ noise level.
The method presented in Liu et al. (2021) was employed to estimate various properties of the outflows, such as the mass, momentum, and energy.Therein, they derived the H 2 column density of the outflow lobes assuming that the gas is in local thermodynamic equilibrium (LTE) and the excitation temperature is 30 K (15 K) for 13 CO (HCO + ).The mass of an outflow lobe can be then estimated based on its size and the H 2 column density.The other physical properties of the outflows are obtained based on the mass, velocity, and size of the lobes (see Appendix A of Liu et al. 2021).In our work, we use the ratios of column densities of N H2 /N13 CO ≈ 5×10 5 and N H2 /N HCO + ≈ 10 8 (Turner et al. 1997;Simon et al. 2001).We also use the inclination of the outflow, i.e., the angle between the long axis of the outflow and the line of sight, to correct the parameters of the outflow.Since the inclination of an outflow cannot be reliably determined, we adopt 57.• 3 as a proxy in this work, which is the mean value for a random distribution (see details in Bontemps et al. 1996;Li et al. 2018).We note that the calculation of the outflow's physical parameters (see Table 4), involve many assumptions, such as the excitation temperature and the inclination.Therefore, we only consider the uncertainties in the outflow lobe region and distance when estimating the errors of the physical parameters.The mass, momentum, and dynamical timescale of the outflow are a few M ⊙ , a few tens of M ⊙ km s −1 , and ∼ 10 5 yr, respectively.Based on Yang et al. (2018), the typical mass of a low-mass outflow is about 0.1 ∼ 1 M ⊙ and the typical mass loss rate is about 10 −7 ∼ 10 −6 M ⊙ yr −1 .These values are much smaller than those found for the outflow in clump E.Moreover, all physical parameters of the outflow in clump E align with the typical values of high-mass outflows as summarized in table 1 of Yang et al. (2018).Therefore, it is unlikely that the outflow in clump E is driven by a single low-mass YSO.We do not find any high-mass YSOs in this region based on the current data, yet there are several low-mass YSOs embedded in the vicinity of the outflow region (see panel (a) of Figures 12 and 13).Therefore, the outflow might be driven by a deeply embedded protostar that has not been detected yet.The contour levels are shown from 40% to 90% with steps of 10% of the peak intensity of each lobe.The black diamonds and cyan pluses represent Class I and Class II YSOs, respectively.The magenta, blue, and red stars are the emission peaks of C 18 O, and the blue lobe and red lobe, respectively.(b) P-V diagram along the white arrow in panel (a).The black contour levels are shown from 10% to 90% with steps of 10% of the peak value.The blue and red dashed lines are the velocity ranges of the blue and red lobes, respectively.(c) The blue spectrum presents 13 CO at the peak position of the blue emission of the 13 CO outflow.The black spectrum is C 18 O at the emission peak of C 18 O.The blue shading of the spectrum indicates the blue line wing velocity of 13 CO.The feature at 8.6 km s −1 is likely unrelated to the outflow, so emission with a velocity below 8.6 km s −1 is not included in the blue wing of the outflow profile.(d) The red spectrum is 13 CO at the peak position of the red emission of the 13 CO outflow.The black spectrum is C 18 O at the emission peak of C 18 O.The red shading of the spectrum indicates the red line wing velocity of 13 CO.

DISCUSSION
In this section, we discuss the physical properties of the bubble (Section 4.1), and the energy cascade of S Mon region (Section 4.2).In Section 4.3, we use the Class II YSOs to investigate the kinematic characteristics of the bubble.

Physical Properties of the Bubble
In the following, we discuss the physical properties of the bubble.Using the shell mass (M shell ) and expansion velocity (V exp ) of the bubble, we can estimate the shell's momentum (P shell ), kinetic energy (E shell ), and kinetic timescale (t kinetic ).These quantities are determined as P shell = M shell V exp , E shell = 0.5M shell V 2 exp , and t kinetic = R bub /V exp , respectively.For the bubble we find that P shell is 4 600 ± 1 300 M ⊙ km s −1 , E shell is (1.5 ± 0.4) × 10 47 erg, and t kinetic is (3.3 ± 0.3) × 10 5 yr.Considering that the expansion velocity serves as a lower limit, the uncertainties of these parameters only come from the mass.
We find 13 B-type stars and a massive binary system, 15 Mon, located in the S Mon region using SIMBAD 1 .Most B-type stars are located in the shell of the bubble rather than at its center (see Figure 19 in Appendix C), and only two B-type stars are located at the center of the bubble.We employ Equation (2) from Arce et al. (2011) to estimate the mass-loss rate ( ṁw ), which can be calculated by: where v w is the wind velocity, and τ w is the wind timescale (i.e., the duration of wind activity).Here, the kinetic timescale (t kinetic ) is employed to estimate the mass-loss rate of the driving source roughly.Given the typical wind velocity of O-or B-type stars is about 1 000-2 000 km s −1 (Chen et al. 2013), the mass-loss rate would be 10 −6 − 10 −5 M ⊙ yr −1 .The value is much larger than the typical mass-loss rate of B-type stars (10 −8 − 10 −11 M ⊙ yr −1 ) but close to that of O-type stars (10 −6 − 10 −7 M ⊙ yr −1 , Chen et al. 2013).Therefore, it is improbable that the B-type stars are the main driving source of the bubble.Although 15 Mon is not precisely at the center of the bubble, it stands out as the closest high-mass star to the bubble's center, apart from the two B-type stars.Therefore, 15 Mon might be the main driving source of the bubble.The dynamic age of the H II region can be used to judge whether the bubble is triggered by 15 Mon.Assuming that the H II region expands in a homogeneous medium, its dynamical age, t dyn is (see the model described by Dyson & Williams 1980): where R HII is the radius of the H II region (Xu et al. 2017), for which we adopt an inner radius of 0.7 ± 0.1 pc.c s is the isothermal sound speed of the ionized gas, and a value of 10 km s −1 is adopted.The radius of the associated Strömgren sphere (R s ) reads: where α B = 2.6 × 10 −13 cm 3 s −1 is the hydrogen recombination coefficient to all levels above the ground level, n 0 is the initial number density of the gas, and Q Ly is the ionizing luminosity.We only consider the number of atoms that can be ionized, so n 0 = M shell 4/3πR 3 bub ≈ (2.6-6.3)×10 4 cm −3 .Since the actual ionized medium within the bubble is likely less massive as compared to the whole shell mass determined from the broad CO ring, the number of atoms that can be ionized might be overestimated.For an O7V star, the ionizing luminosity (Q Ly ) is 10 47.62 s −1 (Martins et al. 2005).Based on these calculations, the obtained dynamical age is (3.5-8.9)×10 5 yr.Even if the secondary member of 15 Mon, an O9.5V star, is also taken into account, the dynamical timescale is only reduced by less than 0.1%.We conclude that the kinetic timescale of the bubble is comparable to the dynamical age of the H II region.Montillaud et al. (2019) estimated the age of 15 Mon using two methods.The first method, based on SED fitting, led to a median age of 1.5 − 2 Myr, with a dispersion between 0.2 and 3 Myr.The second method was based on isochrone fitting in color-magnitude diagrams and led to a median age of 3 Myr, with a dispersion between 1 and 6 Myr.Our estimation on the dynamical timescale of the H II region may be less than the age of 15 Mon.Hence, the bubble is likely to be driven by the expanding H II region of 15 Mon.As the projected position of 15 Mon is not located at the center of the bubble, we presume that 15 Mon may have moved away from its once central position in the bubble.

Energy Cascade in S Mon
Stellar winds inject energy and momentum into the ambient clouds that may help sustain turbulence and disrupt their surroundings (Nakamura & Li 2007;Arce et al. 2011).To assess whether winds can sustain the turbulence in S Mon, we compare the wind energy injection rate into the cloud, Ėw , with the cloud turbulent dissipation rate, L turb,cloud .Based on the 13 CO emission, we can derive the 3D turbulent velocity dispersion (σ 3d ) with an approximation of σ 3d = √ 3σ 1d (Li et al. 2018), where σ 1d is the 1D turbulent velocity dispersion along the line of sight and we consider the mass (M cloud ) of S Mon.The obtained σ 3d is 1.7 ± 0.2 km s −1 and M cloud is 2 200 ± 500 M ⊙ .The cloud turbulent energy, E turb,cloud , is about (6.3 ± 2.5) ×10 46 erg (see the calculation in Appendix B).The cloud turbulent dissipation rate can be calculated as: where t diss is the cloud turbulent dissipation time.We roughly estimate t diss via (McKee & Ostriker 2007): where the cloud diameter is d cloud ∼ (4.0 ± 0.4) pc and σ 1d ∼ (1.0 ± 0.1) km s −1 .Based on Equations ( 5)-( 6), we obtain a cloud turbulent dissipation time of t diss ∼(2.0 ± 0.4) ×10 6 yr, and a cloud turbulent dissipation rate of ∼ (0.5-1.7) ×10 33 erg s −1 .
Since the bubble might be driven by the expanding H II region of 15 Mon, the dynamical age (t dyn ) of the bubble is more suitable than the kinemtaic timescale (t kinetic ) to estimate the mass-loss rate ( ṁw ).Therefore, we reestimated this value, which is approximately (1.9-8.6)×10−6 M ⊙ yr −1 , and the wind energy injection rate can be estimated using the following equation from McKee (1989): The derived wind injection rate is about (4.5-8.1)×10 33erg s −1 , which is larger than the cloud turbulent dissipation rate.This indicates that strong winds serve as a significant energy origin for S Mon, consisting of sustaining the observed turbulence.This result agrees with the findings of other studies (Arce et al. 2011;Li et al. 2015;Feddersen et al. 2018).
In the following, we investigate the turbulent energy in this region.The total wind energy input into the interstellar medium is about E w = 0.5 ṁw (t)v 2 w dt ≃ ∆M v 2 w /2 ≃ (0.7-1.7) ×10 50 erg (Lamers & Cassinelli 1999).Both the wind energy of 15 Mon and the kinetic energy of the bubble (E shell ∼ (1.5 ± 0.4) ×10 47 erg) are higher than the turbulent energy of the cloud (E turb,cloud ∼ (6.3 ± 2.5) ×10 46 erg), thus both of them are help to sustain the turbulence in the S Mon region.Moreover, the wind energy probably plays a crucial role in driving the bubble, since it is about three orders of magnitude higher than the kinetic energy of the bubble.Despite the substantial energy carried away by ionizing photons from the star throughout its lifetime, only 0.01-0.1 percent of this energy is converted into kinetic energy in the shell, with the majority being lost as radiation (Geen et al. 2015).Walch et al. (2012) conducted hydrodynamic simulations to explore the dynamical effects of a single O7 star.They determined that the total radiative energy injected by the ionizing source reaches 7 × 10 51 erg within approximately 1 Myr.A portion of this injected energy is then transformed into kinetic energy, amounting to around 3 × 10 48 erg, which surpasses the kinetic energy we estimated in this work (E shell ∼ (1.5 ± 0.4) ×10 47 erg).Thus, radiative energy contributes to driving the bubble.
Next, we studied the turbulent energy of the outflow of clump E. The total kinetic energy of the outflow in clump E ( E lobe ∼ 4.4 × 10 45 erg) is higher than its turbulent energy (E turb,clump ∼ (5.2 ± 1.1) ×10 44 erg), and it accounts for E lobe /E turb,cloud ∼ 4 − 10% of the turbulent energy of the cloud.This suggests that the outflow in clump E helps to maintain the turbulence in this clump and introduces additional energy to surrounding gas, thereby affecting the turbulence of the local environment (Arce et al. 2010;Li et al. 2018).This result is also consistent with the conclusion made by Li et al. (2020) that the outflow is sufficient to maintain turbulence on a scale of 0.1-0.4pc.Moreover, we also estimate the gravitational binding energy of the S Mon region.The detailed calculation is presented in Appendix B. We find that the total wind energy of 15 Mon is over two orders of magnitude higher than the gravitational binding energy of the cloud, i.e., (2.1 ± 1.0)×10 47 erg, indicating that the wind can disrupt the cloud.Based on the above calculations, we find that half of the gravitational binding energy of the cloud is approximately equal to the kinetic energy of the bubble, and much higher than the total kinetic energy of the outflow in clump E. According to the virial theorem (Nakamura & Li 2014), the bubble has the potential to disrupt the cloud, while the outflow is unable to disrupt the cloud.
We also investigated the gravitational binding energy of clump E (E grav,clump ∼(2.1 ± 0.2)×10 44 erg), which is lower than the total kinetic energy of the outflow, indicating that the outflow contributes to disrupting the clump.The results for the outflow found here are consistent with the conclusions made by Li et al. (2020) in that the outflow activity could potentially disperse material away from the parent cloud on a scale of 0.2 − 1.0 pc.
In conclusion, we find the bubble produced by 15 Mon helps to maintain the turbulence in S Mon and is potentially disrupting the cloud, and while the outflow from clump E only adds a small fraction of the overall turbulence and energy of the system, it still may be able to perturb the local gas.

Bubble from the YSO Perspective
Since some YSOs may be too young to have moved away from their birth sites, they often have the same kinematic properties as their parent clouds (Tobin et al. 2015;Hacar et al. 2016;Großschedl et al. 2021).Thus, we can utilize the motions of YSOs associated with a cloud to represent the motions of the cloud.According to the correlation between YSOs and gas (see Section 3.1 for details), the Class II YSOs, whose parallax accuracies are better than 10% and projected locations are approximately aligned with the cloud were initially selected to investigate the kinematic characteristics of the bubble.
A Gaussian function was used to estimate the average proper motion of the YSO sample and to remove outliers (beyond a 3σ threshold).A total of 38 YSOs were ultimately selected to investigate the kinematic features of the bubble on the plane of the sky. Figure 14 (a) shows the distribution of proper motion of the selected YSOs.The mean proper motion of these sources is ( μl , μb ) = (2.65,-3.17) mas yr −1 , and the 1σ uncertainty range of the Gaussian fitting is (0.14, 0.17) mas yr −1 .
The intrinsic movement of the Class II YSOs on the plane of the sky was calculated by subtracting the mean proper motions.Then the tangential velocities in l and b directions (in units of km s −1 ) are calculated as follows: where κ = 4.74 km yr s −1 , ϖ is the parallax in units of mas, and µ l and µ b are the proper motion along the l and b directions, respectively, in units of mas yr −1 .Figure 14 (b) shows the projected distributions and tangential velocities of the Class II YSOs.The Class II YSOs are mainly located in the molecular gas surrounding the bubble, and their motions are relatively random on the plane of the sky. Figure 15 shows P-V diagrams of the selected YSOs.There is no clear evidence of expansion in the figure .In order to determine whether the YSOs are expanding, we regard the expansion velocity as the radial outward component of the intrinsic velocity from the center of the bubble.The median expansion velocity of the Class II YSOs is found to be 0.0 ± 0.1 km s −1 .Therefore, the YSOs associated with the gas may not exhibit expansion on the plane of the sky.

SUMMARY
We observed four molecular lines ( 12 CO, 13 CO, C 18 O, and HCO + J = 1-0), and collected YSOs in S Mon to investigate the bubble, an outflow identified in one of the clumps, as well as the energy cascade in the S Mon region.This includes the dynamical timescale, possible driving sources, and feedback of the outflow and bubble.The main results are summarized as follows.
1 The Class II YSOs are closely associated with the molecular gas in terms of spatial distribution and kinematic properties.According to the Gaia DR3 parallaxes of the Class II YSOs, we estimate the distance to S Mon to be ∼ 722 ± 9 pc.
2 We discover a molecular bubble in the S Mon region which is consistent with the Herschel H 2 column density maps.The molecular bubble has a radius of ∼ 1.1 ± 0.1 pc, a mass of ∼ 1 400 ± 400 M ⊙ , a momentum of ∼ 4 600 ± 1 300 M ⊙ km s −1 , and kinetic energy of ∼ (1.5 ± 0.4) × 10 47 erg.The dynamical timescale of the H II region of 15 Mon is (3.5-8.9)×10 5 yr, which is likely to be the main driving source of the bubble.
3 We detect a molecular outflow in the shell of the bubble, which indicates that star formation activity is ongoing around the bubble.The outflow has a mass of a few M ⊙ , a momentum of a few tens of M ⊙ km s −1 , kinetic energy of ∼ 10 45 erg, a dynamical timescale of ∼ 10 5 yr, and a mechanical luminosity of a few times 10 −1 L ⊙ .
4 The wind energy helps to sustain the turbulence in the S Mon region and drive the bubble.The bubble can also provide sufficient energy input to sustain the cloud turbulence.The measured outflow in clump E is enough to maintain the turbulence of the clump E. Outflows (like the one detected in clump E) could generally contribute to the turbulence of a cloud, while in S Mon further outflows are not yet confirmed.where T mb,12 is the peak main-beam temperature of 12 CO.
The optical depth can be derived using the following equation (Garden et al. 1991): The relation N H2 /N13 CO ≈ 5 × 10 5 (Simon et al. 2001) is used to estimate the H 2 column density.The mass, M , can be determined by: M = µm H N (H 2 )S, (B4) where µ = 2.72 is the mean molecular weight, m H is the mass of the hydrogen atom (Garden et al. 1991), and S is the projected 2D area.The turbulent energy, E turb , can be given approximately by: where σ 3d is the 3D turbulent velocity dispersion, which can be calculated by: where ∆V FWHM is the 1D FWHM velocity dispersion based on typical 13 CO.
The gravitational binding energy, E grav , can be calculated by: where R is the radius.The excitation temperatures of the eight clumps are between 21.3 and 33.2 K. Since these clumps are regularly located within the bubble, a mean excitation temperature of about 26.9 K is adopted as that of the bubble.The mean column density of the bubble is ∼ 2.5 × 10 20 cm −2 and its mass is ∼ 1 400 M ⊙ .In this work, we consider the uncertainties in gas temperature, shell region, and distance for the shell mass estimation.The total uncertainty caused by these factors is about ±400 M ⊙ .
Figure 18 shows the projected area of the cloud.The red lines denoted its boundary, determined by pixels whose main-beam brightness temperature is larger than 3 × rms in at least three successive channels.The velocity range adopted here is between 7 and 15 km s −1 and the mass of the cloud is derived as 2 200 ± 500 M ⊙ .The uncertainties in the excitation temperature and distance are considered in the estimation of the cloud mass.σ 1d is obtained from Gaussian fitting of the average 13 CO line, which is about 1.0 ± 0.1 km s −1 .The turbulent energy and gravitational binding energy of the cloud are estimated by Equation (B5) and Equation (B7), respectively, which correspond to (6.3 ± 2.5)×10 46 erg and (2.1 ± 1.0)×10 47 erg.
The clump's column density (N H2 ), mass (M clump ), turbulent energy (E turb,clump ), and gravitational binding energy (E grav,clump ) were calculated using the same equations as applied to the cloud.The uncertainties in the distance and clump region were considered in the calculation.
CDF of Gaia distance

Figure 4 .
Figure 4. (a) Projected distribution of Class II YSOs with parallax uncertainties smaller than 10%.The yellow pluses are the Class II YSOs used for the distance estimation.The blue pluses show YSOs with LSR velocities smaller than 7 km s −1 .The magenta pluses show YSOs which do not coincide with the gas distribution.The red triangle represents the projected position of 15 Mon.(b) CDF of Gaia distances.The red vertical line indicates the median value of 722 pc, with a standard error of 9 pc denoted by the gray shadow.The upper and lower percentiles within 1σ around the median are labeled by the dashed gray lines.

Figure 5 .Figure 6 .
Figure 5. Channel maps of 12 CO emission.The red triangle represents the star 15 Mon.The blue dashed circle shows the bubble obtained by fitting the gas, which has a radius and thickness of ∼ 5.4 and ∼ 4.2 arcmin, respectively.The magenta circle indicates the beam size of PMODLH.The velocity range for each map is labeled in the upper right of the map.
Figure 7. (a) Integrated map of 12 CO emission(blue contours) superposed on the Herschel H2 column density map (gray background, Nony et al. 2021), and the green contour lines are the emission of C 18 O.The integrated velocity range is 7.5 − 12.5 km s −1 for both 12 CO and C 18 O, and the contour levels are shown from 30% to 90% with steps of 10% of the peak integrated intensity of the corresponding emission.The orange dashed circle shows the bubble.Clumps identified based on the C 18 O integrated map are labeled A to H. (b) P-V diagram of 12 CO along the arrow in panel (a), where the contour levels are shown from 10% to 90% with steps of 10% of the peak value.The red dashed lines show the expansion velocity range from the visual inspection, i.e., 8.2 − 11.3 km s −1 .(c) Azimuthally averaged radial intensity profile of the 12 CO bubble (black curve), where a Gaussian fit (red curve) to the intensity profile was performed to estimate the radius (peak) and thickness (Full Width at Half Maximum, FWHM).(d) Integrated map of HCO + (magenta contours).The other labels are the same as in panel (a).(e) P-V diagram of HCO + along the arrow in panel (d), and the expansion velocity range is 8.4 − 11.8 km s −1 .The other features are the same as in panel (b).(f) Azimuthally averaged radial intensity profile of the HCO + bubble, where the other labels and symbols are the same as in panel (c).
Figure 8.(a) Herschel H2 column density map, showing a larger area to the southwest of the investigated region (Nony et al. 2021).The red circle depicts the bubble we identified, and the black box shows the region housing the small arc-like structure.(b) Integrated map of 13 CO emission (green contours) in the boxed region in panel (a).The integrated velocity range is 9.5 − 11.0 km s −1 , and the contour levels are shown from 30% to 90% with steps of 10% of the peak integrated intensity of 13 CO emission.(c) P-V diagram of 13 CO along the arrow in panel (b), where the contour levels are shown from 10% to 90% with steps of 10% of the peak integrated intensity of the corresponding emission.

Figure 9 .Figure 10 .
Figure 9. Spectra of 12 CO, 13 CO, C 18 O, and HCO + of the C 18 O peak position of clumps A to H.The red lines represent the sum of the multiple Gaussian fitting components.

Figure 11 .
Figure 11.Distribution of the velocity dispersion (σ).The green circles indicate the bubble region.
Source name.(2) Excitation temperature of the clump, which is obtained by the 12 CO spectral line from the C 18 O peak emission of each clump.(3) Diameter of the clump.(4) Integrated intensity of 13 CO emission.(5) H 2 density of the clump.(6) Mass of the clump.(7) Turbulent energy of the clump.(8) Gravitational energy of the clump.
Figure12.Outflow in clump E identified by 13 CO emission.(a) Integrated C 18 O map with 13 CO blue and red lobe contours.The contour levels are shown from 40% to 90% with steps of 10% of the peak intensity of each lobe.The black diamonds and cyan pluses represent Class I and Class II YSOs, respectively.The magenta, blue, and red stars are the emission peaks of C 18 O, and the blue lobe and red lobe, respectively.(b) P-V diagram along the white arrow in panel (a).The black contour levels are shown from 10% to 90% with steps of 10% of the peak value.The blue and red dashed lines are the velocity ranges of the blue and red lobes, respectively.(c) The blue spectrum presents 13 CO at the peak position of the blue emission of the 13 CO outflow.The black spectrum is C 18 O at the emission peak of C 18 O.The blue shading of the spectrum indicates the blue line wing velocity of 13 CO.The feature at 8.6 km s −1 is likely unrelated to the outflow, so emission with a velocity below 8.6 km s −1 is not included in the blue wing of the outflow profile.(d) The red spectrum is 13 CO at the peak position of the red emission of the 13 CO outflow.The black spectrum is C 18 O at the emission peak of C 18 O.The red shading of the spectrum indicates the red line wing velocity of 13 CO.

Figure 13 .
Figure13.Outflow in clump E identified by HCO + emission.The description of each panel is the same as that in Figure12, except the contour levels in panel (a) are shown from 50% to 90% with steps of 10% of the peak intensity of each lobe.

Figure 14 .Figure 15 .
Figure 14.(a) Distribution of the proper motions of the YSOs, where the average proper motion of these sources is (2.65,-3.17)mas yr −1 , and the 1σ uncertainty range of the Gaussian fitting is (0.14, 0.17) mas yr −1 .(b) Projected distribution and tangential velocity of the YSOs.The magenta arrows depict the velocity for each YSO with respect to the mean velocity of the YSO sample.The red asterisk (0, 0) is the center of the bubble.

Figure 16 .
Figure 16.Color-color diagrams showing the classified YSOs.The green diamonds, red triangles, and magenta crosses present Class I, Class II, Class III YSOs, respectively.

Figure 17 .
Figure 17.Comparison of the VLSR values in the APOGEE and Tobin et al. (2015) samples.The black line represents where VLSR is the same in both catalogs.The red dashed line indicates the mean value of the difference between the APOGEE and Tobin et al. (2015) VLSR values, and the gray filled region depicts the corresponding 1σ region.

Figure 18 .Figure 19 .
Figure 18.Intensity map of 13 CO for the velocity range 7-15 km s −1 .The red lines depict the boundary of the cloud.

Table 1 .
Observational Parameters of the Molecular Lines

Table 2 .
Observed Parameters of Each Clump

Table 3 .
Physical Parameters of the Clumps

Table 4 .
Physical Properties of the Outflow of Clump E

Table 5 .
Description of the YSO SamplesNote-This table is available in its entirety in machine-readable form.