Star Formation in the H ii Region Sh2-87: Evidence of Global Hierarchical Collapse

We present a detailed study of the Sh2-87 H ii region using a multiwavelength data set in optical to radio bands. A Herschel column density map revealed that the host cloud is filamentary in nature, and together they formed a central dense hub. The extinction map generated using near-infrared photometric data also signifies the nonuniform distribution of the cloud and reveals its filamentary nature. We estimated a sizable variable extinction over the region up to A V = 34.4 mag, with an average value of A V = 3.4 mag. Using the various infrared color–color criteria, we identified 13 Class I and 202 Class II young stellar objects (YSOs) and 22 Hα-emitting sources toward this region. Further analysis showed that the cluster is mainly composed of low-mass YSOs with a typical age of ∼3 Myr having masses in the range of 0.1–6.0 M ⊙. The identified evolved YSOs (i.e., Class II YSOs) are primarily distributed along the filaments and in the outer parts of the cloud, while the recent star formation, inferred by the presence of Class I YSOs, ionized gas, and star-forming clumps, is observed in the hub region. The overall star formation scenario in the Sh2-87 region resembles the global hierarchical collapse model of star formation, where younger massive star formation activity is expected at the central hub along with the distribution of evolved low-mass YSOs in the filaments and the outer parts of the cloud.


Introduction
Stars are formed within dense cores of the molecular clouds, and these sites are obscured with immense interstellar gas and dust (Lada & Adams 1992).The interstellar medium (ISM) often appears as a network of filamentary structures (Schneider & Elmegreen 1979;André et al. 2010André et al. , 2014)).The universality of filamentary structures in the ISM was established with the advent of observations from Herschel Space Observatory (André et al. 2010).Several facts revealed that these filaments seem to play an important role in the star formation process (Myers 2009;Dale & Bonnell 2011;Schneider et al. 2012).In star-forming clouds, these filaments often converge to form a common junction of filaments known as a hub-filament system (HFS) with a comparatively higher column density ((N(H 2 )) ∼ 10 22 cm −2 ) but have a low aspect ratio (Myers 2009;Peretto et al. 2014;Chen et al. 2019;Kumar et al. 2020).Such HFSs are regarded as one of the potential sites for the formation of massive stars (8 M e ) and young star clusters (Henshaw et al. 2014;Dewangan et al. 2020;Zhou et al. 2023).Hence, the HFSs serve as the best laboratory to examine the early conditions of star formation (e.g., Dewangan et al. 2017;Kumar et al. 2020;Bhadari et al. 2022;Chung et al. 2023).
Different theories exist to explain the formation of such HFSs.The prevalent HFS is mainly formed either via the inertial-inflow scenario (Padoan et al. 2020) or by a global hierarchical collapse (GHC; Vázquez-Semadeni et al. 2017, 2019).Even though the detailed mechanism of HFS formation in these two theories is entirely different, both mechanisms highlight the longitudinal stream of the materials into the hub along the filaments.According to the inertialinflow model, most of the final stellar mass is directed toward the accreting star by the random velocity field from large scales, but it is uninfluenced by the stellar gravity in its path toward the star (Padoan et al. 2020).As per the inertial-inflow model, a large-scale turbulent compression occurs, where the filaments serve as the location of the intersection of post-shock layers (Zhou et al. 2023).On the other hand, the GHC model predicts that fragmentation takes place in sufficiently dense filaments, and they provide the dense massive core with the blend of gas and stars (Vázquez-Semadeni et al. 2019).According to Gómez & Vázquez-Semadeni (2014), in the case of GHC, the filaments construct part of the large-scale gravitational collapse and guide the flow of gas into the cores.In contrast to the inertial-inflow scenario, the GHC model predicts the cloud scale flow of the matter along the filaments dominated by gravitational energy.Here the collapse is regarded as hierarchical because it comprises a number of small-scale collapses within the bounds of larger-scale ones (Vázquez-Semadeni et al. 2017).In recent statistical and observational multiwavelength studies, several authors (Zhou et al. 2022;Liu et al. 2023;Ma et al. 2023;Zhou et al. 2023) recognized the gravitational GHC model as the mode of star formation.To appraise the existing models of HFS, the study of the physical environment and molecular gas associated with the star-forming region is of primary importance.Thus, using a wealth of multiwavelength data we explored the Sh2-87 (hereafter S87) H II region region that appeared with an embedded HFS.Our study aims to examine the detailed star formation scenario in this region and to investigate if the formation mechanism resemblances to any of the models discussed above.
The S87 region ( ( ) ) is recognized as an optical H II nebula of size 10′ by Sharpless (1959) with a faint and diffused nebulosity.The region is a member of the Vulpecula OB1 association (Vul OB1; Billot et al. 2010), along with other H II regions such as Sh2-86 and Sh2-88.Multiple distance estimations are available for this region.On the basis of radio observations, Bally & Predmore (1983) estimated the distance to the S87 region as 2.0 kpc.From the Galactic rotation curve, Clemens (1985) ascertained the distance of S87 as 2.1 kpc.Using photometric studies of the bright sources associated with S87, Racine (1968) determined E(B − V ) = 1.2 ± 0.1 and a distance of 2.3 ± 0.1 kpc.Crampton et al. (1978) showed that although the star BD +24°3880 is less central than BD +24°3873, it serves as the main exciting source of the region (V = 10.80 mag), and they estimated the distance as 2.3 kpc.From the photometric and spectroscopic studies, Forbes (1989) identified star BD +24°3880 as having spectral class B0.5V and also mentioned this as the ionizing source; they took the color excess E(B − V ) = 0.81 and slightly overestimated the distance as 2.6 kpc.In the present work, we adopted an average distance of 2.3 kpc, which was also previously adopted in Billot et al. (2010).
Several observational studies have been performed on this object in the past.For example, using radio wavelength observations, Felli et al. (1978) detected an ultracompact (UC) H II region in S87.Later, from the Very Large Array (VLA) observation, Bally & Predmore (1983) recognized S87 as a compact H II region in a state of photoionization equilibrium.Henkel et al. (1986) detected two highly variable 22 GHz water maser continua in the S87.Barsony (1989) proclaimed that S87 is a massive pre-main-sequence (PMS) object still embedded in the host molecular cloud, but is crumbling up its surroundings with strong stellar winds; hence a massive biconical outflow is associated with S87.With the aid of near-infrared (NIR) imaging of S87, Chen et al. (2003) showed the presence of two bright nebulae in the eastern and western parts of this region, denoted by S87E and S87W, respectively.Using the submillimeter observations at 850 μm, Xue & Wu (2008) detected three submillimeter clumps, SMM1, SMM2, and SMM3, which lie along an axis from southwest to northwest within S87.They also reported an outflow activity associated with S87E and proposed that the clump is at a very young stage (3-4 Myr).All these observational facts reveal that the region S87 is an active site of star formation.
Figure 1 shows a three-color-composite image of the region constructed using Wide-field Infrared Survey Explorer (WISE) bands (red: 22 μm; green: 12 μm; blue: 3.4 μm).This midinfrared (MIR) three-color image clearly depicts the diffused nebulosity associated with this region.The distribution of the radio continuum emission is also shown in the figure .A UC H II region was previously detected in this region by Bally & Predmore (1983) and Kurtz et al. (1994).Even though the region has been studied from different perspectives and strong evidence of recent star formation activity has been found, no study has been conducted so far to probe the embedded filaments and their role in the star formation process associated with this H II region.This particular aspect is the primary aim of this study, where we explored the kinematic and physical properties of the molecular cloud, identified an HFS and embedded young stellar population, and compared the overall scenario with the recently proposed models for the formation of an HFS.This paper is organized as follows.In Section 2, we presented the archival data sets utilized in this work.Sections 3 and 4 present the results and the discussion.The findings of the work are summarized with a conclusion in Section 5.In addition, a multiwavelength catalog for the young stellar objects (YSOs) is presented in the Appendix.

Optical Data: Pan-STARRS 1
We utilized the optical/NIR data from the Panoramic Survey Telescope and Rapid Response System 1 (Pan-STARRS1, or PS1; Chambers et al. 2016).The PS1 system carried out the observations with the five filters g P1 , r P1 , i P1 , z P1 , and y P1 distributed over wavelengths ranging from 400 to 1000 nm and centered on 481, 617, 752, 866, and 962 nm, respectively (Tonry et al. 2012).

IPHAS Data
The INT/WFC Photometric Hα Survey (IPHAS; Drew et al. 2005) encompasses a 1800 deg2 area of sky in the northern Galactic plane, and it uses the narrowband Hα filter (656.8 nm) and broadband Sloan r (625.4 nm) and i (774.3 nm) filters for observations.We acquired the photometric data for r, i, and Hα bands from the IPHAS data release 2 (DR2) catalog (Barentsen et al. 2014).We selected the sources having brightness 13 < r < 20 and photometric uncertainty of <0.1 mag in DR2.

Astrometric Catalog
The photogeometric distances of the point sources were obtained from Bailer-Jones et al. (2021) estimated using the astrometric data from the Gaia Early Data Release 3 (Gaia EDR3; Gaia Collaboration et al. 2021).These photogeometric distances were estimated on the basis of the parallax, color, and magnitudes of each source.The UC H II region detected by Bally & Predmore (1983) and Kurtz et al. (1994) is also marked by a cyan hexagon.

NIR Point-source Catalogs
The NIR J-, H-, and K-band photometric data were obtained from the sixth data release of the UKIRT Infrared Deep Sky Survey (UKIDSS DR6PLUS) of the Galactic Plane Survey (GPS; Lawrence et al. 2007).Magnitudes for a few saturated sources in the UKIDSS catalog were obtained from the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006).The angular resolutions of UKIDSS and 2MASS images are 0 8 and 2 0, respectively.For the selection of reliable sources from the UKIDSS GPS catalog, we followed the criteria based on the Structured Query Language (SQL6; Lucas et al. 2008;Dewangan et al. 2015).We picked those sources with a magnitude error of 0.1 or less in each band and having a signalto-noise ratio (S/N)  10.We have discarded unreliable sources such as nonstellar sources and saturated sources using the flag listed in the catalog.The selected sources from UKIDSS are fainter than J = 12.5 mag, H = 11.6 mag, and K = 10.2 mag, and for the sources brighter than these limits, we supplemented the sources with 2MASS photometry.

MIR Point-source Catalogs
We obtained the MIR photometric magnitudes of the point objects from the Spitzer Galactic Legacy Infrared Mid-Plane Survey Extraordinaire (GLIMPSE; Benjamin et al. 2003), which have imaged the region with the Infrared Array Camera (IRAC; Fazio et al. 2004) in the four MIR bands centered at 3.6, 4.5, 5.8, and 8.0 μm, and the corresponding angular resolution of the images is about 2″.For good-quality photometric data, we used the sources with the uncertainty σ < 0.2 mag in all of the four IRAC bands from the GLIMPSE I and II and 3D surveys.
In addition, we utilized the photometric MIR photometric magnitudes of the point sources from WISE (Wright et al. 2010) to identify sources that might be saturated/undetected in the IRAC bands.WISE has scanned the sky in the 3.4, 4.6, 12, and 22 μm infrared bands and has an angular resolution of 6″ for the first three bands and 12″ for the last band.We have taken the MIR photometric magnitudes of the point sources from the first three bands ([3.4], [4.6], [12]) of the ALLWISE catalog.To avoid contaminated sources, we utilized the different conditions based on the S/N, χ2 , and uncertainty as mentioned in Section 3 of Koenig & Leisawitz (2014).

Submillimeter Data
The CO molecular line data were obtained from the James Clerk Maxwell Telescope (JCMT; Buckle et al. 2009) archive. 1n this work, we used 12 CO (3-2) data (Program ID: M14AU09) obtained using the Heterodyne Array Receiver Programme (Buckle et al. 2009) having a spatial resolution of 14″.The third axis of the spectral cube is in barycentric frequency in units of GHz; to convert this frequency axis into a radio velocity axis with the standard of rest set to LSRK, we have used the Starlink software's CONVERT and KAPPA packages (Currie et al. 2014).Finally, the spectral map is regridded into 0.5 km s −1 channels and used for further analysis.The CO (J = 3-2) transition traces relatively warmer (∼33 K) and denser gas (2.1 × 10 4 cm −3 ) that is associated directly with star formation (Wilson et al. 2012).
For the dust continuum emission in our target region, we also utilized 850 μm continuum data observed from the Submillimetre Common User Bolometer Array 2 (SCUBA-2; Dempsey et al. 2013) installed on JCMT.The mean rms noise level and beam FWHM of SCUBA-2 at 850 μm are ∼2.9 mJy beam −1 and 13″, respectively.
We also utilized the publicly available 2 Herschel H 2 column density map generated from the 160, 250, 350, and 500 μm wavelength bands.The map has a spatial resolution of ∼12″ (Molinari et al. 2010), and it is generated following the Bayesian point-process procedure (PPMAP; Marsh et al. 2015).

Radio Continuum Data
We obtained the radio continuum map at 1.4 GHz for our target region from the NRAO VLA Sky Survey (NVSS; Condon et al. 1998) archive.The sensitivity and the resolution of this radio map are ∼0.45mJy beam −1 and ∼45″, respectively.We also used the radio continuum data at 5 GHz from the Coordinated Radio and Infrared Survey for High-Mass Star Formation (CORNISH; Hoare et al. 2012).The CORNISH data have a sensitivity of ∼0.3 mJy beam −1 and an angular resolution of 1 5.

Compilation of Multiwavelength Catalog
To generate the final multiwavelength catalog of the point sources, we cross-matched all the data sets mentioned above.We utilized the software TOPCAT (Taylor 2005) to match the coordinates of the stars from different catalogs for a matching radius of 1″.In a few cases, multiple matches were found, and for those we selected the nearest matching sources.Table 1 in the Appendix represents the catalog of the YSOs selected from the MIR and NIR data, whereas Table 2 in the Appendix is the photometric catalog for the Hα emission sources selected from IPHAS data.In catalogs, distances of the objects were obtained from Bailer-Jones et al. (2021).

Identification of Filaments and the Hub
The Herschel column density (N(H 2 )) map of the region is presented in Figure 2(a).Based on the cloud morphology seen in the column density map, we select a 20′ × 20′ area centering on the S87 UC H II region for our analysis.Filaments were identified in the N(H 2 ) map with the aid of the Python-based FILLFINDER 3 algorithm (Koch & Rosolowsky 2015).FILL-FINDER can recognize all of the probable filaments in the input map by the technique of determining their corresponding skeletons through the medial axis transform method.We traced multiple filaments in this region, and they are marked in Figure 2(b).For identification of these filaments we used a global threshold of N(H 2 ) ∼ 3 × 10 21 cm −2 for true emission.As can be seen in Figure 2(b), these skeletons trace the dense gas structures in the column density map.The typical length of these filaments is a few parsecs.The central hub region is defined with a critical H 2 column density of 10 22 cm −2 (Myers 2009).The identified filaments clearly seem to emerge from the central hub region, and overall they reveal the existence of hub-filament morphology as reported in several previous studies (e.g., Myers 2009;Schneider et al. 2012;Baug et al. 2015;Dewangan et al. 2015).

Distribution and Dynamics of Molecular Gas
The region is located at a local standard of rest velocity (V lsr ∼ 23.6 km s −1 ), and the molecular gas of this region is traced at a velocity range from 20.5 to 25.5 km s −1 .To examine the large-scale molecular cloud structure, we generated the integrated intensity map (moment 0 map) using 12 CO (J = 3-2) line data for a velocity range of 20.5-25.5 km s −1 .The corresponding map is shown in Figure 3(a).The mean integrated intensity for this cloud is estimated as 19.1 K km s −1 .The central hub region appears strong in CO emission with integrated intensity 80 K km s −1 .For the velocity distribution of the molecular cloud, we also constructed the velocity distribution map (moment 1 map) of the region (see Figure 3(b)).This map signifies the measure of the intensity-weighted average velocity of the emitting gas.A strong variation of velocity can be noted in the velocity distribution map, specifically where the filaments are identified.The velocities of the gas toward the filaments typically have larger values (24.0 km s −1 ) compared to the central hub (∼22.5 km s −1 ).This particular scenario signifies the gradient in the velocity of molecular gas along the filaments.For further investigation of the gas dynamics, we constructed the velocity dispersion map (derived from moment 2 map) of the region (see Figure 3(c)).A strong velocity dispersion (1.5 km s −1 ) signifying the presence of turbulent gas can be noted toward the central hub of the S87 complex.

Extinction Map
Young star clusters are rich with dust that diminishes the background stellar contamination.Besides, the nonuniform distribution of molecular clouds and dust further makes it more cumbersome to separate the background stars from the young cluster population (Gutermuth et al. 2008).To filtrate these field stars from the young stars embedded in the clusters, we produced the extinction map for our target region.NIR colors are extremely useful for the direct estimation of the extinction of individual stars, arising from interstellar and circumstellar material (Straw & Hyland 1989).As the infrared color excess is proportional to the dust column density, the extinction map generated from the NIR data of the background main-sequence (MS) stars can precisely unveil the distribution of the dust in the cloud.For the mapping of our present region S87, we utilized the NIR UKIDSS catalog.As extinction becomes lesser at longer wavelengths, deeper clouds can be better probed in H and K bands.
We determined the A K values from the (H − K ) colors abiding by the grid method of Gutermuth et al. (2005).Concisely, we can describe the process in the following way.We have divided the region into small uniform grids of size 5″ × 5″, considered the 20 nearest sources from the center of each grid for the estimation of the mean and standard deviation of (H − K ) colors of each grid, and discarded the sources with (H − K ) colors deviating over 3σ from the mean value.To avoid the contribution of the young embedded stars with circumstellar disks, we took only the nonexcess infrared sources.We obtained the value of A K of each grid from the mean (H − K ) colors by the reddening law Flaherty et al. (2007), where we took the intrinsic colors of the stars as (H − K ) int = 0.2 (Allen et al. 2008).We finally converted the A K value to A V by using the extinction ratio A V /A K = 0.112 from Rieke & Lebofsky (1985).Note that this method underestimates the extinction, specifically in the densest parts of the clouds, where fewer numbers of background stars are typically detected.To compensate for this issue, we generated a hybrid extinction map using the NIR data and the JCMT SCUBA-2 850 μm image following the method devised by Gutermuth et al. (2005).The A K /F 850 value was computed from 850 μm dust continuum fluxes using the formula A K /F 850 = 0.272 as reported in Gutermuth et al. (2005).The generated extinction map from dust continuum images was then resampled to the same grid size as the previously derived NIR extinction map.To combine the extinction map, we replace the A K values in the NIR extinction map with the A K values in the millimeter-band extinction map for those grids where millimeter-band extinction is 30% higher than the value derived using NIR data.This hybrid extinction map is used in further analyses.
Figure 4 shows the derived extinction map of the S87 region, and the nonuniform distribution of dust (hence, gas) can be noted in the extinction map.The extinction varies in the wide range from A V = 0.1 to 34.4 mag, with an average extinction of A V = 3.4 mag.This extinction map has a resolution of 5″, and sensitivity is down to A V = 34.4mag.The H 2 column density contours obtained from the Herschel column density (N(H 2 )) map are also overlaid on the map.The extinction map is strongly correlated with the column density map.

Identification of the Young Population
For examination of the overall star formation activity in a star-forming region, it is important to identify the YSOs of the region.As the embedded cluster members are often obscured in the optical bands owing to high extinction, IR data are useful for the identification and characterization of those young cluster members.For identification of the young members in the S87 region, we constructed NIR and MIR color-color diagrams (CCDs) from the magnitudes of point sources reported in UKIDSS, 2MASS, IRAC, and WISE catalogs for a ¢ ´¢ 20 20 area around the center of the cluster.Detailed identification methods are discussed in the following paragraphs.As classical T Tauri stars (CTTSs) appear with distinct For further identification of possible YSOs from the ALLWISE catalog, we have followed the method mentioned in Koenig & Leisawitz (2014).Note that, despite the poor spatial resolution, WISE band magnitudes might be useful in detecting Class I YSOs, as it includes a longer wavelength band (i.e., 12 μm in our case) than Spitzer-IRAC.Before applying the YSO classification scheme, we removed the extragalactic contaminations (star formation galaxies and AGNs).From the (W1-W2)/(W2-W3) CCD (Figure 5(b)) we obtained one Class I and six Class II sources.All the identified Class II sources are already detected in the IRAC CCD, but the Class I source is an entirely new detection.2. MIR and NIR CCD: Due to the nebulosity in the IRAC 5.8 and 8.0 μm, several sources may remain undetected in these two bands, while those might be detected in the IRAC 3.6 and 4.5 μm and NIR K bands.For this set of sources, we applied the phase 2 classification scheme of Gutermuth et al. (2009).Following this method, we first eliminated the field-star population from the dereddened colors ([[3.6]-[4.5]])0 and [K − [3.6]] 0 of each star using our already-generated K-band extinction map (Section 3.3.1).For each source, we took the extinction value of the nearest grid and the color excess ratio was used from Indebetouw et al. (2005).Applying the different conditions in the dereddened color space, we first discarded all extragalactic contamination, and from the classification conditions we identified 2 additional Class I and 42 additional Class II sources as shown in Figure 5(c).3. NIR CCD: NIR observations are crucial for the investigation of large YSO populations, as they can probe the surrounding circumstellar environments of the YSOs.In addition to the MIR data, we also used the infrared J, H , and K photometric data from the UKIDSS and 2MASS catalog for the inclusion of new YSOs.The UKIDSS data were first converted into the 2MASS system according to the transformation equation taken from Wegg et al. (2015), and then they were further converted to the CIT system following the relations of Carpenter (2001).The (J − H)/(H − K ) CCD for our target region is shown in Figure 5(d).To recognize the new YSOs on the basis of their position in the (J − H)/(H − K ) color-color space, we followed the method devised in Lada & Adams (1992), Ojha et al. (2004), andPanja et al. (2022).
To separate out the YSOs from the field-star population, we also analyzed a control field ( ) a = of the same area ( ¢ ´¢ 20 20) and identical photometric depth.From the comparison of the color-color distribution between this control field and our target region, we settled the (J − H) and (H − K ) color cutoff at 1.5 and 0.9 mag, respectively, as most of the stars in the control field have colors below this limit.Hence, the T and P regions of our CCD are almost free from field-star contamination beyond this limit.Using this method, we have 133 candidate YSOs, 10 belonging to the P region (i.e., Class I) and 123 belonging to the T region (i.e., Class II).
For a unified catalog of YSOs, we cross-matched the YSOs catalogs generated using the abovementioned methodology.If a single source was identified in multiple schemes, then for the class of the YSO priority was given to MIR CCD, then MIR-NIR CCD, and finally NIR CCD.Accordingly, we identified 13 Class I and 202 Class II sources.

Identification of Hα Emission Stars
IPHAS (r − i/r − Hα) CCD is useful to identify the Hα emission sources.Most of the Hα emitters, as observed by the IPHAS, may belong to the group of CTTSs, and their Hα emission arises as a result of the presence of hot, infalling gas that is accreted onto the star from a circumstellar disk (Barentsen et al. 2014).The other stellar objects, such as the massive stars (e.g., Be stars; Wolf-Rayet stars), evolved intermediate-mass stars (e.g., Mira variables, unresolved planetary nebulae), and interacting binaries (e.g., cataclysmic variables, symbiotic stars), also show significant Hα emission (Barentsen et al. 2011).
For the identification of Hα emission sources, we followed the method of Barentsen et al. (2014).Figure 6 shows the (r − i/r − Hα) CCD toward our target region S87.Using the unreddened MS emission-line strength of equivalent width (EW) = −10 Å as the CTTS threshold, we identified a total of 22 Hα emission-line stars.Additionally, the selected sources are located above the dashed blue line at the level of 3σ, i.e., the distance between the selected objects and the dashed blue line is larger than three times the average uncertainty in their (r − H α) color (Barentsen et al. 2011).Among the 22 Hαemitting objects, 7 sources are already classified as Class II sources and 1 source is classified as Class I in our YSO classification schemes.The photometric catalog of the Hαemitting sources is given in Table 2 in the Appendix.

Age of the YSOs
The optical color-magnitude diagram (CMD) with PMS isochrones for different ages serves as an authentic tool to estimate the age distribution of the YSOs.On the basis of the specific location of YSOs in the optical CMD, the quantitative ages of the YSOs are generally determined.This method works with an assumption of the negligible contribution of emission from the envelopes and disks in the optical bands.For the age distribution of the YSOs in our target region S87, we utilized PS1 photometry.Figure 7(a) shows the distribution of the YSOs and Hα sources in the CMD (g P1 − y P1 ) versus g P1 with the same symbols as used in the previous CCDs.The zero-age MS (ZAMS) and PMS isochrones for different ages 0.1, 0.5, 1, 3, 5, 10, and 20 Myr for solar metallicity are adapted from Marigo et al. (2017).All isochrones are modified for the cluster distance of 2.3 kpc (taken from Billot et al. 2010), the average visual reddening A V ∼ 3.4 mag (Section 3.3.1),and E(g P1 − y P1 ) = 2.4 mag (Schlafly et al. 2016;Green et al. 2018), and the extinction vectors for the PS1 bands were obtained from Green et al. (2018).Concomitantly, to study the mass distribution of the YSOs, we also plotted the grid of evolutionary tracks for different masses in the range of 0.5-4.0M e , shown by the dotted curved lines.The figure shows a notable age spread in the YSO distribution from 0.1 to  [3.6]] 0 CCD generated using cross-matched IRAC and 2MASS catalogs.The reddening vector corresponding to A K = 2 mag is also shown.(d) NIR (J − H)/(H − K ) CCD of the point sources in the S87 region.The green and red solid lines represent the locus of the unreddened dwarfs and giants adopted from Cohen et al. (1981), and the black solid line is the locus of the CTTSs (Meyer et al. 1997).Three parallel reddening vectors are shown with three black dotted lines, drawn from the base of the dwarf, the tip of the giant branches, and the top of the CTTS loci, respectively.The black arrow represents the reddening vector A K = 5 mag.20 Myr, whereas a majority of the YSOs reclined around 3 Myr and a few stars are located below 0.1 Myr isochrone.The phenomena of age spread are also observed in other young clusters (Ojha et al. 2011;Panwar et al. 2017) and arise possibly as a result of the circumstellar extinction, differential reddening, variability, and photometric uncertainty (Jose et al. 2017).However, in this CMD, the isochrones are well separated, which indicates the smaller effect of binarity in the estimation of ages of the cluster (Panwar et al. 2017).Finally, from this CMD, the average age of the YSOs is estimated as 3 Myr.However, the shortage of good photometric data for major YSOs makes this plot not a very efficient technique for the assertion of the age distribution of this region.

Mass Spectrum of the YSOs
To investigate the masses of the candidate YSOs toward our target region S87, we utilized the NIR (J − H)/J CMD.As the circumstellar disk of YSOs contributes to the longer wavelengths, such as H and K bands, than the J band, the excess emission can cause the overestimation of the masses (Ojha et al. 2004;Panja et al. 2020).Hence, the (J − H)/J CMD is the best choice for the mass spectrum study.Figure 7(b) exhibits the distribution of the selected YSOs in the (J − H)/J color-magnitude spaces.The black solid line is the loci of ZAMS acquired from Girardi et al. (2002), and the blue, red, and magenta solid lines signify different PMS isochrones of ages 1, 3, and 20 Myr, respectively, adopted from Marigo et al. (2017) and adjusted for the cluster distance of 2.3 kpc.To scrutinize the mass spectrum of the YSOs, we considered the average age of 3 Myr (see Section 3.3.4).The dashed parallel slanted lines are the reddening vectors for the different masses drawn from the 3 Myr isochrone.
Based on the comparison of the spread of the YSOs with the theoretical evolutionary isochrones, the average mass of the YSOs is estimated.Here the YSOs indicate a broad spread of the masses ranging from 0.1 to 6.0 M e .It is evident from the figure that the low-mass ending parts of the different isochrones are in close proximity; hence, a change of average age of ∼1-2 Myr will not affect the masses on a large scale.It is also noted that some of the massive stars are also present, and they are placed in the high-extinction (A V ∼ 10-30 mag) zone.Even for the massive stars, this method will give more uncertain results because, in place of the 3 Myr isochrone, if the masses have been estimated from the ZAMS isochrone, the stellar mass will change drastically.In addition, the colors of the YSOs have a huge variation that may arise as a result of the assembled effect of variable extinction, but they might not represent true (J − H) colors of the YSOs (Ojha et al. 2011).In addition, this method of estimation of stellar masses has limitations owing to the uncertainty of the ages and different distances of the different objects (Panja et al. 2020).Furthermore, this method also has other restrictions, arising from the use of several PMS isochrones, the existence of variable extinction, binary components, the circumstellar envelope, etc.Thus, the estimated masses of YSOs could be assumed as representative values and provide an average stellar mass of the cluster members.

Discussion
As is discussed in the previous sections, this region is associated with multiple filaments of length a few parsecs merging to a common overdense hub having H 2 column density >10 22 cm −2 .This is the first report of the identification of a hub-filament structure toward the S87 region.The hub is also associated with strong ionizing emission depicting the presence of a massive star(s).The velocity distribution and velocity dispersion maps of CO line data show evidence of velocity gradient along the filaments.Note that the presence of velocity gradients along the filaments is typically regarded as an indication of gas flow along the filaments (Kirk et al. 2013;Nakamura et al. 2014;Chen et al. 2019;Zhou et al. 2022).In their particle hydrodynamics numerical simulation studies, Gómez & Vázquez-Semadeni (2014) showed that these filamentary structures make river-like arrangements by which molecular gas and materials accelerate from the surrounding cloud environment toward the dense core and also nourish the star-forming cores and protoclusters in the central hub region.Filaments are typically regarded as the channels for the funnelling of gas in the direction of the central potential hub due to the effect of gravitational attraction (Gómez & Vázquez-Semadeni 2014;Vázquez-Semadeni et al. 2019;Ma et al. 2023).
In order to study the global star formation scenario, we conducted the analysis of the stellar members of the target region in ¢ ´¢ 20 20 and identified candidate YSOs.An extinction map generated from the color excess of the background stars reveals the distribution of the dust, and the corresponding extinction map is strongly correlated to the column density map.In addition, the region also shows an ongoing star formation activity as revealed by the identification of YSOs with an average age of 3 Myr and masses ranging from 0.1 to 6.0 M e .The spatial distribution of the YSOs helps    Xue & Wu (2008).The SCUBA contours (green) at 0.35, 1, 2, 5, 7, and 10 mJy arcsec −2 are also overlaid in this image, while the background rms noise level is at 0.07 mJy arcsec −2 .The presence of ionized gas is also marked in white contours of CORNISH 5 GHz data, drawn at the [5σ, 10σ, 20σ] levels (where σ = 0.26 mJy beam −1 ).
estimated the ratio to vary between 0.13 and 0.54.The Class I population generally corresponds to the earlier phases of the star formation in comparison to the Class II population; hence, the higher value of the ratio Class I/II implies a young cluster (Chavarría et al. 2008;Beerer et al. 2010).We calculated the ratio of Class I/II, which is ∼0.3 for the central hub region of S87, and for the outside of the hub, the ratio is negligible.Thus, for S87, this ratio of Class I/ II suggests that the hub region is in the early stages of star formation and is younger compared to the outer parts of the cloud.
From the submillimeter observations, Xue & Wu (2008) identified three clumps SMM1, SMM2, and SMM3 having masses 210, 140, and 110 M e and temperatures 41, 21, and 24 K, respectively.The SCUBA 850 μm contours were overlaid in the zoomed-in view of the HFS, and we marked these three clumps SMM1, SMM2, and SMM3 in Figure 8(b).These three clumps are located in the hub region.The region SMM1 is also associated with the UC H II , previously reported by Felli et al. (1978).High-resolution radio continuum contours obtained from the CORNISH survey are also overplotted on the image, and these contours have a strong overlap with the SMM1.Xue & Wu (2008) also regarded SMM1 and SMM3 as active massive and intermediate-mass star-forming sites, whereas massive and cold SMM2 showed no sign of ongoing star formation activity.
According to the GHC model, the stellar populations should be distributed as per their age gradients, i.e., the oldest stars are located away from the central part of the cluster, whereas the youngest stars should mainly reside at the central parts of the star-forming site.Our observations revealed a similar scenario.A wealth of Class II YSOs (having typical age of ∼13 Myr; Evans et al. 2009) are primarily distributed in the outer part of the cloud and in the filaments, almost without the presence of younger Class I YSOs.The recent star formation activity, inferred from the presence of the UC H II region (typical age of 0.3 Myr; Motte et al. 2018), Class I YSOs (typical age of 0.44 Myr; Evans et al. 2009), and submillimeter clumps, is noted in the central hub region.Note that identified YSOs have a mass range of 0.1-6.0M e (see Section 3.3.5).Thus, the YSOs distributed all over the regions including the filaments and the hub are primarily composed of low-and intermediatemass stars.The only signature of massive star formation (i.e., SMM1 and UC H II region) is noted in the central hub.These observational signatures are in accordance with the observational fact predicted in the GHC scenario.Thus, it is most likely that the GHC (Vázquez-Semadeni et al. 2017, 2019) model is operational for the formation of HFS and also for the overall star formation activity in the S87 region.

Conclusions
We performed a multiwavelength study of the S87 H II region to explore the star formation activity in the region.The major findings of the study are listed below.
1. From the Herschel column density map, we identified the HFS associated with the region S87.Velocity distribution and velocity dispersion maps of the molecular line data show a substantial velocity gradient along the filaments, and the highest velocity dispersion is noted in the hub, signifying its turbulent nature.2. The extinction map generated from the NIR data shows a nonuniform distribution of gas in the region.The extinction map closely follows the column density map of the region and depicts the presence of an HFS. 3. We identified 13 Class I, 202 Class II, and 22 Hα-emitting sources in this region.The estimated average age of YSOs is 3 Myr, with masses in the range of 0.1-6.0M e .4. Evolved Class II YSOs are mostly located in the outer parts of the host cloud and are also distributed along the filaments.The younger star formation activity, along with the signature of massive star formation, is mainly noticed in the hub region.Overall, our comprehensive multiwavelength analysis revealed that the GHC model of star formation is operating in the S87 region.(This table is available in its entirety in machine-readable form.)

Figure 2 .
Figure 2. (a) Herschel column density map of the region for a 30′ × 30′ area centering on S87.The dashed black box indicates the target area in this study.(b) A zoomed-in view (20′ × 20′) of the column density map for the targeted area.Filaments (marked in red) identified using the FILLFINDER algorithm are overlaid on the map.The central hub is marked with a yellow contour at the level of 10 22 cm −2 .A scale bar referring to a spatial scale of 2 pc (assuming a distance of 2.3 kpc to the source) is also marked in both panels.

Figure 3 .
Figure 3. (a) Integrated intensity, (b) velocity distribution, and (c) velocity dispersion maps of the region.A scale bar referring to a spatial scale of 2 pc is also marked in the lower right corners of all the panels.The identified filaments are marked in black.

Figure 4 .
Figure 4.The extinction map of the region generated using the (H − K ) colors of the point sources and the SCUBA 850 μm map.The color bar shows the variation of A V .The H 2 column density contours at [3.5, 4, 6, 10, 20, and 50] × 10 21 cm −2 are also overlaid on the map.

Figure 5 .
Figure 5. (a) MIR ([3.6]−[5.8]/[4.5]−[8.0])CCD.The red and blue squares indicate Class I and Class II sources, respectively, in all the panels.The black arrow is the reddening vector A K = 5 mag.(b) AllWISE (W1−W2)/(W2−W3) CCD of the point sources.(c) The dereddened [[3.6]-[4.5]]0 /[K −[3.6]] 0 CCD generated using cross-matched IRAC and 2MASS catalogs.The reddening vector corresponding to A K = 2 mag is also shown.(d) NIR (J − H)/(H − K ) CCD of the point sources in the S87 region.The green and red solid lines represent the locus of the unreddened dwarfs and giants adopted fromCohen et al. (1981), and the black solid line is the locus of the CTTSs(Meyer et al. 1997).Three parallel reddening vectors are shown with three black dotted lines, drawn from the base of the dwarf, the tip of the giant branches, and the top of the CTTS loci, respectively.The black arrow represents the reddening vector A K = 5 mag.

Figure 6 .
Figure6.The (r − i/r − Hα) CCD of the sources detected in the IPHAS catalog.The solid black and solid green nearly vertical lines correspond to the trend for an unreddened Rayleigh-Jeans continuum and an unreddened optically thick disk continuum(Barentsen et al. 2011), respectively.The solid and dashed blue lines imply the location of the unreddened MS and the most probable location of the unreddened MS stars that have Hα emission-line strengths of EW = −10 Å.The dotted black lines represent the predicted lines of constant emission EW.The green diamonds represent the selected Hα emission sources, whereas the gray dots are the sources having no significant Hα emission.

Figure 7 .
Figure 7. (a) Optical CMD for the age distribution of YSOs and Hα-emitting sources.The PMS isochrones for ages 0.1, 0.5, 1, 3, 5, 10, and 20 Myr, as well as the ZAMS adapted from Marigo et al. (2017), are overplotted.In addition to that, different evolutionary tracks for different masses are also plotted.The object symbols are the same as the previous IR CCD and Figure 6.(b) The NIR CMD displays the mass spectrum of the detected YSOs and Hα emission sources.The black, blue, red, and magenta solid lines are the locus of the ZAMS (Girardi et al. 2002) and PMS isochrones (Marigo et al. 2017) of ages 1, 3, and 20 Myr, respectively.The reddening vectors for different mass values are shown by the parallel black dashed lines, and other symbols for the objects are the same as in Figure 7(a).

Figure 8 .
Figure 8.(a) The spatial distribution of YSOs on the integrated intensity map.The Class I, Class II, and Hα emission sources are marked by red circles, green diamonds, and blue triangles, respectively.The yellow contour corresponding to H 2 column density of 10 22 cm −2 is also overlaid on the map to depict the central hub region.(b) The zoomed-in view (4 2 × 4 2) of the central hub region overlaid with clumps SMM1, SMM2, and SMM3 identified byXue & Wu (2008).The SCUBA contours (green) at 0.35, 1, 2, 5, 7, and 10 mJy arcsec −2 are also overlaid in this image, while the background rms noise level is at 0.07 mJy arcsec −2 .The presence of ionized gas is also marked in white contours of CORNISH 5 GHz data, drawn at the [5σ, 10σ, 20σ] levels (where σ = 0.26 mJy beam −1 ).
Note.The entire catalog is available in the electronic version.(Thistable is available in its entirety in machine-readable form.)11

Table 2
Photometric Catalog of Hα Emitters toward S87 IINote.The entire catalog is available in the electronic version.