Unveiling the Cosmic Cradle: clustering and massive star formation in the enigmatic Galactic bubble N59

In this paper, we have conducted an investigation focused on a segment of the $Spitzer$ mid-infrared bubble N59, specifically referred to as R1 within our study. Situated in the inner Galactic plane, this region stands out for its hosting of five 6.7 GHz methanol masers, as well as numerous compact H II regions, massive clumps, filaments, and prominent bright rims. As 6.7 GHz masers are closely linked to the initial phases of high-mass star formation, exploring regions that exhibit a high abundance of these maser detections provides an opportunity to investigate relatively young massive star-forming sites. To characterize the R1 region comprehensively, we utilize multi-wavelength (archival) data from optical to radio wavelengths, together with $^{13}$CO and C$^{18}$O data. Utilizing the $Gaia$ DR3 data, we estimate the distance towards the bubble to be $4.66 \pm 0.70$ kpc. By combining near-infrared (NIR) and mid-infrared (MIR) data, we identify 12 Class I and 8 Class II sources within R1. Furthermore, spectral energy distribution (SED) analysis of selected sources reveals the presence of four embedded high-mass sources with masses ranging from 8.70-14.20 M$_\odot$. We also identified several O and B-type stars from radio continuum analysis. Our molecular study uncovers two distinct molecular clouds in the region, which, although spatially close, occupy different regions in velocity space. We also find indications of a potential hub-filament system fostering star formation within the confines of R1. Finally, we propose that the feedback from the H II regions has led to the formation of prominent Bright Rimmed Clouds (BRC) within our region of interest.


INTRODUCTION
Massive stars (> 8 M ⊙ ), characterized by their substantial mass and intense luminosity, hold a pivotal role in various astrophysical phenomena.Their significance spans several domains, from contributing to the synthesis of heavy elements and enriching the interstellar medium (ISM) (Zinnecker & Yorke 2007;Motte et al. 2018), to influencing the dynamics of their host galaxies through their radiation and stellar winds, thereby shaping gas motions and fostering the birth of new stars.Despite their important role, the study of massive star formation presents challenges due to their rarity, rapid evolution, and their concealment within dusty cocoons during their early stages.
One promising avenue for probing massive star-forming regions involves the detection of 6.7 GHz methanol masers, as these are closely associated with the initial phases of massive star formation (Breen et al. 2013;Paulson & Pandian 2020).In this study, we focus our attention on the specific star-forming region known as MIR bubble N59.This region stands out for hosting approximately eight 6.7 GHz methanol masers, alongside numerous compact H ii regions, massive clumps, filaments, and prominent bright rims.Analyzing the characteristics of this region, therefore, allows us to glean valuable ★ E-mail: sonu.paulson@tifr.res.ininsights into the interplay between the processes involved in the early phase of high-mass star formation.
The Galactic mid-infrared (MIR) bubble N59, initially documented by Churchwell et al. (2006), is situated at coordinates (l,b)= (33.071 • ,-0.075 • ).Deharveng et al. (2010) offer a concise description of bubble N59, highlighting the presence of multiple bright condensations and radio sources identified from the catalogues of Becker et al. (1994) and Kurtz & Hofner (2005).They also report the identification of two 6.7 GHz methanol masers within the confines of the bubble, referencing compilations from Pestalozzi et al. (2005) and Xu et al. (2008).Notably, these masers are situated in the vicinity of bright dust condensations that are interacting with H ii regions, implying that N59 holds significant potential as a candidate for investigating the mechanisms driving triggered massive star formation.Hattori et al. (2016) have also targeted bubble N59 as part of their sample of 180 IR bubbles for a statistical study of mid-and far-infrared properties of Spitzer Galactic bubbles.They find N59 to be a closed bubble that hosts polycyclic aromatic hydrocarbons (PAH) with luminosities exceeding 10 5 L ⊙ .Subsequently, Hanaoka et al. (2019) conducted a follow-up study involving 247 IR bubbles, including N59, incorporating data from the Herschel survey.The investigation revealed that the fractional PAH luminosity (i.e.L PAH /L TIR ) of N59 was approximately 0.17.Based on the L TIR value (Martins et al. 2005), presence of O-type stars was proposed as one of the reasons for the low L PAH /L TIR , as the intense UV luminosities from such stars would enhance L TIR and expedite the photodissociation of PAH molecules.Thus, N59 turns out to be harbouring several bright stars, including massive stars.While bubble N59 has been featured in several statistical investigations concentrating on dust properties, it has not yet been subjected to a comprehensive multiwavelength analysis.Given the intriguing findings in previous research, such as distinctive dust distribution, the presence of 6.7 GHz methanol masers near prominent dust condensations, and indications of potential triggered star formation, this region merits an in-depth examination.
Fig. 1(a) displays a two-color composite image of N59, showcasing the dust emission at 8.0 µm (green) and 24.0 µm (red) from Spitzer-MIPSGAL.The bubble stands out prominently within this composite map, spanning an area of ∼30 ′ ×30 ′ .It is worth highlighting that a significant concentration of radio sources, ATLASGAL clumps (Schuller et al. 2009), and methanol masers is primarily situated within the cyan square measuring 10' X 10' (hereafter referred to as 'R1').
An enlarged view of R1 is provided in Fig. 1(b), which is a threecolor-composite map created using the Herschel data, where the red, green and blue components correspond to 250 µm, 160 µm and 70 µm, respectively, for region R1.The map also features the position of Galactic filament (marked using the green diamond) reported by Li et al. (2016), who discerned it by employing the DisPerSE algorithm (Sousbie 2011) on the ATLASGAL survey data.The magenta crosses represent 6.7 GHz Class II methanol maser detections.Additionally, the map reveals the presence of multiple bright condensations that coincide with 870 µm ATLASGAL clumps.Numerous radio sources (the blue plus symbols; see Section 2.6) are concentrated within the bubble, primarily around the brightest condensations.The 6.7 GHz methanol masers in R1 were identified using the online tool, MaserDB1 , provided by Ladeyschikov et al. (2019).As the majority of ATLASGAL sources, H ii regions, and methanol maser detections are concentrated in R1, our study focuses extensively on this specific area.

ARCHIVAL DATA
Table 1 presents a summary of the archival data used in this paper.The subsequent portion of this Section offers a concise description of each data sets.

Gaia DR3
The third Gaia data release (Gaia DR3) contains the positions, parallaxes, proper motions, and apparent brightness in the G-band (0.33 − 1.05 µm) for about 1.8 billion stars, based on 34 months of mission observations (Gaia Collaboration, Vallenari et al. 2023).Incorporated within Gaia DR3 are measurements of the apparent brightness in the G  (0.33 − 0.68 µm) and G RP (0.63 − 1.05 µm) bands, providing a wealth of broad-band color information.We used the Gaia data to retrieve the proper motion of stars as well as to determine the distance of the complex.The data was downloaded from the Gaia archive2 .

Near-Infrared Data from UKIDSS
The near-infrared (NIR) J (1.25 µm), H (1.65 µm) and K (2.16 µm) data are retrieved from the UKIRT Infrared Deep Sky Survey (UKIDSS) DR10PLUS Galactic Plane Survey (GPS; Lawrence et al. 2007).The UKIDSS data is obtained using the Wide Field Camera (WFCAM) on the United Kingdom Infrared Telescope.The UKIDSS data is publicly available in WFCAM science archive.Following the selection procedure described in Lucas et al. (2008) and Dewangan et al. (2015), we downloaded the NIR data using the Structured Query Language (see Appendix A).In order to ensure reliable photometry, we discarded the sources that have photometric uncertainty > 0.1 magnitude in each band.

Mid-Infrared Data from Spitzer-IRAC
We obtained the archival mid-infrared (MIR) data using the Spitzer Space Telescope observations under the Galactic Legacy Infrared Mid-Plane Survey Extraordinaire (GLIMPSE) program (Benjamin et al. 2003;Churchwell et al. 2009).GLIMPSE is a legacy science program that mapped the inner Galactic plane using the Infrared Array Camera (IRAC) across four distinct bands: 3.6 (I1), 4.5 (I2), 5.8 (I3), and 8.0 (I4) microns.The angular resolution in the four bands is less than 2 ′′ (Fazio et al. 2004).The GLIMPSE point sources were retrieved from the highly reliable GLIMPSE 1 Spring '07 archive3 .

Far Infrared and Sub-millimeter Data
The far-Infrared (FIR) continuum images are retrieved from the Hershel Space Observatory data archives.We made use of the processed level2_5 images at 70, 160 µm (PACS) as well as 250, 350 and 500 µm (SPIRE), observed as a part of the Herschel infrared Galactic plane survey (Hi-GAL; Molinari et al. 2010a).The maps at different wave bands have different data units, resolution and plate scales.The 70 µm, 160 µm, 250 µm, 350 µm, and 500 µm images have pixel scales of 3.2 ′′ , 6.4 ′′ , 6 ′′ , 10 ′′ , and 14 ′′ , respectively, with resolutions varying from ∼ 5.5 ′′ − 36 ′′ .While the PACS images obtained had the surface brightness unit of Jy pixel −1 , the SPIRE images were in the units of MJy sr −1 .The sub-millimeter continuum map at 870 µm (beam size ∼ 19.2 ′′ ) was retrieved from the ATLASGAL archival survey (Schuller et al. 2009).The FIR and sub-millimeter data are employed to analyze the prevailing physical conditions within the region.

Radio Continuum Data
For obtaining the radio continuum image of the region, we have used the radio interferometer surveys of the Galactic Plane, namely CORNISH (5 GHz; Hoare et al. 2012) Condon et al. 1998).The angular resolution for CORNISH, MAGPIS and NVSS surveys are 1.5 ′′ , 6.0 ′′ and 45 ′′ respectively.While CORNISH and MAGPIS concentrate on compact sources (at high resolution and low resolution respectively), NVSS excels in capturing extended emissions (Thwala et al. 2019).We also obtained archival VLA data at 8.49 GHz for a small part of our R1 region from the NRAO VLA Archive Survey (NVAS) Images Page 5 (Obs Date : 1997 Nov 22, Version : 2007 July 10).

Molecular CO Data
The J= 3−2 line of the 13 CO data is obtained from the CO Heterodyne Inner Milky Way Plane Survey (CHIMPS; Rigby et al. 2016).This high-resolution spectral survey has been carried out using the Heterodyne Array Receiver Program on the 15 m James Clerk Maxwell 5 https://www.vla.nrao.edu/astro/nvas/Telescope (JCMT) in Hawaii.CHIMPS has an angular resolution of 15 ′′ in 0.5 km s −1 velocity channels.

Notes on Distances
An accurate determination of the distance to the star-forming region R1 is essential for estimating parameters like the size and mass of the associated dust condensations.Deharveng et al. (2010) report a kinematic distance of 5.6 kpc for N59, derived primarily from the velocity measurements of the ionized gas using radio recombination lines.This calculation assumes circular rotation around the Galactic center.However, the bubble lies within the "Solar Circle" (i.e.≲ 8 kpc) and therefore carries a distance ambiguity due to two potential solutions (e.g.Kolpak et al. 2003 and references therein).This can be avoided by using the distances derived from parallax.
Thanks to the new Gaia data release (DR3; Vallenari et al. 2023), we now have precise parallax measurements for over 1.8 billion stars, which can be leveraged to derive robust distance estimates.The inversion of parallax to distance is, however, non-trivial, mainly due to the positivity constraint of distance and the non-linear transformation (Luri et al. 2018).In fact, reliable distances cannot be obtained for the vast majority of stars through this inversion (see the discussion by Bailer-Jones et al. 2018;Bailer-Jones et al. 2021).As such, the most appropriate way to implement this conversion is through a probabilistic approach.In this regard, Bailer-Jones et al. ( 2021) present a promising method that uses parallax to provide the full likelihood distribution of the distance.These distances have even been shown to carry smaller uncertainties than kinematic distances, especially for the stars situated at distances less than 5 kpc (Karska et al. 2022).Considering these advantages, we chose to utilize the stellar distances provided by Bailer-Jones et al. (2021).
To this end, we considered all stars belonging to region R1 (i.e, we are not applying the clustering algorithm on any pre identified cluster), with parallax uncertainties below 10 per cent, and computed their membership probabilities using the Unsupervised Photometric Membership Assignment in Stellar Clusters algorithm (UPMASK). 6he algorithm has the advantage of being non-parametric and unsupervised, that is, no a priori selection of field stars is required to serve as a comparison model -which is generally the case in many of the previous cluster membership calculators.Based on the values thus obtained, we selected the stars with membership probabilities ≥ 80 per cent, and used them to infer the mean distance to R1: 4.66 ± 0.70 kpc.More details on calculating membership probabilities using PyUPMASK are provided in Appendix B

Young Stellar Population
The evolutionary state of a young stellar object (YSO) is often signposted by its infrared color excess -an indicator of ionized stellar wind or circumstellar dust (Deng et al. 2022).The infrared colorcolor (CC) and color-magnitude (CM) diagrams display distinct and specific regions occupied by various object types, indicating that each type is well-defined.Additionally, the position of a young stellar object in CC/CM space can provide some valuable information about its evolutionary stage.We utilize the GLIMPSE 1 and UKIDSS data to investigate the infrared excess sources present in R1.However, the selection of YSOs based on excess emission in the infrared is easily contaminated by external polycyclic aromatic hydrocarbon (PAH)/star-forming galaxies, active galactic nuclei (AGNs), or shockexcited extended sources, as they often mimic YSO colors (Gutermuth et al. 2009a).Hence, we eliminated these contaminants prior to classifying the YSOs, using the methods described in Gutermuth et al. (2009a).The following steps were adopted to identify and classify YSOs: (i) IRAC I1-I2-I3-I4 catalog: Initially, the YSOs were identified using their MIR magnitudes.For this purpose, we obtained the Spitzer sources from the GLIMPSE 1 catalog within a search radius of 9 ′ centered at (l,b) = (33.068• , -0.044 • ).We ensured that the errors in each individual band are ≤ 0.15 mag.This generated a comprehensive catalog of 2242 sources.We employed the Phase 1 method described in Gutermuth et al. (2009a, Appendix A.1)  (ii) H-K-I1-I2 catalog: There have been cases where there were non-detections at 5.8 and 8.0 µm, but good quality detections in the NIR bands.We used a combination of H, K (from UKIDSS survey data), and Spitzer-GLIMPSE I1, I2 data, for classifying the YSOs in such scenarios.The H, K, and I1, I2 sources were cross-matched with a 0.6 ′′ matching radius, resulting in a catalog containing 4079 sources.Using the classification scheme described in Gutermuth et al.The dereddened colors were obtained using the color excess ratio (Flaherty et al. 2007).We obtained an additional 6 Class 0/I sources using this method (see Fig. 2(c)) with one of them being associated with region R1.
In summary, within region R1, we detected 11 Class I and 8 Class II sources exclusively using the IRAC catalog, and an additional Class I source by combining data from both NIR and IRAC catalogs.Yasui et al. (2008).The dot-dashed black lines are the reddening vectors, drawn using the reddening laws from Rieke & Lebofsky (1985).Three regions-"F," "T," and "P"-mark the location of different classes of YSOs, with sources in each region colored in orange, yellow and brown, respectively.The horizontal gray line has been drawn at J -H = 0.89.All points are in the Mauna Kea Observatory (MKO) system.This brings the total count of identified YSOs in R1 to 12 Class I and 8 Class II sources.The locations of these identified YSOs are shown in Fig. 3.

Spectral Energy Distribution of YSOs
We conducted spectral energy distribution (SED) modeling for a subset of 12 Class I and 8 Class II YSOs, to estimate their physical parameters.To accomplish this, we employed the grid of YSO models developed by Robitaille et al. (2006) and utilized the online SED fitting tool created by Robitaille et al. (2007).The underlying model includes a pre-main sequence (PMS) star surrounded by a flared accretion disk, which is enveloped by a rotationally flattened envelope featuring cavities carved out by a bipolar outflow.The SED models were computed using the radiation transfer code of Whitney et al. (2003,2004) in a 14-dimensional parameter space, encompassing properties of the central source, infalling envelope, and disk.In the fitting process, a total of 200,000 SED models were generated, and each fitting was characterized by a  2 parameter, comparing the available SED models with the data.The range for distance and interstellar visual extinction (A V ) were considered as free parameters, with a distance range of 4.5 to 5.2 kpc for our sources.
To quantify the extinction, we employ a method that involves measuring the color excess in the infrared wavelengths, as described by Lada et al. (1994).In our study, we utilized the NIR J -H/H -K CC diagram of the region R1 to identify sources displaying reddening.We specifically focused on sources within the "F" region (refer to Fig. 4) that were not initially classified as YSOs.We only considered those sources situated above the J -H = 0.89 threshold as a measure to reduce potential contamination.To estimate the average visual extinction in this region, we adapted a method similar to the Near-IR-Color-Excess technique proposed by Lada et al. (1994), Ojha et al. (2004) and Kainulainen et al. (2007).Initially, we calculated the color excess (H -K) for each selected sources by dereddening them along the reddening vectors until they aligned with the approximate straight-line locus representing dwarf stars.Subsequently, A V was computed for each source using the reddening laws established by Rieke & Lebofsky (1985).The A V values lie in the range of 2.5-9 mag.Furthermore, we plot the K/H-K colour-magnitude diagram (CM-D) for the entire bubble (refer Fig. 5).The almost vertical solid lines represent the zero-age main sequence (ZAMS) loci at a distance of 4.66 kpc reddened by A V = 0, 1, 10, 20 and 30 mag.Slanting lines indicate the reddening vectors for the marked spectral types.The 59 YSOs obtained (see sect.3.2) are depicted as tiny star symbols, with those belonging to region R1 highlighted in blue.Notably, a majority of our YSOs are distributed between A V = 1 and 10 mag.Consequently, we delineate the interstellar visual extinction (A  ) range as 1-10 mag for SED fitting, accounting for uncertainties in our analysis.
Given the extensive number of SED models encompassing a broad parameter space, the fitting process for each source necessitates an increased number of data points and sufficient coverage across the entire wavelength range.Therefore, we selected sources with photometry available in at least all four IRAC bands as a criterion for SED modeling.It is important to note that photometry from Herschel images was not utilized due to the absence of YSO counterparts at those wavelengths, primarily due to poor resolution.
In addition to providing the best-fitting model, the SED fitting tool generates a set of well-fitted models ranked based on their  2 values, which serve as a measure of their relative goodness-of-fit.Following a similar approach as Robitaille et al. (2007), we considered only those models that met the following criterion for further analysis: min < 3 (per data point). (1) As described in Robitaille et al. (2007), it is worth noting that the criterion used for model selection, based on visual examination of SED plots, lacks rigid mathematical justification.Implementing a more stringent criterion could potentially lead to over-interpretation of the results.To determine the physical parameters, a weighted mean value and standard deviation were calculated for each parameter using the subset of models that satisfied equation ( 1).The inverse of the corresponding  2 value was utilized as the weight for each model, following a similar approach as Grave & Kumar (2009).

Mass and Age spectrum
While we conducted SED fitting to all the recognized YSOs within the region R1, a subset of only 14 YSOs satisfied the

Clustering analysis
To investigate the star formation history of R1, it is imperative to identify distinct clusters or groups of YSOs.To accomplish this and assess the spatial distribution of YSOs, we employ the Minimum Spanning Tree (MST) method, which is predicated on the idea that every point data set can be described by a unique, filamentary network (Prim 1957;Barrow et al. 1985).It involves connecting the sources with branches such that their cumulative length is minimised, and the resultant tree is devoid of loops.Then, a cluster is defined as any subset of this network with (at least ) members whose associated branches are shorter than a specified cutoff length ( c ).This enables cluster identification based solely on the spatial distribution of sources, without relying on kinematic information.It has to be noted that MST is particularly effective with a larger number of data points.Therefore, we initially applied the method to the entire bubble N59, for a sample of 59 YSOs (the YSOs obtained in Section 3.2 through (i) and (ii)) and subsequently narrowed our focus to examine the clustering patterns specifically within region R1.While there is not a strict minimum requirement for data points in constructing an MST, its significance lies in how well the data is interconnected.Notably, MST has been successfully applied to YSO samples, even with less than 50 data points (e.g., Gutermuth et al. 2009b;Sharma et al. 2016;Verma et al. 2023).To gauge the reliability of the MST obtained, we examined the histogram illustrating the relationship between MST branch lengths and MST branch numbers for the YSOs (refer to Fig. 7).The histogram reveals a peak at small spacings and a relatively long tail toward large spacings.Such peaked distance distributions often suggest significant subregions above a relatively uniform, elevated surface density.By implementing an MST length threshold (cutoff length,  c ), we can isolate sources closer than this threshold.
Though this method has been successfully applied to numerous nearby star-forming regions (Schmeja & Klessen 2006;Allison et al. 2009;Parker et al. 2012;Saral et al. 2015;Beuret et al. 2017), we note that there is no universally prescribed criterion for determining the appropriate  and  c .Nevertheless, it is common practice to inform  c based on clustering of branch lengths (e.g.Koenig et al. 2008;Gutermuth et al. 2009b;Saral et al. 2015;Pandey et al. 2020), and we adopt the same for our analysis.We remind the reader that the cut-off length is not an extremely rigid value and may vary slightly as the data statistics increase.However, the key takeaway from this analysis remains the identification of a noticeable level of clustering in our specified region R1.
To this end, we first generate the cumulative distribution of the branch length, .The left panel of Fig. 8 shows that the distribution comprises two quasi-linear zones separated by an intermediate zone.
The two quasi-linear zones span 0 ′′ <  < 90 ′′ and 134 ′′ <  < 175 ′′ .The quasi-linear zone at smaller lengths and the intermediate zone essentially represent two levels of clustering, one at a smaller scale than the other.We derive linear fits to these zones (solid brown lines) and take the length corresponding to their intersection as the cut-off -i.e. c = 114 ′′ , which corresponds to 2.59 pc for a distance of 4.66 kpc.
Next, we use this cut-off to identify the small-scale clusters in the tree7 shown in the right panel of Fig. 8, where the YSOs (in blue) are overlaid on the NIR surface density map of the region.The solid and dotted lines are the branches, where the red lines have lengths  <  c .We follow Saral et al. (2015) and use  = 7, considering that this threshold is appropriate for avoiding false identifications stemming from random associations.Thus, the small-scale clusters are those overdensities in the red network that satisfy this threshold.These are demarcated using the dashed-cyan convex-hulls8 (based on the Qhull algorithm; Barber et al. 1996).We term these as the 'core' regions.
Recall that the cumulative  distribution suggests two levels of clustering.This is because the core regions are, in fact, substructures embedded in even larger regions that manifest as the intermediate zone in the left panel of Fig. 8.These regions are traced by the brown branches in the right panel of Fig. 8, corresponding to  < 134 ′′the lower bound of the quasi-linear zone at the higher lengths (left panel; Fig. 8).We define these 'active' regions as overdensities in the brown network (again assuming  = 7), depicted as dashed-yellow convex hulls.
As illustrated in the right panel of Fig. 8, we identify two core regions, labeled as C1 and C2, within R1.Core C1 is found to coin- cide precisely with a prominent dust condensation within the region, highlighting its significance in the context of star formation.Dust, a key player in the star formation process, fulfills a pivotal role by safeguarding molecular clouds from the harmful effects of ultraviolet (UV) radiation while facilitating the essential cooling mechanisms (Williams 2005).Remarkably, three out of the four massive YSOs identified through SED analysis are found within C1 and C2, suggesting that these regions are highly conducive to vigorous star formation.This is further buttressed by the presence of massive AT-LASGAL clumps within these "cores".Furthermore, core C1 stands out with the presence of a 6.7 GHz methanol maser, further reinforcing the notion of active star formation underway within these cores.A more detailed discussion regarding the significance of ATLASGAL clumps and 6.7 GHz masers can be found in Section 4.1.

H II Regions
Fig. 9 shows the compact H II regions in the region from NVSS (red contours), CORNISH, and MAGPIS.Sources from CORNISH and/or MAGPIS have been marked in green crosses and labelled r1-r9.Table 3 shows the details of their detection.In Figs 9(a) and 9(b) -which show the 24m and 70m emission from Spitzer and Herschel, respectively -the contours for NVSS 1.4 GHz emission (red contours) indicate three distinct H II regions, enclosed in dashed boxes labelled B1, B2, and B3.A zoomed-in view of the regions B2 and B3 is shown in Figs 9(c), 9(d), and 9(e) for a better delineation of the features.For the region B1, VLA archival image was available at 8.49 GHz, and is shown in Fig. 10.Emission at 24m and 70m is mainly due to heated dust, such as that by the star formation (Calzetti 2013;Helfand et al. 2006b), and is thus helpful in separating thermal (like H II regions) from nonthermal (like extragalactic origin) sources.We now discuss the H II regions B1-B3 in detail.
For region B1, it can be seen that while the NVSS contours encompassing the "r3-r4" sources have associated thermal dust emission at 24m and 70m, there is no such emission at the location of contours associated with "r5-r6" sources.Combined with the fact that the sources r5 and r6 have been marked as Radio Galaxies in CORNISH catalogue (see Table 3), it can be safely assumed that they are not associated with the N59 region.While the NVSS contours (due to low resolution) subsume both r3 and r4 sources, the background images clearly indicate them to be two separate sources.High-resolution 8.49 GHz image from VLA (Fig. 10) also shows r3 and r4 to be distinct H II regions.Furthermore, in higher resolution 8.49 GHz image, the source r4 appears to be resolved into two distinct compact H II regions (beamsize here is ∼ 9 ′′ ), labelled r4-N and r4-S in Fig. 10.The integrated flux density for r3, r4-N, and r4-S was calculated using the AIPS software (task "JMFIT"), and values returned were used to calculate their spectral types (see Section 3.4.2).
Moving further north, we find weak 1.4 GHz radio continuum emission in the region marked by box B2.It was found to be relatively faint and there are no CORNISH/MAGPIS sources associated with this region.The contours seem to be bounded by 24/70m emission towards their western side (see the zoomed-in view in Fig. 9(c)).This could indicate some embedded source in the west, the ionizing radiation from which is expanding as a possible blister/champagne flow towards the east.
In Figs 9(a) and 9(b), contours in B3 region are associated with the brightest thermal dust emission.Three MAGPIS 6 cm sources were found to be present near the peak of the 1.4 GHz emission, named r7, r8, and r9 here.The zoomed-in view in Figs 9(d) and 9(e) does show a bright source at the location of the peak, but it is not coincident with any of the MAGPIS sources.Coupled with the fact that the   r7-r9 have no detections in MAGPIS 20 cm catalog and CORNISH catalog, it is probable that they are spurious peaks.
Finally, we note the existence of two more CORNISH/MAGPIS sources, named r1 and r2 in Fig. 9.Both these sources are only present in MAGPIS 6 cm catalog with no counterparts in either MAGPIS 20 cm or CORNISH catalog.Both the sources have a high probability of being a sidelobe according to the catalog.In the light of these facts, the source r1 might merit a conclusion of being spurious, but for r2, we find that it has a corresponding source at both 24 and 70m.

Spectral Types
In this section, we calculate the Lyman continuum luminosity for the sources likely to be H II regions.Apropos the discussion in Section 3.4.1,this includes r2, r3, r4-N, r4-S, and NVSS sources in regions B2 and B3.The following expression from Moran (1983) was used for the same : where N  is the Lyman continuum luminosity in photons s −1 , S is the flux density in Jy, D is the distance in kpc,  is the frequency in GHz, and T  is the electron temperature.While for r2, we used the MAGPIS 6 cm (integrated) flux for the above calculation; for r3, r4-N, and r4-S, we use the 8.49 GHz integrated flux obtained using the AIPS software ("JMFIT" task).For the NVSS sources in B2 and B3 boxes, we use the integrated flux from the NVSS catalog.The sources r2, r3, r4-N, r4-S, B2 NVSS source, and B3 NVSS source were found to be of spectral types B1, O9.5, B0.5, B0.5, B0,5, and B0, on a comparison of N c with the tabulated values of Panagia (1973, assuming ZAMS).The results are given in Table 3. Hence there appears to be a rich cluster of O and B type stars in the region.

Molecular Cloud Morphology
We use the J=3-2 transition of 13 CO and C 18 O isotopologues to examine the molecular structure in the region.This transition has a critical density > 10 4 cm −3 (Rigby et al. 2016), with C 18 O(3-2) tracing the densest structures in the region.We first examine the integrated intensity map (moment-0) of the region in 13 CO(3-2) transition in Fig. 11.To avoid confusing artefacts with real features, this image is made using a masked cube of only those regions which show detection above 3 ( being the noise level of the cube) limit in the position-position-velocity (ppv) space.This was achieved by  identifying clumps in the ppv space above the requisite thresholdusing the clumpfind algorithm -and then masking the cube to only show those regions.The masked cube was collapsed in the velocity range ∼ [65, 110] km s −1 (as the molecular emission was found to be in this range only for our region of interest) to make the moment-0 map shown in Fig. 11.We also obtained the linewidth image from the masked cube in this velocity range (discussed in Section 3.5.1).
The image has been overlaid with contours showing local peaks (p1-p13), with spectra extracted therein in 13 CO(3-2) (blue) and C 18 O(3-2) (red) transitions.ATLASGAL sources have also been marked in plus symbols, and seem to be associated with ∼p1/p2, p3, ∼p5, p7, and p11.We examine the 13 CO(3-2) spectra now.The first thing to notice in the extracted spectra is that many of them have complex shapes unlike what a simple gaussian would indicate.For example, while p1, p2, p5, and p8 have broad peaks; p3 and p4 have wide wings towards the blue-shifted velocities; and p11 shows a sharp self-absorbtion feature near the peak.The velocity peaks at p11, and p12 positions show a marked contrast to other regions.While the main velocity peak for p1-p10 regions lies at ∼95-110 km s −1 , the same for p11 and p12 lies at ∼70-80 km s −1 .Among the regions which have velocity peak in ∼95-110 km s −1 range, the peaks p4 and p5 have slightly red-shifted velocities as compared to p1, p2, and p3; with similar trend for p6, p7, and p8 as compared to p9 and p10 (i.e.slightly red-shifted velocity).
It can be noticed that while p1, p2, p3, p9, and p10 have peaks in the range ∼98-102 km s −1 ; p4, p5, p6, and p7 have the peak in ∼100-104 km s −1 range; and p8 has the same in ∼104-108 km s −1 range.At p13, both the velocity components can be seen, albeit the ∼100 km s −1 the component is much stronger.This suggests two different molecular clouds in the region, which spatially might appear as part of the same complex, but are differentiated in velocity space.
To examine the morphology of the two velocity components, we constructed channel maps (using the native cube) for the velocity ranges [65,85] km s −1 and [95,110] km s −1 (Fig. 12).The channel maps clearly reveal two molecular cloud regions.It is worthwhile to note that there is negligible overlap between the spatial locations of  the two molecular cloud structures.In fact, in Fig. 12(b), there is a noticeable "hole" at the location of the blue-shifted velocity cloud.The 13 CO(3-2) integrated intensity peaks p11 and p12 are associated with the blue-shifted molecular cloud.
For comparison, the 13 CO(3-2) linewidth map and the Spitzer 8m emission map for the region have been shown in Figs 13(b) and 13(c), respectively.
Here we notice that while the structures are contiguous in positionvelocity (pv) space, there is considerable fluctuation along the ve-locity axis for most of the pv maps, as one traverses along the line.These fluctuations are most prominent along the directions L1, L2, and L3.For L1, while there is spread in velocity along both sides of ∼100 km s −1 velocity; for L2 and L3, the contiguous structure is blueshifted and redshifted, respectively.Normally, such a shape would imply velocity gradients, and thus a longitudinal flow of matter along the structure.Along L2, the clump at ∼2.6 ′ corresponds to p3 position (also see Fig. 11), which is also associated with an ATLASGAL clump.Along L3, there is clumping at ∼0.4 ′ , 3.2 ′ , and 4.4 ′ , associated with positions p8, p7 (also an ATLASGAL clump), and p6.Since in the linewidth map (Fig. 13(b)), there is higher velocity dispersion at p3 (along L2), thus it is likely that the flow of matter is along the L2 molecular filament towards the gravitational well of p3.Similarly for L3, the relatively higher velocity dispersion at p7 and p6 would suggest longitudinal flow along the direction of the vector, though only p7 is associated with an ATLASGAL clump here.The structures along the vectors L4 and L5 appear to be elongated filamentary structures, and seem to be joining L3 near peaks p6 and p7.The gradient along L4 and L5 is much less than the case for L2 and L3.While the 13 CO(3-2) molecular emission along L5 seems fragmented, it appears as a prominent contiguous structure in Spitzer 8m emission, as well as at 24 and 70m (see Fig. 9).Finally, for L6, the pv slice shows two distinct clouds, as was also discussed earlier (see Fig. 12).The redshifted cloud (∼100 km s −1 ) appears to have a uniform velocity throughout, while the blueshifted cloud has clumps associated with positions p12 (at ∼2 ′ ) and p11 (at ∼3 ′ ).There is a small velocity gradient up to p12, and p11 shows a double-peaked structure indicating self-absorbtion, as was also seen in the spectrum (see Fig. 11).p11 also has an associated ATLASGAL clump and shows a large dispersion in Fig. 13(b).

Implications and Insights from 6.7 GHz masers
According to Deharveng et al. (2010), bubble N59 emerges as a faint structure spanning a diameter of approximately 20 parsecs.Our observations highlight a distinct concentration of star-forming activity primarily situated in the north eastern region of the bubble, marked as R1.Worth highlighting is the presence of five 6.7 GHz methanol masers within R1's bounds (refer to Fig. 14).Intriguingly, two of these maser detections align precisely with the 870 µm dust continuum clumps, exhibiting luminosities exceeding 10 4 L ⊙ .This discovery is consistent with the conclusions drawn by Bourke et al. (2005), suggesting a lower luminosity threshold of approximately 10 3 L ⊙ for a source to host 6.7 GHz methanol masers.Our photometric analysis identifies a massive Class I YSO (M∼8.7 M ⊙ ; hereafter HM1) in the proximity of maser M2.However, an angular offset of approximately 13 ′′ between their positions raises questions about the precise association between the maser and the identified YSO.As previously mentioned, the presence of 6.7 GHz masers typically indicates the early stages of massive star formation.This inference is reinforced by the observed dust temperature values within the maser hosting regions.For example, Fig. 15 displays the maps illustrating the column density and temperature distribution of R1 (from ViaLactea project, see Section 2.5).
Notably, the column density values within this region are found to be approximately on the order of 10 22 cm −2 .The top panel of Fig. 15 clearly depicts that the majority of YSOs are situated within areas characterized by significantly higher column densities (>10 22 cm −2 ).Similar column density values exceeding 10 22 cm −2 have been reported in numerous instances of massive star-forming regions, particularly during their very early stages of star formation (e.g.Rygl et al. 2010, Gieser et al. 2021 andreferences therein).
It is noteworthy that while not all peaks in the column density map coincide with peaks in the temperature map, the locations of methanol masers are consistently associated with regions exhibiting both high column density and elevated temperature (with temperatures > 18 K; Paulson & Pandian ( 2020)).Of particular interest is the column density and dust temperature observed towards HM1.In this context, the recorded values are approximately 3.6 × 10 22 cm −2 and 23.7 K, respectively.Similar dust temperature and column density values have been reported towards other massive star forming sites (Yu et al. 2019;Bally et al. 2010).Dust temperature serves as a dependable indicator of a source's evolutionary stage, with maser clouds generally displaying higher temperatures than less developed Infrared Dark Clouds but consistently lower temperatures than more evolved sources such as H ii regions (Giannetti et al. 2013) In the case of M1 and M2, as derived from Herschel dust temperature maps, the average dust temperature is approximately 21 K.This finding strongly suggests that these sites are positioned in the early stages of massive star formation, characterized by conditions conducive to the birth of massive stars (Paulson & Pandian 2020;Urquhart et al. 2014).
Furthermore, it is noteworthy that the ATLASGAL clump masses exhibit a range from approximately 6×10 2 to 10 3 M ⊙ .The mean clump mass within this sample stands at approximately 1300 M ⊙ , a sufficient mass to give rise to a star cluster that includes at least one massive star with a stellar mass of ≥10 M ⊙ (Urquhart et al. 2014).Particularly intriguing is the affiliation of both maser sources, M1 and M2, with the core region C1 (Refer Fig. 8), implying that this region may host massive stars or be actively involved in ongoing massive star formation processes.Notably, targeted ALMA archival data towards methanol maser M1 exists, and the analysis of the data holds promise for revealing compelling underlying physical phenomena.

Star formation scenarios in R1
The brightest condensation of bubble N59 is positioned towards its northern edge, as depicted in box B3 of Fig. 9.In B2 (refer to Fig. 9(c)), there is a potential indication of a blister or champagnelike flow structure of the H ii region.These structures commonly become evident when a massive star takes shape near the periphery of a molecular cloud.In such scenarios, the ionized gas either diffuses into the inter-cloud medium at the cloud's edge or follows a density gradient within the star's birth cloud.
It is possible that the blue-shifted tail seen at p6 (see Fig. 11 and the L3 pv map in Fig. 13) could indicate the molecular material pushed towards the observer by the champagne flow.
The high-mass star formation in B3 (refer Fig. 11) could be acting as a sink of material along vector L2.This is substantiated by the prominent velocity gradient observed along L2.Similarly, notable velocity gradients are also discernible in the vicinity of p6 and p7 along L3, L4, and L5.These gradients could suggest the directed flow of gas toward these sources.This phenomenon of gas flow along filamentary structures towards sources mirrors the longitudinal flow observed in filaments towards central hubs within the hub-filament system.The association of p7 with an ATLASGAL clump stands out, with its column density estimated to be approximately 10 22 cm −2 (Urquhart et al. 2022).Such regions with comparable column densities have been classified as "hubs" within the context of hub-filament systems, as per existing literature (Myers 2009;Kumar et al. 2020).Consequently, vectors L3, L4, and L5 seem to be integral components of a hub-filament system, wherein p6 and p7 are linked with hubs.This association implies that the region related to B2 could serve as an illustrative case for exploring diverse star formation models.However, delving deeper into understanding the molecular cloud evolution as it gives rise to stars requires further investigation.This involves employing high-resolution studies encompassing molecular transitions of different species, examinations of magnetic fields, and similar approaches to comprehensively grasp the intricate processes at play.
Particular interest is drawn to p11, highlighted in Fig. 11, due to its distinctive attributes.Notably, this cloud exhibits a primary velocity peak at approximately 75 km/s, in stark contrast to the peak velocity of 100 km/s observed in regions p1-p10.An interesting characteristic is the red-asymmetric line profile evident in 13 CO, which contrasts with the single-peaked profile detected in C 18 O, specifically near the central dip of the 13 CO profile.This "red profile," where the redshifted peak of a double-peaked profile is more pronounced for the optically thick line, indicates outflow motions, caused by the absorption of colder blueshifted gas in front of the hot core (e.g.Myers 1998;Liu et al. 2011, and the references therein).In their study of the high-mass star-forming complex G9.62+0.19,Liu et al. (2011) note that around younger cores, the outflow is more active and colder than the more developed ultra-compact H ii (UC H ii) regions.This results in a more pronounced "red profile."In the more mature UC H ii regions, where outflows are weaker, the surrounding gas is thermalized, and the temperature gradient toward the central star is more likely to cause a "blue profile", leading to a higher blue excess compared to the early stages of UC H ii or UC H ii precursors.Considering that the 6.7 GHz methanol maser signals the hot molecular core phase (and the very early stage of UC H ii region -e.g.Walsh et al. 1998;Minier et al. 2001;Ellingsen 2006) of highmass star formation, and given that site p11 hosts a 6.7 GHz methanol maser, the presence of the "red profile" at p11 is well-supported.
Additionally, the most massive star pinpointed in R1, an O9.5 star, is situated in proximity to p11 and exhibits detectable 6.7 GHz maser emission nearby.This locale presents itself as a promising candidate for subsequent high-resolution observations aimed at probing the underlying mechanisms of massive star formation.Further scrutiny reveals that while p11 and p12 exhibit velocities around 70-80 km/s, p13, situated in close proximity to these regions, displays a peak velocity of approximately 100 km/s.This disparity introduces an intriguing scenario of two distinct molecular clouds that, despite their spatial proximity, are differentiated in velocity space.While the possibility of star formation through cloud-cloud collapse might have been conceivable, the notable velocity difference exceeding 15 km/s between the cloud velocities necessitates the dismissal of this scenario.
R1 presents another captivating element with the existence of a bright-rimmed cloud (BRC)-like structure, found towards the site p6 in B2 (refer Figs 9(c) and 11).What makes this feature intriguing is its configuration, where the "head" of the BRC is oriented towards the bright dust condensation (Refer bottom panel of Fig. 15) located in the northern region of R1.Furthermore, B2 exhibits a blister-like morphology, with its ionized region expanding towards the east.This expansion suggests the potential presence of embedded massive star formation in the western part of B2, contributing to the formation of B2's H ii region.It is plausible that this embedded star formation in B2 is triggered by the adjacent H ii region of B3.This connection to ionized emission is significant because BRCs are known to originate from the influence of such ionization, often leading to triggered star formation (Choudhury et al. 2010;Fukuda et al. 2013;Panwar et al. 2020).Additionally, the detection of NVSS emission in B2 indicates the possibility of ongoing robust star formation, potentially indicating a scenario of radiation-driven implosion.Alternatively, these illuminated bright rims may serve as the delineations of regulardensity filaments shaped by ionizing radiation, all converging towards the central area.

General overview of star formation in and around bubble N59
Finally, we want to discuss the overarching bubble morphology of bubble N59 and the associated star formation processes within its vicinity.It is noteworthy to highlight that a considerable number of YSOs identified in Section 3.2 are located on the border of the bubble, as depicted in Fig. 8 (blue dots in right panel).These YSOs are likely influenced by the N59 bubble itself.Three prominent processes have been considered in the literature regarding star formation at the periphery of the H ii regions: collect and collapse (CC; Elmegreen & Lada 1977), radiation-driven implosion (RDI; Sandford et al. 1982), and cloud-cloud collision (CCC) (e.g., Habe & Ohta 1992;Haworth et al. 2015;Torii et al. 2017).In the CC process, a compressed layer of high-density neutral material forms between the ionization front and shock front, leading to regularly spaced fragments along the molecular ring or shell.RDI involves shocks driving into pre-existing density structures, forming stars characterized by cometary globules or optically bright rims.CCC can result in the rapid accumulation of cloud mass into a small volume, forming massive star-forming clumps within the shock interface.Bubble N59 exhibits a broken ring/shell-like structure of gas and dust surrounding the cluster of massive stars, with several YSOs at the edges, suggesting the possible impact of O and B type massive stars located in the R1 region or other areas within the bubble.Notably, Hanaoka et al. (2019) reported a total infrared (TIR) luminosity, log (L TIR /L ⊙ ) ∼ 6.6 , which may not solely originate from the massive stars within the R1 region.This high luminosity could be attributed to a cluster of massive stars located towards the centre of bubble N59.Additionally, the overall morphology indicates that the expansion of an older bubble might have triggered the star formation, leading to the development of the more recent sub-region R1.This further supports the potential presence of an older massive star population at the center of bubble N59.Although, no such massive star clusters were identified at the center of the bubble.
However, it's crucial to note that the massive stars responsible for the development of the bubble may not necessarily lie at the center of the bubble.Hanaoka et al. (2019) proposed a scenario where, for bubbles in the inner Galactic region, the heating sources are positioned nearly opposite to the broken sector directions.This offset could be attributed to the easier expansion of H ii regions toward the lower density direction, causing the lower density side of the shell to break.Another explanation by Hattori et al. (2016) links the offset of heating sources to the formation of massive stars on the boundary of collided clouds, resulting in a broken structure.Interestingly, bubble N59 hosts two Galactic filaments (Li et al. 2016) and a massive protostar near the broken edge (at Galactic coordinates  = 32.996• ;  = 0.042 • ) as catalogued by the RMS9 survey (Urquhart et al. 2016) at its edges.Hanaoka et al. (2019) suggest that if the CCC process is essential for broken bubbles in inner Galactic regions, it is expected that their massive stars would be substantially offset to the opposite side of the broken sector direction.
Further in-depth investigation is required to elucidate the exact mechanism by which bubble N59 formed, which is currently beyond the scope of this paper.In their statistical analysis of 65 IR bubbles, Deharveng et al. (2010) found that over a quarter of the bubbles triggered the formation of massive objects, either through collect and collapse or via the compression of pre-existing clouds (CCC).There- located in R1 within the bubble N59.We take its mean (4.66 kpc) and standard deviation (0.70 kpc) as R1's distance and its uncertainty.

Figure 2 .Figure 3 .
Figure 2. The Class I and Class II YSOs (filled black circles) identified using CC diagrams.Panels (a) and (b) represent Class I and Class II sources, respectively, obtained using the Phase 1 scheme (Gutermuth et al. 2009a).Panel (c) shows Class I and Class II sources in the NIR/IRAC CC plot with areas demarcated (by black dashed and solid lines) as per the Phase 2 scheme.The filled black circles mark the additionally identified Class 0/I sources.Grey dots within both shaded and unshaded regions across panels represent unclassified sources, and sources categorized as contaminants.
for YSO identification.The YSOs are classified into Class I (protostars with circumstellar disks and infalling envelops) and Class II (pre-main sequence stars with optically thick disks), based on their position in the [4.5]-[5.8]Vs [3.6]-[4.5]and [4.5]-[8.0]Vs [3.6]-[5.8]CC diagram, respectively.We obtained a total of 22 Class I and 31 Class II sources using this method, which are shown in Fig. 2 (a) and (b).Among this, 11 Class I and 8 Class II sources fall in R1.

Figure 4 .
Figure 4. J -H/H -K CC diagram for the region.The black curve marks the locus of the dwarfs fromHewitt et al. (2006) The dotted blue line is the Classical T Tauri star (CTTS) locus fromYasui et al. (2008).The dot-dashed black lines are the reddening vectors, drawn using the reddening laws fromRieke & Lebofsky (1985).Three regions-"F," "T," and "P"-mark the location of different classes of YSOs, with sources in each region colored in orange, yellow and brown, respectively.The horizontal gray line has been drawn at J -H = 0.89.All points are in the Mauna Kea Observatory (MKO) system.

Figure 5 .
Figure 5. K/H -K NIR CM-D for bubble N59.Nearly vertical solid lines are the loci of ZAMS stars reddened by A V = 0, 1, 10, 20 and 30 mag.Parallel, slanting lines indicate the reddening vectors for respective spectral types.The small star symbols represent the total YSOs in bubble N59, with those belonging to region R1 colored blue.All points are in the Mauna Kea Observatory (MKO) system.
2 criterion as defined by Equation 1.This subset encompasses 7 Class I and 7 Class II sources.The visual representation in Fig. 6 depicts these YSOs, displaying their respective mass and age properties, overlaid onto the PACS 70 µm image.The orange and green filled circles represent Class I and Class II sources, respectively.The main parameters obtained from the SED analysis are given in Table 2.The clump masses of these sources range from 2.77 -14.20 M ⊙ , and their ages span from 2.13 ×10 4 − 7.46 ×10 5 years.Notably, four of our Class I sources exhibit high mass values, measuring at 14.20, 8.74, and 8.74 and 8.70 M ⊙ respectively.It is worth mentioning that the majority of YSOs observed in region R1 are characterized by relatively young ages, all of them falling below 1 Myr.

Figure 6 .
Figure 6.The mass and age (in units of M ⊙ and year, respectively) of the YSOs in region R1 overlaid on PACS 70 µm image, are displayed in the top and bottom panels.The age is displayed in the log10 scale.Class I and Class II sources are marked in orange and green, respectively.

Figure 7 .
Figure 7. Histogram of the MST branch lengths used for critical length analyses of the YSOs.

Figure 8 .
Figure 8. Cluster identification using MST.The left panel shows the total number of sources that are connected by branches with lengths, , shorter than the value on the -axis.The steeper brown line is the linear fit to the data below 90 ′′ , and the other brown line is a similar fit to the data beyond 134 ′′ (vertical dashed-pink line).The cut-off length ( c ) is opted as the branch length at the intersection of these lines (vertical dashed-black line).The image in the right panel shows the MST overlaid on Spitzer 8.0 µm map.The stars -marked in blue -are connected by solid red ( <  c ), solid light-brown ( c <  < 134 ′′ ) and dotted white curves ( > 134 ′′ ), that together form the branches of the tree.The dashed-yellow and dashed-cyan polygons are convex hulls demarcating the active and core regions, respectively.The two core regions are shaded in blue and are labelled as C1 and C2.The green square shows the region R1.

Figure 10 .
Figure 10.VLA 8.49 GHz image of the B1 box region in Fig. 9.The symbols and red contours are same as in Fig. 9. Blue contours are VLA 8.49 GHz image contours at 0.005 and 0.015 Jy/beam.The radio beam is shown in bottom left corner.Here we can see the r4 source resolved into N and S components.

Figure 15 .
Figure 15.Top and bottom panels illustrate column density and temperature maps for region R1.The column density map is in the unit of 10 20 cm −2 .Class I and Class II YSOs are marked by orange and green filled circles, respectively, while labels M1-M5 indicate the positions of 6.7 GHz methanol masers.

Figure B1 .
Figure B1.The left panel shows the face-on stellar overdensity map for the Milky Way based on the Gaia DR3 (Poggio et al. 2021), where the overdensity decreases from red to blue.The Saggitarius and Scutum spiral arm models from Reid et al. (2019) are shown as blue and red curves, respectively.The Sun and the Galactic centre are shown as the yellow and black crosses, respectively.The Gaia stars located within R1 in projection are shown as grey points.The right panel shows the probability density function of the trigonometric distances of the stars with clustering probability above 0.8.The blue and red vertical lines mark the distances where the Saggitarius and Scutum arm intersect our stellar trail in the left panel.The black curve is the best-fit multi-gaussian model for this distribution; the three Gaussians constituting the black model are shown in blue, red, and green.

Table 1 .
Catalog of various surveys endorsed for this study, covering a wide range of wavelengths from near-infrared to radio frequencies.

Table 2 .
Main parameters from SED analysis.

Table 3 .
H II Regions