# The link between gas and stars in the S254-S258 star-forming region

D. A. Ladeyshchikov1, M. S. Kirsanova2, A. M. Sobolev1, M. Thomasson3, V.
Ossenkopf-Okada4, M. Juvela5, S. A. Khaibrakhmanov1, E. A. Popova1
1 Ural Federal University, 51 Lenin Str., Ekaterinburg 620051, Russia
2 Institute of Astronomy, Russian Academy of Sciences, 48 Pyatnitskaya Str.,
Moscow 119017, Russia
3 Department of Space, Earth and Environment, Chalmers University of
Technology, Onsala Space Observatory, SE-43992 Onsala, Sweden
4 I. Physik. Institut, University of Cologne, D-50937 Cologne, Germany
5 University of Helsinki, Department of Physics, PL 64 (Gustaf Hällströmin
katu 2), Finland<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The paper aims to study relation between the distributions of the young
stellar objects (YSOs) of different ages and the gas-dust constituents of the
S254-S258 star-formation complex. This is necessary to study the time
evolution of the YSO distribution with respect to the gas and dust compounds
which are responsible for the birth of the young stars. For this purpose we
use correlation analysis between different gas, dust and YSOs tracers. We
compared the large-scale CO, HCO+, near-IR extinction, and far-IR Herschel
maps with the density of YSOs of the different evolutionary Classes. The
direct correlation analysis between these maps was used together with the
wavelet-based spatial correlation analysis. This analysis reveals a much
tighter correlation of the gas-dust tracers with the distribution of Class I
YSOs than with that of Class II YSOs. We argue that Class I YSOs which were
initially born in the central bright cluster S255-IR (both N and S parts)
during their evolution to Class II stage ($\sim$2 Myr) had enough time to
travel through the whole S254-S258 star-formation region. Given that the
region contains several isolated YSO clusters, the evolutionary link between
these clusters and the bright central S255-IR (N and S) cluster can be
considered. Despite the complexity of the YSO cluster formation in the non-
uniform medium, the clusters of Class II YSOs in the S254-258 star-formation
region can contain objects born in the different locations of the complex.

###### keywords:

stars: protostars – ISM: clouds – H ii regions – dust, extinction

††pubyear: 2021††pagerange: The link between gas and stars in the S254-S258
star-forming region–C

##  1 Introduction

S254-S258 is a region that shows signs of sequential star formation triggered
by expanding H II regions (see, e.g., Bieging et al. (2007); Minier et al.
(2007); Chavarría et al. (2008); Ojha et al. (2011)). It contains two bright
optical H II regions S255 and S257 with angular sizes about 3′, an extended H
II region S254 with size about 10′ and two compact regions S256 and S258 with
sizes about 1′. The region contains 10 young star clusters that were
identified in the work of Chavarría et al. (2008) (see Figure 12 therein):
S255-2 & S255N, S256, S258, G192.54-0.15, G192.75-0.00, G192.75-0.08,
G192.63-0.00, G192.69-0.25, G192.65-0.08, G192.55-0.01. At present, only the
brightest molecular condensations with embedded clusters S255N, S255-IR, and
S255S were mapped with high-density gas tracing molecular lines (Cyganowski et
al., 2007; Wang et al., 2011; Zinchenko et al., 2012; Zinchenko et al., 2015;
Zemlyanukha et al., 2018). The expansion of the H II regions S255 and S257
(observed recently by e. g. Buslaeva et al., 2021) and compressing the
material "sandwiched" between these H II regions may explain a large amount of
high-density gas and active star formation in the S255-IR region (Minier et
al., 2007). Other molecular condensations with embedded clusters have not been
studied yet in the high-density gas tracers. Most of them reside in the
primordial gas (i.e. in the gas which has not been compressed by expanding H
II regions) of the S254-S258 star formation region mapped by Bieging et al.
(2007); Chavarría et al. (2008). Two of the condensations (S256-south and
S258) are located towards the S256 and S258 H II regions, while others,
including G172.75-0.00, G192.75-0.08, and G192.69-0.25, are located on the
periphery of the region.

The presence of the high-density gas in the position of stellar clusters can
be checked by observations of e.g. CS and HCO+ molecular line emission (Gaches
et al., 2015). The CS(2–1) line in the brightest condensations S255N and
S255-IR was firstly detected by Zinchenko et al. (1994) and was observed
several times in the recent years (Zinchenko et al., 2009; Wang et al., 2011;
Zinchenko et al., 2015) together with the HCO+ line (Zinchenko et al., 2009;
Zinchenko et al., 2015; Shirley et al., 2013). The aim of the present study is
to map the CS(2–1) and HCO+(1–0) emission of the surrounding gas, including
mapping of isolated star clusters in order to study its spatial structure and
kinematics. We also observed optically thin lines C18O(1–0), C34S(2–1), and
H13CO+(1–0) toward selected directions to reveal kinematics of the inner parts
of the clumps with high extinction. We additionally compare different gas and
dust tracers to the distribution of the young stellar objects (YSOs) density
to reveal the connection between young star clusters and the surrounding gas
using the correlation analysis method. Class I and class II YSOs represent
different evolutionary stages of the objects with protoplanetary disks. So,
these studies provide constraints on the distribution of exoplanet-hosting
stars.

##  2 Observations and processing

###  2.1 High-density molecules observations

We use the observations of HCO+, CS and C18O molecular line emission taken in
February 2009 with the Onsala Space Observatory 20m radio telescope (PI –
Maria Kirsanova). The SIS receiver with the frequency range 85-116 MHz was
used for the observations. Extended mapping was done in CS(2–1) and HCO+(1–0)
lines toward the molecular ridge with a size of about $6$′$\times 15$′. Fig. 1
shows the maps of the HCO+(1–0) and CS(2–1) integrated emission together with
spectra in several selected directions. Additionally to these maps we observed
optically thin lines C34S(2–1), C18O(1–0), and H13CO+(1–0) toward ten selected
directions, chosen to trace the inner and outer part of the central bright
S255-IR cluster. These directions are marked by crosses in Figure 1.
Frequencies of the observed lines are given in Table 1. The total amount of
the telescope time was 162 hours. Observations were taken with 40 MHz
bandwidth and 1600 channels using frequency switching mode with a throw of 20
MHz. The spectral resolution corresponds to 0.07 km s-1. The full width at
half-maximum (FWHM) of the telescope beam ranges from $34\arcsec$ to
$43\arcsec$ and is listed in Table 1. We used the beam size as spacing for the
extended maps.

(a) CS (2–1) integrated intensity (b) HCO+ (1–0) integrated intensity

Figure 1: Integrated intensity maps of the high density tracers in S254-258. Red circles represent visible radius of H II regions from DSS red images. Spectra for several positions are shown in the sub-figures. (a) HCO+ (1–0) map, contours go from 0.6 K km s-1 to 5.6 K km s-1 in steps of 0.45 K km s-1. (b) CS(2–1) map, contours go from 0.56 K km s-1 to 9.13 K km s-1 in steps of 0.78 K km s-1.  Table 1: Characteristics of molecular line observations in the S255-258 region. $\theta_{\mathrm{FWHM}}$ is the beam size of the telescope at the specified frequency. Molecule | Line | Frequency | $\theta_{\mathrm{FWHM}}$ | Observation type
---|---|---|---|---
|  | (MHz) | (arcsec) |
CS | (2–1) | 97980.9 | 38.53 | Mapping
HCO+ | (1–0) | 89188.5 | 42.32 | Mapping
C34S | (2–1) | 96412.9 | 39.15 | Points
C18O | (1–0) | 109782.1 | 34.38 | Points
H13CO+ | (1–0) | 86754.2 | 43.51 | Points

The CLASS program from the GILDAS software
package111http://www.iram.fr/IRAMFR/GILDAS was used to process the obtained
spectra. All spectra were converted to the Tmb scale using a primary beam
efficiency. The beam efficiency value varies from 0.35 to 0.55, depending on
the source elevation at the time of observation. We removed a 3rd order
baseline from the raw data to obtain the continuum-subtracted spectra. The
system temperature varies from 200 to 600 K, but in few cases exceeded 800 K.
The spectra were averaged using the inverse system temperature as the
weighting function. As we performed frequency-switching mode observations, we
used the FOLD command of CLASS to deconvolve the frequency-switched spectra.

To estimate column density of CS and HCO+ molecules, their excitation
temperature and optical depth needs to be known. This requires the knowledge
of at least two lines. In our case, we have maps in HCO+(1–0) and CS(2–1), but
H13CO+(1–0) and C34S(2–1) data are available only at six positions. Therefore,
we can only calculate the column densities in the directions of the available
C34S(2–1) and H13CO+(1–0) spectra.

We obtained the line peak optical depth values ($\tau_{0}$) using the ratios
of the peak main-beam temperatures of HCO+(1–0) and CS(2–1) lines and
corresponding isotopes ($\mathrm{H^{13}CO^{+}}$(1–0), C34S(2–1)). The same
equation was used for calculation of the line-integrated optical depths
($\int\tau$) using the integrated main-beam temperatures:

$\frac{T_{\mathrm{mb[isotope]}}}{T_{\mathrm{mb[main]}}}=\frac{1-e^{-\tau_{0}/R}}{1-e^{-\tau_{0}}};\,\,\frac{\int
T_{\mathrm{mb[isotope]}}}{\int
T_{\mathrm{mb[main]}}}=\frac{1-e^{-\int\tau/R}}{1-e^{-\int\tau}},$ (1)

where the isotopic ratio $R$ is assumed
$R(^{34}\mathrm{S}/^{32}\mathrm{S})=22$ and
$R(^{12}\mathrm{C}/^{13}\mathrm{C})=80$. These values are consistent with the
12C/13C ratio of 79.2 (Milam et al., 2005) toward S255 region at the
galactocentric distance of 9.75 kpc.

The excitation temperature ($T_{\mathrm{ex}}$) was estimated using the line
peaks intensities and the optical depths ($\tau_{0}$) from the equation of
radiative transfer:

$T_{\mathrm{mb}}=(J_{v}(T_{\mathrm{ex}})-J_{v}(T_{\mathrm{bg}}))(1-\exp(-\tau_{0}))$
(2)

where $J_{v}(T)=T_{0}/(\exp(T_{0}/T)-1)$ and $T_{\mathrm{bg}}=2.7~{}$K using
the coefficient $T_{0}=h\nu/k$. The used rest frequencies for different lines
are listed in the Table 1.

The column densities of CS and HCO+ were then calculated using the approach
described by Mangum & Shirley (2015).

$\displaystyle N=\frac{\int\tau}{1-\exp(\int\tau)}\frac{\rm
3h}{8\pi^{3}S\mu^{2}}\,\frac{Q_{\mathrm{rot}}}{g_{\mathrm{J}}}\,\frac{\exp(\frac{E_{\rm
u}}{{\rm k}T_{\mathrm{ex}}}))}{\exp(\frac{\rm h\nu}{{\rm
k}/T_{\mathrm{ex}}})-1}\,\frac{1}{f}\int\tau\mathrm{d}v$ (3)

where $S=J^{2}/(J[2J+1])$, $g_{J}=2J+1$ is statistical weight,
$Q_{\mathrm{rot}}\simeq{\rm k}T_{\mathrm{ex}}/{\rm h}/B_{\mathrm{rot}}+1/3$ is
partition function for linear molecules. We assume a filling factor $f=1$. The
constants $B_{\mathrm{rot}}$, $E_{\mathrm{u}}$, and $\mu$ were taken form the
Cologne Database for Molecular Spectroscopy (Müller et al., 2001) and given in
Table 2.

Table 2: Constants that was used for column density calculation from the molecular line observations. Line | T0 | Brot | Eu | $\mu$
---|---|---|---|---
| K | GHz | K | 10-18 esu
CO(1–0) | 5.532 | 57.6 | 5.53 | 0.11
CO(2–1) | 11.064 | 57.6 | 16.60 | 0.11
13CO(1–0) | 5.288 | 55.1 | 5.29 | 0.11
13CO(2–1) | 10.577 | 55.1 | 15.87 | 0.11
CS(2–1) | 4.702 | 24.5 | 7.10 | 1.96
C34S(2–1) | 4.627 | 24.1 | 6.3 | 1.96
HCO+(1–0) | 4.280 | 44.6 | 4.28 | 3.90
H13CO+(1–0) | 4.163 | 43.3 | 4.16 | 3.90

###  2.2 CO data

We used 13CO(1–0) and 13CO(2–1) data to study the whole molecular cloud
associated with S254-258. The estimation of CO column density was done using
the same technique that is described in the Section 2.1. For CO molecule we
have $S=J^{2}/([2J+1])$, as described in the Mangum & Shirley (2015). The data
on 12CO(1–0) and 13CO(1–0) line emission is taken from Chavarría et al.
(2008), data on 12CO(2–1) and 13CO(2–1) lines from Bieging et al. (2009). The
optical depth of CO lines (both peak and integrated) were estimated from
12CO/13CO line ratio (peak and integrated) using equation (1) and the
abundance ratio $R=80$.

The excitation temperature is determined from optically thick 12CO line. For
$\tau\gg 1$ the radiative transfer equation (2) can be inverted to provide

$\displaystyle T_{\rm ex}=T_{0}/\ln\left(1+\frac{T_{0}}{T_{\rm
B}^{12}+T_{0}/(e^{T_{0}/T_{\mathrm{bg}}}-1)}\right)$ (4)

where $T_{\rm B}^{12}$ is the peak 12CO brightness temperature.

We considered the use of gas and dust clump-extraction algorithms like
clumpfind and gaussclump with available CO data (Chavarría et al., 2008;
Bieging et al., 2009). However, they would be limited by the relatively large
velocity resolution of the available data ($1.2$ km s-1) that does not allow
to distinguish kinematically different molecular clumps.

####  2.2.1 Difference between CO(1-0) and CO(2-1) data

Comparing the obtained values of line optical depth and excitation
temperatures for the (1-0) and (2-1) transitions, we have found that
excitation temperatures for both transitions are quite close (see Table 5),
but the (2-1) transition has systematically lower excitation temperatures in
comparison to the (1-0) transition. That can be due to radiation transfer
effects as explained in details in Bernes (1979) (see Figure 1 therein) – the
(2-1) transition traces more outer layers of the molecular cloud, less of the
volume. For example, in Figure 1 of Bernes (1979) the (2-1) excitation
temperature drops from 19 K in the cloud center to 12 K in the outer layers,
while the excitation temperature (1-0) drops only from 19 K to 18 K. Thus, the
excitation temperature of the (1-0) transition is much more close to the
average molecular cloud’s kinetic temperature than the excitation temperature
of the (2-1) transition. The obtained values of optical depth for (2-1) and
(1-0) reveal that the optical depth of (2-1) transition is 1.5-3 times larger
than (1-0) transition consistent with the fact that the (2-1) line rather
traces the surface. Due to lower optical depth of the (1-0) transition it
allows for a more accurate determination of excitation temperature being less
sensitive to radiation transfer effects. Therefore we use only the (1-0)
transition to further study the molecular gas column densities.

We calculate the column density of CO molecules with Equation (3), where we
use the excitation temperature and the line-integrated optical depth of the
(1-0) line. The constants of 13CO molecule that were used for column density
calculation are presented in the Table 2. We convert N(CO) to N(H2) using
abundance [CO]/[H2] = $8\times 10^{-5}$ (Blake et al., 1987). The distance to
the S254-S258 region is assumed $1.59\pm 0.07$ kpc (Rygl et al., 2010), as
measured by trigonometric parallax of 6.7 GHz methanol maser in the S255
region. We used this distance for the mass estimates in the current and
subsequent sections.

###  2.3 Extinction maps

Extinction maps were produced using the NICEST method described by Lombardi
(2009). We used the implementation of Juvela & Montillaud (2016a).

For constructing the extinction maps we used the data from the UKIRT Infrared
Deep Sky Survey (UKIDSS, Lawrence et al., 2007) together with data from the
2MASS survey (Skrutskie et al., 2006). The UKIDSS survey instrument is WFCAM
(Casali et al., 2007) mounted on the UK Infrared Telescope (UKIRT) in Hawaii.
The data from the WFCAM camera have a high quality – resolution of 0.4 arcsec,
photometric accuracy of $\sigma_{J}<0.02^{m}$ and a deep limiting magnitude up
to $19.9^{m}$ in the $J$ band. The S254-S258 star formation region lies in the
Galactic Plane Survey (GPS) area of UKIDSS. We constructed the extinction map
using combined UKIDSS and 2MASS data centered on the S255N region with a
radius of 1.5 degrees.

Due to the inaccuracy of the UKIDSS catalog for bright stars, we replaced the
UKIDSS photometry by 2MASS for stars brighter than 12m for each band using the
TOPCAT software (Taylor, 2005). Photometry of 2MASS was corrected to the
UKIDSS photometric standard using equations (6-8) from Hodgkin et al. (2009).

We filtered the UKIDSS catalogue using the following criteria:

  1. 1.

Included sources are only those with a high probability of being a star
(column “pstar” more than $0.8$ from the WFCAM catalogue);

  2. 2.

At least one combination of two colours (J-H, J-K or H-K) available for each
source;

  3. 3.

Sources that did not have at least two bands with magnitude uncertainties
below 1m were rejected;

The NICEST method is based on the comparison between observed colours between
the J, H and K bands and intrinsic colours of stars, corresponding to zero
extinction.

The NICEST maps were obtained at the resolutions (FWHM) of 1.0-3.0′. The pixel
size was set to FWHM/3. We used an extinction curve with
$E(J-H)/A_{V}=0.10167$ and $E(H-K)/A_{V}=0.06347$ corresponding to the
standard $R_{V}=3.1$ extinction curve (Cardelli et al., 1989) with
$A_{J}/A_{V}=0.26$. We convert the extinction map to the gas column density
using the approach described in the Appendix A. The used conversion factor is
$N(\mathrm{H})/A_{V}=1.87\pm 0.13\times
10^{21}~{}(\mathrm{cm}^{-2}\mathrm{mag^{-1})}$. The NICEST AJ extinction maps
at the best available resolution (1.0 arcmin) in the units of gas column
density are presented in Fig. 2.

The zero level of the NICER extinction map was set from a low-column-density
region that is centered at RA${}_{\mathrm{J2000}}=06^{h}15^{m}24.9^{s}$,
Dec${}_{\mathrm{J2000}}=18\degr 33\arcmin 45\arcsec$ and has a size of
$40\arcmin\times 20$′.

###  2.4 Herschel data

We study dust temperature and column density using the 70-500 µm Hi-GAL data
(Molinari et al., 2016). We used two different data sets: the ViaLactea
results (Marsh et al., 2017) based on the 70-500 µm data and our own fit using
only 160-500 µm data. When comparing the ViaLactea total column density maps
with the extinction maps we found that the values of the gas column density in
the ViaLactea map is lower than the values from the extinction by a factor of
3–5, therefore we tried an independent fit. We believe that the discrepancy
could be due to the missing Planck correction in ViaLactea map for the
Herschel-PACS data making them inconsistent with the Herschel-SPIRE data in
particular for the shortest wavelengths where the PLANCK data cannot be easily
extrapolated. The detailed correction method for Herschel-PACS data using the
Planck data can be found in Lombardi et al. (2014); Abreu-Vicente et al.
(2017) for large-scale maps and in Sadavoy et al. (2018) for small-scale maps.
However, in this work we excluded the 70-100 µm Herschel flux and performed a
grey-body SED fit with a dust emissivity $\kappa\propto\lambda^{-1.8}$ only to
the 160-500 µm range. The resulting fit gives column densities and masses of
the clumps that agree with CO-based values and the extinction maps, however it
is not sensitive to warm dust ($>35$K).

Therefore we used the ViaLactea results to create the contour levels for
clumps and study the dust temperature of clumps, while the 160-500 µm fit was
used to estimate the column density and masses of the clumps and analyze the
correlations between dust and YSOs density maps.

Using the ViaLactea map instead of the 160-500 µm map gives almost no
difference in the correlation analysis (Sections 4.1 and 4.2): the correlation
coefficients differ by less than 2% between these maps. It only affects on the
absolute values of the masses, column densities (Section 3.2), and
temperatures (Section 3.4). We present both ViaLactea and 160-500 µm fit
results in the corresponding sections.

###  2.5 Masses and error estimation

We estimate each clump’s mass using the integration of the clump area in the
mass-in-pixel maps. The column density maps from each of the gas and dust
tracers were converted to the maps of the mass-in-pixel using the following
equation:

$\displaystyle
M=N(\mathrm{H_{2}})\mu_{\mathrm{H}_{2}}m_{\mathrm{H}}\mathrm{A}$ (5)

where $\mu_{\mathrm{H}_{2}}=2.8$ is the mean molecular weight of the
interstellar medium (Kauffmann et al., 2008), $m_{H}$ is the mass of an
hydrogen atom, $\rm A$ is the area of all pixels of a clump in cm2.

The errors on masses, column densities and temperatures are given in the
Tables 4-5. They were estimated using uncertainties package of Python and root
mean square (RMS) maps of the following tracers: 12CO(1-0), 13CO(1-0),
Herschel (ViaLactea) 70-500 µm, Herschel 160-500 µm, and AJ extinction. The
dust temperature uncertainty map for ViaLactea (Marsh et al., 2017) is not
available. Thus we do not include the errors in the Herschel 70-500 dust
temperature estimation.

###  2.6 Young stellar objects density maps

Using catalogs from Chavarría et al. (2008), we constructed density maps of
class I and class II YSOs in order to study an evolutionary link between gas,
dust and formation of YSOs. We created the map of the YSO density using the
Kernel Density Estimate (KDE; Scott, 2015) with a Gaussian kernel, implemented
in the KDEpy222https://kdepy.readthedocs.io/en/latest/index.html package of
Python to make a direct comparison between the YSO density map and different
gas and dust tracers. We tested different kernel widths in the range 45-180″.
The kernel width was selected as 60″ as a balanced value between smoothness
and discretization. The resulting YSO densities are shown in Fig. 6. The
Wavelet-based cross-correlation (WWCC) analysis (Arshakian & Ossenkopf, 2016)
provides a more detailed method of correlation analysis between the YSO
density and gas and dust tracers (CO, HCO+, extinction and Herschel maps) on
different spatial scales.

Before executing the correlation analysis, we convolved and regrid the studied
maps to a uniform beam and uniform step size. The beam size of 60″ matches the
YSO density kernel size. After the convolution, the telescope beam has no
impact on the correlation analysis, and one can compare all scales above that
beam limit. The convolution and regridding were done using Python packages
astropy.convolution and reproject, respectively.

Figure 2: Panels 1-5: Maps of hydrogen column density obtained from different
tracers of gas and dust in the S254-S258 region, including two Herschel maps
from ViaLactea database (Marsh et al., 2017) and from this paper, 13CO(1–0)
map (Chavarría et al., 2008), Near-IR extinction map (see Section 2.3 for
details). In the last panel, the HCO+(1–0) intensity map is shown from the
Onsala observations (this paper). All maps except HCO+(1–0) have the same
units and scale of column density starting from 1$\times 10^{20}$ cm-2 and
goes up to 1$\times 10^{23}$ cm-2. For the HCO+(1–0) map the intensity scale
is shown in the lower right corner. In each panel, molecular clumps associated
with star clusters are shown with a white and cyan thick line. The edges of
molecular clumps marked with white borders are ViaLactea column density of
$5.8\times 10^{21}$ cm-2 and cyan border are ViaLactea column density of
8.0$\times 10^{21}$ cm-2. Red circles in each panel represent the visible
radius of H II regions from DSS-R survey.

##  3 Results

###  3.1 Gas and dust clumps

The structure of the S254-S258 star formation region seen by Herschel (see
Figure 2) maps reveals several bright gas and dust condensations, including
S258, G192.63+0.00, S255N, S256-south, G192.75-0.08 and G192.69-0.25.

We investigate only those molecular clumps that are associated with the star
clusters from Chavarría et al. (2008) (see Figure 12 therein). The borders of
clumps were defined based on the ViaLactea column density map using thresholds
of $5.8\times 10^{21}$ cm-2 for S258, G192.63-0.00, G192.75-0.08, G192.69-0.25
and of $8.0\times 10^{21}$ cm-2 for the S255N and S256-south regions. We used
the two different thresholds for clusters to make the area of bright clumps
associated with YSOs clusters (S255N, S256-south) similar to other clumps.
Otherwise, the threshold of $5.8\times 10^{21}$ cm-2 will occupy the
significant fraction of the S254-S258 star formation region and makes it
impossible to distinguish between the gas related to YSOs clusters S255N and
S256-south (Chavarría et al., 2008). Due to the absence of bright dust
emission in the Herschel bands in G192.75-0.00 and G192.54-0.15 clusters, we
used contours of a YSO density at five stars per square parsec from Chavarría
et al. (2008) as the border for these clusters.

In the original work of Chavarría et al. (2008) authors identify two YSOs
clusters in the direction between S255 and S257 H II regions: S255N&S255-2 and
G192.65-0.08. However, our multi-wavelength analysis (see Figure 2) displays
no difference between these clusters in any of the considered gas-dust
tracers. Thus in further analysis we combine these YSOs clusters into a single
cluster named S255N.

The cluster G192.55-0.01 located to the north of the S255N region and contains
only class II YSOs, but displays no emission in any of the available gas and
dust tracers. Moreover, there is no pronounced border between the cluster and
its surrounding class II YSOs. Thus we did not include it to our list.

The clusters S256-south, G192.63-0.00 and S258 are distinguished from the
central region S255N in the Herschel map. Thus we consider them independently.
However, CO and extinction maps display a smooth distribution of gas and dust
toward S258, G192.63-0.00, S255N and S256-south. Thus these clusters may be
considered as different parts of the main gas ridge in the S254-S258 star
formation region.

Clusters G192.75-0.08 and G192.69-0.25 have both class I and II YSOs and have
pronounced peaks of emission in Herschel and CO maps. Thus they were included
in our list.

Cluster G192.75-0.00 contains only class II YSOs and has no pronounced
emission in the Herschel map. However, it has pronounced peaks in the CO and
extinction maps. Thus we include it to our list, but use the YSO density
contour for the mass estimation instead of Herschel column density contours.

Table 3: Parameters of the observed lines at the different points obtained from the Gaussian fitting and physical conditions in the gas. The pointing center have coordinates $\alpha_{\mathrm{J2000}}=6^{\mathrm{h}}13^{\mathrm{m}}4\aas@@fstack{s}1$, $\delta_{\mathrm{J2000}}=17\degr 58\arcmin 44\arcsec$. Values in the brackets are the formal errors from the class program. Point with asterisks mark (*) have offset of 8.7″ in Right Ascension in different lines. $\tau_{0}$ is the optical depth at the line center channel, $\int\tau\mathrm{d}v$ is line-integrated optical depth, $\mathrm{T_{ex}}$ is the excitation temperature and $N$ is molecular column density.  Line | Point | $\Delta\alpha$ | $\Delta\delta$ | $\mathrm{V_{LSR}}$ | $\mathrm{\Delta V}$ | $\mathrm{T_{mb}^{peak}}$ | $\int\mathrm{T_{mb}}$ | $\tau_{0}$ | $\int\tau\mathrm{d}v$ | $\mathrm{T_{ex}}$ | $N$
---|---|---|---|---|---|---|---|---|---|---|---
|  | arcsec | arcsec | km s-1 | km s-1 | K | K km s-1 |  |  | K | 1014 cm-2
CS | 1 | -32.8 | -4 | 7.41 (0.026) | 2.16 (0.063) | 1.76 | 4.06 (0.1) | 1.13 | 0.58 | 5.63 | 0.10
C34S |  | -32.8 | -4 | 7.35 (0.076) | 1.73 (0.202) | 0.13 | 0.24 (0.023) |  |  |  |
HCO+ |  | -32.8 | -4 | 7.48 (0.019) | 2.88 (0.043) | 2.51 | 7.71 (0.1) | 5.24 | 4.54 | 5.49 | 0.52
H13CO+ |  | -32.8 | -4 | 7.22 (0.081) | 2.52 (0.220) | 0.16 | 0.43 (0.03) |  |  |  |
C18O |  | -32.8 | -4 | 6.87 (0.031) | 2.03 (0.071) | 1.12 | 2.43 (0.075) |  |  |  |
HCO+ | 2 | -111.2 | 116 | 8.22 (0.029) | 3.08 (0.071) | 3.25 | 10.68 (0.2) | 5.05 | 6.63 | 6.27 | 1.37
H13CO+ |  | -111.2 | 116 | 8.80 (0.117) | 4.05 (0.277) | 0.2 | 0.85 (0.05) |  |  |  |
CS | 3 | -112.8 | -4 | 6.97 (0.009) | 2.22 (0.020) | 6.11 | 14.45 (0.1) | 3.19 | 2.37 | 9.53 | 1.74
C34S |  | -112.8 | -4 | 7.10 (0.016) | 1.78 (0.037) | 0.86 | 1.63 (0.03) |  |  |  |
HCO+ |  | -112.8 | -4 | 6.98 (0.012) | 2.96 (0.027) | 4.33 | 13.68 (0.1) | 13.01 | 7.99 | 7.37 | 2.62
H13CO+ |  | -112.8 | -4 | 6.67 (0.026) | 1.86 (0.066) | 0.65 | 1.30 (0.04) |  |  |  |
C18O |  | -112.8 | -4 | 6.74 (0.032) | 2.15 (0.078) | 1.96 | 4.49 (0.13) |  |  |  |
CS | 4* | -152.8 | 36 | 7.29 (0.007) | 2.62 (0.017) | 9.98 | 27.85 (0.15) | 5.14 | 3.61 | 13.25 | 6.31
C34S |  | -144.1 | 36 | 7.29 (0.016) | 1.94 (0.041) | 2.09 | 4.33 (0.075) |  |  |  |
HCO+ |  | -152.8 | 36 | 7.63 (0.019) | 3.86 (0.046) | 4.88 | 20.05 (0.2) |  |  |  |
H13CO+ |  | -144.1 | 36 | 7.06 (0.083) | 3.07 (0.180) | 0.46 | 1.49 (0.08) |  |  |  |
CS | 5 | -151.2 | 96 | 8.53 (0.012) | 3.23 (0.029) | 7.36 | 25.28 (0.2) | 4.59 | 3.67 | 10.61 | 4.57
C34S |  | -151.2 | 96 | 8.43 (0.021) | 2.67 (0.049) | 1.4 | 3.99 (0.07) |  |  |  |
HCO+ |  | -151.2 | 96 | 8.86 (0.013) | 3.11 (0.032) | 6.79 | 22.46 (0.2) | 7.53 | 7.77 | 9.88 | 4.18
H13CO+ |  | -151.2 | 96 | 8.76 (0.061) | 3.20 (0.149) | 0.61 | 2.08 (0.08) |  |  |  |
HCO+ | 6 | -191.2 | 116 | 8.79 (0.035) | 2.81 (0.089) | 2.79 | 8.37 (0.2) |  |  |  |
Table 4: Masses and average column densities estimation of molecular clumps in
S254-S258 region from different gas tracers. The value after the $\pm$ symbol
is the $1\sigma$ RMS of the parameter estimate.

Clump mass, M⊙ |
---|---
| Herschel | $A_{\mathrm{J}}$ | CO
Name | ViaLactea | 160-500 µm |  | (1–0)
S255N |  1286 $\pm$ 5.5 |  2813 $\pm$ 13.4 |  1797 $\pm$ 21.2 |  2638 $\pm$ 49.8
S256-south |  1399 $\pm$ 3 |  2158 $\pm$ 8.3 |  2442 $\pm$ 22.5 |  2669 $\pm$ 50.7
G192.63-0.00 |  432 $\pm$ 1.6 |  621 $\pm$ 4.1 |  697 $\pm$ 10.8 |  886 $\pm$ 31
S258 |  218 $\pm$ 1.1 |  399 $\pm$ 3.7 |  486 $\pm$ 9.5 |  524 $\pm$ 24.6
G192.75-0.08 |  339 $\pm$ 1.4 |  603 $\pm$ 6.6 |  732 $\pm$ 12 |  575 $\pm$ 32.9
G192.69-0.25 |  117 $\pm$ 0.8 |  219 $\pm$ 3.4 |  159 $\pm$ 5.3 |  117 $\pm$ 9.3
G192.75-0.00 |  34 $\pm$ 0.5 |  92 $\pm$ 2.5 |  106 $\pm$ 3.9 |  71 $\pm$ 7.5

N, $\times 10^{21}$ cm-2 |
---|---
| Herschel | $A_{\mathrm{J}}$ | CO
Name | ViaLactea | 160-500 µm |  | (1–0)
S255N |  19.34 $\pm$ 0.02 |  39.96 $\pm$ 0.18 |  26.35 $\pm$ 0.2 |  39.63 $\pm$ 0.69
S256-south |  12.12 $\pm$ 0.01 |  17.56 $\pm$ 0.07 |  19.41 $\pm$ 0.1 |  23 $\pm$ 0.42
G192.63-0.00 |  8.13 $\pm$ 0.02 |  11.11 $\pm$ 0.07 |  11.44 $\pm$ 0.2 |  16.64 $\pm$ 0.57
S258 |  8.31 $\pm$ 0.03 |  14.31 $\pm$ 0.13 |  16.35 $\pm$ 0.3 |  19.69 $\pm$ 0.9
G192.75-0.08 |  8.69 $\pm$ 0.02 |  14.5 $\pm$ 0.16 |  18.12 $\pm$ 0.2 |  14.92 $\pm$ 0.84
G192.69-0.25 |  7.23 $\pm$ 0.03 |  12.81 $\pm$ 0.2 |  8.44 $\pm$ 0.3 |  7.13 $\pm$ 0.56
G192.75-0.00 |  3.28 $\pm$ 0.04 |  8.23 $\pm$ 0.22 |  7.96 $\pm$ 0.3 |  6.69 $\pm$ 0.7

Figure 3: Masses (upper panel) and average column densities (lower panel) of
molecular clumps, associated with star clusters in S254-S258 region.

###  3.2 Masses and column densities of the clumps

Figure 4: Dust and CO excitation temperatures of molecular clumps, associated
with star clusters in S254-S258 region. Designations are the same as in Table
5.

We divide all clusters into two main types: main molecular cloud and isolated
clumps. Clusters S258, G192.63-0.00, S255N and S256-south have gas bridges in
CO and the extinction maps. Thus they may be considered as different parts of
the main molecular cloud in the S254-S258 region. Clusters G192.75-0.00,
G192.69-0.25 and G192.75-0.00 are located at some distance to the main gas
ridge and display no gas bridges with the main molecular cloud in any of
available tracers (13CO(1–0), 13CO(2–1), extinction, Herschel). Thus these
clusters may be considered as isolated molecular clumps.

Table 4 and Fig. 3 present masses and column density estimates for the gas
clumps, including errors. We found that, except for S255N, mass estimates from
the CO, extinction and Herschel maps are close to each other. The CO map gives
the highest column-density and mass estimates in more massive clumps (S255N,
S256-south, G192.63-0.00).

The highest value of column density is found towards the S255N clump –
$N_{\mathrm{H_{2}}}=4.0\pm 0.7\times 10^{22}$ cm-2 according to the CO and
$4.0\pm 0.2\times 10^{22}$ cm-2 according to the 160-500 µm Herschel data.
This value have the same order of column density of $\sim 5.3\times 10^{22}$
cm-2 from the C18O spectral line observations of Minier et al. (2005). Mass
estimation for S255N clump is $2638\pm 50$ $M_{\odot}$ according to CO and
$2813\pm 13$ $M_{\odot}$ according to Herschel data, however extinction map
gives only $1797\pm 21$ $M_{\odot}$. Undersampling of the extinction map and
lack of background stars in S255N may cause this discrepancy.

The most massive clumps from all considered here are S255N and S256-south
($2638\pm 49.8$ and $2669\pm 50.7$ $M_{\mathrm{\odot}}$ respectively,
according to 13CO(1–0) data). According to the Herschel map, they contain 71%
of the mass of our selected clumps and 27% of the mass for the whole studied
region (21′$\times$15′) containing all other clumps and inter-clump medium. In
contrast, the area of these clumps is the only 12% of the studied region.

Table 5: Dust and CO excitation temperature estimation of molecular clumps in the S254-S258 region. 70-500 Herschel Tdust values refer to the mean dust temperatures obtained using the 70-500 µm fit from Marsh et al. (2017). 160-500 Herschel Tdust values were obtained from the fit described in Section 2.4. CO excitation temperatures for (1-0) transition were obtained using FCRAO data (Chavarría et al., 2008) and SMT data for (2-1) transition (Bieging et al., 2009). The value after the $\pm$ symbol is the $1\sigma$ error of the parameter estimate.  |  Tdust, K |  Tex, K
---|---|---
| Herschel | CO
Name | 70-500 | 160-500 µm | (1–0) | (2–1)
S255N | 20.11 |  20.18 $\pm$ 0.07 |  26.21 $\pm$ 0.05 |  25.2 $\pm$ 0.03
S256-south | 15.53 |  16.97 $\pm$ 0.02 |  20.6 $\pm$ 0.04 |  18.9 $\pm$ 0.02
G192.63-0.00 | 17.55 |  19.88 $\pm$ 0.04 |  22.47 $\pm$ 0.06 |  20.7 $\pm$ 0.03
S258 | 19.18 |  20.74 $\pm$ 0.07 |  18.88 $\pm$ 0.08 |  17.4 $\pm$ 0.05
G192.75-0.08 | 13.78 |  14.39 $\pm$ 0.04 |  12.43 $\pm$ 0.06 |  11 $\pm$ 0.04
G192.69-0.25 | 13.67 |  14.2 $\pm$ 0.06 |  12.24 $\pm$ 0.11 |  11 $\pm$ 0.06
G192.75-0.00 | 14.93 |  15.3 $\pm$ 0.12 |  13.35 $\pm$ 0.11 |  12.5 $\pm$ 0.07

###  3.3 High-density gas

Table 3 presents the observed parameters of C34S(2–1) and H13CO+(1–0) lines
together with the CS(2–1), HCO+(1–0) in similar offsets to compare the line
intensities. . Results of the column density estimation are given in the Table
3. Analyzing the values of CS column densities, we found them in the range
$N_{\mathrm{CS}}=(0.10-6.31)\times 10^{14}$ cm-2 for the outer and inner parts
of S255N cluster, respectively. Assuming relative abundance of CS to H2 in the
range 10-9-10-8 we estimate the maximum beam-averaged gas column density in
the range $N_{\mathrm{H_{2}}}=0.63-6.31\times 10^{23}$ cm-2. For comparison,
Zinchenko et al. (2009) estimated the gas column densities toward S255N as
$N_{\mathrm{H_{2}}}=2.0\times 10^{23}$ cm-2.

The available high-density tracer observations confirm the detection of
HCO+(1–0) line toward S258, S256, S255, S257 H II regions and G192.75-0.08,
G192.63-0.00, S256-south star cluster. The CS(2-1) line is detected toward
S258, S256, S255 H II regions and G192.75-0.08, G192.63-0.00 star cluster.

We confirm the presence of high-density gas traced by HCO+ line emission in
the inter-clump bridge between the S255N and S256-south star clusters.The
typical intensity of HCO+(1–0) line in the inter-clump bridge is 1.8 K, while
$\sigma$ is 0.3 K, leading to a detection at 6$\sigma$ level. The presence of
high-density gas in the inter-cluster medium suggest that these clusters are
physically associated and share the same gas and dust material. We also
searched for the inter-clump bridge between S255 and S258 H II regions, but we
detected no HCO+ emission at the level $\sigma\sim 0.3~{}$K. That may lead to
the conclusion that these H II regions do not share the same gas and dust
material and maybe are not physically associated. We also note that the 13CO
column density map (see Figure 2) reveals diffuse emission in the inter-clump
medium between the S255 and S258 H II regions. However, this emission may be
due to the large extent of the gas along the line of sight that leads to high
values of column density.

There are other meaningful differences between various maps in Fig. 2. First,
the peak of the S255N region in the HCO+, Herschel maps is shifted to the
north in comparison to the same peak in the 13CO and extinction maps. The
shift may be related to a different distribution of high-density and
intermediate-density gas inside the S255N cluster.

We used the available observations of C18O molecule in two points to check the
optical depth of 13CO line used for mass estimations (see Section 3.2). In
Point 1 (see Table 3) the peak of the C18O(1–0) line is 1.12 K, while the
13CO(1–0) line has a main-beam temperature of 7.2 K. From the 12CO(1–0) and
13CO(1–0) line ratios and equation (1) we obtain $\tau_{0}=30$ for the
12CO(1-0) line. The abundance ratio $R=80$ then gives $\tau_{0}=0.38$ for the
13CO(1–0) line. Using the ratio between 13CO(1–0) and C18O(1–0) lines in the
equation (1) we obtain the ratio of optical depth for 13CO/C18O $R=7.5$. The
same calculations were done for Point 2, where we obtain $R=8.3$. The
resulting values match the 13CO/C18O abundance ratio R$\sim$7.3 found by
Schöier et al. (2002).

###  3.4 Temperature distribution

Figure 5: ViaLactea dust temperature map of S254-S258 region. Contours are
star clusters associated with high-density gas, same as in the Figure 2. Blue
circles are H II regions, same as in the Figure 2.

The dust temperature estimation from the ViaLactea maps gives 17-20 K for
S255N, G192.63-0.00, S258 and 14-16 K for S256-south, G192.75-0.08,
G192.69-0.25 and G192.75-0.00. In the work of Zinchenko et al. (1997) NH3
molecular line observations of S255N region reveal gas kinetic temperatures
23.0 K and 24.7 K at the S255 core peaks. Their positions is (80″,0″) and
(40″,–40″) in respect to the center field coordinates
RA${}_{\mathrm{J2000}}=06^{h}12^{m}53.3^{s}$, Dec${}_{\mathrm{J2000}}=17\degr
59\arcmin 22\arcsec$. These positions are close to point 5 (distance $\approx$
25″) and 6 (distance $\approx$ 8″) from the Table 3. The temperatures drop to
15-20 K at the edges of the ammonia emitting regions. Therefore, gas and dust
temperatures agree with each other in these dense molecular clumps.

The dust temperature ($T_{\rm d}$) obtained from the Herschel data are shown
in Fig. 5 and the mean values are presented in Table 5 and Figure 4 together
with excitation temperatures of the 13CO(1–0) and 13CO(2–1) lines. We found
that $T_{\rm d}$ is lower than the CO excitation temperature (both (2–1) and
(1–0)) in more massive clumps ($>600M_{\odot}$). In contrast, the $T_{\rm d}$
value is higher than the CO excitation temperature (again both (2–1) and
(1–0)) in less massive clumps. We note that the more massive clumps are
associated with the developed H II regions, and the different behaviour can be
explained by the influence of the H II regions in the heating of the
interstellar medium around these regions. Photon-dominated regions (PDR)
around H II regions always show a gas temperature exceeding the dust
temperature close to the UV sources while the dust temperature exceeds the gas
temperature at larger distances (see e.g. Röllig et al., 2013).

We note that the mean the dust temperature of S256-south (15.5 K) and S258
regions (19.2) are significantly different (see Figure 5 and Table 5), while
the mean HCO+(1–0) line intensity is almost similar for both of them – 1.6 K
for S258 and S256-south regions (see Figure 2). The physical interpretation is
that the dust temperature map measures total energy input from heating
sources, while HCO+ emission shows PDR material, i.e., traces the heating
sources’ UV radiation. We conclude that the heating source in the S258 H II
region produces more thermal emission than in the S256 H II region. According
to Russeil et al. (2007), the spectral type of the heating source is B3V for
the S258 and B2.5V for the S256. Although the difference in spectral class is
minor, the difference in the radiation spectrum is significant. The radiation
fields’ difference may be associated with multiple heating sources in the S258
H II region. The YSOs cluster that ignites the S258 H II region may harbor
several stars that do not emit intense UV radiation and heat the surrounding
gas. This radiation may result in a higher dust temperature in the S258 region
than in the S256-south region. However, more observational data are needed on
the S258 and the S256-south region to clarify the situation.

Figure 6: Input maps that were used for correlation analysis. YSOs density maps are named as YSO I and YSO II according to YSOs class and have arbitrary units. Herschel, 13CO(1–0) and extinction maps have units of column density and the same scale and limits, shown at the bottom right corner. Map of HCO+(1–0) line emission have Tmb units. All maps are log-scaled by the base of 100. White contours follow ViaLactea column density thresholds for selected YSOs clusters (same as in the Figure 2). Red circles represent the visible radius of H II regions. Blue crosses in YSOs density map mark the position of individual YSO.  Table 6: Correlation coefficients and power law indices (slopes) between YSOs class I and II density and different tracers of gas and dust convolved to the beam size of 60 arcsec. Values in the brackets are the 95% confidence intervals of the correlation coefficient – no error propagation was used here and the scatter is purely statistical. “Without S255N” region refers to the S254-S258 region not including the central bright S255N cluster (center at J2000 position 6:12:54.2, +17:58:56; size 3.75$\times$5.54 arcmin). The “HCO+ mask” region refers to the region masked with the available coverage map of HCO+ line – see blue polygon in the Figure 6.  gas and dust tracer | Region | Class I YSOs | Class II YSOs
---|---|---|---
|  |  $r$ [$r_{\mathrm{low}}$, $r_{\mathrm{high}}$] | slope | $\chi^{2}$ | $r$ | slope | $\chi^{2}$
Herschel 160-500 micron | Full | 0.75 [ 0.64, 0.83 ] |  1.02 $\pm$ 0.47 | 0.04 | 0.42 [ 0.29, 0.54 ] |  0.31 $\pm$ 0.24 | 0.05
| Without S255N | 0.47 [ 0.25, 0.64 ] |  0.43 $\pm$ 0.67 | 0.02 | 0.52 [ 0.40, 0.63 ] |  0.33 $\pm$ 0.24 | 0.03
|  HCO+ mask | 0.69 [ 0.54, 0.80 ] |  0.97 $\pm$ 0.57 | 0.05 | 0.25 [ 0.07, 0.42 ] |  0.18 $\pm$ 0.36 | 0.03
13CO(1–0) | Full | 0.60 [ 0.44, 0.72 ] |  1.44 $\pm$ 0.49 | 0.19 | 0.41 [ 0.27, 0.53 ] |  0.74 $\pm$ 0.24 | 0.30
| Without S255N | 0.32 [ 0.08, 0.52 ] |  0.87 $\pm$ 0.70 | 0.22 | 0.44 [ 0.30, 0.56 ] |  0.78 $\pm$ 0.25 | 0.29
|  HCO+ mask | 0.52 [ 0.31, 0.67 ] |  1.17 $\pm$ 0.57 | 0.18 | 0.22 [ 0.03, 0.39 ] |  0.42 $\pm$ 0.36 | 0.26
AJ extinction | Full | 0.69 [ 0.56, 0.79 ] |  0.68 $\pm$ 0.47 | 0.03 | 0.57 [ 0.46, 0.66 ] |  0.51 $\pm$ 0.24 | 0.06
| Without S255N | 0.52 [ 0.32, 0.68 ] |  0.57 $\pm$ 0.67 | 0.03 | 0.61 [ 0.50, 0.70 ] |  0.54 $\pm$ 0.24 | 0.05
|  HCO+ mask | 0.64 [ 0.47, 0.76 ] |  0.52 $\pm$ 0.57 | 0.02 | 0.50 [ 0.34, 0.63 ] |  0.29 $\pm$ 0.36 | 0.02
HCO+(1–0) |  HCO+ mask | 0.44 [ 0.22, 0.61 ] |  0.51 $\pm$ 0.57 | 0.05 | -0.62 [-0.72, -0.48 ] |  -0.42 $\pm$ 0.36 | 0.02

##  4 Correlation between gas and YSOs

Fig. 6 presents YSO density map for both class I and class II YSOs, created
using KDE method (see details in Sec. 2.6). From visual inspection of Fig. 6 ,
one can notice that YSOs of both classes (I and II, see panels 1-2 on the Fig.
6) are associated with the regions of high H2 column density (see panels 4-6
on Fig. 6). We used two independent methods to quantify differences in the
distribution of the of class I and II YSOs and their relation to various gas
and dust tracers. The first one is the power law fit to the YSOs density vs
gas column density, similar to one in Gutermuth et al. (2011). The second one
is correlation analysis on the different spatial scales using Fourier
transformation with the WWCC method Arshakian & Ossenkopf (2016). The
following gas and dust tracers were used to study their relationships to YSOs
density: extinction, 13CO(1–0), Herschel 160-500 µm and HCO+(1–0).

###  4.1 A power law fit to the pixel-to-pixel map of YSOs density and gas
column density

The work of Gutermuth et al. (2011) reveals a power law dependence between
YSOs density and gas column density. The power law is found in scatter plots
of YSOs density against gas and dust column densities. The slopes are found to
be 2.67 and 1.87, while the correlation coefficients are 0.87 and 0.83 for
MonR2 and Ophiuchus molecular clouds, respectively. In the work of Bieging et
al. (2016), those authors report a power law dependence of the mean YSO
surface density on total gas surface density with the slope of 1.63. We check
these relations using the available gas and dust tracers of the S254-S258
region. The column density maps obtained in this paper for each tracer were
used to plot the dependence of gas column density to the class I and II YSOs
density in the direction of each YSO (similar to Figure 9 of Gutermuth et al.
(2011)).

Results of the power law fits and direct correlation analysis for class I and
class II YSOs are presented in the Table 6. We convolved each map to the beam
size of 60″ before going to calculations. The beam size was chosen as maximum
beam size from the avaliable maps (extinction map). We also tested other beam
sizes for calculations. The calculations for 90″ beam size resulted in the
relative increase (8-12%) of correlation coefficient ($r$) in comparison to
60″ calculations. For calculation of $r$ we used the log-based values of gas
column density, main-beam temperature and YSOs density. The calculations were
done for several cases:

  1. 1.

We calculated correlation coefficients and power law slopes for the whole
studied region.

  2. 2.

We excluded the brightest part of S255N cluster to study the surrounding gas
around the S255N region.

  3. 3.

We masked all studied maps using available HCO+(1–0) map coverage to compare
the HCO+ map with other gas tracers.

Analysis of the correlation coefficients between class I YSOs and various gas
tracers in a whole studied region (i) reveal a high degree of correlation
($r=0.60-0.75$) with following maps: Herschel, 13CO(1–0) and extinction. The
HCO+ map has a lower value of $r=0.44$. However, the lower value of $r$ may be
related to the different unit scale of HCO+ map – main-beam temperature
instead of column density for other gas tracers. Using the direct values of
main beam temperature instead of log-based values increase the value of $r$
from $0.44_{-0.22}^{+0.17}$ to $0.62_{-0.17}^{+0.13}$ for HCO+ map.

Similar to Gutermuth et al. (2011), we found a power law dependence between
YSOs density and column density of various gas tracers. The average slope of
power law between YSOs of class I and various gas tracers is 0.9, while slope
values range from 0.51 (HCO+ map) to 1.44 (CO map). The found values are
generally lower than those found for MonR2 (2.67) and Ophiuchus clouds (1.87)
in Gutermuth et al. (2011) and the slope of 1.63 for the S235 region found in
Bieging et al. (2016).

We further studied the behaviour of correlation coefficients on masked
regions. Exclusion of the central bright S255N region (ii) resulted in a
change of correlation coefficients in Herschel, and 13CO(1–0) maps – the
values are decreased from 0.6-0.7 to 0.3-0.4. We conclude that a high degree
of correlation between those maps and class I YSOs is associated with the
bright S255N region, rich in class I YSOs. Other regions in those maps display
relatively high values of $r=0.5-0.6$. However, for the extinction map, the
exclusion of the S255N region does not affect of correlation coefficient
significantly - the value of $r$ is decreased from 0.69 to 0.52. We associate
this with S256-south, S258 and G192.75-0.08 regions that all have high column
densities ($N(\mathrm{H_{2}})\sim 14-19\times 10^{21}$ cm-2, see Table 4) and
rich in class I YSOs (see Figure 6).

We note that HCO+ map has a limited coverage, thus the results for HCO+ would
be much stronger and are subject to change if the mapping of the whole region
becomes available. In order to compare the correlation coefficients of HCO+
map with other maps we have to apply the same HCO+ coverage mask (iii) to all
other maps. The region with available HCO+ data is marked as “HCO+ mask” in
the Table 6. The masking resulted in decrease of correlation coefficient for
all studied maps from 0.6-0.75 to 0.52-0.69. Thus we conclude that correlation
coefficient between HCO+ map and class I YSOs density map is decreased due to
limited coverage of the HCO+ map. This is also confirmed in the spatial
correlation analysis in the “HCO-mask” region (see Section 4.2), where values
of $r$ are almost similar for all gas and dust tracers within the “HCO-mask”
region.

The correlation coefficients and power law slopes between class II YSOs and
different gas and dust tracers are lower than for the class I YSOs. The values
of $r$ range from 0.41 (CO map) to 0.57 (extinction map). The slope for power
law is shallower than the similar slope for class I YSOs and have a mean value
of 0.52 for those tracers that have $r>0$. The exclusion of the S255N region
results in an increase in correlation coefficients from 0.29-0.57 to
0.38-0.61. We associate the increase with the small amount of class II YSOs in
the S255N cluster (see Figure 6), but high values of column density toward
this cluster (see Table 4), the maximum for the whole S254-S258 region.
Excluding such a region removes the main difference between maps of YSOs class
II and gas tracers.

Another notable feature is the negative value of $r=-0.62$ for HCO+ map and
class II YSOs density. If one uses the direct values of HCO+(1–0) main-beam
temperature instead of log-based values, this results in
$r=-0.50^{+0.15}_{-0.13}$. Thus the HCO+ map reveals a statistically
significant ($p<0.001$) negative value of correlation coefficient
$r=-0.58^{+0.14}_{-0.12}$ with class II YSOs density. Visual inspection of
HCO+ and class II YSOs density maps confirms the found negative relation: the
peaks of HCO+ emission do not coincide with maxima in class II YSO density
(see Figure 6).

In the S255N region the peak of YSOs class II density has an offset of $\sim$3
arcmin to the north and to the south from the HCO+ emission peak. In the
S256-south region the peak of class II YSOs density have an offset of
$\sim$2.5 arcmin to the east from the HCO+ emission peak. Only the S258 region
displays no offset between YSOs class II and HCO+ emission peak. The
G192.75-0.08 cluster is covered by line observations but is not associated
with strong HCO+ emission. An average main-beam temperature of HCO+ line in
G192.75-0.08 is 0.7 K, while in the S258 and S255N regions it has values of
$\sim 2$ K and $\sim 5$ K, respectively.

We conclude that class II YSOs density is anti-correlated to the high-density
gas traced by HCO+ map. At the same time, the class II YSOs density map
reveals a weak positive correlation with other tracers of gas and dust
including extinction, 13CO(1–0) and Herschel maps with the value of
$r=0.41-0.57$. In contrast, class I YSOs have much more pronounced positive
correlation with all studied maps with the values of $r$ range from 0.60 to
0.75.

###  4.2 Correlation on different spatial lengths

Direct comparison between YSOs density and gas column density does not give
information about the spatial scale at which correlation becomes more
pronounced. To reveal this spatial scale and measure the dependence of
correlation coefficients on spatial scale quantitatively, we use the WWCC
method (Arshakian & Ossenkopf, 2016). This method gives the possibility to
find correlation coefficients between two images on different scales using
Fourier decomposition of input maps to the different scales and subsequent
calculations of correlation coefficients on different scales. For
calculations, we use three cases, similar to those described in the Section
4.1: the whole S254-258 region, the region without the S255N cluster and the
region masked with the HCO+ map coverage. Excluding the S255N region allows
one to study the surrounding gas around the main cluster. Masking all tracers
by the HCO+ map coverage allow us to compare the HCO+(1–0) map with all other
tracers.

To compare the different scales of YSOs density with gas and dust tracers, we
used the YSO density map with a gaussian kernel width of 60 arcsec (see
details in the Section 2.6). In order to distinguish between actual
correlations and the artifact of comparing the mapped structures with the YSO
density maps with increasingly large kernels we run the “null” test. We create
the random distribution of YSOs in the studied region with the same number of
sources as actual YSO catalog and create the 60 arcsec density map to test the
behaviour of spatial correlation analysis. One can then compare the test
density map with the extinction, CO and Herschel maps. The results of the
spatial correlation analysis together with the “null” test using the WWCC
method are presented in the Figure 7.

(a) Class I YSOs (Whole region) (b) Class II YSOs (Whole region)

(c) Class I YSOs (Without S255N) (d) Class II YSOs (Without S255N)

(e) Class I YSOs (HCO+ mask) (f) Class II YSOs (HCO+ mask)

Figure 7: Result of spatial correlation analysis between YSOs density and gas
and dust tracers for the whole region (panel A,B), for the whole region
without S255N cluster (panel C,D), and the HCO+ map masked region (panel E,F).
The panels A, C, D display correlation coefficients with the class I YSOs,
panels B, D, F display correlation coefficients with class II YSOs. The
vertical red line in each plot displays all studied maps’ beam size (60
arcsec). The gray-filled curve outline the error bars. Values of $p$ are
results of the two-dimensional Kolmogorov-Smirnov (KS) test between WWCC
analysis of the actual YSOs distribution and $N$ randomly distributed sources.
One chose the lowest value of $p$ for all KS-tests with the different number
of random sources ($N=50,100,150,250,500$).

####  4.2.1 Class I YSOs

In the panel A of the Figure 7 we present the correlation between class I YSOs
density and different gas and dust tracers for the whole studied region. At
small scales (up to $\sim$2 arcmin) CO and Herschel maps tracers reveal
similar values of $r=0.5-0.6$. The extinction map at small scales have smaller
values of $r$. The lower correlation between YSO density and extinction at
small scales may be related to the larger relative noise of the extinction
measurements. As we convolve other maps to the 60 arcsec resolution of the
extinction map, they have a much smoother appearance. The uncertainty of the
extinction measurements also increases towards high column densities, as the
higher cloud opacity decreases the density of observed stars. This could also
bias the extinction measurements, as stars are preferentially detected in
those parts of the averaging beam where the column density values are the
lowest (Juvela & Montillaud, 2016b). The other reason may be the shift of the
extinction peak in the S255N region to the south compared to a similar peak in
the CO and Herschel maps. The class I YSOs density map has a peak position
similar to CO and Herschel maps, but the extinction map reveals the peak
shifted to the south. For larger scales (> 2.0 arcmin) all maps have almost
similar values of $r$ within the error bars.

In the panels C and E, the correlation analysis within the masked regions is
presented. In the “no S255N region” (panel C) the values of $r$ for extinction
map is increased in comparison to CO and Herschel maps. Thus for the gas
surrounding the S255N region the extinction map better represent the
distribution of class I YSOs in comparsion to CO and Herschel maps. The direct
pixel-by-pixel correlation analysis from the Section 4 reveal a similar
conclusion – the extinction map in the “no S255N region” have the largest
value of $r=0.52$ among other tracers. From the correlation analysis we found
that this difference related to scales from 1 to 4 arcmin.

From the analysis of Figure 7 (panels A,C) we conclude that for class I YSOs
there is the ‘saturation point’ of $\sim 2-2.5$ arcmin. After reaching of this
scale the correlation coefficients between column density and class I YSOs
density do not increase significantly. We discuss the implication of this
scale in the Section 5. The scale of $\sim 2-2.5$ arcmin is confirmed in all
three studied cases: full region (panel A), excluding the S255N region (panel
C) and the “HCO+ mask” region (panel E). In contrast, for class II YSOs the
distribution of correlation coefficients is much smoother and contains no
saturation point (see panels B,D of Figure 7). Thus it is not obvious at which
spatial scale the correlation becomes more pronounced.

Comparing the correlation coefficients between class I and II YSOs density and
gas-dust tracers (panels A-D) we note that the extinction map has higher
correlation coefficients with YSOs density of both class I and II at large
scales (> 1.5-2 arcmin) in comparison to all other tracers. That is especially
pronounced in the panels B, C and D of the Figure 7. Comparing the YSOs class
I and II density and extinction map visually (see Figure 6), we found that all
regions that have peaks of class I and II YSOs density (S255N, S256-south,
G192.75-0.08, S258) similarly have the peak of extinction. Thus we conclude
that extinction maps is a good tracer of both less evolved class I YSOs and
more evolved class II YSOs.

The panels E and F of the Figure 7 present the result of the spatial
correlation analysis under the “HCO+ mask” region. From the analysis of the
panel E we conclude that HCO+(1–0) map have values of $r$ with class I YSOs
density similar to other tracers of gas and dust within the statistical error.
This may give an idea of distribution of HCO+(1–0) line emission in the whole
unmasked region – it should generally follow other gas and dust tracers,
especially at large scales (> 1.5 arcmin). We also note that at very large
scales (4-8 arcmin) values of $r$ are decreased from 0.8 to 0.7. This decrease
may be associated with the influence of class I YSOs that are outside of the
HCO+ mask coverage. Much stronger decrease is pronounced between class II YSOs
and various gas and dust tracers (see panel f of Figure 7). The number of
class II YSOs that fall outside of HCO+ mask coverage region is very large,
thus we may see even negative values of $r$.

####  4.2.2 Class II YSOs

The correlation analysis between Class II YSOs and gas-dust tracers is
presented in Panels B and D of Figure 6. The correlation of Class II YSOs is
much less pronounced with all studied gas and dust tracers when comparing to
the Class I YSOs (panels A, C). Even from visual inspection of Figure 6, one
can see that the more evolved Class II YSOs are located at some distance from
the regions of high gas column density in S255N and S256-south regions.

If we consider the full studied region (panel B), the best correlation with
class II YSOs is found in extinction map. The situation is not changed
significantly if we exclude the S255N region and study the gas surrounding
S255N (panel d) – extinction maps have the highest $r$ compared to other maps.
However, the absolute values of $r$ for class II YSOs are lower than for class
I YSOs. We note that for the CO map, the values of $r>0.5$ for class II YSOs
become pronounced only at quite large scales (> 4 arcmins) and only when
considering the mask excluding the S255N region (panel D). For the extinction
map, the value of $r>0.5$ reached at the smaller scale of 2 arcmins when
considering the mask excluding the S255N region (panel D). As the S255N region
contain a small number of class II YSOs, but column density is significant,
excluding that region results in much higher values of $r$ for all tracers.
For class I YSOs values of $r=0.5-0.6$ are pronounced at much smaller scales
(see panel C) starting from 1.5 arcmins in the unmasked region. We conclude
that class II YSOs are located at a significant distance from the high column
density regions and have only a weak relation to the maps considered here.

Correlation analysis of various gas tracers with class II YSOs in HCO+ mask
(panel F of the Figure 7) reveal very low and even negative values of $r$.
This may occur because significant part of class II YSOs reside in the regions
outside the HCO+ mask region. Thus we conclude that it is difficult to study
the spatial correlation between class II YSOs and gas and dust tracers in the
HCO+(1–0) masked region. However, direct pixel-by-pixel analysis may be used
in this case, as it has no influence of class II YSOs that are outside of the
mask.

From panels B and D, it is notable that increasing the spatial scale between
class II YSOs and gas-dust tracers results in continuously increasing
correlation coefficients. We associate this increase with the class II YSOs’
proper motions – during their evolution, the class II YSOs have enough time to
fly away from the high-density region. Thus we can see the increase in
correlation coefficients when increasing the spatial scale. However, increased
correlation coefficients at large scales may be associated with map similarity
due to convolution to large kernel size. In the next section, we will use the
test maps to check the difference between actual YSOs and random sources
distribution while comparing them with gas-dust tracer using the WWCC method.

####  4.2.3 Comparsion with test maps

To check the plausibility of comparing actual YSOs density and gas-dust
tracers, we run several tests with the randomly distributed sources. These
maps were analyzed using the WWCC method together with the gas-dust tracers
(see Appendix B). The two-dimensional Kolmogorov-Smirnov (KS) test was used to
check the difference between the actual YSOs correlation with gas-dust tracers
and correlations with the $N$ random sources density map. The different values
of $N$ have been investigated: 50, 100, 150, 250, and 500 sources.

The resulting $p$-values are presented in the legend of Figure 7 and results
of WWCC analysis using the random sources maps presented in Figure 9 of
Appendix B. The KS test reveals a statistically significant difference between
correlation coefficients of gas-dust tracers with class I/II YSOs and randomly
distributed sources. While using the extinction map for WWCC analysis, the
difference between class I and class II YSOs and random sources distribution
is statistically significant in full and “no S255N” regions, i.e., the null
hypothesis about the same distribution can be rejected at a 3$\sigma$ level (p
< 0.001). In the “HCO+ mask” region, the difference between Class I YSOs and
random is also significant at $3\sigma$ level (p < 0.001). The exceptional
case is class II YSOs and map with $N=$150 sources in the “HCO+” mask region –
the KS test gives a $p$-value of 0.0018. Thus one cannot reject a null
hypothesis about the same distribution at a 3$\sigma$ level ($p$ > 0.001).

The same behavior is found for the CO and Herschel maps: both class I and
class II YSOs density maps are statistically different from random sources
distribution except for the “HCO+” mask region, where class II YSOs density
map correlations is similar to the random map of 100 sources ($p$ = 0.14). We
conclude that both class I and class II YSOs display a statistically
significant correlation with gas-dust tracers, not related to the effect of
comparing large-scale structures. However, the “HCO+” mask region has
insufficient spatial coverage to display statistically significant differences
between actual YSOs position and random sources.

Figure 8: Distribution of class I (left panel) and class II (right panel) YSOs
in the S254-S258 star formation region. YSOs are marked with blue crosses. The
black ellipse in the left panel represents the relative span of class I YSOs
in the S255N cluster, while the red ellipse on the right panel represents the
span of class II YSOs using the expansion velocity of 1.8 km s-1 found from
the dispersion of HCO+ line. The black ellipse in the right panel displays
class II YSOs’ span using an expansion velocity of 1.3 km/s found in Bieging
et al. (2016) for the S235-AB region.

##  5 History of the star formation in the S255 region

From the direct comparison between gas and dust column densities and YSOs
density (see Section 4.1), it was shown that the density of class I YSOs has a
positive power law relation to the column density with an average slope of
$\sim 1.2$ and a correlation coefficient of $\sim 0.6-0.7$ for all studied gas
and dust tracers. In comparison, the slope for class II YSOs is $\sim 0.65$
and the correlation coefficient is only $\sim 0.3-0.4$. This difference may be
related to the strong association between young stars and the primordial gas
where these stars were born, traced by the different gas and dust tracers.

More evolved class II YSOs display less pronounced power law relations and
lower values of correlation coefficients due to the greater age of class II
stars, as they have managed to move away from the primordial gas. At small
scales (< 4 arcmin) the correlation between class II YSOs and column-density
tracers is poor ($r=0.2-0.4$) because the stars have already left their natal
clumps, but at larger scales (> 4 arcmin) some correlation remains
($r=0.5-0.6$) because the YSOs have not yet moved far enough.

From the spatial correlation analysis, we found that at $\sim 2-2.5$ arcmin,
there is the ‘saturation point’ where correlation coefficients between class I
YSOs and various gas and dust tracers began to grow much more slowly. At the
assumed distance of 1.59 kpc, this spatial scale corresponds to $\sim$0.9-1.1
pc. From the typical duration of class I YSOs stage (0.5 Myr according to
Heiderman & Evans (2015)), the average YSOs travel speed is then $\sim$1.8-2.3
km s-1. The work of Bieging et al. (2016) argue that the spatial distribution
of YSOs of Class 0/I compared with Class II is consistent with the young stars
having the same velocity dispersion as that which we observe in the gas for
the highest column density regions. From the 13CO(2-1) line in the direction
of the S255N region, we found the median line dispersion of 1.4 km s-1 with
the maximum value of 1.86 km s-1. A similar analysis for HCO+(1–0) line give
median dispersion of 1.3 km s-1 and maximum dispersion of 1.7 km s-1. Thus the
travel speed of $\sim$1.8 km s-1 is consistent with the 13CO(2-1) velocity
dispersion in the direction of high column density regions. However, the
actual YSO movements may be slower than the gas dispersion and the reaching of
the spatial scale turning point. Thus we assume 1.8 km s-1 as an upper limit.

We analyzed YSOs evolutionary stages duration in S255 region similar to
Bieging et al. (2016). We enclose the Class I YSOs in the S255N region in the
ellipse with the major and minor axes of 2.97$\times$1.39′. At the distance of
1.59 kpc, this corresponds to 1.38$\times$0.64 pc. Using a travel speed of 1.8
km s-1 and duration of class II YSO stage of 2 Myr (Heiderman & Evans, 2015)
we conclude that due to random YSO motions, the ellipse of class II YSOs
should be increased by 3.67 pc in both directions that lead to dimensions of
4.32$\times$5.05 pc. Figure 8 reveal the initial ellipse on the left panel and
expanded ellipse using the expansion velocity of 1.8 km s-1 in the right
panel. We additionally put the ellipse for class II YSO using the expansion
speed of 1.3 km s-1, similar to the one used in Bieging et al. (2016). From
the analysis, we found that the almost all YSOs clusters observed in the
S254-S258 region are enclosed in the expanded ellipse of 1.8 km s-1. The
situation is not changed significantly if one used the speed of 1.3 km s-1, as
found by Bieging et al. (2016) for S235-AB region.

In the classical picture of star formation, the gas is concentrated into a
single region, stars forming in situ, and the resulting cluster then
dispersing. In the S254-S258 star-formation region, the absolute peak of class
I YSOs density was found toward the S255N cluster, which also has the absolute
maximum of gas column density. Thus we can assume that this cluster may be the
starting point for many YSOs that already have left its birthplace. Starting
from the S255N, YSOs have enough time to move through the whole S254-S258
region. Given that the whole star formation region contains several isolated
YSOs clusters, the evolutionary link between the S255N cluster and the other
YSOs clusters can be considered. A possible example of YSOs movement is seen
in the S255N and S255N2 clusters. As we observe many class I YSOs and lack of
class II YSOs in the S255N region, we conclude that YSOs are quickly flying
away from the S255N region. However, YSOs are dispersed not in an isotropic
way – most of them have moved toward the S255N2 region.

In the S255N2 region many of class II YSOs are observed, but gas and dust
emission was not detected in all considered tracers, including CO, 13CO, CS,
HCO+, near-IR extinction and Herschel maps. The class II YSOs is much more
scattered in S255N2 region in comparison to other regions rich in class II
YSOs (S256-south, S258, G192.75-0.08, G192.75-0.0). Given that these regions
harbour gas with high column density we suggest that interaction between gas
and stars can create better conditions for stability of star clusters in
comparison to cluster that is just moved away from it’s birth material
(S255N2) – clusters formed in deeper potential wells should disperse more
slowly. However, the stability of YSOs cluster is not only about the gas mass
but also about the mass in the cluster itself. The differences in the cluster
structure could also be related to different dynamics during their formation.

However, with the available data, we can’t resolve why Class II YSOs are
moving primarily to the north of the S255IR region and how other YSOs clusters
(S256-south, G192.75-0.08, G192.75-0.00, S258, G192.63-0.00) are formed. It is
not possible to determine whether star formation in the gas clumps is going
continuously or episodically. Ultimately, Class II YSOs formed earlier in the
same gas condensation and moved away while class I sources are formed more
recently and stay close to positions of their formation. Moreover, the
duration of class II YSOs stage ($\sim$2 Myr), and the gas velocity dispersion
of 1.8 km s-1 allow one to travel through the whole S254-S258 region.

The HCO+ map reveals the regions of primordial gas where stars are being born.
After young stars evolve they are moved to the class I and we can see
relatively high values of $r=0.6-0.8$ between class I YSOs density and gas
column density traced by extinction, CO and HCO+ maps at scales 1.5-3 arcmin.
When YSOs enter the class II stage, they have already moved away from from the
primordial gas, and we can find a good correlation ($r=0.4-0.6$) between class
II YSOs and gas column density only at large scales (>4 arcmin), while small
scales (< 4 arcmin) have low values of $r=0.1-0.4$. Thus the CO and extinction
maps at large scales resemble the "old" star formation history seen in class
II objects. In contrast, HCO+, CO and extinction map at smaller scales (<1.5
arcmin) resembles the recent star formation history seen in class I YSOs.

Concluding the view on the star formation in the S254-S258 region, the
following picture comes out. During the evolution the class II YSOs have
enough time to move thought the whole S254-S258 region. Thus the clusters of
Class II YSOs in the S254-258 star-formation region can contain objects born
in the different locations of the complex. The younger class I YSOs are found
near it’s primordial gas, thus we observe high correlations between class I
and high column density regions. The more evolved class II YSOs are not
directly related to gas column density at small scales. However, we found an
continuously increasing correlation coefficient ($r=0.4-0.8$) between gas
column density and class II YSOs density at large scales (4-8 arcmins). We
associate this increase with the proper motions of class II YSOs – during its
evolution ($\sim$2 Myr) the class II YSOs have enough time to fly away to the
maximum distance of 3.67 pc that corresponds to $\sim 8$ arcmin at the
distance of 1.59 kpc.

##  6 Conclusions

We have studied the star-forming region S254-S258, which includes several H II
regions and star clusters, in the different density tracers of molecular gas
(CO, 13CO, HCO+ and CS line emission) and dust (near-IR extinction and
Herschel maps). The following conclusions are obtained:

  1. 1.

High-density gas was detected in the inter-clump bridge between S255N and
S256-south regions. Thus these regions may have an evolutionary link and trace
different parts of the same molecular cloud.

  2. 2.

Direct pixel-to-pixel correlation analysis revealed that less evolved class I
YSOs are well-correlated with all studied gas and dust tracers with a
correlation coefficient of $r=0.6-0.7$. On the contrary, class II YSOs have a
larger spread and the lower correlation coefficients with the gas and dust
tracers: $r=0.4-0.6$.

  3. 3.

From the spatial correlation analysis it was found that the spatial size at
which the correlations coefficients between gas column density and YSOs reach
maximum values is $\sim$2 arcmins for class I YSOs. That corresponds to the
scale of $\sim 1.4$ pc and mean YSOs velocity of $\sim$1.8 km s-1 while
assuming the distance of 1.59 kpc and duration of class I YSO stage of 0.5
Myr.

  4. 4.

For class II YSOs, the spatial correlation coefficient increase continuously
with the increase of the scale without a saturation point. We associate this
increase with proper motions of class II YSOs. According to the Kolmogorov-
Smirnov test, this increase is statistically different from the correlation
between gas-dust tracers and a random sources distribution.

  5. 5.

In the direction of S255N region the velocity dispersion in the high column
density regions is 1.8 km s-1 according to 13CO(2-1) line emission and 1.7 km
s-1 according to HCO+ line emission. This velocity is consistent with the
velocity found from the spatial correlation analysis.

  6. 6.

Assuming the motion of YSOs with the speed of 1.8 km s-1, starting from the
central S255N reigon, class II YSOs have enough time ($\sim$2 Myr) to reach
the whole S254-S258 region, including isolated YSOs clusters.

  7. 7.

We found a non-isotropic distribution of class II YSOs around S255N that is
shifted to the north and absence of class II YSO toward the S255N cluster.

  8. 8.

The spatial correlation analysis lead to the following conclusions: the HCO+
map reveals the regions of primordial gas where stars are being born. Class I
YSOs have relatively high values of correlation coefficient $r=0.6-0.7$ with
HCO+, CO and extinction map at smaller scales (<2.0 arcmin), thus they
resembles the recent star formation history. In contrast, the CO and
extinction maps at large scales (> 4 arcmins) resemble the "old" star
formation history seen in class II objects.

## Acknowledgment

The work of M. S. K in the Section 2.1 was supported by the RFBR grant
20-02-00643. The work by D. A. L. in Section 3 was supported by the Russian
Ministry of Science and Higher Education, FEUZ-2020-0030. The work of S. A. K.
of correlation analysis in Section 4 was supported by the Russian Science
Foundation grant 19-72-10012. The work of A. M. S in Section 5 was supported
by the Large Scientific Project of the Russian Ministry of Science and Higher
Education “Theoretical and experimental studies of the formation and evolution
of extrasolar planetary systems and characteristics of exoplanets” (No.
075-15-2020-780, contract 780-10). V. O. was supported by the Collaborative
Research Centre 956, sub-project C1, funded by the Deutsche
Forschungsgemeinschaft (DFG), project ID 184018867. The authors acknowledges
support from Onsala Space Observatory for the provisioning of its
facilities/observational support. The Onsala Space Observatory national
research infrastructure is funded through Swedish Research Council grant No
2017-00648.

This work use data products from UKIRT Infrared Deep Sky Survey (UKIDSS)
archive. The UKIDSS project is defined in Lawrence et al. (2007). UKIDSS uses
the UKIRT Wide Field Camera (Casali et al., 2007). The photometric system is
described in Hewett et al. (2006), and the calibration is described in Hodgkin
et al. (2009). The pipeline processing and science archive are described in
Hambly et al. (2008).

## Availability of data

The data underlying this article are available in Figshare, at
https://dx.doi.org/10.6084/m9.figshare.c.5325083

## References

  * Abreu-Vicente et al. (2017) Abreu-Vicente J., Stutz A., Henning T., Keto E., Ballesteros-Paredes J., Robitaille T., 2017, A&A, 604, A65
  * Arshakian & Ossenkopf (2016) Arshakian T. G., Ossenkopf V., 2016, A&A, 585, A98
  * Bernes (1979) Bernes C., 1979, A&A, 73, 67
  * Bieging et al. (2007) Bieging J. H., Peters W. L., Vilaro B. V., Schlottman K., Kulesa C., 2007, in Elmegreen B. G., Palous J., eds, IAU Symposium Vol. 237, Triggered Star Formation in a Turbulent ISM. pp 396–396, doi:10.1017/S1743921307001846
  * Bieging et al. (2009) Bieging J. H., Peters W. L., Vila Vilaro B., Schlottman K., Kulesa C., 2009, AJ, 138, 975
  * Bieging et al. (2016) Bieging J. H., Patel S., Peters W. L., Toth L. V., Marton G., Zahorecz S., 2016, ApJS, 226, 13
  * Blake et al. (1987) Blake G. A., Sutton E. C., Masson C. R., Phillips T. G., 1987, ApJ, 315, 621
  * Bohlin et al. (1978) Bohlin R. C., Savage B. D., Drake J. F., 1978, ApJ, 224, 132
  * Buslaeva et al. (2021) Buslaeva A. I., Kirsanova M. S., Punanova A. F., 2021, Ethynyl around the HII regions S255 and S257 (arXiv:2104.07355)
  * Cardelli et al. (1989) Cardelli J. A., Clayton G. C., Mathis J. S., 1989, ApJ, 345, 245
  * Casali et al. (2007) Casali M., Adamson A., Alves de Oliveira C., 2007, A&A, 467, 777
  * Chavarría et al. (2008) Chavarría L. A., Allen L. E., Hora J. L., Brunt C. M., Fazio G. G., 2008, ApJ, 682, 445
  * Cyganowski et al. (2007) Cyganowski C. J., Brogan C. L., Hunter T. R., 2007, AJ, 134, 346
  * Diplas & Savage (1994) Diplas A., Savage B. D., 1994, ApJS, 93, 211
  * Gaches et al. (2015) Gaches B. A. L., Offner S. S. R., Rosolowsky E. W., Bisbas T. G., 2015, ApJ, 799, 235
  * Gudennavar et al. (2012) Gudennavar S. B., Bubbly S. G., Preethi K., Murthy J., 2012, ApJS, 199, 8
  * Gutermuth et al. (2011) Gutermuth R. A., Pipher J. L., Megeath S. T., Myers P. C., Allen L. E., Allen T. S., 2011, ApJ, 739, 84
  * Güver & Özel (2009) Güver T., Özel F., 2009, MNRAS, 400, 2050
  * Hambly et al. (2008) Hambly N. C., et al., 2008, MNRAS, 384, 637
  * Heiderman & Evans (2015) Heiderman A., Evans Neal J. I., 2015, ApJ, 806, 231
  * Hewett et al. (2006) Hewett P. C., Warren S. J., Leggett S. K., Hodgkin S. T., 2006, MNRAS, 367, 454
  * Hodgkin et al. (2009) Hodgkin S. T., Irwin M. J., Hewett P. C., Warren S. J., 2009, MNRAS, 394, 675
  * Juvela & Montillaud (2016a) Juvela M., Montillaud J., 2016a, A&A, 585, A78
  * Juvela & Montillaud (2016b) Juvela M., Montillaud J., 2016b, A&A, 585, A78
  * Kauffmann et al. (2008) Kauffmann J., Bertoldi F., Bourke T. L., Evans II N. J., Lee C. W., 2008, A&A, 487, 993
  * Lawrence et al. (2007) Lawrence A., Warren S. J., Almaini O., 2007, MNRAS, 379, 1599
  * Lombardi (2009) Lombardi M., 2009, A&A, 493, 735
  * Lombardi et al. (2014) Lombardi M., Bouy H., Alves J., Lada C. J., 2014, A&A, 566, A45
  * Mangum & Shirley (2015) Mangum J. G., Shirley Y. L., 2015, PASP, 127, 266
  * Marsh et al. (2017) Marsh K. A., et al., 2017, MNRAS, 471, 2730
  * Milam et al. (2005) Milam S. N., Savage C., Brewster M. A., Ziurys L. M., Wyckoff S., 2005, ApJ, 634, 1126
  * Minier et al. (2005) Minier V., Burton M. G., Hill T., Pestalozzi M. R., Purcell C. R., Garay G., Walsh A. J., Longmore S., 2005, A&A, 429, 945
  * Minier et al. (2007) Minier V., Peretto N., Longmore S. N., Burton M. G., Cesaroni R., Goddi C., Pestalozzi M. R., André P., 2007, in Elmegreen B. G., Palous J., eds, IAU Symposium Vol. 237, Triggered Star Formation in a Turbulent ISM. pp 160–164, doi:10.1017/S1743921307001391
  * Molinari et al. (2016) Molinari S., et al., 2016, A&A, 591, A149
  * Müller et al. (2001) Müller H. S. P., Thorwirth S., Roth D. A., Winnewisser G., 2001, A&A, 370, L49
  * Ojha et al. (2011) Ojha D. K., et al., 2011, ApJ, 738, 156
  * Predehl & Schmitt (1995) Predehl P., Schmitt J. H. M. M., 1995, A&A, 500, 459
  * Röllig et al. (2013) Röllig M., Szczerba R., Ossenkopf V., Glück C., 2013, A&A, 549, A85
  * Russeil et al. (2007) Russeil D., Adami C., Georgelin Y. M., 2007, A&A, 470, 161
  * Rygl et al. (2010) Rygl K. L. J., Brunthaler A., Reid M. J., Menten K. M., van Langevelde H. J., Xu Y., 2010, A&A, 511, A2
  * Sadavoy et al. (2018) Sadavoy S. I., et al., 2018, ApJ, 852, 102
  * Sault et al. (1995) Sault R. J., Teuben P. J., Wright M. C. H., 1995, in Shaw R. A., Payne H. E., Hayes J. J. E., eds, Astronomical Society of the Pacific Conference Series Vol. 77, Astronomical Data Analysis Software and Systems IV. p. 433
  * Schöier et al. (2002) Schöier F. L., Ryde N., Olofsson H., 2002, A&A, 391, 577
  * Scott (2015) Scott D. W., 2015, Multivariate Density Estimation: Theory, Practice, and Visualization
  * Shirley et al. (2013) Shirley Y. L., et al., 2013, ApJS, 209, 2
  * Skrutskie et al. (2006) Skrutskie M. F., Cutri R. M., Stiening R., 2006, AJ, 131, 1163
  * Taylor (2005) Taylor M. B., 2005, in Shopbell P., Britton M., Ebert R., eds, Astronomical Society of the Pacific Conference Series Vol. 347, Astronomical Data Analysis Software and Systems XIV. p. 29
  * Wang et al. (2011) Wang Y., et al., 2011, A&A, 527, A32
  * Zemlyanukha et al. (2018) Zemlyanukha P. M., Zinchenko I. I., Salii S. V., Ryabukhina O. L., Liu S.-Y., 2018, Astronomy Reports, 62, 326
  * Zhu et al. (2017) Zhu H., Tian W., Li A., Zhang M., 2017, MNRAS, 471, 3494
  * Zinchenko et al. (1994) Zinchenko I., Forsstroem V., Lapinov A., Mattila K., 1994, A&A, 288, 601
  * Zinchenko et al. (1997) Zinchenko I., Henning T., Schreyer K., 1997, A&AS, 124, 385
  * Zinchenko et al. (2009) Zinchenko I., Caselli P., Pirogov L., 2009, MNRAS, 395, 2234
  * Zinchenko et al. (2012) Zinchenko I., Liu S. Y., Su Y. N., Kurtz S., Ojha D. K., Samal M. R., Ghosh S. K., 2012, ApJ, 755, 177
  * Zinchenko et al. (2015) Zinchenko I., et al., 2015, ApJ, 810, 10

##  Appendix A Conversion between the extinction and the gas column density

The extinction maps described in the Section 2.3 were used to estimate the
column density in S254-S258 region. For the conversion between NIR (AJ) and
visual (AV) extinction we adopted the coefficients from Cardelli et al.
(1989): $A_{V}=3.54\times A_{J}=5.26\times A_{H}=8.77\times A_{K}$. To compute
the hydrogen column density from the visual extinction we use the standard
conversion of Bohlin et al. (1978): $N(\mathrm{H})/E_{\mathrm{B-V}}=5.8\times
10^{21}~{}(\mathrm{cm}^{-2}\mathrm{mag^{-1})}$. With the typical ratio of
total to selective extinction $R_{V}=A_{V}/E_{\mathrm{B-V}}\simeq 3.1$
(Cardelli et al., 1989) this gives a hydrogen column density:
$N(\mathrm{H})/A_{V}=1.87\pm 0.13\times
10^{21}~{}(\mathrm{cm}^{-2}\mathrm{mag^{-1})}$

However, in calculation of the gas column density from the extinction, the
conversion factor between $N(H)$ and $A_{V}$ include some uncertainty.
Different values of this ratio were obtained using far-UV extinction
observations (Diplas & Savage, 1994), X-ray sources observations (Predehl &
Schmitt, 1995; Güver & Özel, 2009; Zhu et al., 2017) and interstellar
absorption-line observations (Gudennavar et al., 2012). In the work of
Gudennavar et al. (2012) a survey of $\sim$3000 stars in the Galaxy were used
to determine the conversion factor for total hydrogen (in molecular and atomic
form) and extinction. The resulting value of $N(\mathrm{H})/A_{V}=1.96\pm
0.06\times 10^{21}$ cm-2mag-1 (assuming $R_{V}=3.1$) is close to the canonical
value discussed above. The recent work of Zhu et al. (2017) gives estimation
of $N(\mathrm{H})/A_{V}=2.08\pm 0.02\times 10^{21}$ cm-2mag-1.

##  Appendix B Spatial correlation analysis tests

In this section we examine the behavior of the WWCC analysis between gas-dust
tracers and a random YSOs density map. This is nessesary to ensure that found
correlation coefficients between actual YSOs and gas-dust tracers are not
related to the properties of the WWCC method. In particular, we examine the
increase of correlation between large-scale convolutions of the maps that
actually have no relation to each other.

The random YSOs maps are presented in Figure 9 together with CO and extinction
maps that we used in the WWCC method testing. The following number of random
YSOs were examined: 500, 250, 150, 100, and 50. We generated the catalog of
random YSOs only once, and each time the limit on the source number was
applied. Each YSOs density map was created using a kernel size of 60 arcsec,
similar to the CO and extinction map used in the analysis.

In Figure 10 the results of WWCC analysis are presented. For 50 and 100 YSOs
maps the test reveals slow decrease of the correlation coefficients at large
scales (> 4 arcmin). This decrease occurs in full region and HCO-mask region.
Maps with larger number of YSOs (150-500) display slow increase of correlation
coefficients at large scales (> 4 arcmin). The increase is visible in the full
region and in the region without S255N.

We conclude that for the whole region and the “without S255N” region at the
scales under four arcmin, the correlation coefficients between gas-dust
tracers and random YSOs maps have values below 0.1. However, at large scales
(up to 8 arcmins), it can increase to 0.2-0.3 or even decrease to negative
values if one used a small number of YSOs (50-100). For “HCO+-mask region” the
variation of the correlation coefficients is significant. Thus it is
impossible to make conclusions on the spatial correlations using such a
limited masked region.

Figure 9: Extinction, CO(1-0) and random sources density maps used for WWCC
method testing.

(a) Random sources test (Whole region) (b) Random sources test (Whole region)

(c) Random sources test (Without S255N) (d) Random sources test (Without
S255N)

(e) Random sources test (HCO+ mask) (f) Random sources test (HCO+ mask)

Figure 10: Result of spatial correlation analysis test random distribution of
YSOs density and gas and dust tracers for the whole region (panel a), for the
whole region without S255N cluster (panel b), and the HCO+ map masked region
(panel c). The plot on the left side display correlation coefficients between
randomly distributed sources and the CO map, right-sided plots display
correlation coefficients between random sources and extinction map. The
vertical red line in each plot displays all studied maps’ beam size (60
arcsec).

##  Appendix C Peak velocity and line widths maps of the high-density tracers.

In the Figure 11 we present the maps of peak velocities and line widths (FWHM)
for HCO+(1-0) and CS(2-1) lines. The values are obtained using a three-point
quadratic fit from the moment task of the miriad software package (Sault et
al., 1995).

Figure 11: Maps of line peak velocities and linewidths (FWHM) for HCO+(1-0)
and CS(2-1) lines. Upper panels refer to HCO+(1-0) line and lower panels refer
to CS(2-1) line. Regions are the same as in the Figure 2.