The PDR structure and kinematics around the compact H ii regions S235 A and S235 C with [C ii], [¹³C ii], [O i], and HCO⁺ line profiles Free

ABSTRACT

1 INTRODUCTION

Massive O and B stars reveal the places of their birth by powerful ultraviolet (UV) radiation. The stellar UV photons illuminate and heat the surrounding gas up to thousands of Kelvin. UV photons from the massive stars ionize hydrogen atoms and lead to the formation of H ii regions (e.g. Strömgren 1939; Kahn 1954, and many others). The H ii regions continue to grow in size, driven by the difference of the thermal pressure between the hot ionized and the cold neutral gas (see Mathews 1965; Spitzer 1978, and many others). Shock waves accompany the expansion and move ahead of the ionization front into the neutral surrounding material. The shock compresses the neutral gas and collects it in a dense neutral shell (see simulations by Tenorio-Tagle 1979; Beltrametti, Tenorio-Tagle & Yorke 1982; Franco, Tenorio-Tagle & Bodenheimer 1990; Garcia-Segura & Franco 1996; Hosokawa & Inutsuka 2006; Krumholz, Stone & Gardiner 2007; Kirsanova, Wiebe & Sobolev 2009; Raga, Cantó & Rodríguez 2012; Tremblin et al. 2014; Bisbas et al. 2015; Williams et al. 2018). A photodissociation region (PDR) is formed between the hot ionized and cold molecular material compressed by the shock (e.g. Tielens & Hollenbach 1985a, b; Tielens et al. 1993; Sternberg & Dalgarno 1995; Hollenbach & Tielens 1999; Le Petit et al. 2006; Röllig et al. 2007; Bron et al. 2018). In spite of the long research history, direct observations of the expansion of the compressed PDR material are still missing, while it has been suggested in works by Gómez et al. (2010), Deharveng & Zavagno (2011), Pilleri et al. (2012), Torres-Flores et al. (2013), Schneider et al. (2018), Mookerjea et al. (2019), Pabst et al. (2019), and Sitnik et al. (2019).

The PDRs around H ii regions often appear as ring-like or arc structures in 8-μm Spitzer and 12-μm WISE images due to emission from bending/stretching modes of polycyclic aromatic hydrocarbons (PAHs) excited by the UV radiation, and also in Herschel images at longer wavelengths due to thermal emission from heated dust grains (e.g. Deharveng, Zavagno & Caplan 2005; Churchwell et al. 2006, 2007; Watson et al. 2009; Deharveng et al. 2010; Anderson et al. 2012; Topchieva et al. 2017). Several thousand PDRs have been identified thanks to the space-based surveys, but the PAH or dust emission does not carry any information about their velocity structure. Moreover, the dust is not directly coupled to the gas, and, in general, gas and dust temperatures differ. The gas cooling occurs predominantly through fine-structure lines of [O i] at 63 μm (hot and dense gas) and [C ii] at 158 μm (warm and less dense gas with temperature ≤200 K and density of the order of 100–1000 cm⁻³; see e. g. Tielens & Hollenbach 1985a; Hollenbach & Tielens 1999; Röllig et al. 2006). Only observations of these spectral lines, obtained with high spectral resolution, allow us to trace gas kinematics and the physical conditions of the gas in the emitting region.

The neutral molecular shells around extended H ii regions are non-uniform and clumpy, so it is difficult, if ever possible, to detect the expansion of H ii regions as a whole using position–velocity (pv) diagrams (e. g. Dirienzo et al. 2012; Anderson et al. 2015; Deharveng et al. 2015; Treviño-Morales et al. 2016; Zhang et al. 2016; Anderson et al. 2019; Kirsanova et al. 2019). Common manifestations of inhomogeneous structure near H ii regions are ‘pillars’, which result when the clumpy medium is carved by the UV radiation of hot stars (see e. g. Smith et al. 2010; Flagey et al. 2011; Benaglia et al. 2013; Panwar et al. 2019). Moreover, according to theoretical expectations, the expansion velocity (V_exp) of extended H ii regions embedded in molecular clouds is ∼1 km s⁻¹ with an age about 0.5–1 Myr (see simulations in Hosokawa (see simulations in Hosokawa & Inutsuka 2006; Zavagno et al. 2007; Kirsanova et al. 2009; Kim, Kim & Ostriker 2016; Akimkin et al. 2015, 2017), which is less than the typical velocity dispersions observed in molecular gas in high-mass star-forming regions. In this case, implicit approaches may help. For example, multiwavelength observations allowed Lopez et al. (2014) to estimate the gas pressure exerted on the surrounding neutral gas shells from the ionized gas and conclude that thermal pressure from the ionized gas dominates the expansion of the shells, rather than pressure from stellar wind or from dust grains accelerated by radiation pressure from the ionizing stars. Direct estimation of V_exp from the pv diagrams can also support this conclusion, as was recently shown by Pabst et al. (2019) for a wind-blown bubble.

The best candidates to search for observational signatures of the expansion are embedded compact H ii regions. The simulations mentioned above predict that V_exp is higher for younger compact H ii regions. Moreover, small compact H ii regions might have geometries that are closer to spherical because, at small radii, they have a lower chance of interacting with density fluctuations in the surrounding ISM, resulting in a more uniform expansion in all directions.

The aim of the present work is to study the geometry and gas kinematics in the PDRs around two compact H ii regions, and test whether their observed properties match the predictions of spherical models. We present observations of the [C ii] and [O i] emission from PDRs, complemented by data of molecular line emission from the surrounding natal molecular cloud. We use the [¹³C ii] line to measure the optical depth of the [C ii] emission to study why the line profiles of the latter line are double-peaked: due to high optical depth or gas kinematics. Recent studies (e.g. Graf et al. 2012; Ossenkopf et al. 2015; Mookerjea et al. 2018, 2019; Okada et al. 2019; Guevara et al. 2020) have shown that the [C ii] lines can be optically thick in dense PDRs.

Anderson et al. (2019) examined the radial velocities of the [C ii] line in the extended H ii region Sh2-235 (S235 hereafter). The velocities coincide with the ¹³CO(2–1) and carbon radio recombination lines (RRLs), but differ from the hydrogen RRLs and the fine-structure [N ii] at 205 μm by up to 4 km s⁻¹. The ionized gas in S235 escapes from the parent molecular cloud and surrounding PDR to the observer because of the absence of foreground dense material. There are two compact H ii regions, S235 A and S235 C, in the same molecular cloud, also mapped in [C ii] by Anderson et al. (2019). The velocity difference between the hydrogen RRLs and the ¹³CO(2–1) lines appears to be insignificant there, so we expect that these compact H ii regions might have geometries that are closer to spherical.

2 TARGET REGIONS

The compact H ii regions S235 A and S235 C are located 8 arcmin south of S235 [(see fig. 1 in the work by Anderson et al. (2019) to find a large-scale view S235 and the surrounding regions)]. All of the H ii regions in the area are embedded in the same giant molecular cloud, G174+2.5 (Evans et al. 1981; Heyer, Carpenter & Ladd 1996; Kirsanova et al. 2014; Bieging et al. 2016; Ladeyschikov et al. 2016). For the distance to the H ii regions, we adopt a value of 1.6 kpc from the maser parallax measurements by Burns et al. (2015), who obtained |$1.56^{+0.09}_{-0.08}$| kpc. Consequently, 10 arcsec in our observations translate into a physical scale of 0.08 pc. We also checked the possibility of using Gaia DR2 (Gaia Collaboration et al. 2016, 2018) parallax measurements, in combination with the catalog of inferred distances to Gaia sources by Bailer-Jones et al. (2018) and found that the astrometric solutions for the ionizing stars in S235 A and C are poor. For example, the renormalized unit weight error (RUWE; see Gaia technical note GAIA-C3-TN-LU-LL-124-01), which characterizes the quality of the astrometric data, is greater than 1.3 for both sources. Due to the relative faintness of stars in the region (G = 18–20 mag), the corresponding large uncertainties in the Gaia parallaxes do not allow for a reliable distance determination.

Using infrared spectroscopy, Thompson, Thronson & Campbell (1983) found that the ionizing star in S235 A (which is also known as IRS 3 source in the work by Evans et al. 1981) has spectral type of B0.5, and showed that the star is nearly on the zero-age main sequence. They measured a visual extinction A_V ≈ 10 mag in S235 A. Thompson et al. (1983), and later Felli et al. (2004), determined the luminosity of the ionizing star in S235 A to be L_* = 1.1 × 10⁴ L_⊙. For the ionizing star in S235 C, Thompson et al. (1983) found L_* = (0.8–1.2) × 10⁴ L_⊙ and a spectral type from B3 to B0.3. According to the calibration tables from Smith, Norris & Crowther (2002) and Pecaut & Mamajek (2013), these luminosities correspond to stellar effective temperatures in the range T_eff = 25000–30000 K.

There are four bright point-source objects in the 3.6-μm Spitzer image in the area around S235 A and S235 C. Their coordinates are given in Table 1. The sources S235 A^⋆ and S235 C^⋆ are the ionizing stars of the corresponding H ii regions (see also Felli et al. 1997; Chavarría et al. 2014; Bieging et al. 2016). S235 B is a reflection nebula near a Herbig Be-type star S235B^⋆ (Boley et al. 2009). S235 A-2^⋆ is projected on a PDR to the south-east of S235 A^⋆ and visible in K-band continuum emission (e.g. Felli et al. 1997; Klein et al. 2005).

Glushkov, Denisyuk & Karyagina (1975) determined the electron density (n_e) in S235 A using optical spectroscopy of [S ii] lines and found n_e at selected positions of the order of 10³–10⁴ cm⁻³, as well as a non-uniform distribution of A_V between 6 and 12 mag. Israel & Felli (1978) determined size (S) and n_e of the H ii regions using aperture synthesis observations at 5.0 and 1.4 GHz. They found S = 0.28 and 0.34 pc for S235 A and S235 C, and n_e = 10³ and 4 × 10² cm⁻³ for S235 A and S235 C, respectively. Felli et al. (1997) obtained an image of Br γ emission in S235 A and found a brightness distribution consisting of two spots, with the ionizing star situated between them. Visual inspection of the 3.6-μm image, shown in Fig. 1, gives a projected distance of 0.17–0.20 pc from the ionizing stars in S235 A and S235 C, respectively, to the bright arcs delineating the PDRs from the east and south-east side, and a maximum of about 0.3 pc for the west side of S235 C.

Figure 1.

Left-hand panel: the 3.6-μm Spitzer image of the region (only pixels with S/N ratio higher than 3σ; Fazio & Megeath 2004) with the positions where [¹³C ii] line was observed shown by the white empty circles. Right-hand panel: integrated intensity of the [C ii] 158-μm emission. White dashed lines show cuts for the pv diagrams where the numbers give the angles measured west of south, with the ionizing star at the origin. The integrated intensity of the HCO⁺(3–2) emission is overlaid in black contours at 2, 6, 10, and 14 K km s⁻¹. Point sources from Table 1 are shown by black stars on both panels, labelled by black symbols on the right-hand panel.

Open in new tabDownload slide

3 OBSERVATIONAL DATA

3.1 SOFIA C ii, ¹³C ii, and O i observations

SOFIA/upGREAT¹ mapping observations of the [C ii] 158-μm line were done in 2016–2017 and are described by Anderson et al. (2019). The typical noise level of that data is 5 K at 0.385 km s⁻¹ resolution on a T_mb scale. This noise level only allowed for a tentative detection of the [¹³C ii] line. For a reliable detection, deep integrations were performed towards three positions in guaranteed time on 2018 November 21 and 23, Project-ID 83_0624. For these integrations, we chose positions where the [C ii] line showed some double peak structure, to find out whether that structure is reproduced in [¹³C ii]. The center of the upGREAT Low-Frequency Array (LFA) was pointed towards |$\alpha _1({\rm J2000.0}) = 05^{\rm h}40^{\rm m}53{_{.}^{\rm s}}17$|⁠, δ₁(J2000.0) = +35°42^′11|${_{.}^{\prime\prime}}$|03, and |$\alpha _2({\rm J2000.0}) = 05^{\rm h}40^{\rm m}42{_{.}^{\rm s}}14$|⁠, δ₂(J2000.0) = +35°41^′56|${_{.}^{\prime\prime}}$|99 in S235 A and towards the peak of the [C ii] emission at |$\alpha ({\rm J2000.0}) = 05^{\rm h}40^{\rm m}52{_{.}^{\rm s}}09$| and δ(J2000.0) = +35°38^′37|${_{.}^{\prime\prime}}$|03 in S235 C. The signal was integrated for 115 s at each position. While the detection of [¹³C ii] was the main goal of these observations, the [O i] 63-μm line was observed simultaneously with the High-Frequency Array (HFA) as a side product. Because of the different footprints of the two arrays on the sky, only the central positions match, so that we can compare the [O i] data with the [¹³C ii] data only for the central pixels.

The data were calibrated by the standard GREAT pipeline (Guan et al. 2012) to convert the observed counts into the main beam temperature (T_mb) using a main beam efficiency of η_mb = 0.68 for [¹³C ii] and η_mb = 0.7 for [O i]. For [¹³C ii], both polarizations were averaged. We subtracted linear baselines and spectrally resampled the data to 0.192 km s⁻¹ channel width for [¹³C ii] and 0.14 km s⁻¹ for [O i]. The spatial resolution of the [¹³C ii] data is 14.6 arcsec, while the spatial resolution of the [O i] data is 5.9 arcsec. The sensitivity varies slightly between the individual pixels of the arrays, leading to noise levels of 0.25 K for the best pixels of the LFA and 0.45 K for the worst pixel and noise levels of 0.5 K for the central pixel of the HFA. The [C ii] line was always detected simultaneously with [¹³C ii] with the same noise level at the same positions.

3.2 SMT HCO⁺(3–2) observations

We use previously unpublished observations of the HCO⁺(3–2) emission at 267 GHz made with the 10-m Heinrich Hertz Submillimeter Telescope taken on 2009 April 21 (S235 A and B region) and May 5 (S235C). The observations were made in the on-the-fly (OTF) observing mode, where the telescope scans continuously in a raster pattern with spectra sampled every 0.1 s. For these maps, the scan rate was 5 arcsec s⁻¹, and the spectra were binned by four samples, for an effective sample spacing of 2 arcsec. The scan direction was along Right Ascension and the row spacing was 5 arcsec. The telescope resolution is 28 arcsec (full width at half-maximum), so the maps are well sampled in both coordinates. The map size for the S235AB region was 3.5×3.5 arcmin² centered at |$\alpha ({\rm J2000.0}) = 05^{\rm h}40^{\rm m}51{_{.}^{\rm s}}94$| and δ(J2000.0) = +35°41^′35|${_{.}^{\prime\prime}}$|49. The map of the S235C region was 3×3 arcmin² centered at |$\alpha ({\rm J2000.0}) = 05^{\rm h}40^{\rm m}52{_{.}^{\rm s}}5$| and δ(J2000.0) = +35°38^′24|${_{.}^{\prime\prime}}$|⁠. The receiver was the sideband-separating mixer system as described in Bieging et al. (2016), with sideband separations of about 15 dB, which provides excellent intensity calibration using the standard method (see Bieging et al. 2016). Typical system temperatures were 280–300 K (SSB). The spectrometer was a filterbank with 128 channels of 250-kHz filter width and channel spacing, giving a velocity resolution for the HCO⁺ J = 3–2 line of 0.281 km s⁻¹ covering 36 km s⁻¹ total bandwidth. Typical rms noise levels in each spectrometer channel were about 0.2 K.

3.3 Archival CO and ¹³CO data

We use CO(2–1) and ¹³CO(2–1) data from Bieging et al. (2016) to calculate the CO column density in the same directions where we have [C ii] and [¹³C ii] data. The spatial resolution of the CO data is 38 arcsec and the spectral resolution of the data is 0.3 km s⁻¹.

For the velocity channels where ¹³CO(2–1) line intensity is not higher than 3σ of the noise, we use 〈T_ex〉 averaged over the optically thick channels and equation (5) to calculate N_CO. The uncertainties of N_CO were estimated using standard error propagation formulas.

3.4 Archival continuum data

In order to determine the dust column density from the continuum spectral energy distributions, we collected data from two different instruments. First, archival maps of continuum emission from the Bolocam Galactic Plane Survey (BGPS) were downloaded.² We use the SHARC-II maps at 350 μm (Merello et al. 2015) and the BGPS v2.1 data at 1.1 mm (Ginsburg et al. 2013) without additional manual reduction. The spatial resolution of the BGPS data is 8.5 and 33 arcsec for 350 μm and 1.1 mm, respectively. The SHARC-II camera has 2.59×0.97 arcmin² field of view but is not sensitive to extended emission on scales larger than ≈1 arcmin. For S235 A and S235 C, this corresponds to about 0.5 pc.

Secondly, we use SCUBA-2 data at 450 and 850 μm downloaded from the JCMT Science Archive,³ taken as part of the SCUBA-2 Ambitious Sky Survey (SASSy; Thompson et al. 2007; Nettke et al. 2017). The spatial resolution of the data is ≈9 and 14 arcsec, respectively. The SCUBA-2 camera covers a 45 arcmin² field of view (Holland et al. 2013). We manually reduced the SCUBA-2 data using the 2018A release of starlink (Currie et al. 2014) and versions 1.6.1 of the smurf (Chapin et al. 2013) and 2.5–8 of the kappa packages through the ORAC-DR (Cavanagh et al. 2008) pipeline.

All continuum data were smoothed to the resolution of 33 arcsec provided by the Bolocam data on a common grid. The maps from SCUBA-2 and SHARC-II showed some large-scale negative flux values with a magnitude above the indicated flux errors. These probably stem from the removal of large-scale emission in the chopping observations. We removed those negative values by adding a constant flux to the maps so that no large-scale negative features above the indicated calibration uncertainty remained. This required adding 1.9 mJy arcsec⁻² at 350 μm, 5.0 mJy arcsec⁻² at 450 μm, and 0.1 mJy arcsec⁻² at 850 μm. We added those offsets to the uncertainties in the fitting. The resulting maps then showed typical fluxes around zero outside of S235 AB and S235 C.

The four maps were fitted by a grey-body SED fit using an emissivity law of |$\kappa = \kappa _{250} (\lambda /250\,\rm{\mu m}{})^{-\beta }$| with κ₂₅₀ = 2.16 × 10⁻²⁵ cm⁻² per hydrogen proton and β = 1.8 (Juvela et al. 2015). To convert the hydrogen column density obtained from the dust emission maps or from CO emission maps into an equivalent visual extinction A_V, measured in magnitudes, we use the standard ratio of 1.87 × 10²¹ H-protons per cm⁻²/A_V from Bohlin, Savage & Drake (1978). Variations of β have little impact on the derived column density as the main constraints stems from the long-wavelength data. The temperature fit was constrained to values of at least 9 K, a typical value outside of the sources. This excludes column density spikes in regions with low fluxes due to noise in the data. The uncertainties in the derived equivalent A_V from the fit fall well below 1 mag for all observed points in S235 A and grow up to equivalent A_V ≈ 2 mag for the two most western points in S235 C.

4 OBSERVATIONAL VIEW OF THE PDRs

4.1 Spatial distribution of the line emission

The [C ii] 158-μm emission is detected in both H ii regions as already reported by Anderson et al. (2019). The map of the [C ii] 158-μm integrated line emission is shown in Fig. 1. Peak of the [C ii] emission is shifted by about 10 arcsec to the north-west of the ionizing star in S235 A, corresponding to the physical distance of 0.08 pc, while the emission peak coincides with the ionizing star position in S235 C. Both S235 A and S235 C have elliptical shapes in the [C ii] images, with an axial ratio ≈0.8. S235 C also shows a diffuse tail extending to the north-west. The orientation and size of the [C ii] ellipsoids are similar to the 3.6 μm Spitzer arc-like images of the objects shown in the Fig. 1. There is no enhancement of the [C ii] emission in the direction of the stellar cluster that is situated between the H ii regions.

The spatial distribution of the HCO⁺(3–2) emission significantly differs from [C ii]: Molecular gas is concentrated in the south and south-eastern part of S235 A and in the south-eastern part of S235 C. There are no peaks of the HCO⁺(3–2) emission in the north and north-western parts of the H ii regions. The HCO⁺(3–2) peak around S235 A is situated between the point sources S235 A-2^⋆ and S235 B^⋆ and partially following the bright shell-like emission at 3.6 μm (Fig. 1a). The estimated angular size of the HCO⁺(3–2) cloud around S235 A, calculated through the contour at 2 K km s⁻¹, is 70×100 arcsec² (0.5 × 0.8 pc) elongated in the north–south direction.

The [¹³C ii] spectra in 14 positions of S235 A and in 7 positions of S235 C are shown in Figs 2 and 3. The majority of the detected [¹³C ii] lines show a single peak in contrast to the [C ii] lines that are mostly double-peaked. The brightest [¹³C ii] line has T_mb=4.7 K in S235 A at position 9 with (Δα, Δδ)  = (–12 arcsec, –15 arcsec), while the integrated [¹³C ii] intensities are the same in this position and near S235 A-2^⋆ at position 6 with offset (Δα, Δδ)  = (0 arcsec, 0 arcsec). We have not detected any [¹³C ii] with a signal-to-noise ratio (S/N) > 3 near S235 B^⋆ at position 7 with (Δα, Δδ)  = (–7 arcsec, –33 arcsec). The brightest [¹³C ii] line in the S235 C PDR has T_mb=1.7 K at position 6 with (Δα, Δδ)  = (–16 arcsec, 27 arcsec). The average integrated and line peak intensities are shown in Table 3. These values will be used in Section 5 as selection criteria for theoretical models. We note that the scatter in the [C ii] and [¹³C ii] line intensities is significantly higher in S235 A. Parameters of Gaussian fits to the [¹³C ii] line for each position are given in Table 2.

Figure 2.

[¹³C ii], [C ii], and HCO⁺(3–2) spectra in S235 A shown by red, black, and blue lines, respectively. The [¹³C ii] spectra are the result of the combined three hyperfine components. The velocity range is from –23 to –11 km s⁻¹ in the plots. Offsets of the spectra in arcsec (labelled by black) are given relative to |$\alpha _1({\rm J2000.0}) = 05^{\rm h}40^{\rm m}53{_{.}^{\rm s}}17$| and δ₁(J2000.0) = +35°42^′11|${_{.}^{\prime\prime}}$|03 (spectrum 6). The [C ii] spectra are divided by a factor of 17. The red vertical line corresponds to a velocity of –16.6 km s⁻¹ to show the velocity shift of the spectra. The excursion at the blue end of some of the [¹³C ii] spectra is due to contamination by the red wing of the nearby [C ii] line.

Open in new tabDownload slide

Figure 3.

[¹³C ii], [C ii], and HCO⁺(3–2) spectra in S235 C (right-hand side) shown by red, black, and blue lines, respectively. The [¹³C ii] spectra are the result of the combined three hyperfine components. The velocity range is from –23 to –11 km s⁻¹ in the plots. Offsets of the spectra (labelled by black) are given relative to |$\alpha ({\rm J2000.0}) = 05^{\rm h}40^{\rm m}52{_{.}^{\rm s}}09$| and δ(J2000.0) = +35°38^′37|${_{.}^{\prime\prime}}$|03 (spectrum 4). The [C ii] spectra are divided by a factor of 25. The red vertical line corresponds to a velocity of –16.6 km s⁻¹ to show the velocity shift of the spectra. The excursion at the blue end of some of the [¹³C ii] spectra is due to contamination by the red wing of the nearby [C ii] line.

Open in new tabDownload slide

Due to different array configurations at 63 and 158 μm, the offsets of the [O i] spectra are different from the [¹³C ii] offsets except for the central pixel positions. Fig. 4 shows the [O i] spectra, concentrated in a more compact area compared to the [¹³C ii] spectra. The maximum detected intensity of the [O i] emission has T_mb ≈ 45 K south of S235 A^⋆. The maximum of the [O i] emission has T_mb ≈ 13 K in S235 C south-east of the ionizing star. The [O i] lines have double-peaked profiles, which seems to be a common feature in massive star-forming regions. They are typically interpreted as being due to self-absorption (see, e.g. the recent study of the self-absorbed [O i] line in S 106 by Schneider et al. 2018). The average values of integrated and line peak intensities are shown in Table 3. The averaged values of the [C ii] integrated intensities over the area observed in the [O i] line give almost the same numbers in both PDRs. The resulting average ratio of the [O i] to [C ii] integrated intensities is 0.14 ± 0.04 and 0.04 ± 0.02 in S235 A and S235 C, where the difference in the ratio is almost exclusively due to the different [O i] intensity. Since the [O i] emission is a good tracer of the gas number density, we conclude that these PDRs have different densities, with the density of neutral gas being higher in S235 A. The HCO⁺(3–2) emission, which has a higher critical density compared to the [C ii] line, is strong towards the south-east of both PDRs, with the highest peak around S235 A. In contrast to the [C ii] emission peaks, the peaks on the HCO⁺(3–2) emission map appear in the densest part of PDRs, in the same region where we find bright arcs at 3.6 μm (Fig. 1).

Figure 4.

The [O i] spectra in S235 A (left-hand panel) and S235 C (right-hand panel) superimposed on the Spitzer 3.6-μm image. The velocity range is from –24 to –9 km s⁻¹ in the plots.

Open in new tabDownload slide

4.2 Kinematics of the PDRs

We use the [C ii] position–velocity (pv) diagrams to study gas kinematics in the PDRs. The pv diagrams are computed along the major and minor axes of the elongated [C ii] image of S235 A and in the same directions for S235 C, as shown in Fig. 1 by the dashed lines. The pv diagrams are shown in Fig. 5. The [C ii] lines have double-peaked profiles with the redshifted component around −15 < V_lsr < −14km s⁻¹ and the blueshifted component around −18 < V_lsr < −17km s⁻¹ in both H ii regions. The double-peaked profiles have a dip at V_lsr = −16 km s⁻¹ corresponding to the velocity peak of hydrogen RRL emission; see Anderson et al. (2019). In contrast, the pv diagrams of HCO⁺(3–2) emission show single-peaked line profiles with the peaks falling between the [C ii] peaks. A truncated arc-like shape is visible at the high-density regions in the south-east area of the regions. In the 135° pv diagram through S235 A, the position where the [C ii] emission is sharply confined coincides with the direction of the point source S235 A-2^⋆. The source is projected on the bright arc-like structure visible on the Spitzer images in Fig. 1 (see also images in Dewangan & Anandarao 2011).

Figure 5.

pv diagrams of the [C ii] (colour) and HCO⁺(3–2) (white contours) emission for S235 A (top panels) and S235 C (bottom panels). Dashed vertical lines correspond to positions of the point sources listed in Table 1, from the beginning of the pv diagrams.

Open in new tabDownload slide

The velocity difference between the peaks of the double-peaked [C ii] line and the peak of the HCO⁺(3–2) line is ≈2 km s⁻¹, as can be seen from the pv diagrams. A simple combination of the [C ii] and HCO⁺(3–2) pv diagrams seems to indicate an expansion of the [C ii]-emitting gas of the PDR with V_exp ≈ 2 km s⁻¹ into the quiescent surrounding molecular cloud visible in the HCO⁺(3–2) line. However, to study the gas kinematics in the atomic and molecular layers of the PDRs in details, we extracted HCO⁺(3–2) spectra at the same positions where we have [¹³C ii], and show them in Figs 2 and 3. The single-peaked [¹³C ii] line profiles shown in these figures bring us to the conclusion that the [C ii] dips are related to self-absorption rather than to the expansion of the front and back neutral walls of the PDRs.

The brightest [¹³C ii] line in S235 A has a symmetric single-peaked profile at a velocity V_lsr = −16.6 km s⁻¹. This value of V_lsr is shown by the red vertical line on each spectrum in Fig. 2 to emphasize the [¹³C ii] velocity gradient across the H ii regions. There is an ≈1 km s⁻¹ global west-–east velocity gradient visible in the spectra of the [¹³C ii] lines in S235 C (at positions 1 and 7, the gradient is clearly visible) and a south-west–north-east gradient of a similar, but more uncertain, magnitude in S235 A (positions 3 and 11). The [¹³C ii] spectra in the other directions of S235 A have a complex shape: flat-top in position 6, redshifted and red-skewed in position 3, blue-skewed broad profiles to the north-west of the S235 A^⋆ in positions 8, 11, and 13. Properties of the [C ii] line emission in the directions observed in the [¹³C ii] line are given in Table 4. It is seen that the |$\tau _{\rm C^+} \gt 1$| in all directions where we detect [¹³C ii] at a level higher than 3σ. Average T_ex of the [C ii] and [¹³C ii] line emission is from 40 to about 90 K. We conclude that both H ii regions are embedded in dense and warm PDRs and we observed them through the layer optically thick in the [C ii] line.

The detected [O i] lines have a self-absorption dip, which appears at almost the same velocity (–16 km s⁻¹) in the [C ii] lines. The self absorption is more significant in the [O i] spectra than in the [C ii] spectra in both PDRs (see also averaged over the observed positions spectra in Fig. A1 in Appendix A). The velocity difference between the two line peaks is the same as in the [C ii] lines: ≈4 km s⁻¹. The double-peaked [O i] lines have a blueshifted component brighter than the redshifted component in almost all positions observed in S235 A and S235 C, except for four directions to the south and east of the point source S235 A-2^⋆. Lines with enhanced blue emission are often attributed to infall in centrally condensed sources (Dickel & Auer 1994), but can also be due to inhomogeneous density distributions. The region with a brighter redshifted (blueshifted) component of the [O i] lines coincides with the region of the redshifted (blueshifted) [¹³C ii] lines in S235 A. This indicates that the effect traces the global velocity gradient in the PDRs instead of local expansion or contraction. The self-absorption dip has the same V_lsr for the all spectra in S235 A and S235 C, which implies that a layer of cool gas is situated in front of the warm PDRs, possessing a velocity gradient formed around the H ii regions.

The HCO⁺(3–2) lines have a velocity gradient with similar direction and sign as the [¹³C ii] lines in both PDRs. Consequently, this gradient is indeed related to the properties of the parent molecular cloud where the H ii regions were formed. We fit the HCO⁺(3–2) and [¹³C ii] line shapes with single-peaked Gaussian profiles (see Table 2), and find that the [¹³C ii] line peaks are redshifted relative to HCO⁺(3–2) by about 1 km s⁻¹ at positions 3, 6, and 9 south-east of S235 A, and also at positions 1 and 3 east of S235 C. We interpret this velocity difference as relative motion of the HCO⁺(3–2)-emitting dense molecular material, surrounding the [C ii]-emitting PDR, with ≈ 1 km s⁻¹ in the south-east parts of the PDRs. We find a similar HCO⁺(3–2)-[¹³C ii] velocity shift in the averaged spectra of S235 A in Fig. A1, and an excess of redshifted [¹³C ii] emission compared to HCO⁺(3–2) in S235 C, visible in positions 1 and 3 in Fig. 3, as well as in the averaged spectra.

The lack of a velocity difference between [¹³C ii] and HCO⁺(3–2) in the western parts of both PDRs may be due to the absence of high-density molecular material. The [C ii]- and HCO⁺(3–2)-emitting gas, apparently, flows freely (i.e. with no dense molecular shell formed by shock wave ahead of the ionization front) from the PDRs into the surrounding low-density medium at these locations. The red wings of the [C ii] line visible in positions 8, 11, and 13 in S235 A and in positions 6 and 7 in S235 C indicate that expansion into the low-density medium occurs in the opposite from the observed direction.

The value of V_exp ≈ 1 km s⁻¹ found from the analysis of the [¹³C ii] and HCO⁺(3–2) spectra is about two times less than the value that would be estimated from the [C ii] and HCO⁺(3–2) pv diagrams. Below, in Section 5, we compare the observed gas kinematics with predictions of a numerical model of expanding H ii regions, and show simulated pv diagrams.

4.3 Column densities of C-bearing species

Column density maps of CO and HCO⁺ molecules are shown in Fig. 6. Both maps reveal a concentration of molecular gas to the south and south-east of S235 A in between the point sources S235 A-2^⋆ and S235B^⋆, and to the south-east of S235 C. While the regions with the highest N_CO and N_HCO⁺ around S235 A and S235B^⋆ coincide in the plane of the sky, there is a decrease in the intensity of N_CO near the secondary peak of N_HCO⁺ to the south of S235B^⋆ (around α ≈ 05^h40^m54^s and δ ≈ +35°40^m30^s).

Figure 6.

Left-hand panel: maps of N_CO (colour) and N_HCO⁺ (contours). The contours correspond to 1.0, 2.5, 4.0, and 5.5 × 10¹² cm⁻². Directions of the [¹³C ii] observations are shown by circles. Red circles stand for positions with |$N_{\rm C^+} \gt N_{\rm CO}$|⁠, white circles show the opposite relation. Middle panel: dust column density in equivalent visual extinction A_V scale. The column density maps are given at the original spatial resolution of the observations (38 arcsec for |$N_{\rm CO}$| in the left-hand panel, and 33 arcsec for the middle and right-hand panels). Right-hand panel: dust temperature. Point sources from Table 1 are shown by black stars in the left-hand panel and by red stars in the middle and right-hand panels.

Open in new tabDownload slide

The dust temperature (T_dust) and column density (N_dust) maps are also shown in Fig. 6. The dust column density structure roughly follows the spatial distribution of N_HCO⁺. The T_dust value is enhanced up to about 24 K close to the embedded sources, but the temperature peaks do not agree with the column density peaks but are rather offset. For S235 A, we find the temperature peak south of the main ionizing star, probably due to the contribution from the other embedded sources (e. g. S235 A-2^⋆, S235 B^⋆); in S235 C, the temperature peak is offset to the west, probably due to the density gradient allowing for a more efficient heating of the eastern part of the source. The dust temperature agrees with excitation temperature of CO molecules found by Bieging et al. (2016). These molecules are tracers of the general distribution of gas with relatively low critical density. We will note about higher values of the gas temperatures obtained with ammonia lines in Section 7.

We compare the C⁺ and CO column densities in the directions where the [¹³C ii] was observed in Table 4. Using the dust continuum emission maps described in Section 3.4 and assuming that all carbon is in C⁺ and CO species with the total carbon abundance 1.2 × 10⁻⁴ (see Section 5), we recalculate hydrogen column density in terms of equivalent A_V value (⁠|$N_{\rm H}^{\rm C}$|⁠) and compare it with the dust column density (⁠|$N_{\rm H}^{\rm d}$|⁠) in Table 4, illustrated in Fig. 7. While in the majority of the directions both column densities are in agreement, we find several directions where |$N_{\rm H}^{\rm d} \gt (1.5\!-\!2) N_{\rm H}^{\rm C}$|⁠. These directions surround the dense molecular clump around S235 A from the south, and also to the north-west, as well as in one position to the east of S235 C.

Figure 7.

Top panel: carbon and hydrogen column densities in terms of equivalent A_V. Uncertainty of the hydrogen column density is within 1 and 2 mag in S235 A and S235 C, respectively (see Section 3). Bottom panel: N_C⁺ and N_CO.

Open in new tabDownload slide

To explain this disagreement, we checked our assumptions in the derivation of column density. First, using a map of the optical depth in ¹³CO(2–1) from Bieging et al. (2016), we find |$\tau _{\rm ^{13}CO(2-1)} \approx 1$| in the molecular cloud around S235 A and C. This approach underestimates the N_CO value, which partially explains the missing carbon, at least in positions 4, 7, and 10 in S235 A. Comparison of the panels of Fig. 6 with the N_CO and |$N_{\rm H}^{\rm d}$| distributions shows a relative deficiency of CO in the gas phase south of S235 B^⋆, where the peaks of the |$N_{\rm H}^{\rm d}$| distribution are not followed by the N_CO peaks. Because CO and HCO⁺ are chemically closely related, they should freeze out on to the dust grains simultaneously (see e.g. Bergin & Langer 1997). Thus, we associate the decreasing CO column density to the south of S235 B with the optical depth of the ¹³CO(2–1) line, rather than with freeze-out on to the dust grains, because a secondary peak of N_HCO⁺ is seen at this location.

A second factor contributing to the total carbon column are variations of T_ex for the [C ii] line along the line of sight. We can estimate them by comparing the maximum of the channel T_ex with 〈T_ex〉 in equation (4). The maximum channel T_ex would raise |$N_{\rm C^+}$| only within the interval given in Table 4 included in the upper error bars in Fig. 7.

Looking at S235 A from the north-west to the south-east we distinguish two different types of environments: the PDR where C⁺ is the dominant C-bearing species and the molecular cloud where CO is the dominant C-bearing species. Toward S235 A^⋆, S235 A-2^⋆ and in several directions to the west of the source, the C⁺ column density exceeds the CO column density by a significant factor, |$N_{\rm C^+} \approx (2\!-\!6) \times N_{\rm CO}$|⁠. The directions with |$N_{\rm C^+} \gt N_{\rm CO}$| are shown with red circles in the left-hand panel of Fig. 6. The point sources are situated on the border between the PDR and surrounding molecular cloud where the PDR is expanding; see Fig. 6. The column of carbon locked in C⁺ is higher than in CO in five of the seven observed positions in S235 C. The values of N_CO and N_HCO⁺ are highest in the directions of the dense molecular clumps where |$N_{\rm CO} \gt N_{\rm C^+}$|⁠.

If we average the ratio of the column density N_CO + N_C⁺ to the hydrogen column density computed from the dust emission, we find a gas-phase carbon abundance 7.7 × 10⁻⁵ in S235 A and 1.1 × 10⁻⁴ in S235 C, which is in agreement with the total carbon abundance value used above. Some carbon in the directions with |$N_{\rm H}^{\rm d} \gt N_{\rm H}^{\rm C}$| could be in the gas-phase atomic form, in high-J CO lines, or locked in dust grains. In the directions with |$N_{\rm H}^{\rm d} \gt N_{\rm H}^{\rm C}$| around S235 A, we find 10 ≤ T_dust ≤ 20 K. Taking into account that these positions are situated in the surrounding cold molecular cloud, it is natural to expect that some carbon is depleted from the gas phase due to freeze-out on to the grain surfaces. A rich young stellar cluster is situated between the H ii regions (see Section 7), but it does not provide enough heating to make the dust warmer than 10–15 K.

5 NUMERICAL MODELLING

As the H ii regions S235 A and S235 C are deeply embedded in the molecular cloud, we cannot use published predictions from standard plane-parallel PDR models without foreground absorption (e.g. Kaufman et al. 1999). In order to understand how the observational view by SOFIA is matched by a ’classical’ model of expanding H ii regions, described analytically and numerically by Spitzer (1978), Elmegreen & Lada (1977), Hosokawa & Inutsuka (2006), Raga et al. (2012), Kirsanova et al. (2009), and Bisbas et al. (2015), we make simulations with the MARION model (Kirsanova et al. 2009; Akimkin et al. 2015). The gas-phase chemical network from Röllig et al. (2007) (mainly based on the UMIST99 ratefile; Le Teuff, Millar & Markwick 2000) allows us to reproduce the C⁺, CO, and HCO⁺ abundances in PDRs (see Kirsanova & Wiebe 2019), so we use this network together with the ionization of the atomic species and corresponding recombinations in the H ii regions. The cross-sections for most of photoreactions are taken from the Leiden data base of ‘Photodissociation and photoionization of astrophysically relevant molecules’ (Heays, Bosman & van Dishoeck 2017). We implement the formation of H₂ on grain surfaces and accretion and desorption processes of other neutral species, but other chemical reactions on dust surfaces are not considered in the calculations to save computation time. The rates of accretion and desorption processes are based on the work by Hasegawa & Herbst (1993), with updated desorption energies from Garrod, Wakelam & Herbst (2007). We use the ‘high-metallicity’ initial elemental abundances based on the ‘EA2’ set from Wakelam & Herbst (2008). The initial conditions are cold, molecular, and solid. This means that we start with all carbon and oxygen in CO and H₂O on dust surfaces. Fifty two chemical species are included: H, H⁺, H₂, H|$_2^+$|⁠, H|$_3^+$|⁠, O, O⁺, O⁺⁺, OH⁺, OH, O₂, O₂:d, O|$_2^+$|⁠, H₂O, H₂O⁺, H₃O⁺, C, C⁺, C⁺⁺, CH, CH⁺, CH₂, CH|$_2^+$|⁠, CH₃, CH|$_3^+$|⁠, CH₄, CH|$_4^+$|⁠, CH|$_5^+$|⁠, CO, CO:d, CO⁺, HCO⁺, He, He⁺, S, S⁺, S⁺⁺, Si, Si⁺, H:d, H₂:d, O:d, OH:d, H₂O:d, C:d, CH:d, CH₂:d, CH₃:d, CH₄:d, S:d, Si:d, and e⁻, where the postfix ‘:d’ indicates species on dust grain ices. Chemical species containing Si and S are included only to obtain correct gas temperature in ionized region. We recognize that the chemical network is far from complete, but our choise was motivated by limited computation time. The heating and cooling processes included in the model are listed in Akimkin et al. (2015) and Kirsanova & Wiebe (2019).

The dust grains in the model are represented by silicate and carbonaceous dust grains with 24 size bins for each type. The initial dust-to-gas mass ratio is 0.0088 with the size distribution corresponding to model number 16 from Weingartner & Draine (2001). This model has R_V = 5.5 instead of standard interstellar R_V = 3.1, and we use it because this reddening law implies lower relative amount of small dust grains and, therefore, lower dust extinction and thicker PDR (see below). Values of R_V > 5 were found in PDRs such as the Orion Bar and NGC 7023 (see Marconi et al. 1998; Witt et al. 2006, respectively). The dust model implies that the carbon abundance locked in dust refractory cores is 3 × 10⁻⁵. Therefore, the total carbon abundance in the model is 1.5 × 10⁻⁴. The Si abundance locked in the refractory cores corresponds to solar value of 3.63 × 10⁻⁵. The grain cores are not destroyed in our simulations. The size of the computational domain is 0.5 pc, divided into 550 radial cells. The dust temperature is computed self-consistently from radiative processes and interaction with the gas. The gas temperature is also computed self-consistently, except for the initial moment, when it is set up to the value corresponding to cold molecular clouds.

In the beginning of the simulation, it is assumed that a massive star is embedded in a molecular cloud with uniform gas distribution. At t = 0, the star starts to ionize and heat the surrounding material. We scanned a range of effective temperatures of the ionizing stars (25 000–31 000 K) and of initial gas number densities (5 × 10³–10⁵ cm⁻³) to reproduce the observed properties of the H ii regions and surrounding PDRs. In this work, we present the simulations in this parameter range that give the best agreement with the observations. Initial conditions of gas and dust in these simulations are summarized in Table 5.

In order to simulate the profiles of molecular and atomic line emission for the model physical parameters and chemical abundances calculated with MARION, we apply the spherically symmetric line radiative transfer code simline by Ossenkopf, Trojan & Stutzki (2001). The main partner for the collisional excitation are electrons in the H ii region and hydrogen atoms and molecules in the PDR. We treat this in the radiative transfer computation by switching from electrons to hydrogen when the gas temperature drops below T_gas ≤ 500 K. For HCO⁺, the electron excitation in the H ii region is ignored as there is not significant molecular abundance there. Coefficients for collisional excitation of C⁺ (Sternberg & Dalgarno 1995, table 8), atomic O (Jaquet et al. 1992), and HCO⁺ (Flower 1999) line emission with molecular hydrogen were used. The code uses two parameters for the turbulence description: the width of the turbulent velocity dispersion and the turbulence correlation length. They were set to 2.4 km s⁻¹ (based on typical [C ii] and HCO⁺(3–2) line widths; see below) and 0.1 pc (Xie, Allen & Langer 1995; Ossenkopf et al. 2001), respectively.

6 MODEL PREDICTIONS FOR S235 A AND S235 C

As MARION is a chemo-dynamic model, we can trace the evolution of the physical parameters of the H ii regions and select the moment of the best agreement between the simulated and observed values of the parameters. The following properties need to be matched by the model:

• Radii of the H ii regions are from 0.1 to 0.3 pc.

• n_e = 1000 and 500 cm⁻³ within a factor of 2 for S235 A and S235 C, respectively.

• The profiles of the [¹³C ii] and HCO⁺(3–2) lines are single-peaked.

• The average W([¹³C ii]), W([O i]), and the [¹³C ii] line peak emission is close to the observed values (see Tab. 3) within a factor of 2.

• Expansion velocity is ≈1 km s⁻¹.

• |$3 \le \tau _{\rm C^+} \le 10$| and |$3\le \tau _{\rm C^+} \le 5$| for S235 A and S235 C, respectively.

We used the output from MARION every 5000 to 10 000 yr of model time and calculated the resulting pv diagrams. We compared the model peak line intensity, integrated intensity and the line profiles with the observed values. We find that models with T_eff = 27000 K and initial n_gas = 5 × 10⁴ and 5 × 10³ cm⁻³ provide a reasonable agreement for S235 A and S235 C, respectively. The model age of S235 A is 6 × 10⁴ and of S235 C is 3 × 10⁴ yr. Below we compare the observed view of the H ii regions and their PDRs and the simulated models.

6.1 Physical properties of H ii regions

The spatial distributions of the physical conditions and chemical abundances for both H ii regions are shown in Fig. 8. The values of n_e ≈500 and 260 cm⁻³ for S235 A and S235 C, respectively. Each H ii region is surrounded by a dense molecular shell, with maximum number densities of |$n_{\rm H_2}=2\times 10^5$| and 2 × 10⁴ cm⁻³ for S235 A and S235 C, respectively. The shells are created by the shock waves that spread ahead of the ionization fronts on the border of the H ii regions. The width of the compressed shells is less than 0.1 pc in both models. The expansion velocities are about 1 and 3 km s⁻¹ for S235 A and S235 C, respectively.

Figure 8.

Physical conditions in the models of S235 A (left-hand panel) and S235 C (right-hand panel): number densities of species (three top panels of plots), gas and dust temperature, and also gas velocity (bottom panel).

Open in new tabDownload slide

6.2 Chemical composition

The temperature contrast between the PDR and the surrounding molecular region is higher in the S235 A model compared to the model for S235 C due to relatively low gas density in the latter model. The central source of UV radiation heats the gas and dust over a larger distance in the S235 C model than in the case of S235 A. Carbon is in the form of C⁺ inside of the H ii regions, and also in the neutral dense shells, but the outer border of the shell consist of molecular CO in the model for S235 A. The surrounding gas contains CO in the form of grain ice in the model of S235 A, but mostly consists of C⁺ in the model for S235 C. The abundance of atomic carbon is smaller than the abundances of CO or C⁺ by up to several orders of magnitude in both models, except for a very thin layer in the model for S235 A, where the abundance of atomic C is only five times lower than CO and three times lower than C⁺. Oxygen is in atomic form in the dense shell in both PDRs and also in the undisturbed gas in S235 C. The value of T_gas changes from 20 to 300 K in the PDR of the S235 A model (20 < T_dust < 100 K there), but falls below 50 K only at a distance of 0.45 pc in the model for S235 C (corresponding T_dust < 30 K.)

6.3 Averaged intensities

The radial distributions of the integrated line intensities are shown in the top panel of Fig. 9. The selected models predict that the [¹³C ii] line intensities are brighter on the periphery of the PDRs than in the direction on the ionizing star in agreement with the observations in Figs 2 and 3. The simulated peak values of W([C ii]) are 75 and 50 per cent from the observed in S235 A and S235 C, respectively. Our simulated radiation temperatures of the [¹³C ii] lines about 0.5–1 and 0.3–0.7 K for both PDRs respectively, which is in agreement within factor of 2 with the average values of T_mb from Table 3. Integrated intensities of the simulated [O i] lines are in good agreement with the average values from Table 3. The integrated intensities of the HCO⁺(3–2) lines in both models are lower than those observed by about two and three orders of magnitude in S235A and S235C, respectively.

Figure 9.

Simulated radial integrated intensity distributions (top panels) and pv diagrams in S235 A (left-hand panels) and S235 C (right-hand panels). The HCO⁺(3–2) pv diagrams are shown by white contours on the top of the [C ii] pv diagrams. The contours are given every 14 per cent from the normalized HCO⁺(3–2) peak.

Open in new tabDownload slide

In order to check whether collisions with electrons can solve the discrepancy by enhancing the HCO⁺(3–2) line emission, we performed test described in Appendix C (supplementary material available online). We found that the collisional excitation by electrons can increase the peak HCO⁺(3–2) intensity by a factor of 2 and 5 in S235 A and S235 C, respectively. However, this does not solve the problem of the low intensity of the HCO⁺(3–2) line in the models. A match of the observed intensity requires a significantly enhanced column of HCO+. We have tested this in two ways. We re-simulated the HCO⁺(3–2) lines with an enhanced value of molecular hydrogen number density, but an unchanged value of HCO⁺ relative abundance compared to Fig. 8 (i. e. the column density was enhanced) and found a simulated radiation temperature of the HCO⁺(3–2) line of several Kelvins. Namely, the peak of the HCO⁺(3–2) line was 2 and 6 K for a column density enhancement by a factor of 10 and 100, respectively, in the model for S235 A. These values are in agreement with the observations, as the maximum peak of HCO⁺(3–2) intensity is 4K in S235 A (Fig. 2). If we re-scale the relative abundance of HCO⁺ accordingly, while holding N_HCO⁺ constant, we obtain a peak intensity of 0.4 and 1.6 K for a number density enhancement by a factor of 10 and 100, respectively. The latter value is also within a factor of 2.5 of the peak observational value. The first approach combined a column density enhancement with a density enhancement similar to the effect from the additional collisions with electrons. The relatively small difference between the two cases confirms the conclusion that the molecular excitation has a small effect only. It cannot explain the discrepancy between the observations and the model predictions.

6.4 [C ii] and [O i] optical depths

The simulated profiles of the [C ii] emission in both selected models are single-peaked. The value of |$\tau _{\rm C^+}$| does not exceed 0.4 and 2.0 in the models for S235 A and C, respectively. In spite of the significant optical depth in the S235 C model, the amount of cold gas (T_gas ≈ 30 K) with abundant C⁺ is too small to create double-peaked line profiles. In the model of S235 A, we find only slightly red-skewed [C ii] line profiles, instead of the observed double-peaked profiles. Apparently, the low optical depth of the [C ii] emission in the S235 A model is due to the absence of C⁺ beyond the dense shell in the undisturbed molecular gas, where all the carbon is in the form of solid-phase CO because the value of T_dust in the model is too low to evaporate molecules into the gas phase.

To simulate how the |$\tau _{\rm C^+}$| grows with the width of C⁺-layer in PDR, we increased relative abundance of C⁺ in front of the expanding envelope up to maximum possible value and preserved all other physical values unchanged. Two test models were considered, there the C⁺-layer is twice and three times as thick as in original model shown in Fig. 8. Calculated pv diagrams for considered test models are shown in Fig. 10. The [C ii] pv digrams become significantly skewed in the twice thick C⁺-layer and double-peaked in the three times thick C⁺-layer, respectively. The peak-to-peak velocity difference ≈4 km s⁻¹ between the blue and red components is defined not by thermal but by turbulent velocity dispersion in the model. The maximum channel value of |$\tau _{\rm C^+}$| becomes 1.6 and 4.5 for the twice and three times thick models, respectively, in agreement with the observations presented above. We note that in these test calculations, the [¹³C ii] line also becomes double-peaked due to significant optical depth that is not in agreement with observations. Therefore, the self-absorption features observed in S235 A and S235 C might not be formed in the PDRs, where relative abundance of C⁺ is high, but outside of them where the C⁺ abundance is lower but the absorbing layer is thicker.

Figure 10.

Simulated pv diagrams for the S235 A model with enhanced width of the C⁺-layer up to 0.2 pc (original width multiplied by 2, left-hand panel) and 0.25 pc (original width multiplied by 3, right-hand panel).

Open in new tabDownload slide

The simulated [O i] lines have double-peaked profiles due to the significant optical depth (about 3–4 in the velocity channels corresponding to the self-absorption feature).

6.5 Expansion

In the simulated pv diagrams, we find single-peaked [¹³C ii] line profiles in both of the selected models, despite the expansion of the H ii regions. For the model for S235 A, the expansion velocity is lower than the width of the turbulent velocity dispersion assumed in the radiation transfer calculations with SimLine, therefore we do not even see the line wings on the pv diagram in the direction to the ionizing star. For the S235 C model, we have significant amount of stationary, undisturbed gas where C⁺ is the dominant C-bearing species. Broad [¹³C ii] lines with extended wings in the direction of the centre of the pv diagram indicate the expansion in this model. The HCO⁺(3–2) pv diagrams demonstrate the same behaviour in the selected models [but see Appendix B (supplementary material available online) to find more about the [¹³C ii] and HCO⁺(3–2) pv diagrams].

The peak-to-peak velocity difference between the near and far neutral walls on the simulated [O i] line profiles (2 and 6 km s⁻¹ in the models for S235 A and S235 C) corresponds to the expansion with V_exp = 1 and 3 km s⁻¹, respectively. The simulated V_exp value is in agreement with the observations of S235 A, but exceeds the values observed in S235 C for a factor of 3.

As has been shown in various studies (e. g. Troitsky, Lapinov & Zamozdra 2004; Pavlyuchenkov et al. 2008), moderately optically thick lines can be used to discriminate between expansion (double-peaked profiles with a brighter redshifted component), collapse (double-peaked profiles with a brighter blueshifted component), and other effects (e.g. rotation). The [O i] lines are the best tracers of the PDR expansion in the model due to their double-peaked profiles, with a bright redshifted component, caused by a significant optical depth. We note that the simulated [O i] profiles in the expanding PDRs are double-peaked, although the other PDR tracers, namely the [¹³C ii] and [C ii] lines, can be single-peaked.

In the observations, we found the same peak-to-peak velocity difference on the [C ii] and [O i] line profiles up to 4 km s⁻¹. As shown above, the [C ii] self-absorption feature is explained by the layer of cool absorbing C⁺ on the front of the PDRs. We suggest that the broad and deep [O i] self-absorption feature is related to the same cold absorbing layer and, unfortunately, hide the internal gas kinematics in the PDRs. As recently shown by Guevara et al. (2020), the dense and cold foreground C⁺ gas leads to the self-absorption effects in several objects. While they note that the origin of the cold gas is unknown, this effect may prevent the detection of expansion and/or contraction in the PDRs. Therefore, further theoretical studies of reliable PDR expansion tracers among other lines, which do not suffer from self-absorption in the cold foreground, are required.

While we reproduce the average integrated intensities of the observed [¹³C ii], [C ii], and [O i] lines, we do not reproduce the detailed line shapes of [C ii] and [O i]. The model predicts [O i] (double-peaked) and [C ii] (single-peaked) line profiles with brighter redshifted components; however, the observed profiles are blueshifted in the majority of the positions. From the comparison of the optically-thin [¹³C ii] line profiles with the double-peaked self-absorbed [C ii] and [O i] profiles, we know about the east–west velocity gradient, and also about the gas density gradient from the maps of the CO and HCO⁺(3–2) emission. Therefore, the agreement between the observed and simulated line profiles is not perfect. Change from the redshifted to the blueshifted asymmetry on the observed line profiles might be related with the gas velocity gradient in the PDRs. Foreground absorbing cold material forms the self-absorption feature at ≈−16 km s⁻¹ on the observed profiles.

6.6 Summary on the simulated line profiles

The [O i] emission is not only a reliable tracer of the gas kinematics, but also of the gas density in PDRs. The [O i] line at 63 μm was the most useful in our model fitting procedure. Both the observed and simulated [O i] line profiles have the largest peak-to-peak velocity difference (see also averaged spectra in Fig. A1). While the HCO⁺ lines are more sensitive to the gas number density, we can not obtain the brightness temperatures of several K which are observed in S235 A and S235 C due to insufficient number density of molecular hydrogen. The HCO⁺(3–2) line becomes a reliable kinematic tracer if V_exp is comparable, or exceeds, the turbulent velocity dispersion (models with n_init ≤ 10⁴ cm⁻³). The [C ii] and [¹³C ii] lines show no self-absorption in the considered models due to insufficient column density of cold C⁺ around the PDRs.

While we can reproduce the observed physical conditions and the integrated line intensities with these models (as these properties are determined by UV photons input from the ionizing stars and by the initial number density of the surrounding medium), the gas kinematics and the detailed line profiles strongly depend on the geometry of the regions.

7 DISCUSSION

In the previous sections, we considered S235 A and S235 C just as H ii regions surrounded by PDRs and molecular clouds with no star formation activity, which can affect the properties of the gas and dust. However, the vicinity of these H ii regions is characterized by active star formation. Felli et al. (1997) found a highly obscured young stellar cluster between S235 A and B, as well as a highly variable H₂O maser. Felli et al. (2004) observed various molecular transitions as well as mm and sub-mm continuum emission, and found dense molecular cores with new YSOs responsible for driving the maser activity. Aperture photometry of IRAC point sources done by Felli et al. (2006) showed that more-evolved Class I and II YSOs surround the PDR of S325 A, and there is a younger intermediate-mass protostar embedded in the dense molecular gas and located further from the PDR than the YSOs. With spatially resolved VLA observations, they found a burst of active star formation in the region of S235 A and B, where different evolutionary phases co-exist: from proto-stellar cores to H ii regions. Felli et al. (2006) compared the Spitzer IRAC images of S235 A and B with the spatial distribution of the HCO⁺(1–0) emission and proposed that the H ii region S235 A and the surrounding molecular cloud might be interacting. They also proposed that dynamic action of S235 A may have triggered a second generation of star formation in the surrounding molecular gas. Kirsanova et al. (2008) found that the region of S235 A, S235 B, and S235 C contains the densest and richest cluster of young stars among the whole giant molecular cloud G174+2.5 (see also Dewangan & Anandarao 2011, who identified 55 Class I, 134 Class II, and 51 Class III objects around S235 A and B using Spitzer photometry). We suggest that ionization of CO-bearing molecular gas around the H ii regions by nearby young stellar sources might provide cool C⁺, which is required to form the self-absorbed [C ii] line profiles and also provide additional HCO⁺. One of the most promising examples of this is the poorly studied Class I YSO S235 A-2^⋆ identified by Dewangan & Anandarao (2011) (their source number 5 in table C5). Its contribution to the heating of the gas is unknown, although its proximity to the peaks in density and dust temperature suggest that it might provide a significant contribution.

If we remove the assumption of spherical symmetry, the velocity difference between the optically-thin [¹³C ii] and HCO⁺(3–2) emission in both PDRs can be interpreted in terms of expansion of the PDRs into the surrounding molecular gas. As [¹³C ii] is redshifted relative to the HCO⁺(3–2) line, and the total column in the front wall is smaller than in the rear wall (compare A_V = 10 mag for the foreground gas from Thompson et al. (1983) and the equivalent A_V values from Fig. 7), the PDRs expand into the front molecular HCO⁺-emitting walls of the H ii regions. The absence of the signature of gas motion from the back walls is related to the high column density of the background gas or to the geometry of the PDRs. In S235 A and S235 C, the transitions between diffuse and dense gas occurs in an obviously highly inhomogeneous medium. Moreover, we observe a transition between the CO-rich and CO-frozen molecular gas in S235 A. This compact H ii region appears as a good example of a PDR with a snow line, proposed by Goicoechea et al. (2016), while the H ii region is embedded and not visible face on (e.g. like the Orion Bar PDR).

The model of S235 A does not reproduce the observed high C⁺ column density, while the relatively low gas density in the model of S235 C allows us to obtain the moderately optically thick C⁺ environment of the H ii region, due to large width of C⁺-emitting layer. We propose that there is a highly inhomogeneous and clumpy medium around the H ii regions, which is unresolved in the [C ii] observations. The clumpy structure may be composed, for example, of evaporating globules. These globules could allow the photons from the ionizing star to penetrate deeper into the surrounding cloud and increase the C⁺ column. Signatures of the photon penetration through the clumpy medium were found in, for example, the extended H ii region RCW 120 by Anderson et al. (2010, 2015) and Kirsanova et al. (2019), instead of uniform expansion, and also at the very early stage ultracompact H ii region in Mon R2 by Treviño-Morales et al. (2016). A uniform spherically symmetric H ii region, considered in the classic theoretical studies (see Section 1), has not been found so far.

Interestingly, not only S235 A and S235 C have higher columns in the back wall, but also other H ii regions situated in the giant molecular cloud G174+2.5. Anderson et al. (2019) and Kirsanova et al. (2020) found the same for the extended H ii region S235. Ladeyschikov et al. (2015) determined that the H ii region Sh2-233 is also situated on the near side of the molecular cloud. Visual inspection of the CO emission maps from Bieging et al. (2016) shows that two other H ii regions from G174+2.5, Sh2-231, and Sh2-232, apparently already dispersed their front and back neutral walls. Thus, all four H ii regions in G174+2.5 that are (partly) embedded into molecular cloud have less material on the front than on the rear side.

Besides geometry, our chemical model contains several assumptions that might affect the [¹³C ii] and [C ii] pv diagrams. First, the ’high-metallicity’ initial elemental abundances from Wakelam & Herbst (2008), used in this work for the total gas-phase and grain mantle carbon budget, might overestimate the amount of carbon locked in refractory grain cores and, consequently, underestimate the amount of carbon available for chemical reactions in the model. The carbon abundance, used as a reference in this elemental abundance set, is based on observations of the [C ii] line at λ = 2325 Å towards ζ Oph: x(C) = 1.32 ± 0.32 × 10⁻⁴, while similar observations of χ Per give x(C) = 2.45 × 10⁻⁴ (Cardelli et al. 1993; Savage & Sembach 1996). Sofia et al. (2004) estimated x(C) = 1.62 ± 0.39 × 10⁻⁴ for ζ Oph from [C ii] absorption measurements, and concluded that the diffuse line-of-sight gas-phase value of x(C) might be as high as 2 × 10⁻⁴. Both stars are situated within 250 pc of the Sun, while the closest H ii region, the Orion Nebula, is at about 400 pc (Menten et al. 2007; Großschedl et al. 2019; Kounkel & Covey 2019; Kuhn et al. 2019). The value of x(C) measured towards ζ Oph is successfully used in astrochemical models of the Orion region (e. g. recent work by Cuadrado et al. 2017). However, Rubin et al. (1991) and Osterbrock, Tran & Veilleux (1992) found higher values of x(C) up to x(C) = 2.7 × 10⁻⁴ in the Orion Nebula. For the S106 region, which is about 1.3 kpc from the Sun, Schneider et al. (2018) used x(C) = 2.34 × 10⁻⁴ based also on analysis of the Orion region by Simón-Díaz & Stasińska (2011). A gas-phase elemental abundance x(C) = 8 × 10⁻⁵, significantly lower than the values mentioned above, was found by Fuente et al. (2019) for the nearby dark cloud TMC 1. We see that variations of the gas-phase x(C) in the simulations can be justified observationally in broad limits, and can be tuned to improve the agreement between the observations and simulations. Using higher elemental carbon abundance in our calculations with the ’ISM’ set from cloudy (Ferland et al. 2013) allows us to reproduce the [C ii] double-peaked line profile in the model for S235 C without addition of the external layer of cold C⁺; see Appendix D (supplementary material available online).

Another uncertainty in the initial conditions of the model lies in choosing low gas and dust temperature, and also fixing all carbon in the form of solid CO. Some carbon in the grey positions in Fig. 6 might indeed be hidden from the gas phase on the grain surfaces due to the low dust temperature. Therefore, our choice of solid CO as initial conditions for the model is, at least partially, supported by the observations. Kirsanova et al. (2014) determined the value of T_gas in the direction of S235 A and the young stellar cluster south of it. They found values up to 60 K in S235 A and 20 ≤ T_gas ≤ 32 K towards the cluster. Similar values were recently found by Burns et al. (2019). Kirsanova et al. (2014) also found a gas density higher than 10⁴ cm⁻³ in the molecular cloud around S235 A and in the direction of the embedded stellar cluster. The value of T_gas is higher in the direction of the young stellar cluster than T_dust by up to 10–20 K, and the difference is even higher in S235 A, despite the high density. A similar situation was found, for example, in the direction of S140 PDR by Koumpia et al. (2015). Our model of the expanding structure ‘H ii region + PDR + molecular cloud’ for S235 A can explain such a discrepancy between the T_gas and T_dust values for the gas with density ≥10⁴ cm⁻³; see Fig. 8. We note that various methods for the temperature determination have their own specific features (e.g. Schneider et al. 2016).

As we have shown above, an additional external layer of cold C⁺ would produce the self-absorption gap in the model, but we do not try to reproduce this layer using, for example, external UV radiation in the model, because it requires additional parameters. The [C ii] emission in our models for S235 A and S235 C appears from the PDRs (in agreement with previous theoretical studies by Kaufman, Wolfire & Hollenbach 2006), but the self absorption can be created in the undisturbed gas with lower T_ex, the density and temperature of which represent the cold neutral medium (CNM; e.g. Wolfire et al. 2003). The CNM in the external envelope of the S235 C model may consist of C⁺ and H₂ instead of atomic hydrogen, as suggested in recent theoretical calculations of the thermal and chemical steady state by Bialy & Sternberg (2019). We do not start from steady state abundances, although we emphasize that different initial conditions might change the pv diagrams considered in this study, i.e. increase or decrease the [C ii] self-absorption.

8 CONCLUSIONS

In this work, we presented observations of the PDRs around the compact H ii regions S235 A and S235 C in the [C ii], [¹³C ii] lines at 158 μm, [O i] at 63 μm, and HCO⁺(3–2) at 267 GHz, tracing different layers around the volume of the ionized hydrogen. We summarize our main results as follows:

• The double-peaked [C ii] and [O i] line profiles are due to high optical depth (⁠|$\tau _{\rm C^+}$| is up to 10 in several directions in S235 A), not directly tracing the expansion of the front and back walls of the PDRs. The foreground absorbing material has V_lsr ≈ −16 km s⁻¹, which remains almost unchanged over the observed area. We find an expanding motion of the [C ii]-emitting layer into the molecular layer with a velocity up to 1 km s⁻¹ in both PDRs comparing velocities of the [¹³C ii] and HCO⁺(3–2) lines.

• Gas components neither visible in the [C ii] lines, nor in low-J CO lines, may contribute to the total column in several points across S235 A. Due to the low dust temperature in those directions, we expect CO depletion due to the CO freeze out on to dust grains. Apparently, S235 A represents a good example of an H ii region surrounded by a PDR, molecular cloud, and a snow line, where molecules are absorbed from the gas phase into the dust grains. The H ii regions are situated on the border between the dense molecular and more diffuse atomic gas.

• We are able to reproduce the physical parameters of the H ii regions and integrated intensities of the [¹³C ii], [C ii], and [O i] lines from the PDRs with a spherically symmetric chemo-dynamical model as these properties are defined by spectral type of the ionizing stars and initial gas number density. However, the model does not reproduce the double-peaked [C ii] line profiles related with the self-absorption and high optical depth of the line. To be in agreement with the observations, we need to enhance C⁺ column density outside of the PDRs. We propose that the PDRs are surrounded by cold neutral medium composed of evaporating globules. While the [C ii] line emission originates in the model PDRs, the cold neutral medium is responsible for the self-absorption features in the line profiles.

• The [O i] line profiles are the best tracers of the expanding motion in the considered PDR models in comparison with the [C ii], [¹³C ii], and HCO⁺(3–2) lines. However, the absorbing cold neutral medium might distort the [O i] line profiles originating from the PDRs and contaminate the observed profiles.

SUPPORTING INFORMATION

Figure S1. Simulated pv diagrams in the model with initial gas number density |$n_{\rm init} = 5 \times 10^4\, {\rm cm}^{-3}$| and T_eff = 31 000, 29 000, 27 000, and 25 000 K.

Figure S2. Simulated pv diagrams in the model with T_eff = 29 000 and initial gas number density |$n_{\rm init} = 10^5, 5 \times 10^4,\, 10^4, 5 \times 10^3\ {\rm cm}^{-3}$|⁠.

Figure S3. Top: specific emissivity per column of HCO⁺ molecules for transitions with J ≤ 4.

Figure S4. Simulated pv diagrams in the model with elemental abundances taken from the ‘ISM’ set form cloudy.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

We thank the anonymous referee for critique and suggestions that led to the improvement of this paper.

This work is based on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy (SOFIA). SOFIA is jointly operated by the Universities Space Research Association, Inc. (USRA), under NASA contract NAS2-97001, and the Deutsches SOFIA Institut (DSI) under DLR contract 50 OK 0901 to the University of Stuttgart.

MSK and PAB were supported by Russian Science Foundation, research project 18-72-10132, during their work with the observational data. Theoretical simulations by MSK and YNP were supported by the Russian Foundation for Basic Research, research project 18-32-20049. AMS was supported by the Ministry of Science and Education, FEUZ-2020-0030.

VO was supported by the Collaborative Research Centre 956, sub-project C1, funded by the Deutsche Forschungsgemeinschaft (DFG), project ID 184018867. NS acknowledges support by the Agence National de Recherche (ANR/France) and the Deutsche Forschungsgemeinschaft (DFG/Germany) through the project ‘GENESIS’ (ANR-16-CE92-0035- 01/DFG1591/2-1) and from the BMBF, Projekt Number 50OR1714 (MOBS – MOdellierung von Beobachtungsdaten SOFIA).

The James Clerk Maxwell Telescope is operated by the East Asian Observatory on behalf of The National Astronomical Observatory of Japan, Academia Sinica Institute of Astronomy and Astrophysics; the Korea Astronomy and Space Science Institute, and Center for Astronomical Mega-Science (as well as the National Key R&D Program of China with No. 2017YFA0402700). Additional funding support is provided by the Science and Technology Facilities Council of the United Kingdom and participating universities in the United Kingdom and Canada. We used the data from SASSY project (Program ID MJLSY02).

This research has made use of NASA’s Astrophysics Data System Bibliographic Services; SIMBAD data base, operated at CDS, Strasbourg, France (Wenger et al. 2000); Aladin web page (Bonnarel et al. 2000); astropy, a community-developed core python package for astronomy (Astropy Collaboration et al. 2013); and aplpy, an open-source plotting package for python (http://aplpy.github.com; Robitaille & Bressert 2012).

DATA AVAILABILITY

The data underlying this paper will be shared on reasonable request to the corresponding author.

Footnotes

REFERENCES