© ESO 2018

1 Introduction

Massive stars affect their environment in various ways, thereby shaping the evolution of galaxies and leading to chemical enrichment of the interstellar medium (Zinnecker & Yorke 2007). Due to fast timescales of high-mass star formation, their whole evolution takes place whilst they are still deeply embedded in their natal clump. Indeed, the high-mass star or cluster is hidden at optical and even near-infrared wavelengths. For these reasons, it is more challenging to make observations of high-mass star formation than those of low-mass star formation, and also they are rare and located at much larger distances (Schuller et al. 2009; Urquhart et al. 2018). For these reasons, the feedback between young OB stars and their molecular clouds is still not well understood.

Young massive stars are hot and thus create small regions of hot ionized gas around them that expands into their local environment and interacts with the surrounding molecular gas; these are known as compact H II regions. The H II region phase is the last stage of the massive star formation process. Studying this phase is not only important due to the exciting stars’ influence on their environment, but it is also key to understand how high-mass stars obtain masses above 10 M_⊙ (Churchwell et al. 2010). Indeed, the youngest and most compact H II regions (hyper- and ultra-compact) are often associated with various star formation signposts such as molecular outflows, infall motions, and methanol and water masers, which themselves are pumped by radiation and shocks (Urquhart et al. 2013b, 2015). These signposts suggest that the mass assembly process is still active.

The gas in H II regions surrounding newly born OB stars has velocity dispersions of 25−35 km s⁻¹ determined by a mix of thermal (~20 km s⁻¹) and turbulent motions (Wilson et al. 2009). The kinematics in the ionized gas can be investigated using radio recombination lines (RRLs), particularly toward sources that are deeply embedded in dense molecular clumps. In these dense clumps, observations of Lyman α, Hα, and other ultraviolet, optical, or near-infrared recombination lines are strongly attenuated due to high levels of visual extinction (up to hundreds of mag.) while RRL emission can still escape out of the clumps. Therefore, RRLs are an excellent tool for analyzing the distribution and kinematics of the ionized gas associated with H II regions.

Interestingly, strong maser emission has been observed from submillimeter RRLs toward only a few sources (Thum et al. 1995; Martín-Pintado 2002; Sánchez Contreras et al. 2017). Therefore, it is interesting to investigate how widespread RRL maser emission is towards H II regions.

Submillimeter-RRL observations provide kinematics of H II regions and probe the Lyman continuum photon production rate, which can be used as a measure of star formation rates (SFRs) in external galaxies (Scoville & Murchikova 2013; Bendo et al. 2017). These rates have been compared with those measured by other SFR tracers such as ultraviolet continuum emission, optical/near-infrared recombination lines, mid-/far-infrared continuum emission, and radio continuum emission and have been found to give consistent results (e.g., NGC5253; Bendo et al. 2017).

In our previous study of hydrogen RRLs at millimeter wavelengths (Kim et al. 2017), we identified 178 mm-RRL sources toward 976 compact dust clumps selected from the Atacama Pathfinder EXperiment (APEX) Telescope Large Area Survey of the Galaxy (ATLASGAL) Compact Source Catalogues and GaussClump Source Catalogue (CSCs; Contreras et al. 2013; Urquhart et al. 2014a, and GCSC; Csengeri et al. 2014). These mm-RRLs were found to be associated with embedded H II regions identified from cm-wavelength radio continuum surveys (e.g., Urquhart et al. 2007, 2013b) and were used to derive the properties of the ionized gas, which in turn were compared with the properties of their natal molecular clumps. The properties of the mm-RRLs were found to be consistent with properties derived from the 6 cm radio continuum.

To further characterize the properties of these H II regions, we have re-observed many of these sources targeting RRLs at submillimeter wavelengths, which not only provide brighter flux than mm-RRLs but also allow us to investigate the RRL excitation and whether maser emission is influencing the intensity of RRLs. We have, therefore, carried out submillimeter RRL (submm-RRL) observations using the APEX 12 m telescope to investigate further the H II regions detected in the course of our mm-RRL survey.

The observations and data reduction are explained in Sect. 2. The general properties of the detected submm-RRLs, including the detection rates, and a comparison between submm-RRL and mm-RRL are presented in Sect. 3. The associations with H II regions and molecular clumps are described in Sect. 4. The derivation of the ionizing Lyman photon production rate from the submm-RRL data is described in Sect. 5. Several sources with peculiar line profiles are discussed in Sect. 6. Lastly, we present a summary of our main results in Sect. 7.

2 Observation and data reduction

2.1 Source selection and observational setup

The sources for this targeted survey with the APEX 12 m diameter submillimeter telescope (Güsten et al. 2006) have been selected based only on their peak mm-RRL intensity as measured from our surveys with the Institut de Radioastronomie Millimétrique (IRAM) 30 m and the Mopra 22 m telescopes (Kim et al. 2017). For the selection, the intensity threshold was set to 0.2 K of the peak mm-RRL intensity. This threshold was chosen by comparing the intensity of stacked mm-RRL profiles¹ to H26α RRL data that were observed by APEX as part of a different molecular line survey toward the ATLASGAL “Top 100” sample (see Giannetti et al. 2014; König et al. 2017 for details). In total, 104 ATLASGAL sources were selected from the mm-RRL catalog (Kim et al. 2017), and the submm-RRL observations with the APEX 12 m telescope were carried out primarily in 2015 (Project IDs: M0025-95 and M0018-96) with some additional observations being made in 2016.

The APEX observations used the dual frequency First Light APEX Submillimeter Heterodyne (FLASH) receiver (Klein et al. 2014) for the H25α, H27α, H28α, and H35β transitions and the APEX-1 (HET230; Vassilev et al. 2008) receiver for the H29α and H30α transitions. The observations mainly covered the H25α, H28α, and H35β transitions. However, the H27α, H29α, and H30α transitions were also observed to obtain better line profiles toward some sources for which the other RRL profiles were found to be contaminated by emission from molecular lines. The spectra of the observed submm-RRL transitions, therefore, vary slightly from source to source. In Table 1 we provide details of the transitions observed toward each source and indicate where they are detected. In Table 2 we give the observed frequency, the oscillator strength of each transition, and the full width at half maximum (FWHM) beam widths of the APEX beam at the frequency of the transition, the conversion factor between K and Jy, and the telescope’s main beam efficiency, η_MB ². The beam sizes given in Table 2 are calculated using ${\Theta }_{\text{FWHM}}=7.8^{″}\times \left(800/\nu \right)$³, where ν is a rest frequency of the transition in GHz.

2.2 Datareduction and identification of blended molecular lines

In the submillimeter wavelength regime, the spectra of hot molecular cores, that is, warm, dense regions around embedded massive young stellar objects, are dominated by many lines from complex molecules. As a result, molecular lines can overlap with the submm-RRLs, making the detection and line fitting more complicated, particularly when trying to distinguish between potential RRL maser emission, which can have spectral features as narrow as typical molecular lines (Thum et al. 1995; Martín-Pintado 2002). Therefore, it is necessary to identify any molecular line contribution to observed narrow features on a broad profile.

To identify narrow line features, we used WEEDS within the Continuum and Line Analysis Single-dish Software (CLASS) of the GILDAS package⁴; CLASS is also used for the RRL data reduction. Figure 1 shows an example of identified emission lines of the CH₃OCHO molecule and a submm-RRL (H25α) toward AGAL034.258+00.154. We fit the CH₃OCHO molecular emission assuming local thermal equilibrium (LTE) and that the lines are exposed to the same physical conditions (single value of excitation temperature, column density, linewidth, and peak velocity to all transitions of this molecule within observed frequency bands). The fits to the CH₃OCHO line data obtained from WEEDS and the H25α line are shown in red and green, respectively. The broad profile of the H25α line is consistent with that of other submm-RRL transitions observed towards this source with regard to peak velocity, intensity, and linewidth. In addition, we find that, while the H26α transition is not blended with molecular lines nearby in frequency, it is contaminated in several cases by attenuated strong emission from a CS (7−6) transition (ν = 342 883 MHz) originating from the other side-band of the receiver.

All narrow features could be identified as blended complex lines from molecules. The observed broad features originate from hydrogen submm-RRLs. We, therefore, only use the line parameters determined from fits to the broad velocity features for the analysis in this paper. A polynomial baseline of order one to three was fitted to the line-free channels in a 300 km s⁻¹ wide velocity range and subtracted from the spectrum of each transition. The average root mean square (rms) values for each transition are given in Table 2; these were determined from the line-free channels, which have a velocity resolution of 2.1 km s⁻¹. In some cases the submm-RRLs have been smoothed to a velocity resolution of 4.1 km s⁻¹ to increase the signal-to-noise ratio (S/N). We find that the majority of the submm-RRLs can be fitted with a single Gaussian component. However, for a handful of sources, the profiles have significant red or blue shifted wings and these were fitted with two components, a narrower component that fits the central portion of the line and a much broader component that fits the high-velocity wings. We will discuss those sources with two Gaussian RRL components in Sect. 6. Furthermore, we re-reduced all data of individual and stacked mm-RRLs presented in Kim et al. (2017) using the same velocity resolution as for the submm-RRLs to facilitate comparisons of the mm- and submm-RRLs (see e.g., Fig. 2).

3 Results

3.1 Detection rate

Since the intensity of RRLs increases with the frequency of the RRL (Gordon & Sorochenko 2002), high detection rates of submm-RRLs were expected from the observed mm-RRL sources. Table 3 gives the detection rate for each submm-RRL transition. The detection rates are lower than 100% because our observations are less sensitive than the mm-observations; the non-detections toward six sources might be due to poorer flux sensitivity (for H28α, rms = 0.50 Jy compared to the rms values of the Mopra and IRAM observations, which were 0.20 and 0.05 Jy, respectively). Furthermore, the detections of submm-RRLs towards 11 clumps are ambiguous due to blending from molecular lines that make it difficult to separate the molecular and submm-RRLs profiles; we considered these sources as non-detections, which are indicated in Table 1.

Furthermore, if the blending is not too severe as in the case of the H25α of AGAL034.258+00.154 (Fig. 1), the molecular emission and RRL can be reliably separated, and the detection of the submm-RRL is considered valid. The Hnα and H35β submm-RRL transitions (S∕N ≥ 3σ) have been detected towards a large fraction of the observed sources; these are detected towards 93 (89%) and 34 (40%) clumps, respectively. The detection rate of the H35β compared with the Hnα transitions is consistent with the results of the mm-RRL survey at a similar S/N level.

Fig. 1 Example of the identification of molecular lines blended with submm-RRL emission (black line) in AGAL034.258+00.154. The red profile shows a combined profile of the modeled profiles of CH3 OCHO lines using WEEDS and a Gaussian profile of the H25α RRL. The latter is also shown as the green dashed line. Other emission lines are unidentified molecular lines.

3.2 Properties of submm-RRLs

Figure 2 displays examples of detected submm-RRL spectra (upper panels) and the corresponding mm-RRL spectra (lower panels; Kim et al. 2017) toward AGAL037.874−00.399 and AGAL332.156−0.449. The profiles of the submm- and mm-RRLs show good agreement in peak velocity, linewidth, and profile shape; however, we note a significant difference in the intensity, with the submm-RRLs being approximately a factor of two brighter. Figure 3 shows a comparison of the peak flux of the H28α line to the peak fluxes measured for the other submm-RRL transitions (i.e., n = 25, 26, 27, 29, and 30) (upper panel), and also for the H35β transition (lower panel). We chose the H28α transition as the main transition for these comparisons because this transition is mainly covered by all the observations and is detected with a reasonable S/N level toward the observed sources.

Since we observed multiple submm-RRL transitions toward the majority of our sources, we can investigate whether there is a significant flux difference for different quantum numbers. Overall the peak fluxes determined from the different submm-RRLs agree with each other within the uncertainties (as shown in the upper panel of Fig. 3): the slope determined from a linear bivariate correlated errors and intrinsic scatter bisector (BCES bisector)⁵ fit to the data is 1.00 ± 0.04. In addition to the bestfit, the Pearson correlation coefficient (r) for the whole data set is 0.95 with p−value ≪ 0.0001, revealing that the RRL fluxes are very highly correlated. Although some sources show deviations from the primary trend, their S/N tends to be poor. The consistency of the submm-RRL fluxes indicates that there is no extreme density gradient toward the H II regions on angular scales of a few tens of arc-seconds, which would lead to differences in the measured fluxes (Dupree & Goldberg 1970). By comparing the linewidths of the various submm-RRLs similarly as we have just done for the fluxes, we also find a moderate correlation (r = 0.53 with p−value ≪ 0.0001). The correlation coefficient is considerable but the uncertainty in the linewidths tends to be significantly higher than for the peak fluxes.

The lower panel of Fig. 3 shows the comparison of peak fluxes of the H28α and H35β transitions. As we have seen, the detection rate of the H35β line is lower than that of the H28α line (Table 3) since the H35β lines are about six times fainter than the H28α lines. The slope from BCES bisector fitting is 0.95 and the correlation coefficient (r) of 0.99 with a small p−value ≪ 0.0001 reveals a significant correlation between both transitions. If maser excitation contributed to the observed fluxes, a low ratio of H35β/H28α fluxes is expected to be quite similar to the ratio measured in the submillimeter wavelength regime toward MWC349 (H35β/H28α ~ 0.078). The peak flux ratios of H28α and H35β in the lower panel of Fig. 3 are, however, consistent with the predicted LTE value. The average value of H35β/H28α toward each source, 0.31, is only slightly above the LTE prediction value of 0.28 (Dupree & Goldberg 1970; Thum et al. 1995) but close to the value of 0.32 measured from typical H II regions in the Galaxy (Thum et al. 1995).

In our previous analysis of the mm-RRLs (Kim et al. 2017), we identified two sources in which the observed fluxes could have been the result of maser emission; these were AGAL034.258+00.154 and AGAL043.166+00.011. They were identified by their low Hβ/Hα ratios, which is suggestive of non-LTE conditions due to an increase in the Hα intensity from stimulated maser amplification (as seen in MWC349A; Thum et al. 1995). The submm-RRL profiles for these two sources reveal no features that would indicate the presence of maser emission similar to those seen towards MWC349A at submillimeter wavelengths and their profiles are consistent with pure thermal emission. As discussed in the previous paragraph, the ratio of H35β/H28α of AGAL043.166+00.011 (cyan circle) shown in the lower panel of Fig. 3 does not show a significant offset from the LTE value. In the case of AGAL034.258+00.154, it was impossible to determine this ratio due to contamination of the H35β line. The Hnα line profiles of the source; however, already showed no maser emission. We can, therefore, rule out, based on the new submillimeter data, that AGAL034.258+00.154 and AGAL043.166+00.011 are associated with RRL maser emission.

Fig. 2 Examples of detected submm-RRLs (upper panels) of AGAL037.874−00.399 and AGAL332.156−00.449 and their corresponding mm-RRLs (lower panels) from Kim et al. (2017). The vertical purple dashed lines indicate the systemic velocity measured with H13CO+ (1−0) lines. The different colors present different submm/mm-RRL transitions detected in this paper and Kim et al. (2017). In flux unit, 1 Jy corresponds to ~0.025 K for the APEX 12 m, 0.169 K for the IRAM 30 m, and 0.045 K for the Mopra 22 m data.

Fig. 3 Upper panel: peak flux of the H28α transition versus peak flux of other transitions; H25α (black), H26α (red), H27α (cyan), H29α (orange), and H30α (bright green). All the plotted sources are detected in the H28α transition, and thus there are a few sources missing, which are only detected in other lines. Lower panel: peak flux of the H28α transition versus peak flux of H35β. In both plots, equal fluxes are indicated by the black dashed line, and blue lines indicate the best-fit determined by BCES bisector fits to all data points. The red line in the lower panel presents the predicted LTE value.

3.3 Properties of submm-RRLs in comparison to mm-RRLs

Table 5 compares the linewidth ratios of detected Hnα lines (i.e., n = 25, 26, 27, 28, 30, 39, 40, 41) to the H42α line (in this table we only consider detections with a S∕N > 8 σ). The low quantum numbers correspond to higher frequencies. Here we chose the H42α transition as the main transition to investigate intrinsic linewidths caused by non-thermal motions.

The median ratios are very close to unity for the whole range of transitions (H41α to H25α) and are all within one standard deviation. The similarity in linewidth implies that the different transitions are probing the same volume of gas under the same physical conditions. We note that there are some relative deviations from unity (see maximum and minimum ratios in Table 5). However, only a few particular sources seem to be affected, which might be related to micro-turbulence in the inner H II regions (Jaffe & Martín-Pintado 1999).

Since the linewidths and intensities are comparable in the whole range of submm-RRLs, it is feasible to stack the transitions of adjacent quantum numbers to improve their S/N towards each source. However, even by stacking the RRLs we still found some sources with poor S/N levels. We therefore plotted the uncertainties in the linewidth and peak intensity to determine reliable thresholds that can be used to select high S/N sources. From the distribution of the uncertainties, we chose the S/N threshold to be larger than 8σ for both the linewidths and peak fluxes of mm/submm-RRLs to allow for a reliable comparison with the stacked 3 mm data.

Figure 4 compares the linewidth and peak flux for the stacked mm-RRLs and submm-RRLs. In the upper panel of Fig. 4, the distribution of high-reliability linewidths of mm-RRLs and submm-RRLs reveals a good correlation (Pearson correlation coefficient r = 0.7 with p−value ≪ 0.0001). A linear least square BCES bisector fit to all of the data points is close to the line of equality. The median value and standard deviation of all the linewidths of mm-RRLs are 29.6 km s⁻¹ and 5.3 km s⁻¹, and those values for submm-RRLs are 28.7 km s⁻¹ and 6.2 km s⁻¹. Those standard deviations of mm-RRLs and submm-RRLs are just a factor of two-three more than their velocity resolution (~2 km s⁻¹). In addition, the sources with a good S/N (black circles) are even more tightly clustered around the line of equality. This linewidth comparison and the result in Table 5 are consistent with the hypothesis that both the mm-RRLs and submm-RRLs are probing the same ionized gas and that we can discard any significant influence of pressure broadening on the fitted parameters, which would have affected the mm-RRLs stronger than the submm-RRLs. In general, the pressure broadening is negligible at (sub)millimeter wavelengths except that toward some HII regions withhigh electron density the broadening effect still can affect broadening on lines even at millimeter wavelengths.

In the lower panel of Fig. 4, we compare the peak fluxes of mm- and submm-RRLs. The linear least-squares BCES bisector fit to all of the data reveals the presence of an offset from the line of equality, with the submm-RRL fluxes being significantly higher. The shift in the distribution from the line of equality is approximately a factor of two from mm (86−100 GHz) to submm (231−397 GHz) wavelengths. This shift agrees with the estimated difference in line intensities of Fig. 3 of Peters et al. (2012), assuming a H II region with an electron density of n_e = 5 × 10⁵ cm⁻³ and a temperature of T_e = 10⁴ K. This difference is due to the increase in the line to continuum ratio with increasing frequency expected for RRLs. Furthermore, the RRL intensity computation by Peters et al. (2012) only considers the line emissivity and does not include stimulated emissions. The good agreement we have found with the theoretical predictions is consistent with the detected mm-RRLs and submm-RRLs being emitted from optically thin H II regions. Nevertheless, there are some sources where the fluxes deviate significantly from the observed trend possibly indicating the presence of some H II regions with different physical conditions.

Fig. 4 Upper panel: linewidth comparison of stacked mm- and submm-RRLs. The black filled circles indicate sources with a peak flux with S/N ≥ 8 σ of both RRLs. Lower panel: peak flux comparison of stacked mm- and submm-RRL. The black symbols are the same as in the upper panel. The black dashed lines indicate the locus of equal linewidths and fluxes. The blue lines show the best fit from BCES bisector fitting.

4 Association with molecular clumps

Comparing the peak velocities of the ionized gas determined from the submm- and mm-RRLs and those of the molecular gas determined from the H¹³ CO⁺ (1−0) transition (Kim et al. 2017), we find peak velocity differences ($v_{\text{H}}^{13}{\text{CO}}^{+}-v_{\text{RRL}}$) of them to be in good agreement with a standard deviation of 4.8 km s⁻¹. This is consistent with the hypothesis that the RRLs are associated with young H II regions that are either embedded or still associated with their natal clouds. Thus, the properties of these H II regions may be related to the properties of their molecular clumps. We have, therefore, tried to identify correlations between the ionized and molecular gas using properties of radio emission from the Co-Ordinated Radio “N” Infrared Survey for High-mass star formation (CORNISH; Hoare et al. 2012; Purcell et al. 2013) and Red Midcourse Space Experiment (MSX; Lumsden et al. 2013) 5 GHz radio continuum surveys (Urquhart et al. 2007, 2009) and the properties of dust clumps determined from submm-dust emission of the ATLASGAL survey (Schuller et al. 2009; König et al. 2017; Urquhart et al. 2018).

In Table 6, we provide a summary of the properties of molecular clumps and H II regions toward a small sample of selected sources (the full table is available at the CDS). The sources selected for this table are those that are associated with non-Gaussian RRL profiles, which are discussed in Sect. 6. The peak velocities and linewidths for the first and second Gaussian components were obtained from fits to the stacked mm-RRLs. The distance, dust temperature (T_dust), bolometric luminosities (L_bol) of embedded objects, clump mass (M_clump), and L_bol∕M_clump ratio are taken from Urquhart et al. (2018) and the Lyman photon flux is determined from the 6 cm radio continuum sources correlated with the mm-RRLs (Kim et al. 2017).

Figure 5 shows the relation between the Lyman continuum photon flux emitted by the embedded massive star(s) and the bolometric luminosity generated by the whole cluster of protostellar objects embedded in the molecular clumps. The distribution of bolometric luminosities is shifted above the track of the zero-age main sequence (ZAMS) luminosities for the spectral type of a single star (Panagia 1973), probably because the bolometric luminosity measured is the contribution of a whole cluster of embedded objects (Urquhart et al. 2013b). Only some sources with Lyman photo fluxes equivalent to B-type stars lie below the ZAMS values (red dotted line). Such excess photon fluxes of the B-type stars agree well with the results of several previous studies (e.g., Sánchez-Monge et al. 2013; Urquhart et al. 2013b). The authors showed that the distributions of young B-type H II regions are consistent with the predicted blackbody with the same radius and temperature as a ZAMS star rather than the prediction of Lyman continuum emission (red dotted line).

Except for a few outliers, we find a significant correlation between the two parameters (Spearman correlation coefficient excluding the three outliers, ρ = 0.8 with p−value ≪ 0.0001). The correlation coefficient and the distribution (the BCES bisector fit with a blue line) show a clear correlation, but both bolometric luminosity and Lyman photon flux have a significant distance dependence. We have estimated a partial Spearman correlation (r_{AB, C}) test of bolometric luminosity and Lyman photon flux to eliminate their distance dependence using Eq. (4) of Urquhart et al. (2013a). We obtained a partial Spearman correlation coefficient of 0.5 and Student’s t−value of 3.8. These results of the partial Spearman correlation test also show there is a reasonable correlation between the two parameters, rejecting the null hypothesis and being independent of distance. The good correlation between the Lyman flux and thebolometric luminosity would suggest that the luminosity of the clumps is dominated by the most massive stars in the forming protocluster (cf. Urquhart et al. 2013b). The sources that are best fitted with two Gaussian components, indicated by the purple squares, tend to be the most luminous in the sample except for one (AGAL034.258+00.154). By looking at the dust temperature, we find they tend to increase with the bolometric luminosity and the Lyman photon flux increases, although there are some deviations from this trend. This is consistent with the feedback (radiation, outflows, and strong winds) from the embedded H II regions increasing the temperature in surrounding molecular material. The increase in the dust temperatureis likely to be linked to the evolutionary status of the central objects and as their luminosities increase so do the temperatures of the natal clumps.

In Fig. 5, three outliers (AGAL010.168−00.362, AGAL331. 546−00.067, and AGAL345.408−00.952) significantly deviate to the left from the linear fit. These all have higher bolometric luminosities compared to the Lyman photon flux than the other sources, and we also note that two have a significantly higher dust temperatures (T_dust ≥ 36.5 K; AGAL010.168−00.362 and AGAL345.408−00.952). The cause of the high dust temperatures toward the outliers is likely to be that sizes of the radio continuum sources are much smaller than the aperture used for the infrared and submm photometry and subsequent spectral energy distribution (SED) analysis (see König et al. 2017 for details of the photometry and fitting). Indeed, inspection of the Infrared Array Camera (IRAC) three-color composite images (3.6, 4.6, and 8 μm) of the Galactic Legacy Infrared Mid-Plane Survey Extraordinaire (GLIMPSE) reveals that all three are associated with extended mid-infraredemission (see e.g., Fig. 6). It is also clear from these images that the bolometric luminosity is arising from an extended star formation complex while the radio emission is associated with only a small part of the complex.

Figure 7 presents comparisons of three parameters, the luminosity of the submm-RRLs, bolometric luminosity, and Lyman photon flux toward the submm-RRL detected clumps. In addition to the strong correlation seen in Fig. 5, the luminosity of the submm-RRLs shows pronounced correlations with the bolometric luminosity and Lyman photon flux. The Pearson (r) correlation coefficient gives a high correlation (0.9) with a small p-value ≪ 0.0001. The three parameters increase altogether. It means that bright and hot stars emit stronger UV radiation and thus increase brightness of submm-RRLs and furthermore bolometric luminosity.

The good correlations between Lyman photon flux, submm-RRL, and bolometric luminosities and the high dust temperature (T_dust > 24 K) suggest thatthese submm-RRL sources are already quite evolved. We have plotted the clump mass as a function of the bolometric luminosity for all the ATLASGAL clumps (Urquhart et al. 2018) which contain non-detected and unobserved sources (black dots), and submm-RRL sources (green circles) in Fig. 8. It is evident from this plot that the submm-RRL sources in this survey are a part of the most luminous sources in the whole inner-Galactic plane mapped by ATLASGAL. Figure 9 shows a histogram of the bolometric luminosity to clump mass ratio (L_bol∕M_clump) for the full ATLASGAL sample with the distribution of the submm-RRL associated clumps over-plotted in green. The L_bol∕M_clump ratio is often used as a diagnostic for evolution (e.g., Eden et al. 2013) and shows that the submm-RRL sources are some of the most evolved young stellar objects in the Galaxy. The mm- and submm-RRLs analysis presented here and in Kim et al. (2017) have, therefore, parameterized a sample of some of the most luminous and evolved compact objects in the Galaxy.

Fig. 5 Bolometric luminosity (Lbol) as a function of Lyman continuum photon flux of 6 cm radio continuum emission integrated over the beam of the 3 mm observations. The colors indicate dust temperatures in the molecular clumps. The blue line shows the result of a fit to the datapoints using BCES bisector fitting. The purple squares superposed on the circles are particular sources with two Gaussian component RRLs (see Table 6 and Sect. 6). The red dotted line and red stars indicate bolometric luminosity and Lyman continuum photon flux corresponding to each spectral type (O4 − B0.5) for ZAMS stars from Panagia (1973).

Fig. 6 GLIMPSE IRAC three-color composite image (blue; 3.6 μm, green; 4.6 μm, and red; 8 μm) of one (AGAL345.408−00.952) of the three outliers in Fig. 5. The gray contours represent the 870 μm dust continuum emission from the ATLASGAL survey. The white circle indicates the size of Mopra beam (FWHM, 36″) and the red circle presents a RMS radio continuum source (Urquhart et al. 2007).

Fig. 7 Luminosity of the submm-RRLs versus Lbol of embedded central objects. The color bar indicates the Lyman photon flux of ionizing stars associated with the dust clumps and the submm-RRLs. The best linear fit is presented with the blue line.

Fig. 8 Bolometric luminosity (Lbol) as a functionof Mclump toward the full ATLASGAL sample (Urquhart et al. 2018) that indicates both non-detected and unobserved sources by this submm-RRL survey. The bright green circles indicate sources with a submm-RRL detection.

5 Photoionizing photon production rate, Q

If we assume the detected submm-RRL lines probe most of the photoionizing stars in clumps within the beam of APEX, it is possible to measure the photoionizing photon production rate, Q. In principle, Lyman photon flux and Q measure thesame quantity: the number of emitted photons per second. The photoionizing photon production rate, Q, can be obtained from the submm-RRL flux using Eq. (1) given by Scoville & Murchikova (2013) and Bendo et al. (2017), $$\begin{array}{ccc}\frac{Q(\text{H}n\alpha )}{{\text{s}}^{-1}}\hfill & =\hfill & 3.99\times {10}^{24}\left(\frac{{\alpha }_B}{{\text{cm}}^3{\text{s}}^{-1}}\right)\hfill \\ & & \times {\left(\frac{{ϵ}_{\nu }}{\text{erg}{\text{s}}^{-1}{\text{cm}}^{-3}}\right)}^{-1}\left(\frac{\nu }{\text{GHz}}\right){\left(\frac{D}{\text{kpc}}\right)}^2\left(\frac{{\displaystyle \int F_{\nu }}\text{d}v}{\text{Jy}\text{km}{\text{s}}^{-1}}\right),\hfill \end{array}$$(1)

where α_B is the total recombination coefficient of the sum over captures to all levels above the ground level in Case B (Osterbrock 1989). The effective emissivity, ϵ_ν, is at a frequency ν of observed RRL transition. These α_B and ϵ_ν parameters vary with the electron density and temperature.

We estimated the Q value for theH28α data since this transition is detected toward the majority of the observed sources. The α_B and ϵ_ν for our calculations are taken from the published values by Storey & Hummer (1995). Figure 10 shows the relationship (black circles) between the Lyman photon flux measured by 5 GHz radio continuum emission and the estimated Q(H28α) using the α_B and ϵ_ν based on the chosen T_e = 10⁴ K and n_e = 10⁴ cm⁻³. The cyan and red dashed and dotted lines present variations of the fits resulting from Q values computed using different electron densities and temperatures (T_e = 10⁴ K and n_e = 10⁴ cm⁻³). The variations are not significant in ranges of electron density (10³ −10⁵ cm⁻³) and temperature (7.5 × 10³−1.25 × 10⁴ K). There is a slight shift between the Lyman photon flux and the photoionizing rate Q(H28α) indicating that the single dish recombination line measurement picks up more photoionizing photons than the interferometric continuum measurements, which might resolve out emission. Overall, these two measurements show good agreement, and this confirms that the H28α lines are thermally excited. Furthermore, it confirms that the Q values measured by mm/submm RRLs are useful to measure SFR on extragalactic scales.

Fig. 9 Histogram of Lbol∕Mclump for the full ATLASGAL sources (gray) and the submm-RRL sources (bright green).

Fig. 10 Photoionzing photon rate (Q) measured by the submm-RRLs (this paper) as a function of Lyman photon flux estimated by 5 GHz radio continuum emission from CORNISH and RMS surveys (Urquhart et al. 2007; Purcell et al. 2013). Those two parameters have the same unit (s−1), the number of photons per time. The equality of both axes is indicated with a black dashed line. The color dotted and dashed lines present the best linear fits to the data points calculated for different electron temperature and density in Eq. (1). The black data points are estimated with Te = 104 K and ne = 104 cm−3.

6 Sources of interest

In Sect. 2.2 we mentioned six sources where the RRL spectral profile deviates significantly from the expected Gaussian line shape due to the presence of excess emission seen at higher velocities. We show the stacked mm-RRL spectral profiles for these sources in Fig. 11. In all of these cases, fitting a single Gaussian profile to the RRL results in large residuals. However, the large residuals can be significantly reduced using a two-component fit to these data. One component is centered at the systemic velocity of the sources and a second, broader component is fitted to the emission from the higher velocity gas. The parameters of the Gaussian fits of the two components are given in Table 6. In this section, we have a more detailed look at these sources to understand the nature of these additional high-velocity components. Although the non-Gaussian profiles are observed in both the submm- and mm-RRLs, we have chosen to fit the mm-RRLs as they tend to have a higher S/N and provide therefore more reliable results.

To facilitate the discussion of these sources we complement the spectra presented in Fig. 11 with three-color mid-infrared images that show the local environment of these regions in Fig. 12. It is interesting to note that all of the sources where a high velocity component has been identified are located in the first quadrant. However, these sources were all targeted in the observing campaign with the IRAM 30 m telescope, resulting in higher sensitivity observations than the observations of southern targets with the Mopra 22 m telescope, which therefore might have missed the weaker, broader components in many other sources. We also note that three of the sources show red-shifted emission, two show blue-shifted emission, and one appears to have both blue and red-shifted components. The higher-velocity gas is likely due to either an expanding flow away from the H II regions (Tenorio-Tagle 1979; Tenorio-Tagle et al. 1979) or bow-shock generated by wind-blowing high-mass stars moving at a supersonic velocity (van Buren et al. 1990; van Buren & Mac Low 1992) or possibly the combination of both phenomena.

6.1 AGAL012.804−00.199

This H II region is part of the W33 star-forming complex, and so we have adopted a distance of 2.6 kpc to this clump (Immer et al. 2012). The radio emission consists of two distinct sources, a bipolar structure located towards the center of the clump and a more compact, roughly circular source located to the north. All of the radio emission is correlated with bright mid-infrared emission seen in the GLIMPSE image (see the upper left panel of Fig. 12). The position of the IRAM beam is slightly offset from the peak of the bipolar radio emission, but the brightest lobe of bipolar structure falls within the FWHM beam size. However, emission from the weaker bipolar lobe and the compact source to the north are also likely to be picked up in the outer part of the beam.

The mm-RRL spectrum is shown in the upper left panel of Fig. 11; this reveals the presence of a red wing, which might be associated with an ionized flow perhaps arising from one of the bipolar lobes. The systemic velocity of the clump is 35.4 km s⁻¹, which is broadly consistent with the velocity of the strongest RRL component and so it is the weaker and broader component that is associated with the faster-moving gas. Bieging et al. (1978) also reported different velocity RRL components at low-frequency (1.7, 5, and 8.6 GHz) using the Effelsberg 100 m telescope. They suggested that the components with different velocities are associated with flows of ionized gas, which is consistent with our interpretation.

6.2 AGAL029.954−00.016

This source is associated with W43 South at a distance of 5.2 kpc (Reid et al. 2014) and has a systemic velocity of 97.2 km s⁻¹. The IRAM beam is centered on the peak of the dust emission, and the radio emission is slightly offset to the west; however, it is fully covered by the beam. The radio emission has a cometary morphology with a bright head pointing towards the center of the clump and a tail-like structure that extends away from the center of the clump into the lower-column density region. The radio emission is correlated both spatially and morphologically with the strong mid-infrared emission (see the upper middle panel of Fig. 12). The clump has a mass of ~ 4 × 10³ M_⊙ and a dust temperature of 35.5 K and so feedback from the H II region is having a significant impact on the physical properties of its host clump.

The RRL profile is shown in the upper-middle panel of Fig. 11 and reveals the presence of a high-velocity blue-shifted component; the peak of this component is offset by ~15–20 km s⁻¹ from the systemic velocity of the H II region and the host clump. The high-velocity blue-shifted component corresponds to the velocity range of a Brγ line (around 70 ~90 km s⁻¹) at the position of the head (Martín-Hernández et al. 2003). The linewidth of a Brγ toward the head of the H II region is also found to be ~20 km s⁻¹ broader than that from the tail, but their Brγ observations could not explain it with either bow-shock or champagne flow models. In addition, they also showed the existence of a velocity gradient toward the tail (Martín-Hernández et al. 2003).

Zhu et al. (2008), however, showed that [Ne II] line emission and its velocity maps toward this source well agreed with a model of stellar-wind pressure-driven H II regions with a viewing angle of 135°. In this model, the cometary morphology of the H II region is caused by either relatively slow motion of the star through the molecular cloud or a density gradient in the cloud (Zhu et al. 2008; Arthur & Hoare 2006). According to the model, emission lines from the ionized gas of a cometary head show a broad linewidth with a negative velocity offset concerning the velocity of molecular gas (Zhu et al. 2008). In fact, our broad RRL also has a negative velocity relative to the systemic velocity although it is not clear whether the broad RRL comes from the head of this source. Besides, there is a hot core in front of the UC H II region, and HCO⁺ and SiO molecular lines (Maxia et al. 2001) show the line densities are higher toward the head than the tail. Therefore, our broader RRL component is likely to be associated with accelerated ionized gas toward the tail along the shell created by the relatively slow motion of the star through the the interstellar medium.

Fig. 11 Green and yellow filled areas show two Gaussian fitted profiles of the Hnα transition. The back dashed line is a residual noise signal of the two Gaussian profiles. On the other hand, the red dashed line indicates a residual noise from a one-component Gaussian fit to the data. The vertical cyan dashed line shows the systemic velocity of the clump measured by H13 CO+ (1−0).

Fig. 12 GLIMPSE IRAC three-color composite images (blue; 3.6 μm, green; 4.6 μm, and red; 8 μm). The gray contours trace the 870 μm dust continuum emission from the ATLASGAL survey. The red contours show the 5 GHz radio continuum emission from the CORNISH survey. The cyan dashed circles are centered on the pointing positions of the IRAM observation with FWHM of 29″.

6.3 AGAL034.258+00.154

This ATLASGAL clump is associated with the G34.26+0.15 Complex. The radio emission consists of three radio continuum components; an extended source that has a cometary morphology and two compact regions to the west and north-west; these are referred to as C, A, and B, respectively (following Sewiło et al. 2011). The clump has a mass of ~1500 M_⊙ and dust temperature of ~30 K. The mid-infrared image (upper right panel of Fig. 12) reveals that this region is located at the edge of a larger H II region located to the south-west (region of extended mid-infrared emission). Besides, the morphologies of radio and mid-IR emissions are not so well correlated with each other toward this source. The reason is likely to be extinction (Hoare et al. 2007).

The RRL profile shows a strong blue-shifted wing and is fitted with two components with similar intensities but different linewidths. The velocity of the narrowest component corresponds to the systemic velocity of the clump, but it is offset by ~15 km s⁻¹ from the velocity of the broader component. There is a possibility that some of the RRL emission is coming from the larger H II region located to the south-west given that it falls within the outer part of the IRAM beam. This complex has been observed by Sewiło et al. (2011) at higher resolution in continuum and the H53α RRL. Their study separated the emission from the three radio sources and determined the velocity of A and C to be ~50 km s⁻¹ and that of B to be ~65 km s⁻¹. All of these are broadly consistent with the systemic velocity of the clump and the narrower RRL component. However, we note that source C has a very broad linewidth (49 km s⁻¹) and so may be responsible for the excess seen in the IRAM RRL spectra. Jaffe & Martín-Pintado (1999) also observed 3-mm hydrogen RRLs towards this source with the IRAM 30 m telescope and also reported the presence of a similar board linewidth high-velocity component, which is associated with source C.

It is currently unclear from the present data whether the blue-shifted wing seen in the IRAM spectra is the result of contamination from the nearby extended H II region or a combination of the three radio sources identified; however, the latter seems more likely. Furthermore, the cometary shape of source C has been well explained by the bow-shock model created by wind-blowing massive stars moving at a supersonic velocity with respect to the molecular gas (van Buren et al. 1990).

6.4 W49A: AGAL043.164−00.029 and AGAL043.166+00.011

These two H II regions are located in the W49A star-forming complex. It is located at a distance of ~11 kpc (Zhang et al. 2013) and is one of the two most active star formation regions in the Galaxy (Urquhart et al. 2014b). The region AGAL043.164−00.029 has a clump mass of 3 × 10⁴ M_⊙, luminosity of 1.6 × 10⁶ L_⊙ and systemic velocity of 14.3 km s⁻¹ while AGAL043.166+00.011 has a clump mass of 11 × 10⁴ M_⊙, luminosity of 8 × 10⁶ L_⊙ and systemic velocity of 2.9 km s⁻¹.

6.4.1 AGAL043.164−00.029

The radio emission associated with AGAL043.164−00.029 (lower left panel of Fig. 12) reveals a cometary morphology that is coincident with strong mid-infrared emission; this is known as W49A South. The tail of the cometary structure extends to the west away from the peak submillimeter emission. In front of the head of the cometary H II region, there is another very weak and compact H II region; this is referred to as W49A South-1 (De Pree et al. 1997). The clump looks relatively isolated, so all of the RRLs likely arise from the H II regions. The velocity of the strongest RRL component is correlated with the systemic velocity of the clump. The second component is red-shifted with a velocity 12 km s⁻¹ larger than the systemic velocity.

This red-shifted component is likely to be associated with high-velocity ionized gas that forms the tail, and that is expanding into the lower-density region away from the head of the cometary structure; the head itself is coincident with the highest column density partof the clump. According to the RRL results of De Pree et al. (1997), however, the velocity gradient over the cometary source and a broad linewidth toward the head supports the bow-shock model and a mass-loaded stellar wind in a molecular cloud with a density gradient, instead of the champagne flow model (De Pree et al. 1997; Zhu et al. 2008). Nevertheless, their RRL velocity resolution was too poor to allow for a detailed investigation of the velocity structure of this source.

6.4.2 AGAL043.166+00.011

The region AGAL043.166+00.011 is the most massive and active clump in the Galaxy; it is known as W49 North. It is perhaps the only example of a young massive protocluster outside of the Galactic center region (Urquhart et al. 2018). The radio emission reveals a cluster of compact sources and this has previously been identified as being associated with the highest density of embedded H II regions in the Galaxy (De Pree et al. 1997, 2004; Urquhart et al. 2013b). There is a good correlation between the compact sources seen in the mid-infrared image (lower middle panel of Fig. 12), however, there are few compact radio sources that do not have an infrared counterpart, perhaps due to high extinction.

The RRL spectrum reveals the presence of a relatively modest red-shifted wing that is offset by ~30 km s⁻¹ from the systemic velocity of the clump. Jaffe & Martín-Pintado (1999) also reported a broad mm-RRL component with an offset velocity from the systemic velocity towards this source. Given the large number of compact radio sources and the poor sensitivity to extended radio emission (due to the large heliocentric distance) it is not possible to investigate the nature of the high-velocity ionized gas associated with this region. Nevertheless, the 7-mm observations by De Pree et al. (2004) showed that the W49 North region has several broad RRL sources and these may be associated with signatures of ionized disk winds. However, it is not clear whether the broader component of our RRL is related to the broad RRLs of the 7-mm observations due to different spatial and flux sensitivities between both observations.

6.4.3 AGAL045.121+00.131

This clump is located at a distance of 8 kpc (Reid et al. 2014) and has a systemic velocity of 58.7 km s⁻¹. The clump has a mass of ~7000 M_⊙ and a luminosity of 1 × 10⁶ L_⊙. The radio emission appears to have a bipolar structure extending to the north and south. The radio morphology is mirrored in the mid-infrared image (lower right panel of Fig. 12). We note that discrete compact radio sources do appear in the radio emission instead of one bipolar structure and it is unclear if these indicate the presence of individual H II regions or localized enhancements in the radio emission where the ionization front is impacting on a dense structure in the molecular gas.

The RRL profile shows evidence for both red- and blue-shifted wings, although these are rather modest. The two components fitted to theRRL have similar velocities, both of which are correlated with the systemic velocity of the clump. The broader component is almost twice the linewidth of the narrow component and traces the high-velocity gas. The morphological distribution of the CORNISH continuum emission is correlated with the integrated [Ne II] line emission map presented by Zhu et al. (2008), but the [Ne II] observations did not fully resolve them. The 15.0–92.9 kms⁻¹ velocity channel maps of the [Ne II] emission to the north of this emission correspond to the most extended radio continuum and the mid-infrared source (see the lower right panel of Fig. 12). On the other hand, the southern radio continuum emission structures correspond to the [Ne II] emission seen in the 36.2–71.6 kms⁻¹ velocity channels. Comparing the radio continuum emission with the [Ne II] velocity channel maps may indicate that the two Gaussian components do not arise from a bipolar flow, but may instead be the result of discrete radio sources sharing a similar peak velocity. High-resolution observations are required to obtain better understanding of the structure and kinematics of this region.

6.5 Identification of an ALMA sample

Despite many previous radio continuum observations, the morphology of many H II regions still cannot be fully explained. Studying both the kinematics of ionized gas and the geometry of H II regions provides some insight into their nature. However, the morphology and RRL profiles change with the orientation of our viewing angle and therefore detailed comparison of the observationswith current theoretical models (Zhu et al. 2015; Steggles et al. 2017) are also required to help us understand the different morphologies observed.

Given that most compact H II regions detected have sizes of a few arc-second, high-resolution observations are required. Detailed studies have been previously conducted at cm-wavelengths using the Very Large Array (VLA) and AustraliaTelescope Compact Array (ATCA; Urquhart et al. 2007, 2009; Purcell et al. 2013); however, cm-RRLs can be significantly affected by pressure broadening and the H II regions can often be optically thick resulting in poorly constrained physical parameters. The Atacama Large Millimeter/submillimeter Array (ALMA) not only opens up the possibility to observe mm- and submm-RRLs at a similar resolution to the cm-RRLs observed by the VLA, but has several additional advantages: 1) the line to continuum ratio increases with frequency and so mm- and submm-RRLs are significantly brighter than at cm-wavelengths; 2) pressure broadening is reduced and is in most cases negligible; 3) the continuum emission of the H II regions is optically thin; 4) there are many thermal molecular transitions that allow the ionized and molecular gas to be traced simultaneously and at a similar resolution.

As previously shown, our sample consists of sources of the most massive clumps and most luminous H II regions in the Galaxy, and with the radio continuum and mm- and submm-RRL data we can identify a well-characterized sample for follow-up with ALMA. Our selection criteria are:

• 1.

  The H II regions should be resolved at high-resolution cm-continuum wavelengths (~2–3′′) but at resolutions smaller than the ALMA primary beam at 850 μm (i.e., <18′′). This ensures that we can investigate the kinematic structure of ionized gas over many beams and the source can be observed in a single pointed observation at the highest frequency.

• 2.

  The mm- and submm-RRLs should be bright and relatively isolated from contamination from thermal lines, which can be a significant issue at submillimeter wavelengths.

• 3.

  The H II regions should be relatively isolated from nearby evolved H II regions, which might contribute emission through the side lobes and make imaging difficult.

Applying these criteria to the sample discussed in this paper we have identified 17 candidates for further detailed studies with ALMA (see Table 7 and Fig. B.1). We excluded several sources that satisfied criteria 2 and 3 above since good quality radio continuum data are not available for them. In Table 7 we provide the coordinates of 6 cm radio sources instead of the ATLASGAL dust clumps. Our list of potential ALMA candidates covers radio continuum sources with various radio properties and morphologies. The cometary morphology is the most common in the list, but there is some uncertainty in their morphological classification due to limited sensitivity to angular scales.

These selected sources have clean RRLs profiles at both submm/mm-wavelengths, and thus both wavelengths can be used for observations. Given that the H II regions are optically thin, this means that both mm- and submm-observations are equally able to probe the properties of the ionized gas. As seen in Table 7, their submm-RRLs show brighter peak fluxes than their mm-RRLs. In fact, the line to continuum ratio at submm-wavelengths is a factor of two better at mm-wavelengths (see the lower panel of Fig. 4, and Figs. 3 and 4 in Peters et al. 2012). Observations at submm-wavelengths will also have a factor of three higher resolution compared to mm-observations for the same configuration, however, the primary beam is also three times smaller, and so mm-observations may be more suitable for more extended H II regions. The choice of what frequency at which to observe a particular source can, therefore, be tailored to ensure the region can be covered in a single pointing and that the same physical scales are probed (i.e., mm-wavelengths for near by and extended H II regions and submm-wavelengths for more distant and compact H II regions).

7 Summary and conclusions

We present the results of a set of targeted submillimeter RRL observations of 104 compact H II regions previously identified by Kim et al. (2017). We have observed the H25α, H27α, H28α, H29α, H30α, and H35β RRL transitions. The detection rates for the Hnα and Hnβ transitions are 89% (93) and 40% (34), respectively. These have been compared with the mm-recombination lines reported by Kim et al. (2017) to investigate the physical properties of the ionized gas and its local environment.

We find that the submm-RRLs are approximately a factor of two brighter than the mm-RRLs, which is consistent with LTE, and that the RRL line emission from the H II regions is optically thin. We also find a strong correlation between the velocities of the molecular gas and the ionized gas of the H II region still being embedded in their host cloud.

We find a significant correlation between the bolometric luminosity of the embedded clusters in their dust clumps and the ionizing photon flux produced by the stars. This indicates that the stars of the H II regions provide the majority of the clumps’ luminosity. We also find a trend for increasing dust temperature with the bolometric luminosity consistent with the feedback from the embedded H II regions on its natal environment. Besides, the luminosity of the submm-RRLs presents robust correlations with both bolometric luminosity and Lyman photon flux. Comparing the physical properties of the H II regions and their host clumps to the general population of dust clumps located in the inner-Galaxy, we find the clumps with the submm-RRL detections to be some of the most massive and luminous in the Galaxy.

We used the submm-RRL flux to measure the photoionizing production rate, Q, which corresponds to the Lyman photon flux. The measured Q(H28α) values for H28α lines are consistent with the Lyman photon flux measured from 5 GHz radio continuum emission. This confirms the Q as a tool for measuring SFRs using submm-RRL observations in external galaxies, without the problems of dust obscuration and influences of evolved star and synchrotron emission on continuum fluxes, encountered by other methods to measure SFRs.

A single Gaussian component well describes the profiles of the majority of RRLs. However, we have identified six H II regions where the spectral profiles show the presence of red and/or blue-shifted wings. These wings reveal the presence of a high-velocity component in the ionized gas. We investigate these sources in detail and find that the data towards four H II regions are consistent with the presence of high-velocity flows associated with the head of cometary or bipolar H II regions. The other two clumps are associated with clusters of compact H II regions, and so the observed RRL profiles may result from the blending of the emission from several H II regions. Higher-resolution observations are required to study the kinematics and physical properties of all of the regions in detail and to disentangle their morphology and long lifetime (about 10⁵ years).

Finally, we identify a sample of submm-RRL sources that would be suitable for high-resolution follow-up observations with ALMA. For the potential ALMA candidates, 17 H II regions with mm/submm-RRL detections are selected with the following criteria: appropriate sizes for the ALMA beam; bright mm/submm-RRLs with no contamination by molecular lines; relatively isolated location of the H II regions from nearby evolved H II regions.

Acknowledgements

We would like to thank the referee for their constructive comments and suggestions that have helped to improve thispaper. This paper is based on data acquired with the Atacama Pathfinder EXperiment (APEX). APEX is a collaboration between the Max Planck Institute for Radioastronomy, the European Southern Observatory, and the Onsala Space Observatory. This work was partly carried out within the Collaborative Research Council 956, sub-project A6, funded by the Deutsche Forschungsgemeinschaft (DFG). This document was produced using the Overleaf web application, which can be found at www.overleaf.com. Won-Ju Kim was supported for this research through a stipend from the International Max Planck Research School (IMPRS) for Astronomy and Astrophysics at the Universities of Bonn and Cologne. The ATLASGAL project is a collaboration between the Max-Planck-Gesellschaft, the European Southern Observatory (ESO) and the Universidad de Chile. It includes projects E-181.C-0885, E-078.F-9040(A), M-079.C-9501(A), M-081.C-9501(A) plus Chilean data.

Appendix A Continued table of list of observed sources

Appendix B Potential ALMA candidates

Fig. B.1 GLIMPSE IRAC three-color composite image (blue: 3.6 μm, green: 4.6 μm, and red: 8 μm). Gray contour represents 870 μm dust continuum emission from the ATLASGAL survey. Cyan and blue contours show 5 GHz radio continuum emission from CORNISH and RMS surveys. The arbitrary beam of 20″ is indicated with a bright green hatched circle in the left-bottom corner.

References

All Tables

All Figures

	Fig. 1 Example of the identification of molecular lines blended with submm-RRL emission (black line) in AGAL034.258+00.154. The red profile shows a combined profile of the modeled profiles of CH3 OCHO lines using WEEDS and a Gaussian profile of the H25α RRL. The latter is also shown as the green dashed line. Other emission lines are unidentified molecular lines.
In the text

Fig. 2 Examples of detected submm-RRLs (upper panels) of AGAL037.874−00.399 and AGAL332.156−00.449 and their corresponding mm-RRLs (lower panels) from Kim et al. (2017). The vertical purple dashed lines indicate the systemic velocity measured with H13CO+ (1−0) lines. The different colors present different submm/mm-RRL transitions detected in this paper and Kim et al. (2017). In flux unit, 1 Jy corresponds to ~0.025 K for the APEX 12 m, 0.169 K for the IRAM 30 m, and 0.045 K for the Mopra 22 m data.

In the text

Fig. 3 Upper panel: peak flux of the H28α transition versus peak flux of other transitions; H25α (black), H26α (red), H27α (cyan), H29α (orange), and H30α (bright green). All the plotted sources are detected in the H28α transition, and thus there are a few sources missing, which are only detected in other lines. Lower panel: peak flux of the H28α transition versus peak flux of H35β. In both plots, equal fluxes are indicated by the black dashed line, and blue lines indicate the best-fit determined by BCES bisector fits to all data points. The red line in the lower panel presents the predicted LTE value.

In the text

Fig. 4 Upper panel: linewidth comparison of stacked mm- and submm-RRLs. The black filled circles indicate sources with a peak flux with S/N ≥ 8 σ of both RRLs. Lower panel: peak flux comparison of stacked mm- and submm-RRL. The black symbols are the same as in the upper panel. The black dashed lines indicate the locus of equal linewidths and fluxes. The blue lines show the best fit from BCES bisector fitting.

In the text

Fig. 5 Bolometric luminosity (Lbol) as a function of Lyman continuum photon flux of 6 cm radio continuum emission integrated over the beam of the 3 mm observations. The colors indicate dust temperatures in the molecular clumps. The blue line shows the result of a fit to the datapoints using BCES bisector fitting. The purple squares superposed on the circles are particular sources with two Gaussian component RRLs (see Table 6 and Sect. 6). The red dotted line and red stars indicate bolometric luminosity and Lyman continuum photon flux corresponding to each spectral type (O4 − B0.5) for ZAMS stars from Panagia (1973).

In the text

	Fig. 6 GLIMPSE IRAC three-color composite image (blue; 3.6 μm, green; 4.6 μm, and red; 8 μm) of one (AGAL345.408−00.952) of the three outliers in Fig. 5. The gray contours represent the 870 μm dust continuum emission from the ATLASGAL survey. The white circle indicates the size of Mopra beam (FWHM, 36″) and the red circle presents a RMS radio continuum source (Urquhart et al. 2007).
In the text

	Fig. 7 Luminosity of the submm-RRLs versus Lbol of embedded central objects. The color bar indicates the Lyman photon flux of ionizing stars associated with the dust clumps and the submm-RRLs. The best linear fit is presented with the blue line.
In the text

	Fig. 8 Bolometric luminosity (Lbol) as a functionof Mclump toward the full ATLASGAL sample (Urquhart et al. 2018) that indicates both non-detected and unobserved sources by this submm-RRL survey. The bright green circles indicate sources with a submm-RRL detection.
In the text

	Fig. 9 Histogram of Lbol∕Mclump for the full ATLASGAL sources (gray) and the submm-RRL sources (bright green).
In the text

Fig. 10 Photoionzing photon rate (Q) measured by the submm-RRLs (this paper) as a function of Lyman photon flux estimated by 5 GHz radio continuum emission from CORNISH and RMS surveys (Urquhart et al. 2007; Purcell et al. 2013). Those two parameters have the same unit (s−1), the number of photons per time. The equality of both axes is indicated with a black dashed line. The color dotted and dashed lines present the best linear fits to the data points calculated for different electron temperature and density in Eq. (1). The black data points are estimated with Te = 104 K and ne = 104 cm−3.

In the text

	Fig. 11 Green and yellow filled areas show two Gaussian fitted profiles of the Hnα transition. The back dashed line is a residual noise signal of the two Gaussian profiles. On the other hand, the red dashed line indicates a residual noise from a one-component Gaussian fit to the data. The vertical cyan dashed line shows the systemic velocity of the clump measured by H13 CO+ (1−0).
In the text

	Fig. 12 GLIMPSE IRAC three-color composite images (blue; 3.6 μm, green; 4.6 μm, and red; 8 μm). The gray contours trace the 870 μm dust continuum emission from the ATLASGAL survey. The red contours show the 5 GHz radio continuum emission from the CORNISH survey. The cyan dashed circles are centered on the pointing positions of the IRAM observation with FWHM of 29″.
In the text

	Fig. B.1 GLIMPSE IRAC three-color composite image (blue: 3.6 μm, green: 4.6 μm, and red: 8 μm). Gray contour represents 870 μm dust continuum emission from the ATLASGAL survey. Cyan and blue contours show 5 GHz radio continuum emission from CORNISH and RMS surveys. The arbitrary beam of 20″ is indicated with a bright green hatched circle in the left-bottom corner.
In the text