© The Authors 2026

1 Introduction

High-mass stars (>8 M_⊙) are fundamental drivers of galactic evolution and dynamics through their strong interaction with the surrounding interstellar medium (ISM). These stars play a crucial role in regulating the Galaxy’s energy balance through feedback mechanisms such as strong stellar winds, the intense ionizing ultraviolet (UV) radiation throughout their lifetimes, and the powerful supernova explosions at the end of their life cycles. The UV photons emitted by high-mass stars have enough energy to ionize the surrounding neutral hydrogen, resulting in the formation of H II regions. Because these stars have short lifetimes – only a few million years (see Motte et al. 2018, and reference therein) – and are deeply embedded in their natal molecular clouds, H II regions observed via their radio free-free emission or radio recombination lines (RRLs) serve as strong indicators of recent high-mass star formation (HMSF) activity in the Galaxy.

As H II regions evolve, they expand under the combined influence of internal pressure, stellar radiation, and winds from central OB stars. This expansion can compress the surrounding molecular gas, potentially triggering secondary star formation along the ionization fronts (Deharveng et al. 2010; Elmegreen 2011; Thompson et al. 2012). Young embedded H II regions profoundly affect their environments through feedback mechanisms such as shocks, outflows, and radiation, thereby shaping the conditions for subsequent star formation (Elmegreen & Lada 1977; Deharveng et al. 2003; Elmegreen 2011). However, the extent to which these regions modify the initial conditions or trigger new episodes of HMSF within their parental molecular clouds remains uncertain. A comprehensive understanding of these processes requires accurate measurements of the molecular gas temperature and density near expanding H II regions, along with the physical characteristics of the associated ionized gas and the relationships between them.

Previous extensive radio continuum and recombination-line surveys have unveiled a large number of Galactic H II regions (e.g., Reifenstein et al. 1970; Altenhoff et al. 1979; Lockman 1989; Kuchar & Clark 1997; Urquhart et al. 2007, 2009, 2013; Anderson et al. 2011, 2014; Kalcheva et al. 2018; Gao et al. 2019). Recently, the GLObal view on STAR Formation in the Milky Way survey (GLOSTAR¹; Medina et al. 2019; Brunthaler et al. 2021) cataloged a large population of new and known H II regions using highly sensitive Karl G. Jansky Very Large Array (VLA) continuum data in the D (Medina et al. 2019, 2024) and B configuration (Dzib et al. 2023; Yang et al. 2023) at spatial resolutions of 18″ and 1″, respectively. Khan et al. (2024, hereafter Paper I), compiled a catalog of Galactic H II regions detected with GLOSTAR RRLs data observed with the VLA in the D configuration and analyzed their ionized gas properties. This study revealed that approximately 50% of these H II regions are associated with sources identified in the Atacama Pathfinder EXperiment (APEX) Telescope Large Area Survey of the GALaxy (ATLASGAL; Schuller et al. 2009), which provided a comprehensive view of HMSF at 870 μm – a wavelength particularly sensitive to cold dust emission tracing dense star-forming clumps. While the ATLASGAL survey determined properties of the dust, the gas density and temperature of the clump play a crucial role in regulating their chemistry and star formation activity as well as in influencing the stellar initial mass function. A precise measurement of the dense gas density and temperature in the vicinity of the H II region are essential to deepening our understanding of how H II regions impact their parental molecular clouds and surrounding environments.

H₂CO, which has been ubiquitously detected in the ISM, is primarily formed via the hydrogenation of CO on dust grain surfaces (Simons et al. 2020) and is subsequently released into the gas phase through processes such as UV irradiation or shocks. As a result, it is commonly associated with star-forming and H II regions (Downes et al. 1980; Henkel et al. 1983; Du et al. 2011; Ginsburg et al. 2011; Gong et al. 2023), where its abundance remains relatively stable throughout different stages of the star formation cycle (Mangum et al. 1990; Tang et al. 2018b). Due to its large dipole moment (2.33 Debye; Fabricant et al. 1977), H₂CO exhibits rotational transitions across the centimeter and submillimeter wavelength range (Mangum & Wootten 1993). Spectroscopically, the rotational states of this asymmetric and prolate molecule are described by the total angular momentum and its projection quantum numbers, J_KaKc. Furthermore, due to the nuclear spin of the hydrogen nuclei, H₂CO exists in two spin isomers: ortho (K_a = odd) and para (K_a = even). However, one of the most valuable properties of H₂CO for astrophysical studies is its characteristic “K-doublet” splitting, which arises from the coupling between vibrational and rotational angular momentum states. In its rotational spectrum, this results in closely spaced pairs of spectral lines (or “K-doublets”) with opposite parity, rather than a single transition.

With advancements in receiver technology and increased access to submillimeter windows, astronomers have expanded the utility of H₂CO as a diagnostic tool for characterizing the physical conditions in molecular gas, even using its higher lying submillimeter transitions (Tang et al. 2018a,b, 2021, and references therein). These properties make H₂CO one of the very few molecular species that act as both a thermometer and densitometer for dense gas, particularly in environments where other tracers such as NH₃ or CH₃CN may be less accessible. For these reasons observations of the “K-doublet” transitions of H₂CO serve as a great tool to study the physical properties of the molecular medium and allow observations with a variety of radio telescopes.

To summarize, detailed study of a large sample of embedded H II regions and surrounding molecular clouds is required to probe the underlying physical processes, the physical properties of the ambient dense molecular cloud into which they are expanding, and the effect of HMSF on the parent molecular clouds. To measure the kinetic temperature and gas density of massive star-forming clumps identified in the GLOSTAR H II region and ATLASGAL surveys, here we analyze the rotational transitions of p–H₂CO (J = 3–2 and 4–3). Using these results, we aim to (a) explore possible correlations between the characteristics of molecular gas, ionized gas, and dust in these regions; (b) examine the temperature and density structure of molecular gas surrounding H II regions; and (c) investigate whether the pressure of the surrounding medium is sufficient to confine H II regions. This paper is structured as follows: In Sect. 2, we describe the source sample studied, the selection criteria we used, the H₂CO observations, and data reduction process. In Sect. 3, we outline the methods used to analyze and model the observed spectral line profiles of the p–H₂CO transitions studied. Section 4 presents the results, and this is followed by a discussion in Sect. 5. Finally the main conclusions are summarized in Sect. 6.

2 Sample, observations, and data reduction

2.1 Source selection

We identified 244 Galactic H II regions in Paper I using detected stacked RRLs (H98α–H114α) within the GLOSTAR survey coverage (−2° ≤ ℓ ≤ 60° & |b| ≤ 1° and 76° ≤ ℓ ≤ 83° & −1° ≤ b ≤ 2°). These authors derived physical properties of identified H II regions such as the electron temperature, electron density, and emission measure to characterize the associated ionized gas. By comparing the H II region catalog with dust continuum data, it was found that nearly half of the identified H II regions are spatially associated with dust clumps from the ATLASGAL survey. Urquhart et al. (2018, 2022) previously classified these clumps as H II regions and provided key dust properties such as the bolometric luminosities (L_bol), the clump masses (M_clump), and the dust temperatures (T_dust). Such a combined dataset offers insights into both the ionized gas and dust components associated with H II regions. This study extends the analysis in Paper I by investigating in addition to the ionized gas and dust, the properties of the molecular gas surrounding H II regions and exploring their interplay.

The targets for this study were selected based on the following criteria: (1) the ionized gas properties of the H II regions have already been derived (as cataloged in Paper I); (2) the target has an association with 870 μm continuum emission taken from ATLASGAL catalog, with dust properties derived in Urquhart et al. (2018, 2022); (3) the deconvolved full width half maximum (FWHM) of the H II region (from Paper I) is comparable to or smaller than that of the corresponding ATLASGAL clump; and (4) the target has well-characterized kinematics. We excluded H II regions located toward the Galactic center (−2° ≤ ℓ ≤ 2°) and the Cygnus X complex due to significant uncertainties in their kinematic distances.

These criteria ensure that the targeted H II regions are likely to remain embedded within their natal molecular clouds, enabling an investigation of the physical link between the ionized and molecular components during massive star formation. Based on these criteria, we selected 61 H II regions from Paper I for the p–H₂CO spectral line observations. Figure A.1 presents three–color mid–infrared composite images (8,24, and 70 μm) with contours marking the GLOSTAR 5 GHz radio continuum emission, highlighting the selected H II regions in this study.

Among them, only G30.720–0.082 is unresolved in the GLOSTAR RRL data at 25″ resolution (see Fig. A.2). Six sources (G7.177+0.088, G10.440+0.011, G24.470+0.490, G25.159+0.060, G43.169+0.008, and G43.149+0.013) exhibit extended radio emission surrounding compact bright cores, resulting in sizes larger than their associated ATLASGAL clumps (see Fig. A.2). We retained these sources in the sample to preserve diversity in the population. The selected regions span Galactic longitudes from 3° ≤ ℓ ≤ 53°, encompassing several prominent high-mass star-forming complexes, including W31, W33, W43, and W51.

Figure 1 presents the Galactic distribution of the sources in this study, overlaid on an artist’s impression of the Milky Way². The sources span distances from 1.8 to 20.4 kpc, with values for each H II region and its associated clump taken from Paper I (and references therein). The distribution reveals two distinct populations: one clustered near 4 kpc and another around 12 kpc. Figure A.3 shows the distributions of the angular and physical sizes of the selected H II regions. The angular sizes are concentrated, with a median of 15.6", reflecting our selection criterion based on source angular size. The physical size distribution, with a median of 0.52 pc, is more scattered, due to scatter in distance. This scatter in distance introduces no significant bias in the H II region properties but rather ensures a diverse and representative sample in terms of size and spectral type of central OB stars of H II region as discussed in Appendix B.

Fig. 1 Galactic distribution of the source sample. The locations of the targets are denoted by red stars on an artistic rendering of the Milky Way. The distances to the observed H II regions and their associated clumps are taken from Paper I (and references therein). The spiral arms are marked following the description presented in Reid et al. (2019), and the individual spiral arms are colored as follows: the 3 kpc in magenta, the Norma-Outer Arm in red, the Scutum-Centaurus Arm in blue, the Sagittarius-Carina Arm in green, the Local Arm in cyan, and the Perseus Arm is in yellow.

2.2 APEX 12 m observations of H₂CO

The single pointing observations were carried out from March to April 2023 using the Atacama Pathfinder EXperiment (APEX; Güsten et al. 2006)³ 12 m submillimeter telescope located in the Chajnantor plateau, Chile. Table 1 summarizes the key spectroscopic parameters of the observed p–H₂CO transitions alongside the receiver setups used and the detection rates within the sample. All three transitions of p–H₂CO within the J=3–2 states near 218 GHz were observed using the dual sideband nFLASH230 receiver. At this frequency, the observations were carried out with a half-power beam width of ~27″. In addition to these p–H₂CO transitions, the receiver was tuned such that it also covered the H30α line at 231 GHz in its upper sideband. Total integration times between 7 min and 16 min were spent in this setup. The nFLASH230 receiver was connected to a MPIfR-built high-capacity fast Fourier transform spectrometer (FFTS; Klein et al. 2012) backend with two sidebands. Each sideband has two spectral windows with 4 GHz bandwidth, providing both orthogonal polarizations and resulting in a total bandwidth of 8 GHz. This provided a spectral resolution of 61 kHz, resulting in a velocity resolution of 0.08 km s⁻¹ at 218 GHz.

The three transitions of p–H₂CO within the J=4–3 levels were observed with the SEPIA345 receiver (Meledin et al. 2022) with a beam size of ~20″ and integration times between 8 min to 20 mins. Similar to the nFLASH230 receiver, the SEPIA345 receiver is a dual polarization, dual sideband receiver which covers a total frequency range from 272 GHz to 376 GHz. Each sideband (and polarization) of the SEPIA345 receiver output is recorded using two FFTS spectrometer units, each of them sampling 4 GHz, with a spectral resolution of up to 64 k channels per sideband, this results in velocity resolutions of 0.04 km s⁻¹. Both sets of observations were performed in position-switching mode with off positions offset from the on position of the sources by (±600″, ±600″) in (ℓ, b).

The data reduction was performed using the CLASS software from the GILDAS package⁴. The spectra were converted from antenna temperature (T_A) scales to main-beam brightness temperature (T_MB) scales using a factor of η_F/η_MB, where η_F is the forward efficiency and η_MB⁵ is the telescope main beam efficiency. In our observations, η_F was set to 0.95 and η_MB = 0.79 ± 0.07 for the data collected using both nFLASH230 and SEPIA354 (based on observations of Mars). To enhance the signal-to-noise ratios (S/N) in individual channels for both sets of observations, the data were smoothed to a uniform channel spacing of 0.6 km s⁻¹. The typical line width of H₂CO transitions toward H II regions tends to be on the order of a few kilometer per second, so the channel smoothing has no adverse impact on our results.

3 Analysis

In this section, we describe the multicomponent Gaussian fitting of the p–H₂CO spectral lines observed in this study. We then use the resulting line-fitting parameters to constrain the non-local thermodynamic equilibrium (non-LTE) analysis and derive the kinetic temperature (T_kin), gas density n(H₂), and the p–H₂CO column density N(p–H₂CO).

3.1 Spectral line fitting

We modeled the observed spectral line features using Gaussian profiles that were fit using the MODEL function within the LMFIT⁶ package in Python. As output, it provided the velocity-integrated line intensity (I = ∫ T_mbdv), the local standard of rest (LSR) velocity (v_LSR), the FWHM of the Gaussian line, and the peak brightness temperature on main-beam temperature scales (T_peak).

We began our analysis by fitting single-component Gaussian to the p–H₂CO spectra. However, this approach yields well constrained physical parameters for only 25 out of 61 targets (40%). To improve the modeling, we performed multicomponent Gaussian fits to better model the observed p–H₂CO line profiles. As the brightest transition, p–H₂CO (3_0,3–2_0,2) was detected toward all sources and served as the reference for identifying the number of velocity components and their centroid velocities. These centroid velocities were then fixed for all other p–H₂CO transitions to ensure consistency across lines. Each p–H₂CO (3_0,3–2_0,2) spectrum was initially fit with one Gaussian component, if the residuals are larger than 3σ, a second or even a third component was added, until the standard deviation in the residual is minimized. We verified all fits interactively using the LINE/MINIMIZE task in CLASS/GILDAS. To verify that the H₂CO lines exhibit multiple components and are not affected by optical depth effects such as self-absorption, we compared their spectra with those of the optically thin C¹⁸O line. We did not find any evidence of self-absorption in our sample. We found no preferential occurrence of multiple components in sources located at either near or far distances.

Among the 61 observed sources, the p–H₂CO spectra toward 34 of them displayed multiple Gaussian components, resulting in a total of 105 identified components in the entire sample. Since all six transitions are expected to trace gas with similar kinetic temperature and dynamic motion, we expect them to have similar line widths. However, this need not always be true especially when fitting additional components arising from surrounding envelopes or outflows. Therefore, the line widths were allowed to vary within 1 σ of the value constrained by the line width of the 3_0,3–2_0,2 line of p–H₂CO, whereas the velocity centroids of the components were fixed. Any component with peak brightness temperatures below 3σ level of the rms noise at 0.6 km s⁻¹ or line width less than 3×Δv (1.8 km s⁻¹), is flagged before further analysis. An example of the p–H₂CO multicomponent Gaussian fitting is shown in Fig. 2 (left panel) toward one of the targets in our study, AGAL045.121+00.131 (also listed as G045.122+0.131 in the GLOSTAR H II region catalog).

If a source exhibits multiple velocity components, then the component associated with the source at a p–H₂CO velocity close to the systemic velocity of the corresponding ATLASGAL clump – determined from multiple molecular line surveys and reported by Urquhart et al. (2018, see their Table 1) – is considered the main component and is hereafter referred to as “H₂CO Main.” All secondary components are then considered as sub-components and is hereafter referred to as “H₂CO Sub.” These secondary components may either be associated with the source itself via outflows/inflows, envelops, or may arise from distinct molecular clumps covered in our beam. However, owing to the limited spatial resolution, the exact origin of the “H₂CO Sub” components remains uncertain.

The derived line parameters for p–H₂CO are then used as input to constrain the column densities and the gas densities. However, as we observed two different sets of transitions (J = 3–2 and J = 4–3) at different frequencies (218 GHz and 291 GHz) using two different receivers (nFLASH230 and SEPIA345), with varying HPBWs 27.6″ and 20.7″, the observations resulted in slightly different spatial coverages. Therefore, to ensure uniformity, we further correct for beam filling effects. Previous studies of H₂CO (Immer et al. 2014; Tang et al. 2018a,b) have demonstrated a tight correlation between the integrated intensities of the p–H₂CO transitions and flux density of the dust emission at 870 μm. This suggests that the dense gas traced by H₂CO is associated with dust emission at 870 μm. Based on these findings, we assume that the size of the H₂CO emitting area toward the targets of this study are equivalent to the FWHM source sizes of 870 μm dust emission, as derived by Csengeri et al. (2014). Following Tang et al. (2018b), we corrected for beam dilution as follows, T_MB/η_bf where η_bf is the beam filling factor, given by, ${\eta }_{\mathrm{b}\mathrm{f}}={\theta }_{\mathrm{S}}^2/({\theta }_{\mathrm{S}}^2+{\theta }_{\text{beam}}^2)$. Here, θ_S and θ_beam denote the source and beam size, respectively. The corrected p–H₂CO line intensities with the respective beam filling factors used and other line fitting parameters are summarized in Tables C.1 and C.2.

Fig. 2 Spectra of the different p–H2CO and RRL transitions studied in this work toward AGAL045.121+00.131. Left panels: spectra of all six of the observed p–H2CO transitions (black) along with the multicomponent Gaussian fitting (green and magenta) and residual of the fitting (dashed blue). The total fit is displayed in red and the gray shaded region demarcates the 3σ level of the rms noise. The vertical dashed gray lines mark the centroid velocities of each fit Gaussian component. Right panel: spectra of the stacked GLOSTAR ⟨Hnα⟩ RRL in black alongside the Gaussian fit in red. Also displayed is the 5.3 GHz radio continuum emission map obtained with the GLOSTAR VLA D configuration toward this source and overlaid with white contours is the ATLASGAL 870 μm. The white filled circle shows the 25″ GLOSTAR beam.

3.2 Non-LTE modeling with Pyradex

We carried out non-LTE radiative transfer calculations to derive the physical properties of the clumps, including the gas kinetic temperature, T_kin, spatial density, n(H₂), and column density of p–H₂CO, N(p–H₂CO). We employed the Python wrapper, PyRADEX, to interface with the non-LTE radiative transfer code RADEX (van der Tak et al. 2007), in combination with the emcee⁷ (Foreman-Mackey et al. 2019) package in Python, which implements a Markov chain Monte Carlo (MCMC) algorithm as shown in Yang et al. (2017) and Christensen et al. (in prep.). To reproduce the observed integrated line intensities of the varied p–H₂CO transitions (as determined in the previous Sect. 3.1), the PyRADEX models were run with constraints from the observed line widths using collisional rate coefficients computed for collisions between p–H₂CO and ortho-, and para-H₂ (Wiesenfeld & Faure 2013) and assuming the ortho-to-para ratio (OPR) of 3. For each velocity component toward a given source, the line width was fixed to the average value across the different p–H₂CO transitions probed. This approach minimizes variations and allows us to simultaneously reproduce the physical conditions traced by all six transitions – a valid constraint, since all transitions are expected to originate from the same gas. We further assumed a background temperature, T_bkg, of 2.73 K. We initialized the MCMC process using the curve_fit function from the SciPy⁸ package in Python, then explored the parameter space within n(H₂) = 10–10¹⁰ cm⁻³, T_kin = 10–300 K and N(p–H₂CO) = 10¹⁰–10¹⁷ cm⁻². For each component, we used 400 walkers, discarding the first 100 steps as burn-in and allowing the final 1000 steps to converge on the best-fit integrated line intensity (see also Christensen et al., in prep.). Using this method, we were able to successfully constrain the physical parameters toward 83 of the 105 targets studied in this work. Figure A.4 illustrates the results of the PyRADEX and MCMC computation, displaying the results in the form of a corner plot with 1D histograms of the posterior distributions across the explored parameter space. Table C.3 lists the derived gas kinetic temperatures, volume densities, and column densities of p–H₂CO. Figure 3 displays the distribution of the sources in the n(H₂) − T_kin parameter space. The “H₂CO Main” components of sources with single and multiple Gaussian profiles show a similar distribution, while the extreme values are associated with the “H₂CO Sub” components. The only exceptions are G31.412+0.308, which exhibits the highest T_kin, and G10.462+0.033, which shows the highest n(H₂).

Fig. 3 Distribution of the sources in the n(H2)–Tkin plane. The blue and red colors respectively represent sources with single and multiple component Gaussian profiles. Filled circles denote components classified as “H2CO Main,” while stars indicate those identified as “H2CO Sub.”

4 Results

In the following section, we present the results obtained from the joint PyRADEX and MCMC analyses described in Sect. 3.2. These results include the derived kinetic temperature, gas density, p–H₂CO column density, and the fractional abundance of p–H₂CO, X(p–H₂CO), which we compare with values from previous studies. We also examine the optical depth of p–H₂CO and assess its impact on our analysis (Appendix D).

Figure 4 shows the distributions of gas kinetic temperatures, gas densities, and p–H₂CO column densities, jointly constrained from the PyRADEX and MCMC analyses discussed above while the results toward each source are listed in Table C.3. In the following sections, we discuss these results in detail.

4.1 Gas densities

As shown by the distributions in Fig. 4 for our sample, the derived n(H₂) values range between 8.32 × 10³ cm⁻³ and 1.05 × 10⁷ cm⁻³ with a median⁹ value of (3.7 ± 0.35) × 10⁵ cm⁻³. Our results are consistent with previous studies of starforming regions (Mangum & Wootten 1993; McCauley et al. 2011; Mangum et al. 2013), including the brightest ATLASGAL clumps associated with H II regions, where densities of 6.4 × 10⁵–8.1 × 10⁶ cm⁻³ were reported, with a median value of 1.25 × 10⁵ cm⁻³ (Tang et al. 2018b). Furthermore, a comparison of our results with n(H₂) values determined from 870 μm dust emission (Urquhart et al. 2022) reveals systematically higher densities within our sample (Fig. 5), indicating that these p–H₂CO transitions trace denser gas than the dust continuum. This holds true for both the systemic velocity components and additional components arising from the envelopes or other expanding shells of these H II regions. Figure 4 a shows that the median gas density of the “H₂CO Sub” component similar to that of the “H₂CO Main.”

4.2 Kinetic temperature

As with n(H₂) and N(p–H₂CO), the PyRADEX+MCMC models yielded a range of gas temperatures across the sources in our sample, from 33.7 K to >265 K¹⁰, with a median value of 66 ± 3 K. In contrast, the study by Tang et al. (2018b), which examined ATLASGAL dust clumps at various evolutionary stages, reports a median kinetic temperature of 67 K across their entire sample, increasing to 104 K for sources associated with only H II regions. To assess whether the difference in kinetic temperature derived for the H II region sample by Tang et al. (2018b) and this work are statistically significant, we applied the Anderson-Darling (AD) test to the two distributions. The resulting p-value of approximately 0.001 indicates a statistical difference between the kinetic temperature distributions, providing evidence that T_kin do not arise from the same parent population. This discrepancy in kinetic temperatures may be due to the fact that we decompose the observed spectral features into multiple components likely originating from different environments, while Tang et al. (2018b) carry out their analysis by considering only the integrated line intensities across the entire observed features. The latter may lead to higher line intensity ratios in the analysis conducted by Tang et al. (2018b), resulting in higher overall gas temperatures. Furthermore, we find the values of T_kin derived for the “H₂CO Sub” components (98.4 K) to be higher than that for the “H₂CO Main” components (57.8 K). This suggests that the “H₂CO Sub” components may originate from inflowing/outflowing gas surrounding the star, or from shock-heated molecular gas. Moreover, the lack of correlation between the kinetic temperature of the “H₂CO Sub” component and its velocity difference with both the “H₂CO Main” component and the RRL emission further supports a shock-heated origin for the “H₂CO Sub” components.

The T_kin values derived from the PyRADEX+MCMC analysis can be verified using the integrated intensity ratios between specific p–H₂CO transitions, namely, R₁ = I(3₂₁–2₂₀)/I(3₀₃–2₀₂), and R₂ = I(4₂₂–3₂₁)/I(4₀₄–3₀₃) – which provide an estimate of the T_kin assuming conditions of LTE. This method is applicable if the observed lines are optically thin (discussed further in Appendix D) and originate from high-density regions. Following the approach outlined in Appendix A of Mangum & Wootten (1993), we computed the T_kin using these line ratios: $$T_{\mathrm{L}\mathrm{T}\mathrm{E}}(R_1)=\frac{47.1}{\ln \left(0.556\frac{I(3_{03}-2_{02})}{I(3_{21}-2_{20})}\right)}\mathrm{K},$$(1)

and $$T_{\mathrm{L}\mathrm{T}\mathrm{E}}(R_2)=\frac{47.2}{\ln \left(0.750\frac{I(4_{04}-3_{03})}{I(4_{22}-3_{21})}\right)}\mathrm{K}.$$(2)

In Fig. 6, we refer the temperatures derived using Eq. (1) and Eq. (2) as LTE temperatures, and find the results from both line ratios to be consistent with one another. If our assumption that the p–H₂CO emission is optically thin is valid, then the kinetic temperatures obtained by this method have an uncertainty of ≲30% (Mangum & Wootten 1993). Comparisons between T_kin values derived from LTE and PyRADEX+MCMC non-LTE calculations show overall consistency, although some deviations appear at higher temperatures (Fig. 6). More interestingly, the derived values of the kinematic temperature are significantly less scattered for the “H₂CO Main” velocity components than the “H₂CO Sub” component. Overall, these findings reaffirm the diagnostic power of H₂CO an excellent molecular thermometer.

Fig. 4 Distributions of the pyRADEX+MCMC derived (a) n(H2), (b) Tkin, (c) N(p–H2CO), and (d) p–H2CO fractional abundance, X(p–H2CO) toward the different velocity components identified in our sample, as traced by the p–H2CO. The blue and red histograms represent the distributions for the “H2CO Main” and “H2CO Sub” components, respectively, while the corresponding dashed lines indicate the median values of each distribution.

4.3 Fractional abundance of p–H₂CO

In this section, we compute the fractional abundance of p–H₂CO with respect to H₂ (X(p–H₂CO) = N(p–H₂CO)/N(H₂)), where N(H₂) refers to the H₂ column density derived from the 870 μm continuum emission by Urquhart et al. (2018). The PyRADEX derived N(p–H₂CO) values lie between 4 × 10¹² cm⁻² and 5 × 10¹⁵ cm⁻² with a median value of 1.0 ± 0.1 × 10¹⁴ cm⁻², while N(H₂) was estimated from the peak flux densities of the clumps using $$N({\mathrm{H}}_2)=\frac{{\mathrm{D}}^2{F}_{870}R}{B_{870}(T_{\mathrm{d}}){\theta }_{\text{beam}}{\kappa }_{870}{\mu }_{{\mathrm{H}}_2}{m}_{\mathrm{H}}},$$(3)

where D is the distance to the source; F₈₇₀ is the peak flux density taken from Csengeri et al. (2014), θ_beam; the beam size of ATLASGAL survey (= ${19}_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}2$; Csengeri et al. 2014), μ_H2 is the mean molecular weight per hydrogen molecule assumed as 2.8, and m_H the mass of the hydrogen atom. Further, R is the gas-to-dust ratio assumed to be 100; B₈₇₀(T_d) is the intensity of the blackbody at 870 μm at the dust temperature, T_d, values for which are taken from Urquhart et al. (2018); and κ₈₇₀ is the dust opacity at 870 μm taken to be 1.85 cm² g⁻¹, an average of all dust models presented in Ossenkopf & Henning (1994). This yields N(H₂) values ranging from 6 × 10²⁰ cm⁻² to 8 × 10²² cm⁻² and X(p–H₂CO) values from 2.7 × 10⁻¹⁰ to 4.7 × 10⁻⁷ (see Fig. 4d), with a median value of 7.2 × 10⁻⁹. These results are consistent with values reported for other star-forming regions (Ao et al. 2013; Gerner et al. 2014; Tang et al. 2018a,b). Figure 4c shows that the median N(p–H₂CO) of the “H₂CO Sub” component is higher than that of the “H₂CO Main”, whereas the median X(p–H₂CO) is similar for both components (see Fig. 4d). This suggests that H₂CO can be at least as abundant – if not more so – in the envelope/outflow components, underscoring its presence across diverse environments. However, the “H₂CO Sub” components may actually represent unresolved condensates, since their densities match those of the “H₂CO Main” components, suggesting that the secondary features likely comprise of multiple dense clouds at closely spaced velocities unresolved in our single dish observations.

Fig. 5 Comparison between n(H2) probed by 870 μm emission from dust and p–H2CO as derived from our pyRADEX+MCMC models. The dashed black line displays a slope of unit value.

5 Discussion

5.1 Correlation between physical properties of ionized gas, molecular gas, and dust

Variations in the physical properties of clumps and H II regions associated with star formation offer valuable insights into the evolutionary stages and conditions of HMSF and its impact on the surrounding medium. For instance, over the years, the ratio of bolometric luminosity to clump mass has proven to be a reliable tracer of HMSF stages and has been widely used in numerous studies (Molinari et al. 2008, 2016; Giannetti et al. 2017; Urquhart et al. 2022). Similarly, n_e and RRL line widths have been extensively employed to characterize the evolutionary stages of H II regions (e.g., Sewilo et al. 2004; Kurtz 2005). In order to investigate the relationship between the ionized, molecular gas, and dust in H II regions, we compared the Spearman’s correlation between the different physical properties of molecular gas traced by p–H₂CO derived in the previous sections with those of dust, and ionized gas from Paper I. The resulting correlations are presented in Fig. 7. Note that we have flagged the cross correlation for which the p-value is > 0.0013.

No strong correlation was found between n(H₂) and the other ionized, molecular gas and dust properties (masked in Fig. 7). This suggests that, regardless of the specific characteristics of the HMSF clump, their molecular hydrogen densities typically remain around 10^5–6 cm⁻³, in alignment with the findings of Tang et al. (2018b).

Observations of massive clumps using other molecular line tracers such as CO, NH₃, CH₃CN, CH₃CCH, and CH₃OH, H₂CO (Wienen et al. 2012; Giannetti et al. 2017; Tang et al. 2018b) suggest that these clumps are heated by radiation from embedded massive stars. This implies that the T_kin values derived from these molecules should correlate with the L_bol of the clumps. Giannetti et al. (2017) reported such correlations using kinetic temperatures derived from CH₃CCH, and CH₃OH, with correlation coefficients of 0.68 and 0.54, respectively. Similarly, Tang et al. (2018b) investigated the same relationship using H₂CO and found a somewhat weaker correlation, with a coefficient of ~0.5. However, in this study, we find an even weaker correlation between T_kin and L_bol, with a correlation coefficient of 0.32 (masked in Fig. 7 due to p-value >0.0013).

Although the velocities of RRLs and p–H₂CO lines are generally consistent, there is still a velocity offset between the two, with a standard deviation in the velocity offset of ~4 km s⁻¹ and mean of −0.5 km s⁻¹. For an H II region of roughly 0.52 pc in size, the inferred dynamical age is t_dy ~ 1.3 × 10⁵ yr, which is comparable to the typical expansion timescales of H II regions. Hence, a clump with a velocity difference of 4 km s⁻¹ is consistent with a molecular layer being located near the ionization front of the expanding H II region. Assuming that this velocity difference serves as a proxy for the physical separation between the H II region and the molecular clump, the correlation between the physical parameters is recomputed for sources where this velocity difference is less than 4 km s⁻¹. For these closely associated sources, we find a stronger correlation between T_kin and L_bol, with a coefficient of 0.49 (p-value <0.0013). This result suggests that T_kin is more strongly influenced by bolometric luminosity when the molecular clump lies closer to the H II region.

The p–H₂CO column density correlates strongly with the p–H₂CO FWHM and more weakly with the clump bolometric luminosity. Similar but stronger trends were reported in massive star-forming regions by Immer et al. (2014) and Tang et al. (2018b). These correlations suggest that the p–H₂CO abundance increases with the radiation field. This would indicate that motions such as turbulence and shocks can contribute to the increase of H₂CO in the gas phase. A natural guess is that such motions help the desorption of H₂CO from dust grains.

We also find that N(p–H₂CO) correlates with the electron density of H II regions. The electron density depends on the ionizing photon flux from the central massive star(s), as higher UV photon output generates more free electrons. In addition, the electron density shows a strong correlation with the clump hydrogen column density, N(H₂). This indicates that young H II regions with higher electron densities tend to be formed in clumps with higher H₂ column densities. The electron temperature, T_e, shows only a weak correlation with n_e, with a Spearman’s correlation coefficient of 0.39.

Bolometric luminosity, L_bol, exhibits relatively strong correlations with dust, molecular gas, and ionized gas properties (Fig. 7). These results indicate that the central high-mass star significantly shapes the dynamics, chemistry, and physical conditions of its surrounding molecular and dusty environments. Furthermore, N_Lyc, which directly traces the ionizing capability of the central star, also correlates with L_bol.

Fig. 6 Comparison between the kinetic temperatures derived from pyRADEX+MCMC modeling and using the p–H2CO integrated line intensities ratios R1 (top panel) and R2 (middle panel). Bottom panel: comparison of the kinetic temperatures derived from R1 and R2. The LTE temperature uncertainties are from uncertainties in the integrated intensities of the observed p–H2CO lines. The dashed black lines indicate a slope of unity.

Fig. 7 Correlation matrix showing the relationships between various physical properties of molecular gas, ionized gas, and dust clumps. The color bar and annotations indicate the Spearman’s correlation coefficients. We flagged and excluded correlations with p-values > 0.0013 (3σ), thereby removing statistically insignificant relationships.

5.2 Benchmarking the H₂CO thermometer in H II regions

5.2.1 Comparing gas kinetic temperature and line width

The values of T_kin, derived using the non-LTE radiative transfer analysis presented in Sect. 4.2 are found to be consistently warmer than those derived from the dust emission. Additionally, 17 sources display blue- or red-skewed p–H₂CO line profiles, indicating the presence of outflow or infall motions and associated shocked gas. The observed line profiles and derived gas temperatures of H₂CO suggest that its formation is likely driven by shock heating or outflows and is located close to the H II regions within the parent molecular cloud (Ginsburg et al. 2011; Tang et al. 2017, 2018b).

In Fig. 8, we compare the gas kinetic temperature and line width obtained from the different gas tracers toward the H II regions by Giannetti et al. (2017). This shows that p–H₂CO traces a wide range of the gas temperature and the traced temperature is hotter than molecules such as NH₃. This is based on the findings of Wienen et al. (2012), who reported NH₃ rotational temperatures in the range of 10–28 K and kinetic temperatures between 12 and 35 K. However, these authors note that the low-Eu (NH₃ (1,1) and (2,2)) transitions used in their analysis exhibited different line widths, indicating that the observed NH₃ transitions likely do not trace the same gas. Consequently, assuming equal beam-filling factors for these transitions represents only an approximation and may affect the radiative analysis and the subsequently derived NH₃-bearing gas temperatures. This conclusion is further supported by Giannetti et al. (2017), who showed that these low-lying transitions of NH₃ are not an optimal tracers of warm gas. Similarly, p–H₂CO traces hotter regions compared to that traced by CH₃CCH. Giannetti et al. (2017) found two temperature components associated with the CH₃OH and CH₃CN namely, hot and cool components with an average of 21 and 184 K, and 53 and 238 K for H II regions in their sample. The kinetic temperature derived using p–H₂CO lies between the cool and hot component of CH₃OH and CH₃CN (see Fig. 8). Figure 8 (bottom) shows cumulative distribution function (CDF) comparisons of the molecular line width of different gas tracers. Although the Anderson-Darling test reveals that the distribution of line width of p–H₂CO and CH₃CCH, CH₃OH (cool component) shows no significant differences, that of the warm components show significant differences. The mean value of the line widths of p–H₂CO is greater than CH₃CCH and NH₃ and is less than CH₃OH (hot component) and CH₃CN (hot component). In contrast, the line width of p–H₂CO is similar to the line width of CH₃OH (hot component) and CH₃CN (hot component). If the observed line widths of different molecular tracers were determined purely by thermal motion, they would depend directly on the gas temperature and inversely on the molecular mass. In reality, however, these line widths include both thermal and nonthermal contributions, making direct comparisons between tracers more complex.

Fig. 8 Cumulative distribution functions for temperatures (top panel) and line width (bottom panel) for different H II regions. Different colors represent various tracers and components, as indicated in the legend. The NH3 properties are from Wienen et al. (2012). The CH3CCH, CH3OH (hot and cold), and CH3CN (hot and cold) properties are from Giannetti et al. (2017), whereas dust properties are from Urquhart et al. (2018).

5.2.2 Chemical modeling

In this section, we aim to compare the position of H₂CO relative to the different photodissociation region (PDR) tracers and dense gas tracers around H II regions. We used an updated version (v7.0) of the Meudon PDR code (Le Petit et al. 2006)¹¹ to compare our observations with PDR models and investigate the chemical complexity surrounding an H II region for a typical density-temperature structure. Using the Meudon code, we ran models of constant thermal pressure (isobaric) across a 1D PDR in slab geometry, with a fixed pressure of 4.6×10⁷ K cm⁻³, determined by the average n(H₂) and T_kin as traced by H₂CO toward our sample (see Sects. 4.1 and 4.2). Combining the given density and temperature profiles with external sources of energy: the UV radiation field (G₀) and cosmic-ray ionization rate (ζ), the code determines the abundance of each species by solving the chemical reactions under steady-state conditions. In addition to the external radiation fields described above, we included an internal stellar component—representing irradiation from host star(s)—and an isotropic ambient component. The isotropic field incorporates the interstellar radiation field with a value of 0.5 (in Habing units) on both the observer-facing and the backside of the cloud, as well as blackbody components that simulate dust infrared emission and the cosmic microwave background. A minimal external field is used such that illumination on the backside of the PDR is dominated by the stellar component. This isolates the cloud from external influences, focusing only on the effect of the embedded H II region on the surrounding molecular material.

We estimated the far-UV flux, G₀, assuming that a single massive star is responsible for generating the H II region, from the massive stars that dominate the ionizing photon production in the sources studied. In Paper I, we derived the Lyman continuum photon flux for our sample using integrated radio flux measurements, finding a median value of 1.18 × 10⁴⁸ s⁻¹. This flux corresponds to a stellar spectral type of O5.5V–O5V with an effective surface temperature of approximately 41000 K (Martins et al. 2005). We assumed that the distance between the star and the PDR boundary (d₀) is equal to the radius of H II region, which in our sample has a median value of 0.26 pc. This setup yields a stellar far-UV flux at the PDR surface of G₀ ~ 5.3 × 10⁴ in units of the Habing field (1.6 × 10⁻³ erg cm⁻² s⁻¹; Habing 1968). We note that this estimate of G₀ likely represents a lower limit, as the Lyman continuum flux derived from interferometric radio observations may underestimate the total radio flux, thus leading to an underestimation of the ionizing photon rate. In addition, absorption of UV photons by surrounding dust can further reduce the inferred value. Additionally, the relatively low spatial resolution (25″) of the RRL data used in Paper I could have led to an overestimated H II region size. To explore this further, we also ran the PDR models using a smaller star–PDR distance of d₀ = 0.12 pc, which is the 16^th percentile of our sample H II region size, this results in G₀ values of ~2.5 × 10⁵ in units of Habing field. Both sets of models were run for a fixed value of cosmic-ray ionization rate of H₂ molecules, ζ = 5 × 10⁻¹⁷ s⁻¹, that lies between estimations in diffuse and dense gas (Indriolo et al. 2015; Le Petit et al. 2004). We assumed a total visual extinction of A_max = 10 mag for the cloud and adopted the extinction curve of HD 38087 (Fitzpatrick & Massa 1990) with R_V = 5.6. For the modeling, we used an intermediate value for the column density to reddening ratio, N_H/E(B − V) = 1.05 × 10²² cm⁻² mag⁻¹ (Joblin et al. 2018) and a dust-to-gas mass ratio of 0.01.

A schematic representation of the modeled profile across the H II region associated PDR is presented in Fig. 9 while the results of the Meudon PDR code simulations for d₀ = 0.26 pc (left panel) and d₀ = 0.12 pc (right panel) are presented in Fig. A.5. The central and bottom panels plot normalized molecular line abundances, n(x)/n_peak(x), against A_V for PDR tracers (c − C₃H₂, C₂H, CO⁺, C⁺; Rizzo et al. 2005; Kim et al. 2020) and dense, molecular gas tracers (CN, CH₃OH, C¹⁸O, H₂CO; Kim et al. 2020), respectively, alongside the corresponding gas density (n_H⁺, n_H and n_H2) and T_kin profiles. Consistently between the two sets of models, the PDR tracers peak at low values of n(H₂) ~ 8 × 10³ to 5.6 × 10⁵ cm⁻³ (A_V = 2–5 mag, yellow shaded region in Fig. A.5) and T_kin ~ 900–70 K probing the dissociation front and translucent gas layers. However, as A_v and n(H₂) increase, the abundances of dense gas tracers such as C¹⁸O increase, peaking at A_V ~ 8 mag.

Notably, we observed both H₂CO and CH₃OH exhibit enhanced abundances near the PDR interface, where the gas temperature reaches a few hundred kelvin and the density is around 10^4–5 cm⁻³. This abundance, intermediate to that seen by typical PDR and dense gas tracers, is expected and provides insight into the chemistry of these two species. For instance, CH₃OH, the simplest complex organic molecule (COM; Herbst & van Dishoeck 2009) is formed efficiently in cold environments on grain surfaces by the continuous hydrogenation of CO (Fuchs et al. 2009) and by H-abstraction reactions of icy radicals (Álvarez-Barcia et al. 2018; Simons et al. 2020; Santos et al. 2022). CH₃OH is subsequently released into the gas phase via thermal sublimation or nonthermal desorption processes where the presence of elevated temperatures in the form of shocks, aids this release via sputtering. This explains the CH₃OH abundance peak seen at moderate A_V ’s and high temperatures, highlighting that the formation of CH₃OH only proceeds efficiently on grain surfaces and not in the gas phase. On the other hand, although H₂CO is an important precursor for the formation of CH₃OH on grains, unlike CH₃OH, H₂CO can also form in the cold gas-phase through neutral-neutral reactions between CH₂ and O₂ or CH₃and O (Ramal-Olmedo et al. 2021) or ion-neutral reactions between ${\mathrm{C}\mathrm{H}}_3^{+}$ and OH (Woodall et al. 2007). This reinforces the role of H₂CO as an excellent tracer of gas kinetic temperature (Mangum & Wootten 1993). In addition with increasing UV-radiation field the PDR model produces higher abundances of H₂CO consistent with observation result toward massive star forming regions (Immer et al. 2014; Tang et al. 2018b).

Our Meudon PDR models (Fig. A.5) reproduce H₂CO at gas densities consistent with our non-LTE PyRADEX results but at T_kin values of a few times 100 K. While a subset of p–H₂CO velocity components shows high temperature (Fig. 3) the PDR models fail to reproduce significant amount of cooler H₂CO. This could be because we carried out single-pointing observations which represent averaged spectra toward each clump, the multicomponent profiles likely arise from multiple molecular clumps within the beam, each with distinct velocities and physical conditions. Additional sources of uncertainty include the source specific n(H₂)–T_kin profile assumed. Moreover, the broader, higher-temperature components probably trace gas affected by outflows, inflows, or shocks. In the following section, we examine the origin and nature of these “H₂CO Sub” components in detail.

Fig. 9 Schematic depicting a fully evolved H II region in which accretion onto the central star has ceased but the region remains embedded in its surrounding environment. Additionally, it outlines the various input parameters used in the Meudon PDR code and output molecular abundance as presented in Fig. A.5.

5.3 Nature of secondary component (H₂CO Sub)

The p–H₂CO observations reveal multiple velocity components in a large fraction of the sources in our sample. We classify these components as “H₂CO Main” and “H₂CO Sub” based on their centroid velocities, as described in Sect. 3.1. In this section, we investigate the nature and origin of the “H₂CO Sub” components. As discussed in Sect. 5.2, H₂CO forms in regions with layered density and temperature structures, which can produce multiple velocity components in the observed line profiles. The Meudon PDR models reproduce strong H₂CO emission from high-UV PDR (Fig. A.5). Because these layers lie near the surface of H II regions, they are subject to intense UV radiation, gas inflow/outflow, and turbulent motions, which broaden the observed lines.

Figure 10a shows the distribution of the p–H₂CO (3_0,3–2_0,2) FWHM, indicating that the H₂CO Sub components exhibit systematically broader lines than the Main components, consistent with turbulence-driven broadening. Using the T_kin derived in Sect. 4.2, we calculated the thermal and nonthermal velocity dispersions for each component as $${\sigma }_{\mathrm{T}}=\sqrt{\frac{k{T}_{\mathrm{k}\mathrm{i}\mathrm{n}}}{m_{{\mathrm{H}}_2\mathrm{C}\mathrm{O}}}},{\sigma }_{\mathrm{N}\mathrm{T}}=\sqrt{{\left(\frac{\mathrm{\Delta }v}{8\ln 2}\right)}^2-{\sigma }_{\mathrm{T}}^2},$$(4)

where k is the Boltzmann constant, Δv is the p–H₂CO (3_0,3–2_0,2) FWHM, and m_H2CO is the mass of a formaldehyde molecule. The thermal line width is significantly lower than the nonthermal line width for H₂CO with a median of 0.13 and 0.8 km s⁻¹, respectively, indicating that nonthermal processes dominate the line broadening. This result is consistent with Tang et al. (2018b). Figure 10b further shows that the H₂CO Sub components have higher nonthermal line widths (median ~0.92 km s⁻¹) than the Main components (median ~0.70 km s⁻¹). This enhancement is likely due to the fact the H₂CO Sub components arise from turbulence and shock motions near the surfaces of H II regions.

The median sound speed ($a_{\mathrm{s}}=\sqrt{\frac{k{T}_{\mathrm{k}\mathrm{i}\mathrm{n}}}{\mu m_{\mathrm{H}}}}$, where μ = 2.7 is the mean molecular weight for molecular clouds and m_H is the mass of the hydrogen atom) is ~0.45 km s⁻¹ at a temperature of 66 K. This gives the median Mach number (given as M = σ_NT/a_s) of 2.45 and 1.56 for the Sub and Main components, respectively, confirming that the H₂CO Sub components are dominated by supersonic nonthermal motions and trace strongly turbulent gas within our sample.

We examined the relationship between the nonthermal line width of p–H₂CO and the clump bolometric luminosity for the H₂CO transitions. Figure 10c presents this correlation, where the nonthermal line width serves as a proxy for turbulent motion. We fit a power-law relation for three cases: all components combined, H₂CO Main components, and H₂CO Sub components, with the fitting results summarized in Table 2. The nonthermal line width shows a positive correlation with L_bol, indicating that turbulence in molecular clouds increases with luminosity. This trend suggests that more luminous sources are associated with more turbulent molecular environments (Wang et al. 2009).

The relationship between σ_NT of p–H₂CO (3_0,3–2_0,2) and L_bol follows a power law of the form ${\sigma }_{\mathrm{N}\mathrm{T}}\propto L_{\text{bol}}^{0.1-0.16}$ (see Table 2). The slope is steeper for the H₂CO Sub component (0.16) than for the H₂CO Main component (0.10), supporting the interpretation that the H₂CO Sub emission originates from regions more influenced by massive stars, where turbulence is enhanced—likely near the ionization fronts of H II regions.

Fig. 10 Distributions of (a) FWHM and (b) σNT of the p–H2CO (30,3–20,2) transition. The blue and red histograms represent the H2CO Main and H2CO Sub components, respectively, with the corresponding dashed lines marking their median values. Panel c: relation between σNT and Lbol. The black, blue, and red dashed lines indicate the power-law fits for all components combined, the Main components, and the Sub components, respectively.

5.4 Molecular pressure and confinement of H II region

H II region expand under the ionized gas pressure but are impeded by the ambient molecular gas pressure. Figure 11 (left) compares the thermal pressures of the ionized and molecular gas components. For nearly all sources in our sample, the molecular gas pressure exceeds that of the ionized gas. This implies that the molecular gas pressure can slow down the expansion of embedded H II regions. Moreover showing that high molecular gas pressure persists regardless of the velocity offset between RRLs and p–H₂CO transitions, further supporting the argument that molecular pressure can confine the expansion of embedded H II regions.

However, several caveats exist. Firstly, these H II regions must still be embedded within their natal molecular clumps; our observations though limited by spectral resolution, try to fulfill this condition through our selection of H II regions with sizes smaller than their parental ATLASGAL clump (Sect. 2.1). Secondly, the molecular pressures derived from H₂CO should trace gas closer to the ionization front, this is supported by both observations and the Meudon PDR models (Sect. 5.2). A confirmation of the same awaits higher resolution H₂CO and ionized gas mapping. Thirdly, the ionized gas parameters (n_e, T_e) adopted from Paper I represent spatial averages that may underestimate the true pressure at the ionization boundary. Figure 11 (right) compares the ionized gas pressures with those derived from dust emission at 870 μm (Urquhart et al. 2018, 2022), representing the average clump conditions. In this case, the ionized gas pressure often exceeds that of the molecular component. In contrast to the molecular gas, the pressures inferred from dust emission may not be representative of conditions in the vicinity of the expanding H II regions. Resolving these apparent inconsistencies will require high-resolution observations of the ionized and molecular gas interfaces.

Recently, Faerber et al. (2025) analyzed the expansion signatures of 35 H II regions mapped in [CII] 158 μm emission. Among these, 12 H II regions were identified as candidates for expansion, and eight of them were consistent with a spherical expansion model. They also found that thermal pressure and stellar winds are the primary drivers of H II region expansion in their sample. Previous studies, including that of Anderson et al. (2011) and Paper I, reported that roughly 50% of cataloged Galactic H II regions exhibit bubble-like morphologies at 8 μm, based on GLIMPSE images (Benjamin et al. 2003), which trace PDRs. In our sample, 23 sources exhibit bubble morphology including irregular bubbles and of the remaining 19 sources show compact morphology. To explore these processes in more detail, maps of PDRs and dense molecular gas tracers in the vicinity of expanding H II regions are needed, along with radio continuum and RRL observations at optically thin frequencies. Such observations will allow a detailed analysis of both ionized and molecular gas, as well as high-resolution studies of the gas dynamics. This approach will provide stronger constraints on the formation, evolution, and feedback mechanisms of H II regions.

Fig. 11 Comparison between the thermal pressure of ionized and molecular gas. The color bar indicates the velocity difference between the RRL and p–H2CO transitions. Left panel: n(H2) and Tkin derived from p–H2CO in this work. Right panel: n(H2) and Tdust estimated from the 870 μm emission (Urquhart et al. 2018, 2022), representing the average clump pressure. The dashed black lines mark the slope of equal pressure. The plots are set to the same scale for visual comparison.

6 Conclusion and summary

We have investigated the physical properties traced by molecular and ionized gas in a statistically significant sample of H II regions associated with ATLASGAL dust clumps. We used observations carried out with the APEX 12 m submillimeter telescope of multiple transitions of p–H₂CO associated with its J=3–2 and 4–3 transitions toward these sources. In our earlier study (Paper I), we examined the ionized gas properties of H II regions using centimeter–RRLs from the GLOSTAR survey. In this work, we focused on the molecular gas properties of the embedded H II regions and compared them with the ionized gas characteristics. Our main findings are summarized below:

• We constrained the physical properties of massive clumps associated with H II regions by deriving the gas density and kinetic temperature for 52 clumps. We carried out non-LTE radiative transfer modeling of the six p–H₂CO transitions studied, which serve as effective tracers of density and temperature in dense gas, using PyRADEX+MCMC. The resulting T_kin ranges from 33.7 K to 265 K, with a median value of ~66 K and n(H₂) > 10⁵ cm⁻³ across our sample, and it aligns with findings from previous studies;

• We analyzed the relationships among the physical properties of dust, molecular gas, and ionized gas by calculating the Spearman correlation coefficients and p-values for various parameters. Properties that directly depend on the central high-mass star, such as L_bol and N_Lyc, exhibit stronger and more frequent correlations with other parameters. This suggests that the central high-mass star largely governs and influences the physical conditions of the surrounding dust, ionized gas, and molecular gas;

• By comparing H₂CO with other molecular gas tracers, we found that H₂CO traces denser gas located near the interface between the H II region and the surrounding molecular cloud. The Meudon PDR models reproduced H₂CO abundances peaking near the ionized region’s surface and at gas densities similar to those obtained from PyRADEX but at high values of T_kin. Understanding the origins of the cooler H₂CO components requires further investigation with more detailed modelling;

• The nature of the “H₂CO Sub” components reveals that these components are dominated by supersonic nonthermal motions and trace turbulent gas;

• We investigated whether the surrounding molecular gas confines the expansion of H II regions through excess pressure. Our analysis shows that for most of our sources, the pressure of the neutral gas exceeds that of the ionized gas, regardless of the velocity difference between the ionized and molecular components. This finding suggests that the surrounding neutral molecular gas can slow down or hinder the expansion of H II regions due to its higher pressure. However, there are several caveats, and a confirmation of this fining relies on future high resolution mapping of H II regions and their surroundings in RRL and molecular lines. Such observations will soon be achievable with the unprecedented sensitivity and spatial resolution offered by next-generation radio facilities such as the ngVLA and SKA–Mid.

The expansion and feedback of H II regions significantly shape the ISM, regulate star formation rates, and influence galactic evolution. Future studies should increase the sample size of embedded and expanding H II regions to better explore the relationship between ionized and molecular gas. This will improve our understanding of how H II regions evolve over time and impact their surrounding molecular environment. Placing tighter constraints on the connection between ionized and molecular gas across various evolutionary stages ranging from compact to classical diffuse H II regions will be essential for understanding their development and potential role in triggering star formation.

Data availability

Tables C.1, C.2, and C.3 are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/706/A280

Acknowledgements

During the preparation of this article, we suffered the tragic loss of Prof. Dr Karl M. Menten, S.K.’s PhD supervisor and the principal investigator of the GLOSTAR survey. We dedicate this work to his enduring legacy. We are deeply indebted to the countless discussions with him that have inspired generations of radio astronomers. We will always miss his warmth, wisdom, boundless curiosity, generosity of spirit, kindness, and delightful sense of humor. We sincerely thank the referee for their positive, timely, and constructive feedback. We thank MPIfR and APEX staff for operating the observatory and collecting the data (https://www.apex-telescope.org/ns/). This research made use of information from the ATLASGAL database at (http://atlasgal.mpifr-bonn.mpg.de/cgi-bin/ATLASGAL_DATABASE.cgi) supported by the MPIfR, Bonn. A.M.J. acknowledges support by the CRC1601 (SFB 1601 subproject B2) funded by the DFG (German Research Foundation) – 500700252. S.K. and A.M.J. thanks the Max Planck Society for their support. Y.G. was supported by the Ministry of Science and Technology of China under the National Key R&D Program under grant No. 2023YFA1608200, the National Natural Science Foundation of China (NSFC) under grant No. 12427901, and the Strategic Priority Research Program of the Chinese Academy of Sciences under grant No. XDB0800301. S.N. gratefully acknowledges the Collaborative Research Center 1601 (SFB 1601 subproject B1) funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – 500700252. This research has made use of the SIMBAD database and the VizieR catalog, operated at CDS, Strasbourg, France. This research made use of Astropy (http://www.astropy.org), a community developed core Python package for Astronomy (Astropy Collaboration 2022). This document was prepared using the collaborative tool Overleaf available at: https://www.overleaf.com/.

References

Appendix A Additional figures

Fig. A.1 Three-color composite images of the targets of this study. The images combine mid–infrared data from Herschel 70 μm (red), MIPSGAL 24 μm (green), and GLIMPSE 8.0 μm (blue), overlaid with the GLOSTAR 5 GHz radio continuum contours in white. The filled white circle in corner represents the 25" beam of the radio continuum data from Paper I. Continued.

Fig. A.2 Comparison of physical size of H II regions and associated ATLASGAL clump size taken from Csengeri et al. (2014). The dashed black line shows the line of equality. The red circle represents only unresolved H II region, G030.720–0.082, in our sample.

Fig. A.3 Distribution of (a) angular and (b) physical deconvolved FHWM size of H II region. The vertical dashed line shows the median value of the distribution.

Fig. A.4 Example corner plot for AGAL045.121+00.131. To constrain physical properties of the molecular clouds we used the non-LTE RADEX modelling with MCMC algorithm, here we show the corresponding corner plot for each Gaussian component displaying the 1D histogram of the posterior distribution for each parameter (kinetic temperature, Tkin, H2 density, n(H2), and p–H2CO column density, N(p–H2CO)). The histograms represent the PDFs of the Monte Carlo fit parameters, while the blue curves represent the scatter plot of the joint distribution between the two parameters. The orange solid lines indicate the median of the PDFs, representing the most likely fit parameters, and the dotted blue lines correspond to the 1σ confidence intervals. The subtitle of each corner plot shows most likely value of the parameters with upper and lower 1σ error.

Fig. A.5 Meudon PDR model results for a constant thermal pressure of Pth ~ 4.6 × 107 K cm−3. Left panel: Meudon PDR model run with input spectral type of O5V, distance to PDR of d0 = 0.26 pc, and far-UV flux, G0 ~ 5.3 × 104 in units of Habing field. Right panel: Meudon PDR model run with input spectral type of O 5 V, distance to PDR of d0 = 0.12 pc, and far-UV flux, G0 ~ 2.5 × 105 in units of Habing field. The top plots shows the evolution of hydrogen species – H2, H, and H+, as well as the gas temperature as a function of visual extinction AV. The middle and bottom plots present the spatial density profiles of key PDR and dense gas tracer molecules (see main text), where each density has been normalized to its maximum value within the cloud.

Appendix B Effect of distance on source selection

In Sect. 2.1, we outlined the criteria used to select the sources, and Fig. 1 illustrates their distance distribution. To compare the physical conditions of the ionized and dense molecular gas, we adopted ionized gas parameters from Paper I. Because these parameters vary with distance, we examined the distance dependence within our sample. Figure B.1 shows the trends of N_Lyc, n_e, and T_e as a function of distance. While n_e exhibits only a weak dependence and T_e remains largely independent of distance, N_Lyc increases systematically, reflecting its strong distance dependence. This spread in N_Lyc indicates that our sample encompasses H II regions powered by a wide range of spectral type of central OB star. Overall, the distance variation introduces no significant bias in the derived ionized gas properties but rather ensures a diverse and representative sample.

Fig. B.1 Relation between the distances and (a) NLyc, (b) ne, and (c) Te for our sample.

Appendix C Tables

Appendix D Optical depth

Since much of the analysis presented above assumes optically thin conditions, we here examine the validity of this assumption. To assess the influence of optical depth on the derived quantities, we display in Fig. D.1 the PyRADEX+MCMC reproduced optical depths, τ, for the observed p–H₂CO (3_0,3–2_0,2) transitions (which has the lowest value of E_u and the brightest observed transition) across the range of gas densities and p–H₂CO column densities probed but for fixed values of T_kin, at median T_kin of 66 K (Fig. D.1 left) and at a higher T_kin of 200 K (Fig. D.1 right).

The background color scale represents the counts visualized over hexagonally binned distributions across the derived n(H₂)–N(p–H₂CO) space for our sample. The solid black curves mark contours of the optical depth equal to 1 and 5 at T_kin of 66 and 200 K. The blue dashed curves represent the beam filling factor corrected line intensities of the p–H₂CO (3_0,3–2_0,2), corresponding to the brightest (T_MB =8.26 K), median (T_MB =0.94 K), and the lowest (T_MB =0.18 K) detection in our sample. We find that all the observed p–H₂CO transitions are optically thin (τ < 1) at column densities N(p–H₂CO) < 10¹⁴ cm⁻² and spatial densities n(H₂) = 10^3–8 cm⁻³. We also notice that the transition becomes very optically thick (τ > 5), at column densities higher than 5 × 10¹⁵ cm⁻². The spatial density and column density derived using the observed p–H₂CO line intensities, shown as colored hexagonal bins in Fig. D.1, lie predominantly in or near the optically thin region (τ < 1). Among the observed p–H₂CO transitions, the low-lying p–H₂CO (3_0,3–2_0,2) line is the brightest, while the other transitions, which have higher, upper energy levels (E_u > 34 (K); see Table 1), are expected to have much lower optical depths (τ ≪ 1). This indicates that optical depth effects have only a minimal impact on the derived spatial densities and kinetic temperatures from the p–H₂CO transitions. Furthermore, as was found by Tang et al. (2018b), varying the kinetic temperature only weakly affects the optical depth of the p–H₂CO lines as evident from Fig. D.1.

Fig. D.1 Optical depth of p–H2CO (30,3–20,2) as a function of n(H2)-N(p–H2CO) for a fix value of Tkin. The color hexagonally binned distributions display the PyRADEX+MCMC constrained n(H2) and N(p–H2CO) of our sources. The black solid curves show the optical depth of τ = 1 and 5 for p–H2CO (30,3–20,2) transition on top of n(H2) and N(p–H2CO) space, assuming fixed kinetic temperature of 66 K (the median temperature of our sample; left panel) and 200 K (higher temperature; right panel). The blue dashed curves indicate the contour corresponding to the brightest, average, and faintest observed intensities for p–H2CO (30,3–20,2) transition (see the text).

All Tables

All Figures

Fig. 1 Galactic distribution of the source sample. The locations of the targets are denoted by red stars on an artistic rendering of the Milky Way. The distances to the observed H II regions and their associated clumps are taken from Paper I (and references therein). The spiral arms are marked following the description presented in Reid et al. (2019), and the individual spiral arms are colored as follows: the 3 kpc in magenta, the Norma-Outer Arm in red, the Scutum-Centaurus Arm in blue, the Sagittarius-Carina Arm in green, the Local Arm in cyan, and the Perseus Arm is in yellow.

In the text

Fig. 2 Spectra of the different p–H2CO and RRL transitions studied in this work toward AGAL045.121+00.131. Left panels: spectra of all six of the observed p–H2CO transitions (black) along with the multicomponent Gaussian fitting (green and magenta) and residual of the fitting (dashed blue). The total fit is displayed in red and the gray shaded region demarcates the 3σ level of the rms noise. The vertical dashed gray lines mark the centroid velocities of each fit Gaussian component. Right panel: spectra of the stacked GLOSTAR ⟨Hnα⟩ RRL in black alongside the Gaussian fit in red. Also displayed is the 5.3 GHz radio continuum emission map obtained with the GLOSTAR VLA D configuration toward this source and overlaid with white contours is the ATLASGAL 870 μm. The white filled circle shows the 25″ GLOSTAR beam.

In the text

	Fig. 3 Distribution of the sources in the n(H2)–Tkin plane. The blue and red colors respectively represent sources with single and multiple component Gaussian profiles. Filled circles denote components classified as “H2CO Main,” while stars indicate those identified as “H2CO Sub.”
In the text

Fig. 4 Distributions of the pyRADEX+MCMC derived (a) n(H2), (b) Tkin, (c) N(p–H2CO), and (d) p–H2CO fractional abundance, X(p–H2CO) toward the different velocity components identified in our sample, as traced by the p–H2CO. The blue and red histograms represent the distributions for the “H2CO Main” and “H2CO Sub” components, respectively, while the corresponding dashed lines indicate the median values of each distribution.

In the text

	Fig. 5 Comparison between n(H2) probed by 870 μm emission from dust and p–H2CO as derived from our pyRADEX+MCMC models. The dashed black line displays a slope of unit value.
In the text

Fig. 6 Comparison between the kinetic temperatures derived from pyRADEX+MCMC modeling and using the p–H2CO integrated line intensities ratios R1 (top panel) and R2 (middle panel). Bottom panel: comparison of the kinetic temperatures derived from R1 and R2. The LTE temperature uncertainties are from uncertainties in the integrated intensities of the observed p–H2CO lines. The dashed black lines indicate a slope of unity.

In the text

	Fig. 7 Correlation matrix showing the relationships between various physical properties of molecular gas, ionized gas, and dust clumps. The color bar and annotations indicate the Spearman’s correlation coefficients. We flagged and excluded correlations with p-values > 0.0013 (3σ), thereby removing statistically insignificant relationships.
In the text

Fig. 8 Cumulative distribution functions for temperatures (top panel) and line width (bottom panel) for different H II regions. Different colors represent various tracers and components, as indicated in the legend. The NH3 properties are from Wienen et al. (2012). The CH3CCH, CH3OH (hot and cold), and CH3CN (hot and cold) properties are from Giannetti et al. (2017), whereas dust properties are from Urquhart et al. (2018).

In the text

	Fig. 9 Schematic depicting a fully evolved H II region in which accretion onto the central star has ceased but the region remains embedded in its surrounding environment. Additionally, it outlines the various input parameters used in the Meudon PDR code and output molecular abundance as presented in Fig. A.5.
In the text

Fig. 10 Distributions of (a) FWHM and (b) σNT of the p–H2CO (30,3–20,2) transition. The blue and red histograms represent the H2CO Main and H2CO Sub components, respectively, with the corresponding dashed lines marking their median values. Panel c: relation between σNT and Lbol. The black, blue, and red dashed lines indicate the power-law fits for all components combined, the Main components, and the Sub components, respectively.

In the text

Fig. 11 Comparison between the thermal pressure of ionized and molecular gas. The color bar indicates the velocity difference between the RRL and p–H2CO transitions. Left panel: n(H2) and Tkin derived from p–H2CO in this work. Right panel: n(H2) and Tdust estimated from the 870 μm emission (Urquhart et al. 2018, 2022), representing the average clump pressure. The dashed black lines mark the slope of equal pressure. The plots are set to the same scale for visual comparison.

In the text

	Fig. A.1 Three-color composite images of the targets of this study. The images combine mid–infrared data from Herschel 70 μm (red), MIPSGAL 24 μm (green), and GLIMPSE 8.0 μm (blue), overlaid with the GLOSTAR 5 GHz radio continuum contours in white. The filled white circle in corner represents the 25" beam of the radio continuum data from Paper I. Continued.
In the text

	Fig. A.2 Comparison of physical size of H II regions and associated ATLASGAL clump size taken from Csengeri et al. (2014). The dashed black line shows the line of equality. The red circle represents only unresolved H II region, G030.720–0.082, in our sample.
In the text

	Fig. A.3 Distribution of (a) angular and (b) physical deconvolved FHWM size of H II region. The vertical dashed line shows the median value of the distribution.
In the text

Fig. A.4 Example corner plot for AGAL045.121+00.131. To constrain physical properties of the molecular clouds we used the non-LTE RADEX modelling with MCMC algorithm, here we show the corresponding corner plot for each Gaussian component displaying the 1D histogram of the posterior distribution for each parameter (kinetic temperature, Tkin, H2 density, n(H2), and p–H2CO column density, N(p–H2CO)). The histograms represent the PDFs of the Monte Carlo fit parameters, while the blue curves represent the scatter plot of the joint distribution between the two parameters. The orange solid lines indicate the median of the PDFs, representing the most likely fit parameters, and the dotted blue lines correspond to the 1σ confidence intervals. The subtitle of each corner plot shows most likely value of the parameters with upper and lower 1σ error.

In the text

Fig. A.5 Meudon PDR model results for a constant thermal pressure of Pth ~ 4.6 × 107 K cm−3. Left panel: Meudon PDR model run with input spectral type of O5V, distance to PDR of d0 = 0.26 pc, and far-UV flux, G0 ~ 5.3 × 104 in units of Habing field. Right panel: Meudon PDR model run with input spectral type of O 5 V, distance to PDR of d0 = 0.12 pc, and far-UV flux, G0 ~ 2.5 × 105 in units of Habing field. The top plots shows the evolution of hydrogen species – H2, H, and H+, as well as the gas temperature as a function of visual extinction AV. The middle and bottom plots present the spatial density profiles of key PDR and dense gas tracer molecules (see main text), where each density has been normalized to its maximum value within the cloud.

In the text

	Fig. B.1 Relation between the distances and (a) NLyc, (b) ne, and (c) Te for our sample.
In the text

Fig. D.1 Optical depth of p–H2CO (30,3–20,2) as a function of n(H2)-N(p–H2CO) for a fix value of Tkin. The color hexagonally binned distributions display the PyRADEX+MCMC constrained n(H2) and N(p–H2CO) of our sources. The black solid curves show the optical depth of τ = 1 and 5 for p–H2CO (30,3–20,2) transition on top of n(H2) and N(p–H2CO) space, assuming fixed kinetic temperature of 66 K (the median temperature of our sample; left panel) and 200 K (higher temperature; right panel). The blue dashed curves indicate the contour corresponding to the brightest, average, and faintest observed intensities for p–H2CO (30,3–20,2) transition (see the text).

In the text