# Fragmented atomic shell around S187 H ii region and its interaction with
molecular and ionized gas

Petr Zemlyanukha1, Igor I. Zinchenko1, Evgeny Dombek1,2, Lev E. Pirogov1,
Anastasiia Topchieva3, Gilles Joncas4, Lokesh K. Dewangan5, Devendra K. Ojha6
and Swarna K. Ghosh 6
1Institute of Applied Physics of the Russian Academy of Sciences, 46 Ul’yanov
str., Nizhny Novgorod 603950, Russia
2Lobachevsky State University of Nizhny Novgorod, 23 Gagarin Ave, Nizhny
Novgorod 603950, Russia
3Institute of Astronomy, Russian Academy of Sciences, Pyatnitskaya str., 48 ,
Moscow, Russia
4Centre de Recherche en Astrophysique du Québec, Département de Physique, de
génie physique et d’optique,
Université Laval, Québec, QC, G1V OA6, Canada
5Physical Research Laboratory, Navrangpura, Ahmedabad 380009, India
6Department of Astronomy and Astrophysics, Tata Institute of Fundamental
Research, Homi Bhabha Road, Mumbai 400 005, India

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The environment of S187, a nearby H ii region (1.4$\pm$0.3 kpc), is analyzed.
A surrounding shell has been studied in the H i line, molecular lines, and
also in infrared and radio continua. We report the first evidence of a clumpy
HI environment in its photodissociation region. A background radio galaxy
enables the estimation of the properties of cold atomic gas. The estimated
atomic mass fraction of the shell is $\sim$260 M⊙, the median spin temperature
is $\sim$50 K, the shell size is $\sim$4 pc with typical wall width around 0.2
pc. The atomic shell consists of $\sim$100 fragments. The fragment sizes
correlate with mass with a power-law index of 2.39-2.50. The S187 shell has a
complex kinematical structure, including the expanding quasi spherical layer,
molecular envelope, an atomic sub-bubble inside the shell and two dense cores
(S187 SE and S187 NE) at different stages of evolution. The atomic sub-bubble
inside the shell is young, contains a Class II young stellar object and OH
maser in the centre and the associated YSOs in the walls of the bubble. S187
SE and S187 NE have similar masses ($\sim$1200 M⊙ and $\sim$900 M⊙,
respectively). S187 SE is embedded into the atomic shell and has a number of
associated objects including high mass protostars, outflows, maser sources and
other indicators of ongoing star formation. No YSOs inside S187 NE were
detected, but indications of compression and heating by the H ii region exist.

###### keywords:

stars: formation – ISM: clouds – radio lines:ISM - ISM: HII regions - ISM:
evolution - ISM: bubbles

††pubyear: 2021††pagerange: Fragmented atomic shell around S187 H ii region
and its interaction with molecular and ionized gas–24

##  1 Introduction

The expansion of an H ii region is thought to be a factor that induces star
formation in the interstellar medium especially as a trigger of high mass star
formation. Such a trigger can arise from the overpressurization of gas that
bounds the region (Thompson et al., 2004) leading to compression and
fragmentation. The dynamics of this process is usually studied with Ci/Cii/Oi
observations (Requena-Torres et al., 2016; Okada et al., 2019; Bally et al.,
2018). The atomic hydrogen production in a photodissociation region (PDR) is
quite natural and is predicted by the existing models (Kirsanova & Wiebe,
2019; Kirsanova et al., 2020; Sternberg et al., 2014). Still, there are few
observations of atomic shells around H ii regions in the neutral hydrogen line
(Roger et al., 2004; Cichowolski et al., 2009; Joncas et al., 1992). Such
rarity is caused by the challenges of H i observations. The H i emitting zone
near the ionized region is usually garbled by the Galactic atomic gas in the
line of sight. An estimation of the properties of the atomic gas is difficult.
The H i spin temperature is correlated with the optical depth leading to
uncertainties in estimations (Saha et al., 2018). Absorption profiles are
essential for these estimations. Sometimes one is helped when a strong
background continuum source occurs behind the studied object (Nguyen et al.,
2019).

To study the atomic hydrogen gas around a H ii region, high resolution H i
data should be obtained and compared with kinematics, morphology and
properties of other species.

This paper presents such analysis for the S187 (Sh2-187, KR120, LBN
126.7-00.80) environment. This region is well-studied from the optical to the
radio wavelengths. The S187 H ii region (Fig. 1) is driven by a cluster of B0V
stars with a photon flux of 1047 ionising photons s-1 (Joncas et al., 1992)
historically named as S187 B or BD52. The region is fairly young with an
estimated age of $\simeq$5$\times$105 yr (Joncas et al., 1992). The distance
estimate is 1.4$\pm$0.26 kpc from photometry (Russeil et al., 2007) and
$\sim$0.9 kpc from Gaia DR2 data and CO line analysis of associated large-
scale structures (Yan et al., 2021). Since the uncertainty of S187’s GAIA-
based distance estimation is still not given we adopt the 1.4$\pm$0.26 kpc
distance. The H ii region is surrounded by a molecular and atomic gas shell
(Joncas et al., 1992; Arvidsson & Kerton, 2011). A number of young stellar
objects were detected towards the shell (Kang et al., 2017; Richards et al.,
2012). Three sources with maser emission were also detected including two OH
(Engels & Bunzel, 2015) masers and an H2O maser (Valdettaro et al., 2001)
associated with a CH3OH maser (Slysh et al., 1999), indicating ongoing star
formation in the shell. Two compact radio sources were reported as known
radio-galaxies (Snell & Bally, 1986; Richards et al., 2012). They were
misinterpreted as UC H ii regions (Israel, 1977) and have historical names
S187-1(a&b) (NVSS 012258+614815, ABG85 94) and S187-2(a&b) (NVSS
012250+615451).

Figure 1: The structure of the S187 complex. The 1420 MHz GMRT continuum is in
red contours, tracing the ionized gas. Bright red spots are related to radio
galaxies (S187-1(a&b) and S187-2(a&b)). The red contours near S187-1(a&b) are
caused by limited dynamical range. The noise on the top - right corner is
caused by primary beam correction The 12$\mu$m WSSA image is shown in
grayscale tracing the dust in shocked gas. White contours represent SCUBA
850$\mu m$ emission tracing cold dust in molecular material. The water maser
(Sunada et al., 2007) is marked as a triangle and the OH masers from Engels &
Bunzel (2015) are marked as reversed triangles. S187 NIRS 1 is shown as a
circle and S187 H$\alpha$ as the star.

The formation of IRAS 01202+6133 ($M\sim$21 M⊙, $L\sim$5600 L⊙) is argued to
have been triggered by the expansion of the S187 H ii region (Kang & Kerton,
2012). A pre-main-sequence object S187H$\alpha$ was discovered by Zavagno et
al. (1994) with a mass of 25$\pm$12 M⊙ nearby to IRAS 01202+6133.. A 12CO
outflow is associated with S187H$\alpha$ (Zavagno et al., 1994). A molecular
hydrogen outflow driven by T Tauri star S187 NIRS 1 (S187:SCP 2) reported in
Salas et al. (1998). The H2O (Henkel et al., 1986; Han et al., 1995; Palagi et
al., 1993; Valdettaro et al., 2001) and class II CH3OH maser (Slysh et al.,
1999) sources are located nearby.The existence of such objects is evidence of
ongoing star formation in the shell. An OH maser (Engels & Bunzel, 2015) was
also detected to the north of IRAS 01202+6133, but we are not aware of any
special study of the associated YSO. The dense molecular gas properties were
discussed by Zinchenko et al. (2009); Kang & Kerton (2012) (see also
references therein). The molecular part of the shell was observed in 12CO,
C18O and 13CO (1–0) and discussed in Arvidsson & Kerton (2011), with the total
mass estimate of about 7600 M⊙.

A multifrequency study of S187 and its environment was performed by Joncas et
al. (1992), including H i data with 1′ angular resolution and 1.32 km s-1
spectral resolution at 4.9 K sensitivity. According to Joncas et al. (1992),
the H ii region is interacting with the surrounding neutral gas. The ionised
region is surrounded by the shell containing the H i counterpart. The atomic
gas distribution is inhomogeneous, has a systematic velocity gradient and the
total mass is 70 M⊙. The associated YSO population is reported by Kang et al.
(2017); Snell & Bally (1986). Two types of dust populations with grains larger
than 0.004 $\mu$m were detected in infrared emission (Joncas et al., 1992). In
summary, the S187 environment is a good target to study the interaction of the
ionized, atomic and molecular gas and its impact on star formation.

In this study, we present new GMRT HI and continuum observational data of the
S187 region including two molecular cores (S187 NE & S187 SE) in the shell
around the H ii region that illustrate different phases of the interaction
between molecular, atomic, and ionized gases. Section 2 presents the
observational data at 21 cm, radio continuum, molecular lines HCO+(1–0),
CS(2–1), 12CO(1–0), 13CO(1–0), C18O(1–0), and infrared data from WISE. Section
3 describes the observational results derived from these datasets. More
details and the analysis of the observations are presented in Section 4. Also,
we discuss our results in Section 5, and at the end we give the conclusions.

##  2 Observations and archival data

###  2.1 21 cm line and radio continuum

The 21-cm data were obtained using the Giant Metrewave Radio Telescope (GMRT,
project 26_012) during 2014 September 4-5. The legacy GMRT backend was
configured to the 8.138 kHz (1.72 km s-1) resolution and the 4166.7 kHz
bandwidth. 3C48, 0102+584, 0217+738, and 3C147 were observed as calibration
sources. 3C48 was used as a bandpass and fluxscale calibrator, 3C48, 0102+548,
and 0217+738 as phase calibrators. Absorption features were detected towards
3C147, 0217+738, 0102+584. The channels affected by absorption were excluded
during the calibration. Around 30% of the data were manually flagged mostly
due to radio frequency interference in some baselines. The calibration was
performed in CASA (McMullin et al., 2007). A strong compact radio source NVSS
012258+614815 in the source field was used for phase self-calibration. The
continuum image RMS noise was 0.1 mJy/beam at the field centre for Briggs-
weighted image restoring (beam: 2.28$\times$2.1 arcsec, PA -2∘). A UV-tapered
continuum image was also produced to achieve the spatial sensitivity with a
beam of 8.3$\times$7.0 arcsec, PA -10∘ and 0.7 K RMS.

Two spectral cubes were processed with different restoring beams. One with the
Briggs weighted beam (2.7$\times$2.2 arcsec, PA -13∘) and the other with the
uv-tapered naturally weighted (8.11$\times$7.28 arcsec, PA 7.4∘) beam to
achieve sensitivity to the extended structures. We performed uv continuum
subtraction for the uv-tapered data and combined it with the WSRT CGPS H i
line survey (Taylor et al., 2003) using CASA sdintimaging (Rau et al., 2019)
to fill the flux loss. The standard clean procedure was applied with cutoff
limit of 10 mJy/beam. The tapered combined image has a per-channel RMS noise
around 3 K, the Briggs weighted GMRT-only datacube has a sensitivity of 1
mJy/beam (80 K) RMS noise per channel.

###  2.2 Molecular lines

New HCO+(1–0) observations were carried out in 2021 April with the 20-m
telescope of the Onsala Space Observatory (OSO) (project O2020b-02). The 3 mm
dual-polarization receiver with a SIS mixer was used (Belitsky et al., 2015).
A Fast Fourier Transform spectrum analyser with a 76 kHz frequency resolution
was used. The spectrometer gives a velocity resolution of $\sim 0.26$ km s-1
at 89.2 GHz. To obtain high-sensitivity data, the signals from two
polarizations were co-added during the reduction process, which provided an
increase of sensitivity by a factor of $\sqrt{2}$. The observations were
performed in the frequency switching mode. The system noise temperature varied
in the range of $\sim 170-290$ K depending on weather conditions. The main
beam efficiency was $\sim 0.44$ for elevations of $\sim 30^{\circ}$ at which
the source was observed. The full beam width at half maximum (FWHM) of the
OSO-20m telescope at 86 GHz is 43 arcsec. Pointing and focus were checked
regularly by observations of SiO masers. Mapping was done with a 20 arcsec
grid spacing. The line data were processed using the CLASS program from the
GILDAS  package111http://www.iram.fr/IRAMFR/GILDAS, the XS
software222https://www.chalmers.se/en/researchinfrastructure/oso/radio-
astronomy/Pages/software.aspx, and our original routines.

The CS(2–1) line was observed by Joncas. It was obtained on 1987 November 30
using the 14 m telescope of the Five College Radio Astronomy Observatory
(FCRAO) observatory. The data were not previously published. The telescope was
equipped with a cooled 3 mm receiver connected to a digital spectrum expander
as the backend with resolution of 50 kHz and 256 channels. The weather was
clear, and typical system temperature around 400 K. The individual scans were
processed using the IRAF software (Tody, 1986) and mapped using our scripts
based on astropy package (Astropy Collaboration et al., 2013). The spectral
resolution is 0.25 km s-1 and the angular resolution is 45 arcsec. The beam
efficiency was taken from FCRAO GRS as 0.48 (Jackson et al., 2006). We assume
the same flux uncertainty of 15% as in this paper. This data has a larger
coverage (1h23m51s.2 +61∘51$\arcmin$00″9 centre, 10 arcmin square field) than
the CS (2–1) observations reported in Zinchenko et al. (2009).

Observations of the 12CO (1–0), 13CO (1–0) and C18O (1–0) lines were kindly
provided by Kim Ardvisson and Cris Ardvisson and were done with the 14 m FCRAO
telescope, processed by C. Brunt with the spectral resolution of 0.25 km s-1
and angular resolution of 46 arcsec (Arvidsson & Kerton, 2011). We also used
12CO (1–0) data from FCRAO CGPS to study the extended structures around the
S187 region. The HCO+ (1–0) data towards IRAS 01202+6133 are from Zinchenko et
al. (2009). They were obtained at the Onsala observatory. The spectral
resolution is 84 m s-1 with a 0.4 K RMS noise.

###  2.3 IR data

We used the WISE archive data images to study the dust emission in the region.
The 2MASS (H and Ks) (Cutri et al., 2003) and WISE (i.e., W1, W2, and W3)
photometric data we utilized to classify YSOs in the area and Spitzer IRSA,
unWISE (Meisner et al., 2017) and WSSA (Meisner & Finkbeiner, 2014) archive
data for presentation purposes. We also used the SCUBA data archive for a 850
$\mu$m image that was previously described in Arvidsson & Kerton (2011).

##  3 Observational results

We present the 1420 MHz continuum and the 21-cm line data in absorption and
emission. We discuss the compact features first and then concentrate on the
extended structures. We also present molecular gas emission morphology in
addition to those published in Arvidsson & Kerton (2011).

###  3.1 Radio continuum

The structure of the 1420 MHz continuum emission is consistent with previous
studies (Snell & Bally, 1986; Arvidsson & Kerton, 2011). Six compact sources
are detected in the field, five of them are near the S187 H ii region (Fig.
1), four of them are identified in the NVSS catalog. The fluxes are given in
Table 1.

The parameters were estimated using CASA 2D Gaussian feature fit. The
brightest one, NVSS 012258+614815 (S187-1) is located near the S187 H ii
region. It is known as a radio galaxy (Snell & Bally, 1986; Arvidsson &
Kerton, 2011). We resolve it in two separate compact sources and an elongated
extended weak structure. The compact sources are separated by 10 arcsec (Fig.
3), with brightness of 24000$\pm$120K (S187 1a) and 14300$\pm$120K (Fig. 3).
NVSS 012250+61545 (S187 2) is also split into two compact features. A weak
point source with a flux of 0.67$\pm$0.14 mJy is detected towards IRAS
01202+6133. This source is associated with the Class I protostar reported in
Kang & Kerton (2012). A water vapor maser is present near IRAS 01202+6133 as
demonstrated by Zinchenko et al. (1998). This source is a region of high-mass
star formation (Kang & Kerton, 2012). To the best of our knowledge, no radio
continuum emission was previously reported towards IRAS 01202+6133. Thus,
there is no data to estimate the spectral index. The continuum point source
can represent an UC H ii region.

Figure 2: Closeup view of NVSS 012258+614815 (S187-1 a&b) in the GMRT 1420 MHz
continuum. The beam is shown in the lower left corner. Fluxes are in
logarithmic scale.

Figure 3: The H i spectra towards NVSS 012258+614815 (a and b). Emission
spectra are extracted from the CGPS survey inside the $1^{\prime}$ circle
around the galaxy position. The absorption profiles are obtained from the GMRT
data. The brightness values below 50 K are shown using a linear scale those
above 50K with a logarithmic scale.

S187-1a&b are brighter than the extended H ii region requesting high quality
calibration. The UV-tapered continuum image IS shown in Fig. 1 in red. The
emission associated with S187 is not homogeneous. The shape is close to a
ball-like structure overlaid with a curved feature. This feature is bordered
by 850 $\mu$m SCUBA emission associated with the parental cloud of IRAS
01202+6133, as seen in Fig. 1. There are two weak peaks in radio continuum
emission, located near the SCUBA 850 $\mu$m and WISE 12$\mu$m peaks (Fig. 1).
We associate these features with PDRs and call them S187 PDR1 and S187 PDR2,
(coord. in Table 2, see Fig. 4). The SCUBA 850 $\mu$m emission is detected
near the IR peaks associated with PDR1 and PDR2. Weak continuum emission
surrounds the S187 PDR2 from north and south (Fig. 1, Fig 4, center; , Fig 4,
right). Such morphology is an argument for a small, density bounded zone
bordered or surrounded by the leaking ionising emission and gas.

Table 1: Continuum sources properties and parameters estimated using the CASA 2D Gaussian feature fit. Source | Coordinates, J2000 | Flux at 1420 MHz, mJy | Size, deconvolved: major,minor,p.a., arcsec
---|---|---|---
S187-1a | 01h22m58s.11304(.0006) +061∘48′15″.188(.0037) | 386.4(2.5) | 1.27(0.02)x0.38(0.063)x91(2.5)
S187-1b | 01h22m59s.39627(.002) +061∘48′13″.52(.013) | 116.7(0.8) | 1.4(0.08)x0.35(0.2)x112(4.9)
S187-2a | 01h22m51s.346(.0017) +061∘54′48″.22(.015) | 13.7(0.3) | 1.3(0.1)x0.58(0.16)x144(7.4)
S187-2b | 01h22m50s.91(.02) +061∘54′51″.93(.14) | 7.9(0.7) | 4.42(0.4)x3.4(0.4)x131(22)
IRAS 01202+6133 | 01h23m33s.36(.012) +061∘48′48″.09(.15) | 0.67(0.14) | point
NVSS 012355+614404 | 01h23m54s.1(.012) +061∘43′57″.73(.08) | 1.8(0.24) | point
S187 | 01h23m04s.8(.1) +061∘51′41″.2(.4) | 1200(20) | 175.62(3.6)x141.3(3)x162(4)
Table 2: Molecular and infrared sources sizes estimated using the CASA 2D Gaussian feature fit. The integrated C18O intensity image is marked as $\int$C18O. Source | Coordinates, J2000 | Size: maj,min,bpa, arcmin | source
---|---|---|---
S187-NE | 01h23m29s.9(1.37) +061∘54′28″.8(8.9) | 4.1(0.4)x3.65(0.2)x129(7.5) |  $\int$C18O
S187-SW | 01h23m30s.5(1.2) +061∘48′49″.7(7.0) | 4.2(0.4)x3.4(0.3)x93(17) |  $\int$C18O
S187-PDR1 | 01h23m17s.93(.03) +061∘53′10″.34(.2) | 3.12(.01)x2.00(.006)x128(0.29) | 12$\mu$m WISE
S187-PDR2 | 01h22m50s.2(.02) +061∘52′28″.8(.1) | 2.04(.008)x1.27(.005)x69.4(.3) | 12$\mu$m WISE
HI subbuble | 01h23m39s.8(.1) +061∘51′37″(1.5) | 3.9(.1)x4.1(.1)x90 | 12$\mu$m WISE

###  3.2 The H i data

Using the combined GMRT+CGPS data we are able to obtain full-scale high-
resolution datacubes of the atomic hydrogen lines. From the GMRT data only we
obtain high-contrast absorption profiles towards S187-1(a&b). The 21-cm
emission and absorption line profiles are shown in Fig. 3. These profiles
contain at least 8 spectral features. The emission features at [–5;–20] km s-1
are located near the S187 H ii region, at the same velocities as molecular gas
traced by 12CO (discussed later). These velocities argue for an Orion arm
membership (Joncas et al., 1992).

Figure 4:  a) The observed features in the S187 complex. Background image is
the Spitzer 5.8$\mu$m image overlaying the WISE WSSA 12$\mu$m image (in green
and yellow shades). The C18O average intensity in the [-17.5,-12] km/s range
is shown in blue filled contours while black contours represent the 21-cm line
emission in the same velocity range. The contour width is increasing towards
the brighter emission. The beams (46 arcsec, 8 arcsec, and 8 arcsec for the
C18O, H i and 1420 MHz continuum, respectively) are shown in the lower left
corner of the image. the position of S187 H$\alpha$ is marked as a red star. A
white cross marks the HCO+ line intensity peak position. The other marks are
the same as in Fig 1. b) Close-up picture of S187 PDR2. The image scheme is
the same as for a), except the velocity range is [-13,-10] km/s. Light gray
dotted contours represent 1420 MHz continuum emission at a 7K intensity level.
Declination tick marks are placed at +061∘51′24″and +061∘51′36″. c) Close-up
picture of S187 PDR1. The image scheme is the same as for b), except the line
which is HCO+ the velocity range is [-17.-13] km/s since the S187 PDR1 central
velocity is different.

We resolved NVSS 012258+614815 as two compact sources giving the possibility
for highly accurate estimation of the optical depth and column density of the
absorbing gas (see next section). Other compact sources are too weak to detect
absorption between continuum and hydrogen line brightness levels ($\sim$100 K)
in our beam.

The 21-cm line emission from the combined CGPS+GMRT data are shown as channel
maps in Appendix B (Fig. 22). The maps are overplotted with black contours
which show the 12CO (1-0) emission from Arvidsson & Kerton (2011). The
observed structures are in good consistency with those discussed by Joncas et
al. (1992) despite an order of magnitude improvement in resolution (from 60
arcsec to 8 arcsec) and a factor of 3 improvement in sensitivity (from 0.7 to
0.3 mJy/beam) in restored beam. The H i emission towards S187 is not
homogeneous, the channel maps reveal spatially and kinematically connected
structures up to scales of $\sim$0.6 pc. The features that are spatially
related to the H ii region are detected from –20.7 to –3.5 km s-1. At the
edges of the velocity range, the emission is located towards the centre of the
H ii region tracing the parts of the shell walls oriented with the motion
parallel to the line of sight. As shown in Fig. 4, the integrated H i line
emission correlates well with the 5.8 $\mu$m Spitzer emission. The local peaks
of the H i line are slightly shifted from the infrared data peaks outwards of
centre of the H ii region. The emission region has a shape close to a circle
with inhomogeneous borders except for the north-east direction, where the
border is curved inside.

###  3.3 Molecular lines

The shell around S187 has a molecular counterpart (Joncas et al., 1992;
Zinchenko et al., 2009; Arvidsson & Kerton, 2011). It consists of dense
regions with molecular cores and more diffuse, surrounding gas generally
related to the ionized region envelope. We used the CS(2–1) and HCO+(1–0)
observations as tracers of the dense regions. The molecular gas emission in
the CO isotopologue lines are presented AS channel maps (Figures in appendix
B, 23, 22, 24). The 12CO emission lines are detected all over the field. The
13CO emission region is more extended than the H i associated with S187. The
C18O emission is usually optically thin, tracing the column density of the
molecular gas. We associate it with the dense cores in the molecular shell at
different stages of evolution (discussed below). There are two peaks detected
at different velocities, hereafter called S187 SE, the parental cloud
associated with IRAS 01202+6133 and S187 NE. S187 NE spatially and
kinematically correlates with the 12CO peak while S187 SE does not. At the
velocity of S187 SE ($\sim$ –15 km s-1) the 12CO line is mostly homogeneous
(Fig. 22). The bright peaks of the C18O emission strongly anti-correlate with
the H i data, most strongly towards S187 SE. The molecular emission elongates
towards the S187 exciting star cluster, a small peak of the H i emission is
detected towards the centre of the H ii region. No 12$\mu$m WISE feature is
detected towards or near S187 SE. This absence can be a projection effect on
the line of sight meaning that the molecular gas is not related to the S187
shell, but the anti-correlation between H i and C18O suggests that the
molecular gas of S187 SE is embedded into the atomic medium. No H i emission
is detected at the position and velocity of S187 SE C18O and HCO+ peaks. No
cometary 12CO molecular emission is detected around the C18O peak. The peak of
H i gas, at the -11 km s-1 is located towards the molecular emission peak near
the S187-H$\alpha$ object.

The CS line properties (Fig. 5) are consistent with those in Zinchenko et al.
(2009). In addition we detect an extended feature with a line brightness of
$\sim$0.6$\pm$0.2 K and central velocity of –15 km s-1. The feature extends to
the north of IRAS 01202+6133. The CS and HCO+ peaks are shifted from the C18O
peak outwards of S187B. The HCO+ observations towards S187 NE have similar
properties. Two features also exist. The line peak is centred at 1h23m34s.7
+61∘54$\arcmin$38$\arcsec$.8 with maximum intensity of 2.8$\pm$0.2 K and
velocity of –16$\pm$0.3 km s-1 and a $\sim$2 km s-1 line FWHM. (Fig. 4, a))
The extended feature has a $\sim$3 km s-1 FWHM centered at –15$\pm$0.3 km s-1
and peak intensity of 1.7$\pm$0.25 K. The extended feature is located towards
the S187 PDR1. The profile has non-gaussian features. The line profiles are
shown in Fig. 5.

Figure 5: The C18O(1–0), 12CO(1–0), CS(2–1) and HCO+(1–0) spectra towards some
positions of the S187 shell. The 12CO(1–0) intensity is divided by a factor of
7. Vertical lines mark the C18O peak intensity.

##  4 Physical properties estimation

###  4.1 The properties of the H i gas

The presence of two closely-located strong background sources (S187-1a&b)
gives the possibility for accurate measurement of the absorption profile of
the cold neutral hydrogen gas. We used GMRT-only data (beam: 2.7$\times$2.2
arcsec, PA -13∘) for optical depth profile extraction. Using the emission line
from the CGPS survey we are able to estimate the spin temperature of the
neutral hydrogen. We assume that the properties of the background galactic
atomic gas does not vary significantly across the GMRT field.

We estimate the optical depth towards S187-1(a&b) from the two H i absorption
spectra:

$\tau(v)=-\ln\left[\frac{T_{br}^{a}(v)-T_{br}^{b}(v)}{T_{c}^{a}-T_{c}^{b}}\right],$
(1)

where $T^{a}_{br}(v)$ and $T^{b}_{br}(v)$ are the observed line temperatures
towards the S187-1 components $a$ and $b$ shown in Fig. 3. $T^{a}_{c}$ and
$T^{b}_{c}$ are the continuum brightness temperatures of the components
derived from the same spectra. We suppose the same filling factor for both
components separated by 9.5 arcsec (0.05 pc) and the same optical depth for
the absorbing gas towards both sources. The components $a$ and $b$ have nearly
the same size (Table. 1). The sizes estimated as 1$\times$3 ″and 1$\times$2
″(5 & 15 GHz) for $a$ and $b$, respectively (Snell & Bally, 1986). With this
estimation we also exclude calibration errors between different instruments
and missing flux related to the absorption-emission optical depth estimation
(Nguyen et al., 2019), where absorption spectra are acquired by an
interferometer and emission spectra by single dish. The spin temperature is
calculated using emission spectra averaged near S187-1(a&b) from the CGPS
survey observations over a 1′ region excluding the S187-1 area:

$T_{s}=T_{em}/(1-e^{-\tau(v)}),$

where $T_{em}$ is the line temperature from the CGPS survey (see Fig. 3) and
$\tau$ is taken from Eq. 1. The resulting profiles are shown in Fig. 6. The
high uncertainty of the optical depth estimation at -17.5 km s-1 may be caused
by two factors. First, the presence of extended emission at that channel
increases the RMS and thus the uncertainty of the estimation. Second, S187-1 a
and S187-1 b have similar intensities due to the high optical depth which
leads to A smaller value of the numerator in the expression 1 and thus higher
uncertainties. It is worth noting that $T_{s}$ and $\tau$ estimations are the
line of sight ensemble-averaged values. We do not take into account a possible
beam dilution that may increase the spin temperature value.

Figure 6: Optical depth and spin temperature profiles towards S187-1(a&b). The
flat-layer model fit is shown as a black solid line and the individual layers
emission profiles are shown as dashed lines. The velocity of each layer is
shown by a vertical mark.

To estimate the properties of the atomic hydrogen layer in the entire S187
shell we fit a flat-layers model to the emission spectra from the combined
data. The brightness temperature $T_{b}$ emitted by layer $i$ is

$T^{i}_{b}(v)=T^{i}_{s}(1-e^{-\tau_{i}(v)})+T^{i-1}_{b}(v)e^{-\tau_{i}(v)},$
(2)

where

$\tau_{i}(v)=\tau_{i}\times exp(-(v-\Delta v_{i})^{2}/(2\sigma_{i})),$ (3)

is the optical depth profile of layer $i$, $\tau_{i}$ is the peak optical
depth of layer $i$, $\Delta v_{i}$ is the Doppler velocity and $\sigma_{i}$ is
the Doppler width of the optical depth profile of the same layer. The total
number of layers is defined by the number of components in the Gaussian
decomposition of the optical depth profile derived from the absorption spectra
and was equal to 4. The velocity of the layer decreases with $i$. We assume
that the outer layers represent outer parts of the Galaxy.

We analysed the combined CGPS+GMRT data using the algorithm described in
Appendix A which combines iterative least squares fitting whose initial
conditions are estimated by k-Nearest Neighbors (kNN). Based on the analysis
of the spectral profiles in different positions of the map we constrained some
parameters of the different layers. The component at the velocity of $\sim$ 3
km s-1 does not vary significantly and their parameters are fixed. We also
fixed the velocity, the Doppler width and the spin temperature of the
component at $\sim$ -5 km s-1. These components are probably related to other
galactic structures. Summarizing, nine parameters are varied,
($T^{i}_{s}$,$\tau_{i}$,$\Delta v_{i}$,$\sigma_{i}$) for each component and
$\tau_{i}$ for one component.

Using the constraints described above we fitted the model to the different
positions. The results converged in three passes. The optimal number of the
nearest neighbors is estimated as 5.

As described in the next section, the H i gas has the internal kinematics of
an expanding bubble, so the assumption of a Gaussian line width in this object
is far from ideal. We estimate the column density of the atomic layer from the
following equation using the fit results and the kNN regression:

$N_{HI}=1.823\times 10^{18}cm^{-2}\int\frac{\tau
T_{b}(v)}{[1-exp(-\tau(v))]}dv,$ (4)

from Saha et al. (2018). Here $T_{b}$ is the brightness temperature of the
model component and $\tau$ is from Eq. (3). The column density map of the
atomic gas related to S187 is shown in Fig. 7. The atomic gas associated with
S187 is surrounded by ambient gas. We subtract the mean column density of the
ambient gas and calculate the mass of the H i gas related to S187 as
260$\pm$56 M⊙. The ambient gas column density is estimated by the using the
flat-layer model fit in the outer regions where no 13CO and 12 $\mu$m WISE
emission was detected. Using the model fit and assuming the temperature
broadening we estimated the upper limit of the kinetic temperature of the
atomic gas $T_{kin}<1000$ K. In the case of thermalization, the H i kinetic
temperature should be equal to the spin temperature. The median temperature
over the S187 shell temperature is 50 K. The peak of the spin temperature is
located towards S187 PDR1 with a value of 110 K. The FWHM of the PV diagram
size is $\sim$10 arcmin which we adopt as the H i shell size.

Figure 7: The atomic gas column density map towards S187. The white contours
represent continuum 1420 MHz emission. The marks are the same as in Fig 4.

The HI column density is well-correlated with the bright 12 $\mu$m WISE
emission, which is tracing hot dust and polycyclic aromatic hydrocarbons (PAH)
in PDRs, but with significant deviation from a linear regression. Negative
correlation with the 12 $\mu$m WISE emission is present in the
[0,$\sim$0.0015] Jy/sr WISE flux range. The related points are located at the
periphery and probably caused by a fraction of the atomic gas not being
related to S187. The determination coefficient R2=0.26 is relatively low and
caused by non-systematic deviation from the linear regression. This deviation
exceeds the estimated uncertainties and may trace the different excitation and
PAH generation conditions. It is possible that not all atomic gas is produced
by photo-dissociation of the S187 molecular shell.

####  4.1.1 Atomic gas clumps properties

We analysed the H i datacube using the Fellwalker algorithm (Berry, 2015). The
cutoff limit was set to 7 times the RMS and the [-19,-5] km s-1 range is used.
The results are shown in Fig. 8 and in Fig 22. We detect $\sim$100 artefacts
associated with the S187 atomic shell. We assume that the detected artefacts
are the inhomogeneities of the atomic layer or fragments. Most of the
fragments are concentrated at the shell periphery and correlate well with the
4.6 and 12 $\mu$m WISE emission. Groups of small ones are located near the
S187 PDR1, S187 PDR2 and S187 SE borders. 5% of the fragments are smaller than
the double beam size, most of the fragments are resolved.

Figure 8: Spatial distribution of the H i fragments identified by the CUPID
Fellwalker algorithm. The ellipses represent clump sizes (standard deviation
around centroid). The H i peak intensity is plotted in the background. The
marks are the same as in Fig 4.

We estimate masses from the column density ($N_{HI}=1.823\times
10^{18}cm^{-2}\int T_{b}dV)$) and the geometric sizes of the fragments under
the assumption of low optical depth (Wilson et al., 2012). The fraction of the
optical depth concentrated in the fragments compared to into the interfragment
medium is unknown so we consider it as the lower limit mass estimation. The
fragment size along the line of sight is unknown and is assumed to be the mean
value of sizes in the plane of sky. The results are shown in Fig. 9. The mean
mass of the fragment is 1.14 M⊙. The total mass of the fragments estimated as
$\sim$114 M⊙. Alternatively, the upper limit of fragment masses can be
estimated under the assumption of the same optical depth as for the flat-layer
model. This assumption gives nearly the same total mass of the fragments as
for the entire shell ($\sim$240 M⊙). The actual fragment masses should lay
between these values. Typical volume density is around $\sim$5.4$\times$104
cm-3 for optically thin limit. Detected fragment sizes vary from 0.03 pc to
0.23 pc. The size is estimated as standard deviation of individual pixels
calculated at 1.4 kpc distance, weighted by intensity. Our data indicate a
tight (R2=0.94) power-law correlation between mass estimation and fragment
size (Fig. 10). The power-law index for the mass-size relation is $k=2.39\pm
0.06$.

Figure 9:  The properties of the fragments. Upper left: Individual fragment
atomic gas mass histogram. Upper right: The fragment linear sizes histogram.
Lower left: Virial parameters of the fragments. Lower right: the individual
fragment volume density histogram. The filled bins correspond to the optically
thin conditions and the unfilled ones are based on the flat-layer model. The
bin step has a logarithmic base except for the lower left plot.  Figure 10:
The correlations between fragment parameters. Top left: Linear size versus
line width. Top right: Virial parameters versus linear size scatter. Lower
left: Mass versus linear size. Lower right: Volume density versus linear size.
The red pluses represent the values estimated in the optically thin limit and
the black crosses represent the values estimated by using optical depth values
from the flat-layer model fit results. The best fits are shown as red lines
(thin) and dotted black line ($\tau$ from FLM). The same dependencies from
Roman-Duval et al. (2010); Larson (1981) are shown by dash-dotted lines.

###  4.2 The molecular gas properties

We analysed the 13CO and C18O data from Arvidsson & Kerton (2011). The data
has the same spatial and velocity coverage as H i combined data.

The CO isotopologues column density ($N_{\mathrm{CO}}$) is estimated using Eq.
90 from Mangum & Shirley (2015). We assumed the same excitation temperature
for the 12CO, 13CO and C18O lines and estimated it from 12CO data using Eq.
(15.30) from Wilson et al. (2012). No 12CO emission associated with S187 SE is
detected. For the S187 SE core, the excitation temperature is assumed as 9.5 K
using the kinetic temperature estimation from Kang & Kerton (2012). The
13CO/12CO=0.0145 abundance ratio is calculated for the S187 galactocentric
distance of 9.27 kpc taking into account the galactic gradient of this ratio
(Pineda et al., 2013). The CO to H2 abundance ratio (7.7$\times$10-5 for S187)
is calculated using the galactocentric distance (Fontani et al., 2012). We
estimated the shell mass using the 13CO data as 6400 M⊙. Such estimate is
close to the value from Arvidsson & Kerton (2011) derived from the same data.

We estimated virial parameters for the S187 SE and S187 NE cores using
$\alpha_{vir}=1.2\sigma^{2}\sqrt{S_{eff}}/(10^{-3}M)$, where $M$ is the core
mass in $M_{\odot}$, $\sigma$ in km s-1 is calculated using the second moment
of the ${}^{13}CO$ line, $S_{eff}=\int{T_{b}}ds/T_{b}$ is the core area in pc2
calculated from the 13CO $T_{b}$ image, $ds$ is the elementary area at a
distance of 1.4 kpc (Kauffmann et al., 2013). The $\sigma$ is $\sim$1.4 and
$\sim$0.9 km s-1, the mass is $\sim$1200 and $\sim$900 M⊙ for S187 SE and S187
NE, respectively. The $\alpha_{vir}$ is estimated as 2.75 and 1.4 for S187 SE
and S187 NE, respectively. The core masses are close to each other, the
difference between virial parameters is caused by different line widths, that
is probably caused by non-thermal motions in the cores.

###  4.3 Dust properties from the Akari-FIS & WISE data

We used data from the WISE and AKARI telescopes. The analysis is based on the
3.36, 4.6, 12.08, and 22.19 $\mu$m images from the WISE survey (Wright et al.,
2010) and on 65, 90, 140, and 160 $\mu$m from the AKARI survey (Yamamura et
al., 2010). Our aim was to identify different dust fractions with molecular
and atomic spatial-kinematic features in the shell using the colours at
different IR wavelengths. Spizer data are also available towards S187, but it
has poor coverage and quality. There is no Hershel observations for S187. For
the analysis the WISE and AKARI-FIS data are used instead.

Figure 11: Scatter plot of the [3.4] $-$ [4.6] WISE colour vs the C18O zero
moment. Linear regression is estimated by Huber regression from Pedregosa et
al. (2011). Dot colour indicates the 65 $\mu$m to 90 $\mu$m flux ratio value.
Dots outside of the colorbar range and the C18O sensitivity limit are plotted
as semi-transparent and are excluded from fitting.

Prior to the analysis of the IR images of S187 and the calculation of the flux
ratios, we cleaned the archival IR maps of point sources. They were removed
from the images using an automatic cleaning procedure. The search and removal
of point sources were carried out in three subsequent stages depending on the
brightness of the stars. More information can be found in Topchieva et al.
(2017). The data were then convolved with a kernel as in Aniano et al. (2011)
and reduced to the same pixel sizes and flux units.

The maps (Fig. 12) are plotted for different ratios of infrared wavelengths.
At shorter wavelengths, the emission presumably comes from small grains, while
emission in the 70 and 160 $\mu$m bands comes from relatively large grains
(Draine & Li, 2007). The emission at 12.08 and 22.19 $\mu$m may come from both
stochastically heated small and hot large dust grains and PAH (Pavlyuchenkov
et al., 2012) PAH in 12$\mu$m usually prevails (Reach et al., 2017). The 65 to
90 $\mu$m ratio is used as a large-grain temperature tracer (Draine & Li,
2007).

Figure 12: The S187 shell in the IR continuum. The WISE WSSA 12 $\mu$m image
is in red. The H i column density is shown in contours. The marks are the same
as in Fig 4.

Unlike Reach & Heiles (2021), we found that the the 12$\mu$m to 22$\mu$m ratio
varies from point to point from 0.2 to 0.74 (See Fig. 12). The ratio
correlates with the H i column density. The ratio variation can be associated
with PAH/grains differentiation and processing (Reach & Heiles, 2021). The
ratio peaks are slightly shifted from H i column density peaks implying that
PAH and H i are effectively generated on slightly different conditions.

We also analysed the [3.4]-[4.6] WISE colour and found that it is associated
to the C18O integrated intensity. The colour is plotted in Fig. 11. We also
plot the 65 $\mu$m to 90 $\mu$m ratio from AKARI data in colour. The colour
correlate the C18O intensity under the linear regression. The determination
coefficient $R^{2}$=0.7 for the regression line. The 90$\mu$m flux dominating
region is closer to the regression line than the 65$\mu$m dominated region.
The 90 $\mu$m emission is usually related to large dust grains (Draine & Li,
2007) protected by dense gas represented by C18O emission.

###  4.4 Young stellar population

In order to identify infrared-excess sources or young stellar objects (YSOs)
in our selected target area, we used the 2MASS (H and Ks) and the WISE (i.e.,
W1, W2, and W3) photometric data. Koenig & Leisawitz (2014) presented several
schemes to identify YSOs, and two such schemes are utilized in this paper. We
used one scheme based on the H$-$Ks and W1$-$W2 colour conditions, and the
second one using the W1$-$W2 and W2$-$W3 colours. The results are shown in
Fig. 13.

Figure 13: The distribution of the YSOs near S187. The background image
represents the H i averaged emission in the [-17.5,-12] km s-1 range. The blue
contours represent the C18O emission averaged for the same velocity range.

##  5 Discussion

Our data trace almost all aspects of the high mass star formation: the H ii
region with the ionizing cluster, the PDR regions with the photodissociated H
i layer, the population of young stellar objects, and at least two molecular
cores that could be at different stages of evolution. In Section 5.1 we
discussed large-scale features that may be related to the origin of the S187
complex. In Section 5.2 we discussed the general distribution of the gas and
YSOs around and in the H ii region. In Section 5.3 we discussed the inner
spatial-kinematic structure of the shell.

Figure 14: Map of the S187 surroundings. The averaged CGPS H i emission ( [–14
to –17 km s-1] range) is shown in the background. It is overlaid with the 12CO
emission averaged in the –5 to –17 km s-1 interval. The rectangle indicates
the location of the PV diagram shown in Fig. 15. The crosses correspond to 0
and 2.2 degrees position offsets.

Figure 15: The position-velocity diagram of the atomic and molecular gas near
the S187 region. The PV slice is shown on Fig. 15. Black ellipses indicate the
probable SNR locations.

###  5.1 Large-scale features

One of the interesting feature that is probably not related to S187 is the
existence of narrow absorption and emission HI features at $\sim$5 km s-1
(Fig. 3 and 6). They are in principle "forbidden" in this quadrant by the
rotation curve of the Galaxy (Joncas et al., 1992). The existence of this
feature is interesting but its detailed analysis lies beyond the scope of this
paper. Such features are usually explained by local streaming motions near the
Sun caused by the Galactic density wave (Burton, 1971). Another possible
explanation can be large-scale kinematics caused by the expansion of a
supernova remnant (Cichowolski et al., 2014). The picture of the surrounding
medium is shown in Fig. 15 and is taken from the CGPS H i and 12CO (1–0)
surveys (Taylor et al., 2003). A position-velocity diagram IS shown in Fig.
15. Two nearly circular structures traced in H i and 12CO are detected. Such
PV features can be formed by the expansion of energetic bubbles. The feature
at $\sim$5 km s-1 is probably part of such a bubble. The "unusual" velocity
can be explained by the "snow-plow" effect often accompanying such energetic
phenomena.

We propose that these bubbles are related to supernova events. The northern
one is probably caused by SNR G127.1+0.5 (Leahy & Tian, 2006), and the
southern one is probably related to the pulsar B0105+65 and the G126.1-0.8
shell (Cichowolski et al., 2014). We estimate the dynamical age from
$t_{dyn}=\alpha R/V_{exp}$. $\alpha$ = 0.25 for a radiative-driven expansion
and $\alpha$ = 0.6 for a shell created by the action of stellar winds. The
inner radius of the atomic part of the shell is $R$ and $V_{exp}$ is the
expansion velocity (Cichowolski et al., 2014). The northern bubble age is
t$\sim$3.6$\times$105 yr (radiatively driven) or t$\sim$8.8$\times$105 yr
(wind driven expansion) (Cichowolski et al., 2014). The southern bubble age is
t$\sim$3.0$\times$105 yr or t=7.15$\times$105 yr for a radiative or wind
driven expansion respectively, which gives relatively the same age for both
SNR.

Figure 16: The position-velocity diagram of the atomic and molecular gas.
Normalized IR and radio continuum profiles are plotted at the bottom. The
ellipse represents inner border of the radially expanding shell. The PV slice
is shown in Fig. 20 as the red box.

Figure 17: Scatter plot of the H i vs the 12CO emission. The red-blue colour
indicates the channel velocity. The scatter plot is overlayed with A 2D
histogram shown by black contours. The possible correlation within the [–15.1,
–16.9] km s-1 range is shown AS A white line.

Figure 18: The correlation between the 12$\mu$m WISE flux and H i column
density.

By this view S187 region is located at the intersection of the bubbles. the
molecular gas appears only at $\sim-$12.5 km s-1, on one side of the bubbles.

###  5.2 The shell

Using the data discussed in Sect. 3 we traced different parts of the expanding
gas and dust surrounding the H ii region. The overall picture is shown in Fig.
4. The channel maps of the H i emission compared with the CO isotopologue
lines are shown in Figs. 22, 23, 24. Our observations indicate an expansion of
the atomic gas related to the shell. The expanding part has an almost circular
shape. It is centered on the H ii region. The position-velocity diagram is
shown in Fig. 16. The emission peaks can be associated with the shell walls,
which are brighter near molecular emission peaks. The velocity of the
H$109$$\alpha$ line related to the H ii region is $-$14.6 km s-1 (Joncas et
al., 1992) . It is close to the expanding layer centre velocity. The H i
emission related to S187 is detected in the $-$21.5 to $-$8.5 km s-1 range
(Fig. 16). The PDR atomic layer at the centre should expand directly towards
and outwards with respect of the observer. As a result, no H i emission should
be observed at $\sim-$14.6 km s-1 and towards the H ii region. However, the H
i emission is observed towards the zero offset in the PV diagram (see Fig.
16). This emission cannot be explained as line broadening of the Doppler-
shifted wall. The source of the emission can be the atomic gas inside or
outside of the sphere. The first possibility implies that neutral gas
fragments remain within the H ii region, that are not involved in a systematic
expansion. The other possibility implies the presence of neutral atomic gas
dissociated by leaking radiation from the H ii region that is not involved in
systematic motions either.

Figure 19: The position velocity diagrams across different slices of the shell
(Fig 20). a) - the slice through S187 PDR1 (positive offsets are towards
S187B, the S187 exciting cluster). b) the slice along the direction where no
wall is detected. c) PV slice through S187 PDR2, d) PV slice through S187 SE
(S187-IRS). Normalised plots of the continuum emission and integrated 13CO and
C18O emission along the PV slices are plotted at the bottom of the panel.
(Scale factors: 8 mJy/beam, 7 GJy sr-1, 17 K km s-1, 2 K km s-1, 6 K km s-1,
5$\times$1021 cm-2 for 1420 MHz continuum, 12 $\mu$m WSSA WISE, 13CO, C18O,
HCO+ integrated intensity and H i column density, respectively).

The channel-by-channel scatter plots of the 12CO and H i data does not show
linear correlation (see Fig. 18), except for the [–17,–15] km s-1 range. In
this range a sign of negative correlation is observed and shown as a white
line. It may trace the molecular gas dissociation processes. In the remaining
velocities, no sign of linear correlation is observed, probably resulting from
the limited resolution of the molecular data. The shell’s atomic mass estimate
(260$\pm$56 M⊙) is much higher than estimated in Joncas et al. (1992) (70
M⊙,under the assumption of low optical depth from the absorption profile).
Also, the absorption profiles (Fig. 3) indicate deeper absorption than
reported in Joncas et al. (1992). Such difference arises due to the
significant optical depth correction in our work. Without it both consistent.
In addition, our approach is based on the observations of very compact
sources. The atomic gas is highly fragmented, so the absorbing properties may
vary on small scales. To our luck, the detected fragments have the sizes much
larger (Fig. 9) than the S187-1(a&b) sizes and their separation. The estimate
of the shell’s atomic mass is small compared to the molecular gas mass in the
shell (6400 M⊙, this paper; 7600 M⊙ from Arvidsson & Kerton (2011)) suggesting
that most of the molecular material is intact by the dissociating processes.

The H i shell is thinner towards the dense molecular core, S187 NE. The shell
correlates with the PDR traced at 12 $\mu$m (Fig. 18). Kirsanova & Wiebe
(2019) note that a density-bounded PDR leads to the merging of the H2 and CO
dissociation fronts. In S187 NE the 13CO and H i peaks are separated by less
than 25 arcsec (0.2 pc) (Fig. 24), so the shock fronts are barely separated.
The HCO+ and H i emission regions spatially and kinetically overlap (Fig. 19,
a.) supposing the mixing of the fronts. The 12 $\mu$m peak may be shifted from
the H i due to a projection effect as can be seen from the position-velocity
diagram in Fig. 19, a., on the bottom plots. The 1420 MHz continuum is also
plotted representing the ionized region. (S187 PDR1 is traced by H i and WISE
12 $\mu$m and S187 NE is traced by HCO+ and C18O). The PV cut in the direction
where no PDR and molecular wall is observed is shown in Fig. 19, b. The H i
emission is present, but no compact features exist, the H i layer is broader
than towards S187 PDR1, up to $\sim$2.4 pc. Weak C18O and 13CO emissions are
detected indicating a small amount of molecular material in the shell towards
the location of the PV cut. The PV cut through S187 PDR2 is shown in Fig. 19,
c. The C18O emission is weaker than in the S187 PDR1, but stronger than in cut
b). There are WISE 12 $\mu$m and 1420 MHz continuum peaks suggesting a
compression of the gas.

####  5.2.1 Fragmentation

The shell has systematic motions and related kinematics, as discussed in the
previous section. In addition, we observe smaller features identified as
fragments in Sect. 4.1.1 and estimated their physical properties. The slope of
the mass-size relation (Fig. 10) is the same as found for molecular clouds
(e.g. Elmegreen & Falgarone, 1996; Roman-Duval et al., 2010). Kritsuk et al.
(2013) argued that supersonic turbulence fed by large-scale kinetic energy
injection is at the origin of the law. The H i fragment masses are lower than
masses of molecular clumps estimated by Roman-Duval et al. (2010) for the same
size. This observation may result from the diffusion of the atomic gas during
dissociation or transformation to the ionized gas state. The fragments may
also have significant amounts of molecular material that are not completely
dissociated into atomic gas. This sample can be incomplete on smaller scales
where fragments can not be resolved. Such incompleteness will not affect the
observed relations significantly since the detected fragment sizes vary over
an order of magnitude. The peak of the fragment mass function at $\sim 1-3$ M⊙
(Fig. 9) is very close to the peak of the prestellar core mass function in
molecular clouds (e.g. Könyves et al., 2020). Most probably these fragments in
the H i shell represent remnants of the molecular clumps where molecular
hydrogen was dissociated by the UV radiation of the central star cluster.

The velocity dispersion in the fragments increases with the fragment size
(Fig. 10) resembling one of the well-known “Larson’s laws” (Larson, 1981).
However the slope of our relation is steeper and the absolute values of the
velocity dispersion are much higher than found by Larson for molecular clouds
and condensations of the same size. They are also significantly higher than in
massive protostellar cores (e.g. Zinchenko, 2000). The mass-size relation for
the optically thick limit are closer to those presented in Roman-Duval et al.
(2010), especially for larger fragments. Such differences are in agreement
with these fragments being the eroded material dissociating into the
surroundings with higher turbulent or thermal motions than the ordinary
molecular clumps. The erosion may be caused by photoevaporation (Decataldo et
al., 2019). The high velocity dispersion leads to high values of the virial
parameter for most fragments $\alpha\ga 2$ (Figs. 9, 10), which implies that
they are gravitationally unbound (e.g. Kauffmann et al., 2013). There is a
trend of increasing virial parameter with fragment size. The most compact
fragments are relatively dense and may represent gravitationally bound
objects.

Alternatively, such fragments may form in PDR itself by dynamical compression.
The shock can be caused by the external high pressure from regions heated by
far ultra violet radiation. The fragments then dissociate and/or evaporate and
loose mass. The timescales should be close to the estimated kinematical age of
the sub-bubble (Gorti & Hollenbach, 2002). During the expansion of the H ii
zone, the new material should be involved in the process so the variety of the
physical parameters, especially densities and virial parameters should be
observed. The S187 region is fairly young ($\sim 5\times 10^{5}$yr), but no
radial dependence of the fragment properties is observed. It may be caused by
the observational limits or projection effect. Still, the dependence of
density and virial parameters on the fragment size demonstrate a larger
scatter than the mass dependence on size (Fig. 10).

###  5.3 The stages of star formation

####  5.3.1 The sub-bubble

One of the most complex spatial-kinematical feature in the S187 shell is the
bubble-like H i structure shown in Fig. 4 (labelled H i bubble, the
coordinates are given in Table. 2).H i channel maps and PV diagramm reveal the
presence of the bubble-like structure that is not detected in other tracers.
The spectra are shown in Fig. 5. Weak HCO+ and CS lines are detected. The HCO+
line profile resembles the 12CO line profile towards the sub-bubble. The C18O
and CS peaks are located at –13 km s-1 The position-velocity diagram is shown
in Fig. 20, e). The centre of the position-velocity diagram corresponds to the
Class II source position, and the slice direction is shown in Fig. 20, left.
The 21-cm line emission shows the O-like PV feature tracing the expanding H i
layer around the Class II object with the velocity range from –19 km s-1 to –8
km s-1. The peak of the 12 $\mu$m emission is close to the centre of the ‘O’
“hole” in the H i data. The PV structure of the molecular lines is similar to
the H i, but the velocity dispersion is lower than for the 21-cm line tracing
outer parts of the sub-bubble. No continuum emission is detected towards the
Class II object or the sub-bubble at 1.4 GHz in our data or at 8.2 and 4.9 GHz
(Snell & Bally, 1986). The OH maser is reported in Engels & Bunzel (2015) near
the Class II source. Two molecular line components exist in the 13CO, C18O and
HCO+ emission at –15 and –12.5 km s-1 with the approximate same structure as
the atomic gas. The CS integrated emission has a weak bowed structure in the
direction of the boundary of the sub-bubble, as shown in Fig. 21. The H i
layer has higher intensity towards the H ii region. The C18O and 13CO
position-velocity structure correlates well with the O-shape. The bubble-like
structure is shown in Fig. 20. As shown in Joncas et al. (1992) the stellar
wind from the exciting cluster is pointing towards the bubble which may have
triggered the early, relative to the S187 NE and SE cores, evolution. Based on
the formula from Weaver et al. (1977) we estimate the age of the bubble as
6$\times$104 years.

Figure 20:  (Left:) Locations of the PV slices used in Fig. 19 (white boxes,
with corresponding labeling) and in Fig. 16 (red box). The WISE 3.4 $\mu$m
image is in the background. (Right:) PV slice through cut e, the H i bubble.
The scales, legend and colour bars are the same as in Fig. 19.

This sub-bubble structure reveals an interaction of the Class II object with
ambient gas within the context of the H i shell. The location of the Class II
object near the border of the ionized region suggests A triggered origin for
the exciting star. A number of YSOs are located towards the sub-bubble,
especially at the boundary of S187. Some of them have X-ray emission
counterparts reported in the Chandra point source 2XCO catalog (Evans et al.,
2020). The excess of 12 $\mu$m emission and the 3.4 to 4.6 $\mu$m color
detected towards the sub-bubble (Figs. 1 and 21) imply additional production
of PAHs by the exciting star and the dominance of photo dissociation near the
Class-II source. The Class II source and the associated sub-bubble can lead to
the additional stimulation of the YSO formation in the S187 shell environment.

Figure 21: The molecular gas structure in C18O integrated emission (cyan
contours), CS (2–1) integrated emission in red dashed contours and the tapered
1420 MHz continuum in white (levels: 1, 2 mJy/beam) contours over the
background gray-scale image representing [3.4]–[4.6] colour. Point sources are
not removed for presentation purposes. The marks are the same as in Fig. 4

####  5.3.2 The S187 SE core

The molecular core S187 SE is located to the south-east of the ionizing
cluster S187B (Table 2). S187 SE is associated with the bright IR Class-I
source IRAS 01202+6133 and maser sources. The position-velocity diagram across
the core is plotted in Fig. 19, d. We define the cores as kinematically-
coherent gas detected in C18O. The internal structure of this object is quite
rich. The HCO+, CS, HCN and N2H+ lines were observed by Zinchenko et al.
(2009) and revealed complicated gas structure. The extended structures are
similar to those observed in C18O. Therefore, the core contains at least two
dense gas fragments (probably associated with IRAS 01202+6133 and S187
H$\alpha$). Ammonia emission is also detected towards S187 H$\alpha$ (Zavagno
et al., 1994). There is a peak in the SCUBA 850 $\mu$m image located towards
S187 NIRS1 (See Fig. 1) suggesting the embeddance of this object into the S187
SE core. No associated peak in molecular lines was detected towards S187
NIRS1.

We detected compact continuum emission which possibly traces an UC H ii region
around IRAS 01202+6133. The H i line is less intense towards the centre
position at –15 km s-1. The HCO+ and C18O lines have different profiles (Fig.
5) and spatial distributions (Fig. 19, d): the C18O emission peak is located
towards IRAS 01202+6133. The HCO+ peak at –14.3km s-1 is shifted by $\sim$80
arcsec from the C18O (–13.6 km s-1) peak. The 4.6/3.4 $\mu$m ratio peak is
located towards the HCO+ peak, the CS emission peak is also shifted from IRAS
01202+6133 (Fig. 21). We interpret such differentiation as an interaction
between the UC H ii region and the parent molecular core. The column density
traced by C18O associated with the core is centered at IRAS 01202+6133 where
the dense part of the core is perturbed by the UC H ii region. No self-
absorption feature in HCO+(1–0) (E${}_{u}\approx$4.28K) is detected. The self-
absorption was reported in HCO+(3–2) (E${}_{u}\approx$17K) (Kang & Kerton,
2012). In case the gas is thermalized, this would be in agreement with a
kinetic temperature of 9.5 K for the core.

No 12CO emission peak is detected towards S187 SE. The emission can be
optically thick, unlike S187 NE. The emitting region towards S187 NE is more
extended than the C18O. The 12CO-traced regions, probably, represent the lower
density molecular gas around the core. For S187 SE a different environment is
observed: the C18O, 13CO and HCO+ emission regions are surrounded by H i
emission (see Fig. 19, d). The core is embedded into the shell. Thus, the H i
gas should be produced during photo-dissociation. S187 SE is dense
(1$\times$105 cm-3 Zinchenko et al. (2009)), and bright 12 $\mu$m spots (like
S187 PDR1 and S187 PDR2) should appear at the border between PDR and molecular
gas, but none was detected. The molecular material traced by C18O is elongated
to the centre of the H ii region with the peak near the border of the H i
layer. The orientation of the elongation seems close to the plane of the sky
since no significant gradients are detected in the molecular data position-
velocity diagram (Fig. 19, d). The distribution of the radio continuum
emission indicates pocketing of the molecular material in the H ii region too.
The increasing ratio of the 160 to 12 $\mu$m emission towards the CS and C18O
emission region can be a tracer of the cold dust dominance over the core.

The origin of the [3.4]-[4.6] VS C18O correlation (Fig. 11 and 21) in S187 SE
remains unknown. The presence of recombination lines (Peeters et al., 2004)
can increase the 4.6$\mu$m flux but can only be related to the UC H ii zone.
The H2 lines present in the to 4.6 $\mu$m WISE band also can increase the
colour. H2 emission was observed by Salas et al. (1998) but shows a different
morphology than the colour map. The 3.4 $\mu$m emission can be saturated by a
significant optical depth in the dense core at this wavelength (Pavlyuchenkov
et al., 2012) both in PAH, silicates and grains. The optical depth at 4.6
$\mu$m is lower so such reddening can arise. The [3.4]-[4.6] colour values is
close to those presented in Ybarra et al. (2014) ([3.6]-[4.5] Spitzer) colour
implying the reddening of the PDR emission in the S187 SE core.

The kinematical-morphological structure of the core and a number of active
star formation processes are indicators that S187 SE has suffered from
external pressure. The most obvious (and realistic) explanation is the
pressure caused by H ii region expansion. Still, no shock or ionisation fronts
in PDRs that should be traced in 12 $\mu$m or 3.4 to 4.6 $\mu$m ratio are
detected. The core is elongated towards the center of H ii zone, and the
pressure and ionizing flux can be reduced by the geometry of the core, leading
to perturbation instead of dissociation. Strong dust absorption can also
protect molecular material from dissociation under ionizing cluster radiation.

####  5.3.3 The S187 NE core

The other core (S187 NE) has a mass close to S187 SE (coordinates in Table.
2). No PDR traced by 3.6 or 4.6$\mu$m emission peaks observed under the C18O
contours suggesting that the core is mostly intact from the expansion of the H
ii region. The PDR is located between the H ii zone and the molecular gas, as
shown in Fig. 4 (S187 PDR1). The structure of the PDR, H i and molecular
layers are shown in a PV (Fig. 19, a) plot. The transition between the
unperturbed and the affected by H ii zone material is seen, with the HCO+ line
broadening overlapping with the H i emission. Class-I and Class-II sources are
located between the S187 NE core and S187 H ii region, but no IR sources are
detected in S187 NE according to the WISE images.

The molecular spectra of S187 NE are shown in Fig. 5. Unlike the other cores
and PDRs, the HCO+ and 12CO line profiles are close to the C18O ones, at –16.4
km s-1. No CS peak is observed within A 1 K sensitivity. The CS(2–1) and
HCO+(1–0) transitions upper energies are close (7.1 and 4.28 K, respectively
(Shirley, 2015)). The critical density of the CS(2–1) transition
(1.3$\times$105 cm-3) is two times higher than for HCO+(1–0) (6.8$\times$104
cm-3) (Shirley, 2015). The absence of CS may indicate A low density for this
core or chemical processes varying the abundance.

The PV cut in Fig. 19, a indicates an increase of the FWHM of the C18O line
from 1.4 to 3.6 km s-1 towards the ionized region. This effect can trace the
turbulent motions caused by the expansion of the S187 shell. This effect by
similar effects as for H i fragments (See Sect.5.2.1). The profile is far from
Gaussian, as shown in Fig. 5. The HCO+ line emission (Fig. 19, a) has two
components. One is a compact narrow peak shifted by $\sim$60 arcsec from the
C18O peak outwards of S187 B. The wider ($\sim$5km s-1) weak emission
correlates with the H i emission and is shifted towards S187 B by $\sim$80
arcsec. We associate this component with the photo-dissociation processes,
consistent with modelling (Kirsanova et al., 2020).

The 12CO peak is shifted by 40 arcsec from the C18O peak to the centre of the
H ii region (Fig. 22). Since the 12CO emission can trace the kinetic
temperature (Wilson et al., 2012) this effect may be explained by heating of
the core by the PDR. The PDR is brighter towards the core indicating feedback
from the core to the shell. Another possibility is an optical depth effect,
which should be greater towards the center of the core. The free-free
continuum is also intense in this direction simultaneously with the peaks in
the H i line. The blue spot in Fig. 12 is located just between the cold dust
layer (green) and the PDR (red). This spot is well-correlated with the atomic
gas emission in the 21-cm line and the HCO+ extended feature. This excess
could be related to the shocked gas excited in Br$\alpha$ 4.06 $\mu$m line, as
observed in Orion Bar (Verstraete et al., 2001). This region was observed in
the H2 line (Salas et al., 1998) and no associated peak was reported. We also
note that no extended green object is reported in S187 NE.

Some C18O emission peaks were also detected in the northern and western part
of the shell, as seen in Fig. 23. The coordinates are $\alpha$=1h22m59s.3,
$\delta$=+61∘42′28$\aas@@fstack{\prime\prime}$8 and 1h22m41s.1,
+61∘49′54$\aas@@fstack{\prime\prime}$3\. They are not covered by CS and HCO+
observations and are not discussed in this paper.

###  5.4 The PDRs

Our observations trace the structures of dissociation regions under different
physical conditions. The structures are represented in the position-velocity
diagrams (Fig. 19) and in maps (Fig 4, left). The general structure is close
to that presented in Tielens & Hollenbach (1985); Gorti & Hollenbach (2002).
The close-up map of the S187 PDR1 and S187 PDR2 is shown in Fig 4, b) and c).

We observe free-free continuum, 12$\mu$m, H i CO, HCO+ and C18O emission peaks
sequentially from the S187B exciting cluster to the S187 NE core in S187 PDR1.
The spatial-kinematic structure suggests an edge-on configuration of the PDR
(Fig. 4, c), Fig. 19, a)). The HCO+ line is observed in the molecular core and
towards the wall traced by the H i line, as shown in Fig 19, a, in the center
of the picture. The line width increases towards the S187 PDR1. The line
profile is far from Gaussian and has multiple peaks (Fig. 5). No absorption
profile is observed towards the molecular core S187 NE (unlike for S187 SE,
(Kang & Kerton, 2012)) implying the same towards S187 PDR1. No C18O or SCUBA
850$\mu$m emission is detected towards the HCO+ peak in S187 PDR1. The line
centroid velocity is $\sim$15 km/s.

The HCO+ emitting gas is probably located behind the shocks (5.8 $\mu$m
Spitzer) and the atomic gas shell layer (See Fig. 4, c.). This picture does
not match the classical (Sternberg & Dalgarno, 1995) view of the evolution of
PDR region, where atomic gas is located between H ii and the molecular
material. The structures observed at 5.8 $\mu$m look filamentary, unlike in
S187 PDR2. A Class-II source is detected towards the S187 PDR1 HCO+ line peak.
The gas associated with that peak may survive after the shock passage. The
HCO+ abundance should increase under this process. But, there is no direct
evidence of such scheme and A high-resolution multi-line analysis is necessary
to study the state of the molecular gas in S187-PDR1. The difference between
HCO$+$ and H2 dissociation fronts is modeled in Kirsanova & Wiebe (2019). Such
difference may be increased by fragmentation of the gas. The PV diagram of the
H i line is not symmetric relative TO the molecular line peaks. Molecular
lines are brighter at the positive velocities. The H i line peak is slightly
shifted to negative velocity. Since the atomic gas is generated during
dissociation, the effect may be caused by non-uniform irradiation or weakening
by the HCO+ traced gas fragment.

According to PV diagrams (Fig. 19, c)), the S187 PDR2 region is more or less
edge-on with respect to the observer ($\sim 20^{\circ}$). A weak C18O emission
region associated with SCUBA 850 $\mu$m peak is located at -12 km/s (Fig. 19,
c, Fig. 4, b. and Fig. 1) near S187 PDR2. Structure of the 5.8 $\mu$m Spitzer
and the 4.6 $\mu$m WISE emission suggest shielding by molecular material. The
column density map (Fig. 7) shows that the atomic gas is concentrated near the
H ii zone border. The H i PV diagram indicate presence of the border towards
the H ii zone, same as for (Fig 19, b) where no shielding by molecular
material is observed. The PV diagram (Fig. 19, b) and channel maps (Fig.
22,23,24) trace leaking H ii zone emission trailing up to 4 pc from the center
of the H ii zone. From the spatial-kinematic structure of the H i emission we
conclude that the H ii zone is not effectively shielded by molecular gas or
dust fractions.

##  6 Conclusion

We studied the shell around the S187 H ii region in the H i 21 cm line,
molecular lines, and also in infrared and radio continua. The region is mostly
circular, but the emission is not homogeneous. The 1420 MHz continuum emission
peaks are correlated with the PDR traced by WISE and Spitzer data. We detected
an excess of H i emission towards the H ii region tracing the shell. It has a
shape close to circular, the peaks correlate with PDRs. The position-velocity
diagram in the H i line has the O-type shape typical of an expanding atomic
layer in this case at the boundary of the H ii region. The shell has a
molecular component visible in the 12CO (1–0) line. The C18O (1–0) and H i
emission distribution are spatially and kinematically related. The molecular
part of the shell is more extended than the atomic one. The expansion speed is
$\sim$7.5 km s-1. The sandwich-like composition of the ionized-PDR-HI-
molecular gas structure of the H ii region and its environment is defined by
different tracers.

Summary:

  1. 1.

We present the results of high-resolution 8 arcsec (0.06 pc) observations of
the H i line emission towards the S187 H ii region. To our knowledge, this is
the highest resolution available of the H i emission observations of Galactic
regions.

  2. 2.

We resolved the NVSS 012258+614815 source as two separate sources with
individual H i absorption profiles. By combining with the emission spectra we
estimated optical depth profiles. The mass of the atomic gas related to S187
is estimated as 260$\pm$56 M⊙ using the emission and absorption data. The
median spin temperature over the shell is $\sim$50 K. The size of the shell is
$\sim$4 pc. The H i layer has varying thickness from 0.2 pc to 4 pc.

  3. 3.

The distribution of the neutral atomic gas is close to an expanding sphere,
spatially correlating with the WISE 12 $\mu$m emission and traces the
interaction between the H ii region and the molecular cloud.

  4. 4.

The atomic shell is inhomogeneous and contains $\sim$100 fragments of median
mass around $\sim$1.1 M⊙. The total mass estimation is in 114-240M⊙ (depending
on optical depth assumption) which is close to the overall shell mass. The
fragment sizes vary from 0.03 to 0.23 pc and the index of the mass-size power-
law fit is 2.33-2.6 range. This estimation is close to previously reported
values available in the literature for molecular gas clumps and clouds. These
fragments can be the dissociating molecular prestellar cores under UV erosion.
Alternatively, such fragments may be formed in the PDR by impulse heating and
turbulent motions.

  5. 5.

We detect an H i sub-bubble inside the H i shell. A group of young stellar
objects is detected in the direction of the walls, their formation may be
induced by the interaction between the bubble and the shell. The sub-bubble is
young with A kinematical age of $\sim$105 years.

  6. 6.

Two molecular cores (NE and SE) are present in the S187 environment. They have
similar masses ($\sim$1200 M⊙ and $\sim$900 M⊙, respectively) and different
physical conditions. S187 SE contains a number of young stellar objects, an
IRAS source, outflows and a number of other signs of ongoing star formation,
which is an argument of A past interaction with the H ii zone. The spatial-
kinematic relation between the H i and 12CO emission implies ongoing
interaction between the atomic shell and the molecular core, but no shock
signatures or PDR AT the border are observed. The S187 SE virial parameter is
$\sim$2.75. S187 NE is being compressed by the shell and its virial parameter
is $\sim$1.4. No significant star formation activity is detected inside the
core. Three YSOs are present inside the H i shell contact zone. The
distribution of molecular lines suggests heating and the increase of
turbulence near the contact zone. The core periphery located towards the H i
shell and the H ii zone is eroded by the radiation field.

  7. 7.

The WISE [3.4]$-$[4.6] colour correlates with the column density of the
molecular gas. One of the most reasonable explanation is the effect of photon
absorption by grains in the dense molecular gas of S187 SE.

  8. 8.

The H i line radial plots and position-velocity diagrams are different for the
directions where the H ii region is bounded by dense gas and where it is not.
The PV diagram has an O-type shape where the shell is density bounded. In the
regions where no dense structures are detected, the PV structures are
irregular.

## Acknowledgements

We are thankful to Maria Murga, Yaroslav Pavlyuchenkov and Maria Kirsanova for
helpful discussions and consultations. We appreciate by Kim Ardvisson and Cris
Ardvisson for providing CO isotopologues observations data. We thank the
anonymous reviewer for a critical reading of the manuscript and the useful and
thougtful comments and suggestions, which greatly improved the contents of the
paper.

This work was supported by the Russian Science Foundation under project No.
17-12-01256 (the H i line, molecular observations and analysis, physical
parameter estimations). The dust emission analysis was supported by the
Foundation for the Advancement of Theoretical Physics and Mathematics “BASIS”.

The Canadian Galactic Plane Survey (CGPS) is a Canadian project with
international partners. The Dominion Radio Astrophysical Observatory is
operated as a national facility by the National Research Council of Canada.
The Five College Radio Astronomy Observatory CO Survey of the Outer Galaxy was
supported by NSF grant AST 94-20159. The CGPS is supported by a grant from the
Natural Sciences and Engineering Research Council of Canada.".

We acknowledge support from the Onsala Space Observatory for the provisioning
of its facilities/observational support for obtaining the HCO+ data. The
Onsala Space Observatory national research infrastructure is funded through
Swedish Research Council grant No 2017-00648.

We thank the staff of the GMRT who have made these observations possible. DKO
and SKG acknowledge the support of the Department of Atomic Energy, Government
of India, under project Identification No. RTI 4002.

## Data Availability

The raw data underlying this article are available in the GMRT Online Archive
at https://naps.ncra.tifr.res.in/goa/data/search, and can be accessed with
Project Code 26_012.

The processed data underlying this article (combined GMRT and CGPS datacube)
are available in Zenodoo, at doi:10.5281/zenodo.5070398

## References

  * Aniano et al. (2011) Aniano G., Draine B. T., Gordon K. D., Sandstrom K., 2011, PASP, 123, 1218
  * Arvidsson & Kerton (2011) Arvidsson K., Kerton C. R., 2011, AJ, 141, 153
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Bally et al. (2018) Bally J., Chambers E., Guzman V., Keto E., Mookerjea B., Sandell G., Stanke T., Zinnecker H., 2018, AJ, 155, 80
  * Belitsky et al. (2015) Belitsky V., et al., 2015, A&A, 580, A29
  * Berry (2015) Berry D. S., 2015, Astronomy and Computing, 10, 22
  * Burton (1971) Burton W. B., 1971, A&A, 10, 76
  * Cichowolski et al. (2009) Cichowolski S., Romero G. A., Ortega M. E., Cappa C. E., Vasquez J., 2009, MNRAS, 394, 900
  * Cichowolski et al. (2014) Cichowolski S., Pineault S., Gamen R., Arnal E. M., Suad L. A., Ortega M. E., 2014, MNRAS, 438, 1089
  * Cutri et al. (2003) Cutri R. M., et al., 2003, 2MASS All Sky Catalog of point sources.
  * Decataldo et al. (2019) Decataldo D., Pallottini A., Ferrara A., Vallini L., Gallerani S., 2019, MNRAS, 487, 3377
  * Draine & Li (2007) Draine B. T., Li A., 2007, ApJ, 657, 810
  * Elmegreen & Falgarone (1996) Elmegreen B. G., Falgarone E., 1996, ApJ, 471, 816
  * Engels & Bunzel (2015) Engels D., Bunzel F., 2015, A&A, 582, A68
  * Evans et al. (2020) Evans I. N., et al., 2020, in American Astronomical Society Meeting Abstracts #235. p. 154.05
  * Fontani et al. (2012) Fontani F., Giannetti A., Beltrán M. T., Dodson R., Rioja M., Brand J., Caselli P., Cesaroni R., 2012, MNRAS, 423, 2342
  * Gorti & Hollenbach (2002) Gorti U., Hollenbach D., 2002, ApJ, 573, 215
  * Han et al. (1995) Han F., et al., 1995, Publications of the Purple Mountain Observatory, 14, 185
  * Henkel et al. (1986) Henkel C., Haschick A. D., Guesten R., 1986, A&A, 165, 197
  * Israel (1977) Israel F. P., 1977, A&A, 61, 377
  * Jackson et al. (2006) Jackson J. M., et al., 2006, ApJS, 163, 145
  * Joncas et al. (1992) Joncas G., Durand D., Roger R. S., 1992, ApJ, 387, 591
  * Kang & Kerton (2012) Kang S.-J., Kerton C. R., 2012, ApJ, 759, 13
  * Kang et al. (2017) Kang S.-J., Kerton C. R., Choi M., Kang M., 2017, ApJ, 845, 21
  * Kauffmann et al. (2013) Kauffmann J., Pillai T., Goldsmith P. F., 2013, ApJ, 779, 185
  * Kirsanova & Wiebe (2019) Kirsanova M. S., Wiebe D. S., 2019, MNRAS, 486, 2525
  * Kirsanova et al. (2020) Kirsanova M. S., et al., 2020, MNRAS, 497, 2651
  * Koenig & Leisawitz (2014) Koenig X. P., Leisawitz D. T., 2014, ApJ, 791, 131
  * Könyves et al. (2020) Könyves V., et al., 2020, A&A, 635, A34
  * Kritsuk et al. (2013) Kritsuk A. G., Lee C. T., Norman M. L., 2013, MNRAS, 436, 3247
  * Larson (1981) Larson R. B., 1981, MNRAS, 194, 809
  * Leahy & Tian (2006) Leahy D., Tian W., 2006, A&A, 451, 251
  * Mangum & Shirley (2015) Mangum J. G., Shirley Y. L., 2015, PASP, 127, 266
  * McMullin et al. (2007) McMullin J. P., Waters B., Schiebel D., Young W., Golap K., 2007, in Shaw R. A., Hill F., Bell D. J., eds, Astronomical Society of the Pacific Conference Series Vol. 376, Astronomical Data Analysis Software and Systems XVI. p. 127
  * Meisner & Finkbeiner (2014) Meisner A. M., Finkbeiner D. P., 2014, ApJ, 781, 5
  * Meisner et al. (2017) Meisner A. M., Lang D., Schlegel D. J., 2017, AJ, 154, 161
  * Newville et al. (2014) Newville M., Stensitzki T., Allen D. B., Ingargiola A., 2014, LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python, doi:10.5281/zenodo.11813, https://doi.org/10.5281/zenodo.11813
  * Nguyen et al. (2019) Nguyen H., Dawson J. R., Lee M.-Y., Murray C. E., Stanimirović S., Heiles C., Miville-Deschênes M. A., Petzler A., 2019, ApJ, 880, 141
  * Okada et al. (2019) Okada Y., Güsten R., Requena-Torres M. A., Röllig M., Stutzki J., Graf U. U., Hughes A., 2019, A&A, 621, A62
  * Palagi et al. (1993) Palagi F., Cesaroni R., Comoretto G., Felli M., Natale V., 1993, A&AS, 101, 153
  * Pavlyuchenkov et al. (2012) Pavlyuchenkov Y. N., Wiebe D. S., Akimkin V. V., Khramtsova M. S., Henning T., 2012, MNRAS, 421, 2430
  * Pedregosa et al. (2011) Pedregosa F., et al., 2011, Journal of Machine Learning Research, 12, 2825
  * Peeters et al. (2004) Peeters E., Allamandola L. J., Bauschlicher C. W. J., Hudgins D. M., Sandford S. A., Tielens A. G. G. M., 2004, ApJ, 604, 252
  * Pineda et al. (2013) Pineda J. L., Langer W. D., Velusamy T., Goldsmith P. F., 2013, A&A, 554, A103
  * Pirogov & Zemlyanukha (2021) Pirogov L. E., Zemlyanukha P. M., 2021, Astronomy Reports, 65, 82
  * Rau et al. (2019) Rau U., Naik N., Braun T., 2019, AJ, 158, 3
  * Reach & Heiles (2021) Reach W. T., Heiles C., 2021, ApJ, 909, 71
  * Reach et al. (2017) Reach W. T., Bernard J.-P., Jarrett T. H., Heiles C., 2017, ApJ, 851, 119
  * Requena-Torres et al. (2016) Requena-Torres M. A., Israel F. P., Okada Y., Güsten R., Stutzki J., Risacher C., Simon R., Zinnecker H., 2016, A&A, 589, A28
  * Richards et al. (2012) Richards E. E., Lang C. C., Trombley C., Figer D. F., 2012, AJ, 144, 89
  * Roger et al. (2004) Roger R. S., McCutcheon W. H., Purton C. R., Dewdney P. E., 2004, A&A, 425, 553
  * Roman-Duval et al. (2010) Roman-Duval J., Jackson J. M., Heyer M., Rathborne J., Simon R., 2010, ApJ, 723, 492
  * Russeil et al. (2007) Russeil D., Adami C., Georgelin Y. M., 2007, A&A, 470, 161
  * Saha et al. (2018) Saha P., Roy N., Bhattacharya M., 2018, MNRAS, 480, L126
  * Salas et al. (1998) Salas L., Cruz-González I., Porras A., 1998, ApJ, 500, 853
  * Shirley (2015) Shirley Y. L., 2015, PASP, 127, 299
  * Slysh et al. (1999) Slysh V. I., Valtts I. E., Kalenskii S. V., Voronkov M. A., Palagi F., Tofani G., Catarzi M., 1999, A&AS, 134, 115
  * Snell & Bally (1986) Snell R. L., Bally J., 1986, ApJ, 303, 683
  * Speagle (2019) Speagle J. S., 2019, arXiv e-prints, p. arXiv:1909.12313
  * Sternberg & Dalgarno (1995) Sternberg A., Dalgarno A., 1995, ApJS, 99, 565
  * Sternberg et al. (2014) Sternberg A., Le Petit F., Roueff E., Le Bourlot J., 2014, ApJ, 790, 10
  * Sunada et al. (2007) Sunada K., Nakazato T., Ikeda N., Hongo S., Kitamura Y., Yang J., 2007, PASJ, 59, 1185
  * Taylor et al. (2003) Taylor A. R., et al., 2003, AJ, 125, 3145
  * Thompson et al. (2004) Thompson M. A., White G. J., Morgan L. K., Miao J., Fridlund C. V. M., Huldtgren-White M., 2004, A&A, 414, 1017
  * Tielens & Hollenbach (1985) Tielens A. G. G. M., Hollenbach D., 1985, ApJ, 291, 722
  * Tody (1986) Tody D., 1986, in Crawford D. L., ed., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 627, Instrumentation in astronomy VI. p. 733, doi:10.1117/12.968154
  * Topchieva et al. (2017) Topchieva A., Wiebe D., Kirsanova M., Krushinsky V., 2017, in Balega Y. Y., Kudryavtsev D. O., Romanyuk I. I., Yakunin I. A., eds, Astronomical Society of the Pacific Conference Series Vol. 510, Stars: From Collapse to Collapse. p. 98
  * Valdettaro et al. (2001) Valdettaro R., et al., 2001, A&A, 368, 845
  * Verstraete et al. (2001) Verstraete L., et al., 2001, A&A, 372, 981
  * Weaver et al. (1977) Weaver R., McCray R., Castor J., Shapiro P., Moore R., 1977, ApJ, 218, 377
  * Wilson et al. (2012) Wilson T. L., Rohlfs K., Huttemeister S., 2012, Tools of Radio Astronomy, 5th edition
  * Wright et al. (2010) Wright E. L., et al., 2010, AJ, 140, 1868
  * Yamamura et al. (2010) Yamamura I., Makiuti S., Ikeda N., Fukuda Y., Oyabu S., Koga T., White G. J., 2010, VizieR Online Data Catalog, p. II/298
  * Yan et al. (2021) Yan Q.-Z., Yang J., Sun Y., Su Y., Xu Y., Wang H., Zhou X., Wang C., 2021, A&A, 645, A129
  * Ybarra et al. (2014) Ybarra J. E., Tapia M., Román-Zúñiga C. G., Lada E. A., 2014, ApJ, 794, L25
  * Zavagno et al. (1994) Zavagno A., Deharveng L., Caplan J., 1994, in The P. S., Perez M. R., van den Heuvel E. P. J., eds, Astronomical Society of the Pacific Conference Series Vol. 62, The Nature and Evolutionary Status of Herbig Ae/Be Stars. p. 347
  * Zinchenko (2000) Zinchenko I. I., 2000, Astronomy Letters, 26, 802
  * Zinchenko et al. (1998) Zinchenko I., Pirogov L., Toriseva M., 1998, A&AS, 133, 337
  * Zinchenko et al. (2009) Zinchenko I., Caselli P., Pirogov L., 2009, MNRAS, 395, 2234

##  Appendix A k-Nearest Neighbor regression for H i maps fitting

This appendix describes the usage of the k-Nearest Neighbors regression for
the analysis of the spectral maps using the Gaussian or simple radiative
transfer models. The common way is to perform individual point-by-point least
squares fitting of the model parameters. The major problem of this approach is
that the results of the fitting are heavily dependent on the initial
conditions of the fitting procedure. The H i spectra can consist of a number
of peaks represented by individual Gaussians which lead to cases where no
single statistically important solution exists. This effect increases with
decreasing signal-to-noise ratio. The model parameters and it’s confidence
regions can be found using the sampling of the likelihood function by Monte
Carlo Markov Chains (Speagle, 2019), by grid sampling of dimensionally-reduced
space (Pirogov & Zemlyanukha, 2021) or other similar sampling methods. Such
approach is problematic for spectral map analysis OF individual positions. The
computational complexity of sampling is quite high (especially for high
dimension models) which should be repeated for individual positions and
likelihood functions.

In the case of S187 H i data there is a position on the map where the
parameters are known by analysis of the absorption spectra. Such information
can be used for the map analysis. We assume that the physical parameters are
spatially-correlated, and the results of the previous fittings can be used as
the initial parameters for gradient descent (or similar) methods. We define
the residual:

$\chi^{2}=\frac{1}{N_{v}-N_{p}}\sum_{i=1}^{N_{v}}{\frac{(I_{i}^{obs}-I_{i}^{mod})^{2}}{\sigma^{2}}}\hskip
5.69054pt,$ (5)

where $N_{p}$ is the number of parameters in the model, $N_{v}$ is the number
of spectral channels used for fitting. $I_{i}^{obs}$ is the observed intensity
of the spectral channel $i$ and $I_{i}^{mod}$ is the model intensity (from the
Eq. 2). We use lmfit (Newville et al., 2014) procedure and Levenberg–Marquardt
algorithm to find the minimum of the residual function. In general, the
initial conditions for the fitting procedure can be approximated from previous
calculations in different positions of the map. We use the k-Nearest Neighbor
(kNN) regression as the approximation:

$p_{k}=\frac{\sum_{j=1}^{m}p_{j}d(I_{k},I_{j})^{2}}{\sum_{j=1}^{m}d(I_{k},I_{j})^{2}}\hskip
5.69054pt,$ (6)

where $p_{k}$ is the model parameter value, $d(I_{1},I_{2})$ is the distance
between two spectra defined as:

$d(I_{1},I_{2})=\sqrt{{N_{v}^{-1}\sigma_{1}\sigma_{2}}\times{\sum_{j=1}^{N_{v}}(I_{1j}-I_{2j})^{2}}}\hskip
5.69054pt,$ (7)

where m is the number of Nearest Neighbors used for the estimation. We use the
Euclidean metric in neighbors space. This is equal to a weighted average of
the model parameters from $k$ samples. These $k$ samples are chosen by minimal
squared deviation from current spectra. We normalize the distance by the
number of spectral channels and the RMS deviation estimated in line-free
channels. The effect of the number of channels is not important since it is
the same for the map. The distance to noisier data is increased by multiplying
BY $\sigma$. We use the procedures from scikit-learn (Pedregosa et al., 2011)
in our algorithm. Using this approach we reduce the initial conditions problem
to the adequate sample generation. The absorption spectra decomposition
results and the Least Squares fitting results was used to fill the samples. We
repeat map analysis a number of times until the results converge. The
algorithm used to estimate the physical parameters of the map follows this
sequence:

  1. 1.

Estimate the optical depth and spin temperature profile from the absorption
and emission spectra.

  2. 2.

Fit the profile with the flat-layer model.

  3. 3.

Fit the emission spectra to the nearest positions with the initial conditions
from the previous step. Mark the positions as fitted and put the observed
spectra and the fitted results in the dataset.

  4. 4.

Fit the nearest unfitted positions with the initial conditions estimated using
kNN and the dataset. Mark the positions as fitted.

  5. 5.

Repeat the previous step until the full map is fitted.

  6. 6.

Repeat the previous two steps starting from positions near the absorption:
refit every position with new initial conditions from kNN results. Repeat
until the fitted results for individual positions converge.

The aim of such iterations is to fill the training set for kNN regression. We
repeat the whole map analysis until the initial conditions estimated by kNN
for the least-squares fitting will be close to the fitted results. The initial
conditions obtained by averaging the results from the similar spectra also
reduce the computation cost. The insufficient sampling of the training set on
the first iteration is negated by the next iterations. Such an approach is
relatively computationally efficient because the initial conditions are close
to the fitted results and the gradient descent (or similar) methods need a
large number of iterations to fall into a residual minima. The only assumption
is that the same physical properties produce the same spectra. Still, the
interval of the confident estimation in individual fit is important and can be
used for the analysis of the regions where the estimation is reasonable and
where it is not.

##  Appendix B Figures

Figure 22: The channel maps of the S187 environment. The 12CO (1-0) emission
is in black contours. The combined GMRT+CGPS H i channel maps are shown in the
background image. The beams are shown in lower-left corner of the first
subplot. The velocity ranges are labelled at the top-right of every subplot.
The integrated 12CO (1–0) line emission is shown in the last subplot. The
fragments identified by Fellwalker algorithm (See sect. 4.1.1) is shown as
white crosses. Figure 23: The channel maps of the S187 environment. The C18O
(1-0) emission is in black contours. The combined GMRT+CGPS H i channel maps
are shown in the background image. The beams are shown in lower-left corner of
the first subplot. The velocity ranges are labelled at the top-right of every
subplot. The penultimate last subplot represents integrated emission of C18O
(1–0) and the 21-cm line integrated emission image is shown in the last
subplot.  Figure 24: Same as for Fig. 22, but for 13CO (1–0)